Program hw2_1 implicit none * Program to produce a table Legrendre Polynomials and Associated Legrenddre * Functions. integer*4 max_arg ! Maximum number of arguments parameter ( max_arg = 101 ) * MAIN PROGRAM VARIABLES integer*4 ierr ! IOSTAT error in case there is a problem writing ! to screen integer*4 i, n,m,k ! Loop variables integer*4 narg ! Number of arguments to be computed real*8 darg ! Spacing between argument values real*8 eps ! Small change from -1 and +1 real*8 arg(max_arg) ! Arguments to be evaluated real*8 pnm(max_arg) ! Polynomial and associated functions real*8 plm, plgndr ! Function to compute polynomial and assocated function character*80 fmt_arg **** Write the header lines out write(*,120,iostat=ierr) 120 format('TABLE OF LEGENDRE POLYNOMIALS AND ASSOCIATED FUNCTIONS',/, . '------------------------------------------------------') **** Output header line with the arguments to be used i = 0 narg = 9 eps = 0 darg = 2*(1-eps)/(narg-1) do i = 1,9 arg(i) = -1.d0 + eps + (i-1)*darg end do write(fmt_arg, 140) narg 140 format('("N,M\\Arg ",',I3,'F8.2)') write(*,fmt_arg) (arg(i),i=1,narg) if( ierr.ne.0 ) then print *,'IOSTAT Error ',ierr,' writing header' stop 'HW2_1: Error in writing to screen' end if do n = 0,4 do m = 0,n do k = 1, narg c pnm(k) = plm(n,m,arg(k)) pnm(k) = plgndr(n,m, arg(k)) enddo write(*,160) n,m, (pnm(k),k=1,narg) 160 format(i3,1x,i3,1x,50(F8.3)) end do end do write(*,220,iostat=ierr) 220 format(/,'TABLE OF LEGENDRE POLYNOMIALS AND ASSOCIATED FUNCTIONS', . /,'ALTERNATIVE RECURRSION ALGORITHM' . /,'------------------------------------------------------') write(*,fmt_arg) (arg(i),i=1,narg) do n = 0,4 do m = 0,n do k = 1, narg * Other version that has problems at -1 and +1 pnm(k) = plm(n,m,arg(k)) enddo write(*,160) n,m, (pnm(k),k=1,narg) end do end do end CTITLE PLM real*8 function plm(n,m,z) c this function returns associated Legendre functions of c degree n and order m. * The recursion relationship here computes the associated * functions from earlier ones with the same degree. (The * version on the Wolfram sites, keeps the order the same). * The algorithm here has problems with -1 and +1 arguments. * m - order of function * mi - order of current recursive calculation * n - degree of function * ni - degree of current recursive calculation integer*4 m, mi, n, ni * z - argument (ie cos(theta) ) * p_new - work variable * pn(3) - calculation vector * pm(2) - calculation vector real*8 z, p_new, pn(3), pm(2) c Initialize recursive calculation pn(3) = z pn(2) = 1.d0 pn(1) = 0.d0 mi = 0 ni = 1 c.... Calculate Legendre functions of order m do while (mi .lt. m) do while (ni .le. mi+2) c Increase degree by 1 p_new = ((2*ni + 1)*z*pn(3) - (ni + mi)*pn(2))/(ni-mi+1) pn(1) = pn(2) pn(2) = pn(3) pn(3) = p_new ni = ni + 1 end do c Increase order by 1 ni = ni - 2 pm(1) = ((ni+mi+1)*z*pn(1) - . (ni-mi+1)*pn(2))/sqrt(1.d0-z**2) ni = ni + 1 pm(2) = ((ni+mi+1)*z*pn(2) - . (ni-mi+1)*pn(3))/sqrt(1.d0-z**2) c Store Pm into Pn and repeat pn(3) = pm(2) pn(2) = pm(1) pn(1) = 0.d0 mi = mi + 1 end do c.... Advance degree to n if (ni .gt. n) then plm = pn(2) else do while (ni .lt. n) P_new = ((2*ni+1)*z*pn(3) - (ni + mi)*pn(2))/(ni-mi+1) pn(2) = pn(3) pn(3) = p_new ni = ni + 1 end do plm = pn(3) end if return end CTITLE PLGNDR real*8 function plgndr(l,m,x) * Computes Legendre functions based on algorithm given in: * Press et al., Numerical Recipes, Page 182-183, * * Press, W. H.; Flannery, B. P.; Teukolsky, S. A.; and Vetterling, * W. T. Numerical Recipes in FORTRAN: The Art of Scientific Computing, * 2nd ed. Cambridge, England: Cambridge University Press, p. 252, 1992. * Code implemented here is based on the recursion algorithms in * http://mathworld.wolfram.com/LegendrePolynomial.html * Plm here is for the Abramowitz and Stegun 1972, Press et al. 1992 * definitions (geophysical form) * PASSED VARIABLES integer*4 l, m ! Degree and order to be evaluated real*8 x ! Argument. Value -1 to 1 * LOCAL VARIABLES real*8 pmm, ! Value op Pmm (direct expression . somx2, ! sqrt of (1-x^2) . fact, ! Factor to multiply odd integers . pll, ! Nominally Pll but takes on different values ! in recursion formula . pmmp1 ! P(m-1)m and other values in recurrion routine. integer*4 ll, i ! Counters used in recurrions calculation **** Start pmm = 1.d0 if( m.gt.0 ) then * Implements equation 68 of LegendrePolynomial.html. Note that * n!! means product of odd integers from 1 to n somx2 = sqrt((1.d0-x)*(1.d0+x)) ! Form thought to be faster than ! sqrt(1.d0-x**2) fact = 1.d0 do i = 1,m pmm = pmm*fact*somx2 fact = fact + 2.d0 end do end if **** See if degree == order if( l.eq.m ) then plgndr = pmm else * Implements equation 69 of LegendrePolynomial.html. Degree is one * more than order pmmp1 = x*(2*m+1)*pmm if( l.eq.m+1 ) then plgndr = pmmp1 else * Implements equation 65 of * http://mathworld.wolfram.com/LegendrePolynomial.html do ll = m+2,l pll = (x*(2*ll-1)*pmmp1-(ll+m-1)*pmm)/(ll-m) pmm = pmmp1 pmmp1 = pll end do plgndr = pll endif end if return end