#include "macrodef.fh"

subroutine get_RT_parameters(m0o,Pp0o,Pm0o,apo,amo,bpo,bmo,cpo,cmo,gamo)
implicit none
real*8,intent(out) :: m0o,Pp0o,Pm0o,apo,amo,bpo,bmo,cpo,cmo,gamo

real*8,parameter::m0=1.d0,Pp0=1.d0,Pm0=1.d0,ap=1.d0,am=1.d0
real*8,parameter::bp=0.d0,bm=0.d0,cp=0.d0,cm=0.d0
real*8,parameter::gam=0.5d0

m0o = m0
Pp0o = Pp0
Pm0o = Pm0
apo = ap
amo = am
bpo = bp
bmo = bm
cpo = cp
cmo = cm
gamo = gam
end subroutine get_RT_parameters
!!!---------------------------------------------------------------------------------------------
  function boostbhP(P0,gam,a,b,c,gt,gp)  result(gont)
  implicit none

!~~~~~~> Input parameters:

  real*8, intent(in ):: P0,gam,a,b,c,gt,gp

  real*8::gont

  gont = dcosh(gam)+a*dsinh(gam)*dcos(gt)+dsinh(gam)*dsin(gt)*(b*dcos(gp)+c*dsin(gp))

  gont = P0*gont

  end function boostbhP
!!!!-------------------------------------------------------------------------------------------
#if 1
!! RT ID
subroutine get_initial_null2(ex,crho,sigma,XX,g22,g23,g33,sst,Rmin)
implicit none
! argument variables
integer, intent(in ):: ex(1:3),sst
real*8,intent(in) :: Rmin
double precision,intent(in),dimension(ex(1))::crho
double precision,intent(in),dimension(ex(2))::sigma
double precision,intent(in),dimension(ex(3))::XX
real*8,dimension(ex(1),ex(2),ex(3)),intent(out)::g22,g23,g33

double precision,dimension(ex(3))::R
real*8 :: sr,ss,cr,cs
real*8 :: sr2,ss2,cr2,cs2
integer :: i,j,k
real*8 :: ggr,tgrho,tgsigma
real*8 ::x,y,z,gr,gt,gp
real*8,dimension(ex(1),ex(2),ex(3))::P

real*8 :: PI

real*8 :: m0,Pp0,Pm0,ap,am,bp,bm,cp,cm,gam

real*8::boostbhP

call get_RT_parameters(m0,Pp0,Pm0,ap,am,bp,bm,cp,cm,gam)

R = XX*Rmin/(1-XX)

PI = dacos(-1.d0)

  do i=1,ex(1)
  do j=1,ex(2)
  do k=1,ex(3)
     sr = dsin(crho(i))
     ss = dsin(sigma(j))
     cr = dcos(crho(i))
     cs = dcos(sigma(j))
     sr2 = sr*sr
     ss2 = ss*ss
     cr2 = cr*cr
     cs2 = cs*cs

     g22(i,j,k) = 1-sr2*ss2
     g22(i,j,k) = 1/g22(i,j,k)/g22(i,j,k)

     g23(i,j,k) = -sr*cr*ss*cs*g22(i,j,k)
     g33(i,j,k) = cr2*g22(i,j,k)
     g22(i,j,k) = cs2*g22(i,j,k)

! we want g_AB/r^2 instead of g_AB     
!     g22(i,j,k) = R(k)*R(k)*g22(i,j,k)
!     g23(i,j,k) = R(k)*R(k)*g23(i,j,k)
!     g33(i,j,k) = R(k)*R(k)*g33(i,j,k)

! here fake global coordinate is enough  
    gr = 1.d0
    tgrho = dtan(crho(i))
    tgsigma = dtan(sigma(j))
    select case (sst)
    case (0)
      z = gr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
      x = z*tgrho
      y = z*tgsigma
    case (1)
      z = -gr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
      x = z*tgrho
      y = z*tgsigma
    case (2)
      x = gr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
      y = x*tgrho
      z = x*tgsigma
    case (3)
      x = -gr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
      y = x*tgrho
      z = x*tgsigma
    case (4)
      y = gr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
      x = y*tgrho
      z = y*tgsigma
    case (5)
      y = -gr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
      x = y*tgrho
      z = y*tgsigma
    case default
      write(*,*) "get_initial_null2: not recognized sst = ",sst
      return
    end select
    gt = dacos(z/gr)
    gp = datan2(y,x)

    P(i,j,k) = 1/(1/dsqrt(boostbhP(Pp0,gam,ap,bp,cp,gt,gp))+1/dsqrt(boostbhP(Pm0,gam,am,bm,cm,gt,gp)))**2

  enddo
  enddo
  enddo

  g22 = g22/P**2
  g23 = g23/P**2
  g33 = g33/P**2

return

end subroutine get_initial_null2
#else
!! fake RT for test
subroutine get_initial_null2(ex,crho,sigma,XX,g22,g23,g33,sst,Rmin)
implicit none
! argument variables
integer, intent(in ):: ex(1:3),sst
real*8,intent(in) :: Rmin
double precision,intent(in),dimension(ex(1))::crho
double precision,intent(in),dimension(ex(2))::sigma
double precision,intent(in),dimension(ex(3))::XX
real*8,dimension(ex(1),ex(2),ex(3)),intent(out)::g22,g23,g33

double precision,dimension(ex(3))::R
real*8 :: sr,ss,cr,cs
real*8 :: sr2,ss2,cr2,cs2
integer :: i,j,k
real*8 :: ggr,tgrho,tgsigma
real*8 ::x,y,z,gr,gt,gp
real*8,dimension(ex(1),ex(2),ex(3))::P

real*8 :: PI

real*8 :: m0,Pp0,Pm0,ap,am,bp,bm,cp,cm,gam

real*8::boostbhP

call get_RT_parameters(m0,Pp0,Pm0,ap,am,bp,bm,cp,cm,gam)

R = XX*Rmin/(1-XX)

PI = dacos(-1.d0)

  do i=1,ex(1)
  do j=1,ex(2)
  do k=1,ex(3)
     sr = dsin(crho(i))
     ss = dsin(sigma(j))
     cr = dcos(crho(i))
     cs = dcos(sigma(j))
     sr2 = sr*sr
     ss2 = ss*ss
     cr2 = cr*cr
     cs2 = cs*cs

     g22(i,j,k) = 1-sr2*ss2
     g22(i,j,k) = 1/g22(i,j,k)/g22(i,j,k)

     g23(i,j,k) = -sr*cr*ss*cs*g22(i,j,k)
     g33(i,j,k) = cr2*g22(i,j,k)
     g22(i,j,k) = cs2*g22(i,j,k)

! we want g_AB/r^2 instead of g_AB     
!     g22(i,j,k) = R(k)*R(k)*g22(i,j,k)
!     g23(i,j,k) = R(k)*R(k)*g23(i,j,k)
!     g33(i,j,k) = R(k)*R(k)*g33(i,j,k)

! here fake global coordinate is enough  
    gr = 1.d0
    tgrho = dtan(crho(i))
    tgsigma = dtan(sigma(j))
    select case (sst)
    case (0)
      z = gr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
      x = z*tgrho
      y = z*tgsigma
    case (1)
      z = -gr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
      x = z*tgrho
      y = z*tgsigma
    case (2)
      x = gr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
      y = x*tgrho
      z = x*tgsigma
    case (3)
      x = -gr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
      y = x*tgrho
      z = x*tgsigma
    case (4)
      y = gr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
      x = y*tgrho
      z = y*tgsigma
    case (5)
      y = -gr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
      x = y*tgrho
      z = y*tgsigma
    case default
      write(*,*) "get_initial_null2: not recognized sst = ",sst
      return
    end select
    gt = dacos(z/gr)
    gp = datan2(y,x)

    P(i,j,k) = 1/(1/dsqrt(boostbhP(Pp0,gam,ap,bp,cp,gt,gp))+1/dsqrt(boostbhP(Pm0,gam,am,bm,cm,gt,gp)))**2

  enddo
  enddo
  enddo

  g22 = P

return

end subroutine get_initial_null2
#endif
!!------------------------------------------------------------------------------------------------------------
subroutine std_covdiff(rho,sigma,fs,fr,fss,frr,frs,covf)
implicit none
! argument variables
real*8,intent(in) :: rho,sigma,fs,fr,fss,frr,frs
real*8,intent(out):: covf

real*8 :: t1,t2,t3,t4,t5,t6,t7,t8,t11,t12,t13,t15,t16,t19,t20
real*8 :: t27,t28,t29,t32,t33,t34,t38,t39,t51,t54,t55,t58,t59
real*8 :: t62,t71,t72,t88,t90,t91,t92,t93,t94,t95,t97,t98,t99
real*8 :: t100,t104,t107,t108,t109,t112,t113,t117,t118,t121,t128,t132,t133,t136,t137,t140,t141
real*8 :: t144,t152,t153,t154,t155,t160,t166,t169,t172,t175,t178,t181,t187,t199,t204,t205,t208
real*8 :: t209,t216,t217,t223,t226,t227,t243,t250,t256,t267,t276,t284,t287,t290,t301,t303,t306
real*8 :: t307,t310,t313,t314,t316,t319,t323,t326,t329,t338,t346,t356,t359,t368,t371,t376,t377
real*8 :: t380,t385,t387,t391,t394,t398,t401,t404,t407,t412,t415,t420,t427,t450,t451,t456,t459
real*8 :: t486,t487,t511,t516,t522,t532,t537,t546,t575,t586,t591,t595,t599,t295,t298

      t1 = cos(sigma);
      t2 = t1*t1;
      t3 = t2*t2;
      t4 = t3*t2;
      t5 = t4*fss;
      t6 = 2.0*sigma;
      t7 = cos(t6);
      t8 = t7*t7;
      t11 = cos(rho);
      t12 = t11*t11;
      t13 = t12*t11;
      t15 = sin(rho);
      t16 = fr*t15;
      t19 = t12*t12;
      t20 = t19*t13;
      t27 = t19*t11;
      t28 = t27*fr;
      t29 = t15*t8;
      t32 = 2.0*rho;
      t33 = cos(t32);
      t34 = t33*t33;
      t38 = t11*fr;
      t39 = t15*t3;
      t51 = t19*frr;
      t54 = t19*t12;
      t55 = t54*frr;
      t58 = -2.0*t5*t8-8.0*t2*t13*t16-64.0*t3*t20*t16+32.0*t4*t20*t16+4.0*t28*t29 &
            +4.0*t28*t15*t34-4.0*t38*t39+8.0*t3*t13*t16+32.0*t4*t13*t16+8.0*t2*t27*t16 &
            +4.0*t51*t2-2.0*t55*t34;
      t59 = t3*fss;
      t62 = t12*frr;
      t71 = t19*t19;
      t72 = t71*frr;
      t88 = -32.0*t59*t54+2.0*t62*t3-2.0*t55*t8+64.0*t55*t4-2.0*t5*t34-32.0*t72*t2 &
            +64.0*t72*t3-32.0*t72*t4-4.0*t55*t2-4.0*t51*t3-62.0*t55*t3+60.0*t3*t27*t16;
      t90 = sin(t32);
      t91 = sin(t6);
      t92 = t90*t91;
      t93 = t92*frs;
      t94 = t3*t8;
      t95 = t94*t34;
      t97 = t3*t1;
      t98 = t97*fs;
      t99 = sin(sigma);
      t100 = t98*t99;
      t104 = t2*fss;
      t107 = t54*fr;
      t108 = t90*t33;
      t109 = t108*t2;
      t112 = t19*t8;
      t113 = t112*t34;
      t117 = t12*fr;
      t118 = t108*t3;
      t121 = t12*t3;
      t128 = t19*t2;
      t132 = t3*t3;
      t133 = t132*fss;
      t136 = t93*t95-4.0*t100-32.0*t51*t4+2.0*t104*t19+8.0*t107*t109+t93*t113-62.0*t5*t19 &
            -4.0*t117*t118+2.0*t93*t121*t8+32.0*t2*t20*t16+2.0*t93*t128*t8-32.0*t133*t12;
      t137 = t3*t19;
      t140 = t2*t8;
      t141 = t140*t34;
      t144 = t8*t34;
      t152 = t107*t91;
      t153 = t90*t99;
      t154 = t2*t1;
      t155 = t153*t154;
      t160 = t33*t3*t8;
      t166 = frs*t19;
      t169 = t19*fr;
      t172 = frs*t3;
      t175 = t107*t90;
      t178 = -t93*t137*t8-4.0*t55*t141-2.0*t93*t128*t144+2.0*t62*t95+t93*t137*t144+16.0*t152*t155 &
             +4.0*t117*t90*t160+4.0*t107*t108*t8-t92*t166*t8+8.0*t169*t118-t92*t172*t34+4.0*t175*t160;
      t181 = t169*t90;
      t187 = t33*t2*t8;
      t199 = frs*t2;
      t204 = fs*t3*t154;
      t205 = t19*t99;
      t208 = fs*t154;
      t209 = t54*t99;
      t216 = -8.0*t181*t160-4.0*t107*t118-8.0*t175*t187-t93*t137*t34+4.0*t51*t141+2.0*t93*t128*t34 &
             -4.0*t51*t95+2.0*t92*t199*t12-64.0*t204*t205+32.0*t208*t209-64.0*t98*t209+32.0*t204*t209;
      t217 = t99*t8;
      t223 = t1*fs;
      t226 = t4*fs;
      t227 = t91*t7;
      t243 = t12*t99;
      t250 = 4.0*t98*t217+4.0*t98*t99*t34-4.0*t223*t205-4.0*t226*t227-64.0*t4*t27*t16+2.0*t93*t121*t34 &
            -t92*t166*t34-2.0*t93*t121*t144+8.0*t208*t205+8.0*t98*t243+60.0*t98*t205+32.0*t204*t243;
      t256 = t2*t34;
      t267 = t3*t34;
      t276 = -8.0*t208*t243+t92*t172+t92*t166-4.0*t51*t256-4.0*t51*t140+4.0*t51*t94+2.0*t55*t144 &
             -2.0*t62*t94-2.0*t62*t267-2.0*t55*t94+4.0*t55*t140+4.0*t51*t267;
      t284 = fs*t91*t33;
      t287 = t243*t34;
      t290 = t205*t34;
      t295 = fs*t27*t15;
      t298 = t92*t4;
      t301 = t92*t3;
      t303 = fs*t13*t15;
      t306 = t208*t99;
      t307 = t144*t12;
      t310 = t227*t19;
      t313 = t3*fs;
      t314 = t313*t91;
      t316 = t7*t19*t34;
      t319 = 4.0*t55*t256-2.0*t55*t267+2.0*t55*t95-32.0*t137*t284-8.0*t98*t287+4.0*t98*t290 &
            -8.0*t92*t2*t295-8.0*t298*t295-16.0*t301*t303-8.0*t306*t307-4.0*t226*t310-8.0*t314*t316;
      t323 = t217*t12;
      t326 = t227*t12;
      t329 = t226*t91;
      t338 = t7*t12*t34;
      t346 = t2*fs;
      t356 = 8.0*t208*t323+8.0*t226*t326+4.0*t329*t316+8.0*t208*t287+4.0*t226*t227*t34+8.0*t314*t338 &
           -16.0*t92*t199*t54-8.0*t329*t338-4.0*t346*t310+8.0*t313*t310-8.0*t313*t326+4.0*t346*t91*t316;
      t359 = t205*t8;
      t368 = t153*t97;
      t371 = t169*t91;
      t376 = t13*fr;
      t377 = t29*t2;
      t380 = t376*t15;
      t385 = t4*t12;
      t387 = fr*t7*t90;
      t391 = t15*t2*t34;
      t394 = 8.0*t181*t187+4.0*t98*t359-8.0*t169*t109-8.0*t152*t153*t1-8.0*t117*t91*t368+16.0*t371*t368 &
            -8.0*t152*t368+8.0*t376*t377-8.0*t380*t141-16.0*t371*t155-16.0*t385*t387+8.0*t376*t391;
      t398 = t4*t54;
      t401 = t3*t54;
      t404 = t2*t54;
      t407 = t4*t19;
      t412 = t39*t8;
      t415 = t28*t15;
      t420 = t39*t34;
      t427 = -32.0*t137*t387-16.0*t398*t387+32.0*t401*t387-16.0*t404*t387+32.0*t407*t387-8.0*t28*t377 &
             -8.0*t376*t412-4.0*t415*t95+4.0*t28*t412+4.0*t38*t420+4.0*t28*t420-8.0*t28*t391;
      t450 = t2*t12;
      t451 = t450*t34;
      t456 = -8.0*t376*t420+8.0*t415*t141+8.0*t380*t95-4.0*t28*t29*t34+4.0*t38*t412-8.0*t208*t290 &
             +32.0*t92*t172*t54+32.0*t401*t284-4.0*t100*t113+4.0*t223*t290-2.0*t93*t451-16.0*t398*t284;
      t459 = frs*t4;
      t486 = -16.0*t92*t459*t12+32.0*t407*t284-4.0*t98*t217*t34+2.0*t55-16.0*t385*t284-16.0*t404*t284 &
             -2.0*t104*t112-8.0*t208*t359-8.0*t98*t323-t92*t172*t8+4.0*t223*t359+32.0*t92*t459*t19;
      t487 = t19*t34;
      t511 = t140*t12;
      t516 = 4.0*t59*t487+8.0*t306*t113+8.0*t100*t307-2.0*t5*t487-4.0*t223*t99*t113+16.0*t298*t303 &
            -8.0*t298*fs*t11*t15+16.0*t301*t295+2.0*t104*t113+4.0*t59*t307-2.0*t93*t511-4.0*t59*t113;
      t522 = t8*t12;
      t532 = t144*t450;
      t537 = t12*t34;
      t546 = 2.0*t5*t113-4.0*t5*t307-4.0*t59*t522+2.0*t5-31.0*t92*t172*t19-2.0*t92*t199*t19+2.0*t93*t532 &
            -2.0*t5*t112+4.0*t5*t537-4.0*t59*t537+2.0*t5*t144+4.0*t5*t522;
      t575 = 4.0*t59*t112-2.0*t104*t487-4.0*t107*t108-4.0*t38*t15*t95-16.0*t92*t459*t54-2.0*t92*t172*t12 &
            -4.0*t5*t12+4.0*t59*t12+64.0*t5*t54+64.0*t133*t19-32.0*t133*t54-4.0*t59*t19-4.0*t415;
      t586 = t34*t34;
      t591 = t8*t8;
      t595 = 256.0*t137-32.0*t450+32.0*t451+32.0*t511-32.0*t532+1.0-2.0*t34+t586-2.0*t8+4.0*t144 &
            -2.0*t8*t586+t591-2.0*t591*t34+t591*t586;
      covf = -8.0*(t58+t88+t136+t178+t216+t250+t276+t319+t356+t394+t427+t456+t486+t516+t546+t575)/t595;

return

end subroutine std_covdiff
!!------------------------------------------------------------------------------------------------------------
!! input g_AB and Theta_AB are divided by r^2 indeed
!! input g_00 is also divided by r^2 indeed
! the output g00 is K
#if 1
subroutine get_gauge_g00_K(ex,crho,sigma,X,g22,g23,g33, &
                         Theta22,Theta23,Theta33, g00, Rmin)
implicit none
! argument variables
integer, intent(in ):: ex(1:3)
real*8,intent(in) :: Rmin
double precision,intent(in),dimension(ex(1))::crho
double precision,intent(in),dimension(ex(2))::sigma
double precision,intent(in),dimension(ex(3))::X
real*8,dimension(ex(1),ex(2),ex(3)),intent(in)::g22,g23,g33
real*8,dimension(ex(1),ex(2),ex(3)),intent(in)::Theta22,Theta23,Theta33
real*8,dimension(ex(1),ex(2),ex(3)),intent(out)::g00


double precision,dimension(ex(3))::R
real*8,dimension(ex(1),ex(2),ex(3))::det,gup22,gup23,gup33,KK

real*8 :: sr,ss,cr,cs,sr2,ss2,cr2,cs2,tg22,tg23,tg33
real*8 :: fr,fs,frr,fss,frs,covf

integer :: i,j,k

real*8 :: m0,Pp0,Pm0,ap,am,bp,bm,cp,cm,gam

call get_RT_parameters(m0,Pp0,Pm0,ap,am,bp,bm,cp,cm,gam)

R = X*Rmin/(1-X)
det = g22*g33-g23*g23
gup22 = g33/det
gup23 = -g23/det
gup33 = g22/det

  do i=1,ex(1)
  do j=1,ex(2)
  do k=1,ex(3)
     sr = dsin(crho(i))
     ss = dsin(sigma(j))
     cr = dcos(crho(i))
     cs = dcos(sigma(j))
     sr2 = sr*sr
     ss2 = ss*ss
     cr2 = cr*cr
     cs2 = cs*cs

     tg22 = 1-sr2*ss2
     tg22 = 1/tg22/tg22

     tg23 = -sr*cr*ss*cs*tg22
     tg33 = cr2*tg22
     tg22 = cs2*tg22

! ghat/(g/r^4) indeed
     det(i,j,k) = (tg22*tg33-tg23*tg23)/det(i,j,k)
   enddo
   enddo
   enddo

  do i=1,ex(1)
  do j=1,ex(2)
  do k=1,ex(3)
    
     call rderivs_x_point(ex(1),crho,det(:,j,k),fr,i)
     call rderivs_x_point(ex(2),sigma,det(i,:,k),fs,j)
     call rdderivs_xy_point(ex(1),ex(2),crho,sigma,det(:,:,k),frs,i,j)
     call rdderivs_x_point(ex(1),crho,det(:,j,k),frr,i)
     call rdderivs_x_point(ex(2),sigma,det(i,:,k),fss,j)

     call std_covdiff(crho(i),sigma(j),fs,fr,fss,frr,frs,covf)

     KK(i,j,k) = dsqrt(det(i,j,k))*(1-0.25*covf/R(k)**2)
   enddo
   enddo
   enddo

   g00 = KK

   return

end subroutine get_gauge_g00_K
! the input g00 is K
subroutine get_gauge_g00(ex,crho,sigma,X,g22,g23,g33, &
                         Theta22,Theta23,Theta33, g00, Rmin,fp)
implicit none
! argument variables
integer, intent(in ):: ex(1:3),fp
real*8,intent(in) :: Rmin
double precision,intent(in),dimension(ex(1))::crho
double precision,intent(in),dimension(ex(2))::sigma
double precision,intent(in),dimension(ex(3))::X
real*8,dimension(ex(1),ex(2),ex(3)),intent(in)::g22,g23,g33
real*8,dimension(ex(1),ex(2),ex(3)),intent(out)::Theta22,Theta23,Theta33
real*8,dimension(ex(1),ex(2),ex(3)),intent(in)::g00


double precision,dimension(ex(3))::R
real*8,dimension(ex(1),ex(2),ex(3))::det,gup22,gup23,gup33,KK

real*8 :: sr,ss,cr,cs,sr2,ss2,cr2,cs2,tg22,tg23,tg33
real*8 :: fr,fs,frr,fss,frs,covf

integer :: i,j,k

real*8 :: m0,Pp0,Pm0,ap,am,bp,bm,cp,cm,gam

call get_RT_parameters(m0,Pp0,Pm0,ap,am,bp,bm,cp,cm,gam)

R = X*Rmin/(1-X)
det = g22*g33-g23*g23
gup22 = g33/det
gup23 = -g23/det
gup33 = g22/det

  do i=1,ex(1)
  do j=1,ex(2)
  do k=1,ex(3)
     sr = dsin(crho(i))
     ss = dsin(sigma(j))
     cr = dcos(crho(i))
     cs = dcos(sigma(j))
     sr2 = sr*sr
     ss2 = ss*ss
     cr2 = cr*cr
     cs2 = cs*cs

     tg22 = 1-sr2*ss2
     tg22 = 1/tg22/tg22

     tg23 = -sr*cr*ss*cs*tg22
     tg33 = cr2*tg22
     tg22 = cs2*tg22

     Theta22(i,j,k) = tg22/6/m0
     Theta23(i,j,k) = tg23/6/m0
     Theta33(i,j,k) = tg33/6/m0
   enddo
   enddo
   enddo

  KK = g00

  if(fp == 0)then
  k = 1
  do i=1,ex(1)
  do j=1,ex(2)
    
     call rderivs_x_point(ex(1),crho,KK(:,j,k),fr,i)
     call rderivs_x_point(ex(2),sigma,KK(i,:,k),fs,j)
     call rdderivs_xy_point(ex(1),ex(2),crho,sigma,KK(:,:,k),frs,i,j)
     call rdderivs_x_point(ex(1),crho,KK(:,j,k),frr,i)
     call rdderivs_x_point(ex(2),sigma,KK(i,:,k),fss,j)

     call std_covdiff(crho(i),sigma(j),fs,fr,fss,frr,frs,covf)

     Theta22(i,j,k) = covf*Theta22(i,j,k)
     Theta23(i,j,k) = covf*Theta23(i,j,k)
     Theta33(i,j,k) = covf*Theta33(i,j,k)
   enddo
   enddo
   else
  do i=1,ex(1)
  do j=1,ex(2)
  do k=1,ex(3)
    
     call rderivs_x_point(ex(1),crho,KK(:,j,k),fr,i)
     call rderivs_x_point(ex(2),sigma,KK(i,:,k),fs,j)
     call rdderivs_xy_point(ex(1),ex(2),crho,sigma,KK(:,:,k),frs,i,j)
     call rdderivs_x_point(ex(1),crho,KK(:,j,k),frr,i)
     call rdderivs_x_point(ex(2),sigma,KK(i,:,k),fss,j)

     call std_covdiff(crho(i),sigma(j),fs,fr,fss,frr,frs,covf)

     Theta22(i,j,k) = covf*Theta22(i,j,k)
     Theta23(i,j,k) = covf*Theta23(i,j,k)
     Theta33(i,j,k) = covf*Theta33(i,j,k)
   enddo
   enddo
   enddo
   endif

   return

end subroutine get_gauge_g00
#else
subroutine get_gauge_g00_K(ex,crho,sigma,X,g22,g23,g33, &
                         Theta22,Theta23,Theta33, g00, Rmin)
implicit none
! argument variables
integer, intent(in ):: ex(1:3)
real*8,intent(in) :: Rmin
double precision,intent(in),dimension(ex(1))::crho
double precision,intent(in),dimension(ex(2))::sigma
double precision,intent(in),dimension(ex(3))::X
real*8,dimension(ex(1),ex(2),ex(3)),intent(in)::g22,g23,g33
real*8,dimension(ex(1),ex(2),ex(3)),intent(in)::Theta22,Theta23,Theta33
real*8,dimension(ex(1),ex(2),ex(3)),intent(out)::g00


double precision,dimension(ex(3))::R
real*8,dimension(ex(1),ex(2),ex(3))::det,gup22,gup23,gup33,KK

real*8 :: sr,ss,cr,cs,sr2,ss2,cr2,cs2,tg22,tg23,tg33
real*8 :: fr,fs,frr,fss,frs,covf

integer :: i,j,k

real*8 :: m0,Pp0,Pm0,ap,am,bp,bm,cp,cm,gam

call get_RT_parameters(m0,Pp0,Pm0,ap,am,bp,bm,cp,cm,gam)

R = X*Rmin/(1-X)
! g22 is P
det = dlog(g22**2)

  do i=1,ex(1)
  do j=1,ex(2)
  do k=1,ex(3)
    
     call rderivs_x_point(ex(1),crho,det(:,j,k),fr,i)
     call rderivs_x_point(ex(2),sigma,det(i,:,k),fs,j)
     call rdderivs_xy_point(ex(1),ex(2),crho,sigma,det(:,:,k),frs,i,j)
     call rdderivs_x_point(ex(1),crho,det(:,j,k),frr,i)
     call rdderivs_x_point(ex(2),sigma,det(i,:,k),fss,j)

     call std_covdiff(crho(i),sigma(j),fs,fr,fss,frr,frs,covf)

     KK(i,j,k) = covf
   enddo
   enddo
   enddo

   g00 = g22**2*(1+0.5*KK)

   return

end subroutine get_gauge_g00_K
! the input g00 is K
subroutine get_gauge_g00(ex,crho,sigma,X,g22,g23,g33, &
                         Theta22,Theta23,Theta33, g00, Rmin,fp)
implicit none
! argument variables
integer, intent(in ):: ex(1:3),fp
real*8,intent(in) :: Rmin
double precision,intent(in),dimension(ex(1))::crho
double precision,intent(in),dimension(ex(2))::sigma
double precision,intent(in),dimension(ex(3))::X
real*8,dimension(ex(1),ex(2),ex(3)),intent(in)::g22,g23,g33
real*8,dimension(ex(1),ex(2),ex(3)),intent(out)::Theta22,Theta23,Theta33
real*8,dimension(ex(1),ex(2),ex(3)),intent(in)::g00


double precision,dimension(ex(3))::R
real*8,dimension(ex(1),ex(2),ex(3))::det,gup22,gup23,gup33,KK

real*8 :: sr,ss,cr,cs,sr2,ss2,cr2,cs2,tg22,tg23,tg33
real*8 :: fr,fs,frr,fss,frs,covf

integer :: i,j,k

real*8 :: m0,Pp0,Pm0,ap,am,bp,bm,cp,cm,gam

call get_RT_parameters(m0,Pp0,Pm0,ap,am,bp,bm,cp,cm,gam)

R = X*Rmin/(1-X)

  KK = g00

  if(fp == 0)then
  k = 1
  do i=1,ex(1)
  do j=1,ex(2)
    
     call rderivs_x_point(ex(1),crho,KK(:,j,k),fr,i)
     call rderivs_x_point(ex(2),sigma,KK(i,:,k),fs,j)
     call rdderivs_xy_point(ex(1),ex(2),crho,sigma,KK(:,:,k),frs,i,j)
     call rdderivs_x_point(ex(1),crho,KK(:,j,k),frr,i)
     call rdderivs_x_point(ex(2),sigma,KK(i,:,k),fss,j)

     call std_covdiff(crho(i),sigma(j),fs,fr,fss,frr,frs,covf)

     Theta22(i,j,k) = covf
   enddo
   enddo
   else
  do i=1,ex(1)
  do j=1,ex(2)
  do k=1,ex(3)
    
     call rderivs_x_point(ex(1),crho,KK(:,j,k),fr,i)
     call rderivs_x_point(ex(2),sigma,KK(i,:,k),fs,j)
     call rdderivs_xy_point(ex(1),ex(2),crho,sigma,KK(:,:,k),frs,i,j)
     call rdderivs_x_point(ex(1),crho,KK(:,j,k),frr,i)
     call rdderivs_x_point(ex(2),sigma,KK(i,:,k),fss,j)

     call std_covdiff(crho(i),sigma(j),fs,fr,fss,frr,frs,covf)

     Theta22(i,j,k) = covf
   enddo
   enddo
   enddo
   endif

   Theta22 = -Theta22/12/m0*g22**3
   return

end subroutine get_gauge_g00
#endif
!!---------------------------------------------------------------------------
subroutine get_gauge_g00_real(ex,crho,sigma,X,g22,g23,g33, &
                         Theta22,Theta23,Theta33, g00, Rmin)
implicit none
! argument variables
integer, intent(in ):: ex(1:3)
real*8,intent(in) :: Rmin
double precision,intent(in),dimension(ex(1))::crho
double precision,intent(in),dimension(ex(2))::sigma
double precision,intent(in),dimension(ex(3))::X
real*8,dimension(ex(1),ex(2),ex(3)),intent(in)::g22,g23,g33
real*8,dimension(ex(1),ex(2),ex(3)),intent(in)::Theta22,Theta23,Theta33
real*8,dimension(ex(1),ex(2),ex(3)),intent(out)::g00


double precision,dimension(ex(3))::R
real*8,dimension(ex(1),ex(2),ex(3))::det,gup22,gup23,gup33,KK

real*8 :: sr,ss,cr,cs,sr2,ss2,cr2,cs2,tg22,tg23,tg33
real*8 :: fr,fs,frr,fss,frs,covf

integer :: i,j,k

real*8 :: m0,Pp0,Pm0,ap,am,bp,bm,cp,cm,gam

call get_RT_parameters(m0,Pp0,Pm0,ap,am,bp,bm,cp,cm,gam)

R = X*Rmin/(1-X)
det = g22*g33-g23*g23
gup22 = g33/det
gup23 = -g23/det
gup33 = g22/det

  do i=1,ex(1)
  do j=1,ex(2)
  do k=1,ex(3)
     sr = dsin(crho(i))
     ss = dsin(sigma(j))
     cr = dcos(crho(i))
     cs = dcos(sigma(j))
     sr2 = sr*sr
     ss2 = ss*ss
     cr2 = cr*cr
     cs2 = cs*cs

     tg22 = 1-sr2*ss2
     tg22 = 1/tg22/tg22

     tg23 = -sr*cr*ss*cs*tg22
     tg33 = cr2*tg22
     tg22 = cs2*tg22

! ghat/(g/r^4) indeed
     det(i,j,k) = (tg22*tg33-tg23*tg23)/det(i,j,k)
   enddo
   enddo
   enddo

  do i=1,ex(1)
  do j=1,ex(2)
  do k=1,ex(3)
    
     call rderivs_x_point(ex(1),crho,det(:,j,k),fr,i)
     call rderivs_x_point(ex(2),sigma,det(i,:,k),fs,j)
     call rdderivs_xy_point(ex(1),ex(2),crho,sigma,det(:,:,k),frs,i,j)
     call rdderivs_x_point(ex(1),crho,det(:,j,k),frr,i)
     call rdderivs_x_point(ex(2),sigma,det(i,:,k),fss,j)

     call std_covdiff(crho(i),sigma(j),fs,fr,fss,frr,frs,covf)

     KK(i,j,k) = dsqrt(det(i,j,k))*(1-0.25*covf/R(k)**2)

     g00(i,j,k) = 2*m0/R(k)**3-KK(i,j,k)/R(k)**2 &
                 -(gup22(i,j,k)*Theta22(i,j,k)+2*gup23(i,j,k)*Theta23(i,j,k)+gup33(i,j,k)*Theta33(i,j,k))/2/R(k)
   enddo
   enddo
   enddo

   return

end subroutine get_gauge_g00_real
!!------------------------------------------------------------------------------------------------------------
subroutine get_null_boundary2(ex,crho,sigma,X,g22,g23,g33, &
                         g01,p02,p03,g02,g03,Theta22,Theta23,Theta33, Rmin)
implicit none
! argument variables
integer, intent(in ):: ex(1:3)
real*8,intent(in) :: Rmin
double precision,intent(in),dimension(ex(1))::crho
double precision,intent(in),dimension(ex(2))::sigma
double precision,intent(in),dimension(ex(3))::X
real*8,dimension(ex(1),ex(2),ex(3)),intent(inout)::g22,g23,g33
real*8,dimension(ex(1),ex(2),ex(3)),intent(inout)::g01,p02,p03,g02,g03
real*8,dimension(ex(1),ex(2),ex(3)),intent(inout)::Theta22,Theta23,Theta33

#if 1
real*8 :: fact

!fact = X(1)/X(2)*((1-X(2))/(1-X(1)))
!fact = fact**2
! since we used gAB/r^2 instead of gAB, so fact = 1
fact = 1.d0

g22(:,:,1) = g22(:,:,2)*fact
g23(:,:,1) = g23(:,:,2)*fact
g33(:,:,1) = g33(:,:,2)*fact

g01(:,:,1) = -1.d0

p02(:,:,1) = 0.d0
p03(:,:,1) = 0.d0
g02(:,:,1) = 0.d0
g03(:,:,1) = 0.d0

! have done in get_gauge_g00 
!Theta22(:,:,1) = Theta22(:,:,2)*fact
!Theta23(:,:,1) = Theta23(:,:,2)*fact
!Theta33(:,:,1) = Theta33(:,:,2)*fact
#else
g01 = -1
g02 = 0
g03 = 0
#endif
return

end subroutine get_null_boundary2
!!!--------------------------------------------------------------------------------------------------------------
subroutine get_initial_null3(ex,crho,sigma,XX,g22,g23,g33,sst,Rmin)
implicit none
! argument variables
integer, intent(in ):: ex(1:3),sst
real*8,intent(in) :: Rmin
double precision,intent(in),dimension(ex(1))::crho
double precision,intent(in),dimension(ex(2))::sigma
double precision,intent(in),dimension(ex(3))::XX
real*8,dimension(ex(1),ex(2),ex(3)),intent(out)::g22,g23,g33

double precision,dimension(ex(3))::R
real*8 :: sr,ss,cr,cs
real*8 :: sr2,ss2,cr2,cs2
integer :: i,j,k
real*8 :: ggr,tgrho,tgsigma
real*8 ::x,y,z,gr,gt,gp

real*8 :: gxx,gxy,gyy,tc,ts,PI

double complex :: Zslm,II,Jr,ctp
double complex :: swtf,z220

double complex :: beta0,C1,C2,mx,my,mlx,mly
integer :: nu,m

call initial_null_paramter(beta0,C1,C2,nu,m)

  II = dcmplx(0.d0,1.d0)

R = XX*Rmin/(1-XX)

PI = dacos(-1.d0)

  do i=1,ex(1)
  do j=1,ex(2)
  do k=1,ex(3)
     sr = dsin(crho(i))
     ss = dsin(sigma(j))
     cr = dcos(crho(i))
     cs = dcos(sigma(j))
     sr2 = sr*sr
     ss2 = ss*ss
     cr2 = cr*cr
     cs2 = cs*cs

     gxx = 1-sr2*ss2
     gxx = 1/gxx/gxx

     gxy = -sr*cr*ss*cs*gxx
     gyy = cr2*gxx
     gxx = cs2*gxx
! here fake global coordinate is enough  
    gr = 1.d0
    tgrho = dtan(crho(i))
    tgsigma = dtan(sigma(j))
    tc = dsqrt((1.d0-dsin(crho(i))*dsin(sigma(j)))/2.d0)
    ts = dsqrt((1.d0+dsin(crho(i))*dsin(sigma(j)))/2.d0)
    select case (sst)
    case (0)
      z = gr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
      x = z*tgrho
      y = z*tgsigma
    case (1)
      z = -gr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
      x = z*tgrho
      y = z*tgsigma
    case (2)
      x = gr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
      y = x*tgrho
      z = x*tgsigma
    case (3)
      x = -gr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
      y = x*tgrho
      z = x*tgsigma
    case (4)
      y = gr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
      x = y*tgrho
      z = y*tgsigma
    case (5)
      y = -gr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
      x = y*tgrho
      z = y*tgsigma
    case default
      write(*,*) "get_initial_null2: not recognized sst = ",sst
      return
    end select
    gt = dacos(z/gr)
    gp = datan2(y,x)

    swtf = 2.d0*tc*ts*(ts+II*tc)/dcos(sigma(j))
    if(sst==1 .or. sst==3 .or. sst==4) swtf = dconjg(swtf)
    select case (sst)
    case (0,1)
      swtf = swtf/(dcos(gp)+II*dcos(gt)*dsin(gp))*(dcos(gt)**2+dsin(gt)**2*dcos(gp)**2)
    case (2,3)
      swtf = II*swtf*dsin(gt)
    case (4,5)
      swtf = -II*swtf*dsin(gt)
    end select

    z220 = Zslm(2,2,m,gt,gp)*swtf**2

    if(sst==1 .or. sst==3 .or. sst==4)then
      mx = 2*tc*ts*(ts-II*tc)/dcos(sigma(j))
      my = 2*tc*ts*(ts+II*tc)/dcos(crho(i))
    else
      mx = 2*tc*ts*(ts+II*tc)/dcos(sigma(j))
      my = 2*tc*ts*(ts-II*tc)/dcos(crho(i))
    endif
    mlx = gxx*mx+gxy*my
    mly = gxy*mx+gyy*my
    Jr = (2.4d1*beta0+3.d0*II*nu*C1-II*nu**3*C2)/3.6d1+C1/4.d0/R(k)-C2/1.2d1/R(k)**3
    Jr = dsqrt(dble((2-1)*2*(2+1)*(2+2)))*dreal(Jr)*z220

    ctp = Jr*mlx*mlx+dsqrt(1+abs(Jr)**2)*dconjg(mlx)*mlx
    g22(i,j,k) = dreal(ctp)
    ctp = Jr*mlx*mly+dsqrt(1+abs(Jr)**2)*dconjg(mlx)*mly
    g23(i,j,k) = dreal(ctp)
    ctp = Jr*mly*mly+dsqrt(1+abs(Jr)**2)*dconjg(mly)*mly
    g33(i,j,k) = dreal(ctp)

  enddo
  enddo
  enddo

return

end subroutine get_initial_null3
!!!--------------------------------------------------------------------------------------------------------------
subroutine get_g00_with_t(time,ex,crho,sigma,XX,g00,Rmin,sst)
implicit none
! argument variables
integer, intent(in ):: ex(1:3),sst
real*8,intent(in) :: time,Rmin
double precision,intent(in),dimension(ex(1))::crho
double precision,intent(in),dimension(ex(2))::sigma
double precision,intent(in),dimension(ex(3))::XX
real*8,dimension(ex(1),ex(2),ex(3)),intent(out)::g00

double precision,dimension(ex(3))::R
real*8 :: sr,ss,cr,cs
real*8 :: sr2,ss2,cr2,cs2
integer :: i,j,k
real*8 :: ggr,tgrho,tgsigma
real*8 ::x,y,z,gr,gt,gp

real*8 :: tc,ts,PI

double complex :: Zslm,II,Jr,Ur,Wr
double complex :: swtf,z020,z120,z220

double complex :: beta0,C1,C2
integer :: nu,m

call initial_null_paramter(beta0,C1,C2,nu,m)

  II = dcmplx(0.d0,1.d0)

R = XX*Rmin/(1-XX)

PI = dacos(-1.d0)

  do i=1,ex(1)
  do j=1,ex(2)
  do k=1,ex(3)
     sr = dsin(crho(i))
     ss = dsin(sigma(j))
     cr = dcos(crho(i))
     cs = dcos(sigma(j))
     sr2 = sr*sr
     ss2 = ss*ss
     cr2 = cr*cr
     cs2 = cs*cs

! here fake global coordinate is enough  
    gr = 1.d0
    tgrho = dtan(crho(i))
    tgsigma = dtan(sigma(j))
    tc = dsqrt((1.d0-dsin(crho(i))*dsin(sigma(j)))/2.d0)
    ts = dsqrt((1.d0+dsin(crho(i))*dsin(sigma(j)))/2.d0)
    select case (sst)
    case (0)
      z = gr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
      x = z*tgrho
      y = z*tgsigma
    case (1)
      z = -gr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
      x = z*tgrho
      y = z*tgsigma
    case (2)
      x = gr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
      y = x*tgrho
      z = x*tgsigma
    case (3)
      x = -gr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
      y = x*tgrho
      z = x*tgsigma
    case (4)
      y = gr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
      x = y*tgrho
      z = y*tgsigma
    case (5)
      y = -gr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
      x = y*tgrho
      z = y*tgsigma
    case default
      write(*,*) "get_g00_with_t: not recognized sst = ",sst
      return
    end select
    gt = dacos(z/gr)
    gp = datan2(y,x)

    swtf = 2.d0*tc*ts*(ts+II*tc)/dcos(sigma(j))
    if(sst==1 .or. sst==3 .or. sst==4) swtf = dconjg(swtf)
    select case (sst)
    case (0,1)
      swtf = swtf/(dcos(gp)+II*dcos(gt)*dsin(gp))*(dcos(gt)**2+dsin(gt)**2*dcos(gp)**2)
    case (2,3)
      swtf = II*swtf*dsin(gt)
    case (4,5)
      swtf = -II*swtf*dsin(gt)
    end select

    z020 = Zslm(0,2,m,gt,gp)
    z120 = Zslm(1,2,m,gt,gp)*swtf
    z220 = Zslm(2,2,m,gt,gp)*swtf**2

    Jr = (2.4d1*beta0+3.d0*II*nu*C1-II*nu**3*C2)/3.6d1+C1/4.d0/R(k)-C2/1.2d1/R(k)**3
    Ur = (-24*II*nu*beta0+3*nu*nu*C1-nu**4*C2)/36+2*beta0/R(k)+C1/2/R(k)**2+ &
         II*nu*C2/3/R(k)**3+C2/4/R(k)**4
    Wr = (24*II*nu*beta0-2*nu*C1+nu**4*C2)/6+ &
         (3*II*nu*C1-6*beta0-II*nu**3*C2)/3/R(k) - &
         nu**2*C2/R(k)**2+II*nu*C2/R(k)**3+C2/2/R(k)**4

    Jr = Jr*exp(II*nu*time)
    Ur = Ur*exp(II*nu*time)
    Wr = Wr*exp(II*nu*time)

    g00(i,j,k) = 2*(2*(2+1)*dsqrt(dble((2-1)*2*(2+1)*(2+2)))*dreal(Ur)**2* &
                    dreal(Jr)*dreal(z120**2*dconjg(z220))+ &
                    2*(2+1)*dsqrt(1+(2-1)*2*(2+1)*(2+2)*dreal(Jr)**2*abs(z220)**2)* &
                    dreal(Ur)**2*abs(z120)**2)-(1/R(k)**2+dreal(z020*Wr)/R(k))* &
                    exp(2*dreal(z020*beta0*exp(II*nu*time)))

  enddo
  enddo
  enddo

!if(sst==0 .and. crho(1) <-0.9 .and. sigma(1) <-0.9 .and. XX(1)<0.18182)write(*,*)"time = ",time,g00(1,1,1)

return

end subroutine get_g00_with_t
!!------------------------------------------------------------------------------------------------------------
subroutine get_null_boundary3(time,ex,crho,sigma,XX,g22,g23,g33, &
                         g01,p02,p03,g02,g03,Theta22,Theta23,Theta33, Rmin,sst)
implicit none
! argument variables
integer, intent(in ):: ex(1:3),sst
real*8,intent(in) ::time,Rmin
double precision,intent(in),dimension(ex(1))::crho
double precision,intent(in),dimension(ex(2))::sigma
double precision,intent(in),dimension(ex(3))::XX
real*8,dimension(ex(1),ex(2),ex(3)),intent(inout)::g22,g23,g33
real*8,dimension(ex(1),ex(2),ex(3)),intent(inout)::g01,p02,p03,g02,g03
real*8,dimension(ex(1),ex(2),ex(3)),intent(inout)::Theta22,Theta23,Theta33

double precision,dimension(ex(3))::R
real*8 :: sr,ss,cr,cs
real*8 :: sr2,ss2,cr2,cs2
integer :: i,j,k
real*8 :: ggr,tgrho,tgsigma
real*8 ::x,y,z,gr,gt,gp

real*8 :: gxx,gxy,gyy,tc,ts,PI

double complex :: Zslm,II,Jr,ctp,Jrp,Jrt,Ur,Urp,Wr
double complex :: swtf,z020,z120,z220

double complex :: beta0,C1,C2,mx,my,mlx,mly
integer :: nu,m

call initial_null_paramter(beta0,C1,C2,nu,m)

  II = dcmplx(0.d0,1.d0)

R = XX*Rmin/(1-XX)

PI = dacos(-1.d0)

  do i=1,ex(1)
  do j=1,ex(2)
  do k=1,ex(3)
     sr = dsin(crho(i))
     ss = dsin(sigma(j))
     cr = dcos(crho(i))
     cs = dcos(sigma(j))
     sr2 = sr*sr
     ss2 = ss*ss
     cr2 = cr*cr
     cs2 = cs*cs

     gxx = 1-sr2*ss2
     gxx = 1/gxx/gxx

     gxy = -sr*cr*ss*cs*gxx
     gyy = cr2*gxx
     gxx = cs2*gxx
! here fake global coordinate is enough  
    gr = 1.d0
    tgrho = dtan(crho(i))
    tgsigma = dtan(sigma(j))
    tc = dsqrt((1.d0-dsin(crho(i))*dsin(sigma(j)))/2.d0)
    ts = dsqrt((1.d0+dsin(crho(i))*dsin(sigma(j)))/2.d0)
    select case (sst)
    case (0)
      z = gr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
      x = z*tgrho
      y = z*tgsigma
    case (1)
      z = -gr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
      x = z*tgrho
      y = z*tgsigma
    case (2)
      x = gr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
      y = x*tgrho
      z = x*tgsigma
    case (3)
      x = -gr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
      y = x*tgrho
      z = x*tgsigma
    case (4)
      y = gr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
      x = y*tgrho
      z = y*tgsigma
    case (5)
      y = -gr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma)
      x = y*tgrho
      z = y*tgsigma
    case default
      write(*,*) "get_null_boundary3: not recognized sst = ",sst
      return
    end select
    gt = dacos(z/gr)
    gp = datan2(y,x)

    swtf = 2.d0*tc*ts*(ts+II*tc)/dcos(sigma(j))
    if(sst==1 .or. sst==3 .or. sst==4) swtf = dconjg(swtf)
    select case (sst)
    case (0,1)
      swtf = swtf/(dcos(gp)+II*dcos(gt)*dsin(gp))*(dcos(gt)**2+dsin(gt)**2*dcos(gp)**2)
    case (2,3)
      swtf = II*swtf*dsin(gt)
    case (4,5)
      swtf = -II*swtf*dsin(gt)
    end select

    z020 = Zslm(0,2,m,gt,gp)
    z120 = Zslm(1,2,m,gt,gp)*swtf
    z220 = Zslm(2,2,m,gt,gp)*swtf**2

    if(sst==1 .or. sst==3 .or. sst==4)then
      mx = 2*tc*ts*(ts-II*tc)/dcos(sigma(j))
      my = 2*tc*ts*(ts+II*tc)/dcos(crho(i))
    else
      mx = 2*tc*ts*(ts+II*tc)/dcos(sigma(j))
      my = 2*tc*ts*(ts-II*tc)/dcos(crho(i))
    endif
    mlx = gxx*mx+gxy*my
    mly = gxy*mx+gyy*my

    Jr = (2.4d1*beta0+3.d0*II*nu*C1-II*nu**3*C2)/3.6d1+C1/4.d0/R(k)-C2/1.2d1/R(k)**3
! Jrp = d Jr/d X instead of d Jr/d r    
    Jrp = -C1/4.d0/Rmin/XX(k)**2+C2/1.2d1*3/R(k)**2/Rmin/XX(k)**2
    Ur = (-24*II*nu*beta0+3*nu*nu*C1-nu**4*C2)/36+2*beta0/R(k)+C1/2/R(k)**2+ &
         II*nu*C2/3/R(k)**3+C2/4/R(k)**4
    Urp = -2*beta0/Rmin/XX(k)**2-C1/R(k)/Rmin/XX(k)**2- &
         II*nu*C2/R(k)**2/Rmin/XX(k)**2-C2/R(k)**3/Rmin/XX(k)**2
    Wr = (24*II*nu*beta0-2*nu*C1+nu**4*C2)/6+ &
         (3*II*nu*C1-6*beta0-II*nu**3*C2)/3/R(k) - &
         nu**2*C2/R(k)**2+II*nu*C2/R(k)**3+C2/2/R(k)**4

    Jr = Jr*exp(II*nu*time)
    Jrp = Jrp*exp(II*nu*time)
    Jrt = II*nu*Jr*exp(II*nu*time)
    Ur = Ur*exp(II*nu*time)
    Urp = Urp*exp(II*nu*time)
    Wr = Wr*exp(II*nu*time)

    Jr = dsqrt(dble((2-1)*2*(2+1)*(2+2)))*dreal(Jr)*z220
    Jrt = dsqrt(dble((2-1)*2*(2+1)*(2+2)))*dreal(Jrt)*z220
    Jrp = dsqrt(dble((2-1)*2*(2+1)*(2+2)))*dreal(Jrp)*z220

    g01(i,j,k) = -dexp(2*dreal(z020*beta0*exp(II*nu*time)))
#if 1
    g02(i,j,k) = -dsqrt(dble(2*(2+1)))*dreal(Ur)*dreal(mlx*(z120*dconjg(z220)* &
                  dconjg(Jr)+dsqrt(1+abs(Jr)**2)*dconjg(z120)))
    g03(i,j,k) = -dsqrt(dble(2*(2+1)))*dreal(Ur)*dreal(mly*(z120*dconjg(z220)* &
                  dconjg(Jr)+dsqrt(1+abs(Jr)**2)*dconjg(z120)))
#elif 0
    mlx = mlx/swtf
    mly = mly/swtf
    g02(i,j,k) = dreal(mlx)
    g03(i,j,k) = dreal(mly)
    !if(sst==0 .and. crho(1) <-0.9 .and. sigma(1) <-0.9 .and. XX(1)<0.18182)write(*,*)g02(i,j,k),g03(i,j,k)
    !if(crho(1) <-0.9 .and. sigma(1) <-0.9 .and. XX(1)<0.18182)write(*,*)g02(i,j,k),g03(i,j,k)
#else   
    select case (sst)
    case (0,1)
      tc =-dcos(gp)/(dcos(gt)**2*dcos(gp)**2-dcos(gt)**2-dcos(gp)**2)
      ts = dsin(gp)/(1-dsin(gt)**2*dcos(gp)**2)
    case (2,3)
      tc = 0
      ts = dcos(gp)/(dcos(gt)**2*dcos(gp)**2-dcos(gt)**2-dcos(gp)**2)
    case (4,5)
      tc = 0
      ts =-dsin(gp)/(1-dsin(gt)**2*dcos(gp)**2)
    end select
    g02(i,j,k) = gxx*tc+gxy*ts
    g03(i,j,k) = gxy*tc+gyy*ts
    !if(sst==0 .and. crho(1) <-0.9 .and. sigma(1) <-0.9 .and. XX(1)<0.18182)write(*,*)g02(i,j,k),g03(i,j,k)
    if(crho(1) <-0.9 .and. sigma(1) <-0.9 .and. XX(1)<0.18182)write(*,*)g02(i,j,k),g03(i,j,k),sst
    stop
#endif
    p02(i,j,k) = -dsqrt(dble(2*(2+1)))*dreal(Urp)*dreal(mlx*(z120*dconjg(Jr)+ &
                  dsqrt(1+abs(Jr)**2)*dconjg(z120))) &
                 -dsqrt(dble(2*(2+1)))*dreal(Ur)*dreal(mlx*(z120*dconjg(Jrp)+ &
                  abs(Jrp)*abs(Jr)/dsqrt(1+abs(Jr)**2)*dconjg(z120)))
    p03(i,j,k) = -dsqrt(dble(2*(2+1)))*dreal(Urp)*dreal(mly*(z120*dconjg(Jr)+ &
                  dsqrt(1+abs(Jr)**2)*dconjg(z120))) &
                 -dsqrt(dble(2*(2+1)))*dreal(Ur)*dreal(mly*(z120*dconjg(Jrp)+ &
                  abs(Jrp)*abs(Jr)/dsqrt(1+abs(Jr)**2)*dconjg(z120)))

    ctp = dconjg(Jr)*mlx*mlx+dsqrt(1+abs(Jr)**2)*dconjg(mlx)*mlx
    g22(i,j,k) = dreal(ctp)
    ctp = dconjg(Jr)*mlx*mly+dsqrt(1+abs(Jr)**2)*dconjg(mlx)*mly
    g23(i,j,k) = dreal(ctp)
    ctp = dconjg(Jr)*mly*mly+dsqrt(1+abs(Jr)**2)*dconjg(mly)*mly
    g33(i,j,k) = dreal(ctp)

    ctp = dconjg(Jrt)*mlx*mlx+abs(Jr)*abs(Jrt)/dsqrt(1+abs(Jr)**2)*dconjg(mlx)*mlx
    Theta22(i,j,k) = dreal(ctp)
    ctp = dconjg(Jrt)*mlx*mly+abs(Jr)*abs(Jrt)/dsqrt(1+abs(Jr)**2)*dconjg(mlx)*mly
    Theta23(i,j,k) = dreal(ctp)
    ctp = dconjg(Jrt)*mly*mly+abs(Jr)*abs(Jrt)/dsqrt(1+abs(Jr)**2)*dconjg(mly)*mly
    Theta33(i,j,k) = dreal(ctp)

  enddo
  enddo
  enddo

return

end subroutine get_null_boundary3