/*@@
   @file      SourceData.F77
   @date      
   @author    Gabrielle Allen
   @desc 
              Elliptic initial data for wave equation
   @enddesc 
 @@*/

#include "cctk.h"
#include "cctk_arguments.h"
#include "cctk_parameters.h"

      subroutine UniformCharge(CCTK_FARGUMENTS)

c     Find static field for a uniformly charge sphere
c
c     That is, solve Nabla^2 phi = - 4 pi rho
c     where rho = Q/(4/3 * Pi * R^3)
c     where Q is the total charge and R is the sphere radius

      implicit none

      DECLARE_CCTK_FARGUMENTS
      DECLARE_CCTK_PARAMETERS

      CCTK_REAL pi
      CCTK_REAL AbsTol(3), RelTol(3)
      CCTK_REAL charge_factor

      integer iphi,iMcoeff,iNcoeff
      integer i,j,k,ierr


c     Get variable indices for all grid functions

      call CCTK_VarIndex (iMcoeff,  "idscalarwaveelliptic::Mcoeff")
      if (iMcoeff .lt. 0) then
	call CCTK_WARN(0,"Grid variable index for Mcoeff not found")
      end if

      call CCTK_VarIndex (iNcoeff,  "idscalarwaveelliptic::Ncoeff")
      if (iNcoeff .lt. 0) then
	call CCTK_WARN(0,"Grid variable index for Ncoeff not found")
      end if

      call CCTK_VarIndex (iphi,"wavetoy::phi")
      if (iphi .lt. 0) then
	call CCTK_WARN(0,"Grid variable index for iphi not found")
      end if


c     Set up all coefficients and initial guess for solution

      pi = 4.0*atan(1.0)
      charge_factor = 4.0d0*pi*charge*3.0d0/(4.0d0*pi*radius**3)

      do k=1, cctk_lsh(3)
         do j=1, cctk_lsh(2)
            do i=1, cctk_lsh(1)

               phi(i,j,k) = 0.0d0
                     
               Mcoeff(i,j,k) = 0.0d0

               if (r(i,j,k) <= radius) then
                  Ncoeff(i,j,k) = charge_factor
               else
                  Ncoeff(i,j,k) = 0d0
               end if

            end do
         end do
      end do
      

c     Set tolerance for stopping elliptic solve

      AbsTol(1)=1.0d-5
      AbsTol(2)=1.0d-5
      AbsTol(3)=1.0d-5

      RelTol(1)=-1
      RelTol(2)=-1
      RelTol(3)=-1


c     Set any options needed for the elliptic solve

c      call Ell_SetStringKey(ierr, "yes", "EllLinFlat::Bnd::Robin")
c      write (*,*) "SourceData: SetString ",ierr
      call Ell_SetRealKey(ierr, 1.0d0, "EllLinFlat::Bnd::Robin::inf")
      call Ell_SetIntKey (ierr, 1,     "EllLinFlat::Bnd::Robin::falloff")


c     Call elliptic solver to fill out phi

      write (*,*) "charge: Going into elliptic solver"
      call Ell_LinFlatSolver(
     &	   ierr,
     &     cctkGH, 
     &     iphi, 
     &     iMcoeff, iNcoeff, 
     &     AbsTol, RelTol,
     &     "sor")
      write (*,*) "charge: Exit elliptic solver ierr =",ierr


c     Set up last timestep ... assume time symmetry

      do k=1, cctk_lsh(3)
         do j=1, cctk_lsh(2)
            do i=1, cctk_lsh(1)
                     
               phi_old(i,j,k) = phi(i,j,k)

            end do
         end do
      end do


c     Output exact solution if required

      if (output_tmp==1) then
        do k=1, cctk_lsh(3)
           do j=1, cctk_lsh(2)
              do i=1, cctk_lsh(1)
                     
                 if (r(i,j,k) >= radius) then
                    temp(i,j,k) = charge/r(i,j,k)
                 else
                    temp(i,j,k) = charge/(2.0d0*radius**3)*
     &                   (3.0d0*radius**2-r(i,j,k)**2)
                 end if

              end do
           end do
        end do

        call CCTK_OutputVarAsByMethod(ierr,cctkGH,
     &        "idscalarwaveelliptic::temp",
     &        "IOASCII_1D","phi_exact")
      end if      

      return
      end