/*@@
   @file      exa_init.F
   @date      October 1999
   @author    Miguel Alcubierre
   @desc 
      Initial position and covariant momentum of geodesics
      in exact spacetime.
   @enddesc 
@@*/

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

/*@@
  @routine    exa_init
  @date       August 1997
  @author     Miguel Alcubierre
  @desc

     Gives initial position and covariant momentum of geodesics
     in exact spacetime.

     The initial positions of the geodesic are chosen to correspond
     with the geodesics of the numerical spacetime.

     I also choose the initial geodesics to point in the normal
     direction to the hypersurfaces.

  @enddesc 
  @calls
  @calledby   INITIAL3
  @history 

  @endhistory 
@@*/

      subroutine exa_init(CCTK_ARGUMENTS)

      use arrays

      implicit none

      DECLARE_CCTK_PARAMETERS
      DECLARE_CCTK_ARGUMENTS

      integer i


!     *************************
!     ***   PRINT MESSAGE   ***
!     *************************

      call CCTK_INFO('Initializing timelike geodesics on exact spacetime')


!     **************************************
!     ***   ALLOCATE MEMORY FOR ARRAYS   ***
!     **************************************

      allocate(geote(1:lgeos))
      allocate(geoxe(1:lgeos),geoye(1:lgeos),geoze(1:lgeos))
      allocate(geopxe(1:lgeos),geopye(1:lgeos),geopze(1:lgeos))


!     ************************
!     ***   INITIAL TIME   ***
!     ************************

!     Get initial time (for the moment 0, but this should change).

      geote = 0.0d0


!     ******************************************
!     ***   INITIAL POSITIONS OF GEODESICS   ***
!     ******************************************

!     I assume here that the numerical and exact
!     spacetimes start from the same initial slice,
!     and that the initial positions of the geodesics
!     coincide.

      geoxe = geoxn
      geoye = geoyn
      geoze = geozn


!     *****************************************
!     ***   INITIAL MOMENTUM OF GEODESICS   ***
!     *****************************************

!     Notice that the variables I use here correspond to the
!     spatial components of the covariant momentum.  Since the
!     geodesics are moving initially in the normal direction to
!     the hypersurfaces, those components are zero.

      geopxe = 0.0d0
      geopye = 0.0d0
      geopze = 0.0d0


!     ****************************************
!     ***   CALCULATE INVARIANTS I and J   ***
!     ****************************************

      if (findexa_IJ.eq.1) then

         allocate(geoexaI_r(1:lgeos),geoexaI_i(1:lgeos))
         allocate(geoexaJ_r(1:lgeos),geoexaJ_i(1:lgeos))

         call exa_IJ(CCTK_ARGUMENTS)

      end if


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

      return
      end