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

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

#include "EllBase.h"

      subroutine UniformCharge(CCTK_ARGUMENTS)

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_ARGUMENTS
      DECLARE_CCTK_PARAMETERS
      DECLARE_CCTK_FUNCTIONS

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

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

      character*30 fsolver
      character*200 infoline

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.0d0*atan(1.0d0)
      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) .le. radius) then
                  Ncoeff(i,j,k) = charge_factor
               else
                  Ncoeff(i,j,k) = 0.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

      call Ell_SetStrKey (ierr, "yes", "EllLinFlat::Bnd::Robin")
      call Ell_SetRealKey(ierr, 0.0d0, "EllLinFlat::Bnd::Robin::inf")
      call Ell_SetIntKey (ierr, 1,     "EllLinFlat::Bnd::Robin::falloff")

c     Set parameters for specific solvers
      if (CCTK_EQUALS(solver,"sor")) then
        call Ell_SetIntKey(ierr, sor_maxit,"Ell::SORmaxit")
      end if
     
      call CCTK_FortranString(length,solver,fsolver)	
      write(infoline,'("Going into elliptic solver ",A)') fsolver
      call CCTK_INFO(infoline)

c     Call elliptic solver to fill out phi

      call Ell_LinFlatSolver(
     &	   ierr,
     &     cctkGH, 
     &     iphi, 
     &     iMcoeff, iNcoeff, 
     &     AbsTol, RelTol,
     &     fsolver)
	
      if (ierr .eq. ELL_SUCCESS) then
        call CCTK_INFO("Leaving elliptic solver: solve successful")
      else if (ierr .eq. ELL_NOCONVERGENCE) then
	call CCTK_INFO("Leaving elliptic solver: solver failed to converge")
      else if (ierr .eq. ELL_NOSOLVER) then
	call CCTK_INFO("Elliptic solver not found")
      else 
	write(infoline,'("Leaving elliptic solver: solve failed (Error ",I3,")")') ierr
        call CCTK_INFO(infoline)
      end if

c     Set up last timestep ... assume (first order) time symmetry

      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


c     Output exact solution if required

      if (output_tmp .eq. 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) .ge. 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