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

      subroutine set_astar_trk(CCTK_FARGUMENTS)
      implicit none
      DECLARE_CCTK_FARGUMENTS
      DECLARE_CCTK_PARAMETERS
      INTEGER CCTK_Equals

c     local vars
      integer i,j,k
      CCTK_REAL uxx,uxy,uxz,uyy,uyz,uzz


c     get the trace out
      trK = kxx*ivp_uxx + kyy*ivp_uyy + kzz*ivp_uzz +
     .   2.0d0*kxy*ivp_uxy + 2.0d0*kxz*ivp_uxz + 2.0d0*kyz*ivp_uyz

c     take out the trace
      kxx = kxx - gxx * trK / 3.0d0
      kxy = kxy - gxy * trK / 3.0d0
      kxz = kxz - gxz * trK / 3.0d0
      kyy = kyy - gyy * trK / 3.0d0
      kyz = kyz - gyz * trK / 3.0d0
      kzz = kzz - gzz * trK / 3.0d0

c     Now, since York has astar with covarient indices, we need
c     to raise them.  And since Mathematica cannot handle "_"
c     in variable names, we need to do this as a 3d loop with
c     temp variables.  Sigh.

      do k=1,cctk_lsh(3)
      do j=1,cctk_lsh(2)
      do i=1,cctk_lsh(1)
         uxx = ivp_uxx(i,j,k)
         uxy = ivp_uxy(i,j,k)
         uxz = ivp_uxz(i,j,k)
         uyy = ivp_uyy(i,j,k)
         uyz = ivp_uyz(i,j,k)
         uzz = ivp_uzz(i,j,k)

         ivp_astarxx(i,j,k) = 
     &  uxx**2*kxx(i,j,k)+2.d0*uxx*uxy*kxy(i,j,k)+2.d0*uxx*uxz*kxz(i,j,k
     &  )+uxy**2*kyy(i,j,k)+2.d0*uxy*uxz*kyz(i,j,k)+uxz**2*kzz(i,j,k)

         ivp_astarxy(i,j,k) =
     &  uxx*uxy*kxx(i,j,k)+uxy**2*kxy(i,j,k)+uxx*uyy*kxy(i,j,k)+uxy*uxz*
     &  kxz(i,j,k)+uxx*uyz*kxz(i,j,k)+uxy*uyy*kyy(i,j,k)+uxz*uyy*kyz(i,j
     &  ,k)+uxy*uyz*kyz(i,j,k)+uxz*uyz*kzz(i,j,k)

         ivp_astarxz(i,j,k) =
     &  uxx*uxz*kxx(i,j,k)+uxy*uxz*kxy(i,j,k)+uxx*uyz*kxy(i,j,k)+uxz**2*
     &  kxz(i,j,k)+uxx*uzz*kxz(i,j,k)+uxy*uyz*kyy(i,j,k)+uxz*uyz*kyz(i,j
     &  ,k)+uxy*uzz*kyz(i,j,k)+uxz*uzz*kzz(i,j,k)

         ivp_astaryy(i,j,k) =
     &  uxy**2*kxx(i,j,k)+2.d0*uxy*uyy*kxy(i,j,k)+2.d0*uxy*uyz*kxz(i,j,k
     &  )+uyy**2*kyy(i,j,k)+2.d0*uyy*uyz*kyz(i,j,k)+uyz**2*kzz(i,j,k)

         ivp_astaryz(i,j,k) =
     &  uxy*uxz*kxx(i,j,k)+uxz*uyy*kxy(i,j,k)+uxy*uyz*kxy(i,j,k)+uxz*uyz
     &  *kxz(i,j,k)+uxy*uzz*kxz(i,j,k)+uyy*uyz*kyy(i,j,k)+uyz**2*kyz(i,j
     &  ,k)+uyy*uzz*kyz(i,j,k)+uyz*uzz*kzz(i,j,k)

         ivp_astarzz(i,j,k) =
     &  uxz**2*kxx(i,j,k)+2.d0*uxz*uyz*kxy(i,j,k)+2.d0*uxz*uzz*kxz(i,j,k
     &  )+uyz**2*kyy(i,j,k)+2.d0*uyz*uzz*kyz(i,j,k)+uzz**2*kzz(i,j,k)

      enddo
      enddo
      enddo

      if(CCTK_Equals(ivp_trK_conformal_astar,"input").eq.1) then
         if(trK0.eq.1) then
            trK = 0.0d0
         endif
      endif

      return
      end