/*@@
   @file      Utils.F
   @date      Sat Jan 26 02:28:46 2002
   @author    
   @desc 
   Utility functions for other thorns. Calculation of the determinant
   of the spatial metric and the upper metric.
   @enddesc 
 @@*/

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

 /*@@
   @routine    GRHydro_RefinementLevel
   @date       July 2005
   @author     
   @desc 
   Calculates the current refinement level from flesh data
   @enddesc 
   @calls     
   @calledby   
   @history 
 
   @endhistory 
@@*/

subroutine GRHydro_Debug(CCTK_ARGUMENTS)

  implicit none
  DECLARE_CCTK_ARGUMENTS
  integer i,j,k
  integer nx, ny, nz

  GRHydro_reflevel = int(log10(dble(cctk_levfac(1)))/log10(2.0d0))
  
  nx = cctk_lsh(1)
  ny = cctk_lsh(2)
  nz = cctk_lsh(3)

  if (GRHydro_reflevel .lt. 3) return
 
  do i=4,nx
     do j=4,ny
        do k=4,4
           if( r(i,j,k)-4.0d0 .lt. 96.0d0 .and. &
                r(i,j,k)+4.0d0 .gt. 96.0d0) then
              write(6,"(4i4,1P10E15.6)") GRHydro_reflevel, i,j,k,&
                   r(i,j,k),rho(i,j,k),eps(i,j,k),y_e(i,j,k),temperature(i,j,k)
           endif
        enddo
     enddo
  enddo

!  call flush(6)
!  if(GRHydro_reflevel.eq.4.and.temperature(1,1,1).lt.0.1d0) then
!     call CCTK_WARN(0,"stop")
!  endif

end subroutine GRHydro_Debug


subroutine GRHydro_RefinementLevel(CCTK_ARGUMENTS)

  implicit none
  DECLARE_CCTK_ARGUMENTS

  GRHydro_reflevel = int(log10(dble(cctk_levfac(1)))/log10(2.0d0))

end subroutine GRHydro_RefinementLevel


 /*@@
   @routine    SpatialDeterminant
   @date       Sat Jan 26 02:30:23 2002
   @author     
   @desc 
   Calculates the determinant of the spatial metric.
   PLEASE USE THE MACRO SPATIAL_DETERMINANT from now on!!!
   @enddesc 
   @calls     
   @calledby   
   @history 
 
   @endhistory 
@@*/

subroutine SpatialDeterminant(gxx,gxy,gxz,gyy,gyz,gzz,det)

  implicit none
  
  CCTK_REAL :: det
  CCTK_REAL :: gxx,gxy,gxz,gyy,gyz,gzz

  det = -gxz**2*gyy + 2*gxy*gxz*gyz - gxx*gyz**2 - gxy**2*gzz + gxx*gyy*gzz


!!$  Why is this weird order necessary? Search me. It just seemed 
!!$  to make a really odd NaN go away.

!!$  det =  -gxz**2*gyy + 2*gxy*gxz*gyz 
!!$  det = det - gxx*gyz**2 - gxy**2*gzz
!!$  det = det + gxx*gyy*gzz
  
  return
  
end subroutine SpatialDeterminant



/*@@
@routine    UpperMetric
@date       Sat Jan 26 02:32:26 2002
@author     
@desc 
Calculates the upper metric. The determinant is given, not
calculated.
@enddesc 
@calls     
@calledby   
@history 

@endhistory 
@@*/

subroutine UpperMetric(uxx, uxy, uxz, uyy, uyz, uzz, &
     det, gxx, gxy, gxz, gyy, gyz, gzz)
  
  implicit none
  
  CCTK_REAL :: uxx, uxy, uxz, uyy, uyz, uzz, det, &
       gxx, gxy, gxz, gyy, gyz, gzz, invdet
  
  invdet = 1.d0 / det
  uxx=(-gyz**2 + gyy*gzz)*invdet
  uxy=(gxz*gyz - gxy*gzz)*invdet
  uyy=(-gxz**2 + gxx*gzz)*invdet
  uxz=(-gxz*gyy + gxy*gyz)*invdet
  uyz=(gxy*gxz - gxx*gyz)*invdet
  uzz=(-gxy**2 + gxx*gyy)*invdet
  
  return
  
end subroutine UpperMetric

 /*@@
   @routine    SetMetricTemps
   @date       Wed Mar 23 10:01:08 2005
   @author     Ian Hawke
   @desc 
   If using the new EOS store the metric determinant and inverse metric
   to avoid repeated computation
   @enddesc 
   @calls     
   @calledby   
   @history 
 
   @endhistory 

@@*/

subroutine SetMetricTemps(CCTK_ARGUMENTS)

  implicit none

  DECLARE_CCTK_ARGUMENTS
  DECLARE_CCTK_PARAMETERS

  CCTK_INT :: i, j, k, nx, ny, nz

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

  do k = 1, nz 
    do j = 1, ny 
      do i = 1, nx
        
         GRHydro_det(i,j,k) = 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))
         call UpperMetric(&
              GRHydro_uxx(i,j,k),GRHydro_uxy(i,j,k),GRHydro_uxz(i,j,k),&
              GRHydro_uyy(i,j,k),GRHydro_uyz(i,j,k),GRHydro_uzz(i,j,k),&
              GRHydro_det(i,j,k),&
              gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),&
              gyy(i,j,k),gyz(i,j,k),gzz(i,j,k))        

      end do
    end do
  end do
  
end subroutine SetMetricTemps