/*@@
   @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"

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

subroutine GRHydro_RefinementLevel(CCTK_ARGUMENTS)

  implicit none
  DECLARE_CCTK_ARGUMENTS

  GRHydro_reflevel = aint(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.
   @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
  CCTK_REAL :: psi4pt

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

  do k = 1, nz 
    do j = 1, ny 
      do i = 1, nx
        
        if (conformal_state .eq. 0) then
          call SpatialDeterminant(&
               gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),&
               gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),&
               GRHydro_det(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))        
        else
          psi4pt = psi(i,j,k)**4
          call SpatialDeterminant(&
               psi4pt*gxx(i,j,k),psi4pt*gxy(i,j,k),psi4pt*gxz(i,j,k),&
               psi4pt*gyy(i,j,k),psi4pt*gyz(i,j,k),psi4pt*gzz(i,j,k),&
               GRHydro_det(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),&
               psi4pt*gxx(i,j,k),psi4pt*gxy(i,j,k),psi4pt*gxz(i,j,k),&
               psi4pt*gyy(i,j,k),psi4pt*gyz(i,j,k),psi4pt*gzz(i,j,k))
        end if

      end do
    end do
  end do
  
end subroutine SetMetricTemps