the linear wave routines

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

3 files changed, 852 insertions, 0 deletions

diff --git a/src/LinearWaves.F b/src/LinearWaves.F
new file mode 100644
index 0000000..c0dc44e
--- /dev/null
+++ b/src/LinearWaves.F
@@ -0,0 +1,59 @@
+c/*@@
+c  @file      linearWaves.F
+c  @date      
+c  @author    Joan Masso
+c  @desc 
+c             Driver for the analytical (or almost) linear wave models 
+c  @enddesc 
+c@@*/
+
+c Note that including cactus.h will also include linearWaves.h
+#include "cactus.h"
+
+c/*@@
+c  @routine    linearWaves
+c  @date       
+c  @author     Joan Masso 
+c  @desc 
+c  
+c  @enddesc 
+c  @calls     planewaves teukwaves
+c  @calledby   
+c  @history Converted to Cactus4.0. the code and its comments of the
+   original authors are kept. There is NO BM support at this point.
+   Id like to know how Cactus handles this w/o #ifdef !
+c
+c  @endhistory 
+c@@*/
+
+#include "cctk.h"
+#include "declare_arguments.h"
+#include "declare_parameters.h"
+
+      subroutine linearWaves(CCTK_FARGUMENTS)
+      implicit none
+      
+      DECLARE_CCTK_FARGUMENTS
+      DECLARE_PARMETERS
+
+
+c     Call the Teukolsky wave init routine
+      if ((CCTK_Equals(initial_data,'teukwaves').eq.1) then
+         call teukwaves(CCTK_FARGUMENTS)
+
+c     at this point CACTUS 3.2.0 has the following calls
+c     #ifdef BMUTIL call metricderiv
+c     #ifdef BONAMASSO: vxyz=0; call vectorini;
+
+      elseif (CCTK_Equals(initial_data,'planewaves')) then
+         call planewaves(CCTK_FARGUMENTS)
+
+c     at this point CACTUS 3.2.0 has the following calls     
+c     #ifdef BMUTIL call vectorini;
+
+      else 
+         CCTK_Warn(0,"Thorn linearWaves: model initialization failed");
+      endif
+      
+      return
+      end
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
+      
+
+
+
+
+        
diff --git a/src/teukwaves.F b/src/teukwaves.F
new file mode 100644
index 0000000..c7db62a
--- /dev/null
+++ b/src/teukwaves.F
@@ -0,0 +1,569 @@
+#include "cctk.h"
+#include "declare_arguments.h"
+#include "declare_parameters.h"
+
+/*@@
+  @file     teukwaves.F
+  @date     Jan 96
+  @author   Joan Masso + Malcolm Tobias
+  @desc
+            Routines for the teukolsky waves initialization.
+  @enddesc
+
+@@*/
+        
+/*@@
+  @routine  teukwaves
+  @date     Jan 96
+  @author   Malcolm Tobias + Joan Masso 
+  @desc 
+            Generates the analytic solutions of the
+            metric components for linearized quadrupolar (teukolsky) waves.
+            More changes for cactus port: all 3d arrays are converted
+            to point calculations. What a pain!!!
+
+            Three and a half years later....
+            Converted to Cactus4.0
+
+  @enddesc
+@@*/
+
+      subroutine teukwaves(CCTK_FARGUMENTS)
+      implicit none
+
+      DECLARE_CCTK_FARGUMENTS
+      DECLARE_PARAMETERS
+
+      
+      REAL amp,m,ra,pi
+      REAL wave,wave2,wave3,wave4,wave5,wave6,wave7,wave8
+      REAL wavepulse
+      INTEGER ipacket,iparity
+      INTEGER ingoing,outgoing
+      REAL kappa
+c     point values of x.y z,r
+      REAL xp,yp,zp,rp
+
+c     spherical metric
+      REAL grr,grt,grp,gtt,gtp,gpp
+      REAL hrr,hrt,hrp,htt,htp,hpp
+      
+c     special coefficeints needed for teukolsky waves
+c     All these were 3d arrays in old versions of Newage,G,etc..
+c     uffff....
+      REAL teuk_tp,teuk_tn
+      REAL tp2,tn2,fp,fn,fpa,fna,fpb,fnb,fpc,fnc,fpd,fnd,fpe,fne
+      REAL ca,cb,cc,ck,cl,frr,frt,frp,ftt1,ftt2,ftp,fpp1,fpp2
+      REAL cadot,cbdot,ccdot,ckdot,cldot
+      REAL drt,drp,dtt,dtp,dpp,sint,cost,sinp,cosp,tmp
+
+c     from old spheretocart
+      REAL drx,dry,drz,dtx,dty,dtz,dpx,dpy,dpz
+      
+      INTEGER i,j,k
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+      pi    = 3.14159265358979
+      kappa = sqrt(143.0d0/pi)/12288.0d0
+
+      if (wavecenter.ne.0.) then
+         CCTK_Warn(0,"Teukwaves: wavecenter must be zero for time symmetry")
+      endif
+ 
+c     CACTUS 4.0: all REAL/INTEGER parameter are available as variables and 
+c     need not be derived from the database. Nevertheless I will assign 
+c     the parameter variables to the variables of CAC3.0, to allow for 
+c     easy portation (wavepulse stays wavepulse btw). (GERD)
+      
+      am  = amplitude
+      m   = mvalue
+      wave= wavelength
+
+c     to be consistent with Evans & Abrahams, rescale the amplitude
+c     by kappa, because they use that in the packet and
+c     4/3 because they use I not Teuk. F  
+      if (CCTK_Equals(packet,'eppley')) then
+         ipacket = 1
+         write(*,*)'Teuk.Packet= eppley' 
+      elseif (CCTK_Equals(packet,'evans')) then
+         ipacket = 2
+         write(*,*)'Teuk.Packet= evans'
+c     to be consistent with Evans & Abrahams, rescale the amplitude
+c     by kappa, because they use that in the packet and
+c     4/3 because they use I not Teuk. F
+         amp = amp*(4.0d0/3.0d0)*kappa
+      elseif (CCTK_Equals(packet,'square')) then
+         ipacket = 3
+         write(*,*)'Teuk.Packet= square'
+      else
+         write(*,*)'Teukwaves: unsupported packet'
+         STOP
+      endif      
+
+      
+      write(*,*)'amplitude  = ',amp
+      write(*,*)'     m     = ',m
+      write(*,*)'wavecenter = ',ra
+      write(*,*)'wavelength = ',wave
+      
+      if (CCTK_Equals(parity,'even')) then    
+         iparity = 0
+      else (CONTAINS(parity,'odd')) then 
+         iparity = 1
+      endif
+
+       if(CCTK_Equals(wavegoing,'in')) then
+          ingoing = 1
+       endif
+
+      if(CCTK_Equals(wavegoing,'out')) then
+         outgoing = 1
+      endif
+
+      if(CCTK_Equals(wavegoing,'both')) then
+         ingoing  = 1
+         outgoing = 1
+      endif
+
+      if (ipacket.eq.1) then
+         if (ingoing.eq.0.or.outgoing.eq.0) then
+            write(*,*)'WARNING: Epply packet is only non-singular at'
+            write(*,*)'the origin for ingoing-outgoing combination'
+            write(*,*)'of waves.  Expect problems ;-) '
+         endif
+      endif
+      
+      wave2 = wave*wave
+      wave3 = wave2*wave
+      wave4 = wave3*wave
+      wave5 = wave4*wave
+      wave6 = wave5*wave
+      wave7 = wave6*wave
+      wave8 = wave7*wave
+         
+
+c  ****************************************   
+c     initial data for teukolsky waves
+c  ****************************************
+      do k=1,sh(3)
+         do j=1,sh(2)
+            do i=1,sh(1)
+               
+               xp = x(i,j,k)
+               yp = y(i,j,k)
+               zp = z(i,j,k)
+               rp = r(i,j,k)
+               
+               teuk_tp = (time+rp-ra)
+               teuk_tn = (time-rp+ra)
+               
+               tp2 = teuk_tp**2
+               tn2 = teuk_tn**2
+
+c initialize everything to zero
+               fp = 0.0d0
+               fpa = 0.0d0
+               fpb = 0.0d0
+               fpc = 0.0d0
+               fpd = 0.0d0
+               fpe = 0.0d0
+
+               fn = 0.0d0
+               fna = 0.0d0
+               fnb = 0.0d0
+               fnc = 0.0d0
+               fnd = 0.0d0
+               fne = 0.0d0
+               
+c     eppley package --> x exp(-x^2)  
+               if (ipacket.eq.1) then
+c     to keep the analytic solution valid for the eppley packet:   
+                  if (ra.ne.0.and.ingoing.eq.1.and.outgoing.eq.1) then
+                     teuk_tn = (time-rp-ra)
+                     tn2 = teuk_tn**2
+                  endif
+
+                  if (ingoing.eq.1) then
+
+                     tmp = exp(-tp2)
+                     fp = amp * teuk_tp*tmp
+                     fpa = amp * (1 - 2*tp2)*tmp
+                     fpb = amp * teuk_tp*(4*tp2 - 6)*tmp
+                     fpc = amp * (24*tp2 - 8*tp2*tp2 - 6)*tmp
+                     fpd = amp * teuk_tp*(60 - 80*tp2 + 16*tp2*tp2)*tmp
+                     fpe = amp * (
+     &                    (60 - 360*tp2 + 240*tp2*tp2 - 32*tp2*tp2*tp2)*tmp )
+
+                     if (outgoing.eq.0.and.(teuk_tp+wavepulse).le.0) 
+     &                    then
+
+                        fp  = 0.0d0
+                        fpa = 0.0d0
+                        fpb = 0.0d0
+                        fpc = 0.0d0
+                        fpd = 0.0d0
+                        fpe = 0.0d0
+
+                     endif
+
+                  endif
+
+                  if (outgoing.eq.1) then
+
+                     tmp = exp(-tn2)
+                     fn = amp * teuk_tn*tmp
+                     fna = amp * (1 - 2*tn2)*tmp
+                     fnb = amp * teuk_tn*(4*tn2 - 6)*tmp
+                     fnc = amp * (24*tn2 - 8*tn2*tn2 -6)*tmp
+                     fnd = amp * teuk_tn*(60 - 80*tn2 + 16*tn2*tn2)*tmp
+                     fne = amp * (
+     &                    (60 - 360*tn2 + 240*tn2*tn2 - 32*tn2*tn2*tn2)*tmp )
+
+                  endif
+               endif
+
+      
+c     evans package  --> w^4(1-x^2/w^2)^6
+               if(ipacket.eq.2) then
+
+                  if (ingoing.eq.1) then
+
+                     tmp = 1.0d0 - tp2/wave2
+
+                     if (abs(teuk_tp).lt.wave) then
+
+                        fp  =  amp*wave5*tmp**6
+                        fpa =  -12.0d0*amp*wave3*teuk_tp*tmp**5
+                        fpb =  -12.0d0*amp*wave3*(
+     &                       tmp**5 - 10.0d0*tp2/wave2*tmp**4 )
+                        fpc =  120.0d0*amp*wave3*(
+     &                        3.0d0*teuk_tp/wave2*tmp**4 
+     &                       -8.0d0*teuk_tp*tp2/wave4*tmp**3 )
+                        fpd = 360.0d0*amp*wave3*(
+     &                        1.0d0/wave2*tmp**4 
+     &                       -16.0d0*tp2/wave4*tmp**3
+     &                       +16.0d0*tp2*tp2/wave6*tmp**2 )
+                        fpe = 2880.0d0*amp*wave3*(
+     &                       -5.0d0*teuk_tp/wave4*tmp**3
+     &                       +20.0d0*teuk_tp*tp2/wave6*tmp**2
+     &                       -8.0d0*teuk_tp*tp2*tp2/wave8*tmp )
+
+                     else
+                        fp=0.
+                        fpa=0.
+                        fpb=0.
+                        fpc=0.
+                        fpd=0.
+                        fpe=0.
+                     endif
+
+                  endif
+
+                  if (outgoing.eq.1) then
+
+                     tmp = 1.0d0 - tn2/wave2
+         
+                     if (abs(teuk_tn).lt.wave) then
+
+                        fn  =  amp*wave5*tmp**6
+                        fna =  -12.0d0*amp*wave3*teuk_tn*tmp**5
+                        fnb =  -12.0d0*amp*wave3*(
+     &                       tmp**5 - 10.0d0*tn2/wave2*tmp**4 )
+                        fnc =  120.0d0*amp*wave3*(
+     &                        3.0d0*teuk_tn/wave2*tmp**4 
+     &                       -8.0d0*teuk_tn*tn2/wave4*tmp**3 )
+                        fnd = 360.0d0*amp*wave3*(
+     &                        1.0d0/wave2*tmp**4 
+     &                       -16.0d0*tn2/wave4*tmp**3
+     &                       +16.0d0*tn2*tn2/wave6*tmp**2 )
+                        fne = 2880.0d0*amp*wave3*(
+     &                       -5.0d0*teuk_tn/wave4*tmp**3
+     &                       +20.0d0*teuk_tn*tn2/wave6*tmp**2
+     &                       -8.0d0*teuk_tn*tn2*tn2/wave8*tmp )
+
+                     else
+                        fn=0.
+                        fna=0.
+                        fnb=0.
+                        fnc=0.
+                        fnd=0.
+                        fne=0.
+                     endif
+
+                  endif
+
+               endif
+      
+c     square package --> (1-x^2/w^2)^2
+               if(ipacket.eq.3) then 
+
+                  write(*,*)'WARNING: need to calculate fpe,fne'
+                  write(*,*)'WARNING: need conditionals for in out waves'
+
+                  if (abs(teuk_tp).lt.wave) then
+                     fp = amp * (1-(tp2/wave2))**2
+                     fpa = amp * (-4*teuk_tp)*(1-tp2/wave2)/wave2
+                     fpb = amp * (8*tp2/wave4-4*(1-tp2/wave2)/wave2)
+                     fpc = amp * 24*teuk_tp/wave4
+                     fpd = amp * 24/wave4
+                  else
+                     fp=0.
+                     fpa=0.
+                     fpb=0.
+                     fpc=0.
+                     fpd=0.
+                  endif
+         
+                  if (abs(teuk_tn).lt.wave) then
+                     fn = amp * (1-(tn2/wave2))**2
+                     fna = amp * (-4*teuk_tn)*(1-tn2/wave2)/wave2
+                     fnb = amp * (8*tn2/wave4-4*(1-tn2/wave2)/wave2)
+                     fnc = amp * 24*teuk_tn/wave4
+                     fnd = amp * 24/wave4  
+                  else
+                     fn=0.
+                     fna=0.
+                     fnb=0.
+                     fnc=0.
+                     fnd=0.
+                  endif
+               endif
+      
+c     
+c     coefficients
+c     
+               ca = 3*( (fpb-fnb)/rp**3 -3*(fna+fpa)/rp**4 +3*(fp-fn)/rp**5 )
+               cb = -( -(fnc+fpc)/rp**2 +3*(fpb-fnb)/rp**3 -6*(fna+fpa)/rp**4
+     &              +6*(fp-fn)/rp**5 ) 
+               cc = ( (fpd-fnd)/rp -2*(fnc+fpc)/rp**2 +9*(fpb-fnb)/rp**3
+     &              -21*(fna+fpa)/rp**4 +21*(fp-fn)/rp**5 )/4
+               ck = (fpb-fnb)/rp**2-3*(fpa+fna)/rp**3+3*(fp-fn)/rp**4
+               cl = -(fpc+fnc)/rp+2*(fpb-fnb)/rp**2-3*(fpa+fna)/rp**3+
+     &              3*(fp-fn)/rp**4
+               
+c     If we only have an outgoing wave, we should flip the sign
+c     of the above coefficients.  That is because the above coefficients
+c     assume an ingoing and outgoing wave packet of the form:
+c     wave = ingoing - outgoing.  For a single wave, we want ingoing
+c     and outgoing to have the same sign.
+
+               if(ingoing.eq.0.and.outgoing.eq.1) then
+                  ca=-ca
+                  cb=-cb
+                  cc=-cc
+                  ck=-ck
+                  cl=-cl
+               endif
+               
+               sint = (xp**2+yp**2)**0.5/rp
+               cost = zp/rp
+               sinp = yp/(xp**2+yp**2)**0.5
+               cosp = xp/(xp**2+yp**2)**0.5
+               
+c     mvalue
+               if (m.eq.0) then
+                  frr = 2-3*sint**2
+                  frt = -3*sint*cost
+                  frp = 0.
+                  ftt1 = 3*sint**2
+                  ftt2 = -1.
+                  ftp = 0.
+                  fpp1 = -ftt1
+                  fpp2 = 3*sint**2-1
+                  drt = 0.
+                  drp = -4*cost*sint
+                  dtt = 0.
+                  dtp = -sint**2
+                  dpp = 0.
+               elseif (m.eq.-1) then
+                  frr = 2*sint*cost*sinp
+                  frt = (cost**2-sint**2)*sinp 
+                  frp = cost*cosp 
+                  ftt1 = -frr
+                  ftt2 = 0.
+                  ftp = sint*cosp 
+                  fpp1 = -ftt1
+                  fpp2 = ftt1
+                  drt = -2*cost*sinp
+                  drp = -2*(cost**2-sint**2)*cosp
+                  dtt = -sint*sinp
+                  dtp = -cost*sint*cosp
+                  dpp = sint*sinp
+               elseif (m.eq.-2) then
+                  frr = sint**2*2*sinp*cosp
+                  frt = sint*cost*2*sinp*cosp
+                  frp = sint*(1-2*sinp**2)
+                  ftt1 = (1+cost**2)*2*sinp*cosp
+                  ftt2 = -2*sinp*cosp
+                  fpp1 = -ftt1
+                  fpp2 = cost**2*2*sinp*cosp
+                  drt = 8*sint*sinp*cosp
+                  drp = 4*sint*cost*(1-2*sinp**2)
+                  dtt = -4*cost*sinp*cosp
+                  dtp = (2-sint**2)*(2*sinp**2-1)
+                  dpp = 2*cost*2*sinp*cosp
+               elseif (m.eq.1) then
+                  frr = 2*sint*cost*cosp
+                  frt = (cost**2-sint**2)*cosp
+                  frp = -cost*sinp
+                  ftt1 = -2*sint*cost*cosp
+                  ftt2 = 0.
+                  ftp = -sint*sinp
+                  fpp1 = -ftt1
+                  fpp2 = ftt1
+                  drt = -2*cost*cosp
+                  drp = 2*(cost**2-sint**2)*sinp
+                  dtt = -sint*cosp
+                  dtp = cost*sint*sinp
+                  dpp = sint*cosp
+               elseif (m.eq.2) then
+                  frr = sint**2*(1-2*sinp**2)
+                  frt = sint*cost*(1-2*sinp**2)
+                  frp = -sint*2*sinp*cosp
+                  ftt1 = (1+cost**2)*(1-2*sinp**2)
+                  ftt2 = (2*sinp**2-1)
+                  ftp = cost*2*sinp*cosp
+                  fpp1 = -ftt1
+                  fpp2 = cost**2*(1-2*sinp**2)
+                  drt = 4*sint*(1-2*sinp**2)
+                  drp = -4*sint*cost*2*sinp*cosp
+                  dtt = -2*cost*(1-2*sinp**2)
+                  dtp = (2-sint**2)*2*sinp*cosp
+                  dpp = 2*cost*(1-2*sinp**2)
+               else
+                  write(*,*) 'teukwaves: m should be one of -2,-1,0,1,2'
+                  STOP
+               endif
+               
+c     parity
+               if (iparity.eq.0) then    
+                  grr = 1 + ca*frr
+                  grt = cb*frt*rp
+                  grp = cb*frp*(xp**2+yp**2)**0.5
+                  gtt = rp**2*(1 + cc*ftt1 + ca*ftt2)
+                  gtp = (ca - 2*cc)*ftp*rp*(xp**2+yp**2)**0.5
+                  gpp = (xp**2+yp**2)*(1 + cc*fpp1 + ca*fpp2)
+               else
+                  grr = 1.
+                  grt = ck*drt*rp
+                  grp = ck*drp*rp*sint
+                  gtt = (1+cl*dtt)*rp**2
+                  gtp = cl*dtp*rp**2*sint
+                  gpp = (1+cl*dpp)*rp**2*sint**2
+               endif
+
+c     Finally, convert the spherical components to cartesian
+
+c     define derivatives drx = (dr/dx)               
+               drx = xp/rp
+               dry = yp/rp
+               drz = zp/rp
+      
+               dtx = xp*zp/(rp**2*sqrt(xp**2+yp**2))
+               dty = yp*zp/(rp**2*sqrt(xp**2+yp**2))
+               dtz = (zp**2/rp**2-1)/sqrt(xp**2+yp**2)
+               
+               dpx = -yp/(xp**2+yp**2)
+               dpy = xp/(xp**2+yp**2)
+               dpz = 0.
+
+c     cartesian g
+               gxx(i,j,k) = drx**2*grr+2.*drx*dtx*grt
+     &              +2.*drx*dpx*grp+dtx**2*gtt
+     &              +2.*dtx*dpx*gtp+dpx**2*gpp
+               
+               gyy(i,j,k) = dry**2*grr+2.*dry*dty*grt
+     &              +2.*dry*dpy*grp+dty**2*gtt
+     &              +2.*dty*dpy*gtp+dpy**2*gpp
+ 
+               gzz(i,j,k) = drz**2*grr+2.*drz*dtz*grt
+     &              +2.*drz*dpz*grp+dtz**2*gtt
+     &              +2.*dtz*dpz*gtp+dpz**2*gpp
+               
+               gxy(i,j,k) = drx*dry*grr+(drx*dty+dtx*dry)*grt
+     &              +(drx*dpy+dpx*dry)*grp
+     &              +dtx*dty*gtt+(dtx*dpy+dpx*dty)*gtp+dpx*dpy*gpp
+               
+               gxz(i,j,k) = drx*drz*grr+(drx*dtz+dtx*drz)*grt
+     &              +(drx*dpz+dpx*drz)*grp
+     &              +dtx*dtz*gtt+(dtx*dpz+dpx*dtz)*gtp+dpx*dpz*gpp
+               
+               gyz(i,j,k) = dry*drz*grr+(dry*dtz+dty*drz)*grt
+     &              +(dry*dpz+dpy*drz)*grp
+     &              +dty*dtz*gtt+(dty*dpz+dpy*dtz)*gtp+dpy*dpz*gpp
+                
+
+
+c     for the ext. curv. well need...
+               cadot = 3*( (fpc-fnc)/rp**3 -3*(fnb+fpb)/rp**4 
+     &              +3*(fpa-fna)/rp**5 )
+               cbdot = -( -(fnd+fpd)/rp**2 +3*(fpc-fnc)/rp**3 
+     &              -6*(fnb+fpb)/rp**4
+     &              +6*(fpa-fna)/rp**5 ) 
+               ccdot = ( (fpe-fne)/rp -2*(fnd+fpd)/rp**2 +9*(fpc-fnc)/rp**3
+     &              -21*(fnb+fpb)/rp**4 +21*(fpa-fna)/rp**5 )/4
+               ckdot = (fpc-fnc)/rp**2-3*(fpb+fnb)/rp**3+3*(fpa-fna)/rp**4
+               cldot = -(fpd+fnd)/rp+2*(fpc-fnc)/rp**2-3*(fpb+fnb)/rp**3+
+     &              3*(fpa-fna)/rp**4
+
+c     parity
+               if (iparity.eq.0) then    
+                  hrr = -0.5*(cadot*frr)
+                  hrt = -0.5*(cbdot*frt*rp)
+                  hrp = -0.5*(cbdot*frp*(xp**2+yp**2)**0.5)
+                  htt = -0.5*(rp**2*(ccdot*ftt1 + cadot*ftt2))
+                  htp = -0.5*((cadot - 2*ccdot)*ftp*rp*(xp**2+yp**2)**0.5)
+                  hpp = -0.5*((xp**2+yp**2)*(ccdot*fpp1 + cadot*fpp2))
+               else
+                  hrr = 0.
+                  hrt = -0.5*(ckdot*drt*rp)
+                  hrp = -0.5*(ckdot*drp*rp*sint)
+                  htt = -0.5*((cldot*dtt)*rp**2)
+                  htp = -0.5*(cldot*dtp*rp**2*sint)
+                  hpp = -0.5*((cldot*dpp)*rp**2*sint**2)
+               endif
+               
+c     time symmetry
+               hxx(i,j,k) = drx**2*hrr+2.*drx*dtx*hrt
+     &              +2.*drx*dpx*hrp+dtx**2*htt
+     &              +2.*dtx*dpx*htp+dpx**2*hpp
+
+               hyy(i,j,k) = dry**2*hrr+2.*dry*dty*hrt
+     &              +2.*dry*dpy*hrp+dty**2*htt
+     &              +2.*dty*dpy*htp+dpy**2*hpp
+               
+               hzz(i,j,k) = drz**2*hrr+2.*drz*dtz*hrt
+     &              +2.*drz*dpz*hrp+dtz**2*htt
+     &              +2.*dtz*dpz*htp+dpz**2*hpp
+               
+               hxy(i,j,k) = drx*dry*hrr+(drx*dty+dtx*dry)*hrt
+     &              +(drx*dpy+dpx*dry)*hrp
+     &              +dtx*dty*htt+(dtx*dpy+dpx*dty)*htp+dpx*dpy*hpp
+               
+               hxz(i,j,k) = drx*drz*hrr+(drx*dtz+dtx*drz)*hrt
+     &              +(drx*dpz+dpx*drz)*hrp
+     &              +dtx*dtz*htt+(dtx*dpz+dpx*dtz)*htp+dpx*dpz*hpp
+               
+               hyz(i,j,k) = dry*drz*hrr+(dry*dtz+dty*drz)*hrt
+     &              +(dry*dpz+dpy*drz)*hrp
+     &              +dty*dtz*htt+(dty*dpz+dpy*dtz)*htp+dpy*dpz*hpp
+                
+
+               
+c$$$               hxx(i,j,k) = 0.
+c$$$               hxy(i,j,k) = 0.
+c$$$               hxz(i,j,k) = 0.
+c$$$               hyy(i,j,k) = 0.
+c$$$               hyz(i,j,k) = 0.
+c$$$               hzz(i,j,k) = 0.
+
+
+            enddo
+         enddo
+      enddo
+
+ 
+      return
+      end
+