|
|
SUBROUTINE RYLM(D_COLAT,D_LON,I_ORDER,D_YLMVAL)
|
|
|
c
|
|
|
c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
|
|
c
|
|
|
c $Revision: 1.1 $
|
|
|
c
|
|
|
c The initial version of this subroutine was written by RADEX, INC.
|
|
|
c for VAX/VMS systems.
|
|
|
c
|
|
|
c Subsequent revisions for use with POSIX compliant systems
|
|
|
c have been made by KBB at JHU/APL. These revisions have
|
|
|
c been managed using RCS, and a log of these revisions
|
|
|
c follows.
|
|
|
c
|
|
|
c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
|
|
c RCS Revision History
|
|
|
c
|
|
|
c $Log: not supported by cvs2svn $
|
|
|
c Revision 1.3 1996/03/11 19:25:07 baker
|
|
|
c Revisions for 1995 version of AACGM
|
|
|
c
|
|
|
c Revision 1.2 1994/10/14 10:50:45 baker
|
|
|
c Added the SAVE instruction to make the variables static.
|
|
|
c
|
|
|
c Revision 1.1 94/10/12 15:24:21 15:24:21 baker (Kile Baker S1G)
|
|
|
c Initial revision
|
|
|
c
|
|
|
c
|
|
|
c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
|
|
C
|
|
|
C Purpose:
|
|
|
C
|
|
|
C This subroutine computes the array of real spherical harmonic
|
|
|
C function values Y(L,M) for a given colatitude and longitude up
|
|
|
C to order I_ORDER. The function values are stored in the one-
|
|
|
C dimensional array D_YLMVAL(N). The indexing scheme used is as
|
|
|
C follows:
|
|
|
C
|
|
|
C L 0 1 1 1 2 2 2 2 2 3 3 3 3 3 3 3 4
|
|
|
C M 0 -1 0 1 -2 -1 0 1 2 -3 -2 -1 0 1 2 3 -4 etc.
|
|
|
C
|
|
|
C N 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
|
|
|
C
|
|
|
C Input Arguments:
|
|
|
C
|
|
|
C D_COLAT - Double Precision - The colatitude of the
|
|
|
C point for which the spherical harmonic
|
|
|
C Y(L,M) will be calculated
|
|
|
C
|
|
|
C D_LON - Double Precision - The longitude of the
|
|
|
C point for which the spherical harmonic
|
|
|
C Y(L,M) will be calculated
|
|
|
C
|
|
|
C I_ORDER - Integer - The order of the spherical
|
|
|
C harmonic function expansion. The total
|
|
|
C number of terms computed will be
|
|
|
C (I_ORDER + 1) * (I_ORDER + 1)
|
|
|
C
|
|
|
C Output Argument:
|
|
|
C
|
|
|
C D_YLMVAL - Double Precision array of spherical harmonic
|
|
|
C functions at the point (D_COLAT, D_LON)
|
|
|
C
|
|
|
C Local Variables:
|
|
|
C
|
|
|
C Q_FAC - Double Precision Complex number used for
|
|
|
C recursive computation of the spherical
|
|
|
C harmonic function values. Represents the
|
|
|
C longitude dependence in the recursion
|
|
|
C computation.
|
|
|
C
|
|
|
C Q_VAL - Double Precision Complex number used to
|
|
|
C store spherical harmonic function for
|
|
|
C current L and M values for recursive
|
|
|
C computation for the next set of L and M
|
|
|
C values
|
|
|
C
|
|
|
C COS_THETA - Double Precision Cosine of D_COLAT
|
|
|
C SIN_THETA - Double Precision Sine of D_COLAT
|
|
|
C used to compute Spherical Harmonic functions
|
|
|
C
|
|
|
C COS_LON - Double Precision Cosine of Longitude
|
|
|
C SIN_LON - Double Precision Sine of Longitude
|
|
|
C used to compute Spherical Harmonic functions
|
|
|
C
|
|
|
C CA, CB - Double Precision Coefficients used in the
|
|
|
C recursion computations
|
|
|
C
|
|
|
C LA, LB, LC, - Integer - Pointers for the D_YLMVAL array
|
|
|
C LD, LE, LF corresponding to values of Y(L,M) for
|
|
|
C particular values of L, M for use in the
|
|
|
C recursive computation of Y(L.M)
|
|
|
C
|
|
|
C
|
|
|
C Constants: None
|
|
|
C
|
|
|
C Subroutines Required: None
|
|
|
C
|
|
|
C Files Used at Compile Time: None
|
|
|
C
|
|
|
C Files Used at Run Time: None
|
|
|
C
|
|
|
C Databases Accessed: None
|
|
|
C
|
|
|
C Warnings: None
|
|
|
C
|
|
|
C Revision History:
|
|
|
C
|
|
|
C Written by Radex, Inc., 3 Preston Court, Bedford, MA 01730
|
|
|
C
|
|
|
C
|
|
|
c
|
|
|
C
|
|
|
C
|
|
|
C
|
|
|
C Local Variables:
|
|
|
C
|
|
|
C
|
|
|
C .. Scalar Arguments ..
|
|
|
DOUBLE PRECISION D_COLAT,D_LON
|
|
|
INTEGER I_ORDER
|
|
|
C ..
|
|
|
C .. Array Arguments ..
|
|
|
DOUBLE PRECISION D_YLMVAL(*)
|
|
|
C ..
|
|
|
C .. Local Scalars ..
|
|
|
DOUBLE COMPLEX Q_FAC,Q_VAL
|
|
|
DOUBLE PRECISION CA,CB,COS_LON,COS_THETA,FAC,SIN_LON,SIN_THETA
|
|
|
INTEGER L,LA,LB,LC,LD,LE,LF,M
|
|
|
C ..
|
|
|
C .. External Functions ..
|
|
|
C ..
|
|
|
C .. Intrinsic Functions ..
|
|
|
INTRINSIC DCMPLX,DCOS,DIMAG,DSIN,DBLE
|
|
|
C ..
|
|
|
C .. Save statement ..
|
|
|
SAVE
|
|
|
C ..
|
|
|
COS_THETA = DCOS(D_COLAT)
|
|
|
SIN_THETA = DSIN(D_COLAT)
|
|
|
C
|
|
|
COS_LON = DCOS(D_LON)
|
|
|
SIN_LON = DSIN(D_LON)
|
|
|
C
|
|
|
Q_FAC = -SIN_THETA*DCMPLX(COS_LON,SIN_LON)
|
|
|
C
|
|
|
C Generate Zonal Harmonics (Y(L,0), L = 1, I_ORDER)
|
|
|
C using Standard Recursion Relations (Equation 6.8.7, p 247,
|
|
|
C Numerical Recipes in Fortran, 2. edition, by Press. W. et. al.
|
|
|
C Cambridge University Press, 1992) for case for which m = 0
|
|
|
C
|
|
|
C
|
|
|
D_YLMVAL(1) = 1.0D0
|
|
|
D_YLMVAL(3) = COS_THETA
|
|
|
C
|
|
|
DO 10 L = 1,I_ORDER - 1
|
|
|
LA = (L-1)*L + 1
|
|
|
LB = L*(L+1) + 1
|
|
|
LC = (L+1)*(L+2) + 1
|
|
|
C
|
|
|
CA = (2.0D0*L+1.0D0)/(L+1.0D0)
|
|
|
CB = DBLE(L)/(L+1.0D0)
|
|
|
C
|
|
|
D_YLMVAL(LC) = CA*COS_THETA*D_YLMVAL(LB) - CB*D_YLMVAL(LA)
|
|
|
C
|
|
|
10 CONTINUE
|
|
|
C
|
|
|
C Generate Y(L,L) for L = 1 to I_ORDER. Algorithm based upon
|
|
|
C equation (6.8.8) in Press et al, but incorporates longitude
|
|
|
C dependence.
|
|
|
C
|
|
|
Q_VAL = Q_FAC
|
|
|
C
|
|
|
D_YLMVAL(4) = DBLE(Q_VAL)
|
|
|
D_YLMVAL(2) = -DIMAG(Q_VAL)
|
|
|
DO 20 L = 2,I_ORDER
|
|
|
C
|
|
|
Q_VAL = (2.0D0*L-1.0D0)*Q_FAC*Q_VAL
|
|
|
C
|
|
|
LA = L*L + 2*L + 1
|
|
|
LB = L*L + 1
|
|
|
C
|
|
|
D_YLMVAL(LA) = DBLE(Q_VAL)
|
|
|
D_YLMVAL(LB) = -DIMAG(Q_VAL)
|
|
|
C
|
|
|
20 CONTINUE
|
|
|
C
|
|
|
C Generate Y(L+1,L) to (Y(I_ORDER,L) function values.
|
|
|
C Algorithm based upon equation (6.8.9) in Press et al, but
|
|
|
C incorporates longitude dependence.
|
|
|
C
|
|
|
DO 30 L = 2,I_ORDER
|
|
|
C
|
|
|
LA = L*L
|
|
|
LB = L*L - 2*(L-1)
|
|
|
C
|
|
|
LC = L*L + 2*L
|
|
|
LD = L*L + 2
|
|
|
C
|
|
|
FAC = 2.0D0*L - 1.0D0
|
|
|
C
|
|
|
D_YLMVAL(LC) = FAC*COS_THETA*D_YLMVAL(LA)
|
|
|
D_YLMVAL(LD) = FAC*COS_THETA*D_YLMVAL(LB)
|
|
|
C
|
|
|
30 CONTINUE
|
|
|
C
|
|
|
C Generate remaining Y(L+2,L) to YL(I_ORDER,L) function values.
|
|
|
C Algorithm based upon equation (6.8.7) in Press et al, but
|
|
|
C incorporates longitude dependence.
|
|
|
C
|
|
|
DO 50 M = 1,I_ORDER - 2
|
|
|
C
|
|
|
C INITIALIZATION OF POINTERS FOR CURRENT VALUE OF M
|
|
|
C
|
|
|
LA = (M+1)**2
|
|
|
LB = (M+2)**2 - 1
|
|
|
LC = (M+3)**2 - 2
|
|
|
C
|
|
|
LD = LA - 2*M
|
|
|
LE = LB - 2*M
|
|
|
LF = LC - 2*M
|
|
|
C
|
|
|
DO 40 L = M + 2,I_ORDER
|
|
|
C
|
|
|
C
|
|
|
CA = DBLE(2*L-1)/DBLE(L-M)
|
|
|
CB = DBLE(L+M-1)/DBLE(L-M)
|
|
|
C
|
|
|
D_YLMVAL(LC) = CA*COS_THETA*D_YLMVAL(LB) - CB*D_YLMVAL(LA)
|
|
|
C
|
|
|
D_YLMVAL(LF) = CA*COS_THETA*D_YLMVAL(LE) - CB*D_YLMVAL(LD)
|
|
|
C
|
|
|
C UPDATE POINTERS FOR NEXT M VALUE
|
|
|
C
|
|
|
LA = LB
|
|
|
LB = LC
|
|
|
LC = LB + 2*L + 2
|
|
|
C
|
|
|
LD = LA - 2*M
|
|
|
LE = LB - 2*M
|
|
|
LF = LC - 2*M
|
|
|
C
|
|
|
40 CONTINUE
|
|
|
C
|
|
|
50 CONTINUE
|
|
|
C
|
|
|
C
|
|
|
C
|
|
|
RETURN
|
|
|
END
|
|
|
|
