1870 lines
60 KiB
Fortran
1870 lines
60 KiB
Fortran
|
|
|
|
#include "macrodef.fh"
|
|
|
|
subroutine get_initial_nbhs_null(ex,crho,sigma,x,RJ,IJ,omega,sst,Rmin)
|
|
|
|
implicit none
|
|
!~~~~~~> Input parameters:
|
|
|
|
integer,intent(in ):: ex(1:3),sst
|
|
real*8,intent(in ) :: Rmin
|
|
double precision,intent(in),dimension(ex(1))::crho
|
|
double precision,intent(in),dimension(ex(2))::sigma
|
|
double precision,intent(in),dimension(ex(3))::x
|
|
double precision,intent(inout),dimension(ex(1),ex(2),ex(3))::RJ,IJ,omega
|
|
|
|
!~~~~~~> Other variables:
|
|
real*8 :: xe
|
|
real*8,dimension(ex(1),ex(2)) :: RJe,IJe
|
|
integer :: k
|
|
|
|
xe = x(1)
|
|
RJe = RJ(:,:,1)
|
|
IJe = IJ(:,:,1)
|
|
|
|
do k=1,ex(3)
|
|
RJ(:,:,k) = RJe*(1.d0-x(k))*xe/(1-xe)/x(k)
|
|
IJ(:,:,k) = IJe*(1.d0-x(k))*xe/(1-xe)/x(k)
|
|
enddo
|
|
|
|
omega = 1.d0
|
|
|
|
return
|
|
|
|
end subroutine get_initial_nbhs_null
|
|
!-----------------------------------
|
|
!Eq.(10) of CQG 24, S327 (2007)
|
|
!----------------------------------
|
|
function Zslm(s,l,m,the,phi) result(gont)
|
|
implicit none
|
|
integer,intent(in) :: s,l,m
|
|
real*8,intent(in) :: the,phi
|
|
|
|
double complex :: Yslm,gont,II
|
|
|
|
II=dcmplx(0.d0,1.d0)
|
|
|
|
if(m>0)then
|
|
gont = Yslm(s,l,m,the,phi)
|
|
if(m/2*2==m)then
|
|
gont = gont+Yslm(s,l,-m,the,phi)
|
|
else
|
|
gont = gont-Yslm(s,l,-m,the,phi)
|
|
endif
|
|
gont = gont/dsqrt(2.d0)
|
|
elseif(m<0)then
|
|
gont = -Yslm(s,l,-m,the,phi)
|
|
if(m/2*2==m)then
|
|
gont = gont+Yslm(s,l,m,the,phi)
|
|
else
|
|
gont = gont-Yslm(s,l,m,the,phi)
|
|
endif
|
|
gont = II*gont/dsqrt(2.d0)
|
|
else
|
|
gont = Yslm(s,l,m,the,phi)
|
|
endif
|
|
|
|
return
|
|
|
|
end function Zslm
|
|
|
|
!#define SCH
|
|
|
|
#ifdef SCH
|
|
|
|
subroutine get_initial_null(ex,crho,sigma,R,RJ,IJ,sst,Rmin)
|
|
implicit none
|
|
! argument variables
|
|
integer, intent(in ):: ex(1:3),sst
|
|
real*8,intent(in) :: Rmin
|
|
double precision,intent(in),dimension(ex(1))::crho
|
|
double precision,intent(in),dimension(ex(2))::sigma
|
|
double precision,intent(in),dimension(ex(3))::R
|
|
real*8,dimension(ex(1),ex(2),ex(3)),intent(out)::RJ,IJ
|
|
|
|
RJ = 0.d0
|
|
IJ = 0.d0
|
|
|
|
return
|
|
|
|
end subroutine get_initial_null
|
|
!-------------------------
|
|
subroutine get_null_boundary(ex,crho,sigma,R,beta,RQ,IQ,RU,IU,W,RTheta,ITheta, &
|
|
quR1,quR2,quI1,quI2,qlR1,qlR2,qlI1,qlI2, &
|
|
gR,gI, &
|
|
dquR1,dquR2,dquI1,dquI2, &
|
|
bdquR1,bdquR2,bdquI1,bdquI2, &
|
|
dgR,dgI,bdgR,bdgI,T,Rmin,sst)
|
|
|
|
implicit none
|
|
|
|
!~~~~~~% Input parameters:
|
|
integer,intent(in) :: ex(3),sst
|
|
real*8,intent(in) :: T,Rmin
|
|
double precision,intent(in),dimension(ex(1))::crho
|
|
double precision,intent(in),dimension(ex(2))::sigma
|
|
double precision,intent(in),dimension(ex(3))::R
|
|
real*8,dimension(ex(1),ex(2),ex(3)),intent(inout) :: beta,RQ,IQ,RU,IU
|
|
real*8,dimension(ex(1),ex(2),ex(3)),intent(inout) :: W,RTheta,ITheta
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: quR1,quR2,quI1,quI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: qlR1,qlR2,qlI1,qlI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: gR,gI
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: dquR1,dquR2,dquI1,dquI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: bdquR1,bdquR2,bdquI1,bdquI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: dgR,dgI,bdgR,bdgI
|
|
|
|
integer :: k
|
|
k=1
|
|
|
|
beta(:,:,k) = 0.d0
|
|
W(:,:,k) =-2.d0/R(k)**2/Rmin**2*(1.d0-R(k))**2
|
|
|
|
RQ(:,:,k) = 0.d0
|
|
IQ(:,:,k) = 0.d0
|
|
|
|
RTheta(:,:,k) = 0.d0
|
|
ITheta(:,:,k) = 0.d0
|
|
|
|
RU(:,:,k) = 0.d0
|
|
IU(:,:,k) = 0.d0
|
|
|
|
return
|
|
|
|
end subroutine get_null_boundary
|
|
!-------------------------------------------------------------
|
|
subroutine get_exact_null(ex,crho,sigma,R,RJ,IJ,sst,Rmin,T, &
|
|
quR1,quR2,quI1,quI2,qlR1,qlR2,qlI1,qlI2, &
|
|
gR,gI, &
|
|
dquR1,dquR2,dquI1,dquI2, &
|
|
bdquR1,bdquR2,bdquI1,bdquI2, &
|
|
dgR,dgI,bdgR,bdgI)
|
|
implicit none
|
|
! argument variables
|
|
integer, intent(in ):: ex(1:3),sst
|
|
real*8,intent(in) :: Rmin,T
|
|
double precision,intent(in),dimension(ex(1))::crho
|
|
double precision,intent(in),dimension(ex(2))::sigma
|
|
double precision,intent(in),dimension(ex(3))::R
|
|
real*8,dimension(ex(1),ex(2),ex(3)),intent(out)::RJ,IJ
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: quR1,quR2,quI1,quI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: qlR1,qlR2,qlI1,qlI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: gR,gI
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: dquR1,dquR2,dquI1,dquI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: bdquR1,bdquR2,bdquI1,bdquI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: dgR,dgI,bdgR,bdgI
|
|
|
|
RJ = 0.d0
|
|
IJ = 0.d0
|
|
|
|
return
|
|
|
|
end subroutine get_exact_null
|
|
!-------------------------------------------------------------------------------------------
|
|
subroutine get_null_boundary_c(ex,crho,sigma,R,beta,RQ,IQ,RU,IU,W,RTheta,ITheta, &
|
|
quR1,quR2,quI1,quI2,qlR1,qlR2,qlI1,qlI2, &
|
|
gR,gI, &
|
|
dquR1,dquR2,dquI1,dquI2, &
|
|
bdquR1,bdquR2,bdquI1,bdquI2, &
|
|
dgR,dgI,bdgR,bdgI,T,Rmin,sst)
|
|
|
|
implicit none
|
|
|
|
!~~~~~~% Input parameters:
|
|
integer,intent(in) :: ex(3),sst
|
|
real*8,intent(in) :: T,Rmin
|
|
double precision,intent(in),dimension(ex(1))::crho
|
|
double precision,intent(in),dimension(ex(2))::sigma
|
|
double precision,intent(in),dimension(ex(3))::R
|
|
real*8,dimension(ex(1),ex(2),ex(3)),intent(inout) :: beta,RQ,IQ,RU,IU
|
|
real*8,dimension(ex(1),ex(2),ex(3)),intent(inout) :: W,RTheta,ITheta
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: quR1,quR2,quI1,quI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: qlR1,qlR2,qlI1,qlI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: gR,gI
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: dquR1,dquR2,dquI1,dquI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: bdquR1,bdquR2,bdquI1,bdquI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: dgR,dgI,bdgR,bdgI
|
|
|
|
integer :: k
|
|
|
|
do k=1,ex(3)
|
|
|
|
beta(:,:,k) = 0.d0
|
|
W(:,:,k) =-2.d0/R(k)**2/Rmin**2*(1.d0-R(k))**2
|
|
|
|
RQ(:,:,k) = 0.d0
|
|
IQ(:,:,k) = 0.d0
|
|
|
|
RTheta(:,:,k) = 0.d0
|
|
ITheta(:,:,k) = 0.d0
|
|
|
|
RU(:,:,k) = 0.d0
|
|
IU(:,:,k) = 0.d0
|
|
enddo
|
|
|
|
return
|
|
|
|
end subroutine get_null_boundary_c
|
|
|
|
#else
|
|
|
|
#if 0
|
|
! for some trival check
|
|
#if 1
|
|
!-------------------------------------------------------------
|
|
! Linear wave given in CQG 24S327
|
|
!-------------------------------------------------------------
|
|
subroutine get_initial_null(ex,crho,sigma,R,RJ,IJ,sst,Rmin)
|
|
implicit none
|
|
! argument variables
|
|
integer, intent(in ):: ex(1:3),sst
|
|
real*8,intent(in) :: Rmin
|
|
double precision,intent(in),dimension(ex(1))::crho
|
|
double precision,intent(in),dimension(ex(2))::sigma
|
|
double precision,intent(in),dimension(ex(3))::R
|
|
real*8,dimension(ex(1),ex(2),ex(3)),intent(out)::RJ,IJ
|
|
|
|
integer :: i,j,k
|
|
real*8 ::x,y,z,gr,gt,gp,tgrho,tgsigma
|
|
double complex :: Yslm,II,Jr
|
|
|
|
double complex :: beta0,C1,C2
|
|
integer :: nu,m
|
|
|
|
call initial_null_paramter(beta0,C1,C2,nu,m)
|
|
|
|
II = dcmplx(0.d0,1.d0)
|
|
|
|
do i=1,ex(1)
|
|
do j=1,ex(2)
|
|
do k=1,ex(3)
|
|
! here fake global coordinate is enough
|
|
gr = 1.d0
|
|
tgrho = dtan(crho(i))
|
|
tgsigma = dtan(sigma(j))
|
|
select case (sst)
|
|
case (0)
|
|
z = gr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = z*tgrho
|
|
y = z*tgsigma
|
|
case (1)
|
|
z = -gr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = z*tgrho
|
|
y = z*tgsigma
|
|
case (2)
|
|
x = gr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
y = x*tgrho
|
|
z = x*tgsigma
|
|
case (3)
|
|
x = -gr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
y = x*tgrho
|
|
z = x*tgsigma
|
|
case (4)
|
|
y = gr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = y*tgrho
|
|
z = y*tgsigma
|
|
case (5)
|
|
y = -gr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = y*tgrho
|
|
z = y*tgsigma
|
|
case default
|
|
write(*,*) "get_initial_null: not recognized sst = ",sst
|
|
return
|
|
end select
|
|
gt = dacos(z/gr)
|
|
gp = datan2(y,x)
|
|
|
|
gr = (1.d0-R(k))/R(k)/Rmin
|
|
Jr = (2.4d1*beta0+3.d0*II*nu*C1-II*nu**3*C2)/3.6d1+C1/4.d0*gr-C2/1.2d1*gr**3
|
|
gr = dreal(Jr)
|
|
Jr = Yslm(0,2,m,gt,gp)
|
|
RJ(i,j,k) = gr*dreal(Jr)
|
|
IJ(i,j,k) = gr*dimag(Jr)
|
|
|
|
#if 0
|
|
RJ(i,j,k) = 0.25d0*dsqrt(5.d0/3.1415926)*(3/(1.d0+tgrho*tgrho+tgsigma*tgsigma)-1.d0)
|
|
IJ(i,j,k) = 0.d0
|
|
#endif
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
return
|
|
|
|
end subroutine get_initial_null
|
|
#else
|
|
! for check usage
|
|
subroutine get_initial_null(ex,crho,sigma,R,RJ,IJ,sst,Rmin)
|
|
implicit none
|
|
! argument variables
|
|
integer, intent(in ):: ex(1:3),sst
|
|
real*8,intent(in) :: Rmin
|
|
double precision,intent(in),dimension(ex(1))::crho
|
|
double precision,intent(in),dimension(ex(2))::sigma
|
|
double precision,intent(in),dimension(ex(3))::R
|
|
real*8,dimension(ex(1),ex(2),ex(3)),intent(out)::RJ,IJ
|
|
|
|
real*8 :: thetac,thetas,sr,ss,cr,cs,srss,crcs,tcts,tcts2
|
|
real*8 :: sr2,ss2,cr2,cs2,tc2,ts2
|
|
integer :: i,j,k
|
|
real*8 :: ggr,tgrho,tgsigma
|
|
|
|
real*8 :: PI
|
|
|
|
PI = dacos(-1.d0)
|
|
|
|
do i=1,ex(1)
|
|
do j=1,ex(2)
|
|
do k=1,ex(3)
|
|
sr = dsin(crho(i))
|
|
ss = dsin(sigma(j))
|
|
cr = dcos(crho(i))
|
|
cs = dcos(sigma(j))
|
|
srss = sr*ss
|
|
crcs = cr*cs
|
|
sr2 = sr*sr
|
|
ss2 = ss*ss
|
|
cr2 = cr*cr
|
|
cs2 = cs*cs
|
|
thetac = dsqrt((1.d0-srss)/2.d0)
|
|
thetas = dsqrt((1.d0+srss)/2.d0)
|
|
tc2 = thetac*thetac
|
|
ts2 = thetas*thetas
|
|
tcts = thetac*thetas
|
|
tcts2 = tcts*tcts
|
|
! q^Aq^B@_A@_B Y20
|
|
RJ(i,j,k) =-1.5d0*dsqrt(5.d0/PI)*sr*ss*(4.d0*cr2*cs2+cs2+cr2)
|
|
IJ(i,j,k) = 3.d0*dsqrt(5.d0/PI)*thetac*thetas*(cs2-cr2)
|
|
|
|
! @_rho@_rho Y20
|
|
RJ(i,j,k) = 1.5d0*dsqrt(5.d0/PI)*cs2*cs2*(-cr2*cs2+2*cr2*cr2*cs2-2*cr2*cr2-cs2+3*cr2) &
|
|
/(3*cr2*cs2*cs2*cs2-3*cr2*cr2*cr2*cs2*cs2-3*cr2*cr2*cs2*cs2*cs2+ &
|
|
cr2*cr2*cr2*cs2*cs2*cs2-3*cr2*cr2*cs2-cs2*cs2*cs2-cr2*cr2*cr2+ &
|
|
3*cs2*cr2*cr2*cr2+6*cr2*cr2*cs2*cs2-3*cs2*cs2*cr2)
|
|
IJ(i,j,k) = 0.d0
|
|
! q^Aq^B h_AB
|
|
RJ(i,j,k) = 0.d0
|
|
IJ(i,j,k) = 0.d0
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
return
|
|
|
|
end subroutine get_initial_null
|
|
#endif
|
|
#endif
|
|
!======================================================================================
|
|
!-------------------------------------------------------------
|
|
! Linear wave given in CQG 24S327
|
|
!-------------------------------------------------------------
|
|
subroutine get_initial_null(ex,crho,sigma,R,RJ,IJ,sst,Rmin)
|
|
implicit none
|
|
! argument variables
|
|
integer, intent(in ):: ex(1:3),sst
|
|
real*8,intent(in) :: Rmin
|
|
double precision,intent(in),dimension(ex(1))::crho
|
|
double precision,intent(in),dimension(ex(2))::sigma
|
|
double precision,intent(in),dimension(ex(3))::R
|
|
real*8,dimension(ex(1),ex(2),ex(3)),intent(out)::RJ,IJ
|
|
|
|
integer :: i,j,k
|
|
real*8 ::x,y,z,gr,gt,gp,tgrho,tgsigma,tc,ts
|
|
double complex :: Yslm,II,Jr
|
|
|
|
double complex :: beta0,C1,C2
|
|
integer :: nu,m
|
|
|
|
double complex :: swtf,ff
|
|
|
|
call initial_null_paramter(beta0,C1,C2,nu,m)
|
|
|
|
II = dcmplx(0.d0,1.d0)
|
|
|
|
do i=1,ex(1)
|
|
do j=1,ex(2)
|
|
do k=1,ex(3)
|
|
!fake global coordinate is enough here
|
|
gr = 1.d0
|
|
tgrho = dtan(crho(i))
|
|
tgsigma = dtan(sigma(j))
|
|
tc = dsqrt((1.d0-dsin(crho(i))*dsin(sigma(j)))/2.d0)
|
|
ts = dsqrt((1.d0+dsin(crho(i))*dsin(sigma(j)))/2.d0)
|
|
select case (sst)
|
|
case (0)
|
|
z = gr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = z*tgrho
|
|
y = z*tgsigma
|
|
case (1)
|
|
z = -gr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = z*tgrho
|
|
y = z*tgsigma
|
|
case (2)
|
|
x = gr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
y = x*tgrho
|
|
z = x*tgsigma
|
|
case (3)
|
|
x = -gr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
y = x*tgrho
|
|
z = x*tgsigma
|
|
case (4)
|
|
y = gr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = y*tgrho
|
|
z = y*tgsigma
|
|
case (5)
|
|
y = -gr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = y*tgrho
|
|
z = y*tgsigma
|
|
case default
|
|
write(*,*) "get_initial_null: not recognized sst = ",sst
|
|
return
|
|
end select
|
|
gt = dacos(z/gr)
|
|
gp = datan2(y,x)
|
|
swtf = 2.d0*tc*ts*(ts+II*tc)/dcos(sigma(j))
|
|
if(sst==1 .or. sst==3 .or. sst==4) swtf = dconjg(swtf)
|
|
select case (sst)
|
|
case (0,1)
|
|
swtf = swtf/(dcos(gp)+II*dcos(gt)*dsin(gp))*(dcos(gt)**2+dsin(gt)**2*dcos(gp)**2)
|
|
case (2,3)
|
|
swtf = II*swtf*dsin(gt)
|
|
case (4,5)
|
|
swtf = -II*swtf*dsin(gt)
|
|
end select
|
|
|
|
gr = (1.d0-R(k))/R(k)/Rmin
|
|
Jr = (2.4d1*beta0+3.d0*II*nu*C1-II*nu**3*C2)/3.6d1+C1/4.d0*gr-C2/1.2d1*gr**3
|
|
gr = dreal(Jr)
|
|
Jr = Yslm(2,2,m,gt,gp)
|
|
ff = dsqrt(dble((2-1)*2*(2+1)*(2+2)))*gr*Jr*swtf**2
|
|
RJ(i,j,k) = dreal(ff)
|
|
IJ(i,j,k) = dimag(ff)
|
|
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
return
|
|
|
|
end subroutine get_initial_null
|
|
!==============================================================================================
|
|
|
|
#if 0
|
|
! for checking derivs_eth and dderivs_eth
|
|
!-------------------------
|
|
subroutine get_null_boundary(ex,crho,sigma,R,beta,RQ,IQ,RU,IU,W,RTheta,ITheta, &
|
|
quR1,quR2,quI1,quI2,qlR1,qlR2,qlI1,qlI2, &
|
|
gR,gI, &
|
|
dquR1,dquR2,dquI1,dquI2, &
|
|
bdquR1,bdquR2,bdquI1,bdquI2, &
|
|
dgR,dgI,bdgR,bdgI,T,Rmin,sst)
|
|
|
|
implicit none
|
|
|
|
!~~~~~~% Input parameters:
|
|
integer,intent(in) :: ex(3),sst
|
|
real*8,intent(in) :: T,Rmin
|
|
double precision,intent(in),dimension(ex(1))::crho
|
|
double precision,intent(in),dimension(ex(2))::sigma
|
|
double precision,intent(in),dimension(ex(3))::R
|
|
real*8,dimension(ex(1),ex(2),ex(3)),intent(inout) :: beta,RQ,IQ,RU,IU
|
|
real*8,dimension(ex(1),ex(2),ex(3)),intent(inout) :: W,RTheta,ITheta
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: quR1,quR2,quI1,quI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: qlR1,qlR2,qlI1,qlI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: gR,gI
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: dquR1,dquR2,dquI1,dquI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: bdquR1,bdquR2,bdquI1,bdquI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: dgR,dgI,bdgR,bdgI
|
|
|
|
double complex,dimension(ex(1),ex(2)) :: Y20,dY20,ddY20,f
|
|
integer :: i,j,k
|
|
real*8 ::x,y,z,hgr,gt,gp,tgrho,tgsigma
|
|
double complex :: Yslm,II,Jr
|
|
|
|
double complex :: beta0,C1,C2
|
|
integer :: nu,m
|
|
|
|
call initial_null_paramter(beta0,C1,C2,nu,m)
|
|
|
|
II = dcmplx(0.d0,1.d0)
|
|
|
|
k=1
|
|
do i=1,ex(1)
|
|
do j=1,ex(2)
|
|
! fake global coordinate is enough
|
|
hgr = 1.d0
|
|
tgrho = dtan(crho(i))
|
|
tgsigma = dtan(sigma(j))
|
|
select case (sst)
|
|
case (0)
|
|
z = hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = z*tgrho
|
|
y = z*tgsigma
|
|
case (1)
|
|
z = -hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = z*tgrho
|
|
y = z*tgsigma
|
|
case (2)
|
|
x = hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
y = x*tgrho
|
|
z = x*tgsigma
|
|
case (3)
|
|
x = -hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
y = x*tgrho
|
|
z = x*tgsigma
|
|
case (4)
|
|
y = hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = y*tgrho
|
|
z = y*tgsigma
|
|
case (5)
|
|
y = -hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = y*tgrho
|
|
z = y*tgsigma
|
|
case default
|
|
write(*,*) "get_initial_null: not recognized sst = ",sst
|
|
return
|
|
end select
|
|
gt = dacos(z/hgr)
|
|
gp = datan2(y,x)
|
|
|
|
hgr = (1.d0-R(k))/R(k)/Rmin
|
|
Jr = (2.4d1*beta0+3.d0*II*nu*C1-II*nu**3*C2)/3.6d1+C1/4.d0*hgr-C2/1.2d1*hgr**3
|
|
RTheta(i,j,k) = dreal(Jr*nu*(II*dcos(nu*T)-dsin(nu*T)))
|
|
f(i,j) = dreal(Jr*(dcos(nu*T)+II*dsin(nu*T)))
|
|
|
|
Jr = (-2.4d1*II*nu*beta0+3.d0*nu*nu*C1-nu**4*C2)/36.d0+2.d0*beta0*hgr&
|
|
+C1/2.d0*hgr*hgr+II*nu*C2/3.d0*hgr**3+C2/4.d0*hgr**4
|
|
RU(i,j,k) = dreal(Jr*(dcos(nu*T)+II*dsin(nu*T)))
|
|
|
|
! Re(r^2*Ul_,r*exp(i nu T)) of CQG 24S327, (12) indeed
|
|
Jr = -2.d0*beta0-C1*hgr-II*nu*C2*hgr**2-C2*hgr**3
|
|
RQ(i,j,k) = dreal(Jr*(dcos(nu*T)+II*dsin(nu*T)))
|
|
|
|
Jr = (2.4d1*II*nu*beta0-3.d0*nu*nu*C1+nu**4*C2)/6.d0+(3.d0*II*nu*C1-6.d0*beta0-II*nu**3*C2)/3.d0*hgr&
|
|
-nu*nu*C2*hgr*hgr+II*nu*C2*hgr**3+C2/2.d0*hgr**4
|
|
W(i,j,k) = dreal(Jr*(dcos(nu*T)+II*dsin(nu*T)))
|
|
|
|
Y20(i,j) = Yslm(0,2,m,gt,gp)
|
|
|
|
enddo
|
|
enddo
|
|
|
|
call derivs_eth(ex(1:2),crho,sigma,Y20,dY20,0,1, &
|
|
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
|
|
call dderivs_eth(ex(1:2),crho,sigma,Y20,ddY20,0,1,1, &
|
|
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k),&
|
|
dquR1(:,:,k),dquR2(:,:,k),dquI1(:,:,k),dquI2(:,:,k), &
|
|
bdquR1(:,:,k),bdquR2(:,:,k),bdquI1(:,:,k),bdquI2(:,:,k), &
|
|
dgR(:,:,k),dgI(:,:,k),bdgR(:,:,k),bdgI(:,:,k))
|
|
|
|
beta(:,:,k) = dreal(beta0*(dcos(nu*T)+II*dsin(nu*T)))*dreal(Y20)
|
|
W(:,:,k) = W(:,:,k)*dreal(Y20)
|
|
|
|
f = dexp(-2.d0*beta(:,:,k))*(f*ddY20*RQ(:,:,k)*dconjg(dY20)+dsqrt(1.d0+abs(f*ddY20))*RQ(:,:,k)*dY20)
|
|
RQ(:,:,k) = dreal(f)
|
|
IQ(:,:,k) = dimag(f)
|
|
|
|
f = ddY20*RTheta(:,:,k)
|
|
RTheta(:,:,k) = dreal(f)
|
|
ITheta(:,:,k) = dimag(f)
|
|
|
|
f = dY20*RU(:,:,k)
|
|
RU(:,:,k) = dreal(f)
|
|
IU(:,:,k) = dimag(f)
|
|
|
|
return
|
|
|
|
end subroutine get_null_boundary
|
|
#else
|
|
!-------------------------
|
|
subroutine get_null_boundary(ex,crho,sigma,R,beta,RQ,IQ,RU,IU,W,RTheta,ITheta, &
|
|
quR1,quR2,quI1,quI2,qlR1,qlR2,qlI1,qlI2, &
|
|
gR,gI, &
|
|
dquR1,dquR2,dquI1,dquI2, &
|
|
bdquR1,bdquR2,bdquI1,bdquI2, &
|
|
dgR,dgI,bdgR,bdgI,T,Rmin,sst)
|
|
|
|
implicit none
|
|
|
|
!~~~~~~% Input parameters:
|
|
integer,intent(in) :: ex(3),sst
|
|
real*8,intent(in) :: T,Rmin
|
|
double precision,intent(in),dimension(ex(1))::crho
|
|
double precision,intent(in),dimension(ex(2))::sigma
|
|
double precision,intent(in),dimension(ex(3))::R
|
|
real*8,dimension(ex(1),ex(2),ex(3)),intent(inout) :: beta,RQ,IQ,RU,IU
|
|
real*8,dimension(ex(1),ex(2),ex(3)),intent(inout) :: W,RTheta,ITheta
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: quR1,quR2,quI1,quI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: qlR1,qlR2,qlI1,qlI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: gR,gI
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: dquR1,dquR2,dquI1,dquI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: bdquR1,bdquR2,bdquI1,bdquI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: dgR,dgI,bdgR,bdgI
|
|
|
|
integer :: i,j,k
|
|
real*8 ::x,y,z,hgr,gt,gp,tgrho,tgsigma,tc,ts,rf
|
|
double complex :: Yslm,II,Jr,swtf,ff
|
|
|
|
double complex :: beta0,C1,C2
|
|
integer :: nu,m
|
|
|
|
call initial_null_paramter(beta0,C1,C2,nu,m)
|
|
|
|
II = dcmplx(0.d0,1.d0)
|
|
|
|
k=1
|
|
do i=1,ex(1)
|
|
do j=1,ex(2)
|
|
! fake global coordinate is enough here
|
|
hgr = 1.d0
|
|
tgrho = dtan(crho(i))
|
|
tgsigma = dtan(sigma(j))
|
|
tc = dsqrt((1.d0-dsin(crho(i))*dsin(sigma(j)))/2.d0)
|
|
ts = dsqrt((1.d0+dsin(crho(i))*dsin(sigma(j)))/2.d0)
|
|
select case (sst)
|
|
case (0)
|
|
z = hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = z*tgrho
|
|
y = z*tgsigma
|
|
case (1)
|
|
z = -hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = z*tgrho
|
|
y = z*tgsigma
|
|
case (2)
|
|
x = hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
y = x*tgrho
|
|
z = x*tgsigma
|
|
case (3)
|
|
x = -hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
y = x*tgrho
|
|
z = x*tgsigma
|
|
case (4)
|
|
y = hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = y*tgrho
|
|
z = y*tgsigma
|
|
case (5)
|
|
y = -hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = y*tgrho
|
|
z = y*tgsigma
|
|
case default
|
|
write(*,*) "get_null_boundary: not recognized sst = ",sst
|
|
return
|
|
end select
|
|
gt = dacos(z/hgr)
|
|
gp = datan2(y,x)
|
|
swtf = 2.d0*tc*ts*(ts+II*tc)/dcos(sigma(j))
|
|
if(sst==1 .or. sst==3 .or. sst==4) swtf = dconjg(swtf)
|
|
select case (sst)
|
|
case (0,1)
|
|
swtf = swtf/(dcos(gp)+II*dcos(gt)*dsin(gp))*(dcos(gt)**2+dsin(gt)**2*dcos(gp)**2)
|
|
case (2,3)
|
|
swtf = II*swtf*dsin(gt)
|
|
case (4,5)
|
|
swtf = -II*swtf*dsin(gt)
|
|
end select
|
|
|
|
hgr = (1.d0-R(k))/R(k)/Rmin
|
|
Jr = (2.4d1*beta0+3.d0*II*nu*C1-II*nu**3*C2)/3.6d1+C1/4.d0*hgr-C2/1.2d1*hgr**3
|
|
rf = dreal(Jr*nu*(II*dcos(nu*T)-dsin(nu*T)))
|
|
Jr = Yslm(2,2,m,gt,gp)
|
|
ff = dsqrt(dble((2-1)*2*(2+1)*(2+2)))*rf*Jr*swtf**2
|
|
RTheta(i,j,k) = dreal(ff)
|
|
ITheta(i,j,k) = dimag(ff)
|
|
|
|
rf = dreal(Yslm(0,2,m,gt,gp))
|
|
beta(i,j,k) = rf*dreal(beta0*(dcos(nu*T)+II*dsin(nu*T)))
|
|
Jr = (2.4d1*II*nu*beta0-3.d0*nu*nu*C1+nu**4*C2)/6.d0+(3.d0*II*nu*C1-6.d0*beta0-II*nu**3*C2)/3.d0*hgr&
|
|
-nu*nu*C2*hgr*hgr+II*nu*C2*hgr**3+C2/2.d0*hgr**4
|
|
W(i,j,k) = rf*dreal(Jr*(dcos(nu*T)+II*dsin(nu*T)))
|
|
|
|
Jr = (-2.4d1*II*nu*beta0+3.d0*nu*nu*C1-nu**4*C2)/36.d0+2.d0*beta0*hgr&
|
|
+C1/2.d0*hgr*hgr+II*nu*C2/3.d0*hgr**3+C2/4.d0*hgr**4
|
|
rf = dreal(Jr*(dcos(nu*T)+II*dsin(nu*T)))
|
|
Jr = Yslm(1,2,m,gt,gp)
|
|
ff = dsqrt(dble(2*(2+1)))*rf*Jr*swtf
|
|
RU(i,j,k) = dreal(ff)
|
|
IU(i,j,k) = dimag(ff)
|
|
|
|
Jr = -2.d0*beta0-C1*hgr-II*nu*C2*hgr**2-C2*hgr**3
|
|
rf = dreal(Jr*(dcos(nu*T)+II*dsin(nu*T)))
|
|
Jr = Yslm(1,2,m,gt,gp)
|
|
ff = dsqrt(dble(2*(2+1)))*rf*Jr*swtf !! U_,r
|
|
Jr = (2.4d1*beta0+3.d0*II*nu*C1-II*nu**3*C2)/3.6d1+C1/4.d0*hgr-C2/1.2d1*hgr**3
|
|
rf = dreal(Jr*(dcos(nu*T)+II*dsin(nu*T)))
|
|
Jr = Yslm(2,2,m,gt,gp)
|
|
Jr = dsqrt(dble((2-1)*2*(2+1)*(2+2)))*rf*Jr*swtf**2 !! J
|
|
rf = dsqrt(1.d0+abs(Jr)**2) !! K
|
|
ff = dexp(-2.d0*beta(i,j,k))*(Jr*dconjg(ff)+rf*ff)
|
|
RQ(i,j,k) = dreal(ff)
|
|
IQ(i,j,k) = dimag(ff)
|
|
|
|
enddo
|
|
enddo
|
|
|
|
return
|
|
|
|
end subroutine get_null_boundary
|
|
|
|
#endif
|
|
|
|
#if 0
|
|
! for checking dderivs_eth
|
|
!-------------------------------------------------------------
|
|
! Linear wave given in CQG 24S327
|
|
!-------------------------------------------------------------
|
|
subroutine get_exact_null(ex,crho,sigma,R,RJ,IJ,sst,Rmin,T, &
|
|
quR1,quR2,quI1,quI2,qlR1,qlR2,qlI1,qlI2, &
|
|
gR,gI, &
|
|
dquR1,dquR2,dquI1,dquI2, &
|
|
bdquR1,bdquR2,bdquI1,bdquI2, &
|
|
dgR,dgI,bdgR,bdgI)
|
|
implicit none
|
|
! argument variables
|
|
integer, intent(in ):: ex(1:3),sst
|
|
real*8,intent(in) :: Rmin,T
|
|
double precision,intent(in),dimension(ex(1))::crho
|
|
double precision,intent(in),dimension(ex(2))::sigma
|
|
double precision,intent(in),dimension(ex(3))::R
|
|
real*8,dimension(ex(1),ex(2),ex(3)),intent(out)::RJ,IJ
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: quR1,quR2,quI1,quI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: qlR1,qlR2,qlI1,qlI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: gR,gI
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: dquR1,dquR2,dquI1,dquI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: bdquR1,bdquR2,bdquI1,bdquI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: dgR,dgI,bdgR,bdgI
|
|
|
|
double complex,dimension(ex(1),ex(2)) :: Y20,ddY20,f
|
|
integer :: i,j,k
|
|
real*8 ::x,y,z,hgr,gt,gp,tgrho,tgsigma
|
|
double complex :: Yslm,II,Jr
|
|
|
|
double complex :: beta0,C1,C2
|
|
integer :: nu,m
|
|
|
|
call initial_null_paramter(beta0,C1,C2,nu,m)
|
|
|
|
II = dcmplx(0.d0,1.d0)
|
|
|
|
do i=1,ex(1)
|
|
do j=1,ex(2)
|
|
do k=1,ex(3)
|
|
! fake global coordinate is enough here
|
|
hgr = 1.d0
|
|
tgrho = dtan(crho(i))
|
|
tgsigma = dtan(sigma(j))
|
|
select case (sst)
|
|
case (0)
|
|
z = hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = z*tgrho
|
|
y = z*tgsigma
|
|
case (1)
|
|
z = -hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = z*tgrho
|
|
y = z*tgsigma
|
|
case (2)
|
|
x = hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
y = x*tgrho
|
|
z = x*tgsigma
|
|
case (3)
|
|
x = -hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
y = x*tgrho
|
|
z = x*tgsigma
|
|
case (4)
|
|
y = hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = y*tgrho
|
|
z = y*tgsigma
|
|
case (5)
|
|
y = -hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = y*tgrho
|
|
z = y*tgsigma
|
|
case default
|
|
write(*,*) "get_exact_null: not recognized sst = ",sst
|
|
return
|
|
end select
|
|
gt = dacos(z/hgr)
|
|
gp = datan2(y,x)
|
|
|
|
hgr = (1.d0-R(k))/R(k)/Rmin
|
|
Jr = (2.4d1*beta0+3.d0*II*nu*C1-II*nu**3*C2)/3.6d1+C1/4.d0*hgr-C2/1.2d1*hgr**3
|
|
Y20(i,j) = Yslm(0,2,m,gt,gp)
|
|
RJ(i,j,k) = dreal(Jr*(dcos(nu*T)+II*dsin(nu*T)))
|
|
IJ(i,j,k) = 0.d0
|
|
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
k=1
|
|
call dderivs_eth(ex(1:2),crho,sigma,Y20,ddY20,0,1,1, &
|
|
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k),&
|
|
dquR1(:,:,k),dquR2(:,:,k),dquI1(:,:,k),dquI2(:,:,k), &
|
|
bdquR1(:,:,k),bdquR2(:,:,k),bdquI1(:,:,k),bdquI2(:,:,k), &
|
|
dgR(:,:,k),dgI(:,:,k),bdgR(:,:,k),bdgI(:,:,k))
|
|
|
|
do k=1,ex(3)
|
|
f = ddY20*RJ(:,:,k)
|
|
RJ(:,:,k) = dreal(f)
|
|
IJ(:,:,k) = dimag(f)
|
|
enddo
|
|
|
|
return
|
|
|
|
end subroutine get_exact_null
|
|
|
|
#else
|
|
!-------------------------------------------------------------
|
|
! Linear wave given in CQG 24S327
|
|
!-------------------------------------------------------------
|
|
subroutine get_exact_null(ex,crho,sigma,R,RJ,IJ,sst,Rmin,T, &
|
|
quR1,quR2,quI1,quI2,qlR1,qlR2,qlI1,qlI2, &
|
|
gR,gI, &
|
|
dquR1,dquR2,dquI1,dquI2, &
|
|
bdquR1,bdquR2,bdquI1,bdquI2, &
|
|
dgR,dgI,bdgR,bdgI)
|
|
implicit none
|
|
! argument variables
|
|
integer, intent(in ):: ex(1:3),sst
|
|
real*8,intent(in) :: Rmin,T
|
|
double precision,intent(in),dimension(ex(1))::crho
|
|
double precision,intent(in),dimension(ex(2))::sigma
|
|
double precision,intent(in),dimension(ex(3))::R
|
|
real*8,dimension(ex(1),ex(2),ex(3)),intent(out)::RJ,IJ
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: quR1,quR2,quI1,quI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: qlR1,qlR2,qlI1,qlI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: gR,gI
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: dquR1,dquR2,dquI1,dquI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: bdquR1,bdquR2,bdquI1,bdquI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: dgR,dgI,bdgR,bdgI
|
|
|
|
integer :: i,j,k
|
|
real*8 ::x,y,z,hgr,gt,gp,tgrho,tgsigma,tc,ts
|
|
double complex :: Yslm,II,Jr,swtf,ff
|
|
|
|
double complex :: beta0,C1,C2
|
|
integer :: nu,m
|
|
|
|
call initial_null_paramter(beta0,C1,C2,nu,m)
|
|
|
|
II = dcmplx(0.d0,1.d0)
|
|
|
|
do i=1,ex(1)
|
|
do j=1,ex(2)
|
|
do k=1,ex(3)
|
|
! fake global coordinate is enough here
|
|
hgr = 1.d0
|
|
tgrho = dtan(crho(i))
|
|
tgsigma = dtan(sigma(j))
|
|
tc = dsqrt((1.d0-dsin(crho(i))*dsin(sigma(j)))/2.d0)
|
|
ts = dsqrt((1.d0+dsin(crho(i))*dsin(sigma(j)))/2.d0)
|
|
select case (sst)
|
|
case (0)
|
|
z = hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = z*tgrho
|
|
y = z*tgsigma
|
|
case (1)
|
|
z = -hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = z*tgrho
|
|
y = z*tgsigma
|
|
case (2)
|
|
x = hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
y = x*tgrho
|
|
z = x*tgsigma
|
|
case (3)
|
|
x = -hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
y = x*tgrho
|
|
z = x*tgsigma
|
|
case (4)
|
|
y = hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = y*tgrho
|
|
z = y*tgsigma
|
|
case (5)
|
|
y = -hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = y*tgrho
|
|
z = y*tgsigma
|
|
case default
|
|
write(*,*) "get_exact_null: not recognized sst = ",sst
|
|
return
|
|
end select
|
|
gt = dacos(z/hgr)
|
|
gp = datan2(y,x)
|
|
swtf = 2.d0*tc*ts*(ts+II*tc)/dcos(sigma(j))
|
|
if(sst==1 .or. sst==3 .or. sst==4) swtf = dconjg(swtf)
|
|
select case (sst)
|
|
case (0,1)
|
|
swtf = swtf/(dcos(gp)+II*dcos(gt)*dsin(gp))*(dcos(gt)**2+dsin(gt)**2*dcos(gp)**2)
|
|
case (2,3)
|
|
swtf = II*swtf*dsin(gt)
|
|
case (4,5)
|
|
swtf = -II*swtf*dsin(gt)
|
|
end select
|
|
|
|
hgr = (1.d0-R(k))/R(k)/Rmin
|
|
Jr = (2.4d1*beta0+3.d0*II*nu*C1-II*nu**3*C2)/3.6d1+C1/4.d0*hgr-C2/1.2d1*hgr**3
|
|
hgr = dreal(Jr*(dcos(nu*T)+II*dsin(nu*T)))
|
|
Jr = Yslm(2,2,m,gt,gp)
|
|
ff = dsqrt(dble((2-1)*2*(2+1)*(2+2)))*hgr*Jr*swtf**2
|
|
RJ(i,j,k) = dreal(ff)
|
|
IJ(i,j,k) = dimag(ff)
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
return
|
|
|
|
end subroutine get_exact_null
|
|
|
|
#endif
|
|
|
|
#if 0
|
|
! for checking derivs_eth and dderivs_eth
|
|
!-------------------------
|
|
subroutine get_null_boundary_c(ex,crho,sigma,R,beta,RQ,IQ,RU,IU,W,RTheta,ITheta, &
|
|
quR1,quR2,quI1,quI2,qlR1,qlR2,qlI1,qlI2, &
|
|
gR,gI, &
|
|
dquR1,dquR2,dquI1,dquI2, &
|
|
bdquR1,bdquR2,bdquI1,bdquI2, &
|
|
dgR,dgI,bdgR,bdgI,T,Rmin,sst)
|
|
|
|
implicit none
|
|
|
|
!~~~~~~% Input parameters:
|
|
integer,intent(in) :: ex(3),sst
|
|
real*8,intent(in) :: T,Rmin
|
|
double precision,intent(in),dimension(ex(1))::crho
|
|
double precision,intent(in),dimension(ex(2))::sigma
|
|
double precision,intent(in),dimension(ex(3))::R
|
|
real*8,dimension(ex(1),ex(2),ex(3)),intent(inout) :: beta,RQ,IQ,RU,IU
|
|
real*8,dimension(ex(1),ex(2),ex(3)),intent(inout) :: W,RTheta,ITheta
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: quR1,quR2,quI1,quI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: qlR1,qlR2,qlI1,qlI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: gR,gI
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: dquR1,dquR2,dquI1,dquI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: bdquR1,bdquR2,bdquI1,bdquI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: dgR,dgI,bdgR,bdgI
|
|
|
|
double complex,dimension(ex(1),ex(2)) :: Y20,dY20,ddY20,f
|
|
integer :: i,j,k
|
|
real*8 ::x,y,z,hgr,gt,gp,tgrho,tgsigma
|
|
double complex :: Yslm,II,Jr
|
|
|
|
double complex :: beta0,C1,C2
|
|
integer :: nu,m
|
|
|
|
call initial_null_paramter(beta0,C1,C2,nu,m)
|
|
|
|
II = dcmplx(0.d0,1.d0)
|
|
|
|
! write(*,*) abs(II) confirms abs == cabs
|
|
|
|
do k=1,ex(3)
|
|
do i=1,ex(1)
|
|
do j=1,ex(2)
|
|
! fake global coordinate is enough here
|
|
hgr = 1.d0
|
|
tgrho = dtan(crho(i))
|
|
tgsigma = dtan(sigma(j))
|
|
select case (sst)
|
|
case (0)
|
|
z = hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = z*tgrho
|
|
y = z*tgsigma
|
|
case (1)
|
|
z = -hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = z*tgrho
|
|
y = z*tgsigma
|
|
case (2)
|
|
x = hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
y = x*tgrho
|
|
z = x*tgsigma
|
|
case (3)
|
|
x = -hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
y = x*tgrho
|
|
z = x*tgsigma
|
|
case (4)
|
|
y = hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = y*tgrho
|
|
z = y*tgsigma
|
|
case (5)
|
|
y = -hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = y*tgrho
|
|
z = y*tgsigma
|
|
case default
|
|
write(*,*) "get_null_boundary_c: not recognized sst = ",sst
|
|
return
|
|
end select
|
|
gt = dacos(z/hgr)
|
|
gp = datan2(y,x)
|
|
|
|
hgr = (1.d0-R(k))/R(k)/Rmin
|
|
Jr = (2.4d1*beta0+3.d0*II*nu*C1-II*nu**3*C2)/3.6d1+C1/4.d0*hgr-C2/1.2d1*hgr**3
|
|
RTheta(i,j,k) = dreal(Jr*nu*(II*dcos(nu*T)-dsin(nu*T)))
|
|
f(i,j) = dreal(Jr*(dcos(nu*T)+II*dsin(nu*T)))
|
|
|
|
Jr = (-2.4d1*II*nu*beta0+3.d0*nu*nu*C1-nu**4*C2)/36.d0+2.d0*beta0*hgr&
|
|
+C1/2.d0*hgr*hgr+II*nu*C2/3.d0*hgr**3+C2/4.d0*hgr**4
|
|
RU(i,j,k) = dreal(Jr*(dcos(nu*T)+II*dsin(nu*T)))
|
|
|
|
! Re(r^2*Ul_,r*exp(i nu T)) of CQG 24S327, (12) indeed
|
|
Jr = -2.d0*beta0-C1*hgr-II*nu*C2*hgr**2-C2*hgr**3
|
|
RQ(i,j,k) = dreal(Jr*(dcos(nu*T)+II*dsin(nu*T)))
|
|
|
|
Jr = (2.4d1*II*nu*beta0-3.d0*nu*nu*C1+nu**4*C2)/6.d0+(3.d0*II*nu*C1-6.d0*beta0-II*nu**3*C2)/3.d0*hgr&
|
|
-nu*nu*C2*hgr*hgr+II*nu*C2*hgr**3+C2/2.d0*hgr**4
|
|
W(i,j,k) = dreal(Jr*(dcos(nu*T)+II*dsin(nu*T)))
|
|
|
|
Y20(i,j) = Yslm(0,2,m,gt,gp)
|
|
|
|
enddo
|
|
enddo
|
|
|
|
call derivs_eth(ex(1:2),crho,sigma,Y20,dY20,0,1, &
|
|
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
|
|
call dderivs_eth(ex(1:2),crho,sigma,Y20,ddY20,0,1,1, &
|
|
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k),&
|
|
dquR1(:,:,k),dquR2(:,:,k),dquI1(:,:,k),dquI2(:,:,k), &
|
|
bdquR1(:,:,k),bdquR2(:,:,k),bdquI1(:,:,k),bdquI2(:,:,k), &
|
|
dgR(:,:,k),dgI(:,:,k),bdgR(:,:,k),bdgI(:,:,k))
|
|
|
|
beta(:,:,k) = dreal(beta0*cdexp(II*nu*T))*dreal(Y20)
|
|
W(:,:,k) = W(:,:,k)*dreal(Y20)
|
|
|
|
f = dexp(-2.d0*beta(:,:,k))*(f*ddY20*RQ(:,:,k)*dconjg(dY20)+dsqrt(1.d0+abs(f*ddY20))*RQ(:,:,k)*dY20)
|
|
RQ(:,:,k) = dreal(f)
|
|
IQ(:,:,k) = dimag(f)
|
|
|
|
f = ddY20*RTheta(:,:,k)
|
|
RTheta(:,:,k) = dreal(f)
|
|
ITheta(:,:,k) = dimag(f)
|
|
|
|
f = dY20*RU(:,:,k)
|
|
RU(:,:,k) = dreal(f)
|
|
IU(:,:,k) = dimag(f)
|
|
enddo
|
|
|
|
return
|
|
|
|
end subroutine get_null_boundary_c
|
|
|
|
#else
|
|
|
|
!-------------------------
|
|
subroutine get_null_boundary_c(ex,crho,sigma,R,beta,RQ,IQ,RU,IU,W,RTheta,ITheta, &
|
|
quR1,quR2,quI1,quI2,qlR1,qlR2,qlI1,qlI2, &
|
|
gR,gI, &
|
|
dquR1,dquR2,dquI1,dquI2, &
|
|
bdquR1,bdquR2,bdquI1,bdquI2, &
|
|
dgR,dgI,bdgR,bdgI,T,Rmin,sst)
|
|
|
|
implicit none
|
|
|
|
!~~~~~~% Input parameters:
|
|
integer,intent(in) :: ex(3),sst
|
|
real*8,intent(in) :: T,Rmin
|
|
double precision,intent(in),dimension(ex(1))::crho
|
|
double precision,intent(in),dimension(ex(2))::sigma
|
|
double precision,intent(in),dimension(ex(3))::R
|
|
real*8,dimension(ex(1),ex(2),ex(3)),intent(inout) :: beta,RQ,IQ,RU,IU
|
|
real*8,dimension(ex(1),ex(2),ex(3)),intent(inout) :: W,RTheta,ITheta
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: quR1,quR2,quI1,quI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: qlR1,qlR2,qlI1,qlI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: gR,gI
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: dquR1,dquR2,dquI1,dquI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: bdquR1,bdquR2,bdquI1,bdquI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: dgR,dgI,bdgR,bdgI
|
|
|
|
integer :: i,j,k
|
|
real*8 ::x,y,z,hgr,gt,gp,tgrho,tgsigma,tc,ts,rf
|
|
double complex :: Yslm,II,Jr,swtf,ff
|
|
|
|
double complex :: beta0,C1,C2
|
|
integer :: nu,m
|
|
|
|
#if 0
|
|
real*8 :: betax,KK,KKx,Wx
|
|
double complex :: CJ,DCJ,CJx,CJxx,DCJx,CU,CUx,DCU,DCUx,bDCU,bDCUx,CB,DCB,bDCB,CBx
|
|
double complex :: Cnu,Cnux,Ck,fCTheta,fCThetax,Theta_rhs
|
|
|
|
|
|
T=0.25d0
|
|
i=1
|
|
j=1
|
|
k=1
|
|
hgr = 1.d0
|
|
beta(i,j,k) = dreal(beta0*cdexp(II*nu*T))
|
|
CB = beta(i,j,k)
|
|
DCB = CB
|
|
bDCB = CB
|
|
betax = 0.d0
|
|
Jr = (2.4d1*II*nu*beta0-3.d0*nu*nu*C1+nu**4*C2)/6.d0+(3.d0*II*nu*C1-6.d0*beta0-II*nu**3*C2)/3.d0/hgr&
|
|
-nu*nu*C2/hgr/hgr+II*nu*C2/hgr**3+C2/2.d0/hgr**4
|
|
W(i,j,k) = dreal(Jr*cdexp(II*nu*T))
|
|
Jr = -(3.d0*II*nu*C1-6.d0*beta0-II*nu**3*C2)/3.d0/hgr**2&
|
|
+2.d0*nu*nu*C2/hgr**3-3.d0*II*nu*C2/hgr**4-2.d0*C2/hgr**5
|
|
Wx = dreal(Jr*cdexp(II*nu*T))*(Rmin+hgr)**2/Rmin
|
|
Jr = (2.4d1*beta0+3.d0*II*nu*C1-II*nu**3*C2)/3.6d1+C1/4.d0/hgr-C2/1.2d1/hgr**3
|
|
CJ = dreal(Jr*cdexp(II*nu*T))
|
|
fCTheta = dreal(Jr*II*nu*cdexp(II*nu*T))
|
|
DCJ=CJ
|
|
KK = dsqrt(1.d0+cdabs(CJ)**2)
|
|
Jr = -C1/4.d0/hgr**2+C2/4.d0/hgr**4
|
|
rf = dreal(Jr*cdexp(II*nu*T))
|
|
CJx = rf*(Rmin+hgr)**2/Rmin
|
|
fCThetax = dreal(Jr*II*nu*cdexp(II*nu*T))*(Rmin+hgr)**2/Rmin
|
|
Jr = C1/2.d0/hgr**3-C2/hgr**5
|
|
rf = dreal(Jr*cdexp(II*nu*T))
|
|
CJxx = rf*(Rmin+hgr)**4/Rmin**2+2.d0*(Rmin+hgr)/Rmin*CJx
|
|
DCJx = CJx
|
|
KKx = dreal(CJ*dconjg(CJx))/KK
|
|
Jr = (-2.4d1*II*nu*beta0+3.d0*nu*nu*C1-nu**4*C2)/36.d0+2.d0*beta0/hgr&
|
|
+C1/2.d0/hgr/hgr+II*nu*C2/3.d0/hgr**3+C2/4.d0/hgr**4
|
|
CU = dreal(Jr*cdexp(II*nu*T))
|
|
bDCU = CU
|
|
DCU=CU
|
|
Jr = -2.d0*beta0/hgr/hgr-C1/hgr**3-II*nu*C2/hgr**4-C2/hgr**5
|
|
rf = dreal(Jr*cdexp(II*nu*T))
|
|
CUx = rf*(Rmin+hgr)**2/Rmin
|
|
DCUx=CUx
|
|
bDCUx=CUx
|
|
|
|
Cnu = CJ
|
|
Cnux = CJx
|
|
Ck = 0.d0
|
|
hgr = hgr/(Rmin+hgr)
|
|
|
|
|
|
call getndxs(T,crho(i),sigma(j),hgr,beta(i,j,k),KK,CU,bDCU,DCU, &
|
|
CB,DCB,W(i,j,k),CJ,DCJ,bDCB,Cnu,Ck,fCTheta,sst,Rmin)
|
|
call getdxs(T,crho(i),sigma(j),hgr,betax,KKx,CUx,DCUx,bDCUx, &
|
|
Wx,CJx,CJxx,DCJx,Cnux,fCThetax,sst,Rmin)
|
|
! write(*,*) 2.d0*hgr*(1.d0-hgr)*fCThetax-(-(hgr*(1-hgr)*DCUx+2.d0*DCU)+2.d0/hgr/Rmin*(1.d0-hgr)*DCB &
|
|
! +(1.d0-hgr)**3/Rmin*(2.d0*CJx+hgr*CJxx)-2.d0*fCTheta)
|
|
! stop
|
|
write(*,*) fCThetax-Theta_rhs(hgr,Rmin,beta(i,j,k),betax,KK,KKx,CU,CUx,DCUx,bDCU,bDCUx, &
|
|
DCU,CB,DCB,W(i,j,k),Wx,CJ,DCJ, &
|
|
CJx,CJxx,DCJx,bDCB,Cnu,Cnux,Ck,fCTheta)
|
|
stop
|
|
#endif
|
|
call initial_null_paramter(beta0,C1,C2,nu,m)
|
|
|
|
II = dcmplx(0.d0,1.d0)
|
|
|
|
do i=1,ex(1)
|
|
do j=1,ex(2)
|
|
do k=1,ex(3)
|
|
!fake global coordinate is enough here
|
|
hgr = 1.d0
|
|
tgrho = dtan(crho(i))
|
|
tgsigma = dtan(sigma(j))
|
|
tc = dsqrt((1.d0-dsin(crho(i))*dsin(sigma(j)))/2.d0)
|
|
ts = dsqrt((1.d0+dsin(crho(i))*dsin(sigma(j)))/2.d0)
|
|
select case (sst)
|
|
case (0)
|
|
z = hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = z*tgrho
|
|
y = z*tgsigma
|
|
case (1)
|
|
z = -hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = z*tgrho
|
|
y = z*tgsigma
|
|
case (2)
|
|
x = hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
y = x*tgrho
|
|
z = x*tgsigma
|
|
case (3)
|
|
x = -hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
y = x*tgrho
|
|
z = x*tgsigma
|
|
case (4)
|
|
y = hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = y*tgrho
|
|
z = y*tgsigma
|
|
case (5)
|
|
y = -hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = y*tgrho
|
|
z = y*tgsigma
|
|
case default
|
|
write(*,*) "get_null_boundary: not recognized sst = ",sst
|
|
return
|
|
end select
|
|
gt = dacos(z/hgr)
|
|
gp = datan2(y,x)
|
|
swtf = 2.d0*tc*ts*(ts+II*tc)/dcos(sigma(j))
|
|
if(sst==1 .or. sst==3 .or. sst==4) swtf = dconjg(swtf)
|
|
select case (sst)
|
|
case (0,1)
|
|
swtf = swtf/(dcos(gp)+II*dcos(gt)*dsin(gp))*(dcos(gt)**2+dsin(gt)**2*dcos(gp)**2)
|
|
case (2,3)
|
|
swtf = II*swtf*dsin(gt)
|
|
case (4,5)
|
|
swtf = -II*swtf*dsin(gt)
|
|
end select
|
|
|
|
hgr = (1.d0-R(k))/R(k)/Rmin
|
|
Jr = (2.4d1*beta0+3.d0*II*nu*C1-II*nu**3*C2)/3.6d1+C1/4.d0*hgr-C2/1.2d1*hgr**3
|
|
rf = dreal(Jr*nu*II*cdexp(II*nu*T))
|
|
Jr = Yslm(2,2,m,gt,gp)*swtf**2
|
|
ff = dsqrt(dble((2-1)*2*(2+1)*(2+2)))*rf*Jr
|
|
RTheta(i,j,k) = dreal(ff)
|
|
ITheta(i,j,k) = dimag(ff)
|
|
|
|
rf = dreal(Yslm(0,2,m,gt,gp))
|
|
beta(i,j,k) = rf*dreal(beta0*cdexp(II*nu*T))
|
|
Jr = (2.4d1*II*nu*beta0-3.d0*nu*nu*C1+nu**4*C2)/6.d0+(3.d0*II*nu*C1-6.d0*beta0-II*nu**3*C2)/3.d0*hgr&
|
|
-nu*nu*C2*hgr*hgr+II*nu*C2*hgr**3+C2/2.d0*hgr**4
|
|
W(i,j,k) = rf*dreal(Jr*cdexp(II*nu*T))
|
|
|
|
Jr = (-2.4d1*II*nu*beta0+3.d0*nu*nu*C1-nu**4*C2)/36.d0+2.d0*beta0*hgr&
|
|
+C1/2.d0*hgr*hgr+II*nu*C2/3.d0*hgr**3+C2/4.d0*hgr**4
|
|
rf = dreal(Jr*cdexp(II*nu*T))
|
|
Jr = Yslm(1,2,m,gt,gp)*swtf
|
|
ff = dsqrt(dble(2*(2+1)))*rf*Jr
|
|
RU(i,j,k) = dreal(ff)
|
|
IU(i,j,k) = dimag(ff)
|
|
|
|
Jr = -2.d0*beta0-C1*hgr-II*nu*C2*hgr**2-C2*hgr**3
|
|
rf = dreal(Jr*cdexp(II*nu*T))
|
|
Jr = Yslm(1,2,m,gt,gp)*swtf
|
|
ff = dsqrt(dble(2*(2+1)))*rf*Jr !! U_,r
|
|
Jr = (2.4d1*beta0+3.d0*II*nu*C1-II*nu**3*C2)/3.6d1+C1/4.d0*hgr-C2/1.2d1*hgr**3
|
|
rf = dreal(Jr*cdexp(II*nu*T))
|
|
Jr = Yslm(2,2,m,gt,gp)*swtf**2
|
|
Jr = dsqrt(dble((2-1)*2*(2+1)*(2+2)))*rf*Jr !! J
|
|
rf = dsqrt(1.d0+cdabs(Jr)**2) !! K
|
|
ff = dexp(-2.d0*beta(i,j,k))*(Jr*dconjg(ff)+rf*ff)
|
|
RQ(i,j,k) = dreal(ff)
|
|
IQ(i,j,k) = dimag(ff)
|
|
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
return
|
|
|
|
end subroutine get_null_boundary_c
|
|
#endif
|
|
|
|
!==========================================================
|
|
subroutine initial_null_paramter(beta0,C1,C2,nu,m)
|
|
|
|
implicit none
|
|
|
|
double complex,intent(out) :: beta0,C1,C2
|
|
integer,intent(out) :: nu,m
|
|
|
|
nu=1
|
|
m=0
|
|
beta0 = dcmplx(0.d0,1.d-6)
|
|
C1 = dcmplx(3.d-6,0.d0)
|
|
C2 = dcmplx(1.d-6,0.d0)
|
|
|
|
end subroutine initial_null_paramter
|
|
|
|
#if 1
|
|
subroutine get_exact_null_theta(ex,crho,sigma,R,RTheta,ITheta,sst,Rmin,T, &
|
|
quR1,quR2,quI1,quI2,qlR1,qlR2,qlI1,qlI2, &
|
|
gR,gI, &
|
|
dquR1,dquR2,dquI1,dquI2, &
|
|
bdquR1,bdquR2,bdquI1,bdquI2, &
|
|
dgR,dgI,bdgR,bdgI)
|
|
implicit none
|
|
! argument variables
|
|
integer, intent(in ):: ex(1:3),sst
|
|
real*8,intent(in) :: Rmin,T
|
|
double precision,intent(in),dimension(ex(1))::crho
|
|
double precision,intent(in),dimension(ex(2))::sigma
|
|
double precision,intent(in),dimension(ex(3))::R
|
|
real*8,dimension(ex(1),ex(2),ex(3)),intent(out)::RTheta,ITheta
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: quR1,quR2,quI1,quI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: qlR1,qlR2,qlI1,qlI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: gR,gI
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: dquR1,dquR2,dquI1,dquI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: bdquR1,bdquR2,bdquI1,bdquI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: dgR,dgI,bdgR,bdgI
|
|
|
|
integer :: i,j,k
|
|
real*8 ::x,y,z,hgr,gt,gp,tgrho,tgsigma,tc,ts
|
|
double complex :: Yslm,II,Jr,swtf,ff
|
|
|
|
double complex :: beta0,C1,C2
|
|
integer :: nu,m
|
|
|
|
call initial_null_paramter(beta0,C1,C2,nu,m)
|
|
|
|
II = dcmplx(0.d0,1.d0)
|
|
|
|
do i=1,ex(1)
|
|
do j=1,ex(2)
|
|
do k=1,ex(3)
|
|
! fake global coordinate is enough here
|
|
hgr = 1.d0
|
|
tgrho = dtan(crho(i))
|
|
tgsigma = dtan(sigma(j))
|
|
tc = dsqrt((1.d0-dsin(crho(i))*dsin(sigma(j)))/2.d0)
|
|
ts = dsqrt((1.d0+dsin(crho(i))*dsin(sigma(j)))/2.d0)
|
|
select case (sst)
|
|
case (0)
|
|
z = hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = z*tgrho
|
|
y = z*tgsigma
|
|
case (1)
|
|
z = -hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = z*tgrho
|
|
y = z*tgsigma
|
|
case (2)
|
|
x = hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
y = x*tgrho
|
|
z = x*tgsigma
|
|
case (3)
|
|
x = -hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
y = x*tgrho
|
|
z = x*tgsigma
|
|
case (4)
|
|
y = hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = y*tgrho
|
|
z = y*tgsigma
|
|
case (5)
|
|
y = -hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = y*tgrho
|
|
z = y*tgsigma
|
|
case default
|
|
write(*,*) "get_exact_null_theta: not recognized sst = ",sst
|
|
return
|
|
end select
|
|
gt = dacos(z/hgr)
|
|
gp = datan2(y,x)
|
|
swtf = 2.d0*tc*ts*(ts+II*tc)/dcos(sigma(j))
|
|
if(sst==1 .or. sst==3 .or. sst==4) swtf = dconjg(swtf)
|
|
select case (sst)
|
|
case (0,1)
|
|
swtf = swtf/(dcos(gp)+II*dcos(gt)*dsin(gp))*(dcos(gt)**2+dsin(gt)**2*dcos(gp)**2)
|
|
case (2,3)
|
|
swtf = II*swtf*dsin(gt)
|
|
case (4,5)
|
|
swtf = -II*swtf*dsin(gt)
|
|
end select
|
|
|
|
hgr = (1.d0-R(k))/R(k)/Rmin
|
|
Jr = (2.4d1*beta0+3.d0*II*nu*C1-II*nu**3*C2)/3.6d1+C1/4.d0*hgr-C2/1.2d1*hgr**3
|
|
hgr = dreal(Jr*nu*(-dsin(nu*T)+II*dcos(nu*T)))
|
|
Jr = Yslm(2,2,m,gt,gp)
|
|
ff = dsqrt(dble((2-1)*2*(2+1)*(2+2)))*hgr*Jr*swtf**2
|
|
RTheta(i,j,k) = dreal(ff)
|
|
ITheta(i,j,k) = dimag(ff)
|
|
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
return
|
|
|
|
end subroutine get_exact_null_theta
|
|
!------------------------------------------------------------------------------------------------
|
|
subroutine get_exact_null_theta_x(ex,crho,sigma,R,RThetax,IThetax,sst,Rmin,T, &
|
|
quR1,quR2,quI1,quI2,qlR1,qlR2,qlI1,qlI2, &
|
|
gR,gI, &
|
|
dquR1,dquR2,dquI1,dquI2, &
|
|
bdquR1,bdquR2,bdquI1,bdquI2, &
|
|
dgR,dgI,bdgR,bdgI)
|
|
implicit none
|
|
! argument variables
|
|
integer, intent(in ):: ex(1:3),sst
|
|
real*8,intent(in) :: Rmin,T
|
|
double precision,intent(in),dimension(ex(1))::crho
|
|
double precision,intent(in),dimension(ex(2))::sigma
|
|
double precision,intent(in),dimension(ex(3))::R
|
|
real*8,dimension(ex(1),ex(2),ex(3)),intent(out)::RThetax,IThetax
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: quR1,quR2,quI1,quI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: qlR1,qlR2,qlI1,qlI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: gR,gI
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: dquR1,dquR2,dquI1,dquI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: bdquR1,bdquR2,bdquI1,bdquI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: dgR,dgI,bdgR,bdgI
|
|
|
|
integer :: i,j,k
|
|
real*8 ::x,y,z,hgr,gt,gp,tgrho,tgsigma,tc,ts
|
|
double complex :: Yslm,II,Jr,swtf,ff
|
|
|
|
double complex :: beta0,C1,C2
|
|
integer :: nu,m
|
|
|
|
call initial_null_paramter(beta0,C1,C2,nu,m)
|
|
|
|
II = dcmplx(0.d0,1.d0)
|
|
|
|
do i=1,ex(1)
|
|
do j=1,ex(2)
|
|
do k=1,ex(3)
|
|
! fake global coordinate is enough here
|
|
hgr = 1.d0
|
|
tgrho = dtan(crho(i))
|
|
tgsigma = dtan(sigma(j))
|
|
tc = dsqrt((1.d0-dsin(crho(i))*dsin(sigma(j)))/2.d0)
|
|
ts = dsqrt((1.d0+dsin(crho(i))*dsin(sigma(j)))/2.d0)
|
|
select case (sst)
|
|
case (0)
|
|
z = hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = z*tgrho
|
|
y = z*tgsigma
|
|
case (1)
|
|
z = -hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = z*tgrho
|
|
y = z*tgsigma
|
|
case (2)
|
|
x = hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
y = x*tgrho
|
|
z = x*tgsigma
|
|
case (3)
|
|
x = -hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
y = x*tgrho
|
|
z = x*tgsigma
|
|
case (4)
|
|
y = hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = y*tgrho
|
|
z = y*tgsigma
|
|
case (5)
|
|
y = -hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = y*tgrho
|
|
z = y*tgsigma
|
|
case default
|
|
write(*,*) "get_exact_null_theta_x: not recognized sst = ",sst
|
|
return
|
|
end select
|
|
gt = dacos(z/hgr)
|
|
gp = datan2(y,x)
|
|
swtf = 2.d0*tc*ts*(ts+II*tc)/dcos(sigma(j))
|
|
if(sst==1 .or. sst==3 .or. sst==4) swtf = dconjg(swtf)
|
|
select case (sst)
|
|
case (0,1)
|
|
swtf = swtf/(dcos(gp)+II*dcos(gt)*dsin(gp))*(dcos(gt)**2+dsin(gt)**2*dcos(gp)**2)
|
|
case (2,3)
|
|
swtf = II*swtf*dsin(gt)
|
|
case (4,5)
|
|
swtf = -II*swtf*dsin(gt)
|
|
end select
|
|
|
|
hgr = (1.d0-R(k))/R(k)/Rmin
|
|
Jr = C1/4.d0-C2/1.2d1*3*hgr**2
|
|
Jr = -Jr/Rmin/R(k)/R(k)
|
|
hgr = dreal(Jr*nu*(-dsin(nu*T)+II*dcos(nu*T)))
|
|
Jr = Yslm(2,2,m,gt,gp)
|
|
ff = dsqrt(dble((2-1)*2*(2+1)*(2+2)))*hgr*Jr*swtf**2
|
|
RThetax(i,j,k) = dreal(ff)
|
|
IThetax(i,j,k) = dimag(ff)
|
|
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
return
|
|
|
|
end subroutine get_exact_null_theta_x
|
|
!-------------------------
|
|
subroutine get_exact_Jul(ex,crho,sigma,R,RJul,IJul, &
|
|
quR1,quR2,quI1,quI2,qlR1,qlR2,qlI1,qlI2, &
|
|
gR,gI, &
|
|
dquR1,dquR2,dquI1,dquI2, &
|
|
bdquR1,bdquR2,bdquI1,bdquI2, &
|
|
dgR,dgI,bdgR,bdgI,T,Rmin,sst)
|
|
|
|
implicit none
|
|
|
|
!~~~~~~% Input parameters:
|
|
integer,intent(in) :: ex(3),sst
|
|
real*8,intent(in) :: T,Rmin
|
|
double precision,intent(in),dimension(ex(1))::crho
|
|
double precision,intent(in),dimension(ex(2))::sigma
|
|
double precision,intent(in),dimension(ex(3))::R
|
|
real*8,dimension(ex(1),ex(2),ex(3)),intent(inout) :: RJul,IJul
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: quR1,quR2,quI1,quI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: qlR1,qlR2,qlI1,qlI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: gR,gI
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: dquR1,dquR2,dquI1,dquI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: bdquR1,bdquR2,bdquI1,bdquI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: dgR,dgI,bdgR,bdgI
|
|
|
|
integer :: i,j,k
|
|
real*8 ::x,y,z,hgr,gt,gp,tgrho,tgsigma,tc,ts,rf
|
|
double complex :: Yslm,II,Jr,swtf,ff
|
|
|
|
double complex :: beta0,C1,C2
|
|
integer :: nu,m
|
|
|
|
call initial_null_paramter(beta0,C1,C2,nu,m)
|
|
|
|
II = dcmplx(0.d0,1.d0)
|
|
|
|
do i=1,ex(1)
|
|
do j=1,ex(2)
|
|
do k=1,ex(3)
|
|
!fake global coordinate is enough here
|
|
hgr = 1.d0
|
|
tgrho = dtan(crho(i))
|
|
tgsigma = dtan(sigma(j))
|
|
tc = dsqrt((1.d0-dsin(crho(i))*dsin(sigma(j)))/2.d0)
|
|
ts = dsqrt((1.d0+dsin(crho(i))*dsin(sigma(j)))/2.d0)
|
|
select case (sst)
|
|
case (0)
|
|
z = hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = z*tgrho
|
|
y = z*tgsigma
|
|
case (1)
|
|
z = -hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = z*tgrho
|
|
y = z*tgsigma
|
|
case (2)
|
|
x = hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
y = x*tgrho
|
|
z = x*tgsigma
|
|
case (3)
|
|
x = -hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
y = x*tgrho
|
|
z = x*tgsigma
|
|
case (4)
|
|
y = hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = y*tgrho
|
|
z = y*tgsigma
|
|
case (5)
|
|
y = -hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = y*tgrho
|
|
z = y*tgsigma
|
|
case default
|
|
write(*,*) "get_exact_Jul: not recognized sst = ",sst
|
|
return
|
|
end select
|
|
gt = dacos(z/hgr)
|
|
gp = datan2(y,x)
|
|
swtf = 2.d0*tc*ts*(ts+II*tc)/dcos(sigma(j))
|
|
if(sst==1 .or. sst==3 .or. sst==4) swtf = dconjg(swtf)
|
|
select case (sst)
|
|
case (0,1)
|
|
swtf = swtf/(dcos(gp)+II*dcos(gt)*dsin(gp))*(dcos(gt)**2+dsin(gt)**2*dcos(gp)**2)
|
|
case (2,3)
|
|
swtf = II*swtf*dsin(gt)
|
|
case (4,5)
|
|
swtf = -II*swtf*dsin(gt)
|
|
end select
|
|
|
|
hgr = (1.d0-R(k))/R(k)/Rmin
|
|
Jr = C1/4.d0-C2/4.d0*hgr**2
|
|
rf = dreal(Jr*nu*II*cdexp(II*nu*T))
|
|
Jr = Yslm(2,2,m,gt,gp)*swtf**2
|
|
ff = dsqrt(dble((2-1)*2*(2+1)*(2+2)))*rf*Jr
|
|
RJul(i,j,k) = dreal(ff)
|
|
IJul(i,j,k) = dimag(ff)
|
|
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
return
|
|
|
|
end subroutine get_exact_Jul
|
|
!-------------------------
|
|
subroutine get_fake_Ju(ex,crho,sigma,R,RJul,IJul, &
|
|
quR1,quR2,quI1,quI2,qlR1,qlR2,qlI1,qlI2, &
|
|
gR,gI, &
|
|
dquR1,dquR2,dquI1,dquI2, &
|
|
bdquR1,bdquR2,bdquI1,bdquI2, &
|
|
dgR,dgI,bdgR,bdgI,T,Rmin,sst)
|
|
|
|
implicit none
|
|
|
|
!~~~~~~% Input parameters:
|
|
integer,intent(in) :: ex(3),sst
|
|
real*8,intent(in) :: T,Rmin
|
|
double precision,intent(in),dimension(ex(1))::crho
|
|
double precision,intent(in),dimension(ex(2))::sigma
|
|
double precision,intent(in),dimension(ex(3))::R
|
|
real*8,dimension(ex(1),ex(2),ex(3)),intent(inout) :: RJul,IJul
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: quR1,quR2,quI1,quI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: qlR1,qlR2,qlI1,qlI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: gR,gI
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: dquR1,dquR2,dquI1,dquI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: bdquR1,bdquR2,bdquI1,bdquI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: dgR,dgI,bdgR,bdgI
|
|
|
|
integer :: i,j,k
|
|
real*8 ::x,y,z,hgr,gt,gp,tgrho,tgsigma,tc,ts,rf
|
|
double complex :: Yslm,II,Jr,swtf,ff
|
|
|
|
double complex :: beta0,C1,C2
|
|
integer :: nu,m
|
|
|
|
call initial_null_paramter(beta0,C1,C2,nu,m)
|
|
|
|
II = dcmplx(0.d0,1.d0)
|
|
|
|
do i=1,ex(1)
|
|
do j=1,ex(2)
|
|
do k=1,ex(3)
|
|
!fake global coordinate is enough here
|
|
hgr = 1.d0
|
|
tgrho = dtan(crho(i))
|
|
tgsigma = dtan(sigma(j))
|
|
tc = dsqrt((1.d0-dsin(crho(i))*dsin(sigma(j)))/2.d0)
|
|
ts = dsqrt((1.d0+dsin(crho(i))*dsin(sigma(j)))/2.d0)
|
|
select case (sst)
|
|
case (0)
|
|
z = hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = z*tgrho
|
|
y = z*tgsigma
|
|
case (1)
|
|
z = -hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = z*tgrho
|
|
y = z*tgsigma
|
|
case (2)
|
|
x = hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
y = x*tgrho
|
|
z = x*tgsigma
|
|
case (3)
|
|
x = -hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
y = x*tgrho
|
|
z = x*tgsigma
|
|
case (4)
|
|
y = hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = y*tgrho
|
|
z = y*tgsigma
|
|
case (5)
|
|
y = -hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = y*tgrho
|
|
z = y*tgsigma
|
|
case default
|
|
write(*,*) "get_fake_Ju: not recognized sst = ",sst
|
|
return
|
|
end select
|
|
gt = dacos(z/hgr)
|
|
gp = datan2(y,x)
|
|
swtf = 2.d0*tc*ts*(ts+II*tc)/dcos(sigma(j))
|
|
if(sst==1 .or. sst==3 .or. sst==4) swtf = dconjg(swtf)
|
|
select case (sst)
|
|
case (0,1)
|
|
swtf = swtf/(dcos(gp)+II*dcos(gt)*dsin(gp))*(dcos(gt)**2+dsin(gt)**2*dcos(gp)**2)
|
|
case (2,3)
|
|
swtf = II*swtf*dsin(gt)
|
|
case (4,5)
|
|
swtf = -II*swtf*dsin(gt)
|
|
end select
|
|
|
|
hgr = (1.d0-R(k))/R(k)/Rmin
|
|
Jr = (2.4d1*beta0+3.d0*II*nu*C1-II*nu**3*C2)/3.6d1+C1/4.d0*hgr-C2/1.2d1*hgr**3
|
|
rf = dreal(Jr*nu*II*cdexp(II*nu*T))
|
|
ff = dcmplx(rf,0.d0)
|
|
RJul(i,j,k) = dreal(ff)
|
|
IJul(i,j,k) = dimag(ff)
|
|
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
if(sst == 0 .and. crho(1) < -dacos(-1.d0)/4 .and. sigma(1) < -dacos(-1.d0)/4)then
|
|
write(*,*)"T=",T,"exp(i T)=",cdexp(II*nu*T)
|
|
write(*,*)RJul(ex(1)/2,ex(2)/2,ex(3)),RJul(ex(1)/2,ex(2)/2,ex(3)-1),R(2)-R(1)
|
|
endif
|
|
|
|
return
|
|
|
|
end subroutine get_fake_Ju
|
|
!-------------------------
|
|
subroutine get_exact_omegau(ex,crho,sigma,R,omegau, &
|
|
quR1,quR2,quI1,quI2,qlR1,qlR2,qlI1,qlI2, &
|
|
gR,gI, &
|
|
dquR1,dquR2,dquI1,dquI2, &
|
|
bdquR1,bdquR2,bdquI1,bdquI2, &
|
|
dgR,dgI,bdgR,bdgI,T,Rmin,sst)
|
|
|
|
implicit none
|
|
|
|
!~~~~~~% Input parameters:
|
|
integer,intent(in) :: ex(3),sst
|
|
real*8,intent(in) :: T,Rmin
|
|
double precision,intent(in),dimension(ex(1))::crho
|
|
double precision,intent(in),dimension(ex(2))::sigma
|
|
double precision,intent(in),dimension(ex(3))::R
|
|
real*8,dimension(ex(1),ex(2),ex(3)),intent(inout) :: omegau
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: quR1,quR2,quI1,quI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: qlR1,qlR2,qlI1,qlI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: gR,gI
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: dquR1,dquR2,dquI1,dquI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: bdquR1,bdquR2,bdquI1,bdquI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: dgR,dgI,bdgR,bdgI
|
|
|
|
integer :: i,j,k
|
|
real*8 ::x,y,z,hgr,gt,gp,tgrho,tgsigma,tc,ts,rf
|
|
double complex :: Yslm,II,Jr,swtf,ff
|
|
|
|
double complex :: beta0,C1,C2
|
|
integer :: nu,m
|
|
|
|
call initial_null_paramter(beta0,C1,C2,nu,m)
|
|
|
|
II = dcmplx(0.d0,1.d0)
|
|
|
|
do i=1,ex(1)
|
|
do j=1,ex(2)
|
|
do k=1,ex(3)
|
|
!fake global coordinate is enough here
|
|
hgr = 1.d0
|
|
tgrho = dtan(crho(i))
|
|
tgsigma = dtan(sigma(j))
|
|
tc = dsqrt((1.d0-dsin(crho(i))*dsin(sigma(j)))/2.d0)
|
|
ts = dsqrt((1.d0+dsin(crho(i))*dsin(sigma(j)))/2.d0)
|
|
select case (sst)
|
|
case (0)
|
|
z = hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = z*tgrho
|
|
y = z*tgsigma
|
|
case (1)
|
|
z = -hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = z*tgrho
|
|
y = z*tgsigma
|
|
case (2)
|
|
x = hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
y = x*tgrho
|
|
z = x*tgsigma
|
|
case (3)
|
|
x = -hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
y = x*tgrho
|
|
z = x*tgsigma
|
|
case (4)
|
|
y = hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = y*tgrho
|
|
z = y*tgsigma
|
|
case (5)
|
|
y = -hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = y*tgrho
|
|
z = y*tgsigma
|
|
case default
|
|
write(*,*) "get_exact_omegau: not recognized sst = ",sst
|
|
return
|
|
end select
|
|
gt = dacos(z/hgr)
|
|
gp = datan2(y,x)
|
|
swtf = 2.d0*tc*ts*(ts+II*tc)/dcos(sigma(j))
|
|
if(sst==1 .or. sst==3 .or. sst==4) swtf = dconjg(swtf)
|
|
select case (sst)
|
|
case (0,1)
|
|
swtf = swtf/(dcos(gp)+II*dcos(gt)*dsin(gp))*(dcos(gt)**2+dsin(gt)**2*dcos(gp)**2)
|
|
case (2,3)
|
|
swtf = II*swtf*dsin(gt)
|
|
case (4,5)
|
|
swtf = -II*swtf*dsin(gt)
|
|
end select
|
|
|
|
hgr = (1.d0-R(k))/R(k)/Rmin
|
|
Jr = (2.4d1*beta0+3.d0*II*nu*C1-II*nu**3*C2)/3.6d1+C1/4.d0*hgr-C2/1.2d1*hgr**3
|
|
rf = dreal(Jr*nu*II*cdexp(II*nu*T))
|
|
Jr = Yslm(0,2,m,gt,gp)
|
|
ff = -dble(2*(2+1))/2.d0*rf*Jr
|
|
omegau(i,j,k) = dreal(ff)
|
|
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
return
|
|
|
|
end subroutine get_exact_omegau
|
|
!-------------------------
|
|
subroutine get_exact_eth2omega(ex,crho,sigma,R,Reth2omega,Ieth2omega, &
|
|
quR1,quR2,quI1,quI2,qlR1,qlR2,qlI1,qlI2, &
|
|
gR,gI, &
|
|
dquR1,dquR2,dquI1,dquI2, &
|
|
bdquR1,bdquR2,bdquI1,bdquI2, &
|
|
dgR,dgI,bdgR,bdgI,T,Rmin,sst)
|
|
|
|
implicit none
|
|
|
|
!~~~~~~% Input parameters:
|
|
integer,intent(in) :: ex(3),sst
|
|
real*8,intent(in) :: T,Rmin
|
|
double precision,intent(in),dimension(ex(1))::crho
|
|
double precision,intent(in),dimension(ex(2))::sigma
|
|
double precision,intent(in),dimension(ex(3))::R
|
|
real*8,dimension(ex(1),ex(2),ex(3)),intent(inout) :: Reth2omega,Ieth2omega
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: quR1,quR2,quI1,quI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: qlR1,qlR2,qlI1,qlI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: gR,gI
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: dquR1,dquR2,dquI1,dquI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: bdquR1,bdquR2,bdquI1,bdquI2
|
|
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: dgR,dgI,bdgR,bdgI
|
|
|
|
integer :: i,j,k
|
|
real*8 ::x,y,z,hgr,gt,gp,tgrho,tgsigma,tc,ts,rf
|
|
double complex :: Yslm,II,Jr,swtf,ff
|
|
|
|
double complex :: beta0,C1,C2
|
|
integer :: nu,m
|
|
|
|
call initial_null_paramter(beta0,C1,C2,nu,m)
|
|
|
|
II = dcmplx(0.d0,1.d0)
|
|
|
|
do i=1,ex(1)
|
|
do j=1,ex(2)
|
|
do k=1,ex(3)
|
|
!fake global coordinate is enough here
|
|
hgr = 1.d0
|
|
tgrho = dtan(crho(i))
|
|
tgsigma = dtan(sigma(j))
|
|
tc = dsqrt((1.d0-dsin(crho(i))*dsin(sigma(j)))/2.d0)
|
|
ts = dsqrt((1.d0+dsin(crho(i))*dsin(sigma(j)))/2.d0)
|
|
select case (sst)
|
|
case (0)
|
|
z = hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = z*tgrho
|
|
y = z*tgsigma
|
|
case (1)
|
|
z = -hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = z*tgrho
|
|
y = z*tgsigma
|
|
case (2)
|
|
x = hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
y = x*tgrho
|
|
z = x*tgsigma
|
|
case (3)
|
|
x = -hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
y = x*tgrho
|
|
z = x*tgsigma
|
|
case (4)
|
|
y = hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = y*tgrho
|
|
z = y*tgsigma
|
|
case (5)
|
|
y = -hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
|
|
x = y*tgrho
|
|
z = y*tgsigma
|
|
case default
|
|
write(*,*) "get_exact_eth2omega: not recognized sst = ",sst
|
|
return
|
|
end select
|
|
gt = dacos(z/hgr)
|
|
gp = datan2(y,x)
|
|
swtf = 2.d0*tc*ts*(ts+II*tc)/dcos(sigma(j))
|
|
if(sst==1 .or. sst==3 .or. sst==4) swtf = dconjg(swtf)
|
|
select case (sst)
|
|
case (0,1)
|
|
swtf = swtf/(dcos(gp)+II*dcos(gt)*dsin(gp))*(dcos(gt)**2+dsin(gt)**2*dcos(gp)**2)
|
|
case (2,3)
|
|
swtf = II*swtf*dsin(gt)
|
|
case (4,5)
|
|
swtf = -II*swtf*dsin(gt)
|
|
end select
|
|
|
|
hgr = (1.d0-R(k))/R(k)/Rmin
|
|
Jr = (2.4d1*beta0+3.d0*II*nu*C1-II*nu**3*C2)/3.6d1+C1/4.d0*hgr-C2/1.2d1*hgr**3
|
|
rf = dreal(Jr*cdexp(II*nu*T))
|
|
Jr = Yslm(2,2,m,gt,gp)
|
|
ff = -dble(2*(2+1))/2.d0*dsqrt(dble((2-1)*2*(2+1)*(2+2)))*rf*Jr
|
|
Reth2omega(i,j,k) = dreal(ff)
|
|
Ieth2omega(i,j,k) = dimag(ff)
|
|
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
return
|
|
|
|
end subroutine get_exact_eth2omega
|
|
#endif
|
|
#endif
|