!//当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
