c  Mahc_Evolve: Numerical Evolution of the General Relativistic Hydro Equations
c  Copyright (C) 2001 Mark Miller

/*@@
  @file      Mahc_Util.F
  @date      June 2001
  @author    Mark Miller
  @song      "Everybodys Got Something to Hide Except Me and My Monkey"
                                                                   - Beatles
  @desc
     Various Mahc Utilities.
  @enddesc
@@*/
#include "cctk.h"
#include "cctk_Parameters.h"
#include "cctk_Arguments.h"

/*@@
  @routine    mahc_calctmunu
  @date       June 2001
  @author     Mark Miller
  @song       "Sexy Sadie"   - Beatles
  @desc
     Calculate the (raised index) components of the stress energy tensor
     from the input 4-metric (3 metric, lapse, and shift) and the input
     perfect fluid components:
     INPUT:   3-metric components (indices lowered):  gxx,gxy,gxz,gyy,gyz,gzz
              shift (indices raised): betax,betay,betaz
              lapse: alp
              rest mass density scalar: rho
              specific internal energy scalar: eps
              pressure scalar: press
              velocity of fluid (indices raised): velx, vely, velz
      OUTPUT: stress energy tensor components (indices raised): ttt,ttx,
              tty,ttz,txx,txy,txz,tyy,tyz,tzz
  @enddesc
  @calls
  @calledby
  @history
  @endhistory
@@*/

      subroutine mahc_calctmunu(gxx,gxy,gxz,gyy,gyz,gzz,betax,betay,betaz,
     .   alp,rho,eps,press,velx,vely,velz,ttt,ttx,tty,ttz,txx,txy,txz,tyy,
     .   tyz,tzz)
      implicit none
      CCTK_REAL gxx,gxy,gxz,gyy,gyz,gzz,betax,betay,betaz
      CCTK_REAL alp,rho,eps,press,velx,vely,velz,ttt,ttx
      CCTK_REAL tty,ttz,txx,txy,txz,tyy,tyz,tzz

c     Local Vars
      CCTK_REAL det,uxx,uxy,uxz,uyy,uyz,uzz
      CCTK_REAL lorentz_fact,rhoenthalpy
      CCTK_REAL velxlow,velylow,velzlow,v2

c     calculate the inverse 3-metric:
      det = gxx*gyy*gzz + 2.0d0*gxy*gxz*gyz - gxx * (gyz**2) -
     .      gyy * (gxz**2) - gzz * (gxy**2)
      uxx=(-gyz**2 + gyy*gzz)/det
      uxy=(gxz*gyz - gxy*gzz)/det
      uyy=(-gxz**2 + gxx*gzz)/det
      uxz=(-gxz*gyy + gxy*gyz)/det
      uyz=(gxy*gxz - gxx*gyz)/det
      uzz=(-gxy**2 + gxx*gyy)/det

c     calculate the lorentz factor
      velxlow = gxx*velx + gxy*vely + gxz*velz
      velylow = gxy*velx + gyy*vely + gyz*velz
      velzlow = gxz*velx + gyz*vely + gzz*velz
      v2 = velx*velxlow + vely*velylow + velz*velzlow 
      lorentz_fact = 1.0d0/sqrt(1.0d0 - v2)

c     calculate the specific relativistic enthalpy * rho
      rhoenthalpy = rho + rho*eps + press

c     calculate the components of the stress energy tensor
      ttt = (rhoenthalpy*lorentz_fact**2 - press)/(alp**2)

      ttx = rhoenthalpy*(velx - betax/alp)*(lorentz_fact**2)/alp +
     .      press * betax / (alp**2)
      tty = rhoenthalpy*(vely - betay/alp)*(lorentz_fact**2)/alp +
     .      press * betay / (alp**2)
      ttz = rhoenthalpy*(velz - betaz/alp)*(lorentz_fact**2)/alp +
     .      press * betaz / (alp**2)

      txx = rhoenthalpy*(lorentz_fact**2)*
     .      (velx - betax/alp)*(velx - betax/alp) +
     .      press*(uxx - betax*betax/(alp**2))
      txy = rhoenthalpy*(lorentz_fact**2)*
     .      (velx - betax/alp)*(vely - betay/alp) +
     .      press*(uxy - betax*betay/(alp**2))
      txz = rhoenthalpy*(lorentz_fact**2)*
     .      (velx - betax/alp)*(velz - betaz/alp) +
     .      press*(uxz - betax*betaz/(alp**2))
      tyy = rhoenthalpy*(lorentz_fact**2)*
     .      (vely - betay/alp)*(vely - betay/alp) +
     .      press*(uyy - betay*betay/(alp**2))
      tyz = rhoenthalpy*(lorentz_fact**2)*
     .      (vely - betay/alp)*(velz - betaz/alp) +
     .      press*(uyz - betay*betaz/(alp**2))
      tzz = rhoenthalpy*(lorentz_fact**2)*
     .      (velz - betaz/alp)*(velz - betaz/alp) +
     .      press*(uzz - betaz*betaz/(alp**2))

      return
      end