/*@@
   @file      GRHydro_Bondi.F90 
   @date      Wed Jan 13 13:00:49 EST 2010
   @author    Scott C. Noble
   @desc 
   Hydro initial data for the relativistic Bondi solution about 
   a single Schwarzschild black hole. 
   @enddesc 
 @@*/

/*
   Calculates the Bondi solution, or the spherically symmetric hydrostationary 
   solution to a fluid on a static fixed background spacetime. We assume that one can
   calculate a radius "r" from the grid and that with respect to this radial coordinate, 
   the solution satisfies 

   d (\rho u^r) / dr    = 0 

   Assumes that the equation of state is  P = K \rho^\Gamma   and K  is set by 
   the location of the sonic point.     


 -- Implicitly assumes that there is no spin in the geometry as there is no Bondi 
    solution for spinning black holes.  If a spin is specified, a spherically symmetric 
    is still assumed but the 4-velocity is set consistently with the spinning spacetime. 
*/

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

# define M_PI     3.14159265358979323846d0  /* pi */

!!$Newton-Raphson parameters:

#define velx(i,j,k) vel(i,j,k,1)
#define vely(i,j,k) vel(i,j,k,2)
#define velz(i,j,k) vel(i,j,k,3)
#define sx(i,j,k) scon(i,j,k,1)
#define sy(i,j,k) scon(i,j,k,2)
#define sz(i,j,k) scon(i,j,k,3)
#define Bvecx(i,j,k) Bvec(i,j,k,1)
#define Bvecy(i,j,k) Bvec(i,j,k,2)
#define Bvecz(i,j,k) Bvec(i,j,k,3)
#define Bconsx(i,j,k) Bcons(i,j,k,1)
#define Bconsy(i,j,k) Bcons(i,j,k,2)
#define Bconsz(i,j,k) Bcons(i,j,k,3)

subroutine GRHydro_BondiM_Iso(CCTK_ARGUMENTS)

  implicit none
  
  DECLARE_CCTK_ARGUMENTS
  DECLARE_CCTK_PARAMETERS
  DECLARE_CCTK_FUNCTIONS
  
  CCTK_INT :: i, j, k, nx, ny, nz, imin, jb,N_points
  PARAMETER (N_points=2000)
  CCTK_REAL :: M, Msq, Mdot, rs, gam, rmin_bondi, rmax_bondi, cs_sq,cs,vs_sq,vs,rhos,gtemp,hs, Kval, Qdot
  CCTK_REAL :: logrmin,dlogr,rhotmp,utmp,vtmp,rspher
  CCTK_REAL :: r_bondi(N_points), logr_bondi(N_points), rho_bondi(N_points), u_bondi(N_points), v_bondi(N_points)
  CCTK_REAL :: drhodr, det, rhocheck, rhocheck2, riso, rnew, rsch, ucheck
  CCTK_REAL :: uiso, uisocheck, vcheck, ucheck2, vcheck2, xhat,yhat, zhat, xp, yp, zp

  !!$set_bondi_parameters
  M = bondi_central_mass(1)
  Msq  = M*M
  Mdot = mdot_sonicpt_bondi
  rs = r_sonicpt_bondi
  gam = gl_gamma

  write(*,*) 'Bondi_pars:',M,mdot_sonicpt_bondi,r_sonicpt_bondi,gl_gamma

  rmin_bondi = M * bondi_rmin(1)
  rmax_bondi = M * bondi_rmax(1)

  cs_sq  =  M / ( 2.*rs - 3.*M ) 

  if( cs_sq > (gam - 1.)) then 
     cs_sq = gam - 1.
     rs = 0.5 * M * ( 3. + 1./cs_sq ) 
  endif
  
  cs     =  sqrt(cs_sq)
  vs_sq  =  M / ( 2. * rs )  
  vs     =  sqrt(vs_sq)
  rhos   =  Mdot / ( 4. * M_PI * vs * rs * rs ) 
  gtemp  =  gam - 1.
  hs     =  1. / ( 1. - cs_sq / (gam - 1.) )
  Kval      = hs * cs_sq * rhos**(-gtemp) / gam 
  Qdot   = hs * hs * ( 1. - 3. * vs_sq ) 
  
  logrmin = log10(rmin_bondi)
  dlogr = (log10(rmax_bondi) - logrmin)/(1.*(N_points-1))

  write(*,*)'More pars:',cs,vs,rhos,hs,Kval,Qdot,logrmin,dlogr

  rhotmp=1.0d30
  imin=1

  do i=1,N_points   
     logr_bondi(i) = logrmin + dlogr*(i-1)
     r_bondi(i) = 10.**(logr_bondi(i))
     utmp = abs(r_bondi(i) - r_sonicpt_bondi) 
     if (utmp < rhotmp) then 
        rhotmp = utmp 
        imin = i 
     endif
  enddo

  rhotmp = -1.  !!$ start with guess

  do i=imin,N_points
     rspher = r_bondi(i) 
     call find_bondi_solution( rspher, rhotmp, utmp, vtmp, rs, rhos, M, Mdot, Kval, gam, Qdot )
     if(rhotmp < initial_rho_abs_min) then
        rhotmp = initial_rho_abs_min
        utmp = Kval * rhotmp**gl_gamma  / (gl_gamma - 1.)
     endif
     rho_bondi(i) = rhotmp
     u_bondi(i)   = utmp
     v_bondi(i)   = vtmp
  end do
  
  rhotmp = -1.
  
  do i=imin-1,1,-1
     rspher = r_bondi(i)
     call find_bondi_solution( rspher, rhotmp, utmp, vtmp, rs, rhos, M, Mdot, Kval, gam, Qdot )
     if(rhotmp < initial_rho_abs_min) then
        rhotmp = initial_rho_abs_min
        utmp = K * rhotmp**gl_gamma  / (gl_gamma - 1.)
     endif
     rho_bondi(i) = rhotmp
     u_bondi(i)   = utmp
     v_bondi(i)   = vtmp
  enddo
  
  write(*,*)"i=1:",r_bondi(1),rho_bondi(1),u_bondi(1),v_bondi(1)
  write(*,*)"i=100:",r_bondi(100),rho_bondi(100),u_bondi(100),v_bondi(100)
  write(*,*)"i=1000:",r_bondi(1000),rho_bondi(1000),u_bondi(1000),v_bondi(1000)
  write(*,*)"i=1500:",r_bondi(1500),rho_bondi(1500),u_bondi(1500),v_bondi(1500)

!!$    // find the derivative near r=M
  rnew = 2.25 * M
  j = floor ( 0.5 +  (log10(rnew) - logrmin) / dlogr )  
  rhocheck = rho_bondi(j)
  call find_bondi_solution(rnew,rhocheck, ucheck, vcheck, rs, rhos, M, Mdot, Kval, gam, Qdot )
  uisocheck = 4.0*vcheck/3.0
  
  rnew = 0.25 * 3.02**2 * M/1.01
  j = floor( 0.5 +  (log10(rnew) - logrmin) / dlogr ) 
  rhocheck2 = rho_bondi(j)
  call find_bondi_solution( rnew, rhocheck2, ucheck2, vcheck2, rs, rhos, M, Mdot, Kval, gam, Qdot )
  drhodr = 100.0*(rhocheck2-rhocheck)/M 

  write(6,*)'Rhocheck:',rhocheck,rhocheck2,drhodr

  nx = cctk_lsh(1)
  ny = cctk_lsh(2)
  nz = cctk_lsh(3)
  
  do i=1,nx
     do j=1,ny
        do k=1,nz
           
           xp=x(i,j,k)
           yp=y(i,j,k)
           zp=z(i,j,k)
           
           riso = sqrt(xp*xp + yp*yp + zp*zp +1.0e-16)
           xhat = xp/riso
           yhat = yp/riso
           zhat = zp/riso
           
           if(riso < 1.0e-7) then
              gxx(i,j,k) = 1.0e4
              gyy(i,j,k) = 1.0e4
              gzz(i,j,k) = 1.0e4
              gxy(i,j,k) = 0.0
              gxz(i,j,k) = 0.0
              gyz(i,j,k) = 0.0
           endif

           alp(i,j,k) = 1.0/gxx(i,j,k)**2

           if(riso > M) then
              rsch = 0.25 * ( 2.*riso + M)**2 / riso
              jb = floor( 0.5 +  (log10(rsch) - logrmin) / dlogr ) 

              if(jb > N_points)jb = N_points

              rhotmp = rho_bondi(jb)
              call find_bondi_solution( rsch,rhotmp, utmp, vtmp, rs, rhos, M, Mdot, Kval, gam, Qdot) 
              rho(i,j,k) = rhotmp
              uiso = vtmp / (1.0 - M/2.0/riso) / (1.0+ M/2.0/riso)
           else
              if(riso > 0.5*M) then
                 rho(i,j,k) = rhocheck+drhodr*riso*(riso-M)/M
              else 
                 rho(i,j,k) = (rhocheck-drhodr*M/4.0)*(1.-cos(2.*M_PI*riso/M))/2.0
              endif
              utmp = Kval * rho(i,j,k)**( gam ) / (gam - 1.)
              uiso = uisocheck * riso / M
           endif
           eps(i,j,k) = utmp/rhotmp
           
           w_lorentz(i,j,k) = sqrt(1.0+gxx(i,j,k) * uiso**2)
           velx(i,j,k) = -1.0*uiso * xhat / w_lorentz(i,j,k)
           vely(i,j,k) = -1.0*uiso * yhat / w_lorentz(i,j,k)
           velz(i,j,k) = -1.0*uiso * zhat / w_lorentz(i,j,k)
           
           det=SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k))

           Bvecx(i,j,k) = bondi_bmag*M**2*xhat/sqrt(det)/riso**2
           Bvecy(i,j,k) = bondi_bmag*M**2*yhat/sqrt(det)/riso**2
           Bvecz(i,j,k) = bondi_bmag*M**2*zhat/sqrt(det)/riso**2
           
           call Prim2ConGenM(GRHydro_eos_handle,gxx(i,j,k),gxy(i,j,k), &
		  gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k), &
		  det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k), &
		  tau(i,j,k),Bconsx(i,j,k),Bconsy(i,j,k),Bconsz(i,j,k), &
                  rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k), &
		  eps(i,j,k),press(i,j,k), &
                  Bvecx(i,j,k),Bvecy(i,j,k),Bvecz(i,j,k),w_lorentz(i,j,k))
           
        end do
     end do
  end do

  densrhs = 0.d0
  srhs = 0.d0
  taurhs = 0.d0
  Bconsrhs = 0.d0

  return

end subroutine GRHydro_BondiM_Iso