the linear wave routines

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

1 files changed, 224 insertions, 0 deletions

diff --git a/src/planewaves.F b/src/planewaves.F
new file mode 100644
index 0000000..17b8431
--- /dev/null
+++ b/src/planewaves.F
@@ -0,0 +1,224 @@
+/*@@
+   @file      planewaves.f
+   @date      Wed Jan 24 16:51:03 1996
+   @author    Malcolm Tobias + Joan Masso
+   @desc 
+	Routine to initialize linear plane wave spacetime
+   @enddesc 
+   @hdate Tue Mar 16 12:38:04 1999 @hauthor Gerd Lanfermann
+   @hdesc Converted to Cactus4.0. the code and its comments of the
+   original authors are kept. There is NO BM suport at this point.
+   Id like to know how Cactus handles this w/o #ifdef !
+      
+@@*/
+
+ /*@@
+   @routine    planewaves
+   @date       Mon Jan 29 11:57:00 1996
+   @author     Malcolm Tobias + Joan Masso
+   @desc   
+	Routine to initialize a linear plane wave travelling in an 
+	arbitrary direction.  The form of the packet is:  <p>   
+	$ A*exp\left[-(kp_ix^i-\omega_p (time-ra))^2\right]
+        cos(k_ix^i-\omega \ time) $ 
+	<p>
+	where: <p>
+	A = amplitude of the wave  <br>
+	k  = the wave # of the sine wave <br>
+	$\omega$ = the frequency of the sine wave  <br>
+	kp = the wave # of the gaussian modulating the sine wave<br>
+	$\omega_p$ = the frequency of the gaussian <br>
+	ra = the initial position of the packet(s)
+   @enddesc 
+   @calledby linearWaves
+   @calls 
+@@*/
+
+#include "cctk.h"
+#include "declare_arguments.h"
+#include "declare_parameters.h"
+
+      subroutine planewaves(CCTK_FARGUMENTS)
+      implicit none
+      
+      DECLARE_CCTK_FARGUMENTS
+      DECLARE_PARAMETERS
+      
+      INTEGER iin,iout
+      REAL pi
+      REAL amplitude,ra,the,phi
+      REAL wave,wavep
+      REAL kx,ky,kz,w
+      REAL kxp,kyp,kzp,wp
+      
+c     local arrays (scalars! Man, what a sweat!)
+      REAL plus,minus,plusp,minusp,ain,aout
+      
+      INTEGERE i,j,k
+
+      pi = 3.14159265358979d0
+
+      c     Wave characteristics
+
+c     amplitude is declared by DECLARE_PARAMETER, for all other parameters
+c     we stick to the naming in the cactus3.2.0 version of cactus
+
+c     wavecentering
+      ra   = wavecenter
+     
+c     wavelength
+      wave = wavelength   
+
+c     wavepulse
+      wavep= wavepulse     
+
+c     determine whether the wave is ingoing, outgoing, or both
+      iin  = 0
+      iout = 0
+      if(CCTK_Equals(wavesgoing,'in').ne.0) then
+         iin = 1
+      elseif(CCTK_Equals(wavesgoing,'out').ne.0) then
+         iout = 1
+      elseif(CCTK_Equals(wavesgoing,'both').ne.0) then
+         iin = 1
+         iout = 1
+      endif
+
+c     determine the direction for the plane wave to travel
+c     and convert it from degrees to radians
+      the = pi*wavetheta/180.d0
+      phi = pi*wavephi/180.d0
+      
+      write(*,*)'Plane waves'
+      write(*,*)'amplitude  = ',amplitude
+      write(*,*)'wavecenter = ',ra
+      write(*,*)'wavelength = ',wave
+      write(*,*)'wavepulse  = ',wavep
+
+c     precalc
+      kx = 2*pi*sin(the)*cos(phi)/wave
+      ky = 2*pi*sin(the)*sin(phi)/wave
+      kz = 2*pi*cos(the)/wave
+      w = (kx*kx+ky*ky+kz*kz)**0.5
+      
+      kxp = 2*pi*sin(the)*cos(phi)/wavep
+      kyp = 2*pi*sin(the)*sin(phi)/wavep
+      kzp = 2*pi*cos(the)/wavep
+      wp = (kxp*kxp+kyp*kyp+kzp*kzp)**0.5
+
+      c  *************** plane waves ********************
+      do k=1,sh(3)
+         do j=1,sh(2)
+            do i=1,sh(1)		 
+               plus = (kx*x(i,j,k)+ky*y(i,j,k)+kz*z(i,j,k)+w*time)
+               minus = (kx*x(i,j,k)+ky*y(i,j,k)+kz*z(i,j,k)-w*time)
+               
+               plusp =(kxp*x(i,j,k)+kyp*y(i,j,k)+kzp*z(i,j,k)+wp*(time-ra))
+               minusp =(kxp*x(i,j,k)+kyp*y(i,j,k)+kzp*z(i,j,k)-wp*(time-ra))
+               
+               ain = iin*amplitude*cos(plus)*exp(-(plusp)**2)
+               aout = iout*amplitude*cos(minus)*exp(-(minusp)**2)
+	
+c     the metric functions
+        gxx(i,j,k) = cos(the)**2*cos(phi)**2*(1+ain+aout) +
+     &     sin(phi)**2*(1-ain-aout) + sin(the)**2*cos(phi)**2
+
+	gxy(i,j,k) = cos(the)**2*sin(phi)*cos(phi)*(1+ain+aout) +
+     &    sin(phi)*cos(phi)*(1-ain-aout) - sin(the)**2*sin(phi)*cos(phi)
+
+	gxz(i,j,k) = sin(the)*cos(the)*cos(phi)*(1+ain+aout) -
+     &    sin(the)*cos(the)*cos(phi)
+
+	gyy(i,j,k) = cos(the)**2*sin(phi)**2*(1+ain+aout) +
+     &    cos(phi)**2*(1-ain-aout) + sin(the)**2*sin(phi)**2
+
+	gyz(i,j,k) = -sin(the)*cos(the)*sin(phi)*(1+ain+aout) +
+     &    sin(the)*cos(the)*sin(phi)
+
+	gzz(i,j,k) = sin(the)**2*(1+ain+aout) + cos(the)**2
+
+
+c     At this pint, the CAC3.2.0 version has some code for the dxgxx,....
+c     defined, how does CAC4.0 handles this without  #ifdef THORN_BONAMASSO ?
+        
+        c  and finally the extrinsic curvatures
+c       Joan: I removed the division by alp as the initialization
+c       of the lapse may take place later. This data is consistent
+c       for initial lapse equal one.
+	
+	hxx(i,j,k) = amplitude*
+     &  (cos(the)**2*cos(phi)**2*
+     &  (iin*(-w*sin(plus)*exp(-(plusp)**2) -
+     &  2.d0*wp*plusp*cos(plus)*exp(-(plusp)**2)) + 
+     &  iout*(w*sin(minus)*exp(-(minusp)**2) +
+     &  2.d0*wp*minusp*cos(minus)*exp(-(minusp)**2))) - 
+     &  sin(phi)**2*
+     &  (iin*(-w*sin(plus)*exp(-(plusp)**2) -
+     &  2.d0*wp*plusp*cos(plus)*exp(-(plusp)**2)) + 
+     &  iout*(w*sin(minus)*exp(-(minusp)**2) +
+     &  2.d0*wp*minusp*cos(minus)*exp(-(minusp)**2))))/(-2.d0)
+
+	
+	hxy(i,j,k) = amplitude*
+     &  (cos(the)**2*sin(phi)*cos(phi)*
+     &  (iin*(-w*sin(plus)*exp(-(plusp)**2) -
+     &  2.d0*wp*plusp*cos(plus)*exp(-(plusp)**2)) + 
+     &  iout*(w*sin(minus)*exp(-(minusp)**2) +
+     &  2.d0*wp*minusp*cos(minus)*exp(-(minusp)**2))) -
+     &  sin(phi)*cos(phi)*
+     &  (iin*(-w*sin(plus)*exp(-(plusp)**2) -
+     &  2.d0*wp*plusp*cos(plus)*exp(-(plusp)**2)) + 
+     &  iout*(w*sin(minus)*exp(-(minusp)**2) +
+     &  2.d0*wp*minusp*cos(minus)*exp(-(minusp)**2))))/(-2.d0)
+
+	
+	hxz(i,j,k) = amplitude*
+     &  (sin(the)*cos(the)*cos(phi)*
+     &  (iin*(-w*sin(plus)*exp(-(plusp)**2) -
+     &  2.d0*wp*plusp*cos(plus)*exp(-(plusp)**2)) + 
+     &  iout*(w*sin(minus)*exp(-(minusp)**2) +
+     &  2.d0*wp*minusp*cos(minus)*exp(-(minusp)**2))))/(-2.d0)
+
+	
+	hyy(i,j,k) = amplitude*
+     &  (cos(the)**2*sin(phi)**2*
+     &  (iin*(-w*sin(plus)*exp(-(plusp)**2) -
+     &  2.d0*wp*plusp*cos(plus)*exp(-(plusp)**2)) + 
+     &  iout*(w*sin(minus)*exp(-(minusp)**2) +
+     &  2.d0*wp*minusp*cos(minus)*exp(-(minusp)**2))) -
+     &  cos(phi)**2*
+     &  (iin*(-w*sin(plus)*exp(-(plusp)**2) -
+     &  2.d0*wp*plusp*cos(plus)*exp(-(plusp)**2)) + 
+     &  iout*(w*sin(minus)*exp(-(minusp)**2) +
+     &  2.d0*wp*minusp*cos(minus)*exp(-(minusp)**2))))/(-2.d0)
+
+	
+	hyz(i,j,k) = amplitude*
+     &  (-sin(the)*cos(the)*sin(phi)*
+     &  (iin*(-w*sin(plus)*exp(-(plusp)**2) -
+     &  2.d0*wp*plusp*cos(plus)*exp(-(plusp)**2)) + 
+     &  iout*(w*sin(minus)*exp(-(minusp)**2) +
+     &  2.d0*wp*minusp*cos(minus)*exp(-(minusp)**2))))/(-2.d0)
+
+	 
+	hzz(i,j,k) = amplitude*
+     &  (sin(the)**2*
+     &  (iin*(-w*sin(plus)*exp(-(plusp)**2) -
+     &  2.d0*wp*plusp*cos(plus)*exp(-(plusp)**2)) + 
+     &  iout*(w*sin(minus)*exp(-(minusp)**2) +
+     &  2.d0*wp*minusp*cos(minus)*exp(-(minusp)**2))))/(-2.d0)
+	
+
+c     loop over sh ends here:
+          enddo 
+        enddo
+      enddo
+      
+      return
+      end
+      
+
+
+
+
+