连带勒让德多项式

!//当theta<0.1°时,Pn2(cos(theta)) 计算不准!
!//当theta=90°时,Pn, Pn2(cos(theta)) 计算不准!
!// 一般情况可以使用.
!//当theta<0.1°时,Pn2(cos(theta)) 计算不准!
!//当theta=90°时,Pn, Pn2(cos(theta)) 计算不准!
!// 一般情况可以使用.
module LegendreP012
  !// 计算连带勒让德(不含"Condon-Shortley phase")
  implicit none
  contains
    !// 函数LgendrePn1,LgendrePn2需要调用LegndrePn.
    function LegendrePn(theta,degree) result(Pn)
      !// n阶勒让德函数
      !// theat in Unit of Degree, that is degree measure
      !// degree is the degree of Legendre Functions
      implicit none
      real   (kind=8),intent(in) :: theta
      integer(kind=8),intent(in) :: degree
      real   (kind=8) :: Pn
 
      !// temp variable
      integer(kind=8) :: l                  !// for degree loop
      real(kind=8) :: pfirst,psecond,pthird !// p_(k-1),p_k,p_(k+1)
      real(kind=8) :: ct,pi                 !// theta in unit of radian,pi=3.14
      real(kind=8) :: x                     !// cos(ct)
      pi=3.14159265358979d0
      ct=theta/180.0d0*pi
      x=cos(ct)
      l=0
      select case (degree)
      case(0)
        Pn=1.0d0
      case(1)
        Pn=x
      case(2:)
        pfirst=1.0d0
        psecond=x
        !// 递推计算 A,B->C
        l=2
        do
          pthird=(2-1.0d0/l)*x*psecond-(1-1.0d0/l)*pfirst
          if (l==degree) then
            Pn=pthird
            exit
          end if
          pfirst=psecond
          psecond=pthird
          l=l+1
        end do
        case default
          write(*,*) 'Error in "LegendrePn", Degree of Legendre Function should not be less than 0!'
      end select
 
    end function LegendrePn
 
    function LegendrePn1(theta,degree) result(Pn1)
      !// n阶1次勒让德函数
      !// Pn1, the degree is n and the order is 1
      !// theat in Unit of Degree, that is degree measure
      !// degree is the degree of Legendre Functions
      implicit none
      real   (kind=8),intent(in) :: theta
      integer(kind=8),intent(in) :: degree
      real   (kind=8) :: Pn1
 
      !// temp variable
      real(kind=8) :: pfirst,psecond
      real(kind=8) :: ct,pi                 !// theta in unit of radian,pi=3.14
      real(kind=8) :: x,y                   !// cos(ct), sin(ct)
      pi=3.14159265358979d0
      ct=theta/180.0d0*pi
      x=cos(ct)
      y=sin(ct)
 
      select case (degree)
      case(0)
        Pn1=0.0d0
      case(1:)
        pfirst=LegendrePn(theta,degree-1)
        psecond=LegendrePn(theta,degree)
        Pn1=degree/y*(pfirst-x*psecond)
      case default
          write(*,*) 'Error in "LegendrePn1", Degree of Legendre Function should not be less than 0!'
      end select
 
    end function LegendrePn1
 
    function LegendrePn2(theta,degree) result(Pn2)
      !// n阶2次勒让德函数
      !// Pn2, the degree is n and the order is 1
      !// theat in Unit of Degree, that is degree measure
      !// degree is the degree of Legendre Functions
 
      !// When theta is very small, LegendrePn2 will give wrong result!
 
      implicit none
      real   (kind=8),intent(in) :: theta
      integer(kind=8),intent(in) :: degree
      real   (kind=8) :: Pn2
 
      !// temp variable
      real(kind=8) :: pfirst,psecond
      real(kind=8) :: ct,pi                 !// theta in unit of radian,pi=3.14
      real(kind=8) :: x,y,x2,factor1,factor2   !// cos(ct), sin(ct), cos(2*ct)
 
      pi=3.14159265358979d0
      ct=theta/180.0d0*pi
      x=cos(ct)
      y=sin(ct)
      x2=cos(2.0d0*ct)
 
      select case (degree)
      case(0,1)
        Pn2=0.0d0
      case(2:)
          !// use recursion formula to calculate
          pfirst=LegendrePn(theta,degree-1)
          psecond=LegendrePn(theta,degree)
 
          factor1=2*x
          factor2=-((degree+1)*y*y+2*x*x)
          Pn2=degree/(y*y)*(factor1*pfirst+factor2*psecond)
          !// 或者
    !     Pn2=2*degree/(1-x2)*(factor1*pfirst+factor2*psecond)
    !
    !     !// 或者
    !     factor1=2*x/(y*y)
    !     factor2=1.0d0*(1-degree)-2.0d0/(y*y)
    !     Pn2=degree*factor1*pfirst+degree*factor2*psecond
 
      case default
          write(*,*) 'Error in "LegendrePn2", Degree of Legendre Function should not be less than 0!'
      end select
 
    end function LegendrePn2
 
 
    !// 子例行程序CLegendrePn,CLegendrePn1,CLegendrePn2是独立的,可以单独使用.
    subroutine CLegendrePn(choose,ntheta,thetalist,degree,LegendrePn,LegendrePn_sine,LegendrePnd,LegendrePndd)
 
      !// 使用时: 需要3步1),2),3)
      !// 1) call CLegendrePn ('prepare',   2_8,list_angle(1:2),loop_degree,LegendrePnd=pnd)
      !// 2) call CLegendrePn ('keep',      2_8,list_angle(1:2),loop_degree,LegendrePnd=pnd) loop_degree=0/1/2 开始连续调用即可.
      !// 3) call CLegendrePn ('deallocate',2_8,list_angle(1:2),loop_degree,LegendrePnd=pnd)
      !//
      implicit none
      character(len=*),intent(in) :: choose            !// prepare,keep,deallocate
      integer (kind=8),intent(in) :: ntheta            !// length of theta list
      real    (kind=8),intent(in) :: thetalist(ntheta) !// theats in unit of degree, that is degree measure
      integer (kind=8),intent(in) :: degree            !// Legendre的阶次
 
      real   (kind=8),intent(out),optional :: LegendrePn(ntheta),LegendrePn_sine(ntheta),LegendrePnd(ntheta),LegendrePndd(ntheta)
      !// 勒让德函数以及导数: Pn, Pn/sine(theta), dPn(cos(theta))/d(theta), d^2Pn(theta)/d(theta^2) 二阶导数
 
      !// temp variable
      integer(kind=8) :: theta_loop
      real   (kind=8) :: theta
      real   (kind=8) :: Pn(ntheta),Pn_sine(ntheta),Pnd(ntheta),Pndd(ntheta)
      !// 勒让德函数以及导数: Pn, Pn/sine(theta), dPn(cos(theta))/d(theta)
      integer(kind=8) :: l                                   !// for degree loop
      real   (kind=8),allocatable :: pfirst(:),psecond(:)    !// p_(k-1),p_k
      integer(kind=8) :: allstatus1,allstatus2
      real   (kind=8) :: pthird(ntheta)                      !// p_(k+1)
      real   (kind=8) :: ct,pi                               !// theta in unit of radian,pi=3.14
      real   (kind=8) :: xct,yct                             !// cos(ct),sin(ct)
      real   (kind=8) :: xctlist(ntheta),yctlist(ntheta)     !// all cos(ct), sin(ct)
      save :: pfirst,psecond
 
 
      !// 常数初值
      pi=3.14159265358979d0
 
      if (trim(choose)=='prepare' .or. trim(choose)=='PREPARE') then
        !// Legendre递推初值
        if (.not.allocated(pfirst)) allocate(pfirst(ntheta))
        if (.not.allocated(psecond)) allocate(psecond(ntheta))
        pfirst=1.0d0
        do theta_loop=1,ntheta
          theta=thetalist(theta_loop)
          ct=theta/180.0d0*pi
          xct=cos(ct)
          psecond(theta_loop)=xct
        end do
        Pn=0.0d0
        Pn_sine=0.0d0
        Pnd=0.0d0
        Pndd=0.0d0
      else if (trim(choose)=='keep' .or. trim(choose)=='KEEP') then
        !// get xct and yct
        do theta_loop=1,ntheta
          theta=thetalist(theta_loop)
          ct=theta/180.0d0*pi
          xctlist(theta_loop)=cos(ct)
          yctlist(theta_loop)=sin(ct)
        end do
        l=0
        !// get result for different degree
        select case (degree)
        case(0)
          Pn=1.0d0
          Pnd=0.0d0
          Pndd=0.0d0
          do theta_loop=1,ntheta
            yct=yctlist(theta_loop)
            Pn_sine(theta_loop)=1.0d0/yct
          end do
        case(1)
          do theta_loop=1,ntheta
            xct=xctlist(theta_loop)
            yct=yctlist(theta_loop)
            Pn(theta_loop)=xct
            Pn_sine(theta_loop)=xct/yct
            Pnd(theta_loop)=-yct
            Pndd(theta_loop)=-xct
          end do
        case(2:)
          !// 递推计算 A,B->C
          l=degree
          do theta_loop=1,ntheta
            xct=xctlist(theta_loop)
            yct=yctlist(theta_loop)
            pthird(theta_loop)=(2-1.0d0/l)*xct*psecond(theta_loop)-(1-1.0d0/l)*pfirst(theta_loop)
            Pn(theta_loop)=pthird(theta_loop)
            Pn_sine(theta_loop)=pthird(theta_loop)/yct
            Pnd(theta_loop)=l/yct*(xct*pthird(theta_loop)-psecond(theta_loop))
            Pndd(theta_loop)=l/(yct**2)*(-(1+l*yct**2)*pthird(theta_loop)+xct*psecond(theta_loop))
 
            pfirst(theta_loop)=psecond(theta_loop)
            psecond(theta_loop)=pthird(theta_loop)
          end do
        case default
          write(*,*) 'Error in "CLegendrePn", Degree of Legendre Function should not be less than 0!'
        end select
      else if (trim(choose)=='deallocate' .or. trim(choose)=='DEALLOCATE') then
        deallocate(pfirst, stat=allstatus1)
        deallocate(psecond,stat=allstatus2)
        if (allstatus1/=0 .or. allstatus2/=0) then
          write(*,*) 'Error! When deallocate arrays in "CLegendrePn"!'
        end if
      else
        write(*,*) 'Error! The first argument of "CLegendrePn" must be "prepare", or "keep", or "deallocate"!'
      end if
 
      !// 返回值
      if(present(LegendrePn)) LegendrePn=Pn
      if(present(LegendrePn_sine)) LegendrePn_sine=Pn_sine
      if(present(LegendrePnd)) LegendrePnd=Pnd
      if(present(LegendrePndd)) LegendrePndd=Pndd
 
 
    end subroutine CLegendrePn
 
    subroutine CLegendrePn1(choose,ntheta,thetalist,degree,LegendrePn1,LegendrePn1_sine,LegendrePn1d,LegendrePn1dd)
      !Attation: we define Pn1 = -d(Pn)/d(theta), as it used in Dislocation's thory, and it DO NOT have Condon–Shortley phase, not Same as Mathematica!
      !// 使用方法和前面一样.
      !// 再递推,可以从0/1/2开始向后连续(Continue)调用
      implicit none
      character(len=*),intent(in) :: choose             !// prepare,keep,deallocate
      integer (kind=8),intent(in) :: ntheta             !// length of theta list
      real    (kind=8),intent(in) :: thetalist(ntheta)  !// theat in unit of degree, that is degree measure
      integer (kind=8),intent(in) :: degree             !// Legendre的阶次
      real    (kind=8),intent(out),optional :: LegendrePn1(ntheta),LegendrePn1_sine(ntheta),&
                                               LegendrePn1d(ntheta),LegendrePn1dd(ntheta)
      !// 勒让德函数以及导数: Pn1, Pn1/sine(theta), dPn1(cos(theta))/d(theta), d^2Pn1(cos(theta))/d(theta^2)
      !// temp variable
      real   (kind=8) :: theta
      integer(kind=8) :: theta_loop
      real   (kind=8) :: Pn1(ntheta),Pn1_sine(ntheta),Pn1d(ntheta),Pn1dd(ntheta)
      integer(kind=8) :: l                               !// for degree loop
      real   (kind=8),allocatable:: pfirst(:),psecond(:) !// p_(k-1),p_k
      integer(kind=8) :: allstatus1,allstatus2
      real   (kind=8) :: pthird(ntheta)                  !// p_(k+1)
      real   (kind=8) :: ct,pi                           !// theta in unit of radian,pi=3.14
      real   (kind=8) :: xct,yct                         !// cos(ct),sin(ct)
      real   (kind=8) :: xctlist(ntheta),yctlist(ntheta) !// all cos(ct), sin(ct)
      save :: pfirst,psecond
      !// 常数初值
      pi=3.14159265358979d0
 
      if (trim(choose)=='prepare' .or. trim(choose)=='PREPARE') then
        !// Legendre递推初值
        if (.not.allocated(pfirst)) allocate(pfirst(ntheta))
        if (.not.allocated(psecond)) allocate(psecond(ntheta))
        pfirst=1.0d0
        do theta_loop=1,ntheta
          theta=thetalist(theta_loop)
          ct=theta/180.0d0*pi
          xct=cos(ct)
          psecond(theta_loop)=xct
        end do
        Pn1=0.0d0
        Pn1_sine=0.0d0
        Pn1d=0.0d0
        Pn1dd=0.0d0
      else if (trim(choose)=='keep' .or. trim(choose)=='KEEP') then
        !// get xct and yct
        do theta_loop=1,ntheta
          theta=thetalist(theta_loop)
          ct=theta/180.0d0*pi
          xctlist(theta_loop)=cos(ct)
          yctlist(theta_loop)=sin(ct)
        end do
        l=0
        !// get result for different degree
        select case (degree)
        case(0)
          Pn1=0.0d0
          Pn1_sine=0.0d0
          Pn1d=0.0d0
          Pn1dd=0.0d0
        case(1)
          Pn1_sine=1.0d0
          do theta_loop=1,ntheta
            xct=xctlist(theta_loop)
            yct=yctlist(theta_loop)
            Pn1(theta_loop)=(1-xct**2)/yct
            Pn1d(theta_loop)=xct
            Pn1dd(theta_loop)=-yct
          end do
        case(2:)
          !// 递推计算 A,B->C
          l=degree
          do theta_loop=1,ntheta
            xct=xctlist(theta_loop)
            yct=yctlist(theta_loop)
            pthird(theta_loop)=(2-1.0d0/l)*xct*psecond(theta_loop)-(1-1.0d0/l)*pfirst(theta_loop)
            Pn1(theta_loop)=-l/yct*(-psecond(theta_loop)+xct*pthird(theta_loop))
            Pn1_sine(theta_loop)=Pn1(theta_loop)/yct
            Pn1d(theta_loop)=l/(yct**2)*(-xct*psecond(theta_loop)+(1+l*yct**2)*pthird(theta_loop))
            Pn1dd(theta_loop)=l/(yct**3)*( (1-l-l**2 + (1+l+l**2)*(xct**2) )*psecond(theta_loop) + &
                                          (-2+(l**2)*(yct**2))*xct*pthird(theta_loop) )
 
            pfirst(theta_loop)=psecond(theta_loop)
            psecond(theta_loop)=pthird(theta_loop)
 
          end do
        case default
          write(*,*) 'Error in "CLegendrePn1", Degree of Legendre Function should not be less than 0!'
        end select
      else if (trim(choose)=='deallocate' .or. trim(choose)=='DEALLOCATE') then
        deallocate(pfirst, stat=allstatus1)
        deallocate(psecond,stat=allstatus2)
        if (allstatus1/=0 .or. allstatus2/=0) then
          write(*,*) 'Error! When deallocate arrays in "CLegendrePn1"!'
        end if
      else
        write(*,*) 'Error! The first argument of "CLegendrePn1" must be "prepare", or "keep", or "deallocate"!'
      end if
 
      if(present(LegendrePn1)) LegendrePn1=Pn1
      if(present(LegendrePn1_sine)) LegendrePn1_sine=Pn1_sine
      if(present(LegendrePn1d)) LegendrePn1d=Pn1d
      if(present(LegendrePn1dd)) LegendrePn1dd=Pn1dd
    end subroutine CLegendrePn1
 
    subroutine CLegendrePn2(choose,ntheta,thetalist,degree,LegendrePn2,LegendrePn2_sine,LegendrePn2d,LegendrePn2dd)
 
      !// 使用时,必须先调用一次CLegendrePn2(theta,0_8,.true.,Pn2)!
      !// 再递推,可以从0/1/2开始向后连续(Continue)调用
      implicit none
      character(len=*),intent(in) :: choose             !// prepare,keep,deallocate
      integer (kind=8),intent(in) :: ntheta             !// length of theta list
      real    (kind=8),intent(in) :: thetalist(ntheta)  !// theat in unit of degree, that is degree measure
      integer (kind=8),intent(in) :: degree             !// Legendre的阶次
      real    (kind=8),intent(out),optional :: LegendrePn2(ntheta),LegendrePn2_sine(ntheta),&
                                               LegendrePn2d(ntheta),LegendrePn2dd(ntheta)
      !// 勒让德函数以及导数: Pn2, Pn2/sine(theta), dPn2(cos(theta))/d(theta)
      !// temp variable
      real   (kind=8) :: theta
      integer(kind=8) :: theta_loop
      real   (kind=8) :: Pn2(ntheta),Pn2_sine(ntheta),Pn2d(ntheta),Pn2dd(ntheta)
      integer(kind=8) :: l                               !// for degree loop
      real   (kind=8),allocatable:: pfirst(:),psecond(:) !// p_(k-1),p_k
      integer(kind=8) :: allstatus1,allstatus2
      real   (kind=8) :: pthird(ntheta)                  !// p_(k+1)
      real   (kind=8) :: ct,pi                           !// theta in unit of radian,pi=3.14
      real   (kind=8) :: xct,yct                         !// cos(ct),sin(ct)
      real   (kind=8) :: xctlist(ntheta),yctlist(ntheta) !// all cos(ct), sin(ct)
      real   (kind=8) :: factor1,factor2
      save :: pfirst,psecond
      !// 常数初值
      pi=3.14159265358979d0
 
      if (trim(choose)=='prepare' .or. trim(choose)=='PREPARE') then
        !// Legendre递推初值
        if (.not.allocated(pfirst)) allocate(pfirst(ntheta))
        if (.not.allocated(psecond)) allocate(psecond(ntheta))
        pfirst=1.0d0
        do theta_loop=1,ntheta
          theta=thetalist(theta_loop)
          ct=theta/180.0d0*pi
          xct=cos(ct)
          psecond(theta_loop)=xct
        end do
        Pn2=0.0d0
        Pn2_sine=0.0d0
        Pn2d=0.0d0
        Pn2dd=0.0d0
      else if (trim(choose)=='keep' .or. trim(choose)=='KEEP') then
        !// get xct and yct
        do theta_loop=1,ntheta
          theta=thetalist(theta_loop)
          ct=theta/180.0d0*pi
          xctlist(theta_loop)=cos(ct)
          yctlist(theta_loop)=sin(ct)
        end do
        l=0
        !// get result for different degree
        select case (degree)
        case(0,1)
          Pn2=0.0d0
          Pn2_sine=0.0d0
          Pn2d=0.0d0
          Pn2dd=0.0d0
        case(2:)
          !// 递推计算 A,B->C
          l=degree
          do theta_loop=1,ntheta
            xct=xctlist(theta_loop)
            yct=yctlist(theta_loop)
            factor1=2*xct
            factor2=-((l+1)*yct**2+2*xct**2)
            pthird(theta_loop)=(2-1.0d0/l)*xct*psecond(theta_loop)-(1-1.0d0/l)*pfirst(theta_loop)
            Pn2(theta_loop)=l/(yct**2)*(factor1*psecond(theta_loop)+factor2*pthird(theta_loop))
            Pn2_sine(theta_loop)=Pn2(theta_loop)/yct
 
            factor1=(l**2+l+2)*yct**2-4
            factor2=((l-l**2)*yct**2+4)*xct
            Pn2d(theta_loop)=l/(yct**3)*(factor1*psecond(theta_loop)+factor2*pthird(theta_loop))
 
            factor1=2*(1+3*l-l**2) + (8-9*l-l**2)*(yct**2) + 2*(l**2-3*l+5)*(xct**2)
            factor2=-4*(l+1) + (l*(l+1)**2)*(yct**2) + 4*(l-2)*xct**2 - (l**2)*(l-1)*(xct**2)*(yct**2)
            Pn2dd(theta_loop)=l/(yct**4)*(factor1*xct*psecond(theta_loop)+factor2*pthird(theta_loop))
 
            pfirst(theta_loop)=psecond(theta_loop)
            psecond(theta_loop)=pthird(theta_loop)
          end do
        case default
          write(*,*) 'Error in "CLegendrePn2", Degree of Legendre Function should not be less than 0!'
        end select
      else if (trim(choose)=='deallocate' .or. trim(choose)=='DEALLOCATE') then
        deallocate(pfirst, stat=allstatus1)
        deallocate(psecond,stat=allstatus2)
        if (allstatus1/=0 .or. allstatus2/=0) then
          write(*,*) 'Error! When deallocate arrays in "CLegendrePn2"!'
        end if
      end if
      if(present(LegendrePn2)) LegendrePn2=Pn2
      if(present(LegendrePn2_sine)) LegendrePn2_sine=Pn2_sine
      if(present(LegendrePn2d)) LegendrePn2d=Pn2d
      if(present(LegendrePn2dd)) LegendrePn2dd=Pn2dd
 
    end subroutine CLegendrePn2 
 
end module LegendreP012