/*@@
   @file      num_IJ.F
   @date      October 1999
   @author    Miguel Alcubierre
   @desc
      Interpolate the curvature invariants I and J at the
      positions of the geodesics in numerical spacerime.
   @enddesc 
@@*/

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

/*@@
  @routine    num_IJ
  @date       November 1997
  @author     Miguel Alcubierre
  @desc
              Interpolate the curvature invariants I and J
              at the positions of the geodesics in the
              numerical spacetime.

              The call for this routine comes from PsiKadelia.
              The reason for this is that PsiKadelia calculates
              the invariants only for output, i.e. the 3D arrays
              containing them do not persist after exit from
              PsiKadelia.  The interpolation into the geodesic
              positions must therefore be done while still inside
              PsiKadelia.

              Another quirk of this is the fact that PsiKadelia is
              only ever called if one of the quatities it calculates
              requires output.  One must therefore ask for output of
              the invariants at least as slice fields at the same
              time as one needs output from the geodesics.

  @enddesc 
  @calls
  @calledby
  @history 

  @endhistory 
@@*/

      subroutine num_IJ(CCTK_ARGUMENTS)

      use arrays

      implicit none

      integer npoints

      DECLARE_CCTK_PARAMETERS
      DECLARE_CCTK_ARGUMENTS

      CCTK_REAL, dimension (lgeos) :: XS,YS,ZS
      CCTK_REAL, dimension (lgeos) :: G1,G2,G3,G4


!     ***********************************
!     ***   FIND INVARIANTS I and J   ***
!     ***********************************

!     Interpolate invariants.

      npoints = lgeos

      XS = geoxn
      YS = geoyn
      ZS = geozn

#ifdef THORN_PSIKADELIA
c      call getpoints(icurvre3d,XS,YS,ZS,G1,npoints)
c      call getpoints(icurvim3d,XS,YS,ZS,G2,npoints)
c      call getpoints(jcurvre3d,XS,YS,ZS,G3,npoints)
c      call getpoints(jcurvim3d,XS,YS,ZS,G4,npoints)
#endif

      geonumI_r = G1
      geonumI_i = G2
      geonumJ_r = G3
      geonumJ_i = G4


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

      return
      end