Program for evaluating the Legendre function (test2d19.dat).
This can be copied from the end of file test2d19.dat, which is supplied with the CPO2D package.
C Program for evaluating the Legendre function Aug 1998.
C (written quickly, so not sophisticated)
OPEN (UNIT=1,FILE='tempin.dat',STATUS='OLD')
C =input data
OPEN (UNIT=2,FILE='tempres.dat',STATUS='UNKNOWN')
C =results file
pi = 4.*ATAN(1.)
READ (1,*) an , x
C =order, argument (ie cos(phi) )
READ (1,*) m
C =number of summation terms
phi = ACOS(x)
ans = SIN((an + 1.)*phi)
C =1st term in summation
coeff = (an + 1.)/(2.*an + 3.)
C =coefficient of 2nd term
ai = an + 3.
C -for argument of sin of 2nd term
DO i = 2 , m
ans = ans + coeff*SIN(ai*phi)
C numbers for next term:
f3 = real(2*i - 1)
f5 = f3 + 2.
f2 = real(i)
coeff = coeff*f3*(an + f2)/(f2*(2.*an + f5))
ai = ai + 2.
ENDDO
ans = ans*(2.**(an + 2.))/pi
C
C Evaluate the factorial factors:
C (do each one separately, so that they can be checked)
eps = 1.E-12
fln_fac = (an + 0.5)*log(an + 1.0) - (an + 1.) + alog(sqrt(2.*pi))
C -will be log of factorial(an) = gamma(an+1)
term_fln = 1.E10
real_m = 0.
C =summation index m
100 real_m = real_m + 1.
IF (term_fln.LT.eps) GOTO 300
C -ie summation over m finished
fm = 0.5*real_m/((real_m + 1.)*(real_m + 2.))
real_k = 0.
C =summation index k
sum_over_k = 0.
200 real_k = real_k + 1.
term_k = (an + 1.0 + real_k)**(-(real_m + 1.))
sum_over_k = sum_over_k + term_k
term_fln = fm*sum_over_k
IF (fm*term_k.LT.eps) THEN
C -ie summation over k finished
fln_fac = fln_fac + term_fln
GOTO 100
ELSE
GOTO 200
ENDIF
300 f_fac = exp(fln_fac)
C And then same procedure for the double factorial:
fln_fac2 = (an + 1.0)*log(an + 1.5) - (an + 1.5) +
& alog(sqrt(2.*pi))
C -will be log of double factorial (2*an+1)
C = gamma(an+1.5)*2**(an+1)/sqrt(pi)
term_fln = 1.E10
real_m = 0.
110 real_m = real_m + 1.
IF (term_fln.LT.eps) GOTO 310
fm = 0.5*real_m/((real_m + 1.)*(real_m + 2.))
real_k = 0.
sum_over_k = 0.
210 real_k = real_k + 1.
term_k = (an + 1.5 + real_k)**(-(real_m + 1.))
sum_over_k = sum_over_k + term_k
term_fln = fm*sum_over_k
IF (fm*term_k.LT.eps) THEN
fln_fac2 = fln_fac2 + term_fln
GOTO 110
ELSE
GOTO 210
ENDIF
310 f_fac2 = exp(fln_fac2)
f_fac2 = f_fac2*2.**(an+1)/sqrt(pi)
ans = ans*f_fac/f_fac2
WRITE (2,'(''no of terms n cos(phi) result''/
&I8,2X,1P3E14.6)') m , an , x , ans
STOP
END