cosmetic changes + adding back the new profile for GaugeWave that was

lost in the last commit.


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/Exact/trunk@82 e296648e-0e4f-0410-bd07-d597d9acff87

1 files changed, 60 insertions, 64 deletions

diff --git a/src/GaugeWave.F b/src/GaugeWave.F
index b11e6c6..78d6ddd 100644
--- a/src/GaugeWave.F
+++ b/src/GaugeWave.F
@@ -4,7 +4,7 @@
 
 #define Pi (4 * atan(1.d0))
 
-      
+
       subroutine  gauge_wave(
      $     x, y, z, t,
      $     gdtt, gdtx, gdty, gdtz, 
@@ -25,12 +25,17 @@
 
       logical firstcall
 
-      CCTK_REAL a, o, H, half, d, fs
+      CCTK_REAL a, o, H, d, fs
+      CCTK_REAL zero,half,one
 
       data firstcall /.true./
-      save firstcall, a, o, half, d, fs
+      save firstcall, a, o, half, d
 
-      
+C     Numbers.
+
+      zero = 0.0d0
+      half = 0.5d0
+      one  = 1.0d0
 
 C     Get parameters of the exact solution.
 
@@ -38,25 +43,26 @@ C     Get parameters of the exact solution.
 
          a = GW_a
          o = GW_omega
-         half = 1/2d0
-C         d = GW_del / 8d0 / atan(1.d0)
-         d = GW_del * half / Pi
+         d = GW_del
          fs = GW_phase_shift
 
          firstcall = .false.
 
       end if
 
-C     How should the wave look like
+C     How should the wave look like.
+
       if (CCTK_Equals(GW_H,"sin").eq.1) then
-        
+
+        d = GW_del * half / Pi
+
         if (GW_diagonal.ne.0) then
 
-          H = 1 - a * sin((x-y)/sqrt(2.d0)/d - o*t/d)
+          H = one - a * sin((x-y)/sqrt(2.d0)/d - o*t/d)
  
         else
 
-          H = 1 - a * sin((x-o*t)/d - fs )
+          H = one - a * sin((x-o*t)/d - fs)
 
         end if
 
@@ -64,65 +70,55 @@ C     How should the wave look like
 
         H = exp(a*(x/d-o*t))
 
+      elseif (CCTK_Equals(GW_H,"gau").eq.1) then
+
+        H = one - a*dexp(-(x-t)**2/d**2)
+
       end if
 
-C     write metric      
-
-      if (GW_diagonal.ne.0) then
-         gdxx = half + half *H 
-         gdyy = half + half *H 
-         gdxy = half - half *H  
-         gdzz = 1
-         gdzx = 0
-         gdyz = 0
-
-         gdtt = - H
-         gdtx = 0
-         gdty = 0
-         gdtz = 0
-
-C     and upper metric
-
-         guxx = half + half / H 
-         guyy = half + half / H
-         guxy = half - half / H
-         guzz = 1
-         guyz = 0
-         guzx = 0
-      
-         gutt = - 1/H
-         gutx = 0
-         guty = 0
-         gutz = 0
+C     write metric.
+
+      if (GW_diagonal .eq. 1) then
+
+        gdxx = half * H + half
+        gdxy = -half* H + half
+        gdyy = half * H + half
+
+        guxx = half / H + half
+        guxy =-half / H + half
+        guyy = half / H + half
 
       else
 
-         gdxx = H 
-         gdyy = 1
-         gdzz = 1
-         gdxy = 0 
-         gdzx = 0
-         gdyz = 0
-
-         gdtt = - H
-         gdtx = 0
-         gdty = 0
-         gdtz = 0
-
-C     upper metric
-
-         guxx = 1/H
-         guyy = 1
-         guzz = 1
-         guxy = 0
-         guyz = 0
-         guzx = 0
-              
-         gutt = - 1/H
-         gutx = 0
-         guty = 0
-         gutz = 0
+        gdxx = H
+        gdxy = zero
+        gdyy = one
+
+        guxx = one/H
+        guxy = zero
+        guyy = one
+
       end if
+      
+      gdtt = - H
+      gdtx = zero
+      gdty = zero
+      gdtz = zero
+
+      gdzx = zero
+      gdyz = zero
+      gdzz = one
+
+C     and upper metric.
+
+      gutt = - one/H
+      gutx = zero
+      guty = zero
+      gutz = zero
+
+      guyz = zero
+      guzx = zero
+      guzz = one
 
       return
       end