Adding Gaussian solution for the point r=0 (thanks to Horst Beyer for this)

git-svn-id: http://svn.cactuscode.org/arrangements/CactusWave/IDScalarWaveC/trunk@20 071ba4cb-2270-47cd-a45a-3e542e115664

1 files changed, 11 insertions, 5 deletions

diff --git a/src/InitialData.c b/src/InitialData.c
index 1ecc970..c7ac1d3 100644
--- a/src/InitialData.c
+++ b/src/InitialData.c
@@ -98,11 +98,17 @@ void IDScalarWaveC_InitialData(CCTK_ARGUMENTS)
 
           phi[index] = amplitude*exp( - sqr( (R - radius) / sigma ) );
 
-          phi_p[index] = amplitude/2.0*(R-dt)/R*
-	    exp( - sqr( (R - radius - dt)/ sigma ) )
-	    + amplitude/2.0*(R+dt)/R*
-	    exp( - sqr( (R - radius + dt)/ sigma ) )
-	    ;
+ 	  if (R == 0.0) 
+	  {
+	    phi_p[index] = (1.0 - 2.0*dt*dt/sigma)*exp(-dt*dt/sigma);
+	  }
+	  else
+	  {
+	    phi_p[index] = amplitude/2.0*(R-dt)/R*
+	      exp( - sqr( (R - radius - dt)/ sigma ) )
+	      + amplitude/2.0*(R+dt)/R*
+	      exp( - sqr( (R - radius + dt)/ sigma ) );
+	  }
         } 
       }
     }