#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, helper    
  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(*,*) "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
		  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
		  drrho(i,j,k) = 0
		  A(i,j,k) = 0 
		  B(i,j,k) = TOV_Gamma*press(i,j,k)/rho(i,j,k)*exp(2*PHI(i,j,k)-2*LAMBDA(i,j,k)) 
		  C(i,j,k) = 0 
		  Xi(i,j,k)  = (sin(pio*r(i,j,k)/TOV_surface)*cos(pio*r(i,j,k)/TOV_surface)-pio*r(i,j,k)/TOV_surface)*SWTNS_amplitude
		  Xeta(i,j,k) = 0
		  drXi(i,j,k) = pio * SWTNS_amplitude *(cos(2*pio*r(i,j,k)/TOV_surface)-1)/TOV_surface
		  H(i,j,k) = drXi(i,j,k) 
		  drH(i,j,k) = -(2*pio**2*SWTNS_amplitude*sin(2*pio*r(i,j,k)/TOV_surface))/(TOV_surface**2)
		  omega(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 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))
		  drLAMBDA(i,j,k) = (4*pio*r(i,j,k)**3 * rho(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) = -(rho(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)))
		  drrho(i,j,k) = 0 !manual approximation
		  A(i,j,k) = 0 !manual
		  B(i,j,k) = 0 !manual
		  C(i,j,k) = exp(2*PHI(i,j,k)-2*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
		  omega(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))
		  drLAMBDA(i,j,k) = (4*pio*r(i,j,k)**3 * rho(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) = -(rho(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)))
		  drrho(i,j,k) = drpress(i,j,k)*rho(i,j,k)/(TOV_Gamma*press(i,j,k))
		  A(i,j,k) = (exp(2*PHI(i,j,k)-2*LAMBDA(i,j,k))/rho(i,j,k))*(TOV_Gamma*press(i,j,k)*(drLAMBDA(i,j,k) + 3*drPHI(i,j,k)) &
				-(TOV_Gamma-1)*press(i,j,k)*drPHI(i,j,k)+TOV_Gamma*drpress(i,j,k) -2*TOV_Gamma*press(i,j,k)/r(i,j,k) )
		  B(i,j,k) = exp(2*PHI(i,j,k)-2*LAMBDA(i,j,k))*TOV_Gamma*press(i,j,k)/rho(i,j,k)
		  C(i,j,k) = exp(2*PHI(i,j,k)-2*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)  = (sin(pio*r(i,j,k)/TOV_surface)*cos(pio*r(i,j,k)/TOV_surface)-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) = pio * SWTNS_amplitude *(cos(2*pio*r(i,j,k)/TOV_surface)-1)/TOV_surface 
		  H(i,j,k) = drXi(i,j,k) 
		  drH(i,j,k) = -(2*pio**2*SWTNS_amplitude*sin(2*pio*r(i,j,k)/TOV_surface))/(TOV_surface**2)
		  omegasquare(i,j,k) = 0!(-1)*(A(i,j,k)*drXi(i,j,k)/Xi(i,j,k) + B(i,j,k)*drH(i,j,k)/Xi(i,j,k)+C(i,j,k))
		  if(omegasquare(i,j,k) <= 0) then
		  omega(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
		  omega(i,j,k) = sqrt(omegasquare(i,j,k))
		  Pi(i,j,k)  = Xi(i,j,k)*omega(i,j,k)
		  drPi(i,j,k) = drXi(i,j,k)*omega(i,j,k)
		  Xidot(i,j,k) = Xi(i,j,k)*omega(i,j,K)
   		  Hdot(i,j,k) = drXi(i,j,k)*omega(i,j,k)
		  Pidot(i,j,k) = -Xi(i,j,k)*omegasquare(i,j,k)
		  end if
   		  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"
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