SUBROUTINE plm (lmax, mdiff, ctheta, plmn) USE Numeric_Kinds_Module ! computes normalized associated legendre polynomials (plmn) by recursio ! from l = 0 to l = lmax for some fixed mdiff. we use abramowitz and st ! recurrance relation 8.5.3. notice that this relation fixes m and recu ! over l. our method is to use the assoc. legendre poly. with m = mdiff ! and l = mdiff as the starting point in the recurrance. the next term ! the phase is that of condon and shortley, abramowitz and stegun, ! and gradshteyn and ryzhik. ! then each plmn is normalized to 1.0 on (-1.,1.). ! for either positive or negative mdiff, multiplication by exp(1*mdiff*p ! sqrt(2.*pi) gives the spherical harmonics of condon and shortley, rose ! messiah, and jackson. ! this routine assumes that the argument really is the cosine of the ang ! variable uses ! name |dim | use !----------------------------------------------------------------------- ! ctheta r dp cosine of the angle theta, argument of the polynomials ! i i do loop counter ! j i do loop counter ! lmax i the maxium value of l that we will use ! m i the absolute value of mdiff ! mdiff i the difference of the values of m in the calling program ! part* r dp used these to break the recurrance relation into 4 pieces ! plmn r dp array that contains the normalized associated legendre ! polynomials from l=0 to l=lmax for m=mdiff ! prfac1 r dp prefactor for calculating the polynomial for l = m ! prfac2 r dp prefactor for calculating the normalization ! start i starting value of l for the recurrance relation ! stheta r dp the sine of the angle theta IMPLICIT NONE ! this is a vax only statement INTEGER :: i, j, lmax, m, mdiff, start REAL(dp) :: ctheta, part1, part2, part3, part4, plmn (0:lmax), & prfac1, prfac2, stheta m = abs (mdiff) ! check to see if the indicies are proper IF(m.gt.lmax.or.lmax.lt.0)THEN WRITE(6, * ) 'illegal values of indices in plm' STOP ENDIF ! zero out all the terms where l 0 than for m < 0. reference: arfken pp. 568, exercise 12.5.4 ! and abramowitz & segun 8.6.17 respectively. ! initialize prfac1 for m=0. prfac1 = 1.0d0 IF(mdiff.lt.0)THEN do 20 i = 1, m prfac1 = dble (i) * prfac1 20 ENDDO prfac1 = 1.0d0 / prfac1 ELSEIF (mdiff.gt.0)THEN do 30 i = m + 1, 2 * m prfac1 = dble (i) * prfac1 30 ENDDO ! this takes care of the phase of condon and shortley prfac1 = prfac1 * ( ( - 1.0d0) **m) ENDIF ! now use the prefactor to calculate plmn for m = l stheta = dsqrt (1.0d0 - ctheta * ctheta) IF(stheta.eq.0.0d0)THEN IF(m.eq.0)THEN plmn (0) = 1.0d0 ELSE ! reference: arfken 12.94 plmn (m) = 0.0d0 ENDIF ELSE plmn (m) = prfac1 * (stheta**m) / (2.0d0**m) ENDIF ! set up for calculating plmn ! if m=0 we are calculating legendre polynomials and we must ! provide the plmn(1) and start recurring at plmn(2). IF(m.eq.0.and.lmax.ne.0)THEN plmn (1) = ctheta start = 2 ! WRITE(6,*)'ctheta=',ctheta ELSE ! otherwise if m not equal to 0, start at m+1 start = m + 1 ENDIF ! calculate plmn via the recurrance relation up to lmax do 50 i = start, lmax part1 = 2.0d0 * dble (i) * ctheta * plmn (i - 1) part2 = (dble (i + mdiff) ) * plmn (i - 2) part3 = ctheta * plmn (i - 1) part4 = plmn (i - 2) plmn (i) = (part1 - part2 - (part3 - part4) ) / (dble (i - & mdiff) ) 50 ENDDO ! normalization; as from abramowitz & segun pp. 332 or ! arfken 12.148 do 60 i = m, lmax prfac2 = 1.0d0 ! compute the prefactor (prfac2); notice that if m=0 this do loop ! is skipped do 65 j = i - m + 1, i + m prfac2 = prfac2 * dble (j) 65 ENDDO IF(mdiff.le.0)THEN plmn (i) = dsqrt (prfac2 * (2.0d0 * dble (i) + 1.0d0) & / 2.0d0) * plmn (i) ELSEIF (mdiff.gt.0)THEN plmn (i) = dsqrt ( (2.0d0 * dble (i) + 1.0d0) / (2.0d0 * & prfac2) ) * plmn (i) ENDIF 60 ENDDO RETURN END SUBROUTINE plm