#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)
		grid_r(i,j,k) = r(i,j,k)
		param_dr(i,j,k) = 1
		param_dx(i,j,k) = 1
		param_dy(i,j,k) = 1
		param_dz(i,j,k) = 1
	       end do
	end do
  end do
!$OMP END PARALLEL DO

do k=1,cctk_lsh(3)
  	do j=1,cctk_lsh(2)
  	      do i=1,cctk_lsh(1)
		if( Seve_in_cube(i,j,k) == 1 ) then !manually inside the cube
			if (grid_r(i,j,k) == 0) then !at r= 0 every direction is equally r direction
			dxdr(i,j,k) = 1/3
			dydr(i,j,k) = 1/3
			dzdr(i,j,k) = 1/3
			else !projection in r direction
			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
	        else if( Seve_in_cube(i,j,k) == 0 ) then !using iJ = dxda for spherical grid
		dxdr(i,j,k) = iJ13(i,j,k)
		dydr(i,j,k) = iJ23(i,j,k)
		dzdr(i,j,k) = iJ33(i,j,k)
		end if	
	   end do
	end do
 end do

	
write(*,*) "Coordinate transformation factors for spatial derivatives have been prepared."

!$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)
		!grid_r = r(i,j,k)
		  if(grid_r(i,j,k) == 0) then
		  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 !-TOV_Gamma*press(i,j,k)/(rho(i,j,k)*rprec)*exp(2*PHI(i,j,k)-2*LAMBDA(i,j,k)) manual approximation for r=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 !exp(2*PHI(i,j,k)-2*LAMBDA(i,j,k))*(4-8*pio*exp(2*LAMBDA(i,j,k))*press(i,j,k) )manual approximation for r=0
		  else if (grid_r(i,j,k) >= (TOV_surface - 4*rprec)) then !surface interval + outer region
		  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) )
		  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) )
		  end if
	       end do
	end do
  end do
!$OMP END PARALLEL DO
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)
		  !grid_r = r(i,j,k)
		  if(grid_r(i,j,k) == 0) then
		  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)
		  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 (grid_r(i,j,k) <= TOV_surface .AND. grid_r(i,j,k) >0 ) then
		  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)/(r(i,j,k)**2) * exp(PHI(i,j,k))
		  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)
   		  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 (grid_r(i,j,k) > TOV_surface) then
		  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
               	  end if
	       end do
	end do
  end do
!$OMP END PARALLEL DO

write(*,*)"calculated initial fluid displacemnts, initial amplitude of Xi is: Xi = " , SWTNS_amplitude
 
  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, PHI, param_dx, J11, J21, J31, &
                                  J12, J22, J32, J13, J23, J33, -1_ik)
             call globalDiff_gv (cctkGH, 1, PHI, param_dy, J11, J21, J31, &
                                  J12, J22, J32, J13, J23, J33, -1_ik)
             call globalDiff_gv (cctkGH, 2, 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

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