From 145184eb4a4321d132384d9f3f4ed3aad31968e3 Mon Sep 17 00:00:00 2001 From: lanfer Date: Tue, 16 Mar 1999 13:58:06 +0000 Subject: the linear wave routines git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/IDLinearWaves/trunk@3 5c0f84ea-6048-4d6e-bfa4-55cd5f2e0dd7 --- src/planewaves.F | 224 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 224 insertions(+) create mode 100644 src/planewaves.F (limited to 'src/planewaves.F') 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:

+ $ A*exp\left[-(kp_ix^i-\omega_p (time-ra))^2\right] + cos(k_ix^i-\omega \ time) $ +

+ where:

+ A = amplitude of the wave
+ k = the wave # of the sine wave
+ $\omega$ = the frequency of the sine wave
+ kp = the wave # of the gaussian modulating the sine wave
+ $\omega_p$ = the frequency of the gaussian
+ 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 + + + + + + -- cgit v1.2.3