/*@@
   @file      exa_integrate.F
   @date      October 1999
   @author    Miguel Alcubierre
   @desc 
      Integrate geodesics in exact spacetime.
   @enddesc 
@@*/

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

/*@@
  @routine    exa_integrate
  @date       August 1997
  @author     Miguel Alcubierre
  @desc
              This routine integrates geodesics in the exact
              spacetime using a fourth order Runge-Kutta. The geodesics
              are integrated using proper time as the parameter.

              It should be called after the main code has updated
              the geometry (including lapse and shift), and also
              after the routine "timegeodesic_num_integrate" has
              advanced the proper time.

  @enddesc 
  @calls
  @calledby
  @history

  @endhistory
@@*/

      subroutine exa_integrate(CCTK_ARGUMENTS)

      use arrays

      implicit none

      DECLARE_CCTK_PARAMETERS
      DECLARE_CCTK_ARGUMENTS

      integer i

      CCTK_REAL ts,xs,ys,zs,dtt,dtx,dty,dtz
      CCTK_REAL pxs,pys,pzs,dtpx,dtpy,dtpz
      CCTK_REAL kt1,kt2,kt3,kt4
      CCTK_REAL kx1,kx2,kx3,kx4
      CCTK_REAL ky1,ky2,ky3,ky4
      CCTK_REAL kz1,kz2,kz3,kz4
      CCTK_REAL kpx1,kpx2,kpx3,kpx4
      CCTK_REAL kpy1,kpy2,kpy3,kpy4
      CCTK_REAL kpz1,kpz2,kpz3,kpz4
      CCTK_REAL dtau
      CCTK_REAL half,third,sixth

!     Description of variables:
!
!     i,j,k            Counters.
!
!     ts               Auxilliary time.
!     xs,ys,zs         Auxilliary position.
!
!     dtx,dty,dtz      Derivatives of position with respect to
!                      proper time.
!
!     pxs,pys,pzs      Auxilliary momentum.
!
!     dtpx,dtpy,dtpz   Derivatives of momentum with respect to
!                      proper time.
!
!     kx#,ky#,kz#      Runge-Kutta contributions to the change of position.
!
!     kpx#,kpy#,kpz#   Runge-Kutta contributions to the change of momentum.
!
!     dtau             Increment in proper time.


!     *************************************
!     ***   IF  iteration = 0  RETURN   ***
!     *************************************

!     Because of the calling structure of CACTUS, this routine
!     is called even after time step 0.  This should not happen
!     because the geodesics will not be sincronized properly
!     any more.  To avoid it, I make the routine return here
!     if  iteration = 0.

      if (cctk_iteration.eq.0) return


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

      half = 0.5d0

      third = 1.0d0/3.0d0
      sixth = 1.0d0/6.0d0


!     *************************************
!     ***   LOOP OVER THE GRID POINTS   ***
!     *************************************

      do i=1,lgeos

!        **********************
!        ***   FIND  dtau   ***
!        **********************

         dtau = geodtau(i)


!        *******************************************
!        ***   FIRST RUNGE-KUTTA CONTRIBUTIONS   ***
!        *******************************************

         ts = geote(i)

         xs = geoxe(i)
         ys = geoye(i)
         zs = geoze(i)

         pxs = geopxe(i)
         pys = geopye(i)
         pzs = geopze(i)

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

         kt1 = dtau*dtt

         kx1 = dtau*dtx
         ky1 = dtau*dty
         kz1 = dtau*dtz

         kpx1 = dtau*dtpx
         kpy1 = dtau*dtpy
         kpz1 = dtau*dtpz


!        ********************************************
!        ***   SECOND RUNGE-KUTTA CONTRIBUTIONS   ***
!        ********************************************

         ts = geote(i) + half*kt1

         xs = geoxe(i) + half*kx1
         ys = geoye(i) + half*ky1
         zs = geoze(i) + half*kz1

         pxs = geopxe(i) + half*kpx1
         pys = geopye(i) + half*kpy1
         pzs = geopze(i) + half*kpz1

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

         kt2 = dtau*dtt

         kx2 = dtau*dtx
         ky2 = dtau*dty
         kz2 = dtau*dtz

         kpx2 = dtau*dtpx
         kpy2 = dtau*dtpy
         kpz2 = dtau*dtpz


!        *******************************************
!        ***   THIRD RUNGE-KUTTA CONTRIBUTIONS   ***
!        *******************************************

         ts = geote(i) + half*kt2

         xs = geoxe(i) + half*kx2
         ys = geoye(i) + half*ky2
         zs = geoze(i) + half*kz2

         pxs = geopxe(i) + half*kpx2
         pys = geopye(i) + half*kpy2
         pzs = geopze(i) + half*kpz2

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

         kt3 = dtau*dtt

         kx3 = dtau*dtx
         ky3 = dtau*dty
         kz3 = dtau*dtz

         kpx3 = dtau*dtpx
         kpy3 = dtau*dtpy
         kpz3 = dtau*dtpz


!        ********************************************
!        ***   FOURTH RUNGE-KUTTA CONTRIBUTIONS   ***
!        ********************************************

         ts = geote(i) + kt3

         xs = geoxe(i) + kx3
         ys = geoye(i) + ky3
         zs = geoze(i) + kz3

         pxs = geopxe(i) + kpx3
         pys = geopye(i) + kpy3
         pzs = geopze(i) + kpz3

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

         kt4 = dtau*dtt

         kx4 = dtau*dtx
         ky4 = dtau*dty
         kz4 = dtau*dtz

         kpx4 = dtau*dtpx
         kpy4 = dtau*dtpy
         kpz4 = dtau*dtpz


!        **************************************************
!        ***   UPDATE THE TIME, POSITION AND MOMENTUM   ***
!        **************************************************

!        Time.

         geote(i) = geote(i) + sixth*(kt1 + kt4)
     &            + third*(kt2 + kt3)

!        Position.

         geoxe(i) = geoxe(i) + sixth*(kx1 + kx4)
     &            + third*(kx2 + kx3)
         geoye(i) = geoye(i) + sixth*(ky1 + ky4)
     &            + third*(ky2 + ky3)
         geoze(i) = geoze(i) + sixth*(kz1 + kz4)
     &            + third*(kz2 + kz3)

!        Momentum.

         geopxe(i) = geopxe(i) + sixth*(kpx1 + kpx4)
     &             + third*(kpx2 + kpx3)
         geopye(i) = geopye(i) + sixth*(kpy1 + kpy4)
     &             + third*(kpy2 + kpy3)
         geopze(i) = geopze(i) + sixth*(kpz1 + kpz4)
     &             + third*(kpz2 + kpz3)

      end do


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

      if (findexa_IJ.eq.1) then
         call exa_IJ(CCTK_ARGUMENTS)
      end if


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

      return
      end