/*@@
   @file      exa_dt.F
   @date      October 1999
   @author    Miguel Alcubierre
   @desc 
      Calculates time derivatives of the position and momentum
      in terms of proper time.
   @enddesc 
@@*/

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

/*@@
  @routine    exa_dt
  @date       August 1997
  @author     Miguel Alcubierre
  @desc
              This subroutine calculates time derivatives of the
              position and momentum in terms of proper time.

  @enddesc 
  @calls
  @calledby
  @history

  @endhistory
@@*/

      subroutine exa_dt(CCTK_ARGUMENTS,ts,xs,ys,zs,pxs,pys,pzs,
     .           dtt,dtx,dty,dtz,dtpx,dtpy,dtpz)

      implicit none

      DECLARE_CCTK_PARAMETERS
      DECLARE_CCTK_ARGUMENTS

      integer l,m,p,q

      CCTK_REAL xs,ys,zs,ts,pxs,pys,pzs
      CCTK_REAL dtt,dtx,dty,dtz,dtpx,dtpy,dtpz

      CCTK_REAL A,P0,IDET
      CCTK_REAL one,two
      CCTK_REAL aux1,aux2

      CCTK_REAL, dimension (3)     :: B,DA,PS
      CCTK_REAL, dimension (3,3)   :: GL,GU,DB,KL
      CCTK_REAL, dimension (3,3,3) :: DG

!     Description of variables:
!
!     l,m,p,q         Counters.
!
!     ts              Coordinate time.
!
!     xs,ys,zs        Position.
!
!     pxs,pys,pzs     Covariant momentum.
!
!     dtt             Derivative of coordinate time with respect to
!                     proper time.
!
!     dtx,dty,dtz     Derivatives of position with respect to proper time.
!
!     dtpx,dtpy,dtpz  Derivatives of momentum with respect to proper time.
!
!     A               Lapse.
!
!     B(i)            Shift.
!
!     GU(i,j)         Spatial metric.
!
!     GL(i,j)         Inverse spatial metric.
!
!     KL(i,j)         Extrinsic curvature
!
!     IDET            1/det(g)
!
!     DA(i)           d a / a
!                      i
!                        j
!     DB(i,j)         d b / 2
!                      i
!
!     DG(i,j,k)       d g  / 2
!                      i kj
!
!     PS(i)           Covariant momentum.
!

!     **************************
!     ***   DEFINE NUMBERS   ***
!     **************************

      one  = 1.0d0
      two  = 2.0d0


!     **************************************************
!     ***   FIND EXACT VALUES OF GEOMETRIC OBJECTS   ***
!     **************************************************

!     Here I add conditionals for the preprocessor.
!     If thorn_exact does not exist, the code should
!     still compile since it is always possible that
!     one does not want to integrate the geodesics on
!     the exact spacetime.  However, if the thorn is not
!     there and one still calls this routine, the code
!     should give a warning and stop cleanly.

#ifdef THORN_EXACT
      call exactdata(xs,ys,zs,ts,
     &   GL(1,1),GL(2,2),GL(3,3),GL(1,2),GL(2,3),GL(1,3),
     &   KL(1,1),KL(2,2),KL(3,3),KL(1,2),KL(2,3),KL(1,3),
     &   DG(1,1,1),DG(1,2,2),DG(1,3,3),DG(1,1,2),DG(1,2,3),DG(1,1,3),
     &   DG(2,1,1),DG(2,2,2),DG(2,3,3),DG(2,1,2),DG(2,2,3),DG(2,1,3),
     &   DG(3,1,1),DG(3,2,2),DG(3,3,3),DG(3,1,2),DG(3,2,3),DG(3,1,3),
     &   A,DA(1),DA(2),DA(3),B(1),B(2),B(3),
     &   DB(1,1),DB(1,2),DB(1,3),
     &   DB(2,1),DB(2,2),DB(2,3),
     &   DB(3,1),DB(3,2),DB(3,3))
#else
        GL = 0
        KL = 0
        DG = 0
        A = 0
        B = 0
        DA = 0
        DB = 0
#endif

      GL(2,1) = GL(1,2)
      GL(3,1) = GL(1,3)
      GL(3,2) = GL(2,3)

      KL(2,1) = KL(1,2)
      KL(3,1) = KL(1,3)
      KL(3,2) = KL(2,3)

      do l=1,3
         DG(l,2,1) = DG(l,1,2)
         DG(l,3,1) = DG(l,1,3)
         DG(l,3,2) = DG(l,2,3)
      end do


!     *******************************
!     ***   FIND INVERSE METRIC   ***
!     *******************************

!     Find the inverse of the determinant of the 3-metric.

      IDET = one/(GL(1,1)*GL(2,2)*GL(3,3) + two*GL(1,2)*GL(1,3)*GL(2,3)
     &     - (GL(1,1)*GL(2,3)**2 + GL(2,2)*GL(1,3)**2
     &     + GL(3,3)*GL(1,2)**2))

!     Find inverse metric components.

      GU(1,1) = IDET*(GL(2,2)*GL(3,3) - GL(2,3)**2)
      GU(2,2) = IDET*(GL(1,1)*GL(3,3) - GL(1,3)**2)
      GU(3,3) = IDET*(GL(1,1)*GL(2,2) - GL(1,2)**2)

      GU(1,2) = IDET*(GL(1,3)*GL(2,3) - GL(3,3)*GL(1,2))
      GU(1,3) = IDET*(GL(1,2)*GL(2,3) - GL(2,2)*GL(1,3))
      GU(2,3) = IDET*(GL(1,2)*GL(1,3) - GL(1,1)*GL(2,3))

      GU(2,1) = GU(1,2)
      GU(3,1) = GU(1,3)
      GU(3,2) = GU(2,3)


!     ********************
!     ***   FIND  P0   ***
!     ********************

      P0 = dsqrt(one + GU(1,1)*pxs**2 + GU(2,2)*pys**2 + GU(3,3)*pzs**2
     &   + two*(GU(1,2)*pxs*pys + GU(1,3)*pxs*pzs + GU(2,3)*pys*pzs))/A


!     ******************************
!     ***   DERIVATIVE OF TIME   ***
!     ******************************

      dtt = P0


!     ***********************************
!     ***   DERIVATIVES OF POSITION   ***
!     ***********************************

      dtx = - P0*B(1) + (GU(1,1)*pxs + GU(1,2)*pys + GU(1,3)*pzs)
      dty = - P0*B(2) + (GU(1,2)*pxs + GU(2,2)*pys + GU(2,3)*pzs)
      dtz = - P0*B(3) + (GU(1,3)*pxs + GU(2,3)*pys + GU(3,3)*pzs)


!     ***********************************
!     ***   DERIVATIVES OF POSITION   ***
!     ***********************************

!     Find contributions from derivatives of the lapse and shift
!     to the derivatives of the momentum.

      aux1 = - (A*P0)**2
      aux2 = two*P0

      dtpx = aux1*DA(1) + aux2*(pxs*DB(1,1) + pys*DB(1,2) + pzs*DB(1,3))
      dtpy = aux1*DA(2) + aux2*(pxs*DB(2,1) + pys*DB(2,2) + pzs*DB(2,3))
      dtpz = aux1*DA(3) + aux2*(pxs*DB(3,1) + pys*DB(3,2) + pzs*DB(3,3))

!     Add contributions from derivatives of the metric.

      PS(1) = pxs
      PS(2) = pys
      PS(3) = pzs

      do l=1,3
         do m=1,3
            do p=1,3
               do q=1,3

                  dtpx = dtpx + GU(l,p)*GU(m,q)*PS(p)*PS(q)*DG(1,l,m)
                  dtpy = dtpy + GU(l,p)*GU(m,q)*PS(p)*PS(q)*DG(2,l,m)
                  dtpz = dtpz + GU(l,p)*GU(m,q)*PS(p)*PS(q)*DG(3,l,m)

               end do
            end do
         end do
      end do


!     ***************
!     ***   END   ***
!     ***************

      return
      end