Files
AMSS-NCKU/AMSS_NCKU_source/NullEvol.f90
2026-01-13 15:01:15 +08:00

4027 lines
153 KiB
Fortran

#include "macrodef.fh"
!#define OLD
! 0: rk4, 1: Adams-Moulton
#define RKorAM 0
function beta_rhs(xx,CJx,Kx) result(gont)
implicit none
double complex,intent(in) :: CJx
real*8,intent(in) :: xx,Kx
real*8 :: gont
gont = xx*(1.d0-xx)/8.d0*(dreal(CJx*dconjg(CJx))-Kx*Kx)
return
end function beta_rhs
function Q_rhs(xx,CJ,CJx,DCJx,KK,Ck,Ckx,Cnux,KKx,CBx,Cnu,DCJ,CB,CQ) result(gont)
implicit none
double complex,intent(in) :: CJ,CJx,DCJx,Ck,Ckx,Cnux,CBx,Cnu,DCJ,CB,CQ
real*8,intent(in) :: xx,KK,KKx
double complex :: gont
gont = -KK*(Ckx+Cnux)+Cnu*KKx+CJ*dconjg(Ckx)+2.d0*CBx &
+dconjg(Cnu)*CJx+dconjg(CJ)*DCJx-dconjg(Ck)*CJx &
+(dconjg(Cnu)*(CJx-CJ*CJ*dconjg(CJx)) &
+DCJ*(dconjg(CJx)-dconjg(CJ*CJ)*CJx)/2.d0/KK/KK) &
-2.d0*(2.d0*CB+CQ)/xx/(1.d0-xx)
return
end function Q_rhs
function U_rhs(xx,Rmin,beta,KK,CQ,CJ) result(gont)
implicit none
double complex,intent(in) :: CQ,CJ
real*8,intent(in) :: xx,Rmin,beta,KK
double complex :: gont
#if 1
gont = dexp(2.d0*beta)/Rmin/xx/xx*(KK*CQ-CJ*dconjg(CQ))
#else
gont = CQ/Rmin/xx/xx
#endif
#if 0
if(cdabs(gont)>1)then
write(*,*)beta,KK,CQ,CJ
stop
endif
#endif
return
end function U_rhs
function W_rhs(xx,Rmin,beta,KK,DCB,CB,CJ,Cnu,Ck,W, &
CQ,bDCk,bDCnu,bDCB,bDCU,bDCUx,DCJ) result(gont)
implicit none
double complex,intent(in) :: DCB,CB,CJ,Cnu,Ck,CQ,bDCk,bDCnu,bDCB,bDCU,bDCUx,DCJ
real*8,intent(in) :: xx,Rmin,beta,KK,W
real*8 :: Ric,gont
Ric = dreal(2.d0*KK+bDCnu-bDCk+(DCJ*dconjg(DCJ)-Cnu*dconjg(Cnu))/4.d0/KK)
gont = dreal(dexp(2.d0*beta)*(Ric/2.d0-KK*(bDCB+CB*dconjg(CB))+dconjg(CJ)*(bDCB+CB*CB) &
+(Cnu-Ck)*dconjg(CB))-1.d0+2.d0*Rmin*xx/(1.d0-xx)*(bDCU-W) &
+Rmin*xx*xx/2.d0*bDCUx-dexp(2.d0*beta)/4.d0* &
(KK*KK-CJ*dconjg(CJ))*(KK*dconjg(CQ)-dconjg(CJ)*CQ)*CQ)
gont = gont/Rmin/xx/xx
return
end function W_rhs
function Theta_rhs(xx,Rmin,beta,betax,KK,KKx,CU,CUx,DCUx,bDCU,bDCUx,DCU,CB,DCB,W,Wx,CJ,DCJ,CJx,CJxx, &
DCJx,bDCB,Cnu,Cnux,Ck,Theta) result(gont)
implicit none
double complex,intent(in) :: CU,CUx,DCUx,bDCU,bDCUx,DCU,CB,DCB,CJ,DCJ,CJx,CJxx,DCJx
double complex,intent(in) :: bDCB,Cnu,Cnux,Ck,Theta
real*8,intent(in) :: xx,Rmin,beta,betax,KK,KKx,W,Wx
double complex :: JH,II,gont
real*8 :: V,Vx,Pu
II = dcmplx(0.d0,1.d0)
V = xx*Rmin/(1.d0-xx)*(1.d0+xx*Rmin/(1.d0-xx)*W)
Vx = Rmin/(1.d0-xx)**2+2.d0*xx*Rmin*Rmin/(1.d0-xx)**3*W+xx*xx*Rmin*Rmin/(1.d0-xx)**2*Wx
Pu = 2.d0*xx*(1.d0-xx)/KK*dreal(Theta*(dconjg(CJx)*KK-dconjg(CJ)*KKx))
JH = (1.d0-xx)/xx/Rmin*dexp(2.d0*beta)*(-KK*DCJ*dconjg(CB)+ &
(KK*Cnu+(KK*KK-1.d0)*DCJ-2.d0*KK*Ck)*CB+CJ* &
((2.d0*Ck-Cnu)*dconjg(CB)-2.d0*KK*(bDCB+CB*dconjg(CB))+ &
2.d0*dreal((Cnu-Ck)*dconjg(CB)+dconjg(CJ)*(DCB+CB*CB)))) &
+0.5d0*Rmin*xx**3*(1.d0-xx)*dexp(-2.d0*beta)* &
((KK*CUx+CJ*dconjg(CUx))**2- &
CJ*dreal(dconjg(CUx)*(KK*CUx+CJ*dconjg(CUx)))) &
-0.5d0*(Cnu*(xx*(1.d0-xx)*CUx+2.d0*CU)+DCJ*(xx*(1.d0-xx)*dconjg(CUx)+ &
2.d0*dconjg(CU)))+CJ*II*dimag(xx*(1.d0-xx)*bDCUx+2.d0*bDCU) &
-xx*(1.d0-xx)*CJx*dreal(bDCU) &
+xx*(1.d0-xx)*(dconjg(CU)*DCJ+CU*Cnu)*II*dimag(CJ*dconjg(CJx)) &
-xx*(1.d0-xx)*(dconjg(CU)*DCJx+CU*Cnux) &
-2.d0*xx*(1.d0-xx)*(CJ*KKx-KK*CJx)*(dreal(dconjg(CU)*Ck)+ &
II*dimag(KK*bDCU-dconjg(CJ)*DCU)) &
-8.d0*CJ*((1.d0-xx)**2/Rmin+xx*(1.d0-xx)*W)*betax
gont = -KK*(xx*(1-xx)*DCUx+2.d0*DCU)+2.d0*(1.d0-xx)/xx/Rmin*dexp(2.d0*beta)*(DCB+CB*CB) &
-(xx*(1.d0-xx)*Wx+W)*CJ+JH+CJ*Pu-2.d0*Theta &
-(1.d0-xx)*(1.d0-xx)/xx/xx/Rmin/Rmin*V*(CJ+xx*(1.d0-xx)*CJx) &
+(1.d0-xx)*(1.d0-xx)*(1.d0-xx)/xx/Rmin/Rmin*Vx*(CJ+xx*(1.d0-xx)*CJx) &
+(1.d0-xx)**4/xx/Rmin/Rmin*V*(2.d0*CJx+xx*CJxx)
#if 0
gont = -(xx*(1-xx)*DCUx+2.d0*DCU)+2.d0*(1.d0-xx)/xx/Rmin*DCB &
-2.d0*Theta &
+(1.d0-xx)**3/Rmin*(2.d0*CJx+xx*CJxx)
#endif
gont = gont/2.d0/xx/(1.d0-xx)
return
end function Theta_rhs
!///////////////////////////////////////////////////////////////////////////////////////////////////////////////
subroutine fake_Theta_rhs(lx,X,rhs,Theta)
implicit none
integer,intent(in) :: lx
double complex,dimension(lx),intent(in) :: Theta
double complex,dimension(lx),intent(out) :: rhs
real*8,dimension(lx),intent(in) :: X
call cderivs_x(lx,X,Theta,rhs)
return
end subroutine fake_Theta_rhs
!///////////////////////////////////////////////////////////////////////////////////////////////////////////////
! try other guy's old method
function Theta_rhs_o(xx,Rmin,beta,betax,KK,KKx,CU,CUx,DCUx,bDCU,bDCUx,DCU,CB,DCB,W,Wx,CJ,DCJ,CJx,CJxx, &
DCJx,bDCB,Cnu,Cnux,Ck,Theta) result(gont)
implicit none
double complex,intent(in) :: CU,CUx,DCUx,bDCU,bDCUx,DCU,CB,DCB,CJ,DCJ,CJx,CJxx,DCJx
double complex,intent(in) :: Cnu,Cnux,Ck,bDCB,Theta
real*8,intent(in) :: xx,Rmin,beta,betax,KK,KKx,W,Wx
double complex :: JH,II,gont
real*8 :: V,Vx,Pu
II = dcmplx(0.d0,1.d0)
V = xx*Rmin/(1.d0-xx)*(1.d0+xx*Rmin/(1.d0-xx)*W)
Vx = Rmin/(1.d0-xx)**2+2.d0*xx*Rmin*Rmin/(1.d0-xx)**3*W+xx*xx*Rmin*Rmin/(1.d0-xx)**2*Wx
Pu = 2.d0*xx*(1.d0-xx)/KK*dreal(Theta*(dconjg(CJx)*KK-dconjg(CJ)*KKx))
JH = (1.d0-xx)/xx/Rmin*dexp(2.d0*beta)*(-KK*DCJ*dconjg(CB)+ &
(KK*Cnu+(KK*KK-1.d0)*DCJ-2.d0*KK*Ck)*CB+CJ* &
((2.d0*Ck-Cnu)*dconjg(CB)-2.d0*KK*(bDCB+CB*dconjg(CB))+ &
2.d0*dreal((Cnu-Ck)*dconjg(CB)+dconjg(CJ)*(DCB+CB*CB)))) &
+0.5d0*Rmin*xx**3*(1.d0-xx)*dexp(-2.d0*beta)* &
((KK*CUx+CJ*dconjg(CUx))**2- &
CJ*dreal(dconjg(CUx)*(KK*CUx+CJ*dconjg(CUx)))) &
-0.5d0*(Cnu*(xx*(1.d0-xx)*CUx+2.d0*CU)+DCJ*(xx*(1.d0-xx)*dconjg(CUx)+ &
2.d0*dconjg(CU)))+CJ*II*dimag(xx*(1.d0-xx)*bDCUx+2.d0*bDCU) &
-xx*(1.d0-xx)*CJx*dreal(bDCU) &
+xx*(1.d0-xx)*(dconjg(CU)*DCJ+CU*Cnu)*II*dimag(CJ*dconjg(CJx)) &
-xx*(1.d0-xx)*(dconjg(CU)*DCJx+CU*Cnux) &
-2.d0*xx*(1.d0-xx)*(CJ*KKx-KK*CJx)*(dreal(dconjg(CU)*Ck)+ &
II*dimag(KK*bDCU-dconjg(CJ)*DCU)) &
-8.d0*CJ*((1.d0-xx)**2/Rmin+xx*(1.d0-xx)*W)*betax
gont = -KK*(xx*(1-xx)*DCUx+2.d0*DCU)+2.d0*(1.d0-xx)/xx/Rmin*dexp(2.d0*beta)*(DCB+CB*CB) &
-(xx*(1.d0-xx)*Wx+W)*CJ+JH+CJ*Pu &
-(1.d0-xx)*(1.d0-xx)/xx/xx/Rmin/Rmin*V*(CJ+xx*(1.d0-xx)*CJx) &
+(1.d0-xx)*(1.d0-xx)*(1.d0-xx)/xx/Rmin/Rmin*Vx*(CJ+xx*(1.d0-xx)*CJx) &
+(1.d0-xx)**4/xx/Rmin/Rmin*V*(2.d0*CJx+xx*CJxx)
return
end function Theta_rhs_o
#if (RKorAM == 0)
!--------------------------------------------------------------------
! this R is indeed x
function NullEvol_Theta_o(ex,crho,sigma,R,RJ,IJ,RU,IU,beta,RB,IB, &
Rnu,Inu,Rk,Ik,RTheta,ITheta,W,Rmin, &
qlR1,qlR2,qlI1,qlI2,quR1,quR2,quI1,quI2,gR,gI) result(gont)
implicit none
integer,intent(in ):: ex(1:3)
real*8,intent(in),dimension(ex(1))::crho
real*8,intent(in),dimension(ex(2))::sigma
real*8,intent(in),dimension(ex(3))::R
real*8, dimension(ex(1),ex(2),ex(3)),intent(in) :: beta,W
real*8, dimension(ex(1),ex(2),ex(3)),intent(in) :: RJ,IJ,RU,IU,RB,IB,Rnu,Inu,Rk,Ik
real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: RTheta,ITheta
real*8,intent(in) :: Rmin
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: qlR1,qlR2,qlI1,qlI2
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: quR1,quR2,quI1,quI2
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: gR,gI
! gont = 0: success; gont = 1: something wrong
integer::gont
double complex,dimension(ex(1),ex(2),ex(3)) :: CU,DCU,bDCU,CB,DCB,bDCB,CJ,DCJ
double complex :: CTheta0,CTheta,CTheta1,RHS
integer :: i,j,k,RK4
double complex,dimension(ex(3)) :: Cnu,Ck,HCnu,HCk,HCU,HDCU,CUx,HCUx,DCUx,HDCUx,HbDCU,bDCUx,HbDCUx
double complex,dimension(ex(3)) :: Cnux,HCnux,HCJ,HDCJ,CJx,HCJx,CJxx,HCJxx,DCJx,HDCJx,HCB,HDCB,HbDCB
real*8,dimension(ex(3)) :: KK,KKx,HKK,HKKx,Hbeta,betax,Hbetax,HW,Wx,HWx
double complex :: Theta_rhs_o
real*8 :: dR
!!! sanity check
dR = sum(RJ)+sum(IJ)+sum(RU)+sum(IU)+sum(beta)+sum(RB)+sum(IB) + &
sum(Rnu)+sum(Inu)+sum(Rk)+sum(Ik)+sum(RTheta)+sum(ITheta)
if(dR.ne.dR) then
if(sum(RJ).ne.sum(RJ))write(*,*)"NullEvol_Theta: find NaN in RJ"
if(sum(IJ).ne.sum(IJ))write(*,*)"NullEvol_Theta: find NaN in IJ"
if(sum(RU).ne.sum(RU))write(*,*)"NullEvol_Theta: find NaN in RU"
if(sum(IU).ne.sum(IU))write(*,*)"NullEvol_Theta: find NaN in IU"
if(sum(beta).ne.sum(beta))write(*,*)"NullEvol_Theta: find NaN in beta"
if(sum(RB).ne.sum(RB))write(*,*)"NullEvol_Theta: find NaN in RB"
if(sum(IB).ne.sum(IB))write(*,*)"NullEvol_Theta: find NaN in IB"
if(sum(Rnu).ne.sum(Rnu))write(*,*)"NullEvol_Theta: find NaN in Rnu"
if(sum(Inu).ne.sum(Inu))write(*,*)"NullEvol_Theta: find NaN in Inu"
if(sum(Rk).ne.sum(Rk))write(*,*)"NullEvol_Theta: find NaN in Rk"
if(sum(Ik).ne.sum(Ik))write(*,*)"NullEvol_Theta: find NaN in Ik"
if(sum(RTheta).ne.sum(RTheta))write(*,*)"NullEvol_Theta: find NaN in RTheta"
if(sum(ITheta).ne.sum(ITheta))write(*,*)"NullEvol_Theta: find NaN in ITheta"
gont = 1
return
endif
dR = R(2) - R(1)
CU = dcmplx(RU,IU)
CB = dcmplx(RB,IB)
CJ = dcmplx(RJ,IJ)
do k=1,ex(3)
call derivs_eth(ex(1:2),crho,sigma,CU(:,:,k),DCU(:,:,k),1,1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
call derivs_eth(ex(1:2),crho,sigma,CU(:,:,k),bDCU(:,:,k),1,-1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
call derivs_eth(ex(1:2),crho,sigma,CB(:,:,k),DCB(:,:,k),1,1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
call derivs_eth(ex(1:2),crho,sigma,CB(:,:,k),bDCB(:,:,k),1,-1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
call derivs_eth(ex(1:2),crho,sigma,CJ(:,:,k),DCJ(:,:,k),2,1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
enddo
do j=1,ex(2)
do i=1,ex(1)
CTheta0 = dcmplx(RTheta(i,j,1),ITheta(i,j,1))
Cnu = dcmplx(Rnu(i,j,:),Inu(i,j,:))
Ck = dcmplx(Rk(i,j,:),Ik(i,j,:))
call cget_half_x(ex(3),CB(i,j,:),HCB)
call cget_half_x(ex(3),DCB(i,j,:),HDCB)
call cget_half_x(ex(3),bDCB(i,j,:),HbDCB)
call cget_half_x(ex(3),Cnu,HCnu)
call cderivs_x(ex(3),R,Cnu,Cnux)
call cget_half_x(ex(3),Cnux,HCnux)
call cget_half_x(ex(3),Ck,HCk)
call rget_half_x(ex(3),beta(i,j,:),Hbeta)
call rderivs_x(ex(3),R,beta(i,j,:),betax)
call rget_half_x(ex(3),betax,Hbetax)
KK = dsqrt(1.d0+RJ(i,j,:)*RJ(i,j,:)+IJ(i,j,:)*IJ(i,j,:))
call rget_half_x(ex(3),KK,HKK)
call rderivs_x(ex(3),R,KK,KKx)
call rget_half_x(ex(3),KKx,HKKx)
call rderivs_x(ex(3),R,W,Wx)
call rget_half_x(ex(3),Wx,HWx)
call rget_half_x(ex(3),W(i,j,:),HW)
call cget_half_x(ex(3),CU(i,j,:),HCU)
call cderivs_x(ex(3),R,DCU(i,j,:),DCUx)
call cderivs_x(ex(3),R,CU(i,j,:),CUx)
call cget_half_x(ex(3),DCUx,HDCUx)
call cget_half_x(ex(3),CUx,HCUx)
call cget_half_x(ex(3),DCU(i,j,:),HDCU)
call cderivs_x(ex(3),R,bDCU(i,j,:),bDCUx)
call cget_half_x(ex(3),bDCUx,HbDCUx)
call cget_half_x(ex(3),bDCU(i,j,:),HbDCU)
call cderivs_x(ex(3),R,CJ(i,j,:),CJx)
call cdderivs_x(ex(3),R,CJ(i,j,:),CJxx)
call cget_half_x(ex(3),CJx,HCJx)
call cget_half_x(ex(3),CJxx,HCJxx)
call cderivs_x(ex(3),R,DCJ(i,j,:),DCJx)
call cget_half_x(ex(3),DCJx,HDCJx)
do k=1,ex(3)-1
RHS = Theta_rhs_o(R(k)+dR/2.d0,Rmin,Hbeta(k),betax(k),HKK(k),KKx(k),HCU(k),CUx(k),DCUx(k),HbDCU(k),bDCUx(k), &
HDCU(k),HCB(k),HDCB(k),HW(k),Wx(k),HCJ(k),HDCJ(k), &
CJx(k),CJxx(k),DCJx(k),HbDCB(k),HCnu(k),Cnux(k),HCk(k),CTheta0)
CTheta1 = RHS-(1-2.d0*(R(k)+dR/2.d0)*(1.d0-R(k)-dR/2.d0)/dR)*CTheta0
CTheta1 = CTheta1/(1+2.d0*(R(k)+dR/2.d0)*(1.d0-R(k)-dR/2.d0)/dR)
CTheta0 = CTheta1
RTheta(i,j,k+1) = dreal(CTheta0)
ITheta(i,j,k+1) = dimag(CTheta0)
enddo
enddo
enddo
gont = 0
return
end function NullEvol_Theta_o
!------------------------------------------------------------------------------
! this R is indeed x
function NullEvol_beta(ex,crho,sigma,R,RJ,IJ,beta,KKx,HKKx) result(gont)
implicit none
integer,intent(in ):: ex(1:3)
real*8,intent(in),dimension(ex(1))::crho
real*8,intent(in),dimension(ex(2))::sigma
real*8,intent(in),dimension(ex(3))::R
real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: beta
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: RJ,IJ,KKx,HKKx
! gont = 0: success; gont = 1: something wrong
integer::gont
real*8 :: dR
double complex, dimension(ex(3)):: CJ,CJx,HCJx
real*8 :: betah0,betah1,betah,rhs
integer :: i,j,k,RK4
real*8 :: beta_rhs
!!! sanity check
dR = sum(RJ)+sum(IJ)+sum(beta)+sum(KKx)+sum(HKKx)
if(dR.ne.dR) then
if(sum(RJ).ne.sum(RJ))write(*,*)"NullEvol_beta: find NaN in RJ"
if(sum(IJ).ne.sum(IJ))write(*,*)"NullEvol_beta: find NaN in IJ"
if(sum(beta).ne.sum(beta))write(*,*)"NullEvol_beta: find NaN in beta"
if(sum(KKx).ne.sum(KKx))write(*,*)"NullEvol_beta: find NaN in KKx"
if(sum(HKKx).ne.sum(HKKx))write(*,*)"NullEvol_beta: find NaN in HKKx"
gont = 1
return
endif
dR = R(2) - R(1)
do j=1,ex(2)
do i=1,ex(1)
betah0 = beta(i,j,1)
CJ = dcmplx(RJ(i,j,:),IJ(i,j,:))
call cderivs_x(ex(3),R,CJ,CJx)
call cget_half_x(ex(3),CJx,HCJx)
#ifdef OLD
do k = 1,ex(3)-1
! note our CJx(ex(3)) = (CJ(ex(3))-CJ(ex(3)-1))/dR
! note our KKx(ex(3)) = (KK(ex(3))-KK(ex(3)-1))/dR
rhs = beta_rhs(R(k)+dR/2.d0,CJx(k+1),KKx(i,j,k+1))
beta(i,j,k+1) = beta(i,j,k) + rhs*dR
enddo
#else
do k=1,ex(3)-1
RK4 = 0
rhs = beta_rhs(R(k),CJx(k),KKx(i,j,k))
call rungekutta4_scalar(dR,betah0,betah,rhs,RK4)
RK4 = 1
betah1 = beta_rhs(R(k)+dR/2.d0,HCJx(k),HKKx(i,j,k))
call rungekutta4_scalar(dR,betah0,betah1,rhs,RK4)
call rswap(betah,betah1)
RK4 = 2
betah1 = beta_rhs(R(k)+dR/2.d0,HCJx(k),HKKx(i,j,k))
call rungekutta4_scalar(dR,betah0,betah1,rhs,RK4)
call rswap(betah,betah1)
RK4 = 3
betah1 = beta_rhs(R(k+1),CJx(k+1),KKx(i,j,k+1))
call rungekutta4_scalar(dR,betah0,betah1,rhs,RK4)
call rswap(betah0,betah1)
beta(i,j,k+1) = betah0
enddo
! above k takes ex(3)-1 then do not need this closing step
#if 1
! closing step
k = ex(3)-1
! note our CJx(ex(3)) = (CJ(ex(3))-CJ(ex(3)-1))/dR
! note our KKx(ex(3)) = (KK(ex(3))-KK(ex(3)-1))/dR
rhs = beta_rhs(R(k)+dR/2.d0,CJx(k+1),KKx(i,j,k+1))
beta(i,j,k+1) = beta(i,j,k) + rhs*dR
#endif
#endif
enddo
enddo
gont = 0
return
end function NullEvol_beta
!------------------------------------------------------------------------------
! this R is indeed x
function NullEvol_Q(ex,crho,sigma,R,RJ,IJ,Rk,Ik,Rnu,Inu,RB,IB,RQ,IQ,KK,Hkk,KKx,HKKx, &
qlR1,qlR2,qlI1,qlI2,quR1,quR2,quI1,quI2,gR,gI) result(gont)
implicit none
integer,intent(in ):: ex(1:3)
real*8,intent(in),dimension(ex(1))::crho
real*8,intent(in),dimension(ex(2))::sigma
real*8,intent(in),dimension(ex(3))::R
real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: RQ,IQ
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: RJ,IJ,KK,Hkk,KKx,HKKx
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: Rk,Ik,Rnu,Inu,RB,IB
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: qlR1,qlR2,qlI1,qlI2
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: quR1,quR2,quI1,quI2
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: gR,gI
! gont = 0: success; gont = 1: something wrong
integer::gont
real*8 :: xx,dR
double complex :: CQ0,CQ,CQ1,RHS
double complex,dimension(ex(3)) :: CJx,HCJx,DCJx,HDCJx,Ck,Ckx,HCkx,Cnu,Cnux,HCnux,CB,CBx,HCBx
double complex,dimension(ex(3)) :: HCJ,HCk,HCnu,HCB,HDCJ
double complex, dimension(ex(1),ex(2),ex(3)) :: CJ,DCJ
integer :: i,j,k,RK4
double complex :: Q_rhs
!!! sanity check
dR = sum(RJ)+sum(IJ)+sum(RQ)+sum(IQ) &
+sum(RK)+sum(IK)+sum(Rnu)+sum(Inu)+sum(RB)+sum(IB) &
+sum(KK)+sum(HKK)+sum(KKx)+sum(HKKx)
if(dR.ne.dR) then
if(sum(RJ).ne.sum(RJ))write(*,*)"NullEvol_Q: find NaN in RJ"
if(sum(IJ).ne.sum(IJ))write(*,*)"NullEvol_Q: find NaN in IJ"
if(sum(RQ).ne.sum(RQ))write(*,*)"NullEvol_Q: find NaN in RQ"
if(sum(IQ).ne.sum(IQ))write(*,*)"NullEvol_Q: find NaN in IQ"
if(sum(RK).ne.sum(RK))write(*,*)"NullEvol_Q: find NaN in RK"
if(sum(IK).ne.sum(IK))write(*,*)"NullEvol_Q: find NaN in IK"
if(sum(Rnu).ne.sum(Rnu))write(*,*)"NullEvol_Q: find NaN in Rnu"
if(sum(Inu).ne.sum(Inu))write(*,*)"NullEvol_Q: find NaN in Inu"
if(sum(RB).ne.sum(RB))write(*,*)"NullEvol_Q: find NaN in RB"
if(sum(IB).ne.sum(IB))write(*,*)"NullEvol_Q: find NaN in IB"
if(sum(KK).ne.sum(KK))write(*,*)"NullEvol_Q: find NaN in KK"
if(sum(HKK).ne.sum(HKK))write(*,*)"NullEvol_Q: find NaN in HKK"
if(sum(KKx).ne.sum(KKx))write(*,*)"NullEvol_Q: find NaN in KKx"
if(sum(HKKx).ne.sum(HKKx))write(*,*)"NullEvol_Q: find NaN in HKKx"
gont = 1
return
endif
dR = R(2) - R(1)
CJ = dcmplx(RJ,IJ)
do k=1,ex(3)
call derivs_eth(ex(1:2),crho,sigma,CJ(:,:,k),DCJ(:,:,k),2,1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
enddo
do j=1,ex(2)
do i=1,ex(1)
CQ0 = dcmplx(RQ(i,j,1),IQ(i,j,1))
call cderivs_x(ex(3),R,CJ(i,j,:),CJx)
call cget_half_x(ex(3),CJx,HCJx)
call cget_half_x(ex(3),CJ,HCJ)
call cderivs_x(ex(3),R,DCJ(i,j,:),DCJx)
call cget_half_x(ex(3),DCJx,HDCJx)
call cget_half_x(ex(3),DCJ,HDCJ)
Ck = dcmplx(Rk(i,j,:),Ik(i,j,:))
call cderivs_x(ex(3),R,Ck,Ckx)
call cget_half_x(ex(3),Ckx,HCkx)
call cget_half_x(ex(3),Ck,HCk)
Cnu = dcmplx(Rnu(i,j,:),Inu(i,j,:))
call cderivs_x(ex(3),R,Cnu,Cnux)
call cget_half_x(ex(3),Cnux,HCnux)
call cget_half_x(ex(3),Cnu,HCnu)
CB = dcmplx(RB(i,j,:),IB(i,j,:))
call cderivs_x(ex(3),R,CB,CBx)
call cget_half_x(ex(3),CBx,HCBx)
call cget_half_x(ex(3),CB,HCB)
#ifdef OLD
do k = 1,ex(3)-1
xx = R(k)+dR/2.d0
! note our HCJ(ex(3)-1) = (CJ(ex(3))+CJ(ex(3)-1))/2
! note our KKx(ex(3)) = (KK(ex(3))-KK(ex(3)-1))/dR
RHS = Q_rhs(xx,HCJ(k),CJx(k+1),DCJx(k+1),HKK(i,j,k),HCk(k),Ckx(k+1),Cnux(k+1),KKx(i,j,k+1),CBx(k+1),HCnu(k),HDCJ(k),HCB(k),0)
RHS = RHS+CQ0*(1.d0/dR-1.d0/xx/(1.d0-xx))
CQ0 = RHS/(1.d0/dR+1.d0/xx/(1.d0-xx))
RQ(i,j,k+1) = dreal(CQ0)
IQ(i,j,k+1) = dimag(CQ0)
enddo
#else
do k=1,ex(3)-2
RK4 = 0
RHS = Q_rhs(R(k),CJ(i,j,k),CJx(k),DCJx(k),KK(i,j,k),Ck(k),Ckx(k),Cnux(k),KKx(i,j,k),CBx(k),Cnu(k),DCJ(i,j,k),CB(k),CQ0)
call rungekutta4_cplxscalar(dR,CQ0,CQ,RHS,RK4)
RK4 = 1
CQ1 = Q_rhs(R(k)+dR/2.d0,HCJ(k),HCJx(k),HDCJx(k),HKK(i,j,k),HCk(k),HCkx(k),HCnux(k),HKKx(i,j,k), &
HCBx(k),HCnu(k),HDCJ(k),HCB(k),CQ)
call rungekutta4_cplxscalar(dR,CQ0,CQ1,RHS,RK4)
call cswap(CQ,CQ1)
RK4 = 2
CQ1 = Q_rhs(R(k)+dR/2.d0,HCJ(k),HCJx(k),HDCJx(k),HKK(i,j,k),HCk(k),HCkx(k),HCnux(k),HKKx(i,j,k), &
HCBx(k),HCnu(k),HDCJ(k),HCB(k),CQ)
call rungekutta4_cplxscalar(dR,CQ0,CQ1,RHS,RK4)
call cswap(CQ,CQ1)
RK4 = 3
CQ1 = Q_rhs(R(k+1),CJ(i,j,k+1),CJx(k+1),DCJx(k+1),KK(i,j,k+1),Ck(k+1),Ckx(k+1),Cnux(k+1),KKx(i,j,k+1), &
CBx(k+1),Cnu(k+1),DCJ(i,j,k+1),CB(k+1),CQ)
call rungekutta4_cplxscalar(dR,CQ0,CQ1,RHS,RK4)
call cswap(CQ0,CQ1)
RQ(i,j,k+1) = dreal(CQ0)
IQ(i,j,k+1) = dimag(CQ0)
enddo
#if 0
k = ex(3)
CQ0 = -2*CB(k)
RQ(i,j,k+1) = dreal(CQ0)
IQ(i,j,k+1) = dimag(CQ0)
#else
! closing step
k = ex(3)-1
CQ0 = dcmplx(RQ(i,j,k),IQ(i,j,k))
xx = R(k)+dR/2.d0
! note our HCJ(ex(3)-1) = (CJ(ex(3))+CJ(ex(3)-1))/2
! note our KKx(ex(3)) = (KK(ex(3))-KK(ex(3)-1))/dR
RHS = Q_rhs(xx,HCJ(k),CJx(k+1),DCJx(k+1),HKK(i,j,k), &
HCk(k),Ckx(k+1),Cnux(k+1),KKx(i,j,k+1),CBx(k+1),HCnu(k),HDCJ(k),HCB(k),dcmplx(0.d0,0.d0))
RHS = RHS+CQ0*(1.d0/dR-1.d0/xx/(1.d0-xx))
CQ0 = RHS/(1.d0/dR+1.d0/xx/(1.d0-xx))
RQ(i,j,k+1) = dreal(CQ0)
IQ(i,j,k+1) = dimag(CQ0)
#endif
#endif
enddo
enddo
gont = 0
return
end function NullEvol_Q
!--------------------------------------------------------------------
! this R is indeed x
function NullEvol_U(ex,crho,sigma,R,RJ,IJ,RQ,IQ,KK,HKK,beta,RU,IU, &
Rmin) result(gont)
implicit none
integer,intent(in ):: ex(1:3)
real*8,intent(in),dimension(ex(1))::crho
real*8,intent(in),dimension(ex(2))::sigma
real*8,intent(in),dimension(ex(3))::R
real*8, dimension(ex(1),ex(2),ex(3)),intent(in) :: RJ,IJ,RQ,IQ,beta,KK,HKK
real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: RU,IU
real*8,intent(in) :: Rmin
! gont = 0: success; gont = 1: something wrong
integer::gont
real*8 :: dR
double complex :: CU0,CU,CU1,RHS
integer :: i,j,k,RK4
double complex,dimension(ex(3)) :: CJ,CQ,HCJ,HCQ
real*8,dimension(ex(3)) :: Hbeta
double complex :: U_rhs
!!! sanity check
dR = sum(RJ)+sum(IJ)+sum(RQ)+sum(IQ)+sum(beta)+sum(RU)+sum(IU)+sum(KK)+sum(HKK)
if(dR.ne.dR) then
if(sum(RJ).ne.sum(RJ))write(*,*)"NullEvol_U: find NaN in RJ"
if(sum(IJ).ne.sum(IJ))write(*,*)"NullEvol_U: find NaN in IJ"
if(sum(RQ).ne.sum(RQ))write(*,*)"NullEvol_U: find NaN in RQ"
if(sum(IQ).ne.sum(IQ))write(*,*)"NullEvol_U: find NaN in IQ"
if(sum(beta).ne.sum(beta))write(*,*)"NullEvol_U: find NaN in beta"
if(sum(RU).ne.sum(RU))write(*,*)"NullEvol_U: find NaN in RU"
if(sum(IU).ne.sum(IU))write(*,*)"NullEvol_U: find NaN in IU"
if(sum(KK).ne.sum(KK))write(*,*)"NullEvol_U: find NaN in KK"
if(sum(HKK).ne.sum(HKK))write(*,*)"NullEvol_U: find NaN in HKK"
gont = 1
return
endif
dR = R(2) - R(1)
do j=1,ex(2)
do i=1,ex(1)
CU0 = dcmplx(RU(i,j,1),IU(i,j,1))
CJ = dcmplx(RJ(i,j,:),IJ(i,j,:))
CQ = dcmplx(RQ(i,j,:),IQ(i,j,:))
call cget_half_x(ex(3),CJ,HCJ)
call cget_half_x(ex(3),CQ,HCQ)
call rget_half_x(ex(3),beta(i,j,:),Hbeta)
#ifdef OLD
do k = 1,ex(3)-1
! note our HCJ(ex(3)-1) = (CJ(ex(3))+CJ(ex(3)-1))/2
! note our KKx(ex(3)) = (KK(ex(3))-KK(ex(3)-1))/dR
RHS = U_rhs(R(k)+dR/2,Rmin,Hbeta(k),HKK(i,j,k),HCQ(k),HCJ(k))
CU0 = CU0+RHS*dR
RU(i,j,k+1) = dreal(CU0)
IU(i,j,k+1) = dimag(CU0)
enddo
#else
do k=1,ex(3)-2
RK4 = 0
RHS = U_rhs(R(k),Rmin,beta(i,j,k),KK(i,j,k),CQ(k),CJ(k))
call rungekutta4_cplxscalar(dR,CU0,CU,RHS,RK4)
RK4 = 1
CU1 = U_rhs(R(k)+dR/2.d0,Rmin,Hbeta(k),HKK(i,j,k),HCQ(k),HCJ(k))
call rungekutta4_cplxscalar(dR,CU0,CU1,RHS,RK4)
call cswap(CU,CU1)
RK4 = 2
CU1 = U_rhs(R(k)+dR/2.d0,Rmin,Hbeta(k),HKK(i,j,k),HCQ(k),HCJ(k))
call rungekutta4_cplxscalar(dR,CU0,CU1,RHS,RK4)
call cswap(CU,CU1)
RK4 = 3
CU1 = U_rhs(R(k+1),Rmin,beta(i,j,k+1),KK(i,j,k+1),CQ(k+1),CJ(k+1))
call rungekutta4_cplxscalar(dR,CU0,CU1,RHS,RK4)
call cswap(CU0,CU1)
RU(i,j,k+1) = dreal(CU0)
IU(i,j,k+1) = dimag(CU0)
enddo
! above k takes ex(3)-1 then do not need closing step
#if 1
! closing step
k = ex(3)-1
CU0 = dcmplx(RU(i,j,k),IU(i,j,k))
! note our HCJ(ex(3)-1) = (CJ(ex(3))+CJ(ex(3)-1))/2
! note our KKx(ex(3)) = (KK(ex(3))-KK(ex(3)-1))/dR
RHS = U_rhs(R(k)+dR/2,Rmin,Hbeta(k),HKK(i,j,k),HCQ(k),HCJ(k))
CU0 = CU0+RHS*dR
RU(i,j,k+1) = dreal(CU0)
IU(i,j,k+1) = dimag(CU0)
#endif
#endif
enddo
enddo
gont = 0
return
end function NullEvol_U
!----------------------------------------------------------------------------------------
! this R is indeed x
function NullEvol_W(ex,crho,sigma,R,RJ,IJ,RB,IB,Rnu,Inu,Rk,Ik, &
RU,IU,RQ,IQ,W,beta,KK,HKK,Rmin, &
qlR1,qlR2,qlI1,qlI2,quR1,quR2,quI1,quI2,gR,gI) result(gont)
implicit none
integer,intent(in ):: ex(1:3)
real*8,intent(in),dimension(ex(1))::crho
real*8,intent(in),dimension(ex(2))::sigma
real*8,intent(in),dimension(ex(3))::R
real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: W
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: RJ,IJ,RB,IB
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: Rnu,Inu,Rk,Ik
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: RU,IU,RQ,IQ,beta,KK,HKK
real*8,intent(in ) :: Rmin
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: qlR1,qlR2,qlI1,qlI2
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: quR1,quR2,quI1,quI2
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: gR,gI
! gont = 0: success; gont = 1: something wrong
integer::gont
real*8 :: dR
real*8, dimension(ex(3)) :: Hbeta
double complex, dimension(ex(1),ex(2),ex(3)) :: CU,DCU,bDCU
double complex, dimension(ex(1),ex(2),ex(3)) :: CB,DCB,bDCB,CJ,DCJ,Cnu,bDCnu,Ck,bDCk
double complex, dimension(ex(3)) :: HCB,HDCB,HbDCB,HCJ,HDCJ,HCnu,HbDCnu,HCk,HbDCk
double complex, dimension(ex(3)) :: HbDCU,bDCUx,HbDCUx,CQ,HCQ
real*8 :: Wh0,Wh1,Wh,rhs
integer :: i,j,k,RK4
real*8 :: xx,W_rhs
!!! sanity check
dR = sum(RJ)+sum(IJ)+sum(beta)+sum(RB)+sum(IB)+sum(Rnu)+sum(Inu) &
+sum(Rk)+sum(Ik)+sum(W)+sum(RU)+sum(IU)+sum(RQ)+sum(IQ)&
+sum(KK)+sum(HKK)
if(dR.ne.dR) then
if(sum(RJ).ne.sum(RJ))write(*,*)"NullEvol_W: find NaN in RJ"
if(sum(IJ).ne.sum(IJ))write(*,*)"NullEvol_W: find NaN in IJ"
if(sum(beta).ne.sum(beta))write(*,*)"NullEvol_W: find NaN in beta"
if(sum(RB).ne.sum(RB))write(*,*)"NullEvol_W: find NaN in RB"
if(sum(IB).ne.sum(IB))write(*,*)"NullEvol_W: find NaN in IB"
if(sum(Rnu).ne.sum(Rnu))write(*,*)"NullEvol_W: find NaN in Rnu"
if(sum(Inu).ne.sum(Inu))write(*,*)"NullEvol_W: find NaN in Inu"
if(sum(Rk).ne.sum(Rk))write(*,*)"NullEvol_W: find NaN in Rk"
if(sum(Ik).ne.sum(Ik))write(*,*)"NullEvol_W: find NaN in Ik"
if(sum(RU).ne.sum(RU))write(*,*)"NullEvol_W: find NaN in RU"
if(sum(IU).ne.sum(IU))write(*,*)"NullEvol_W: find NaN in IU"
if(sum(RQ).ne.sum(RQ))write(*,*)"NullEvol_W: find NaN in RQ"
if(sum(IQ).ne.sum(IQ))write(*,*)"NullEvol_W: find NaN in IQ"
if(sum(W).ne.sum(W))write(*,*)"NullEvol_W: find NaN in W"
if(sum(KK).ne.sum(KK))write(*,*)"NullEvol_W: find NaN in KK"
if(sum(HKK).ne.sum(HKK))write(*,*)"NullEvol_W: find NaN in HKK"
gont = 1
return
endif
dR = R(2) - R(1)
CB = dcmplx(RB,IB)
CU = dcmplx(RU,IU)
Ck = dcmplx(Rk,Ik)
Cnu = dcmplx(Rnu,Inu)
CJ = dcmplx(RJ,IJ)
do k=1,ex(3)
call derivs_eth(ex(1:2),crho,sigma,CJ(:,:,k),DCJ(:,:,k),2,1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
call derivs_eth(ex(1:2),crho,sigma,CB(:,:,k),DCB(:,:,k),1,1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
call derivs_eth(ex(1:2),crho,sigma,CB(:,:,k),bDCB(:,:,k),1,-1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
call derivs_eth(ex(1:2),crho,sigma,CU(:,:,k),bDCU(:,:,k),1,-1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
call derivs_eth(ex(1:2),crho,sigma,Ck(:,:,k),bDCk(:,:,k),1,-1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
call derivs_eth(ex(1:2),crho,sigma,Cnu(:,:,k),bDCnu(:,:,k),1,-1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
enddo
do j=1,ex(2)
do i=1,ex(1)
Wh0 = W(i,j,1)
call rget_half_x(ex(3),beta(i,j,:),Hbeta)
call cderivs_x(ex(3),R,bDCU,bDCUx)
call cget_half_x(ex(3),bDCUx,HbDCUx)
call cget_half_x(ex(3),bDCU,HbDCU)
call cget_half_x(ex(3),DCJ,HDCJ)
call cget_half_x(ex(3),DCB,HDCB)
call cget_half_x(ex(3),bDCB,HbDCB)
call cget_half_x(ex(3),CB,HCB)
call cget_half_x(ex(3),CJ,HCJ)
call cget_half_x(ex(3),Cnu,HCnu)
call cget_half_x(ex(3),Ck,HCk)
CQ = dcmplx(RQ(i,j,:),IQ(i,j,:))
call cget_half_x(ex(3),CQ,HCQ)
call cget_half_x(ex(3),bDCk,HbDCk)
call cget_half_x(ex(3),bDCnu,HbDCnu)
#ifdef OLD
do k = 1,ex(3)-1
xx = R(k)+dR/2
! note our HCJ(ex(3)-1) = (CJ(ex(3))+CJ(ex(3)-1))/2
! note our KKx(ex(3)) = (KK(ex(3))-KK(ex(3)-1))/dR
rhs = W_rhs(xx,Rmin,Hbeta(k),HKK(i,j,k),HDCB(k),HCB(k),HCJ(k),HCnu(k),HCk(k),0, &
HCQ(k),HbDCk(k),HbDCnu(k),HbDCB(k),HbDCU(k),bDCUx(k+1),HDCJ(k))
rhs = rhs+Wh0*(1.d0/dR-1.d0/xx/(1.d0-xx))
W(i,j,k+1) = rhs/(1.d0/dR+1.d0/xx/(1.d0-xx))
enddo
#else
do k=1,ex(3)-2
RK4 = 0
rhs = W_rhs(R(k),Rmin,beta(i,j,k),KK(i,j,k),DCB(i,j,k),CB(i,j,k),CJ(i,j,k),Cnu(i,j,k),Ck(i,j,k),Wh0, &
CQ(k),bDCk(i,j,k),bDCnu(i,j,k),bDCB(i,j,k),bDCU(i,j,k),bDCUx(k),DCJ(i,j,k))
call rungekutta4_scalar(dR,Wh0,Wh,rhs,RK4)
RK4 = 1
Wh1 = W_rhs(R(k)+dR/2.d0,Rmin,Hbeta(k),HKK(i,j,k),HDCB(k),HCB(k),HCJ(k),HCnu(k),HCk(k),Wh, &
HCQ(k),HbDCk(k),HbDCnu(k),HbDCB(k),HbDCU(k),HbDCUx(k),HDCJ(k))
call rungekutta4_scalar(dR,Wh0,Wh1,rhs,RK4)
call rswap(Wh,Wh1)
RK4 = 2
Wh1 = W_rhs(R(k)+dR/2.d0,Rmin,Hbeta(k),HKK(i,j,k),HDCB(k),HCB(k),HCJ(k),HCnu(k),HCk(k),Wh, &
HCQ(k),HbDCk(k),HbDCnu(k),HbDCB(k),HbDCU(k),HbDCUx(k),HDCJ(k))
call rungekutta4_scalar(dR,Wh0,Wh1,rhs,RK4)
call rswap(Wh,Wh1)
RK4 = 3
Wh1 = W_rhs(R(k+1),Rmin,beta(i,j,k+1),KK(i,j,k+1),DCB(i,j,k+1),CB(i,j,k+1),CJ(i,j,k+1),Cnu(i,j,k+1),Ck(i,j,k+1),Wh, &
CQ(k+1),bDCk(i,j,k+1),bDCnu(i,j,k+1),bDCB(i,j,k+1),bDCU(i,j,k+1),bDCUx(k+1),DCJ(i,j,k+1))
call rungekutta4_scalar(dR,Wh0,Wh1,rhs,RK4)
call rswap(Wh0,Wh1)
W(i,j,k+1) = Wh0
enddo
#if 0
k = ex(3)
W(i,j,k) = dreal(bDCU(i,j,k))
#else
! closing step
k = ex(3)-1
Wh0 = W(i,j,k)
xx = R(k)+dR/2
! note our HCJ(ex(3)-1) = (CJ(ex(3))+CJ(ex(3)-1))/2
! note our KKx(ex(3)) = (KK(ex(3))-KK(ex(3)-1))/dR
rhs = W_rhs(xx,Rmin,Hbeta(k),HKK(i,j,k),HDCB(k),HCB(k),HCJ(k),HCnu(k),HCk(k),0.d0, &
HCQ(k),HbDCk(k),HbDCnu(k),HbDCB(k),HbDCU(k),bDCUx(k+1),HDCJ(k))
rhs = rhs+Wh0*(1.d0/dR-1.d0/xx/(1.d0-xx))
W(i,j,k+1) = rhs/(1.d0/dR+1.d0/xx/(1.d0-xx))
#endif
#endif
enddo
enddo
gont = 0
return
end function NullEvol_W
!-----------------------------------------------------------------------------------------------
! given exact Theta_x
! this R is indeed x
function NullEvol_Theta_givenx(ex,crho,sigma,R,RJ,IJ,RU,IU,beta,RB,IB, &
Rnu,Inu,Rk,Ik,RTheta,ITheta,W,KK,HKK,KKx,HKKx,Rmin, &
qlR1,qlR2,qlI1,qlI2,quR1,quR2,quI1,quI2,gR,gI, &
dquR1,dquR2,dquI1,dquI2, &
bdquR1,bdquR2,bdquI1,bdquI2, &
dgR,dgI,bdgR,bdgI,T,sst) result(gont)
implicit none
integer,intent(in ):: ex(1:3),sst
real*8,intent(in),dimension(ex(1))::crho
real*8,intent(in),dimension(ex(2))::sigma
real*8,intent(in),dimension(ex(3))::R
real*8, dimension(ex(1),ex(2),ex(3)),intent(in) :: beta,W,KK,HKK,KKx,HKKx
real*8, dimension(ex(1),ex(2),ex(3)),intent(in) :: RJ,IJ,RU,IU,RB,IB,Rnu,Inu,Rk,Ik
real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: RTheta,ITheta
real*8,intent(in) :: Rmin,T
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: qlR1,qlR2,qlI1,qlI2
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: quR1,quR2,quI1,quI2
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: 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
! gont = 0: success; gont = 1: something wrong
integer::gont
real*8,dimension(ex(3))::HR
real*8,dimension(ex(1),ex(2),ex(3)) :: RThetax,IThetax,HRThetax,HIThetax
double complex,dimension(ex(3)) :: fRHS,HfRHS
real*8 :: xx,dR
integer :: i,j,k,RK4
double complex :: CTheta0,CTheta,CTheta1,RHS
integer,parameter :: ks=1
!!! sanity check
dR = sum(RJ)+sum(IJ)+sum(RU)+sum(IU)+sum(beta)+sum(RB)+sum(IB) + &
sum(Rnu)+sum(Inu)+sum(Rk)+sum(Ik)+sum(RTheta)+sum(ITheta) + &
sum(KK)+sum(HKK)+sum(KKx)+sum(HKKx)
if(dR.ne.dR) then
if(sum(RJ).ne.sum(RJ))write(*,*)"NullEvol_Theta: find NaN in RJ"
if(sum(IJ).ne.sum(IJ))write(*,*)"NullEvol_Theta: find NaN in IJ"
if(sum(RU).ne.sum(RU))write(*,*)"NullEvol_Theta: find NaN in RU"
if(sum(IU).ne.sum(IU))write(*,*)"NullEvol_Theta: find NaN in IU"
if(sum(beta).ne.sum(beta))write(*,*)"NullEvol_Theta: find NaN in beta"
if(sum(RB).ne.sum(RB))write(*,*)"NullEvol_Theta: find NaN in RB"
if(sum(IB).ne.sum(IB))write(*,*)"NullEvol_Theta: find NaN in IB"
if(sum(Rnu).ne.sum(Rnu))write(*,*)"NullEvol_Theta: find NaN in Rnu"
if(sum(Inu).ne.sum(Inu))write(*,*)"NullEvol_Theta: find NaN in Inu"
if(sum(Rk).ne.sum(Rk))write(*,*)"NullEvol_Theta: find NaN in Rk"
if(sum(Ik).ne.sum(Ik))write(*,*)"NullEvol_Theta: find NaN in Ik"
if(sum(RTheta).ne.sum(RTheta))write(*,*)"NullEvol_Theta: find NaN in RTheta"
if(sum(ITheta).ne.sum(ITheta))write(*,*)"NullEvol_Theta: find NaN in ITheta"
if(sum(KK).ne.sum(KK))write(*,*)"NullEvol_Theta: find NaN in KK"
if(sum(HKK).ne.sum(HKK))write(*,*)"NullEvol_Theta: find NaN in HKK"
if(sum(KKx).ne.sum(KKx))write(*,*)"NullEvol_Theta: find NaN in KKx"
if(sum(HKKx).ne.sum(HKKx))write(*,*)"NullEvol_Theta: find NaN in HKKx"
gont = 1
return
endif
dR = R(2) - R(1)
HR = R+dR/2
call 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)
call get_exact_null_theta_x(ex,crho,sigma,HR,HRThetax,HIThetax,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)
do j=1,ex(2)
do i=1,ex(1)
CTheta0 = dcmplx(RTheta(i,j,ks),ITheta(i,j,ks))
fRHS = dcmplx(RThetax(i,j,:),IThetax(i,j,:))
HfRHS = dcmplx(HRThetax(i,j,:),HIThetax(i,j,:))
! call cget_half_x(ex(3),fRHS,HfRHS)
do k=ks,ex(3)-1
RK4 = 0
RHS = fRHS(k)
call rungekutta4_cplxscalar(dR,CTheta0,CTheta,RHS,RK4)
RK4 = 1
CTheta1 = HfRHS(k)
call rungekutta4_cplxscalar(dR,CTheta0,CTheta1,RHS,RK4)
call cswap(CTheta,CTheta1)
RK4 = 2
CTheta1 = HfRHS(k)
call rungekutta4_cplxscalar(dR,CTheta0,CTheta1,RHS,RK4)
call cswap(CTheta,CTheta1)
RK4 = 3
CTheta1 = fRHS(k+1)
call rungekutta4_cplxscalar(dR,CTheta0,CTheta1,RHS,RK4)
call cswap(CTheta0,CTheta1)
RTheta(i,j,k+1) = dreal(CTheta0)
ITheta(i,j,k+1) = dimag(CTheta0)
enddo
#if 0
! closing step
k = ex(3)-1
RHS = fRHS(k)
CTheta0 = dcmplx(RTheta(i,j,k),ITheta(i,j,k))+RHS*dR
RTheta(i,j,k+1) = dreal(CTheta0)
ITheta(i,j,k+1) = dimag(CTheta0)
#endif
enddo
enddo
gont = 0
return
end function NullEvol_Theta_givenx
!-----------------------------------------------------------------------------------------------
#if 1
! real evolve
! for eth_x, eth first, _x later
! this R is indeed x
function NullEvol_Theta(ex,crho,sigma,R,RJ,IJ,RU,IU,beta,RB,IB, &
Rnu,Inu,Rk,Ik,RTheta,ITheta,W,KK,HKK,KKx,HKKx,Rmin, &
qlR1,qlR2,qlI1,qlI2,quR1,quR2,quI1,quI2,gR,gI) result(gont)
implicit none
integer,intent(in ):: ex(1:3)
real*8,intent(in),dimension(ex(1))::crho
real*8,intent(in),dimension(ex(2))::sigma
real*8,intent(in),dimension(ex(3))::R
real*8, dimension(ex(1),ex(2),ex(3)),intent(in) :: beta,W,KK,HKK,KKx,HKKx
real*8, dimension(ex(1),ex(2),ex(3)),intent(in) :: RJ,IJ,RU,IU,RB,IB,Rnu,Inu,Rk,Ik
real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: RTheta,ITheta
real*8,intent(in) :: Rmin
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: qlR1,qlR2,qlI1,qlI2
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: quR1,quR2,quI1,quI2
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: gR,gI
! gont = 0: success; gont = 1: something wrong
integer::gont
double complex,dimension(ex(1),ex(2),ex(3)) :: CU,DCU,bDCU,CB,DCB,bDCB,CJ,DCJ
double complex :: CTheta0,CTheta,CTheta1,RHS
integer :: i,j,k,RK4
double complex,dimension(ex(3)) :: Cnu,Ck,HCnu,HCk,HCU,HDCU,CUx,HCUx,DCUx,HDCUx,HbDCU,bDCUx,HbDCUx
double complex,dimension(ex(3)) :: Cnux,HCnux,HCJ,HDCJ,CJx,HCJx,CJxx,HCJxx,DCJx,HDCJx,HCB,HDCB,HbDCB
real*8,dimension(ex(3)) :: Hbeta,betax,Hbetax,HW,Wx,HWx
double complex :: Theta_rhs,Theta_rhs_o
real*8 :: xx,dR
!!! sanity check
dR = sum(RJ)+sum(IJ)+sum(RU)+sum(IU)+sum(beta)+sum(RB)+sum(IB) + &
sum(Rnu)+sum(Inu)+sum(Rk)+sum(Ik)+sum(RTheta)+sum(ITheta) + &
sum(KK)+sum(HKK)+sum(KKx)+sum(HKKx)
if(dR.ne.dR) then
if(sum(RJ).ne.sum(RJ))write(*,*)"NullEvol_Theta: find NaN in RJ"
if(sum(IJ).ne.sum(IJ))write(*,*)"NullEvol_Theta: find NaN in IJ"
if(sum(RU).ne.sum(RU))write(*,*)"NullEvol_Theta: find NaN in RU"
if(sum(IU).ne.sum(IU))write(*,*)"NullEvol_Theta: find NaN in IU"
if(sum(beta).ne.sum(beta))write(*,*)"NullEvol_Theta: find NaN in beta"
if(sum(RB).ne.sum(RB))write(*,*)"NullEvol_Theta: find NaN in RB"
if(sum(IB).ne.sum(IB))write(*,*)"NullEvol_Theta: find NaN in IB"
if(sum(Rnu).ne.sum(Rnu))write(*,*)"NullEvol_Theta: find NaN in Rnu"
if(sum(Inu).ne.sum(Inu))write(*,*)"NullEvol_Theta: find NaN in Inu"
if(sum(Rk).ne.sum(Rk))write(*,*)"NullEvol_Theta: find NaN in Rk"
if(sum(Ik).ne.sum(Ik))write(*,*)"NullEvol_Theta: find NaN in Ik"
if(sum(RTheta).ne.sum(RTheta))write(*,*)"NullEvol_Theta: find NaN in RTheta"
if(sum(ITheta).ne.sum(ITheta))write(*,*)"NullEvol_Theta: find NaN in ITheta"
if(sum(KK).ne.sum(KK))write(*,*)"NullEvol_Theta: find NaN in KK"
if(sum(HKK).ne.sum(HKK))write(*,*)"NullEvol_Theta: find NaN in HKK"
if(sum(KKx).ne.sum(KKx))write(*,*)"NullEvol_Theta: find NaN in KKx"
if(sum(HKKx).ne.sum(HKKx))write(*,*)"NullEvol_Theta: find NaN in HKKx"
gont = 1
return
endif
dR = R(2) - R(1)
CU = dcmplx(RU,IU)
CB = dcmplx(RB,IB)
CJ = dcmplx(RJ,IJ)
do k=1,ex(3)
call derivs_eth(ex(1:2),crho,sigma,CU(:,:,k),DCU(:,:,k),1,1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
call derivs_eth(ex(1:2),crho,sigma,CU(:,:,k),bDCU(:,:,k),1,-1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
call derivs_eth(ex(1:2),crho,sigma,CB(:,:,k),DCB(:,:,k),1,1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
call derivs_eth(ex(1:2),crho,sigma,CB(:,:,k),bDCB(:,:,k),1,-1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
call derivs_eth(ex(1:2),crho,sigma,CJ(:,:,k),DCJ(:,:,k),2,1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
enddo
do j=1,ex(2)
do i=1,ex(1)
CTheta0 = dcmplx(RTheta(i,j,1),ITheta(i,j,1))
Cnu = dcmplx(Rnu(i,j,:),Inu(i,j,:))
Ck = dcmplx(Rk(i,j,:),Ik(i,j,:))
call cget_half_x(ex(3),CB(i,j,:),HCB)
call cget_half_x(ex(3),DCB(i,j,:),HDCB)
call cget_half_x(ex(3),bDCB(i,j,:),HbDCB)
call cget_half_x(ex(3),Cnu,HCnu)
call cderivs_x(ex(3),R,Cnu,Cnux)
call cget_half_x(ex(3),Cnux,HCnux)
call cget_half_x(ex(3),Ck,HCk)
call rget_half_x(ex(3),beta(i,j,:),Hbeta)
call rderivs_x(ex(3),R,beta(i,j,:),betax)
call rget_half_x(ex(3),betax,Hbetax)
call rderivs_x(ex(3),R,W,Wx)
call rget_half_x(ex(3),Wx,HWx)
call rget_half_x(ex(3),W(i,j,:),HW)
call cget_half_x(ex(3),CU(i,j,:),HCU)
call cderivs_x(ex(3),R,DCU(i,j,:),DCUx)
call cderivs_x(ex(3),R,CU(i,j,:),CUx)
call cget_half_x(ex(3),DCUx,HDCUx)
call cget_half_x(ex(3),CUx,HCUx)
call cget_half_x(ex(3),DCU(i,j,:),HDCU)
call cderivs_x(ex(3),R,bDCU(i,j,:),bDCUx)
call cget_half_x(ex(3),bDCUx,HbDCUx)
call cget_half_x(ex(3),bDCU(i,j,:),HbDCU)
call cderivs_x(ex(3),R,CJ(i,j,:),CJx)
call cdderivs_x(ex(3),R,CJ(i,j,:),CJxx)
call cget_half_x(ex(3),CJx,HCJx)
call cget_half_x(ex(3),CJxx,HCJxx)
call cderivs_x(ex(3),R,DCJ(i,j,:),DCJx)
call cget_half_x(ex(3),DCJx,HDCJx)
! old type code: PRD 54, 6153, Eq.(32) etc.
#if 0
! start up part
k = 1
RHS = Theta_rhs_o(R(k)+dR/2.d0,Rmin,Hbeta(k),betax(k),HKK(i,j,k),KKx(i,j,k),HCU(k),CUx(k),DCUx(k),HbDCU(k),bDCUx(k), &
HDCU(k),HCB(k),HDCB(k),HW(k),Wx(k),HCJ(k),HDCJ(k), &
CJx(k),CJxx(k),DCJx(k),HbDCB(k),HCnu(k),Cnux(k),HCk(k),CTheta0)
CTheta1 = RHS-(1-2.d0*(R(k)+dR/2.d0)*(1.d0-R(k)-dR/2.d0)/dR)*CTheta0
CTheta0 = CTheta1/(1+2.d0*(R(k)+dR/2.d0)*(1.d0-R(k)-dR/2.d0)/dR)
RTheta(i,j,k+1) = dreal(CTheta0)
ITheta(i,j,k+1) = dimag(CTheta0)
do k=1,ex(3)-2
RHS = Theta_rhs_o(R(k+1),Rmin,beta(i,j,k+1),betax(k+1),KK(i,j,k+1),KKx(i,j,k+1),CU(i,j,k+1),CUx(k+1),DCUx(k+1),bDCU(i,j,k+1),bDCUx(k+1), &
DCU(i,j,k+1),CB(i,j,k+1),DCB(i,j,k+1),W(i,j,k+1),Wx(k+1),CJ(i,j,k+1),DCJ(i,j,k+1), &
CJx(k+1),CJxx(k+1),DCJx(k+1),bDCB(i,j,k+1),Cnu(k+1),Cnux(k+1),Ck(k+1),CTheta0)
CTheta1 = RHS-(1-R(k+1)*(1.d0-R(k+1))/dR)*(dcmplx(RTheta(i,j,k),ITheta(i,j,k)))
CTheta0 = CTheta1/(1+R(k+1)*(1.d0-R(k+1))/dR)
RTheta(i,j,k+2) = dreal(CTheta0)
ITheta(i,j,k+2) = dimag(CTheta0)
enddo
#endif
#ifdef OLD
do k = 1,ex(3)-1
xx = R(k)+dR/2
! note our HCJ(ex(3)-1) = (CJ(ex(3))+CJ(ex(3)-1))/2
! note our KKx(ex(3)) = (KK(ex(3))-KK(ex(3)-1))/dR
! note our fxx(ex(3)) = (f(ex(3))-2.d0*f(ex(3)-1)+f(ex(3)-2))/dR
RHS = Theta_rhs(xx,Rmin,Hbeta(k),betax(k+1),HKK(i,j,k),KKx(i,j,k+1),HCU(k),CUx(k+1),DCUx(k+1),HbDCU(k),bDCUx(k+1), &
HDCU(k),HCB(k),HDCB(k),HW(k),Wx(k+1),HCJ(k),HDCJ(k), &
CJx(k+1),CJxx(k+1),DCJx(k+1),HbDCB(k),HCnu(k),Cnux(k+1),HCk(k),0)
RHS = RHS+CTheta0*(1.d0/dR-0.5d0/xx/(1.d0-xx))
CTheta0 = RHS/(1.d0/dR+0.5d0/xx/(1.d0-xx))
RTheta(i,j,k+1) = dreal(CTheta0)
ITheta(i,j,k+1) = dimag(CTheta0)
enddo
#else
do k=1,ex(3)-2
RK4 = 0
RHS = Theta_rhs(R(k),Rmin,beta(i,j,k),betax(k),KK(i,j,k),KKx(i,j,k),CU(i,j,k),CUx(k),DCUx(k),bDCU(i,j,k),bDCUx(k), &
DCU(i,j,k),CB(i,j,k),DCB(i,j,k),W(i,j,k),Wx(k),CJ(i,j,k),DCJ(i,j,k), &
CJx(k),CJxx(k),DCJx(k),bDCB(i,j,k),Cnu(k),Cnux(k),Ck(k),CTheta0)
call rungekutta4_cplxscalar(dR,CTheta0,CTheta,RHS,RK4)
RK4 = 1
CTheta1 = Theta_rhs(R(k)+dR/2.d0,Rmin,Hbeta(k),Hbetax(k),HKK(i,j,k),HKKx(i,j,k), &
HCU(k),HCUx(k),HDCUx(k),HbDCU(k),HbDCUx(k), &
HDCU(k),HCB(k),HDCB(k),HW(k),HWx(k),HCJ(k),HDCJ(k), &
HCJx(k),HCJxx(k),HDCJx(k),HbDCB(k),HCnu(k),HCnux(k),HCk(k),CTheta)
call rungekutta4_cplxscalar(dR,CTheta0,CTheta1,RHS,RK4)
call cswap(CTheta,CTheta1)
RK4 = 2
CTheta1 = Theta_rhs(R(k)+dR/2.d0,Rmin,Hbeta(k),Hbetax(k),HKK(i,j,k),HKKx(i,j,k), &
HCU(k),HCUx(k),HDCUx(k),HbDCU(k),HbDCUx(k), &
HDCU(k),HCB(k),HDCB(k),HW(k),HWx(k),HCJ(k),HDCJ(k), &
HCJx(k),HCJxx(k),HDCJx(k),HbDCB(k),HCnu(k),HCnux(k),HCk(k),CTheta)
call rungekutta4_cplxscalar(dR,CTheta0,CTheta1,RHS,RK4)
call cswap(CTheta,CTheta1)
RK4 = 3
CTheta1 = Theta_rhs(R(k+1),Rmin,beta(i,j,k+1),betax(k+1),KK(i,j,k+1),KKx(i,j,k+1), &
CU(i,j,k+1),CUx(k+1),DCUx(k+1),bDCU(i,j,k+1),bDCUx(k+1), &
DCU(i,j,k+1),CB(i,j,k+1),DCB(i,j,k+1),W(i,j,k+1),Wx(k+1),CJ(i,j,k+1),DCJ(i,j,k+1), &
CJx(k+1),CJxx(k+1),DCJx(k+1),bDCB(i,j,k+1),Cnu(k+1),Cnux(k+1),Ck(k+1),CTheta)
call rungekutta4_cplxscalar(dR,CTheta0,CTheta1,RHS,RK4)
call cswap(CTheta0,CTheta1)
RTheta(i,j,k+1) = dreal(CTheta0)
ITheta(i,j,k+1) = dimag(CTheta0)
enddo
#if 0
k = ex(3)
CTheta0 = -KK(i,j,k)*DCU(i,j,k)-(CU(i,j,k)*Cnu(k)+dconjg(CU(i,j,k))*DCJ(i,j,k))/2 &
+CJ(i,j,k)*(bDCU(i,j,k)-dconjg(bDCU(i,j,k)))/2 - W(i,j,k)*CJ(i,j,k)/2
RTheta(i,j,k) = dreal(CTheta0)
ITheta(i,j,k) = dimag(CTheta0)
#else
! closing step
k = ex(3)-1
CTheta0 = dcmplx(RTheta(i,j,k),ITheta(i,j,k))
xx = R(k)+dR/2
! note our HCJ(ex(3)-1) = (CJ(ex(3))+CJ(ex(3)-1))/2
! note our KKx(ex(3)) = (KK(ex(3))-KK(ex(3)-1))/dR
! note our fxx(ex(3)) = (f(ex(3))-2.d0*f(ex(3)-1)+f(ex(3)-2))/dR
RHS = Theta_rhs(xx,Rmin,Hbeta(k),betax(k+1),HKK(i,j,k),KKx(i,j,k+1),HCU(k),CUx(k+1),DCUx(k+1),HbDCU(k),bDCUx(k+1), &
HDCU(k),HCB(k),HDCB(k),HW(k),Wx(k+1),HCJ(k),HDCJ(k), &
CJx(k+1),CJxx(k+1),DCJx(k+1),HbDCB(k),HCnu(k),Cnux(k+1),HCk(k),dcmplx(0.d0,0.d0))
RHS = RHS+CTheta0*(1.d0/dR-0.5d0/xx/(1.d0-xx))
CTheta0 = RHS/(1.d0/dR+0.5d0/xx/(1.d0-xx))
RTheta(i,j,k+1) = dreal(CTheta0)
ITheta(i,j,k+1) = dimag(CTheta0)
#endif
#endif
enddo
enddo
gont = 0
return
end function NullEvol_Theta
!--------------------------------------------------------------------
! check with fake_Theta_rhs
#elif 0
! this R is indeed x
function NullEvol_Theta(ex,crho,sigma,R,RJ,IJ,RU,IU,beta,RB,IB, &
Rnu,Inu,Rk,Ik,RTheta,ITheta,W,KK,HKK,KKx,HKKx,Rmin, &
qlR1,qlR2,qlI1,qlI2,quR1,quR2,quI1,quI2,gR,gI) result(gont)
implicit none
integer,intent(in ):: ex(1:3)
real*8,intent(in),dimension(ex(1))::crho
real*8,intent(in),dimension(ex(2))::sigma
real*8,intent(in),dimension(ex(3))::R
real*8, dimension(ex(1),ex(2),ex(3)),intent(in) :: beta,W,KK,HKK,KKx,HKKx
real*8, dimension(ex(1),ex(2),ex(3)),intent(in) :: RJ,IJ,RU,IU,RB,IB,Rnu,Inu,Rk,Ik
real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: RTheta,ITheta
real*8,intent(in) :: Rmin
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: qlR1,qlR2,qlI1,qlI2
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: quR1,quR2,quI1,quI2
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: gR,gI
! gont = 0: success; gont = 1: something wrong
integer::gont
double complex,dimension(ex(1),ex(2),ex(3)) :: CU,DCU,bDCU,CB,DCB,bDCB,CJ,DCJ
double complex :: CTheta0,CTheta,CTheta1,RHS
integer :: i,j,k,RK4
integer,parameter :: ks=1
double complex,dimension(ex(3)) :: Cnu,Ck,HCnu,HCk,HCU,HDCU,CUx,HCUx,DCUx,HDCUx,HbDCU,bDCUx,HbDCUx
double complex,dimension(ex(3)) :: Cnux,HCnux,HCJ,HDCJ,CJx,HCJx,CJxx,HCJxx,DCJx,HDCJx,HCB,HDCB,HbDCB
real*8,dimension(ex(3)) :: Hbeta,betax,Hbetax,HW,Wx,HWx
double complex :: Theta_rhs,Theta_rhs_o
real*8 :: xx,dR
!!! sanity check
dR = sum(RJ)+sum(IJ)+sum(RU)+sum(IU)+sum(beta)+sum(RB)+sum(IB) + &
sum(Rnu)+sum(Inu)+sum(Rk)+sum(Ik)+sum(RTheta)+sum(ITheta) + &
sum(KK)+sum(HKK)+sum(KKx)+sum(HKKx)
if(dR.ne.dR) then
if(sum(RJ).ne.sum(RJ))write(*,*)"NullEvol_Theta: find NaN in RJ"
if(sum(IJ).ne.sum(IJ))write(*,*)"NullEvol_Theta: find NaN in IJ"
if(sum(RU).ne.sum(RU))write(*,*)"NullEvol_Theta: find NaN in RU"
if(sum(IU).ne.sum(IU))write(*,*)"NullEvol_Theta: find NaN in IU"
if(sum(beta).ne.sum(beta))write(*,*)"NullEvol_Theta: find NaN in beta"
if(sum(RB).ne.sum(RB))write(*,*)"NullEvol_Theta: find NaN in RB"
if(sum(IB).ne.sum(IB))write(*,*)"NullEvol_Theta: find NaN in IB"
if(sum(Rnu).ne.sum(Rnu))write(*,*)"NullEvol_Theta: find NaN in Rnu"
if(sum(Inu).ne.sum(Inu))write(*,*)"NullEvol_Theta: find NaN in Inu"
if(sum(Rk).ne.sum(Rk))write(*,*)"NullEvol_Theta: find NaN in Rk"
if(sum(Ik).ne.sum(Ik))write(*,*)"NullEvol_Theta: find NaN in Ik"
if(sum(RTheta).ne.sum(RTheta))write(*,*)"NullEvol_Theta: find NaN in RTheta"
if(sum(ITheta).ne.sum(ITheta))write(*,*)"NullEvol_Theta: find NaN in ITheta"
if(sum(KK).ne.sum(KK))write(*,*)"NullEvol_Theta: find NaN in KK"
if(sum(HKK).ne.sum(HKK))write(*,*)"NullEvol_Theta: find NaN in HKK"
if(sum(KKx).ne.sum(KKx))write(*,*)"NullEvol_Theta: find NaN in KKx"
if(sum(HKKx).ne.sum(HKKx))write(*,*)"NullEvol_Theta: find NaN in HKKx"
gont = 1
return
endif
dR = R(2) - R(1)
CU = dcmplx(RU,IU)
CB = dcmplx(RB,IB)
CJ = dcmplx(RJ,IJ)
do k=1,ex(3)
call derivs_eth(ex(1:2),crho,sigma,CU(:,:,k),DCU(:,:,k),1,1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
call derivs_eth(ex(1:2),crho,sigma,CU(:,:,k),bDCU(:,:,k),1,-1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
call derivs_eth(ex(1:2),crho,sigma,CB(:,:,k),DCB(:,:,k),1,1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
call derivs_eth(ex(1:2),crho,sigma,CB(:,:,k),bDCB(:,:,k),1,-1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
call derivs_eth(ex(1:2),crho,sigma,CJ(:,:,k),DCJ(:,:,k),2,1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
enddo
do j=1,ex(2)
do i=1,ex(1)
CTheta0 = dcmplx(RTheta(i,j,ks),ITheta(i,j,ks))
Cnu = dcmplx(RTheta(i,j,:),ITheta(i,j,:))
call fake_Theta_rhs(ex(3),R,Ck,Cnu)
call cget_half_x(ex(3),Ck,HCk)
do k=ks,ex(3)-1
RK4 = 0
RHS = Ck(k)
call rungekutta4_cplxscalar(dR,CTheta0,CTheta,RHS,RK4)
RK4 = 1
CTheta1 = HCk(k)
call rungekutta4_cplxscalar(dR,CTheta0,CTheta1,RHS,RK4)
call cswap(CTheta,CTheta1)
RK4 = 2
CTheta1 = HCk(k)
call rungekutta4_cplxscalar(dR,CTheta0,CTheta1,RHS,RK4)
call cswap(CTheta,CTheta1)
RK4 = 3
CTheta1 = Ck(k+1)
call rungekutta4_cplxscalar(dR,CTheta0,CTheta1,RHS,RK4)
call cswap(CTheta0,CTheta1)
RTheta(i,j,k+1) = dreal(CTheta0)
ITheta(i,j,k+1) = dimag(CTheta0)
enddo
enddo
enddo
gont = 0
return
end function NullEvol_Theta
#else
! for eth_x, _x first, eth second
! this R is indeed x
function NullEvol_Theta(ex,crho,sigma,R,RJ,IJ,RU,IU,beta,RB,IB, &
Rnu,Inu,Rk,Ik,RTheta,ITheta,W,KK,HKK,KKx,HKKx,Rmin, &
qlR1,qlR2,qlI1,qlI2,quR1,quR2,quI1,quI2,gR,gI) result(gont)
implicit none
integer,intent(in ):: ex(1:3)
real*8,intent(in),dimension(ex(1))::crho
real*8,intent(in),dimension(ex(2))::sigma
real*8,intent(in),dimension(ex(3))::R
real*8, dimension(ex(1),ex(2),ex(3)),intent(in) :: beta,W,KK,HKK,KKx,HKKx
real*8, dimension(ex(1),ex(2),ex(3)),intent(in) :: RJ,IJ,RU,IU,RB,IB,Rnu,Inu,Rk,Ik
real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: RTheta,ITheta
real*8,intent(in) :: Rmin
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: qlR1,qlR2,qlI1,qlI2
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: quR1,quR2,quI1,quI2
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: gR,gI
! gont = 0: success; gont = 1: something wrong
integer::gont
double complex,dimension(ex(1),ex(2),ex(3)) :: CU,DCU,bDCU,CB,DCB,bDCB,CJ,DCJ
double complex :: CTheta0,CTheta,CTheta1,RHS
integer :: i,j,k,RK4
double complex,dimension(ex(1),ex(2),ex(3)) :: Cnu,Ck,HCnu,HCk,HCU,HDCU,CUx,HCUx,DCUx,HDCUx,HbDCU,bDCUx,HbDCUx
double complex,dimension(ex(1),ex(2),ex(3)) :: Cnux,HCnux,HCJ,HDCJ,CJx,HCJx,CJxx,HCJxx,DCJx,HDCJx,HCB,HDCB,HbDCB
real*8,dimension(ex(1),ex(2),ex(3)) :: Hbeta,betax,Hbetax,HW,Wx,HWx
double complex :: Theta_rhs
real*8 :: xx,dR
!!! sanity check
dR = sum(RJ)+sum(IJ)+sum(RU)+sum(IU)+sum(beta)+sum(RB)+sum(IB) + &
sum(Rnu)+sum(Inu)+sum(Rk)+sum(Ik)+sum(RTheta)+sum(ITheta) + &
sum(KK)+sum(HKK)+sum(KKx)+sum(HKKx)
if(dR.ne.dR) then
if(sum(RJ).ne.sum(RJ))write(*,*)"NullEvol_Theta: find NaN in RJ"
if(sum(IJ).ne.sum(IJ))write(*,*)"NullEvol_Theta: find NaN in IJ"
if(sum(RU).ne.sum(RU))write(*,*)"NullEvol_Theta: find NaN in RU"
if(sum(IU).ne.sum(IU))write(*,*)"NullEvol_Theta: find NaN in IU"
if(sum(beta).ne.sum(beta))write(*,*)"NullEvol_Theta: find NaN in beta"
if(sum(RB).ne.sum(RB))write(*,*)"NullEvol_Theta: find NaN in RB"
if(sum(IB).ne.sum(IB))write(*,*)"NullEvol_Theta: find NaN in IB"
if(sum(Rnu).ne.sum(Rnu))write(*,*)"NullEvol_Theta: find NaN in Rnu"
if(sum(Inu).ne.sum(Inu))write(*,*)"NullEvol_Theta: find NaN in Inu"
if(sum(Rk).ne.sum(Rk))write(*,*)"NullEvol_Theta: find NaN in Rk"
if(sum(Ik).ne.sum(Ik))write(*,*)"NullEvol_Theta: find NaN in Ik"
if(sum(RTheta).ne.sum(RTheta))write(*,*)"NullEvol_Theta: find NaN in RTheta"
if(sum(ITheta).ne.sum(ITheta))write(*,*)"NullEvol_Theta: find NaN in ITheta"
if(sum(KK).ne.sum(KK))write(*,*)"NullEvol_Theta: find NaN in KK"
if(sum(HKK).ne.sum(HKK))write(*,*)"NullEvol_Theta: find NaN in HKK"
if(sum(KKx).ne.sum(KKx))write(*,*)"NullEvol_Theta: find NaN in KKx"
if(sum(HKKx).ne.sum(HKKx))write(*,*)"NullEvol_Theta: find NaN in HKKx"
gont = 1
return
endif
dR = R(2) - R(1)
CU = dcmplx(RU,IU)
CB = dcmplx(RB,IB)
CJ = dcmplx(RJ,IJ)
Cnu = dcmplx(Rnu,Inu)
Ck = dcmplx(Rk,Ik)
do j=1,ex(2)
do i=1,ex(1)
call cderivs_x(ex(3),R,Cnu(i,j,:),Cnux(i,j,:))
call rderivs_x(ex(3),R,beta(i,j,:),betax(i,j,:))
call rderivs_x(ex(3),R,W(i,j,:),Wx(i,j,:))
call cderivs_x(ex(3),R,CU(i,j,:),CUx(i,j,:))
call cderivs_x(ex(3),R,CJ(i,j,:),CJx(i,j,:))
call cdderivs_x(ex(3),R,CJ(i,j,:),CJxx(i,j,:))
enddo
enddo
do k=1,ex(3)
call derivs_eth(ex(1:2),crho,sigma,CU(:,:,k),DCU(:,:,k),1,1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
call derivs_eth(ex(1:2),crho,sigma,CUx(:,:,k),DCUx(:,:,k),1,1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
call derivs_eth(ex(1:2),crho,sigma,CU(:,:,k),bDCU(:,:,k),1,-1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
call derivs_eth(ex(1:2),crho,sigma,CUx(:,:,k),bDCUx(:,:,k),1,-1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
call derivs_eth(ex(1:2),crho,sigma,CB(:,:,k),DCB(:,:,k),1,1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
call derivs_eth(ex(1:2),crho,sigma,CB(:,:,k),bDCB(:,:,k),1,-1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
call derivs_eth(ex(1:2),crho,sigma,CJ(:,:,k),DCJ(:,:,k),2,1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
call derivs_eth(ex(1:2),crho,sigma,CJx(:,:,k),DCJx(:,:,k),2,1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
enddo
do j=1,ex(2)
do i=1,ex(1)
call cget_half_x(ex(3),CB(i,j,:),HCB(i,j,:))
call cget_half_x(ex(3),DCB(i,j,:),HDCB(i,j,:))
call cget_half_x(ex(3),bDCB(i,j,:),HbDCB(i,j,:))
call cget_half_x(ex(3),Cnu(i,j,:),HCnu(i,j,:))
call cget_half_x(ex(3),Cnux(i,j,:),HCnux(i,j,:))
call cget_half_x(ex(3),Ck(i,j,:),HCk(i,j,:))
call rget_half_x(ex(3),beta(i,j,:),Hbeta(i,j,:))
call rget_half_x(ex(3),betax(i,j,:),Hbetax(i,j,:))
call rget_half_x(ex(3),Wx(i,j,:),HWx(i,j,:))
call rget_half_x(ex(3),W(i,j,:),HW(i,j,:))
call cget_half_x(ex(3),CU(i,j,:),HCU(i,j,:))
call cget_half_x(ex(3),DCUx(i,j,:),HDCUx(i,j,:))
call cget_half_x(ex(3),CUx(i,j,:),HCUx(i,j,:))
call cget_half_x(ex(3),DCU(i,j,:),HDCU(i,j,:))
call cget_half_x(ex(3),bDCUx(i,j,:),HbDCUx(i,j,:))
call cget_half_x(ex(3),bDCU(i,j,:),HbDCU(i,j,:))
call cget_half_x(ex(3),CJx(i,j,:),HCJx(i,j,:))
call cget_half_x(ex(3),CJxx(i,j,:),HCJxx(i,j,:))
call cget_half_x(ex(3),DCJx(i,j,:),HDCJx(i,j,:))
enddo
enddo
do j=1,ex(2)
do i=1,ex(1)
CTheta0 = dcmplx(RTheta(i,j,1),ITheta(i,j,1))
do k=1,ex(3)-2
RK4 = 0
RHS = Theta_rhs(R(k),Rmin,beta(i,j,k),betax(i,j,k),KK(i,j,k),KKx(i,j,k),CU(i,j,k),CUx(i,j,k),DCUx(i,j,k),bDCU(i,j,k),bDCUx(i,j,k), &
DCU(i,j,k),CB(i,j,k),DCB(i,j,k),W(i,j,k),Wx(i,j,k),CJ(i,j,k),DCJ(i,j,k), &
CJx(i,j,k),CJxx(i,j,k),DCJx(i,j,k),bDCB(i,j,k),Cnu(i,j,k),Cnux(i,j,k),Ck(i,j,k),CTheta0)
call rungekutta4_cplxscalar(dR,CTheta0,CTheta,RHS,RK4)
RK4 = 1
CTheta1 = Theta_rhs(R(k)+dR/2.d0,Rmin,Hbeta(i,j,k),Hbetax(i,j,k),HKK(i,j,k),HKKx(i,j,k), &
HCU(i,j,k),HCUx(i,j,k),HDCUx(i,j,k),HbDCU(i,j,k),HbDCUx(i,j,k), &
HDCU(i,j,k),HCB(i,j,k),HDCB(i,j,k),HW(i,j,k),HWx(i,j,k),HCJ(i,j,k),HDCJ(i,j,k), &
HCJx(i,j,k),HCJxx(i,j,k),HDCJx(i,j,k),HbDCB(i,j,k),HCnu(i,j,k),HCnux(i,j,k),HCk(i,j,k),CTheta)
call rungekutta4_cplxscalar(dR,CTheta0,CTheta1,RHS,RK4)
call cswap(CTheta,CTheta1)
RK4 = 2
CTheta1 = Theta_rhs(R(k)+dR/2.d0,Rmin,Hbeta(i,j,k),Hbetax(i,j,k),HKK(i,j,k),HKKx(i,j,k), &
HCU(i,j,k),HCUx(i,j,k),HDCUx(i,j,k),HbDCU(i,j,k),HbDCUx(i,j,k), &
HDCU(i,j,k),HCB(i,j,k),HDCB(i,j,k),HW(i,j,k),HWx(i,j,k),HCJ(i,j,k),HDCJ(i,j,k), &
HCJx(i,j,k),HCJxx(i,j,k),HDCJx(i,j,k),HbDCB(i,j,k),HCnu(i,j,k),HCnux(i,j,k),HCk(i,j,k),CTheta)
call rungekutta4_cplxscalar(dR,CTheta0,CTheta1,RHS,RK4)
call cswap(CTheta,CTheta1)
RK4 = 3
CTheta1 = Theta_rhs(R(k+1),Rmin,beta(i,j,k+1),betax(i,j,k+1),KK(i,j,k+1), &
KKx(i,j,k+1),CU(i,j,k+1),CUx(i,j,k+1),DCUx(i,j,k+1),bDCU(i,j,k+1),bDCUx(i,j,k+1), &
DCU(i,j,k+1),CB(i,j,k+1),DCB(i,j,k+1),W(i,j,k+1),Wx(i,j,k+1),CJ(i,j,k+1),DCJ(i,j,k+1), &
CJx(i,j,k+1),CJxx(i,j,k+1),DCJx(i,j,k+1),bDCB(i,j,k+1),Cnu(i,j,k+1),Cnux(i,j,k+1),Ck(i,j,k+1),CTheta)
call rungekutta4_cplxscalar(dR,CTheta0,CTheta1,RHS,RK4)
call cswap(CTheta0,CTheta1)
RTheta(i,j,k+1) = dreal(CTheta0)
ITheta(i,j,k+1) = dimag(CTheta0)
enddo
k = ex(3)
CTheta0 = -KK(i,j,k)*DCU(i,j,k)-(CU(i,j,k)*Cnu(i,j,k)+dconjg(CU(i,j,k))*DCJ(i,j,k))/2 &
+CJ(i,j,k)*(bDCU(i,j,k)-dconjg(bDCU(i,j,k)))/2 - W(i,j,k)*CJ(i,j,k)/2
RTheta(i,j,k) = dreal(CTheta0)
ITheta(i,j,k) = dimag(CTheta0)
enddo
enddo
gont = 0
return
end function NullEvol_Theta
#endif
#elif (RKorAM == 1)
!------------------------------------------------------------------------------
! this R is indeed x
function NullEvol_beta(ex,crho,sigma,R,RJ,IJ,beta,KKx,HKKx) result(gont)
implicit none
integer,intent(in ):: ex(1:3)
real*8,intent(in),dimension(ex(1))::crho
real*8,intent(in),dimension(ex(2))::sigma
real*8,intent(in),dimension(ex(3))::R
real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: beta
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: RJ,IJ,KKx,HKKx
! gont = 0: success; gont = 1: something wrong
integer::gont
real*8 :: dR,beta_rhs
double complex, dimension(ex(3)):: CJ,CJx
real*8, dimension(ex(3)) :: rhs
integer :: i,j,k
real*8,parameter :: F5o12=2.d0/1.2d1,F2o3=2.d0/3.d0,F1o12=1.d0/1.2d1
real*8,parameter :: F3o8=3.d0/8.d0,F19o24=1.9d1/2.4d1,F5o24=5.d0/2.4d1,F1o24=1.d0/2.4d1
!!! sanity check
dR = sum(RJ)+sum(IJ)+sum(beta)+sum(KKx)+sum(HKKx)
if(dR.ne.dR) then
if(sum(RJ).ne.sum(RJ))write(*,*)"NullEvol_beta: find NaN in RJ"
if(sum(IJ).ne.sum(IJ))write(*,*)"NullEvol_beta: find NaN in IJ"
if(sum(beta).ne.sum(beta))write(*,*)"NullEvol_beta: find NaN in beta"
if(sum(KKx).ne.sum(KKx))write(*,*)"NullEvol_beta: find NaN in KKx"
if(sum(HKKx).ne.sum(HKKx))write(*,*)"NullEvol_beta: find NaN in HKKx"
gont = 1
return
endif
dR = R(2) - R(1)
do j=1,ex(2)
do i=1,ex(1)
CJ = dcmplx(RJ(i,j,:),IJ(i,j,:))
#if 0
call cderivs_sw_x(ex(3),R,CJ,CJx)
#else
call cderivs_x(ex(3),R,CJ,CJx)
#endif
do k=1,ex(3)
rhs(k) = beta_rhs(R(k),CJx(k),KKx(i,j,k))
enddo
k = 1
beta(i,j,k+1) = beta(i,j,k) + (rhs(k+1)+rhs(k))*dR/2
k = 2
beta(i,j,k+1) = beta(i,j,k) + (F5o12*rhs(k+1) + F2o3*rhs(k) - F1o12*rhs(k-1))*dR
do k=3,ex(3)-1
beta(i,j,k+1) = beta(i,j,k) + (F3o8*rhs(k+1) + F19o24*rhs(k) - F5o24*rhs(k-1) + F1o24*rhs(k-2))*dR
enddo
enddo
enddo
gont = 0
return
end function NullEvol_beta
!------------------------------------------------------------------------------
! this R is indeed x
function NullEvol_Q(ex,crho,sigma,R,RJ,IJ,Rk,Ik,Rnu,Inu,RB,IB,RQ,IQ,KK,Hkk,KKx,HKKx, &
qlR1,qlR2,qlI1,qlI2,quR1,quR2,quI1,quI2,gR,gI) result(gont)
implicit none
integer,intent(in ):: ex(1:3)
real*8,intent(in),dimension(ex(1))::crho
real*8,intent(in),dimension(ex(2))::sigma
real*8,intent(in),dimension(ex(3))::R
real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: RQ,IQ
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: RJ,IJ,KK,KKx,HKK,HKKx
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: Rk,Ik,Rnu,Inu,RB,IB
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: qlR1,qlR2,qlI1,qlI2
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: quR1,quR2,quI1,quI2
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: gR,gI
! gont = 0: success; gont = 1: something wrong
integer::gont
real*8 :: xx,dR
double complex,dimension(ex(3)) :: CQ,RHS
real*8, dimension(ex(3)) :: gunc
double complex,dimension(ex(3)) :: CJx,DCJx,Ck,Ckx,Cnu,Cnux,CB,CBx
double complex, dimension(ex(1),ex(2),ex(3)) :: CJ,DCJ
integer :: i,j,k
double complex :: ZEO,Q_rhs
real*8,parameter :: F5o12=2.d0/1.2d1,F2o3=2.d0/3.d0,F1o12=1.d0/1.2d1
real*8,parameter :: F3o8=3.d0/8.d0,F19o24=1.9d1/2.4d1,F5o24=5.d0/2.4d1,F1o24=1.d0/2.4d1
!!! sanity check
dR = sum(RJ)+sum(IJ)+sum(RQ)+sum(IQ) &
+sum(RK)+sum(IK)+sum(Rnu)+sum(Inu)+sum(RB)+sum(IB) &
+sum(KK)+sum(KKx)
if(dR.ne.dR) then
if(sum(RJ).ne.sum(RJ))write(*,*)"NullEvol_Q: find NaN in RJ"
if(sum(IJ).ne.sum(IJ))write(*,*)"NullEvol_Q: find NaN in IJ"
if(sum(RQ).ne.sum(RQ))write(*,*)"NullEvol_Q: find NaN in RQ"
if(sum(IQ).ne.sum(IQ))write(*,*)"NullEvol_Q: find NaN in IQ"
if(sum(RK).ne.sum(RK))write(*,*)"NullEvol_Q: find NaN in RK"
if(sum(IK).ne.sum(IK))write(*,*)"NullEvol_Q: find NaN in IK"
if(sum(Rnu).ne.sum(Rnu))write(*,*)"NullEvol_Q: find NaN in Rnu"
if(sum(Inu).ne.sum(Inu))write(*,*)"NullEvol_Q: find NaN in Inu"
if(sum(RB).ne.sum(RB))write(*,*)"NullEvol_Q: find NaN in RB"
if(sum(IB).ne.sum(IB))write(*,*)"NullEvol_Q: find NaN in IB"
if(sum(KK).ne.sum(KK))write(*,*)"NullEvol_Q: find NaN in KK"
if(sum(KKx).ne.sum(KKx))write(*,*)"NullEvol_Q: find NaN in KKx"
gont = 1
return
endif
dR = R(2) - R(1)
ZEO = dcmplx(0.d0,0.d0)
CJ = dcmplx(RJ,IJ)
do k=1,ex(3)
call derivs_eth(ex(1:2),crho,sigma,CJ(:,:,k),DCJ(:,:,k),2,1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
enddo
do j=1,ex(2)
do i=1,ex(1)
CQ(1) = dcmplx(RQ(i,j,1),IQ(i,j,1))
Ck = dcmplx(Rk(i,j,:),Ik(i,j,:))
Cnu = dcmplx(Rnu(i,j,:),Inu(i,j,:))
CB = dcmplx(RB(i,j,:),IB(i,j,:))
#if 0
call cderivs_sw_x(ex(3),R,CJ(i,j,:),CJx)
call cderivs_sw_x(ex(3),R,DCJ(i,j,:),DCJx)
call cderivs_sw_x(ex(3),R,Ck,Ckx)
call cderivs_sw_x(ex(3),R,Cnu,Cnux)
call cderivs_sw_x(ex(3),R,CB,CBx)
#else
call cderivs_x(ex(3),R,CJ(i,j,:),CJx)
call cderivs_x(ex(3),R,DCJ(i,j,:),DCJx)
call cderivs_x(ex(3),R,Ck,Ckx)
call cderivs_x(ex(3),R,Cnu,Cnux)
call cderivs_x(ex(3),R,CB,CBx)
#endif
do k = 1,ex(3)
RHS(k) = Q_rhs(R(k),CJ(i,j,k),CJx(k),DCJx(k),KK(i,j,k),Ck(k),Ckx(k),Cnux(k),KKx(i,j,k),CBx(k),Cnu(k),DCJ(i,j,k),CB(k),ZEO)
gunc(k) = -2/R(k)/(1-R(k))
enddo
k = 1
CQ(k+1) = CQ(k) + (RHS(k+1)+RHS(k)+CQ(k)*gunc(k))*dR/2
CQ(k+1) = CQ(k+1)/(1-0.5*dR*gunc(k+1))
k = 2
CQ(k+1) = CQ(k) + (F5o12*RHS(k+1) + F2o3*(RHS(k)+CQ(k)*gunc(k)) - F1o12*(RHS(k-1)+CQ(k-1)*gunc(k-1)))*dR
CQ(k+1) = CQ(k+1)/(1-F5o12*dR*gunc(k+1))
do k=3,ex(3)-2
CQ(k+1) = CQ(k) + (F3o8*RHS(k+1) + F19o24*(RHS(k)+CQ(k)*gunc(k)) - F5o24*(RHS(k-1)+CQ(k-1)*gunc(k-1)) &
+ F1o24*(RHS(k-2)+CQ(k-2)*gunc(k-2)))*dR
CQ(k+1) = CQ(k+1)/(1-F3o8*dR*gunc(k+1))
enddo
k = ex(3)
CQ(k) = -2*CB(k)
RQ(i,j,:) = dreal(CQ)
IQ(i,j,:) = dimag(CQ)
enddo
enddo
gont = 0
return
end function NullEvol_Q
!--------------------------------------------------------------------
! this R is indeed x
function NullEvol_U(ex,crho,sigma,R,RJ,IJ,RQ,IQ,KK,HKK,beta,RU,IU, &
Rmin) result(gont)
implicit none
integer,intent(in ):: ex(1:3)
real*8,intent(in),dimension(ex(1))::crho
real*8,intent(in),dimension(ex(2))::sigma
real*8,intent(in),dimension(ex(3))::R
real*8, dimension(ex(1),ex(2),ex(3)),intent(in) :: RJ,IJ,RQ,IQ,beta,KK,HKK
real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: RU,IU
real*8,intent(in) :: Rmin
! gont = 0: success; gont = 1: something wrong
integer::gont
real*8 :: dR
double complex,dimension(ex(3)) :: CU0,RHS
integer :: i,j,k
double complex :: U_rhs
double complex,dimension(ex(3)) :: CJ,CQ
real*8,parameter :: F5o12=2.d0/1.2d1,F2o3=2.d0/3.d0,F1o12=1.d0/1.2d1
real*8,parameter :: F3o8=3.d0/8.d0,F19o24=1.9d1/2.4d1,F5o24=5.d0/2.4d1,F1o24=1.d0/2.4d1
!!! sanity check
dR = sum(RJ)+sum(IJ)+sum(RQ)+sum(IQ)+sum(beta)+sum(RU)+sum(IU)+sum(KK)
if(dR.ne.dR) then
if(sum(RJ).ne.sum(RJ))write(*,*)"NullEvol_U: find NaN in RJ"
if(sum(IJ).ne.sum(IJ))write(*,*)"NullEvol_U: find NaN in IJ"
if(sum(RQ).ne.sum(RQ))write(*,*)"NullEvol_U: find NaN in RQ"
if(sum(IQ).ne.sum(IQ))write(*,*)"NullEvol_U: find NaN in IQ"
if(sum(beta).ne.sum(beta))write(*,*)"NullEvol_U: find NaN in beta"
if(sum(RU).ne.sum(RU))write(*,*)"NullEvol_U: find NaN in RU"
if(sum(IU).ne.sum(IU))write(*,*)"NullEvol_U: find NaN in IU"
if(sum(KK).ne.sum(KK))write(*,*)"NullEvol_U: find NaN in KK"
gont = 1
return
endif
dR = R(2) - R(1)
do j=1,ex(2)
do i=1,ex(1)
CU0(1) = dcmplx(RU(i,j,1),IU(i,j,1))
CJ = dcmplx(RJ(i,j,:),IJ(i,j,:))
CQ = dcmplx(RQ(i,j,:),IQ(i,j,:))
do k = 1,ex(3)
RHS(k) = U_rhs(R(k),Rmin,beta(i,j,k),KK(i,j,k),CQ(k),CJ(k))
enddo
k = 1
CU0(k+1) = CU0(k) + (RHS(k+1)+RHS(k))*dR/2
k = 2
CU0(k+1) = CU0(k) + (F5o12*RHS(k+1) + F2o3*RHS(k) - F1o12*RHS(k-1))*dR
do k=3,ex(3)-1
CU0(k+1) = CU0(k) + (F3o8*RHS(k+1) + F19o24*RHS(k) - F5o24*RHS(k-1) &
+ F1o24*RHS(k-2))*dR
enddo
RU(i,j,:) = dreal(CU0)
IU(i,j,:) = dimag(CU0)
enddo
enddo
gont = 0
return
end function NullEvol_U
!----------------------------------------------------------------------------------------
! this R is indeed x
function NullEvol_W(ex,crho,sigma,R,RJ,IJ,RB,IB,Rnu,Inu,Rk,Ik, &
RU,IU,RQ,IQ,W,beta,KK,HKK,Rmin, &
qlR1,qlR2,qlI1,qlI2,quR1,quR2,quI1,quI2,gR,gI) result(gont)
implicit none
integer,intent(in ):: ex(1:3)
real*8,intent(in),dimension(ex(1))::crho
real*8,intent(in),dimension(ex(2))::sigma
real*8,intent(in),dimension(ex(3))::R
real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: W
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: RJ,IJ,RB,IB
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: Rnu,Inu,Rk,Ik
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: RU,IU,RQ,IQ,beta,KK,HKK
real*8,intent(in ) :: Rmin
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: qlR1,qlR2,qlI1,qlI2
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: quR1,quR2,quI1,quI2
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: gR,gI
! gont = 0: success; gont = 1: something wrong
integer::gont
real*8 :: dR
double complex, dimension(ex(1),ex(2),ex(3)) :: CU,DCU,bDCU
double complex, dimension(ex(1),ex(2),ex(3)) :: CB,DCB,bDCB,CJ,DCJ,Cnu,bDCnu,Ck,bDCk
double complex, dimension(ex(3)) :: bDCUx,CQ
integer :: i,j,k
real*8, dimension(ex(3)) :: rhs,gunc
real*8 :: zeo,W_rhs
real*8,parameter :: F5o12=2.d0/1.2d1,F2o3=2.d0/3.d0,F1o12=1.d0/1.2d1
real*8,parameter :: F3o8=3.d0/8.d0,F19o24=1.9d1/2.4d1,F5o24=5.d0/2.4d1,F1o24=1.d0/2.4d1
!!! sanity check
dR = sum(RJ)+sum(IJ)+sum(beta)+sum(RB)+sum(IB)+sum(Rnu)+sum(Inu) &
+sum(Rk)+sum(Ik)+sum(W)+sum(RU)+sum(IU)+sum(RQ)+sum(IQ)&
+sum(KK)
if(dR.ne.dR) then
if(sum(RJ).ne.sum(RJ))write(*,*)"NullEvol_W: find NaN in RJ"
if(sum(IJ).ne.sum(IJ))write(*,*)"NullEvol_W: find NaN in IJ"
if(sum(beta).ne.sum(beta))write(*,*)"NullEvol_W: find NaN in beta"
if(sum(RB).ne.sum(RB))write(*,*)"NullEvol_W: find NaN in RB"
if(sum(IB).ne.sum(IB))write(*,*)"NullEvol_W: find NaN in IB"
if(sum(Rnu).ne.sum(Rnu))write(*,*)"NullEvol_W: find NaN in Rnu"
if(sum(Inu).ne.sum(Inu))write(*,*)"NullEvol_W: find NaN in Inu"
if(sum(Rk).ne.sum(Rk))write(*,*)"NullEvol_W: find NaN in Rk"
if(sum(Ik).ne.sum(Ik))write(*,*)"NullEvol_W: find NaN in Ik"
if(sum(RU).ne.sum(RU))write(*,*)"NullEvol_W: find NaN in RU"
if(sum(IU).ne.sum(IU))write(*,*)"NullEvol_W: find NaN in IU"
if(sum(RQ).ne.sum(RQ))write(*,*)"NullEvol_W: find NaN in RQ"
if(sum(IQ).ne.sum(IQ))write(*,*)"NullEvol_W: find NaN in IQ"
if(sum(W).ne.sum(W))write(*,*)"NullEvol_W: find NaN in W"
if(sum(KK).ne.sum(KK))write(*,*)"NullEvol_W: find NaN in KK"
gont = 1
return
endif
dR = R(2) - R(1)
zeo = 0.d0
CB = dcmplx(RB,IB)
CU = dcmplx(RU,IU)
Ck = dcmplx(Rk,Ik)
Cnu = dcmplx(Rnu,Inu)
CJ = dcmplx(RJ,IJ)
do k=1,ex(3)
call derivs_eth(ex(1:2),crho,sigma,CJ(:,:,k),DCJ(:,:,k),2,1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
call derivs_eth(ex(1:2),crho,sigma,CB(:,:,k),DCB(:,:,k),1,1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
call derivs_eth(ex(1:2),crho,sigma,CB(:,:,k),bDCB(:,:,k),1,-1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
call derivs_eth(ex(1:2),crho,sigma,CU(:,:,k),bDCU(:,:,k),1,-1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
call derivs_eth(ex(1:2),crho,sigma,Ck(:,:,k),bDCk(:,:,k),1,-1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
call derivs_eth(ex(1:2),crho,sigma,Cnu(:,:,k),bDCnu(:,:,k),1,-1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
enddo
do j=1,ex(2)
do i=1,ex(1)
#if 0
call cderivs_sw_x(ex(3),R,bDCU,bDCUx)
#else
call cderivs_x(ex(3),R,bDCU,bDCUx)
#endif
CQ = dcmplx(RQ(i,j,:),IQ(i,j,:))
do k = 1,ex(3)
rhs(k) = W_rhs(R(k),Rmin,beta(i,j,k),KK(i,j,k),DCB(i,j,k),CB(i,j,k),CJ(i,j,k),Cnu(i,j,k),Ck(i,j,k),zeo, &
CQ(k),bDCk(i,j,k),bDCnu(i,j,k),bDCB(i,j,k),bDCU(i,j,k),bDCUx(k),DCJ(i,j,k))
gunc(k) = -2/R(k)/(1-R(k))
enddo
k = 1
W(i,j,k+1) = W(i,j,k) + (rhs(k+1)+rhs(k)+W(i,j,k)*gunc(k))*dR/2
W(i,j,k+1) = W(i,j,k+1)/(1-0.5*dR*gunc(k+1))
k = 2
W(i,j,k+1) = W(i,j,k) + (F5o12*rhs(k+1) + F2o3*(rhs(k)+W(i,j,k)*gunc(k)) - F1o12*(rhs(k-1)+W(i,j,k-1)*gunc(k-1)))*dR
W(i,j,k+1) = W(i,j,k+1)/(1-F5o12*dR*gunc(k+1))
do k=3,ex(3)-2
W(i,j,k+1) = W(i,j,k) + (F3o8*rhs(k+1) + F19o24*(rhs(k)+W(i,j,k)*gunc(k)) - F5o24*(rhs(k-1)+W(i,j,k-1)*gunc(k-1)) &
+ F1o24*(rhs(k-2)+W(i,j,k-2)*gunc(k-2)))*dR
W(i,j,k+1) = W(i,j,k+1)/(1-F3o8*dR*gunc(k+1))
enddo
k = ex(3)
W(i,j,k) = dreal(bDCU(i,j,k))
enddo
enddo
gont = 0
return
end function NullEvol_W
!--------------------------------------------------------------------
! this R is indeed x
function NullEvol_Theta(ex,crho,sigma,R,RJ,IJ,RU,IU,beta,RB,IB, &
Rnu,Inu,Rk,Ik,RTheta,ITheta,W,KK,HKK,KKx,HKKx,Rmin, &
qlR1,qlR2,qlI1,qlI2,quR1,quR2,quI1,quI2,gR,gI) result(gont)
implicit none
integer,intent(in ):: ex(1:3)
real*8,intent(in),dimension(ex(1))::crho
real*8,intent(in),dimension(ex(2))::sigma
real*8,intent(in),dimension(ex(3))::R
real*8, dimension(ex(1),ex(2),ex(3)),intent(in) :: beta,W,KK,KKx,HKK,HKKx
real*8, dimension(ex(1),ex(2),ex(3)),intent(in) :: RJ,IJ,RU,IU,RB,IB,Rnu,Inu,Rk,Ik
real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: RTheta,ITheta
real*8,intent(in) :: Rmin
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: qlR1,qlR2,qlI1,qlI2
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: quR1,quR2,quI1,quI2
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: gR,gI
! gont = 0: success; gont = 1: something wrong
integer::gont
double complex,dimension(ex(1),ex(2),ex(3)) :: CU,DCU,bDCU,CB,DCB,bDCB,CJ,DCJ
double complex,dimension(ex(3)) :: CTheta0,RHS
integer :: i,j,k,RK4
double complex,dimension(ex(3)) :: Cnu,Ck,CUx,DCUx,bDCUx
double complex,dimension(ex(3)) :: Cnux,CJx,CJxx,DCJx
real*8,dimension(ex(3)) :: betax,Wx,gunc
double complex :: Theta_rhs,ZEO
real*8 :: dR
real*8,parameter :: F5o12=2.d0/1.2d1,F2o3=2.d0/3.d0,F1o12=1.d0/1.2d1
real*8,parameter :: F3o8=3.d0/8.d0,F19o24=1.9d1/2.4d1,F5o24=5.d0/2.4d1,F1o24=1.d0/2.4d1
!!! sanity check
dR = sum(RJ)+sum(IJ)+sum(RU)+sum(IU)+sum(beta)+sum(RB)+sum(IB) + &
sum(Rnu)+sum(Inu)+sum(Rk)+sum(Ik)+sum(RTheta)+sum(ITheta) + &
sum(KK)+sum(KKx)+sum(W)
if(dR.ne.dR) then
if(sum(RJ).ne.sum(RJ))write(*,*)"NullEvol_Theta: find NaN in RJ"
if(sum(IJ).ne.sum(IJ))write(*,*)"NullEvol_Theta: find NaN in IJ"
if(sum(RU).ne.sum(RU))write(*,*)"NullEvol_Theta: find NaN in RU"
if(sum(IU).ne.sum(IU))write(*,*)"NullEvol_Theta: find NaN in IU"
if(sum(beta).ne.sum(beta))write(*,*)"NullEvol_Theta: find NaN in beta"
if(sum(RB).ne.sum(RB))write(*,*)"NullEvol_Theta: find NaN in RB"
if(sum(IB).ne.sum(IB))write(*,*)"NullEvol_Theta: find NaN in IB"
if(sum(Rnu).ne.sum(Rnu))write(*,*)"NullEvol_Theta: find NaN in Rnu"
if(sum(Inu).ne.sum(Inu))write(*,*)"NullEvol_Theta: find NaN in Inu"
if(sum(Rk).ne.sum(Rk))write(*,*)"NullEvol_Theta: find NaN in Rk"
if(sum(Ik).ne.sum(Ik))write(*,*)"NullEvol_Theta: find NaN in Ik"
if(sum(RTheta).ne.sum(RTheta))write(*,*)"NullEvol_Theta: find NaN in RTheta"
if(sum(ITheta).ne.sum(ITheta))write(*,*)"NullEvol_Theta: find NaN in ITheta"
if(sum(KK).ne.sum(KK))write(*,*)"NullEvol_Theta: find NaN in KK"
if(sum(KKx).ne.sum(KKx))write(*,*)"NullEvol_Theta: find NaN in KKx"
if(sum(W).ne.sum(W))write(*,*)"NullEvol_Theta: find NaN in W"
gont = 1
return
endif
dR = R(2) - R(1)
ZEO = dcmplx(0.d0,0.d0)
CU = dcmplx(RU,IU)
CB = dcmplx(RB,IB)
CJ = dcmplx(RJ,IJ)
do k=1,ex(3)
call derivs_eth(ex(1:2),crho,sigma,CU(:,:,k),DCU(:,:,k),1,1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
call derivs_eth(ex(1:2),crho,sigma,CU(:,:,k),bDCU(:,:,k),1,-1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
call derivs_eth(ex(1:2),crho,sigma,CB(:,:,k),DCB(:,:,k),1,1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
call derivs_eth(ex(1:2),crho,sigma,CB(:,:,k),bDCB(:,:,k),1,-1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
call derivs_eth(ex(1:2),crho,sigma,CJ(:,:,k),DCJ(:,:,k),2,1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
enddo
do j=1,ex(2)
do i=1,ex(1)
CTheta0(1) = dcmplx(RTheta(i,j,1),ITheta(i,j,1))
Cnu = dcmplx(Rnu(i,j,:),Inu(i,j,:))
Ck = dcmplx(Rk(i,j,:),Ik(i,j,:))
#if 0
call cderivs_sw_x(ex(3),R,Cnu,Cnux)
call rderivs_sw_x(ex(3),R,beta(i,j,:),betax)
call rderivs_sw_x(ex(3),R,W,Wx)
call cderivs_sw_x(ex(3),R,DCU(i,j,:),DCUx)
call cderivs_sw_x(ex(3),R,CU(i,j,:),CUx)
call cderivs_sw_x(ex(3),R,bDCU(i,j,:),bDCUx)
call cderivs_sw_x(ex(3),R,CJ(i,j,:),CJx)
call cdderivs_sw_x(ex(3),R,CJ(i,j,:),CJxx)
call cderivs_sw_x(ex(3),R,DCJ(i,j,:),DCJx)
#else
call cderivs_x(ex(3),R,Cnu,Cnux)
call rderivs_x(ex(3),R,beta(i,j,:),betax)
call rderivs_x(ex(3),R,W,Wx)
call cderivs_x(ex(3),R,DCU(i,j,:),DCUx)
call cderivs_x(ex(3),R,CU(i,j,:),CUx)
call cderivs_x(ex(3),R,bDCU(i,j,:),bDCUx)
call cderivs_x(ex(3),R,CJ(i,j,:),CJx)
call cdderivs_x(ex(3),R,CJ(i,j,:),CJxx)
call cderivs_x(ex(3),R,DCJ(i,j,:),DCJx)
#endif
do k = 1,ex(3)
rhs(k) = Theta_rhs(R(k),Rmin,beta(i,j,k),betax(k),KK(i,j,k),KKx(i,j,k),CU(i,j,k),CUx(k),DCUx(k),bDCU(i,j,k),bDCUx(k), &
DCU(i,j,k),CB(i,j,k),DCB(i,j,k),W(i,j,k),Wx(k),CJ(i,j,k),DCJ(i,j,k), &
CJx(k),CJxx(k),DCJx(k),bDCB(i,j,k),Cnu(k),Cnux(k),Ck(k),ZEO)
gunc(k) = -1/R(k)/(1-R(k))
enddo
k = 1
CTheta0(k+1) = CTheta0(k) + (RHS(k+1)+RHS(k)+CTheta0(k)*gunc(k))*dR/2
CTheta0(k+1) = CTheta0(k+1)/(1-0.5*dR*gunc(k+1))
k = 2
CTheta0(k+1) = CTheta0(k) + (F5o12*RHS(k+1) + F2o3*(RHS(k)+CTheta0(k)*gunc(k)) - F1o12*(RHS(k-1)+CTheta0(k-1)*gunc(k-1)))*dR
CTheta0(k+1) = CTheta0(k+1)/(1-F5o12*dR*gunc(k+1))
do k=3,ex(3)-2
CTheta0(k+1) = CTheta0(k) + (F3o8*RHS(k+1) + F19o24*(RHS(k)+CTheta0(k)*gunc(k)) - F5o24*(RHS(k-1)+CTheta0(k-1)*gunc(k-1)) &
+ F1o24*(RHS(k-2)+CTheta0(k-2)*gunc(k-2)))*dR
CTheta0(k+1) = CTheta0(k+1)/(1-F3o8*dR*gunc(k+1))
enddo
k = ex(3)
CTheta0(k) = -KK(i,j,k)*DCU(i,j,k)-(CU(i,j,k)*Cnu(k)+dconjg(CU(i,j,k))*DCJ(i,j,k))/2 &
+CJ(i,j,k)*(bDCU(i,j,k)-dconjg(bDCU(i,j,k)))/2 - W(i,j,k)*CJ(i,j,k)/2
RTheta(i,j,:) = dreal(CTheta0)
ITheta(i,j,:) = dimag(CTheta0)
enddo
enddo
gont = 0
return
end function NullEvol_Theta
#else
#error "not recognized RKorAM"
#endif
!=====================================================================================================================================
! basic tool routines
subroutine rswap(r1,r2)
implicit none
!~~~~~~% Input parameters:
real*8,intent(inout) :: r1,r2
real*8 :: r
r = r1
r1= r2
r2= r
return
end subroutine rswap
!----
subroutine cswap(r1,r2)
implicit none
!~~~~~~% Input parameters:
double complex,intent(inout) :: r1,r2
double complex :: r
r = r1
r1= r2
r2= r
return
end subroutine cswap
! center type finite difference
!====================================================================================
!----
subroutine rderivs_x(lx,X,f,fx)
implicit none
!~~~~~~% Input parameters:
integer,intent(in) :: lx
real*8,intent(in),dimension(lx) :: X
real*8,intent(in),dimension(lx) :: f
real*8,intent(out),dimension(lx) :: fx
real*8 :: dX
dX = X(2)-X(1)
#ifdef OLD
fx(1:lx-1) = (f(2:lx)-f(1:lx-1))/dX
fx(lx) = (f(lx)-f(lx-1))/dX
#else
#if (ghost_width == 2)
fx(2:lx-1) = (f(3:lx)-f(1:lx-2))/2.d0/dX
fx(1) = (f(2)-f(1))/dX
fx(lx) = (f(lx)-f(lx-1))/dX
#elif (ghost_width == 3)
fx(3:lx-2) = (f(1:lx-4)-8.d0*f(2:lx-3)+8.d0*f(4:lx-1)-f(5:lx))/1.2d1/dX
fx(2) = (f(3)-f(1))/2.d0/dX
fx(lx-1) = (f(lx)-f(lx-2))/2.d0/dX
fx(1) = (f(2)-f(1))/dX
fx(lx) = (f(lx)-f(lx-1))/dX
! fx(1) =-(2.5d1*f(1)-4.8d1*f(2)+3.6d1*f(3)-1.6d1*f(4)+3.d0*f(5))/1.2d1/dX
! fx(2) =-(3.d0*f(1)+1.d1*f(2)-1.8d1*f(3)+6.d0*f(4)-f(5))/1.2d1/dX
#elif (ghost_width == 4)
fx(4:lx-3) = (-f(1:lx-6)+9.d0*f(2:lx-5)-4.5d1*f(3:lx-4)+4.5d1*f(5:lx-2)-9.d0*f(6:lx-1)+f(7:lx))/6.d1/dX
fx(3) = (f(1)-8.d0*f(2)+8.d0*f(4)-f(5))/1.2d1/dX
fx(lx-2) = (f(lx-4)-8.d0*f(lx-3)+8.d0*f(lx-1)-f(lx))/1.2d1/dX
fx(2) = (f(3)-f(1))/2.d0/dX
fx(lx-1) = (f(lx)-f(lx-2))/2.d0/dX
fx(1) = (f(2)-f(1))/dX
fx(lx) = (f(lx)-f(lx-1))/dX
#elif (ghost_width == 5)
fx(5:lx-4) = (3.d0*f(1:lx-8)-3.2d1*f(2:lx-7)+1.68d2*f(3:lx-6)-6.72d2*f(4:lx-5)+ &
6.72d2*f(6:lx-3)-1.68d2*f(7:lx-2)+3.2d1*f(8:lx-1)-3.d0*f(9:lx))/8.4d2/dX
fx(4) = (-f(1)+9.d0*f(2)-4.5d1*f(3)+4.5d1*f(5)-9.d0*f(6)+f(7))/6.d1/dX
fx(lx-3) = (-f(lx-6)+9.d0*f(lx-5)-4.5d1*f(lx-4)+4.5d1*f(lx-2)-9.d0*f(lx-1)+f(lx))/6.d1/dX
fx(3) = (f(1)-8.d0*f(2)+8.d0*f(4)-f(5))/1.2d1/dX
fx(lx-2) = (f(lx-4)-8.d0*f(lx-3)+8.d0*f(lx-1)-f(lx))/1.2d1/dX
fx(2) = (f(3)-f(1))/2.d0/dX
fx(lx-1) = (f(lx)-f(lx-2))/2.d0/dX
fx(1) = (f(2)-f(1))/dX
fx(lx) = (f(lx)-f(lx-1))/dX
#endif
#endif
return
end subroutine rderivs_x
!----
subroutine rderivs_x_point(lx,X,f,fx,k)
implicit none
!~~~~~~% Input parameters:
integer,intent(in) :: lx,k
real*8,intent(in),dimension(lx) :: X
real*8,intent(in),dimension(lx) :: f
real*8,intent(out) :: fx
real*8 :: dX
dX = X(2)-X(1)
#ifdef OLD
if(k .eq. lx)then
fx = (f(lx)-f(lx-1))/dX
else
fx = (f(k+1)-f(k))/dX
endif
#else
#if (ghost_width == 2)
if(k .gt. 1 .and. k .lt. lx) then
fx = (f(k+1)-f(k-1))/2.d0/dX
elseif(k.eq.1) then
fx = (f(2)-f(1))/dX
elseif(k.eq.lx) then
fx = (f(lx)-f(lx-1))/dX
endif
#elif (ghost_width == 3)
if(k .gt. 2 .and. k .lt. lx-1) then
fx = (f(k-2)-8.d0*f(k-1)+8.d0*f(k+1)-f(k+2))/1.2d1/dX
elseif(k.eq.1) then
fx = (f(2)-f(1))/dX
elseif(k.eq.lx) then
fx = (f(lx)-f(lx-1))/dX
elseif(k.eq.2) then
fx = (f(3)-f(1))/2.d0/dX
elseif(k.eq.lx-1) then
fx = (f(lx)-f(lx-2))/2.d0/dX
endif
#elif (ghost_width == 4)
if(k .gt. 3 .and. k .lt. lx-2) then
fx = (-f(k-3)+9.d0*f(k-2)-4.5d1*f(k-1)+4.5d1*f(k+1)-9.d0*f(k+2)+f(k+3))/6.d1/dX
elseif(k.eq.1) then
fx = (f(2)-f(1))/dX
elseif(k.eq.lx) then
fx = (f(lx)-f(lx-1))/dX
elseif(k.eq.2) then
fx = (f(3)-f(1))/2.d0/dX
elseif(k.eq.lx-1) then
fx = (f(lx)-f(lx-2))/2.d0/dX
elseif(k.eq.3) then
fx = (f(1)-8.d0*f(2)+8.d0*f(4)-f(5))/1.2d1/dX
elseif(k.eq.lx-2) then
fx = (f(lx-4)-8.d0*f(lx-3)+8.d0*f(lx-1)-f(lx))/1.2d1/dX
endif
#elif (ghost_width == 5)
if(k .gt. 4 .and. k .lt. lx-3) then
fx = (3.d0*f(k-4)-3.2d1*f(k-3)+1.68d2*f(k-2)-6.72d2*f(k-1)+ &
6.72d2*f(k+1)-1.68d2*f(k+2)+3.2d1*f(k+3)-3.d0*f(k+4))/8.4d2/dX
elseif(k.eq.1) then
fx = (f(2)-f(1))/dX
elseif(k.eq.lx) then
fx = (f(lx)-f(lx-1))/dX
elseif(k.eq.2) then
fx = (f(3)-f(1))/2.d0/dX
elseif(k.eq.lx-1) then
fx = (f(lx)-f(lx-2))/2.d0/dX
elseif(k.eq.3) then
fx = (f(1)-8.d0*f(2)+8.d0*f(4)-f(5))/1.2d1/dX
elseif(k.eq.lx-2) then
fx = (f(lx-4)-8.d0*f(lx-3)+8.d0*f(lx-1)-f(lx))/1.2d1/dX
elseif(k.eq.4) then
fx = (-f(1)+9.d0*f(2)-4.5d1*f(3)+4.5d1*f(5)-9.d0*f(6)+f(7))/6.d1/dX
elseif(k.eq.lx-3) then
fx = (-f(lx-6)+9.d0*f(lx-5)-4.5d1*f(lx-4)+4.5d1*f(lx-2)-9.d0*f(lx-1)+f(lx))/6.d1/dX
endif
#endif
#endif
return
end subroutine rderivs_x_point
!----
subroutine cderivs_x(lx,X,f,fx)
implicit none
!~~~~~~% Input parameters:
integer,intent(in) :: lx
real*8,intent(in),dimension(lx) :: X
double complex,intent(in),dimension(lx) :: f
double complex,intent(out),dimension(lx) :: fx
real*8 :: dX
dX = X(2)-X(1)
#ifdef OLD
fx(1:lx-1) = (f(2:lx)-f(1:lx-1))/dX
fx(lx) = (f(lx)-f(lx-1))/dX
#else
#if (ghost_width == 2)
fx(2:lx-1) = (f(3:lx)-f(1:lx-2))/2.d0/dX
fx(1) = (f(2)-f(1))/dX
fx(lx) = (f(lx)-f(lx-1))/dX
#elif (ghost_width == 3)
fx(3:lx-2) = (f(1:lx-4)-8.d0*f(2:lx-3)+8.d0*f(4:lx-1)-f(5:lx))/1.2d1/dX
fx(2) = (f(3)-f(1))/2.d0/dX
fx(lx-1) = (f(lx)-f(lx-2))/2.d0/dX
fx(1) = (f(2)-f(1))/dX
fx(lx) = (f(lx)-f(lx-1))/dX
! fx(1) =-(2.5d1*f(1)-4.8d1*f(2)+3.6d1*f(3)-1.6d1*f(4)+3.d0*f(5))/1.2d1/dX
! fx(2) =-(3.d0*f(1)+1.d1*f(2)-1.8d1*f(3)+6.d0*f(4)-f(5))/1.2d1/dX
#elif (ghost_width == 4)
fx(4:lx-3) = (-f(1:lx-6)+9.d0*f(2:lx-5)-4.5d1*f(3:lx-4)+4.5d1*f(5:lx-2)-9.d0*f(6:lx-1)+f(7:lx))/6.d1/dX
fx(3) = (f(1)-8.d0*f(2)+8.d0*f(4)-f(5))/1.2d1/dX
fx(lx-2) = (f(lx-4)-8.d0*f(lx-3)+8.d0*f(lx-1)-f(lx))/1.2d1/dX
fx(2) = (f(3)-f(1))/2.d0/dX
fx(lx-1) = (f(lx)-f(lx-2))/2.d0/dX
fx(1) = (f(2)-f(1))/dX
fx(lx) = (f(lx)-f(lx-1))/dX
#elif (ghost_width == 5)
fx(5:lx-4) = (3.d0*f(1:lx-8)-3.2d1*f(2:lx-7)+1.68d2*f(3:lx-6)-6.72d2*f(4:lx-5)+ &
6.72d2*f(6:lx-3)-1.68d2*f(7:lx-2)+3.2d1*f(8:lx-1)-3.d0*f(9:lx))/8.4d2/dX
fx(4) = (-f(1)+9.d0*f(2)-4.5d1*f(3)+4.5d1*f(5)-9.d0*f(6)+f(7))/6.d1/dX
fx(lx-3) = (-f(lx-6)+9.d0*f(lx-5)-4.5d1*f(lx-4)+4.5d1*f(lx-2)-9.d0*f(lx-1)+f(lx))/6.d1/dX
fx(3) = (f(1)-8.d0*f(2)+8.d0*f(4)-f(5))/1.2d1/dX
fx(lx-2) = (f(lx-4)-8.d0*f(lx-3)+8.d0*f(lx-1)-f(lx))/1.2d1/dX
fx(2) = (f(3)-f(1))/2.d0/dX
fx(lx-1) = (f(lx)-f(lx-2))/2.d0/dX
fx(1) = (f(2)-f(1))/dX
fx(lx) = (f(lx)-f(lx-1))/dX
#endif
#endif
return
end subroutine cderivs_x
!----
subroutine cdderivs_x(lx,X,f,fxx)
implicit none
!~~~~~~% Input parameters:
integer,intent(in) :: lx
real*8,intent(in),dimension(lx) :: X
double complex,intent(in),dimension(lx) :: f
double complex,intent(out),dimension(lx) :: fxx
real*8 :: dX
dX = X(2)-X(1)
dX = dX*dX
#ifdef OLD
fxx(1:lx-2) = (f(3:lx)-2.0*f(2:lx-1)+f(1:lx-2))/dX
fxx(lx-1) = (f(lx)-2.0*f(lx-1)+f(lx-2))/dX
fxx(lx ) = (f(lx)-2.0*f(lx-1)+f(lx-2))/dX
#else
#if (ghost_width == 2)
fxx(2:lx-1) = (f(3:lx)-2.d0*f(2:lx-1)+f(1:lx-2))/dX
fxx(1) = (f(3)-2.d0*f(2)+f(1))/dX
fxx(lx) = (f(lx)-2.d0*f(lx-1)+f(lx-2))/dX
#elif (ghost_width == 3)
fxx(3:lx-2) = (-f(1:lx-4)+1.6d1*f(2:lx-3)-3.d1*f(3:lx-2)+1.6d1*f(4:lx-1)-f(5:lx))/1.2d1/dX
fxx(2) = (f(3)-2.d0*f(2)+f(1))/dX
fxx(lx-1) = (f(lx)-2.d0*f(lx-1)+f(lx-2))/dX
fxx(1) = (f(3)-2.d0*f(2)+f(1))/dX
fxx(lx) = (f(lx)-2.d0*f(lx-1)+f(lx-2))/dX
#elif (ghost_width == 4)
fxx(4:lx-3) = (2.d0*f(1:lx-6)-2.7d1*f(2:lx-5)+2.7d2*f(3:lx-4)-4.9d2*f(4:lx-3) &
+2.7d2*f(5:lx-2)-2.7d1*f(6:lx-1)+2.d0*f(7:lx))/1.8d2/dX
fxx(3) = (-f(1)+1.6d1*f(2)-3.d1*f(3)+1.6d1*f(4)-f(5))/1.2d1/dX
fxx(lx-2) = (-f(lx-4)+1.6d1*f(lx-3)-3.d1*f(lx-2)+1.6d1*f(lx-1)-f(lx))/1.2d1/dX
fxx(2) = (f(3)-2.d0*f(2)+f(1))/dX
fxx(lx-1) = (f(lx)-2.d0*f(lx-1)+f(lx-2))/dX
fxx(1) = (f(3)-2.d0*f(2)+f(1))/dX
fxx(lx) = (f(lx)-2.d0*f(lx-1)+f(lx-2))/dX
#elif (ghost_width == 5)
fxx(5:lx-4) = (-9.d0*f(1:lx-8)+1.28d2*f(2:lx-7)-1.008d3*f(3:lx-6)+8.064d3*f(4:lx-5)-1.435d4*f(5:lx-4) &
+8.064d3*f(6:lx-3)-1.008d3*f(7:lx-2)+1.28d2*f(8:lx-1)-9.d0*f(9:lx))/5.04d3/dX
fxx(4) = (2.d0*f(1)-2.7d1*f(2)+2.7d2*f(3)-4.9d2*f(4) &
+2.7d2*f(5)-2.7d1*f(6)+2.d0*f(7))/1.8d2/dX
fxx(lx-3) = (2.d0*f(lx-6)-2.7d1*f(lx-5)+2.7d2*f(lx-4)-4.9d2*f(lx-3) &
+2.7d2*f(lx-2)-2.7d1*f(lx-1)+2.d0*f(lx))/1.8d2/dX
fxx(3) = (-f(1)+1.6d1*f(2)-3.d1*f(3)+1.6d1*f(4)-f(5))/1.2d1/dX
fxx(lx-2) = (-f(lx-4)+1.6d1*f(lx-3)-3.d1*f(lx-2)+1.6d1*f(lx-1)-f(lx))/1.2d1/dX
fxx(2) = (f(3)-2.d0*f(2)+f(1))/dX
fxx(lx-1) = (f(lx)-2.d0*f(lx-1)+f(lx-2))/dX
fxx(1) = (f(3)-2.d0*f(2)+f(1))/dX
fxx(lx) = (f(lx)-2.d0*f(lx-1)+f(lx-2))/dX
#endif
#endif
return
end subroutine cdderivs_x
!----
subroutine rdderivs_x(lx,X,f,fxx)
implicit none
!~~~~~~% Input parameters:
integer,intent(in) :: lx
real*8,intent(in),dimension(lx) :: X
real*8,intent(in),dimension(lx) :: f
real*8,intent(out),dimension(lx) :: fxx
real*8 :: dX
dX = X(2)-X(1)
dX = dX*dX
#ifdef OLD
fxx(1:lx-2) = (f(3:lx)-2.0*f(2:lx-1)+f(1:lx-2))/dX
fxx(lx-1) = (f(lx)-2.0*f(lx-1)+f(lx-2))/dX
fxx(lx ) = (f(lx)-2.0*f(lx-1)+f(lx-2))/dX
#else
#if (ghost_width == 2)
fxx(2:lx-1) = (f(3:lx)-2.d0*f(2:lx-1)+f(1:lx-2))/dX
fxx(1) = (f(3)-2.d0*f(2)+f(1))/dX
fxx(lx) = (f(lx)-2.d0*f(lx-1)+f(lx-2))/dX
#elif (ghost_width == 3)
fxx(3:lx-2) = (-f(1:lx-4)+1.6d1*f(2:lx-3)-3.d1*f(3:lx-2)+1.6d1*f(4:lx-1)-f(5:lx))/1.2d1/dX
fxx(2) = (f(3)-2.d0*f(2)+f(1))/dX
fxx(lx-1) = (f(lx)-2.d0*f(lx-1)+f(lx-2))/dX
fxx(1) = (f(3)-2.d0*f(2)+f(1))/dX
fxx(lx) = (f(lx)-2.d0*f(lx-1)+f(lx-2))/dX
#elif (ghost_width == 4)
fxx(4:lx-3) = (2.d0*f(1:lx-6)-2.7d1*f(2:lx-5)+2.7d2*f(3:lx-4)-4.9d2*f(4:lx-3) &
+2.7d2*f(5:lx-2)-2.7d1*f(6:lx-1)+2.d0*f(7:lx))/1.8d2/dX
fxx(3) = (-f(1)+1.6d1*f(2)-3.d1*f(3)+1.6d1*f(4)-f(5))/1.2d1/dX
fxx(lx-2) = (-f(lx-4)+1.6d1*f(lx-3)-3.d1*f(lx-2)+1.6d1*f(lx-1)-f(lx))/1.2d1/dX
fxx(2) = (f(3)-2.d0*f(2)+f(1))/dX
fxx(lx-1) = (f(lx)-2.d0*f(lx-1)+f(lx-2))/dX
fxx(1) = (f(3)-2.d0*f(2)+f(1))/dX
fxx(lx) = (f(lx)-2.d0*f(lx-1)+f(lx-2))/dX
#elif (ghost_width == 5)
fxx(5:lx-4) = (-9.d0*f(1:lx-8)+1.28d2*f(2:lx-7)-1.008d3*f(3:lx-6)+8.064d3*f(4:lx-5)-1.435d4*f(5:lx-4) &
+8.064d3*f(6:lx-3)-1.008d3*f(7:lx-2)+1.28d2*f(8:lx-1)-9.d0*f(9:lx))/5.04d3/dX
fxx(4) = (2.d0*f(1)-2.7d1*f(2)+2.7d2*f(3)-4.9d2*f(4) &
+2.7d2*f(5)-2.7d1*f(6)+2.d0*f(7))/1.8d2/dX
fxx(lx-3) = (2.d0*f(lx-6)-2.7d1*f(lx-5)+2.7d2*f(lx-4)-4.9d2*f(lx-3) &
+2.7d2*f(lx-2)-2.7d1*f(lx-1)+2.d0*f(lx))/1.8d2/dX
fxx(3) = (-f(1)+1.6d1*f(2)-3.d1*f(3)+1.6d1*f(4)-f(5))/1.2d1/dX
fxx(lx-2) = (-f(lx-4)+1.6d1*f(lx-3)-3.d1*f(lx-2)+1.6d1*f(lx-1)-f(lx))/1.2d1/dX
fxx(2) = (f(3)-2.d0*f(2)+f(1))/dX
fxx(lx-1) = (f(lx)-2.d0*f(lx-1)+f(lx-2))/dX
fxx(1) = (f(3)-2.d0*f(2)+f(1))/dX
fxx(lx) = (f(lx)-2.d0*f(lx-1)+f(lx-2))/dX
#endif
#endif
return
end subroutine rdderivs_x
!----
subroutine rdderivs_x_point(lx,X,f,fxx,k)
implicit none
!~~~~~~% Input parameters:
integer,intent(in) :: lx,k
real*8,intent(in),dimension(lx) :: X
real*8,intent(in),dimension(lx) :: f
real*8,intent(out) :: fxx
real*8 :: dX
dX = X(2)-X(1)
dX = dX*dX
#ifdef OLD
if(k.lt.lx-1) then
fxx = (f(k+2)-2.0*f(k+1)+f(k))/dX
elseif(k.eq.lx-1) then
fxx = (f(lx)-2.0*f(lx-1)+f(lx-2))/dX
elseif(k.eq.lx) then
fxx = (f(lx)-2.0*f(lx-1)+f(lx-2))/dX
endif
#else
#if (ghost_width == 2)
if(k.gt.1 .and. k.lt.lx) then
fxx = (f(k+1)-2.d0*f(k)+f(k-1))/dX
elseif(k.eq.1) then
fxx = (f(3)-2.d0*f(2)+f(1))/dX
elseif(k.eq.lx) then
fxx = (f(lx)-2.d0*f(lx-1)+f(lx-2))/dX
endif
#elif (ghost_width == 3)
if(k.gt.2 .and. k.lt.lx-1) then
fxx = (-f(k-2)+1.6d1*f(k-1)-3.d1*f(k)+1.6d1*f(k+1)-f(k+2))/1.2d1/dX
elseif(k.eq.1) then
fxx = (f(3)-2.d0*f(2)+f(1))/dX
elseif(k.eq.lx) then
fxx = (f(lx)-2.d0*f(lx-1)+f(lx-2))/dX
elseif(k.eq.2) then
fxx = (f(3)-2.d0*f(2)+f(1))/dX
elseif(k.eq.lx-1) then
fxx = (f(lx)-2.d0*f(lx-1)+f(lx-2))/dX
endif
#elif (ghost_width == 4)
if(k.gt.3 .and. k.lt.lx-2)then
fxx = (2.d0*f(k-3)-2.7d1*f(k-2)+2.7d2*f(k-1)-4.9d2*f(k) &
+2.7d2*f(k+1)-2.7d1*f(k+2)+2.d0*f(k+3))/1.8d2/dX
elseif(k.eq.1) then
fxx = (f(3)-2.d0*f(2)+f(1))/dX
elseif(k.eq.lx) then
fxx = (f(lx)-2.d0*f(lx-1)+f(lx-2))/dX
elseif(k.eq.2) then
fxx = (f(3)-2.d0*f(2)+f(1))/dX
elseif(k.eq.lx-1) then
fxx = (f(lx)-2.d0*f(lx-1)+f(lx-2))/dX
elseif(k.eq.3) then
fxx = (-f(1)+1.6d1*f(2)-3.d1*f(3)+1.6d1*f(4)-f(5))/1.2d1/dX
elseif(k.eq.lx-2) then
fxx = (-f(lx-4)+1.6d1*f(lx-3)-3.d1*f(lx-2)+1.6d1*f(lx-1)-f(lx))/1.2d1/dX
endif
#elif (ghost_width == 5)
if(k.gt.4 .and. k.lt.lx-3) then
fxx = (-9.d0*f(k-4)+1.28d2*f(k-3)-1.008d3*f(k-2)+8.064d3*f(k-1)-1.435d4*f(k) &
+8.064d3*f(k+1)-1.008d3*f(k+2)+1.28d2*f(k+3)-9.d0*f(k+4))/5.04d3/dX
elseif(k.eq.1) then
fxx = (f(3)-2.d0*f(2)+f(1))/dX
elseif(k.eq.lx) then
fxx = (f(lx)-2.d0*f(lx-1)+f(lx-2))/dX
elseif(k.eq.2) then
fxx = (f(3)-2.d0*f(2)+f(1))/dX
elseif(k.eq.lx-1) then
fxx = (f(lx)-2.d0*f(lx-1)+f(lx-2))/dX
elseif(k.eq.3) then
fxx = (-f(1)+1.6d1*f(2)-3.d1*f(3)+1.6d1*f(4)-f(5))/1.2d1/dX
elseif(k.eq.lx-2) then
fxx = (-f(lx-4)+1.6d1*f(lx-3)-3.d1*f(lx-2)+1.6d1*f(lx-1)-f(lx))/1.2d1/dX
elseif(k.eq.4) then
fxx = (2.d0*f(1)-2.7d1*f(2)+2.7d2*f(3)-4.9d2*f(4) &
+2.7d2*f(5)-2.7d1*f(6)+2.d0*f(7))/1.8d2/dX
elseif(k.eq.lx-3) then
fxx = (2.d0*f(lx-6)-2.7d1*f(lx-5)+2.7d2*f(lx-4)-4.9d2*f(lx-3) &
+2.7d2*f(lx-2)-2.7d1*f(lx-1)+2.d0*f(lx))/1.8d2/dX
endif
#endif
#endif
return
end subroutine rdderivs_x_point
!----
subroutine rdderivs_xy_point(lx,ly,X,Y,f,fxy,i,j)
implicit none
!~~~~~~% Input parameters:
integer,intent(in) :: lx,ly,i,j
real*8,intent(in),dimension(lx) :: X
real*8,intent(in),dimension(ly) :: Y
real*8,intent(in),dimension(lx,ly) :: f
real*8,intent(out) :: fxy
real*8 :: dX,dY
dX = X(2)-X(1)
dY = Y(2)-Y(1)
!! we only consider inner points
#if (ghost_width == 2)
if(i>1 .and. j>1.and.i<lx .and.j<ly)then
fxy = (f(i-1,j-1)-f(i+1,j-1)-f(i-1,j+1)+f(i+1,j+1))/(4*dX*dY)
else
fxy = 0.d0
endif
#elif (ghost_width == 3)
if(i>2 .and. j>2.and.i<lx-1 .and.j<ly-1)then
fxy = ( (f(i-2,j-2)-8*f(i-1,j-2)+8*f(i+1,j-2)-f(i+2,j-2)) &
-8 *(f(i-2,j-1)-8*f(i-1,j-1)+8*f(i+1,j-1)-f(i+2,j-1)) &
+8 *(f(i-2,j+1)-8*f(i-1,j+1)+8*f(i+1,j+1)-f(i+2,j+1)) &
- (f(i-2,j+2)-8*f(i-1,j+2)+8*f(i+1,j+2)-f(i+2,j+2)))/(144*dX*dY)
else
fxy = 0.d0
endif
#elif (ghost_width == 4)
if(i>3 .and. j>3.and.i<lx-2 .and.j<ly-2)then
fxy = (- (-f(i-3,j-3)+9*f(i-2,j-3)-45*f(i-1,j-3)+45*f(i+1,j-3)-9*f(i+2,j-3)+f(i+3,j-3)) &
+9 *(-f(i-3,j-2)+9*f(i-2,j-2)-45*f(i-1,j-2)+45*f(i+1,j-2)-9*f(i+2,j-2)+f(i+3,j-2)) &
-45*(-f(i-3,j-1)+9*f(i-2,j-1)-45*f(i-1,j-1)+45*f(i+1,j-1)-9*f(i+2,j-1)+f(i+3,j-1)) &
+45*(-f(i-3,j+1)+9*f(i-2,j+1)-45*f(i-1,j+1)+45*f(i+1,j+1)-9*f(i+2,j+1)+f(i+3,j+1)) &
-9 *(-f(i-3,j+2)+9*f(i-2,j+2)-45*f(i-1,j+2)+45*f(i+1,j+2)-9*f(i+2,j+2)+f(i+3,j+2)) &
+ (-f(i-3,j+3)+9*f(i-2,j+3)-45*f(i-1,j+3)+45*f(i+1,j+3)-9*f(i+2,j+3)+f(i+3,j+3)))/(3.6d3*dX*dY)
else
fxy = 0.d0
endif
#elif (ghost_width == 5)
if(i>4 .and. j>4.and.i<lx-3 .and.j<ly-3)then
fxy = ( 3 *( 3*f(i-4,j-4)-32*f(i-3,j-4)+168*f(i-2,j-4)-672*f(i-1,j-4) &
-3*f(i+4,j-4)+32*f(i+3,j-4)-168*f(i+2,j-4)+672*f(i+1,j-4)) &
-32 *( 3*f(i-4,j-3)-32*f(i-3,j-3)+168*f(i-2,j-3)-672*f(i-1,j-3) &
-3*f(i+4,j-3)+32*f(i+3,j-3)-168*f(i+2,j-3)+672*f(i+1,j-3)) &
+168*( 3*f(i-4,j-2)-32*f(i-3,j-2)+168*f(i-2,j-2)-672*f(i-1,j-2) &
-3*f(i+4,j-2)+32*f(i+3,j-2)-168*f(i+2,j-2)+672*f(i+1,j-2)) &
-672*( 3*f(i-4,j-1)-32*f(i-3,j-1)+168*f(i-2,j-1)-672*f(i-1,j-1) &
-3*f(i+4,j-1)+32*f(i+3,j-1)-168*f(i+2,j-1)+672*f(i+1,j-1)) &
+672*( 3*f(i-4,j+1)-32*f(i-3,j+1)+168*f(i-2,j+1)-672*f(i-1,j+1) &
-3*f(i+4,j+1)+32*f(i+3,j+1)-168*f(i+2,j+1)+672*f(i+1,j+1)) &
-168*( 3*f(i-4,j+2)-32*f(i-3,j+2)+168*f(i-2,j+2)-672*f(i-1,j+2) &
-3*f(i+4,j+2)+32*f(i+3,j+2)-168*f(i+2,j+2)+672*f(i+1,j+2)) &
+32 *( 3*f(i-4,j+3)-32*f(i-3,j+3)+168*f(i-2,j+3)-672*f(i-1,j+3) &
-3*f(i+4,j+3)+32*f(i+3,j+3)-168*f(i+2,j+3)+672*f(i+1,j+3)) &
-3 *( 3*f(i-4,j+4)-32*f(i-3,j+4)+168*f(i-2,j+4)-672*f(i-1,j+4) &
-3*f(i+4,j+4)+32*f(i+3,j+4)-168*f(i+2,j+4)+672*f(i+1,j+4)) )/(7.056d5*dX*dY)
else
fxy = 0.d0
endif
#endif
return
end subroutine rdderivs_xy_point
! side-winded type finite difference (arXiv:1208.3891)
!====================================================================================
subroutine rderivs_sw_x(lx,X,f,fx)
implicit none
!~~~~~~% Input parameters:
integer,intent(in) :: lx
real*8,intent(in),dimension(lx) :: X
real*8,intent(in),dimension(lx) :: f
real*8,intent(out),dimension(lx) :: fx
real*8 :: dX
integer :: i
dX = X(2)-X(1)
#if (ghost_width == 2)
!#error
do i=1,lx-3
fx(i) = (-3.d0*f(i)+4.d0*f(i+1)-f(i+2))/2.d0/dX ! forward difference, 2nd-order accuracy
enddo
do i=lx-2,lx-1
fx(i) = (f(i+1)-f(i-1))/2.d0/dX ! centered difference, 2nd-order accuracy
enddo
i=lx
fx(i) = (3.d0*f(i)-4.d0*f(i-1)+f(i-2))/2.d0/dX ! backward difference, 2nd-order accuracy
#elif (ghost_width == 3)
do i=1,lx-5
fx(i) = (-2.5d1/1.2d1*f(i)+4*f(i+1)-3*f(i+2)+4.d0/3.d0*f(i+3)-0.25d0*f(i+4))/dX ! forward difference, 4th-order accuracy
enddo
do i=lx-4,lx-2
fx(i) = (f(i-2)-8.d0*f(i-1)+8.d0*f(i+1)-f(i+2))/1.2d1/dX ! centered difference, 4th-order accuracy
enddo
do i=lx-1,lx
fx(i) = (2.5d1/1.2d1*f(i)-4*f(i-1)+3*f(i-2)-4.d0/3.d0*f(i-3)+0.25d0*f(i-4))/dX ! backward difference, 4th-order accuracy
enddo
#elif (ghost_width == 4)
!#error
do i=1,lx-7
fx(i) = (-147.d0*f(i)+360.d0*f(i+1)-450.d0*f(i+2)+400.d0*f(i+3)-225.d0*f(i+4)+72.d0*f(i+5)-10.d0*f(i+6)) /60.d0/dX ! forward difference, 6th-order accuracy
enddo
do i=lx-6,lx-3
fx(i) = (f(i-3)-9.d0*f(i-2)+45.d0*f(i-1)-45.d0*f(i+1)+9.0*f(i+2)-f(i+3)) /60.d0/dX ! centered difference, 6th-order accuracy
enddo
do i=lx-2,lx
fx(i) = (147.d0*f(i)-360.d0*f(i-1)+450.d0*f(i-2)-400.d0*f(i-3)+225.d0*f(i-4)-72.d0*f(i-5)+10.d0*f(i-6)) /60.d0/dX ! backward difference, 6th-order accuracy
enddo
#elif (ghost_width == 5)
!#error
do i=1,lx-9
fx(i) = (-(761.d0/280.d0)*f(i)+8.d0*f(i+1)-14.d0*f(i+2)+(56.d0/3.d0)*f(i+3)-(35.d0/2.d0)*f(i+4) &
+(56.d0/5.d0)*f(i+5)-(14.d0/3.d0)*f(i+6)+(8.d0/7.d0)*f(i+7)-(1.d0/8.d0)*f(i+8)) /dX
! forward difference, 8th-order accuracy
enddo
do i=lx-8,lx-4
fx(i) = (3.d0*f(i-4)-32.d0*f(i-3)+168.d0*f(i-2)-672.d0*f(i-1)+672.d0*f(i+1)-168.d0*f(i+2)+32.d0*f(i+3)-3.d0*f(i+4)) /840.d0/dX
! centered difference, 8th-order accuracy
enddo
do i=lx-3,lx
fx(i) = ((761.d0/280.d0)*f(i)-8.d0*f(i-1)+14.d0*f(i-2)-(56.d0/3.d0)*f(i-3)+(35.d0/2.d0)*f(i-4) &
-(56.d0/5.d0)*f(i-5)+(14.d0/3.d0)*f(i-6)-(8.d0/7.d0)*f(i-7)+(1.d0/8.d0)*f(i-8)) /dX
! backward difference, 8th-order accuracy
enddo
#endif
return
end subroutine rderivs_sw_x
!----
subroutine cderivs_sw_x(lx,X,f,fx)
implicit none
!~~~~~~% Input parameters:
integer,intent(in) :: lx
real*8,intent(in),dimension(lx) :: X
double complex,intent(in),dimension(lx) :: f
double complex,intent(out),dimension(lx) :: fx
real*8 :: dX
integer :: i
dX = X(2)-X(1)
#if (ghost_width == 2)
!#error
do i=1,lx-3
fx(i) = (-3.d0*f(i)+4.d0*f(i+1)-f(i+2))/2.d0/dX ! forward difference, 2nd-order accuracy
enddo
do i=lx-2,lx-1 ! width = 2
fx(i) = (f(i+1)-f(i-1))/2.d0/dX ! centered difference, 2nd-order accuracy
enddo
i=lx
fx(i) = (3.d0*f(i)-4.d0*f(i-1)+f(i-2))/2.d0/dX ! backward difference, 2nd-order accuracy
#elif (ghost_width == 3)
do i=1,lx-5
fx(i) = (-2.5d1/1.2d1*f(i)+4*f(i+1)-3*f(i+2)+4.d0/3.d0*f(i+3)-0.25d0*f(i+4))/dX ! forward difference, 4th-order accuracy
enddo
do i=lx-4,lx-2 ! width = 3
fx(i) = (f(i-2)-8.d0*f(i-1)+8.d0*f(i+1)-f(i+2))/1.2d1/dX ! centered difference, 4th-order accuracy
enddo
do i=lx-1,lx
fx(i) = (2.5d1/1.2d1*f(i)-4*f(i-1)+3*f(i-2)-4.d0/3.d0*f(i-3)+0.25d0*f(i-4))/dX ! backward difference, 4th-order accuracy
enddo
#elif (ghost_width == 4)
!#error
do i=1,lx-7
fx(i) = (-147.d0*f(i)+360.d0*f(i+1)-450.d0*f(i+2)+400.d0*f(i+3)-225.d0*f(i+4)+72.d0*f(i+5)-10.d0*f(i+6)) /60.d0/dX
! forward difference, 6th-order accuracy
enddo
do i=lx-6,lx-3 ! width = 4
fx(i) = (f(i-3)-9.d0*f(i-2)+45.d0*f(i-1)-45.d0*f(i+1)+9.0*f(i+2)-f(i+3)) /60.d0/dX ! centered difference, 6th-order accuracy
enddo
do i=lx-2,lx
fx(i) = (147.d0*f(i)-360.d0*f(i-1)+450.d0*f(i-2)-400.d0*f(i-3)+225.d0*f(i-4)-72.d0*f(i-5)+10.d0*f(i-6)) /60.d0/dX
! backward difference, 6th-order accuracy
enddo
#elif (ghost_width == 5)
!#error
do i=1,lx-9
fx(i) = (-(761.d0/280.d0)*f(i)+8.d0*f(i+1)-14.d0*f(i+2)+(56.d0/3.d0)*f(i+3)-(35.d0/2.d0)*f(i+4) &
+(56.d0/5.d0)*f(i+5)-(14.d0/3.d0)*f(i+6)+(8.d0/7.d0)*f(i+7)-(1.d0/8.d0)*f(i+8)) /dX
! forward difference, 8th-order accuracy
enddo
do i=lx-8,lx-4 ! width = 5
fx(i) = (3.d0*f(i-4)-32.d0*f(i-3)+168.d0*f(i-2)-672.d0*f(i-1)+672.d0*f(i+1)-168.d0*f(i+2)+32.d0*f(i+3)-3.d0*f(i+4)) /840.d0/dX
! centered difference, 8th-order accuracy
enddo
do i=lx-3,lx
fx(i) = ((761.d0/280.d0)*f(i)-8.d0*f(i-1)+14.d0*f(i-2)-(56.d0/3.d0)*f(i-3)+(35.d0/2.d0)*f(i-4) &
-(56.d0/5.d0)*f(i-5)+(14.d0/3.d0)*f(i-6)-(8.d0/7.d0)*f(i-7)+(1.d0/8.d0)*f(i-8)) /dX
! backward difference, 8th-order accuracy
enddo
#endif
return
end subroutine cderivs_sw_x
!----
subroutine cdderivs_sw_x(lx,X,f,fxx)
implicit none
!~~~~~~% Input parameters:
integer,intent(in) :: lx
real*8,intent(in),dimension(lx) :: X
double complex,intent(in),dimension(lx) :: f
double complex,intent(out),dimension(lx) :: fxx
real*8 :: dX
integer :: i
dX = X(2)-X(1)
dX = dX*dX
#if (ghost_width == 2)
!#error
do i=1,lx-3
fxx(i) = (f(i)-2.d0*f(i+1)+f(i+2))/dX ! forward difference, 2nd-order accuracy
enddo
do i=lx-2,lx-1
fxx(i) = (f(i-1)-2.d0*f(i)+f(i+1))/dX ! centered difference, 2nd-order accuracy
enddo
i=lx
fxx(i) = (f(i)-2.d0*f(i-1)+f(i-2))/dX ! backward difference, 2nd-order accuracy
#elif (ghost_width == 3)
do i=1,lx-5
fxx(i) = ((15.d0/4)*f(i)-(77.d0/6)*f(i+1)+(107.d0/6)*f(i+2)-13.d0*f(i+3)+(61.d0/12.d0)*f(i+4)-(5.d0/6.d0)*f(i+5))/dX ! forward difference, 4th-order accuracy
enddo
do i=lx-4,lx-2 ! width = 3
fxx(i) = (-f(i-2)+16.d0*f(i-1)-30.d0*f(i)+16.d0*f(i+1)-f(i+2))/12.d0/dX ! centered difference, 4th-order accuracy
enddo
do i=lx-1,lx
fxx(i) = ((15.d0/4)*f(i)-(77.d0/6)*f(i-1)+(107.d0/6)*f(i-2)-13.d0*f(i-3)+(61.d0/12.d0)*f(i-4)-(5.d0/6.d0)*f(i-5))/dX ! backward difference, 4th-order accuracy
enddo
#elif (ghost_width == 4)
do i=1,lx-7
fxx(i) = (812.d0*f(i)-3132.d0*f(i+1)+5265.d0*f(i+2)-5080.d0*f(i+3)+2970.d0*f(i+4)-972.d0*f(i+5)+137.d0*f(i+6)) /180.d0/dX ! forward difference, 6th-order accuracy
enddo
do i=lx-6,lx-3 ! width = 4
fxx(i) = (2.d0*f(i-3)-27.d0*f(i-2)+270.d0*f(i-1)-490.d0*f(i)+270.d0*f(i+1)-27.d0*f(i+2)+2.d0*f(i+3)) /180.d0/dX ! centered difference, 6th-order accuracy
enddo
do i=lx-2,lx
fxx(i) = (812.d0*f(i)-3132.d0*f(i-1)+5265.d0*f(i-2)-5080.d0*f(i-3)+2970.d0*f(i-4)-972.d0*f(i-5)+137.d0*f(i-6)) /180.d0/dX ! backward difference, 6th-order accuracy
enddo
#elif (ghost_width == 5)
do i=1,lx-9
fxx(i) = ( (29531.d0/5040.d0)*f(i) - (962.d0/35.d0)*f(i+1) + (621.d0/10.d0)*f(i+2) &
- (4006.d0/45.d0)*f(i+3) + (691.d0/8.d0)*f(i+4) - (282.d0/5.d0)*f(i+5) &
+ (2143.d0/90.d0)*f(i+6) - (206.d0/35.d0)*f(i+7) + (363.d0/560.d0)*f(i+8) ) /dX
! forward difference, 8th-order accuracy
enddo
do i=lx-8,lx-4 ! width = 5
fxx(i) = ( - 9.d0*f(i-4) + 128.d0*f(i-3) - 1008.d0*f(i-2) + 8064.d0*f(i-1) &
- 14350.d0*f(i) + 8064.d0*f(i+1) - 1008.d0*f(i+2) + 128.d0*f(i+3) - 9.d0*f(i+4) ) / 5040.d0/dX
! centered difference, 8th-order accuracy
enddo
do i=lx-3,lx
fxx(i) = ( (29531.d0/5040.d0)*f(i) - (962.d0/35.d0)*f(i-1) + (621.d0/10.d0)*f(i-2) &
- (4006.d0/45.d0)*f(i-3) + (691.d0/8.d0)*f(i-4) - (282.d0/5.d0)*f(i-5) &
+ (2143.d0/90.d0)*f(i-6) - (206.d0/35.d0)*f(i-7) + (363.d0/560.d0)*f(i-8) ) /dX
! backward difference, 8th-order accuracy
enddo
#endif
return
end subroutine cdderivs_sw_x
!--------------------
subroutine rget_half_x(lx,f,hf)
implicit none
!~~~~~~% Input parameters:
integer,intent(in) :: lx
real*8,intent(in),dimension(lx) :: f
real*8,intent(out),dimension(lx) :: hf
#ifdef OLD
hf(1:lx-1) = (f(2:lx)+f(1:lx-1))/2.d0
hf(lx) = hf(lx-1)
#else
hf(lx) = 0.d0
#if (ghost_width == 2)
hf(1:lx-1) = (f(1:lx-1)+f(2:lx))/2.d0
#elif (ghost_width == 3)
hf(2:lx-2) = (-f(1:lx-3)+9.d0*f(2:lx-2)+9.d0*f(3:lx-1)-f(4:lx))/1.6d1
hf(1) = (f(2)+f(1))/2.d0
hf(lx-1) = (f(lx)+f(lx-1))/2.d0
#elif (ghost_width == 4)
hf(3:lx-3) = 3.d0/2.56d2*(f(1:lx-5)+f(6:lx))-2.5d1/2.56d2*(f(2:lx-4)+f(5:lx-1))+7.5d1/1.28d2*(f(3:lx-3)+f(4:lx-2))
hf(2) = (-f(1)+9.d0*f(2)+9.d0*f(3)-f(4))/1.6d1
hf(lx-2) = (-f(lx-1)+9.d0*f(lx-2)+9.d0*f(lx-3)-f(lx-4))/1.6d1
hf(1) = (f(2)+f(1))/2.d0
hf(lx-1) = (f(lx)+f(lx-1))/2.d0
#elif (ghost_width == 5)
hf(4:lx-4) = (-5.d0*(f(1:lx-7)+f(8:lx))+4.9d1*(f(2:lx-6)+f(7:lx-1))-2.45d2*(f(3:lx-5)+f(6:lx-2))+1.225d3*(f(4:lx-4)+f(5:lx-3))) &
/2.048d3
hf(3) = 3.d0/2.56d2*(f(1)+f(6))-2.5d1/2.56d2*(f(2)+f(5))+7.5d1/1.28d2*(f(3)+f(4))
hf(lx-3) = 3.d0/2.56d2*(f(lx)+f(lx-5))-2.5d1/2.56d2*(f(lx-1)+f(lx-4))+7.5d1/1.28d2*(f(lx-2)+f(lx-3))
hf(2) = (-f(1)+9.d0*f(2)+9.d0*f(3)-f(4))/1.6d1
hf(lx-2) = (-f(lx-1)+9.d0*f(lx-2)+9.d0*f(lx-3)-f(lx-4))/1.6d1
hf(1) = (f(2)+f(1))/2.d0
hf(lx-1) = (f(lx)+f(lx-1))/2.d0
#endif
#endif
return
end subroutine rget_half_x
!--------------------
subroutine rget_half_x_point(lx,f,hf,k)
implicit none
!~~~~~~% Input parameters:
integer,intent(in) :: lx,k
real*8,intent(in),dimension(lx) :: f
real*8,intent(out) :: hf
#ifdef OLD
if(k.eq.lx .or. k.eq.lx-1)then
hf = (f(lx)+f(lx-1))/2.d0
else
hf = (f(k+1)+f(k))/2.d0
endif
#else
if(k .eq. lx) hf = 0.d0
#if (ghost_width == 2)
if(k .lt. lx) hf = (f(k+1)+f(k))/2.d0
#elif (ghost_width == 3)
if(k .eq. 1) hf = (f(2)+f(1))/2.d0
if(k .eq. lx-1) hf = (f(lx)+f(lx-1))/2.d0
if(k .gt. 1 .and. k .lt. lx-1) hf = (-f(k-1)+9.d0*f(k)+9.d0*f(k+1)-f(k+2))/1.6d1
#elif (ghost_width == 4)
if(k .eq. 1) hf = (f(2)+f(1))/2.d0
if(k .eq. lx-1) hf = (f(lx)+f(lx-1))/2.d0
if(k .eq. 2) hf = (-f(1)+9.d0*f(2)+9.d0*f(3)-f(4))/1.6d1
if(k .eq. lx-2) hf =(-f(lx-1)+9.d0*f(lx-2)+9.d0*f(lx-3)-f(lx-4))/1.6d1
if(k .gt. 2 .and. k .lt. lx-2) hf = 3.d0/2.56d2*(f(k-2)+f(k+3))-2.5d1/2.56d2*(f(k-1)+f(k+2))+7.5d1/1.28d2*(f(k)+f(k+1))
#elif (ghost_width == 5)
if(k .eq. 1) hf = (f(2)+f(1))/2.d0
if(k .eq. lx-1) hf = (f(lx)+f(lx-1))/2.d0
if(k .eq. 2) hf = (-f(1)+9.d0*f(2)+9.d0*f(3)-f(4))/1.6d1
if(k .eq. lx-2) hf =(-f(lx-1)+9.d0*f(lx-2)+9.d0*f(lx-3)-f(lx-4))/1.6d1
if(k .eq. 3) hf = 3.d0/2.56d2*(f(1)+f(6))-2.5d1/2.56d2*(f(2)+f(5))+7.5d1/1.28d2*(f(3)+f(4))
if(k .eq. lx-3) hf = 3.d0/2.56d2*(f(lx)+f(lx-5))-2.5d1/2.56d2*(f(lx-1)+f(lx-4))+7.5d1/1.28d2*(f(lx-2)+f(lx-3))
if(k .gt. 3 .and. k .lt. lx-3) then
hf = (-5.d0*(f(k-3)+f(k+4))+4.9d1*(f(k-2)+f(k+3))-2.45d2*(f(k-1)+f(k+2))+1.225d3*(f(k)+f(k+1)))/2.048d3
endif
#endif
#endif
return
end subroutine rget_half_x_point
!--------------------
subroutine cget_half_x(lx,f,hf)
implicit none
!~~~~~~% Input parameters:
integer,intent(in) :: lx
double complex,intent(in),dimension(lx) :: f
double complex,intent(out),dimension(lx) :: hf
#ifdef OLD
hf(1:lx-1) = (f(2:lx)+f(1:lx-1))/2.d0
hf(lx) = hf(lx-1)
#else
hf(lx) = 0.d0
#if (ghost_width == 2)
hf(1:lx-1) = (f(1:lx-1)+f(2:lx))/2.d0
#elif (ghost_width == 3)
hf(2:lx-2) = (-f(1:lx-3)+9.d0*f(2:lx-2)+9.d0*f(3:lx-1)-f(4:lx))/1.6d1
hf(1) = (f(2)+f(1))/2.d0
hf(lx-1) = (f(lx)+f(lx-1))/2.d0
#elif (ghost_width == 4)
hf(3:lx-3) = 3.d0/2.56d2*(f(1:lx-5)+f(6:lx))-2.5d1/2.56d2*(f(2:lx-4)+f(5:lx-1))+7.5d1/1.28d2*(f(3:lx-3)+f(4:lx-2))
hf(2) = (-f(1)+9.d0*f(2)+9.d0*f(3)-f(4))/1.6d1
hf(lx-2) = (-f(lx-1)+9.d0*f(lx-2)+9.d0*f(lx-3)-f(lx-4))/1.6d1
hf(1) = (f(2)+f(1))/2.d0
hf(lx-1) = (f(lx)+f(lx-1))/2.d0
#elif (ghost_width == 5)
hf(4:lx-4) = (-5.d0*(f(1:lx-7)+f(8:lx))+4.9d1*(f(2:lx-6)+f(7:lx-1))-2.45d2*(f(3:lx-5)+f(6:lx-2))+1.225d3*(f(4:lx-4)+f(5:lx-3))) &
/2.048d3
hf(3) = 3.d0/2.56d2*(f(1)+f(6))-2.5d1/2.56d2*(f(2)+f(5))+7.5d1/1.28d2*(f(3)+f(4))
hf(lx-3) = 3.d0/2.56d2*(f(lx)+f(lx-5))-2.5d1/2.56d2*(f(lx-1)+f(lx-4))+7.5d1/1.28d2*(f(lx-2)+f(lx-3))
hf(2) = (-f(1)+9.d0*f(2)+9.d0*f(3)-f(4))/1.6d1
hf(lx-2) = (-f(lx-1)+9.d0*f(lx-2)+9.d0*f(lx-3)-f(lx-4))/1.6d1
hf(1) = (f(2)+f(1))/2.d0
hf(lx-1) = (f(lx)+f(lx-1))/2.d0
#endif
#endif
return
end subroutine cget_half_x
!---------
subroutine derivs_eth(ex,X,Y,f,eth_f,spin,e, &
quR1,quR2,quI1,quI2,gR,gI)
implicit none
!~~~~~~% Input parameters:
! spin: spin weight of f; e: eth (1) or eth bar (-1)
integer,intent(in) :: spin,e
integer,intent(in),dimension(2) :: ex
real*8,intent(in),dimension(ex(1)) :: X
real*8,intent(in),dimension(ex(2)) :: Y
double complex,intent(in),dimension(ex(1),ex(2)) :: f
double complex,intent(out),dimension(ex(1),ex(2)) :: eth_f
real*8,intent(in),dimension(ex(1),ex(2)) :: quR1,quR2,quI1,quI2,gR,gI
double complex,dimension(ex(1),ex(2)) :: fx,fy
double complex,dimension(ex(1),ex(2)) :: qu1,qu2,gama
integer :: i,j
!sanity check
if(e.ne.1 .and.e.ne.-1)then
write(*,*) "derivs_eth: bad input e = ",e
return
endif
qu1 = dcmplx(quR1,quI1)
qu2 = dcmplx(quR2,quI2)
gama = dcmplx(gR,gI)
do j=1,ex(2)
call cderivs_x(ex(1),X,f(:,j),fx(:,j))
enddo
do i=1,ex(1)
call cderivs_x(ex(2),Y,f(i,:),fy(i,:))
enddo
if(e.eq.1)then
eth_f = qu1*fx+qu2*fy+spin*gama*f
else
eth_f = dconjg(qu1)*fx+dconjg(qu2)*fy-spin*dconjg(gama)*f
endif
return
end subroutine derivs_eth
!---------
! dqu: q^B @_B q^A
! bdqu: bar{q}^B @_B q^A
! dg: q^B @_B Gamma
! bdg: q^B @_B bar{Gamma}
subroutine dderivs_eth(ex,X,Y,f,eth_f,spin,e1,e2, &
quR1,quR2,quI1,quI2,gR,gI, &
dquR1,dquR2,dquI1,dquI2, &
bdquR1,bdquR2,bdquI1,bdquI2, &
dgR,dgI,bdgR,bdgI)
implicit none
!~~~~~~% Input parameters:
! spin: spin weight of f; e: eth (1) or eth bar (-1)
integer,intent(in) :: spin,e1,e2
integer,intent(in),dimension(2) :: ex
real*8,intent(in),dimension(ex(1)) :: X
real*8,intent(in),dimension(ex(2)) :: Y
double complex,intent(in),dimension(ex(1),ex(2)) :: f
double complex,intent(out),dimension(ex(1),ex(2)) :: eth_f
real*8,intent(in),dimension(ex(1),ex(2)) :: quR1,quR2,quI1,quI2,gR,gI
real*8,intent(in),dimension(ex(1),ex(2)) :: dquR1,dquR2,dquI1,dquI2
real*8,intent(in),dimension(ex(1),ex(2)) :: bdquR1,bdquR2,bdquI1,bdquI2
real*8,intent(in),dimension(ex(1),ex(2)) :: dgR,dgI,bdgR,bdgI
double complex,dimension(ex(1),ex(2)) :: fx,fy,fxx,fxy,fyy
double complex,dimension(ex(1),ex(2)) :: qu1,qu2,dqu1,dqu2,bdqu1,bdqu2,gama,dgama,bdgama
integer :: i,j
!sanity check
if((e1.ne.1 .and. e1.ne.-1).or.(e2.ne.1 .and. e2.ne.-1))then
write(*,*) "dderivs_eth: bad input e1 = ",e1,"e2 = ",e2
return
endif
qu1 = dcmplx(quR1,quI1)
qu2 = dcmplx(quR2,quI2)
gama = dcmplx(gR,gI)
dqu1 = dcmplx(dquR1,dquI1)
dqu2 = dcmplx(dquR2,dquI2)
bdqu1 = dcmplx(bdquR1,bdquI1)
bdqu2 = dcmplx(bdquR2,bdquI2)
dgama = dcmplx(dgR,dgI)
bdgama = dcmplx(bdgR,bdgI)
do j=1,ex(2)
call cderivs_x(ex(1),X,f(:,j),fx(:,j))
call cdderivs_x(ex(1),X,f(:,j),fxx(:,j))
enddo
do i=1,ex(1)
call cderivs_x(ex(2),Y,f(i,:),fy(i,:))
call cdderivs_x(ex(2),Y,f(i,:),fyy(i,:))
enddo
do j=1,ex(2)
call cderivs_x(ex(1),X,fy(:,j),fxy(:,j))
enddo
if(e1.eq.1.and.e2.eq.1)then
! eth_eth
eth_f = qu1*qu1*fxx+2.d0*qu1*qu2*fxy+qu2*qu2*fyy &
+(dqu1+(2*spin+1)*gama*qu1)*fx+(dqu2+(2*spin+1)*gama*qu2)*fy &
+spin*(dgama+(spin+1)*gama*gama)*f
elseif(e1.eq.1.and.e2.eq.-1)then
! eth_ethb
eth_f = qu1*dconjg(qu1)*fxx+qu1*dconjg(qu2)*fxy+qu2*dconjg(qu1)*fxy+qu2*dconjg(qu2)*fyy &
+(dconjg(bdqu1)-spin*dconjg(gama)*qu1+(spin-1)*gama*dconjg(qu1))*fx &
+(dconjg(bdqu2)-spin*dconjg(gama)*qu2+(spin-1)*gama*dconjg(qu2))*fy &
-spin*(bdgama+(spin-1)*gama*dconjg(gama))*f
elseif(e1.eq.-1.and.e2.eq.1)then
! ethb_eth
eth_f = qu1*dconjg(qu1)*fxx+qu1*dconjg(qu2)*fxy+qu2*dconjg(qu1)*fxy+qu2*dconjg(qu2)*fyy &
+(bdqu1+spin*gama*dconjg(qu1)-(spin+1)*dconjg(gama)*qu1)*fx &
+(bdqu2+spin*gama*dconjg(qu2)-(spin+1)*dconjg(gama)*qu2)*fy &
+spin*(dconjg(bdgama)-(spin+1)*gama*dconjg(gama))*f
else
! ethb_ethb
eth_f = dconjg(qu1*qu1)*fxx+2.d0*dconjg(qu1*qu2)*fxy+dconjg(qu2*qu2)*fyy &
+(dconjg(dqu1)-(2*spin-1)*dconjg(gama*qu1))*fx+(dconjg(dqu2)-(2*spin-1)*dconjg(gama*qu2))*fy &
-spin*(dconjg(dgama)-(spin-1)*dconjg(gama*gama))*f
endif
return
end subroutine dderivs_eth
!------------------------------------------------------------------------------------------
! CQG 24S327, Eq.(19)--(22)
!------------------------------------------------------------------------------------------
subroutine setup_dyad(ex,crho,sigma,R, &
quR1,quR2,quI1,quI2, &
qlR1,qlR2,qlI1,qlI2, &
gR,gI, &
dquR1,dquR2,dquI1,dquI2, &
bdquR1,bdquR2,bdquI1,bdquI2, &
dgR,dgI,bdgR,bdgI, &
gx,gy,gz,sst,Rmin)
implicit none
!~~~~~~% Input parameters:
integer,intent(in ):: ex(1:3),sst
real*8,intent(in),dimension(ex(1))::crho
real*8,intent(in),dimension(ex(2))::sigma
real*8,intent(in),dimension(ex(3))::R
real*8,intent(in) :: Rmin
real*8,intent(out),dimension(ex(1),ex(2),ex(3)) :: quR1,quR2,quI1,quI2
real*8,intent(out),dimension(ex(1),ex(2),ex(3)) :: qlR1,qlR2,qlI1,qlI2
real*8,intent(out),dimension(ex(1),ex(2),ex(3)) :: gR,gI
real*8,intent(out),dimension(ex(1),ex(2),ex(3)) :: dquR1,dquR2,dquI1,dquI2
real*8,intent(out),dimension(ex(1),ex(2),ex(3)) :: bdquR1,bdquR2,bdquI1,bdquI2
real*8,intent(out),dimension(ex(1),ex(2),ex(3)) :: dgR,dgI,bdgR,bdgI
real*8,intent(out),dimension(ex(1),ex(2),ex(3)) :: gx,gy,gz
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,dimension(ex(1),ex(2),ex(3)) :: tp
do j=1,ex(2)
do i=1,ex(1)
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
qlR1(i,j,:) = thetac*cs/4.d0/tcts2
qlI1(i,j,:) = thetas*cs/4.d0/tcts2
qlR2(i,j,:) = thetac*cr/4.d0/tcts2
qlI2(i,j,:) =-thetas*cr/4.d0/tcts2
quR1(i,j,:) = 2.d0*tcts*thetas/cs
quI1(i,j,:) = 2.d0*tcts*thetac/cs
quR2(i,j,:) = 2.d0*tcts*thetas/cr
quI2(i,j,:) =-2.d0*tcts*thetac/cr
gR(i,j,:) = (crcs*crcs*(sr+ss)+(cr*cr-cs*cs)*(ss-sr))/4.d0/thetac/crcs
gI(i,j,:) = (crcs*crcs*(sr-ss)+(cs*cs-cr*cr)*(ss+sr))/4.d0/thetas/crcs
dquR1(i,j,:) = 3.d0*cr2*ss*(cs2+cr2*ss2)/2.d0/cr/cs2
dquI1(i,j,:) =-thetac*thetas*sr*(cs2+3.d0*cr2*ss2)/cr/cs2
dquR2(i,j,:) = 3.d0*cs2*sr*(cr2+cs2*sr2)/2.d0/cs/cr2
dquI2(i,j,:) = thetac*thetas*ss*(cr2+3.d0*cs2*sr2)/cs/cr2
bdquR1(i,j,:) = sr*cs2*(cs2+cr2*ss2)/2.d0/cr/cs2
bdquI1(i,j,:) = thetac*thetas*ss*(8.d0*tc2*ts2-3.d0*cs2*sr2-cr2)/cr/cs2
bdquR2(i,j,:) = ss*cr2*(cr2+cs2*sr2)/2.d0/cs/cr2
bdquI2(i,j,:) =-thetac*thetas*sr*(8.d0*tc2*ts2-3.d0*cr2*ss2-cs2)/cs/cr2
ggr = tc2*ts2/cr2/cs2
dgR(i,j,:) =-3.d0*ggr*sr*ss*(cs2+cr2)
dgI(i,j,:) = 6.d0*ggr*thetac*thetas*(cs2-cr2)
bdgR(i,j,:) = ggr*(2.d0*cr2*cs2+cr2+cs2)
bdgI(i,j,:) = 0.d0
do k=1,ex(3)
ggr = R(k)*Rmin/(1.d0-R(k))
tgrho = dtan(crho(i))
tgsigma = dtan(sigma(j))
select case (sst)
case (0)
gz(i,j,k) = ggr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
gx(i,j,k) = gz(i,j,k)*tgrho
gy(i,j,k) = gz(i,j,k)*tgsigma
case (1)
gz(i,j,k) = -ggr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
gx(i,j,k) = gz(i,j,k)*tgrho
gy(i,j,k) = gz(i,j,k)*tgsigma
case (2)
gx(i,j,k) = ggr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
gy(i,j,k) = gx(i,j,k)*tgrho
gz(i,j,k) = gx(i,j,k)*tgsigma
case (3)
gx(i,j,k) = -ggr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
gy(i,j,k) = gx(i,j,k)*tgrho
gz(i,j,k) = gx(i,j,k)*tgsigma
case (4)
gy(i,j,k) = ggr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
gx(i,j,k) = gy(i,j,k)*tgrho
gz(i,j,k) = gy(i,j,k)*tgsigma
case (5)
gy(i,j,k) = -ggr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
gx(i,j,k) = gy(i,j,k)*tgrho
gz(i,j,k) = gy(i,j,k)*tgsigma
case default
write(*,*) "setup_dyad: not recognized sst = ",sst
return
end select
enddo
enddo
enddo
! note the difference between right hand coordinate and left hand coordinate
if(sst==1 .or. sst==3 .or. sst==4)then
quI1 = -quI1
quI2 = -quI2
qlI1 = -qlI1
qlI2 = -qlI2
gI = -gI
dquI1 = -dquI1
dquI2 = -dquI2
bdquI1 = -bdquI1
bdquI2 = -bdquI2
dgI = -dgI
bdgI = -bdgI
endif
return
end subroutine setup_dyad
!---------------------------------------------------------------------------------
! interface to c
!---------------------------------------------------------------------------------
subroutine eth_derivs(ex,X,Y,Rfi,Ifi,Reth_f,Ieth_f,spin,e, &
quR1,quR2,quI1,quI2,gR,gI)
implicit none
!~~~~~~% Input parameters:
! spin: spin weight of f; e: eth (1) or eth bar (-1)
integer,intent(in) :: spin,e
integer,intent(in),dimension(3) :: ex
real*8,intent(in),dimension(ex(1)) :: X
real*8,intent(in),dimension(ex(2)) :: Y
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: Rfi,Ifi
real*8,intent(out),dimension(ex(1),ex(2),ex(3)) :: Reth_f,Ieth_f
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: quR1,quR2,quI1,quI2,gR,gI
double complex,dimension(ex(1),ex(2)) :: f,eth_f,fx,fy
double complex,dimension(ex(1),ex(2)) :: qu1,qu2,gama
integer :: k
do k=1,ex(3)
f = dcmplx(Rfi(:,:,k),Ifi(:,:,k))
call derivs_eth(ex(1:2),X,Y,f,eth_f,spin,e, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
Reth_f(:,:,k) = dreal(eth_f)
Ieth_f(:,:,k) = dimag(eth_f)
enddo
return
end subroutine eth_derivs
!---------------------------------------------------------------------------------
subroutine eth_dderivs(ex,X,Y,Rfi,Ifi,Reth_f,Ieth_f,spin,e1,e2, &
quR1,quR2,quI1,quI2,gR,gI, &
dquR1,dquR2,dquI1,dquI2, &
bdquR1,bdquR2,bdquI1,bdquI2, &
dgR,dgI,bdgR,bdgI)
implicit none
!~~~~~~% Input parameters:
! spin: spin weight of f; e: eth (1) or eth bar (-1)
integer,intent(in) :: spin,e1,e2
integer,intent(in),dimension(3) :: ex
real*8,intent(in),dimension(ex(1)) :: X
real*8,intent(in),dimension(ex(2)) :: Y
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: Rfi,Ifi
real*8,intent(out),dimension(ex(1),ex(2),ex(3)) :: Reth_f,Ieth_f
real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: quR1,quR2,quI1,quI2,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)) :: f,eth_f,fx,fy
double complex,dimension(ex(1),ex(2)) :: qu1,qu2,gama
integer :: k
do k=1,ex(3)
f = dcmplx(Rfi(:,:,k),Ifi(:,:,k))
call dderivs_eth(ex(1:2),X,Y,f,eth_f,spin,e1,e2, &
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))
Reth_f(:,:,k) = dreal(eth_f)
Ieth_f(:,:,k) = dimag(eth_f)
enddo
return
end subroutine eth_dderivs
!---------------------------------------------------------------------------------
! fill symmetric boundary buffer points
!---------------------------------------------------------------------------------
subroutine fill_symmetric_boundarybuffer(ex,crho,sigma,R,drho,dsigma, &
quR1,quR2,quI1,quI2,qlR1,qlR2,qlI1,qlI2, &
Rv,Iv,Symmetry,sst,spin)
implicit none
!~~~~~~% Input parameters:
integer,intent(in),dimension(3) :: ex
integer,intent(in) :: Symmetry,sst,spin
real*8,intent(in),dimension(ex(1))::crho
real*8,intent(in),dimension(ex(2))::sigma
real*8,intent(in),dimension(ex(3))::R
real*8,intent(in) :: drho,dsigma
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(inout),dimension(ex(1),ex(2),ex(3)) :: Rv,Iv
integer :: i,j,k,t
select case (Symmetry)
case (0)
return
case (1)
if((sst==2.or.sst==4).and.dabs(sigma(1)+ghost_width*dsigma) < dsigma/2.d0)then
do k=1,ex(3)
do j=1,ghost_width
do i=1,ex(1)
#ifdef Vertex
#ifdef Cell
#error Both Cell and Vertex are defined
#endif
t = 2*ghost_width+2-j
#endif
#ifdef Cell
#ifdef Vertex
#error Both Cell and Vertex are defined
#endif
t = 2*ghost_width+1-j
#endif
call symmetrymap(quR1(i,j,k),quR2(i,j,k),quI1(i,j,k),quI2(i,j,k),qlR1(i,j,k),qlR2(i,j,k),qlI1(i,j,k),qlI2(i,j,k), &
quR1(i,t,k),quR2(i,t,k),quI1(i,t,k),quI2(i,t,k),qlR1(i,t,k),qlR2(i,t,k),qlI1(i,t,k),qlI2(i,t,k), &
Rv(i,j,k),Iv(i,j,k),Rv(i,t,k),Iv(i,t,k),2,spin)
enddo
enddo
enddo
endif
if((sst==3.or.sst==5).and.dabs(sigma(ex(2))-ghost_width*dsigma) < dsigma/2.d0)then
do k=1,ex(3)
do j=ex(2)-ghost_width+1,ex(2)
do i=1,ex(1)
t = ex(2)-j+1
#ifdef Vertex
#ifdef Cell
#error Both Cell and Vertex are defined
#endif
t = ex(2)-2*ghost_width-1+t
#endif
#ifdef Cell
#ifdef Vertex
#error Both Cell and Vertex are defined
#endif
t = ex(2)-2*ghost_width+t
#endif
call symmetrymap(quR1(i,j,k),quR2(i,j,k),quI1(i,j,k),quI2(i,j,k),qlR1(i,j,k),qlR2(i,j,k),qlI1(i,j,k),qlI2(i,j,k), &
quR1(i,t,k),quR2(i,t,k),quI1(i,t,k),quI2(i,t,k),qlR1(i,t,k),qlR2(i,t,k),qlI1(i,t,k),qlI2(i,t,k), &
Rv(i,j,k),Iv(i,j,k),Rv(i,t,k),Iv(i,t,k),2,spin)
enddo
enddo
enddo
endif
case (2)
if(dabs(crho(1)+ghost_width*drho) < drho/2.d0)then
if(dabs(sigma(1)+ghost_width*dsigma) < dsigma/2.d0)then
do k=1,ex(3)
do j=1,ghost_width
do i=ghost_width+1,ex(1)
#ifdef Vertex
#ifdef Cell
#error Both Cell and Vertex are defined
#endif
t = 2*ghost_width+2-j
#endif
#ifdef Cell
#ifdef Vertex
#error Both Cell and Vertex are defined
#endif
t = 2*ghost_width+1-j
#endif
call symmetrymap(quR1(i,j,k),quR2(i,j,k),quI1(i,j,k),quI2(i,j,k),qlR1(i,j,k),qlR2(i,j,k),qlI1(i,j,k),qlI2(i,j,k), &
quR1(i,t,k),quR2(i,t,k),quI1(i,t,k),quI2(i,t,k),qlR1(i,t,k),qlR2(i,t,k),qlI1(i,t,k),qlI2(i,t,k), &
Rv(i,j,k),Iv(i,j,k),Rv(i,t,k),Iv(i,t,k),2,spin)
enddo
enddo
enddo
endif
do k=1,ex(3)
do j=1,ex(2)
do i=1,ghost_width
#ifdef Vertex
#ifdef Cell
#error Both Cell and Vertex are defined
#endif
t = 2*ghost_width+2-i
#endif
#ifdef Cell
#ifdef Vertex
#error Both Cell and Vertex are defined
#endif
t = 2*ghost_width+1-i
#endif
call symmetrymap(quR1(i,j,k),quR2(i,j,k),quI1(i,j,k),quI2(i,j,k),qlR1(i,j,k),qlR2(i,j,k),qlI1(i,j,k),qlI2(i,j,k), &
quR1(t,j,k),quR2(t,j,k),quI1(t,j,k),quI2(t,j,k),qlR1(t,j,k),qlR2(t,j,k),qlI1(t,j,k),qlI2(t,j,k), &
Rv(i,j,k),Iv(i,j,k),Rv(t,j,k),Iv(t,j,k),1,spin)
enddo
enddo
enddo
else
if(dabs(sigma(1)+ghost_width*dsigma) < dsigma/2.d0)then
do k=1,ex(3)
do j=1,ghost_width
do i=1,ex(1)
#ifdef Vertex
#ifdef Cell
#error Both Cell and Vertex are defined
#endif
t = 2*ghost_width+2-j
#endif
#ifdef Cell
#ifdef Vertex
#error Both Cell and Vertex are defined
#endif
t = 2*ghost_width+1-j
#endif
call symmetrymap(quR1(i,j,k),quR2(i,j,k),quI1(i,j,k),quI2(i,j,k),qlR1(i,j,k),qlR2(i,j,k),qlI1(i,j,k),qlI2(i,j,k), &
quR1(i,t,k),quR2(i,t,k),quI1(i,t,k),quI2(i,t,k),qlR1(i,t,k),qlR2(i,t,k),qlI1(i,t,k),qlI2(i,t,k), &
Rv(i,j,k),Iv(i,j,k),Rv(i,t,k),Iv(i,t,k),2,spin)
enddo
enddo
enddo
endif
endif
end select
return
end subroutine fill_symmetric_boundarybuffer
!-------------------------------------------------------------------------------------------
! note Theta has spin weight 2, but the its determinant does not satisfy our
! assumption, so in principle this routine can not be applied to Theta
! but since we never calculate eth derivative on Theta, it does not matter
!-------------------------------------------------------------------------------------------
subroutine symmetrymap(oquR1,oquR2,oquI1,oquI2,oqlR1,oqlR2,oqlI1,oqlI2, &
iquR1,iquR2,iquI1,iquI2,iqlR1,iqlR2,iqlI1,iqlI2, &
oRv,oIv,iRv,iIv,ind,spin)
implicit none
!~~~~~~% Input parameters:
integer,intent(in) :: ind,spin
real*8,intent(in) :: iquR1,iquR2,iquI1,iquI2,iqlR1,iqlR2,iqlI1,iqlI2
real*8,intent(in) :: iRv,iIv
real*8,intent(in) :: oquR1,oquR2,oquI1,oquI2,oqlR1,oqlR2,oqlI1,oqlI2
real*8,intent(out) :: oRv,oIv
double complex :: iqu1,iqu2,iql1,iql2,oqu1,oqu2,oql1,oql2
real*8 :: h11,h12,h22,KK
double complex :: r
iqu1 = dcmplx(iquR1,iquI1)
iqu2 = dcmplx(iquR2,iquI2)
iql1 = dcmplx(iqlR1,iqlI1)
iql2 = dcmplx(iqlR2,iqlI2)
oqu1 = dcmplx(oquR1,oquI1)
oqu2 = dcmplx(oquR2,oquI2)
oql1 = dcmplx(oqlR1,oqlI1)
oql2 = dcmplx(oqlR2,oqlI2)
! spin weight
select case (spin)
case (-2)
r = dcmplx(iRv,iIv)
#if 1
KK=dsqrt(iRv*iRv+iIv*iIv+1.d0)
h11 = dreal(r*iql1*iql1+KK*dconjg(iql1)*iql1)
h12 =-dreal(r*iql1*iql2+KK*dconjg(iql1)*iql2)
h22 = dreal(r*iql2*iql2+KK*dconjg(iql2)*iql2)
#else
h11 = dreal(r*iql1*iql1)
h12 =-dreal(r*iql1*iql2)
h22 = dreal(r*iql2*iql2)
#endif
r = h11*dconjg(oqu1)*dconjg(oqu1)+2.d0*h12*dconjg(oqu1)*dconjg(oqu2)+h22*dconjg(oqu2)*dconjg(oqu2)
oRv = dreal(r)
oIv = dimag(r)
case (-1)
r = dcmplx(iRv,iIv)
h11 = dreal(r*iql1)
h22 = dreal(r*iql2)
if(ind == 1) h11 = -h11
if(ind == 2) h22 = -h22
r = h11*dconjg(oqu1)+h22*dconjg(oqu2)
oRv = dreal(r)
oIv = dimag(r)
case (0)
oRv = iRv
oIv = iIv
case(1)
r = dcmplx(iRv,iIv)
h11 = dreal(r*dconjg(iql1))
h22 = dreal(r*dconjg(iql2))
if(ind == 1) h11 = -h11
if(ind == 2) h22 = -h22
r = h11*oqu1+h22*oqu2
oRv = dreal(r)
oIv = dimag(r)
case (2)
r = dcmplx(iRv,iIv)
#if 1
KK=dsqrt(iRv*iRv+iIv*iIv+1.d0)
h11 = dreal(r*dconjg(iql1)*dconjg(iql1)+KK*dconjg(iql1)*iql1)
h12 =-dreal(r*dconjg(iql1)*dconjg(iql2)+KK*dconjg(iql1)*iql2)
h22 = dreal(r*dconjg(iql2)*dconjg(iql2)+KK*dconjg(iql2)*iql2)
#else
h11 = dreal(r*dconjg(iql1)*dconjg(iql1))
h12 =-dreal(r*dconjg(iql1)*dconjg(iql2))
h22 = dreal(r*dconjg(iql2)*dconjg(iql2))
#endif
r = h11*oqu1*oqu1+2.d0*h12*oqu1*oqu2+h22*oqu2*oqu2
oRv = dreal(r)
oIv = dimag(r)
case default
write(*,*)"symmetrymap does not yet support spin weight ",spin
end select
return
end subroutine symmetrymap
!-------------------------
subroutine calculate_K(ex,X,Y,Z,RJ,IJ,KK,HKK,KKx,HKKx)
implicit none
!~~~~~~% Input parameters:
integer,intent(in) :: ex(3)
real*8,intent(in) :: X(ex(1)),Y(ex(2)),Z(ex(3))
real*8,dimension(ex(1),ex(2),ex(3)),intent(in) :: RJ,IJ
real*8,dimension(ex(1),ex(2),ex(3)),intent(out) :: KK,HKK,KKx,HKKx
integer :: i,j
KK = dsqrt(1.d0+RJ*RJ+IJ*IJ)
do j=1,ex(2)
do i=1,ex(1)
call rget_half_x(ex(3),KK(i,j,:),HKK(i,j,:))
call rderivs_x(ex(3),Z,KK(i,j,:),KKx(i,j,:))
call rget_half_x(ex(3),KKx(i,j,:),HKKx(i,j,:))
enddo
enddo
return
end subroutine calculate_K
!--------------------------------------------------------------------
! this R is indeed x
function Eq_Theta(ex,crho,sigma,R,RJ,IJ,RU,IU,beta,RB,IB, &
Rnu,Inu,Rk,Ik,RTheta,ITheta,W,Rmin, &
qlR1,qlR2,qlI1,qlI2,quR1,quR2,quI1,quI2,gR,gI) result(gont)
implicit none
integer,intent(in ):: ex(1:3)
real*8,intent(in),dimension(ex(1))::crho
real*8,intent(in),dimension(ex(2))::sigma
real*8,intent(in),dimension(ex(3))::R
real*8, dimension(ex(1),ex(2),ex(3)),intent(in) :: beta,W
real*8, dimension(ex(1),ex(2),ex(3)),intent(in) :: RJ,IJ,RU,IU,RB,IB,Rnu,Inu,Rk,Ik
real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: RTheta,ITheta
real*8,intent(in) :: Rmin
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: qlR1,qlR2,qlI1,qlI2
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: quR1,quR2,quI1,quI2
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: gR,gI
! gont = 0: success; gont = 1: something wrong
integer::gont
double complex,dimension(ex(1),ex(2),ex(3)) :: CU,DCU,bDCU,CB,DCB,bDCB,CJ,DCJ
double complex :: CTheta0,CTheta,CTheta1,RHS
integer :: i,j,k,RK4
double complex,dimension(ex(3)) :: Cnu,Ck,CUx,DCUx,bDCUx
double complex,dimension(ex(3)) :: Cnux,CJx,CJxx,DCJx
double complex,dimension(ex(3)) :: fCTheta,CThetax
real*8,dimension(ex(3)) :: KK,KKx,betax,Wx
double complex :: Theta_rhs,Theta_rhs_o
real*8 :: dR
!!! sanity check
dR = sum(RJ)+sum(IJ)+sum(RU)+sum(IU)+sum(beta)+sum(RB)+sum(IB) + &
sum(Rnu)+sum(Inu)+sum(Rk)+sum(Ik)+sum(RTheta)+sum(ITheta)
if(dR.ne.dR) then
if(sum(RJ).ne.sum(RJ))write(*,*)"NullEvol_Theta: find NaN in RJ"
if(sum(IJ).ne.sum(IJ))write(*,*)"NullEvol_Theta: find NaN in IJ"
if(sum(RU).ne.sum(RU))write(*,*)"NullEvol_Theta: find NaN in RU"
if(sum(IU).ne.sum(IU))write(*,*)"NullEvol_Theta: find NaN in IU"
if(sum(beta).ne.sum(beta))write(*,*)"NullEvol_Theta: find NaN in beta"
if(sum(RB).ne.sum(RB))write(*,*)"NullEvol_Theta: find NaN in RB"
if(sum(IB).ne.sum(IB))write(*,*)"NullEvol_Theta: find NaN in IB"
if(sum(Rnu).ne.sum(Rnu))write(*,*)"NullEvol_Theta: find NaN in Rnu"
if(sum(Inu).ne.sum(Inu))write(*,*)"NullEvol_Theta: find NaN in Inu"
if(sum(Rk).ne.sum(Rk))write(*,*)"NullEvol_Theta: find NaN in Rk"
if(sum(Ik).ne.sum(Ik))write(*,*)"NullEvol_Theta: find NaN in Ik"
if(sum(RTheta).ne.sum(RTheta))write(*,*)"NullEvol_Theta: find NaN in RTheta"
if(sum(ITheta).ne.sum(ITheta))write(*,*)"NullEvol_Theta: find NaN in ITheta"
gont = 1
return
endif
dR = R(2) - R(1)
CU = dcmplx(RU,IU)
CB = dcmplx(RB,IB)
CJ = dcmplx(RJ,IJ)
do k=1,ex(3)
call derivs_eth(ex(1:2),crho,sigma,CU(:,:,k),DCU(:,:,k),1,1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
call derivs_eth(ex(1:2),crho,sigma,CU(:,:,k),bDCU(:,:,k),1,-1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
call derivs_eth(ex(1:2),crho,sigma,CB(:,:,k),DCB(:,:,k),1,1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
call derivs_eth(ex(1:2),crho,sigma,CB(:,:,k),bDCB(:,:,k),1,-1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
call derivs_eth(ex(1:2),crho,sigma,CJ(:,:,k),DCJ(:,:,k),2,1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
enddo
do j=1,ex(2)
do i=1,ex(1)
CTheta0 = dcmplx(RTheta(i,j,1),ITheta(i,j,1))
fCTheta = dcmplx(RTheta(i,j,:),ITheta(i,j,:))
call cderivs_x(ex(3),R,fCTheta,CThetax)
Cnu = dcmplx(Rnu(i,j,:),Inu(i,j,:))
Ck = dcmplx(Rk(i,j,:),Ik(i,j,:))
call cderivs_x(ex(3),R,Cnu,Cnux)
call rderivs_x(ex(3),R,beta(i,j,:),betax)
KK = dsqrt(1.d0+RJ(i,j,:)*RJ(i,j,:)+IJ(i,j,:)*IJ(i,j,:))
call rderivs_x(ex(3),R,KK,KKx)
call rderivs_x(ex(3),R,W,Wx)
call cderivs_x(ex(3),R,DCU(i,j,:),DCUx)
call cderivs_x(ex(3),R,CU(i,j,:),CUx)
call cderivs_x(ex(3),R,bDCU(i,j,:),bDCUx)
call cderivs_x(ex(3),R,CJ(i,j,:),CJx)
call cdderivs_x(ex(3),R,CJ(i,j,:),CJxx)
call cderivs_x(ex(3),R,DCJ(i,j,:),DCJx)
RTheta(i,j,1) = 0.d0
ITheta(i,j,1) = 0.d0
do k=1,ex(3)-1
#if 0
call check_daxiao(beta(i,j,k))
call check_daxiao(CU(i,j,k))
call check_daxiao(bDCU(i,j,k))
call check_daxiao(DCU(i,j,k))
call check_daxiao(CB(i,j,k))
call check_daxiao(DCB(i,j,k))
call check_daxiao(W(i,j,k))
call check_daxiao(CJ(i,j,k))
call check_daxiao(DCJ(i,j,k))
call check_daxiao(bDCB(i,j,k))
call check_daxiao(Cnu(k))
call check_daxiao(Ck(k))
call check_daxiao(fCTheta(k))
call check_daxiao(betax(k))
call check_daxiao(KKx(k))
call check_daxiao(CUx(k))
call check_daxiao(DCUx(k))
call check_daxiao(bDCUx(k))
call check_daxiao(Wx(k))
call check_daxiao(CJx(k))
call check_daxiao(CJxx(k))
call check_daxiao(DCJx(k))
call check_daxiao(Cnux(k))
#endif
RHS = Theta_rhs(R(k),Rmin,beta(i,j,k),betax(k),KK(k),KKx(k),CU(i,j,k),CUx(k),DCUx(k),bDCU(i,j,k),bDCUx(k), &
DCU(i,j,k),CB(i,j,k),DCB(i,j,k),W(i,j,k),Wx(k),CJ(i,j,k),DCJ(i,j,k), &
CJx(k),CJxx(k),DCJx(k),bDCB(i,j,k),Cnu(k),Cnux(k),Ck(k),fCTheta(k))
RHS = RHS - CThetax(k)
RTheta(i,j,k+1) = dreal(RHS)
ITheta(i,j,k+1) = dimag(RHS)
enddo
write(*,*)RTheta(i,j,:)
stop
enddo
enddo
gont = 0
return
end function Eq_Theta
subroutine check_daxiao(f)
implicit none
double complex,intent(in) :: f
real*8,parameter :: eps=1.d-5
if(cdabs(f)>eps) write(*,*) f
return
end subroutine check_daxiao
subroutine check_factor(T,crho,sigma,R,sst,Rmin)
implicit none
integer,intent(in) :: sst
real*8,intent(in) :: T,crho,sigma,R,Rmin
real*8 ::x,y,z,hgr,gt,gp,tgrho,tgsigma,tc,ts,rf
double complex :: Yslm,II,Jr,swtf,ff
hgr = R*Rmin/(1.d0-R)
tgrho = dtan(crho)
tgsigma = dtan(sigma)
tc = dsqrt((1.d0-dsin(crho)*dsin(sigma))/2.d0)
ts = dsqrt((1.d0+dsin(crho)*dsin(sigma))/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)
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
write(*,*) dsqrt(dble((2-1)*2*(2+1)*(2+2)))*Yslm(2,2,0,gt,gp)*swtf**2
return
end subroutine check_factor
subroutine getdxs(T,crho,sigma,R,betax,KKx,CUx,DCUx,bDCUx,Wx,CJx,CJxx,DCJx,Cnux,CThetax,sst,Rmin)
implicit none
integer,intent(in) :: sst
real*8,intent(in) :: T,crho,sigma,R,Rmin
real*8,intent(out) :: betax,KKx,Wx
double complex,intent(out) :: CUx,DCUx,bDCUx,CJx,CJxx,DCJx,Cnux,CThetax
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)
hgr = R*Rmin/(1.d0-R)
tgrho = dtan(crho)
tgsigma = dtan(sigma)
tc = dsqrt((1.d0-dsin(crho)*dsin(sigma))/2.d0)
ts = dsqrt((1.d0+dsin(crho)*dsin(sigma))/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)
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
betax = 0.d0
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(Yslm(0,2,m,gt,gp))*dreal(Jr*cdexp(II*nu*T))*(Rmin+hgr)**2/Rmin
! Wx = dreal(Jr*cdexp(II*nu*T))*(Rmin+hgr)**2/Rmin
KKx = 0.d0
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 = dsqrt(dble(2*(2+1)))*Yslm(1,2,m,gt,gp)*swtf*rf*(Rmin+hgr)**2/Rmin
! CUx = rf*(Rmin+hgr)**2/Rmin
DCUx = dsqrt(dble((2-1)*2*(2+1)*(2+2)))*Yslm(2,2,m,gt,gp)*swtf**2*rf*(Rmin+hgr)**2/Rmin
! DCUx = rf*(Rmin+hgr)**2/Rmin
bDCUx =-dble(2*(2+1))*Yslm(0,2,m,gt,gp)*rf*(Rmin+hgr)**2/Rmin
! bDCUx = rf*(Rmin+hgr)**2/Rmin
Jr = -C1/4.d0/hgr**2+C2/4.d0/hgr**4
rf = dreal(Jr*cdexp(II*nu*T))
CJx = dsqrt(dble((2-1)*2*(2+1)*(2+2)))*Yslm(2,2,m,gt,gp)*swtf**2*rf*(Rmin+hgr)**2/Rmin
! CJx = rf*(Rmin+hgr)**2/Rmin
Cnux =-dble((2-1)*(2+2))*dsqrt(dble(2*(2+1)))*Yslm(1,2,m,gt,gp)*swtf*rf*(Rmin+hgr)**2/Rmin
! Cnux = rf*(Rmin+hgr)**2/Rmin
DCJx = 0.d0
rf = dreal(Jr*II*nu*cdexp(II*nu*T))
CThetax = dsqrt(dble((2-1)*2*(2+1)*(2+2)))*Yslm(2,2,m,gt,gp)*swtf**2*rf*(Rmin+hgr)**2/Rmin
! CThetax = rf*(Rmin+hgr)**2/Rmin
Jr = C1/2.d0/hgr**3-C2/hgr**5
rf = dreal(Jr*cdexp(II*nu*T))
CJxx = dsqrt(dble((2-1)*2*(2+1)*(2+2)))*Yslm(2,2,m,gt,gp)*swtf**2*rf*(Rmin+hgr)**4/Rmin**2+2.d0*(Rmin+hgr)/Rmin*CJx
! CJxx = rf*(Rmin+hgr)**4/Rmin**2+2.d0*(Rmin+hgr)/Rmin*CJx
#if 0
DCUx = DCUx*dsqrt(dble((2-1)*2*(2+1)*(2+2)))*Yslm(2,2,m,gt,gp)*swtf**2
CJx = CJx*dsqrt(dble((2-1)*2*(2+1)*(2+2)))*Yslm(2,2,m,gt,gp)*swtf**2
CJxx = CJxx*dsqrt(dble((2-1)*2*(2+1)*(2+2)))*Yslm(2,2,m,gt,gp)*swtf**2
CThetax = CThetax*dsqrt(dble((2-1)*2*(2+1)*(2+2)))*Yslm(2,2,m,gt,gp)*swtf**2
#endif
return
end subroutine getdxs
subroutine getndxs(T,crho,sigma,R,beta,KK,CU,bDCU,DCU,CB,DCB,W,CJ,DCJ,bDCB,Cnu,Ck,CTheta,sst,Rmin)
implicit none
integer,intent(in) :: sst
real*8,intent(in) :: T,crho,sigma,R,Rmin
real*8,intent(out) :: beta,KK,W
double complex,intent(out) :: CU,bDCU,DCU,CB,DCB,CJ,DCJ,bDCB,Cnu,Ck,CTheta
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)
hgr = R*Rmin/(1.d0-R)
tgrho = dtan(crho)
tgsigma = dtan(sigma)
tc = dsqrt((1.d0-dsin(crho)*dsin(sigma))/2.d0)
ts = dsqrt((1.d0+dsin(crho)*dsin(sigma))/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)
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
beta = dreal(Yslm(0,2,m,gt,gp))*dreal(beta0*cdexp(II*nu*T))
! beta = dreal(beta0*cdexp(II*nu*T))
CB = dsqrt(dble(2*(2+1)))*Yslm(1,2,m,gt,gp)*swtf*dreal(beta0*cdexp(II*nu*T))
! CB = dreal(beta0*cdexp(II*nu*T))
DCB = dsqrt(dble((2-1)*2*(2+1)*(2+2)))*Yslm(2,2,m,gt,gp)*swtf**2*dreal(beta0*cdexp(II*nu*T))
! DCB = dreal(beta0*cdexp(II*nu*T))
bDCB =-dble(2*(2+1))*Yslm(0,2,m,gt,gp)*dreal(beta0*cdexp(II*nu*T))
! bDCB = 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 = dreal(Yslm(0,2,m,gt,gp))*dreal(Jr*cdexp(II*nu*T))
! W = dreal(Jr*cdexp(II*nu*T))
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))
CJ = dsqrt(dble((2-1)*2*(2+1)*(2+2)))*Yslm(2,2,m,gt,gp)*swtf**2*rf
! CJ = rf
DCJ = 0.d0
Cnu =-dsqrt(dble((2+2)*(2-2+1)*(2-1)*2*(2+1)*(2+2)))*Yslm(1,2,m,gt,gp)*swtf*rf
! Cnu = rf
KK = dsqrt(1.d0+cdabs(CJ)**2)
Ck = 0.d0
rf = dreal(Jr*II*nu*cdexp(II*nu*T))
CTheta = dsqrt(dble((2-1)*2*(2+1)*(2+2)))*Yslm(2,2,m,gt,gp)*swtf**2*rf
! CTheta = rf
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))
CU = dsqrt(dble(2*(2+1)))*Yslm(1,2,m,gt,gp)*swtf*rf
! CU = rf
DCU = dsqrt(dble((2-1)*2*(2+1)*(2+2)))*Yslm(2,2,m,gt,gp)*swtf**2*rf
! DCU = rf
bDCU =-dble(2*(2+1))*Yslm(0,2,m,gt,gp)*rf
! bDCU = rf
#if 0
DCU = DCU*dsqrt(dble((2-1)*2*(2+1)*(2+2)))*Yslm(2,2,m,gt,gp)*swtf**2
DCB = DCB*dsqrt(dble((2-1)*2*(2+1)*(2+2)))*Yslm(2,2,m,gt,gp)*swtf**2
CTheta = CTheta*dsqrt(dble((2-1)*2*(2+1)*(2+2)))*Yslm(2,2,m,gt,gp)*swtf**2
#endif
return
end subroutine getndxs
!--------------------------------------------------------------------
! this R is indeed x
function Eq_Theta_2(ex,crho,sigma,R,RJ,IJ,RU,IU,beta,RB,IB, &
Rnu,Inu,Rk,Ik,RTheta,ITheta,W,Rmin, &
qlR1,qlR2,qlI1,qlI2,quR1,quR2,quI1,quI2,gR,gI, &
T,sst) result(gont)
implicit none
integer,intent(in ):: ex(1:3),sst
real*8,intent(in),dimension(ex(1))::crho
real*8,intent(in),dimension(ex(2))::sigma
real*8,intent(in),dimension(ex(3))::R
real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: beta,W
real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: RJ,IJ,RU,IU,RB,IB,Rnu,Inu,Rk,Ik
real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: RTheta,ITheta
real*8,intent(in) :: Rmin,T
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: qlR1,qlR2,qlI1,qlI2
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: quR1,quR2,quI1,quI2
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: gR,gI
! gont = 0: success; gont = 1: something wrong
integer::gont
double complex,dimension(ex(1),ex(2),ex(3)) :: CU,DCU,bDCU,CB,DCB,bDCB,CJ,DCJ
double complex :: CTheta0,CTheta,CTheta1,RHS
integer :: i,j,k,RK4
double complex,dimension(ex(3)) :: Cnu,Ck,HCnu,HCk,HCU,HDCU,CUx,HCUx,DCUx,HDCUx,HbDCU,bDCUx,HbDCUx
double complex,dimension(ex(3)) :: Cnux,HCnux,HCJ,HDCJ,CJx,HCJx,CJxx,HCJxx,DCJx,HDCJx,HCB,HDCB,HbDCB
double complex,dimension(ex(3)) :: fCTheta,CThetax
real*8,dimension(ex(3)) :: KK,KKx,HKK,HKKx,Hbeta,betax,Hbetax,HW,Wx,HWx
double complex :: Theta_rhs,Theta_rhs_o
real*8 :: dR
!!! sanity check
dR = sum(RJ)+sum(IJ)+sum(RU)+sum(IU)+sum(beta)+sum(RB)+sum(IB) + &
sum(Rnu)+sum(Inu)+sum(Rk)+sum(Ik)+sum(RTheta)+sum(ITheta)
if(dR.ne.dR) then
if(sum(RJ).ne.sum(RJ))write(*,*)"NullEvol_Theta: find NaN in RJ"
if(sum(IJ).ne.sum(IJ))write(*,*)"NullEvol_Theta: find NaN in IJ"
if(sum(RU).ne.sum(RU))write(*,*)"NullEvol_Theta: find NaN in RU"
if(sum(IU).ne.sum(IU))write(*,*)"NullEvol_Theta: find NaN in IU"
if(sum(beta).ne.sum(beta))write(*,*)"NullEvol_Theta: find NaN in beta"
if(sum(RB).ne.sum(RB))write(*,*)"NullEvol_Theta: find NaN in RB"
if(sum(IB).ne.sum(IB))write(*,*)"NullEvol_Theta: find NaN in IB"
if(sum(Rnu).ne.sum(Rnu))write(*,*)"NullEvol_Theta: find NaN in Rnu"
if(sum(Inu).ne.sum(Inu))write(*,*)"NullEvol_Theta: find NaN in Inu"
if(sum(Rk).ne.sum(Rk))write(*,*)"NullEvol_Theta: find NaN in Rk"
if(sum(Ik).ne.sum(Ik))write(*,*)"NullEvol_Theta: find NaN in Ik"
if(sum(RTheta).ne.sum(RTheta))write(*,*)"NullEvol_Theta: find NaN in RTheta"
if(sum(ITheta).ne.sum(ITheta))write(*,*)"NullEvol_Theta: find NaN in ITheta"
gont = 1
return
endif
dR = R(2) - R(1)
CU = dcmplx(RU,IU)
CB = dcmplx(RB,IB)
CJ = dcmplx(RJ,IJ)
do k=1,ex(3)
call derivs_eth(ex(1:2),crho,sigma,CU(:,:,k),DCU(:,:,k),1,1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
call derivs_eth(ex(1:2),crho,sigma,CU(:,:,k),bDCU(:,:,k),1,-1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
call derivs_eth(ex(1:2),crho,sigma,CB(:,:,k),DCB(:,:,k),1,1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
call derivs_eth(ex(1:2),crho,sigma,CB(:,:,k),bDCB(:,:,k),1,-1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
call derivs_eth(ex(1:2),crho,sigma,CJ(:,:,k),DCJ(:,:,k),2,1, &
quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k))
enddo
do j=ghost_width+1,ex(2)-ghost_width
do i=ghost_width+1,ex(1)-ghost_width
CTheta0 = dcmplx(RTheta(i,j,1),ITheta(i,j,1))
fCTheta = dcmplx(RTheta(i,j,:),ITheta(i,j,:))
call cderivs_x(ex(3),R,fCTheta,CThetax)
Cnu = dcmplx(Rnu(i,j,:),Inu(i,j,:))
Ck = dcmplx(Rk(i,j,:),Ik(i,j,:))
call cget_half_x(ex(3),CB(i,j,:),HCB)
call cget_half_x(ex(3),DCB(i,j,:),HDCB)
call cget_half_x(ex(3),bDCB(i,j,:),HbDCB)
call cget_half_x(ex(3),Cnu,HCnu)
call cderivs_x(ex(3),R,Cnu,Cnux)
call cget_half_x(ex(3),Cnux,HCnux)
call cget_half_x(ex(3),Ck,HCk)
call rget_half_x(ex(3),beta(i,j,:),Hbeta)
call rderivs_x(ex(3),R,beta(i,j,:),betax)
call rget_half_x(ex(3),betax,Hbetax)
KK = dsqrt(1.d0+RJ(i,j,:)*RJ(i,j,:)+IJ(i,j,:)*IJ(i,j,:))
call rget_half_x(ex(3),KK,HKK)
call rderivs_x(ex(3),R,KK,KKx)
call rget_half_x(ex(3),KKx,HKKx)
call rderivs_x(ex(3),R,W,Wx)
call rget_half_x(ex(3),Wx,HWx)
call rget_half_x(ex(3),W(i,j,:),HW)
call cget_half_x(ex(3),CU(i,j,:),HCU)
call cderivs_x(ex(3),R,DCU(i,j,:),DCUx)
call cderivs_x(ex(3),R,CU(i,j,:),CUx)
call cget_half_x(ex(3),DCUx,HDCUx)
call cget_half_x(ex(3),CUx,HCUx)
call cget_half_x(ex(3),DCU(i,j,:),HDCU)
call cderivs_x(ex(3),R,bDCU(i,j,:),bDCUx)
call cget_half_x(ex(3),bDCUx,HbDCUx)
call cget_half_x(ex(3),bDCU(i,j,:),HbDCU)
call cderivs_x(ex(3),R,CJ(i,j,:),CJx)
call cdderivs_x(ex(3),R,CJ(i,j,:),CJxx)
call cget_half_x(ex(3),CJx,HCJx)
call cget_half_x(ex(3),CJxx,HCJxx)
call cderivs_x(ex(3),R,DCJ(i,j,:),DCJx)
call cget_half_x(ex(3),DCJx,HDCJx)
RTheta(i,j,1) = 0.d0
ITheta(i,j,1) = 0.d0
do k=1,ex(3)-1
! call getndxs(T,crho(i),sigma(j),R(k),beta(i,j,k),KK(k),CU(i,j,k),bDCU(i,j,k),DCU(i,j,k), &
! CB(i,j,k),DCB(i,j,k),W(i,j,k),CJ(i,j,k),DCJ(i,j,k),bDCB(i,j,k),Cnu(k),Ck(k),fCTheta(k),sst,Rmin)
! call getdxs(T,crho(i),sigma(j),R(k),betax(k),KKx(k),CUx(k),DCUx(k),bDCUx(k), &
! Wx(k),CJx(k),CJxx(k),DCJx(k),Cnux(k),CThetax(k),sst,Rmin)
RHS = Theta_rhs(R(k),Rmin,beta(i,j,k),betax(k),KK(k),KKx(k),CU(i,j,k),CUx(k),DCUx(k),bDCU(i,j,k),bDCUx(k), &
DCU(i,j,k),CB(i,j,k),DCB(i,j,k),W(i,j,k),Wx(k),CJ(i,j,k),DCJ(i,j,k), &
CJx(k),CJxx(k),DCJx(k),bDCB(i,j,k),Cnu(k),Cnux(k),Ck(k),fCTheta(k))
RHS = RHS - CThetax(k)
#if 0
if(cdabs(RHS)>1.d-9)then
#if 0
write(*,*)beta(i,j,k),KK(k),CU(i,j,k),bDCU(i,j,k),DCU(i,j,k),CB(i,j,k),DCB(i,j,k)
write(*,*)W(i,j,k),CJ(i,j,k),DCJ(i,j,k),bDCB(i,j,k),Cnu(k),Ck(k),fCTheta(k)
call getndxs(T,crho(i),sigma(j),R(k),beta(i,j,k),KK(k),CU(i,j,k),bDCU(i,j,k),DCU(i,j,k), &
CB(i,j,k),DCB(i,j,k),W(i,j,k),CJ(i,j,k),DCJ(i,j,k),bDCB(i,j,k),Cnu(k),Ck(k),fCTheta(k),sst,Rmin)
write(*,*)"VS"
write(*,*)beta(i,j,k),KK(k),CU(i,j,k),bDCU(i,j,k),DCU(i,j,k),CB(i,j,k),DCB(i,j,k)
write(*,*)W(i,j,k),CJ(i,j,k),DCJ(i,j,k),bDCB(i,j,k),Cnu(k),Ck(k),fCTheta(k)
#endif
write(*,*)betax(k),KKx(k),CUx(k),DCUx(k),bDCUx(k)
write(*,*)Wx(k),CJx(k),CJxx(k),DCJx(k),Cnux(k),CThetax(k)
call getdxs(T,crho(i),sigma(j),R(k),betax(k),KKx(k),CUx(k),DCUx(k),bDCUx(k), &
Wx(k),CJx(k),CJxx(k),DCJx(k),Cnux(k),CThetax(k),sst,Rmin)
write(*,*)"VS"
write(*,*)betax(k),KKx(k),CUx(k),DCUx(k),bDCUx(k)
write(*,*)Wx(k),CJx(k),CJxx(k),DCJx(k),Cnux(k),CThetax(k)
! write(*,*)RHS
! call check_factor(T,crho(i),sigma(j),R(k),sst,Rmin)
stop
endif
#endif
RTheta(i,j,k+1) = dreal(RHS)
ITheta(i,j,k+1) = dimag(RHS)
enddo
enddo
enddo
gont = 0
return
end function Eq_Theta_2