c  Mahc_Init_Data: Various Initial Data Routines for Mahc_Evolve
c  Copyright (C) 2001 Mark Miller

/*@@
  @file      tov_iso.F
  @date      December 1999
  @author    Mark Miller
  @desc
     Set TOV isotropic coordinates as initial data.
  @enddesc
@@*/
#include "cctk.h"
#include "cctk_Parameters.h"
#include "cctk_Arguments.h"
#include "cctk_DefineThorn.h"


/*@@
  @routine    tov_iso
  @date       December 1999
  @author     Mark Miller
  @desc
     Set TOV isotropic coordinates as initial data.
  @enddesc
  @calls
  @calledby
  @history
  @endhistory
@@*/


      subroutine tov_iso(CCTK_FARGUMENTS)
      implicit none
      DECLARE_CCTK_FARGUMENTS
      DECLARE_CCTK_PARAMETERS

c     LOCAL VARS
c     TOV SOLVER PARAMETERS
      CCTK_REAL tov_centden,rmaxtov,tov_eos_k,tov_eos_gamma,tiny
      integer ntov,tov_solver,nx,ny,nz
      integer imaxtov

      integer i,j,k,CCTK_Equals
      CCTK_REAL rr,xx,yy,zz,mval,phival,mgrr,dgrr,r_schwarval
      CCTK_REAL alpha,dalpha,pi
      CCTK_REAL one,zero,tov_ns_atmos_fact

      CCTK_REAL rho_temp,m_temp,phi_temp,eps_temp
      logical tov_brill
      CCTK_REAL tov_brill_lambda_r,tov_brill_lambda_z,tov_brill_amp
      CCTK_REAL tov_brill_r_cent
      CCTK_REAL cos_phi,sin_phi,g_rho_rho,g_rho_z,g_z_z
      CCTK_REAL r2_brill,q_brill
#ifdef CACTUSEOS_EOS_1D_BASE
#include "EOS_1D_Base.inc"
#endif

      nx = cctk_lsh(1)
      ny = cctk_lsh(2)
      nz = cctk_lsh(3)

      one = 1.0d0
      zero = 0.0d0
      tiny = 1.0d-29
      rmaxtov = mahc_init_rmaxtov
      ntov = mahc_init_ntov
      tov_centden = centden
      tov_eos_gamma = eos_gamma
      tov_eos_k = eos_k

c     clear out the memory, since T3E seems to like this
      do i=1,ntov
         presstov(i) = 0.0d0
         epstov(i) = 0.0d0
         rhotov(i) = 0.0d0
         phitov(i) = 0.0d0
         mtov(i) = 0.0d0
         r_schwartov(i) = 0.0d0
         r_isotov(i) = 0.0d0
      enddo

      tov_ns_atmos_fact = ns_atmos_fact

      pi = 4.0d0 * atan(1.0d0)

      imaxtov = 0
      call solvetov_iso(presstov,rhotov,epstov,mtov,phitov,r_schwartov,
     .   r_isotov,ntov,imaxtov,tov_centden,tov_eos_gamma,tov_eos_k,rmaxtov,
     .   mahc_init_eos_handle)

      write(*,*) '*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/'
      write(*,*) '*/*/              MAHC_INIT_DATA                */*/'
      write(*,*) '*/*/             TOV INITIAL DATA               */*/'
      write(*,*) '*/*/          (Isotropic coordinates)           */*/'
      write(*,*) '*/*/             ----------------               */*/'
      write(*,*) '*/*/         Evolving a Neutron Star has        */*/'
      write(*,*) '*/*/            never been so fun!              */*/'
      write(*,*) '*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/'
      write(*,*) '*/*/ tov_iso:  central density is ',tov_centden
      write(*,*) '*/*/           tov_eos_gamma is ',tov_eos_gamma
      write(*,*) '*/*/           tov_eos_k is ',tov_eos_k
      write(*,*) '*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/'

      do k=1,nz
      do j=1,ny
      do i=1,nx

         rr = r(i,j,k) + tiny
         xx = x(i,j,k)
         yy = y(i,j,k)
         zz = z(i,j,k)

         call gettov_iso(press(i,j,k),rho(i,j,k),eps(i,j,k),mval,phival,
     .      r_schwarval,rr,presstov,rhotov,epstov,mtov,phitov,r_schwartov,
     .      r_isotov,ntov,imaxtov,rmaxtov)

         alpha = exp(phival)

c        DO USUAL TOV
c        metric components (r = 0 taken care of by tiny and gettov_iso)
         gxx(i,j,k) = (r_schwarval/rr)**2
         gxy(i,j,k) = 0.0d0
         gxz(i,j,k) = 0.0d0
         gyy(i,j,k) = gxx(i,j,k)
         gyz(i,j,k) = 0.0d0
         gzz(i,j,k) = gxx(i,j,k)

c        extrinsic curvature
         kxx(i,j,k) = 0.0d0
         kxy(i,j,k) = 0.0d0
         kxz(i,j,k) = 0.0d0
         kyy(i,j,k) = 0.0d0
         kyz(i,j,k) = 0.0d0
         kzz(i,j,k) = 0.0d0

c        lapse
         alp(i,j,k) = alpha

c        atmosphere
         if(rho(i,j,k) .lt. tov_centden/tov_ns_atmos_fact) then
            rho(i,j,k) = tov_centden/tov_ns_atmos_fact
         endif

       

      enddo
      enddo
      enddo

      if(.not.(CCTK_EQUALS(mahc_init_eos,'tables'))) then
         eps = (tov_eos_k * rho**(tov_eos_gamma - 1.0d0)) / 
     .        (tov_eos_gamma - 1.0d0 + 1.e-30)
      else
#ifdef CACTUSEOS_EOS_1D_BASE
         do k=1,nz
         do j=1,ny
         do i=1,nx
            eps(i,j,k) = EOS_1D_epsfromrho(mahc_init_eos_handle,rho(i,j,k))
         enddo
         enddo
         enddo

#else
         write(*,*) 'if you want to use tables in mahc_init_data, you'
         write(*,*) 'need to compile in with EOS_1D_Base'
         call CCTK_WARN(0,"tov_iso")
#endif

      endif

c     set evolved quantities

      dens = rho * sqrt( gxx*gyy*gzz +
     .      2.0d0*gxy*gxz*gyz -
     .      gxx * (gyz**2) -
     .      gyy * (gxz**2) -
     .      gzz * (gxy**2) )
      tau = dens * eps
      sx = 0.0d0
      sy = 0.0d0
      sz = 0.0d0

      if(.not.(CCTK_EQUALS(mahc_eos,'tables'))) then
         press = (tov_eos_gamma - 1.0d0) * rho * eps
      else
#ifdef CACTUSEOS_EOS_1D_BASE
         do k=1,nz
         do j=1,ny
         do i=1,nx
            press(i,j,k) = EOS_1D_pressfromrho(mahc_init_eos_handle,rho(i,j,k))
         enddo
         enddo
         enddo

#else
         write(*,*) 'if you want to use tables in mahc_init_data, you'
         write(*,*) 'need to compile in with EOS_1D_Base'
         call CCTK_WARN(0,"tov_iso")
#endif
      endif

      return
      end