1 files changed, 30 insertions, 18 deletions

diff --git a/src/SourceData.F b/src/SourceData.F
index 8bd4222..43bfe34 100644
--- a/src/SourceData.F
+++ b/src/SourceData.F
@@ -1,3 +1,11 @@
+ /*@@
+   @file      SourceData.F77
+   @date      
+   @author    Gabrielle Allen
+   @desc 
+              Elliptic initial data for wave equation
+   @enddesc 
+ @@*/
 
 #include "cctk.h"
 #include "cctk_arguments.h"
@@ -26,7 +34,6 @@ c     where Q is the total charge and R is the sphere radius
       pi = 4.0*atan(1.0)
 
       charge_factor = 4.0d0*pi*charge*3.0d0/(4.0d0*pi*radius**3)
-                  print *,charge_factor
 
       do k=1, cctk_lsh(3)
          do j=1, cctk_lsh(2)
@@ -34,8 +41,10 @@ c     where Q is the total charge and R is the sphere radius
                      
                Mcoeff(i,j,k) = 0.0d0
 
+               print *,r(i,j,k),radius
                if (r(i,j,k) <= radius) then
                   Ncoeff(i,j,k) = charge_factor
+                  print *,"Here"
                else
                   Ncoeff(i,j,k) = 0d0
                end if
@@ -58,14 +67,16 @@ c     where Q is the total charge and R is the sphere radius
       RelTol(2)=-1
       RelTol(3)=-1
 
-c     Call to elliptic solver will go here
+c     Call elliptic solver to fill out phi
+
       call Ell_LinFlatSolver(ierr,cctkGH, 
      &     iphi, 
      &     iMcoeff, iNcoeff, 
      &     AbsTol, RelTol,
      &     "sor")
 
-c     Set up last timestep
+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)
@@ -78,23 +89,24 @@ c     Set up last timestep
 
 
 c     Output exact solution if required
-      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
 
       if (output_tmp==1) then
-         call CCTK_OutputVarAsByMethod(ierr,cctkGH,
+        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