Add a new type of plane waves: these are sinusoidal, as necessary for

the Mexico tests.


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/IDLinearWaves/trunk@89 5c0f84ea-6048-4d6e-bfa4-55cd5f2e0dd7

4 files changed, 137 insertions, 0 deletions

diff --git a/param.ccl b/param.ccl
index 5793b63..36c0959 100644
--- a/param.ccl
+++ b/param.ccl
@@ -7,6 +7,7 @@ EXTENDS KEYWORD initial_data
 "teukwaves"           :: "linear waves initial data- Teukolsky waves"
 "planewaves"          :: "linear waves initial data- plane waves"
 "standing_planewaves" :: "linear waves initial data- standing plane waves"
+"sine_planewaves"     :: "linear waves initial data- sine shaped plane waves"
 }
 
 USES KEYWORD metric_type 
diff --git a/schedule.ccl b/schedule.ccl
index e2b8975..db7c85e 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -5,6 +5,7 @@ if (CCTK_Equals(initial_data,"planewaves"))
   schedule IDLinearWaves_ParamChecker at CCTK_PARAMCHECK
   {
     LANG: C
+    OPTIONS: global
   } "Check that the metric_type is recognised"
 
   schedule IDLinearWaves_PlaneWaves in ADMBase_InitialData
@@ -18,6 +19,7 @@ if (CCTK_Equals(initial_data,"standing_planewaves"))
   schedule IDLinearWaves_ParamChecker at CCTK_PARAMCHECK
   {
     LANG: C
+    OPTIONS: global
   } "Check that the metric_type is recognised"
 
   schedule IDLinearWaves_StandWaves in ADMBase_InitialData
@@ -31,6 +33,7 @@ if (CCTK_Equals(initial_data,"teukwaves"))
   schedule IDLinearWaves_ParamChecker at CCTK_PARAMCHECK
   {
     LANG: C
+    OPTIONS: global
   } "Check that the metric_type is recognised"
 
   schedule IDLinearWaves_TeukWaves in ADMBase_InitialData
@@ -39,3 +42,16 @@ if (CCTK_Equals(initial_data,"teukwaves"))
   } "Construct linear Teukolsky wave initial data"  
 }
 
+if (CCTK_Equals(initial_data,"sine_planewaves")) 
+{
+  schedule IDLinearWaves_ParamChecker at CCTK_PARAMCHECK
+  {
+    LANG: C
+    OPTIONS: global
+  } "Check that the metric_type is recognised"
+
+  schedule IDLinearWaves_SinePlaneWaves in ADMBase_InitialData
+  {
+    LANG: Fortran
+  } "Construct linear plane wave initial data"  
+}
diff --git a/src/make.code.defn b/src/make.code.defn
index ff57d67..33166fc 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -4,6 +4,7 @@
 # Source files in this directory
 SRCS = 	ParamCheck.c\
 	planewaves.F77\
+	sine_planewaves.F77\
 	standing_planewaves.F77\
 	teukwaves.F77
 
diff --git a/src/sine_planewaves.F77 b/src/sine_planewaves.F77
new file mode 100644
index 0000000..689ae6d
--- /dev/null
+++ b/src/sine_planewaves.F77
@@ -0,0 +1,119 @@
+c     $Header$
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Functions.h"
+#include "cctk_Parameters.h"
+
+
+      subroutine IDLinearWaves_SinePlaneWaves (CCTK_ARGUMENTS)
+
+      implicit none
+      
+      DECLARE_CCTK_ARGUMENTS
+      DECLARE_CCTK_FUNCTIONS
+      DECLARE_CCTK_PARAMETERS
+      
+      CCTK_REAL pi
+      CCTK_REAL the,phi
+      CCTK_REAL st,ct,sp,cp
+      CCTK_REAL kx,ky,kz,w
+      CCTK_REAL b,bt
+      
+      INTEGER i,j,k
+      
+      CHARACTER*200 infoline
+      
+      pi = 3.14159265358979d0
+
+c     determine the direction for the plane wave to travel
+c     and convert it from degrees to radians
+      the = pi*wavetheta/180
+      phi = pi*wavephi/180
+
+c     precalc
+      st = sin(the)
+      ct = cos(the)
+      sp = sin(phi)
+      cp = cos(phi)
+      kx = 2*pi*st*cp/wavelength
+      ky = 2*pi*st*sp/wavelength
+      kz = 2*pi*ct/wavelength
+      w = sqrt(kx*kx+ky*ky+kz*kz)
+
+c     too lazy
+      if (wavetheta .ne. 90) then
+         call CCTK_WARN (0, "wavetheta.ne.90 is not implemented")
+      end if
+
+c  *************** plane waves ********************
+      do k=1,cctk_lsh(3)
+         do j=1,cctk_lsh(2)
+            do i=1,cctk_lsh(1)           
+               
+               b = amplitude * sin(kx*x(i,j,k)+ky*y(i,j,k)+kz*z(i,j,k)-w*cctk_time)
+               bt = - amplitude * w * cos(kx*x(i,j,k)+ky*y(i,j,k)+kz*z(i,j,k)-w*cctk_time)
+               
+c     the metric functions
+               gxx(i,j,k) = sp*sp * (1+b) + cp*cp * (1)
+               gxy(i,j,k) = cp*sp * (1-(1+b))
+               gyy(i,j,k) = cp*cp * (1+b) + sp*sp * (1)
+               gxz(i,j,k) = 0
+               gyz(i,j,k) = 0
+               gzz(i,j,k) = 1-b
+        
+c     and the extrinsic curvature
+               
+               kxx(i,j,k) = sp*sp * (bt/2) + cp*cp * (0)
+               kxy(i,j,k) = cp*sp * (-bt/2)
+               kyy(i,j,k) = cp*cp * (bt/2) + sp*sp * (0)
+               kxz(i,j,k) = 0
+               kyz(i,j,k) = 0
+               kzz(i,j,k) = bt/2
+               
+          enddo 
+        enddo
+      enddo
+      
+c     initialize the conformal factor
+      if (CCTK_EQUALS(metric_type, "static conformal")) then 
+         
+         conformal_state = 1
+         
+         if (CCTK_EQUALS(conformal_storage, "factor+derivs")) then
+            
+            conformal_state = 2
+            
+         else if (CCTK_EQUALS(conformal_storage, "factor+derivs+2nd derivs")) then
+            conformal_state = 3
+            
+         end if
+         
+         do k=1,cctk_lsh(3)
+            do j=1,cctk_lsh(2)
+               do i=1,cctk_lsh(1)                
+                  
+                  psi(i,j,k)   = 1d0
+                  
+                  if(conformal_state .gt. 1) then
+                     psix(i,j,k)  = 0d0
+                     psiy(i,j,k)  = 0d0
+                     psiz(i,j,k)  = 0d0    
+                  endif
+                  
+                  if(conformal_state .gt. 2) then
+                     psixy(i,j,k) = 0d0
+                     psixz(i,j,k) = 0d0
+                     psiyz(i,j,k) = 0d0  
+                     psixx(i,j,k) = 0d0
+                     psiyy(i,j,k) = 0d0
+                     psizz(i,j,k) = 0d0
+                  endif
+                  
+               end do
+            end do
+         end do
+         
+      end if
+      
+      end