#include "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Functions.h"
#include "cctk_Parameters.h"

subroutine SWTNS_init (CCTK_ARGUMENTS)
  implicit none
  DECLARE_CCTK_ARGUMENTS
  DECLARE_CCTK_FUNCTIONS
  DECLARE_CCTK_PARAMETERS
  integer :: i, j, k
  !CCTK_REAL :: TOV_R global
  CCTK_REAL :: pio, mua, mub, muc, mud    
  pio = acos(-1.0d0) 
  TOV_surface = TOV_R
  write(*,*)"Check if TOV_R is global..."
  write(*,*)"Radius = ", TOV_R
  write(*,*) "Surface = " , TOV_surface
  write(*,*) "precission limit for grid points = " , rprec
  write(*,*) "pio = " , pio
  write(*,*) "TOV_Gamma = " , TOV_Gamma
 ! write(*,*) "grid boundary = " , gridboundary
  write(*,*)"start initial"


!$OMP PARALLEL DO private(i,j,k)
do k=1,cctk_lsh(3)
  	do j=1,cctk_lsh(2)
  	      do i=1,cctk_lsh(1) !coordinate transformation factors
		  grid_r(i,j,k) = r(i,j,k)
		  if(patchnumber(i,j,k) > 0) then
		  dxdr(i,j,k) = iJ13(i,j,k)
		  dydr(i,j,k) = iJ23(i,j,k)
		  dzdr(i,j,k) = iJ33(i,j,k)
		  else if (patchnumber(i,j,k) == 0 .AND. grid_r(i,j,k) > 0) then
	          dxdr(i,j,k) = x(i,j,k)/sqrt(x(i,j,k)**2+y(i,j,k)**2+z(i,j,k)**2)
		  dydr(i,j,k) = y(i,j,k)/sqrt(x(i,j,k)**2+y(i,j,k)**2+z(i,j,k)**2)
		  dzdr(i,j,k) = z(i,j,k)/sqrt(x(i,j,k)**2+y(i,j,k)**2+z(i,j,k)**2)
		  end if
		  if(grid_r(i,j,k) == 0) then !r=0
		  dxdr(i,j,k) = 0
		  dydr(i,j,k) = 0
		  dzdr(i,j,k) = 0 !filling initial data starts here
		  PHI(i,j,k)  = TOV_PHI(i,j,k)
		  drPHI(i,j,k) = 0
		  LAMBDA(i,j,k) = 0
		  drLAMBDA(i,j,k) = 0
		  drpress(i,j,k) = 0
		  mu(i,j,k) = rho(i,j,k)*(1+eps(i,j,k))
		  Cs_square(i,j,k) = TOV_Gamma*press(i,j,k)/(mu(i,j,k)+press(i,j,k))
		  A(i,j,k) = 0
		  B(i,j,k) = Cs_square(i,j,k)*exp(2*(PHI(i,j,k)-LAMBDA(i,j,k)))  
		  C(i,j,k) = 0
		  Xi(i,j,k)  = (r(i,j,k)-TOV_surface)**2*sin(pio*r(i,j,k)/TOV_surface)*SWTNS_amplitude
		  Xeta(i,j,k) = 0
		  drXi(i,j,k) = SWTNS_amplitude*(2*(r(i,j,k)-TOV_surface)*sin(pio*r(i,j,k)/TOV_surface)&
				+(r(i,j,k)-TOV_surface)**2*cos(pio*r(i,j,k)/TOV_surface)*pio/TOV_surface)
		  H(i,j,k) = drXi(i,j,k) 
		  drH(i,j,k) = SWTNS_amplitude*(2*sin(pio*r(i,j,k)/TOV_surface)-(pio**2*(r(i,j,k)-TOV_surface)**2)/(TOV_surface**2)&
				*sin(pio*r(i,j,k)/TOV_surface)+(4*pio*(r(i,j,k)-TOV_surface))/TOV_surface*cos(pio*r(i,j,k)/TOV_surface))
		  Pi(i,j,k)  = Xi(i,j,k)
		  drPi(i,j,k) = H(i,j,k)
		  Xidot(i,j,k) = Pi(i,j,k)
   		  Hdot(i,j,k) = drPi(i,j,k)
		  Pidot(i,j,k) = A(i,j,k)*H(i,j,k) +B(i,j,k)*drH(i,j,k) + C(i,j,k)*Xi(i,j,k)
		  else if (press(i,j,k) == 0 .OR. rho(i,j,k) == 0) then !outside the star 
		  PHI(i,j,k)  = TOV_PHI(i,j,k)
		  drPHI(i,j,k) = (TOV_mr(i,j,k) + 4*pio*r(i,j,k)**3 * press(i,j,k))/(r(i,j,k)*(r(i,j,k)-2*TOV_mr(i,j,k)))
		  LAMBDA(i,j,k) = -0.5*log(1-2*TOV_mr(i,j,k)/r(i,j,k))
		  mu(i,j,k) = rho(i,j,k)*(1+eps(i,j,k))
		  Cs_square(i,j,k) = 0
		  drLAMBDA(i,j,k) = (4*pio*r(i,j,k)**3 * mu(i,j,k)-TOV_mr(i,j,k))/(r(i,j,k)*(r(i,j,k)-2*TOV_mr(i,j,k)))
		  drpress(i,j,k) = -(mu(i,j,k) + press(i,j,k))*(TOV_mr(i,j,k)+4*pio*r(i,j,k)**3*press(i,j,k))/(r(i,j,k)*(r(i,j,k)-2*TOV_mr(i,j,k)))
		   !manual approximation
		  A(i,j,k) = exp(2*(PHI(i,j,k)-LAMBDA(i,j,k)))*(Cs_square(i,j,k)*(drLAMBDA(i,j,k)+3*drPHI(i,j,k))&
				 - TOV_Gamma*(TOV_mr(i,j,k)+4*pio*r(i,j,k)**3*press(i,j,k))/(r(i,j,k)*(r(i,j,k)-2*TOV_mr(i,j,k)))&
				 -2*Cs_square(i,j,k)/r(i,j,k) )
		  B(i,j,k) = 0 !Cs=0
		  C(i,j,k) = exp(2*(PHI(i,j,k)-LAMBDA(i,j,k)))*(drPHI(i,j,k)**2+4*drPHI(i,j,k)/r(i,j,k)-8*pio*exp(2*LAMBDA(i,j,k))*press(i,j,k) )
		  Xi(i,j,k)  = 0
		  Xeta(i,j,k) = 0
		  drXi(i,j,k) = 0
		  H(i,j,k) = 0
		  drH(i,j,k) = 0
		  Pi(i,j,k)  = 0
		  drPi(i,j,k) = 0
		  Xidot(i,j,k) = 0
   		  Pidot(i,j,k) = 0
		  Hdot(i,j,k) = 0
		  else
		  PHI(i,j,k)  = TOV_PHI(i,j,k)
		  drPHI(i,j,k) = (TOV_mr(i,j,k) + 4*pio*r(i,j,k)**3 * press(i,j,k))/(r(i,j,k)*(r(i,j,k)-2*TOV_mr(i,j,k)))
		  LAMBDA(i,j,k) = -0.5*log(1-2*TOV_mr(i,j,k)/r(i,j,k))
		  mu(i,j,k) = rho(i,j,k)*(1+eps(i,j,k))
		  drpress(i,j,k) = -(mu(i,j,k) + press(i,j,k))*(TOV_mr(i,j,k)+4*pio*r(i,j,k)**3*press(i,j,k))/(r(i,j,k)*(r(i,j,k)-2*TOV_mr(i,j,k)))
		  Cs_square(i,j,k) = TOV_Gamma*press(i,j,k)/(mu(i,j,k)+press(i,j,k))
		  drLAMBDA(i,j,k) = (4*pio*r(i,j,k)**3 * mu(i,j,k)-TOV_mr(i,j,k))/(r(i,j,k)*(r(i,j,k)-2*TOV_mr(i,j,k)))
		  A(i,j,k) = exp(2*(PHI(i,j,k)-LAMBDA(i,j,k)))*(Cs_square(i,j,k)*(drLAMBDA(i,j,k)+3*drPHI(i,j,k))&
				 - TOV_Gamma*(TOV_mr(i,j,k)+4*pio*r(i,j,k)**3*press(i,j,k))/(r(i,j,k)*(r(i,j,k)-2*TOV_mr(i,j,k)))&
				 -2*Cs_square(i,j,k)/r(i,j,k) )
		  B(i,j,k) = Cs_square(i,j,k)*exp(2*(PHI(i,j,k)-LAMBDA(i,j,k)))
		  C(i,j,k) = exp(2*(PHI(i,j,k)-LAMBDA(i,j,k)))*(drPHI(i,j,k)**2+4*drPHI(i,j,k)/r(i,j,k)-8*pio*exp(2*LAMBDA(i,j,k))*press(i,j,k) )
		  Xi(i,j,k)  = (r(i,j,k)-TOV_surface)**2*sin(pio*r(i,j,k)/TOV_surface)*SWTNS_amplitude
		  Xeta(i,j,k) = Xi(i,j,k)*exp(LAMBDA(i,j,k)+PHI(i,j,k))/(r(i,j,k)**2)
		  drXi(i,j,k) = SWTNS_amplitude*(2*(r(i,j,k)-TOV_surface)*sin(pio*r(i,j,k)/TOV_surface)&
				+(r(i,j,k)-TOV_surface)**2*cos(pio*r(i,j,k)/TOV_surface)*pio/TOV_surface)
		  H(i,j,k) = drXi(i,j,k) 
		  drH(i,j,k) = SWTNS_amplitude*(2*sin(pio*r(i,j,k)/TOV_surface)-(pio**2*(r(i,j,k)-TOV_surface)**2)/(TOV_surface**2)&
				*sin(pio*r(i,j,k)/TOV_surface)+(4*pio*(r(i,j,k)-TOV_surface))/TOV_surface*cos(pio*r(i,j,k)/TOV_surface))
		  Pi(i,j,k)  = Xi(i,j,k)
		  drPi(i,j,k) = H(i,j,k)
		  Xidot(i,j,k) = Pi(i,j,k)
   		  Hdot(i,j,k) = drPi(i,j,k)
		  Pidot(i,j,k) = A(i,j,k)*H(i,j,k) +B(i,j,k)*drH(i,j,k) + C(i,j,k)*Xi(i,j,k)
		  end if
	       end do
	end do
  end do
!$OMP END PARALLEL DO

write(*,*) "Coordinate transformation factors for spatial derivatives have been prepared for Thornburg04 and Thornburg04nc."
write(*,*) "Calculated metric potentials, derivatives and coefficients A,B,C"
!$OMP PARALLEL DO private(i,j,k)
do k=1,cctk_lsh(3)
  	do j=1,cctk_lsh(2)
  	      do i=1,cctk_lsh(1)
		if(grid_r(i,j,k) <= smoothrad) then
		A(i,j,k) = exp(2*(PHI(i,j,k)-LAMBDA(i,j,k)))*(Cs_square(i,j,k)*(drLAMBDA(i,j,k)+3*drPHI(i,j,k))&
				 - TOV_Gamma*(TOV_mr(i,j,k)+4*pio*smoothrad**3*press(i,j,k))/(smoothrad*(smoothrad-2*TOV_mr(i,j,k)))&
				 -2*Cs_square(i,j,k)/smoothrad )
		end if
	       end do
	end do
  end do
!$OMP END PARALLEL DO

write(*,*) "smothed out A at r=0 area with a radius of ", smoothrad  

write(*,*)"done initial!"

call SWTNS_drparainit(CCTK_ARGUMENTS)
 	  
end subroutine SWTNS_init 

subroutine SWTNS_drparainit (CCTK_ARGUMENTS)
implicit none
  DECLARE_CCTK_ARGUMENTS
  DECLARE_CCTK_FUNCTIONS
  DECLARE_CCTK_PARAMETERS
  
   CCTK_INT, parameter :: izero = 0
 !  CCTK_INT :: i, j, k
       
   integer, parameter :: ik = kind (izero)
   
   integer :: na, nb, nc
   
   CCTK_REAL, dimension(:,:,:), allocatable &
        :: param_dx, param_dy, param_dz
	   

   na = cctk_lsh(1); nb = cctk_lsh(2); nc = cctk_lsh(3)

 allocate( param_dx(na,nb,nc), &
	   param_dy(na,nb,nc), &
	   param_dz(na,nb,nc) )

             
	!check derivative algorithem with drPHI
	
	     call globalDiff_gv (cctkGH, 0_ik, PHI, param_dx, J11, J21, J31, &
                                  J12, J22, J32, J13, J23, J33, -1_ik)
             call globalDiff_gv (cctkGH, 1_ik, PHI, param_dy, J11, J21, J31, &
                                  J12, J22, J32, J13, J23, J33, -1_ik)
             call globalDiff_gv (cctkGH, 2_ik, PHI, param_dz, J11, J21, J31, &
                                  J12, J22, J32, J13, J23, J33, -1_ik)

		param_dr = dxdr * param_dx + dydr * param_dy + dzdr * param_dz
		deriv_error = drPHI - param_dr

deallocate (param_dx, param_dy,param_dz)
  !write(*,*) "Calculated drP"
end subroutine SWTNS_drparainit