c     -*-Fortran-*-
c     $Header$

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



#define NVAR 1			/* number of equations to solve */
#define NAUX 4			/* number of auxiliary variables */
#define MAXLEVEL 3		/* number of existing multigrid levels */

#define VARIABLES				\
     VAR(0, phi,      "wavetoy::phi")		\
     RES(0, residual, "IDSWTEF::residual")	\
     AUX(0, x,        "grid::x")		\
     AUX(1, y,        "grid::y")		\
     AUX(2, z,        "grid::z")		\
     AUX(3, r,        "grid::r")



      subroutine IDSWTEF_uniform_charge (CCTK_ARGUMENTS)
      implicit none
      DECLARE_CCTK_ARGUMENTS
      DECLARE_CCTK_FUNCTIONS
      DECLARE_CCTK_PARAMETERS
      
      integer sym(3)
      integer i,j,k
      integer ierr
      
c     Set symmetries
      sym(1) = 1
      sym(2) = 1
      sym(3) = 1
      call SetCartSymVN (ierr, cctkGH, sym, "IDSWTEF::l1_phi")
      if (ierr.ne.0) call CCTK_WARN (0, "internal error")
      call SetCartSymVN (ierr, cctkGH, sym, "IDSWTEF::l2_phi")
      if (ierr.ne.0) call CCTK_WARN (0, "internal error")
      call SetCartSymVN (ierr, cctkGH, sym, "IDSWTEF::l3_phi")
      if (ierr.ne.0) call CCTK_WARN (0, "internal error")
      call SetCartSymVN (ierr, cctkGH, sym, "IDSWTEF::residual")
      if (ierr.ne.0) call CCTK_WARN (0, "internal error")
      call SetCartSymVN (ierr, cctkGH, sym, "IDSWTEF::l1_residual")
      if (ierr.ne.0) call CCTK_WARN (0, "internal error")
      call SetCartSymVN (ierr, cctkGH, sym, "IDSWTEF::l2_residual")
      if (ierr.ne.0) call CCTK_WARN (0, "internal error")
      call SetCartSymVN (ierr, cctkGH, sym, "IDSWTEF::l3_residual")
      if (ierr.ne.0) call CCTK_WARN (0, "internal error")
      
c     Initial data for the solver
      do k=1,cctk_lsh(3)
         do j=1,cctk_lsh(2)
            do i=1,cctk_lsh(1)
               phi(i,j,k) = 0
            end do
         end do
      end do
      
c     Call solver
      call IDSWTEF_call_solver (cctkgh)
      
c     Create time-symmetric initial data
      do k=1,cctk_lsh(3)
         do j=1,cctk_lsh(2)
            do i=1,cctk_lsh(1)
               phi_p(i,j,k) = phi(i,j,k)
            end do
         end do
      end do
      
      end
      
      
      
c     Calculate the residual
#include "MG_FARGUMENTS_pre.h"
      subroutine IDSWTEF_calc_residual (MG_FARGUMENTS)
#include "MG_FARGUMENTS_post.h"
      implicit none
#include "DECLARE_MG_FARGUMENTS.h"
      DECLARE_CCTK_FUNCTIONS
      DECLARE_CCTK_PARAMETERS
      
      CCTK_REAL pi
      CCTK_REAL rho		/* charge density */
      CCTK_REAL dx(3)
      CCTK_REAL idx2(3)		/* inverse squared grid spacing */
      
      integer i,j,k
      
      pi = 4 * atan(1.0d0)
      rho = charge / (4.0d0/3.0d0 * pi * radius**3)
      
      dx(1) = CCTK_DELTA_SPACE(1)
      dx(2) = CCTK_DELTA_SPACE(2)
      dx(3) = CCTK_DELTA_SPACE(3)
      
      idx2(1) = 1 / dx(1)**2
      idx2(2) = 1 / dx(2)**2
      idx2(3) = 1 / dx(3)**2
      
c     Calculate residual:
c     res = Laplace phi + 4 pi rho
      do k=1+cctk_nghostzones(3),cctk_lsh(3)-cctk_nghostzones(3)
         do j=1+cctk_nghostzones(2),cctk_lsh(2)-cctk_nghostzones(2)
            do i=1+cctk_nghostzones(1),cctk_lsh(1)-cctk_nghostzones(1)
               residual(i,j,k) =
     $              + (phi(i-1,j,k) - 2*phi(i,j,k) + phi(i+1,j,k)) * idx2(1)
     $              + (phi(i,j-1,k) - 2*phi(i,j,k) + phi(i,j+1,k)) * idx2(2)
     $              + (phi(i,j,k-1) - 2*phi(i,j,k) + phi(i,j,k+1)) * idx2(3)
               if (r(i,j,k).lt.radius) then
                  residual(i,j,k) = residual(i,j,k) + 4 * pi * rho
               end if
            end do
         end do
      end do
      
#include "RELEASE_MG_FARGUMENTS.h"
      end
      
      
      
c     Apply the boundary conditions to the solution
#include "MG_FARGUMENTS_pre.h"
      subroutine IDSWTEF_apply_bounds (MG_FARGUMENTS)
#include "MG_FARGUMENTS_post.h"
      implicit none
#include "DECLARE_MG_FARGUMENTS.h"
      DECLARE_CCTK_FUNCTIONS
      DECLARE_CCTK_PARAMETERS
      
      CCTK_REAL zero, one
      parameter (zero=0, one=1)
      
      integer ierr
      
c     Synchronise
      call CCTK_SyncGroupWithVarI (ierr, cctkGH, phi_index());
      if (ierr.ne.0) call CCTK_WARN (0, "internal error")
      
c     Apply the symmetry boundary conditions
c$$$      call CartSymGN (ierr, cctkGH, "WaveToy::scalarevolve")
c$$$      if (ierr.ne.0) call CCTK_WARN (0, "internal error")
      call MG_CartSymVI (ierr, cctkGH, phi_index());
      if (ierr.ne.0) call CCTK_WARN (0, "internal error")
      
c     Apply the outer boundary conditions
c$$$      call BndRobinGN (ierr, cctkGH, cctk_nghostzones, zero, 1,
c$$$     $     "WaveToy::scalarevolve")
c$$$      if (ierr.ne.0) call CCTK_WARN (0, "internal error")
      call MG_BndRobinVI (ierr, cctkGH, zero, one, 
     $     phi_index(), x_index(), y_index(), z_index(), r_index())
      if (ierr.ne.0) call CCTK_WARN (0, "internal error")
      
#include "RELEASE_MG_FARGUMENTS.h"
      end
      
      
      
c     Apply the boundary conditions to the auxiliary variables
#include "MG_FARGUMENTS_pre.h"
      subroutine IDSWTEF_fixup_aux (MG_FARGUMENTS)
#include "MG_FARGUMENTS_post.h"
      implicit none
#include "DECLARE_MG_FARGUMENTS.h"
      DECLARE_CCTK_FUNCTIONS
      DECLARE_CCTK_PARAMETERS
      
c     There really should be a way for the user to set specific
c     restriction operators for the auxiliary variables.
      
      CCTK_REAL dx(3)
      integer i,j,k
      
      dx(1) = CCTK_DELTA_SPACE(1)
      dx(2) = CCTK_DELTA_SPACE(2)
      dx(3) = CCTK_DELTA_SPACE(3)
      
      do k=1,cctk_lsh(3)
         do j=1,cctk_lsh(2)
            do i=1,cctk_lsh(1)
               x(i,j,k) = cctk_origin_space(1) + (cctk_lbnd(1) + i - 1) * dx(1)
               y(i,j,k) = cctk_origin_space(2) + (cctk_lbnd(2) + j - 1) * dx(2)
               z(i,j,k) = cctk_origin_space(3) + (cctk_lbnd(3) + k - 1) * dx(3)
               r(i,j,k) = sqrt(x(i,j,k)**2 + y(i,j,k)**2 + z(i,j,k)**2)
            end do
         end do
      end do
      
#include "RELEASE_MG_FARGUMENTS.h"
      end