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

subroutine SWTNS_calc_rhs (CCTK_ARGUMENTS)
  implicit none
  DECLARE_CCTK_ARGUMENTS
  DECLARE_CCTK_FUNCTIONS
  DECLARE_CCTK_PARAMETERS
  
   CCTK_INT, parameter :: izero = 0
   CCTK_INT :: i, j, k
   CCTK_REAL :: pio
    
   integer, parameter :: ik = kind (izero)
   
   integer :: na, nb, nc
   
   CCTK_REAL, dimension(:,:,:), allocatable &
        :: TOVold, dxXi, dyXi, dzXi, &
            dxPi, dyPi, dzPi, dxH, dyH, dzH

   na = cctk_lsh(1); nb = cctk_lsh(2); nc = cctk_lsh(3)
 !   write(*,*)"allocate RHS..."
  
    !pio=3.14159265359
    !write(*,*) "Surface = " , TOV_surface
allocate( dxXi(na, nb, nc), &
           dyXi(na, nb, nc), &
           dzXi(na, nb, nc), &
           dxPi(na, nb, nc), &
           dyPi(na, nb, nc), &
           dzPi(na, nb, nc), &
           dxH(na, nb, nc), &
           dyH(na, nb, nc), &
           dzH(na, nb, nc), &
	   TOVold(na,nb,nc)  )

!Updating TOV_surface 
 outerloop: 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) > (TOV_surface -rprec) .AND. grid_r(i,j,k) < (TOV_surface + rprec) ) then
		 TOV_surface = TOV_surface + Xi(i,j,k)
		 !write(*,*) "i = ", i
!		 write(*,*) "j ", j
!		 write(*,*) "k = ", k
!		 write(*,*) "r = ", r(i,j,K)
!		 write(*,*) "Xi(i,j,k) = ", Xi(i,j,k)	
!                write(*,*) "TOV_surface = ", TOV_surface
		exit outerloop
		end if
	      end do
	end do 
  end do outerloop

!Boundary Xi(r=0) = 0
do k = 1, cctk_lsh(3)	   
     do j = 1, cctk_lsh(2)
        do i = 1, cctk_lsh(1)
	   if (grid_r(i,j,k) <= 4*rprec ) then
		Xi(i,j,k) = 0
		Xeta(i,j,k) = 0
	   else
		Xeta(i,j,k) = Xi(i,j,k)*exp(PHI(i,j,k))/(r(i,j,k)**2)
	   end if
   	end do
    end do
end do
	
!drXi
           call globalDiff_gv (cctkGH, 0_ik, Xi, dxXi, J11, J21, J31, &
                                  J12, J22, J32, J13, J23, J33,-1_ik)
           call globalDiff_gv (cctkGH, 1_ik, Xi, dyXi, J11, J21, J31, &
                                  J12, J22, J32, J13, J23, J33, -1_ik)
           call globalDiff_gv (cctkGH, 2_ik, Xi, dzXi, J11, J21, J31, &
                                  J12, J22, J32, J13, J23, J33, -1_ik)
	   
	   drXi = dxXi*dxdr + dyXi*dydr + dzXi*dzdr

!Boundary Xi'(r=R) = 0
do k = 1, cctk_lsh(3)	   
     do j = 1, cctk_lsh(2)
        do i = 1, cctk_lsh(1)
	   if (grid_r(i,j,k) >= (TOV_surface - 4*rprec) .AND. grid_r(i,j,k) <= TOV_surface ) then
		drXi(i,j,k) = 0
		H(i,j,k) = 0
	   end if
   	end do
    end do
end do          
		  
!Ableitung drPi
           call globalDiff_gv (cctkGH, 0_ik, Pi, dxPi, J11, J21, J31, &
                                  J12, J22, J32, J13, J23, J33,-1_ik)
           call globalDiff_gv (cctkGH, 1_ik, Pi, dyPi, J11, J21, J31, &
                                  J12, J22, J32, J13, J23, J33, -1_ik)
           call globalDiff_gv (cctkGH, 2_ik, Pi, dzPi, J11, J21, J31, &
                                  J12, J22, J32, J13, J23, J33, -1_ik)
	  
	   drPi = dxPi*dxdr + dyPi*dydr + dzPi*dzdr
             
!Ableitung drH
           call globalDiff_gv (cctkGH, 0_ik, H, dxH, J11, J21, J31, &
                                  J12, J22, J32, J13, J23, J33,-1_ik)
           call globalDiff_gv (cctkGH, 1_ik, H, dyH, J11, J21, J31, &
                                  J12, J22, J32, J13, J23, J33, -1_ik)
           call globalDiff_gv (cctkGH, 2_ik, H, dzH, J11, J21, J31, &
                                  J12, J22, J32, J13, J23, J33, -1_ik)
           
	  drH = dxH*dxdr + dyH*dydr + dzH*dzdr         

!!! Wave Eq.
         
	  
do k = 1, cctk_lsh(3)	   
     do j = 1, cctk_lsh(2)
        do i = 1, cctk_lsh(1)
	   if (grid_r(i,j,k) <= TOV_surface) then !inside the star
	   	 Xidot(i,j,k) = Pi(i,j,k)
	  	 Hdot(i,j,k) = drPi(i,j,k)
	  	 Pidot(i,j,k) = A(i,j,k)*drXi(i,j,k) + B(i,j,k)*drH(i,j,k) + C(i,j,k)*Xi(i,j,k)
	   else if (grid_r(i,j,k) > TOV_surface) then !outside the star
        	Xi(i,j,k) = 0
		Xeta(i,j,k) = 0
		H(i,j,k) = 0
		Pi(i,j,k) = 0		
		drXi(i,j,k) = 0
		drH(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
	
deallocate ( dxXi, dyXi, dzXi, dxPi, dyPi, dzPi, dxH, dyH, dzH, TOVold )

call SWTNS_drpara(CCTK_ARGUMENTS)
  
end subroutine SWTNS_calc_rhs

subroutine SWTNS_drpara (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, PHI, param_dr, J11, J21, J31, &
                                  J12, J22, J32, iJ13, iJ23, iJ33, -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_drpara