##// END OF EJS Templates
DATABASES » madrigal
RSS
Add restart to nginx container
Add restart to nginx container
jespinoza - - Load All Authors

File last commit:

r0:b84e1135c2c4
r22:1c8c96a4f254 master
Show More
rylm.f
250 lines | 7.2 KiB | text/x-fortran | FortranFixedLexer
/ source / madf / geolib / rylm.f
History | Source | Raw |Copy content |Copy permalink
jespinoza
Initial
 r0 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