subroutine BelPnmdt(pnm,dpt1,dpt2,maxn,t) !计算规格化Pnm及其对t一、二阶导数 !规格化Pnm采用改进的Belikov递推算法 !一、二阶导数非奇异递推算法 implicit none integer::maxn,n,m,k,kk,k1,k2,k3,astat(6) real*8::pnm((maxn+2)**2),dpt1((maxn+2)**2),dpt2((maxn+2)**2) real*8::t,u,a,b,c,d,e real*8::anm,bnm,cnm,da !--------------------------------------------------------------------------- pnm=0.d0;dpt1=0.d0;dpt2=0.d0 u=dsqrt(1-t**2);pnm(1)=1.d0 pnm(2)=dsqrt(3.d0)*t;pnm(3)=dsqrt(3.d0)*u do n=1,maxn a=dsqrt(2.d0*n+1.d0)/dsqrt(2.d0*n-1.d0) b=dsqrt(2.d0*(n-1.d0)*(2.d0*n+1.d0))/dsqrt((2.d0*n-1.d0)*n) kk=n*(n+1)/2+1;k1=n*(n-1)/2+1;k2=n*(n-1)/2+2 pnm(kk)=a*t*pnm(k1)-b*u/2.d0*pnm(k2) do m=1,n kk=n*(n+1)/2+m+1;k1=n*(n-1)/2+m+1 k2=n*(n-1)/2+m+2;k3=n*(n-1)/2+m c=a/n*dsqrt(dble(n**2-m**2));d=a/n/2.d0*dsqrt(dble((n-m)*(n-m-1))) e=a/n/2.d0*dsqrt(dble((n+m)*(n+m-1))) if(m==1)e=e*dsqrt(2.d0) pnm(kk)=c*t*pnm(k1)-d*u*pnm(k2)+e*u*pnm(k3) enddo enddo !计算一、二阶导数 dpt1=0.d0;dpt2=0.d0 dpt1(1)=0.d0;dpt1(2)=-dsqrt(3.d0)*u;dpt1(3)=dsqrt(3.d0)*t dpt2(1)=0.d0;dpt2(2)=-dsqrt(3.d0)*t;dpt2(3)=-dsqrt(3.d0)*u do n=2,maxn kk=n*(n+1)/2+1;k1=n*(n+1)/2+2 !m=0 dpt1(kk)=-dsqrt(dble(n*(n+1))/2.d0)*pnm(k1) kk=n*(n+1)/2+2;k1=n*(n+1)/2+1;k2=n*(n+1)/2+3 !m=1 anm=dsqrt(dble(2*n))*dsqrt(dble(n+1))/2.d0 bnm=-dsqrt(dble(n-1))*dsqrt(dble(n+2))/2.d0 dpt1(kk)=anm*pnm(k1)+bnm*pnm(k2) kk=n*(n+1)/2+1;k1=n*(n+1)/2+3 !m=0 anm=-dsqrt(dble(n*(n+1))/2.d0) bnm=dsqrt(dble(n*(n-1)))*dsqrt(dble((n+1)*(n+2)))/dsqrt(8.d0) dpt2(kk)=anm*pnm(kk)+bnm*pnm(k1) kk=n*(n+1)/2+2;k1=n*(n+1)/2+4 !m=1 anm=-dble(2*n*(n+1)+(n-1)*(n+2))/4.d0 bnm=dsqrt(dble((n-2)*(n-1)))*dsqrt(dble((n+2)*(n+3)))/4.d0 dpt2(kk)=anm*pnm(kk)+bnm*pnm(k1) do m=2,n kk=n*(n+1)/2+m+1;k1=n*(n+1)/2+m;k2=n*(n+1)/2+m+2 anm=dsqrt(dble(n+m))*dsqrt(dble(n-m+1))/2.d0 bnm=-dsqrt(dble(n-m))*dsqrt(dble(n+m+1))/2.d0 dpt1(kk)=anm*pnm(k1)+bnm*pnm(k2) k1=n*(n+1)/2+m-1;k2=n*(n+1)/2+m+3 anm=dsqrt(dble((n-m+1)*(n-m+2)))*dsqrt(dble((n+m-1)*(n+m)))/4.d0 bnm=-dble((n-m+1)*(n+m)+(n-m)*(n+m+1))/4.d0 cnm=dsqrt(dble((n-m-1)*(n-m)))*dsqrt(dble((n+m+1)*(n+m+2)))/4.d0 dpt2(kk)=anm*pnm(k1)+bnm*pnm(kk)+cnm*pnm(k2) enddo enddo 904 return end ! !*************************************************************************** ! subroutine BelPnm(pnm,maxn,t) !计算规格化Pnm !规格化Pnm采用改进的Belikov递推算法 implicit none integer::maxn,n,m,k,kk,k1,k2,k3,astat(6) real*8::pnm((maxn+2)**2) real*8::t,u,a,b,c,d,e real*8::anm,bnm,cnm,da !--------------------------------------------------------------------------- pnm=0.d0;u=dsqrt(1-t**2);pnm(1)=1.d0 pnm(2)=dsqrt(3.d0)*t;pnm(3)=dsqrt(3.d0)*u do n=1,maxn a=dsqrt(2.d0*n+1.d0)/dsqrt(2.d0*n-1.d0) b=dsqrt(2.d0*(n-1.d0)*(2.d0*n+1.d0))/dsqrt((2.d0*n-1.d0)*n) kk=n*(n+1)/2+1;k1=n*(n-1)/2+1;k2=n*(n-1)/2+2 pnm(kk)=a*t*pnm(k1)-b*u/2.d0*pnm(k2) do m=1,n kk=n*(n+1)/2+m+1;k1=n*(n-1)/2+m+1 k2=n*(n-1)/2+m+2;k3=n*(n-1)/2+m c=a/n*dsqrt(dble(n**2-m**2));d=a/n/2.d0*dsqrt(dble((n-m)*(n-m-1))) e=a/n/2.d0*dsqrt(dble((n+m)*(n+m-1))) if(m==1)e=e*dsqrt(2.d0) pnm(kk)=c*t*pnm(k1)-d*u*pnm(k2)+e*u*pnm(k3) enddo enddo 904 return end ! !*************************************************************************** ! subroutine Pndpn_dt(p,dp,n,t) implicit none integer::n,m,k real*8 ::p(n+5),dp(n+5),t !--------------------------------------------------------------------------- ! 计算全部Legendre函数pn及其导数dpn.p(1)=p0,dp(1)=dp0 p(1)=1.d0;p(2)=t;p(3)=1.5d0*t**2-0.5d0;p(4)=2.5d0*t**3-1.5d0*t dp(1)=0.d0;dp(2)=1.d0;dp(3)=3.d0*t;dp(4)=7.5d0*t**2-1.5d0 do m=3,n+1 k=m-1 p(m)=real(2*k-1)/real(k)*t*p(m-1)-real(k-1)/real(k)*p(m-2) dp(m)=t*dp(k)+real(m)*p(k) enddo end ! !****************************************************************************** ! subroutine LegPn_dt2(pn,dp1,dp2,n,t) !计算Pn(t)及其对θ一、二阶导数 implicit none integer::n,k real*8::pn(n),dp1(n),dp2(n),t,u !--------------------------------------------------------------------------- u=dsqrt(1-t**2) pn(1)=t;pn(2)=1.5d0*t**2-0.5d0 dp1(1)=-u;dp1(2)=-3.d0*t*u dp2(1)=-t;dp2(2)=3.d0*(1.d0-2.d0*t**2) do k=3,n pn(k)=dble(2*k-1)/dble(k)*t*pn(k-1)-dble(k-1)*pn(k-2)/dble(k) dp1(k)=dble(2*k-1)*(t*dp1(k-1)-u*pn(k-1))-dble(k-1)*dp1(k-2) dp1(k)=dp1(k)/dble(k) dp2(k)=(t*dp2(k-1)-2.d0*u*dp1(k-1)-t*pn(k-1))*dble(2*k-1) dp2(k)=(dp2(k)-dble(k-1)*dp2(k-2))/dble(k) enddo end