From d4ab2056efeb6149e0d87330e5eb182152a47b10 Mon Sep 17 00:00:00 2001 From: allen Date: Wed, 19 Sep 2001 17:00:21 +0000 Subject: Changing to Fortran 77 git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/IDLinearWaves/trunk@63 5c0f84ea-6048-4d6e-bfa4-55cd5f2e0dd7 --- src/planewaves.F | 256 ------------------------------------------------------- 1 file changed, 256 deletions(-) delete mode 100644 src/planewaves.F (limited to 'src/planewaves.F') diff --git a/src/planewaves.F b/src/planewaves.F deleted file mode 100644 index a0bf788..0000000 --- a/src/planewaves.F +++ /dev/null @@ -1,256 +0,0 @@ -/*@@ - @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 "cctk_Arguments.h" -#include "cctk_Parameters.h" -c Using macro definitions from Einstein -#include "CactusEinstein/Einstein/src/Einstein.h" - - - subroutine planewaves(CCTK_ARGUMENTS) - implicit none - - DECLARE_CCTK_ARGUMENTS - DECLARE_CCTK_PARAMETERS - - INTEGER iin,iout - CCTK_REAL pi - CCTK_REAL ra,the,phi - CCTK_REAL wave,wavep - CCTK_REAL kx,ky,kz,w - CCTK_REAL kxp,kyp,kzp,wp - CHARACTER*200 :: infoline - -c local arrays (scalars! Man, what a sweat!) - CCTK_REAL plus,minus,plusp,minusp,ain,aout - - INTEGER i,j,k - - INTEGER CCTK_Equals - - 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 - - call CCTK_INFO('Plane waves') - write(infoline,'(A13,G12.7)') - & ' amplitude = ',amplitude - call CCTK_INFO(infoline) - write(infoline,'(A14,G12.7)') - & ' wavecenter = ',ra - call CCTK_INFO(infoline) - write(infoline,'(A14,G12.7)') - & ' wavelength = ',wave - call CCTK_INFO(infoline) - write(infoline,'(A14,G12.7)') - & ' wavepulse = ',wavep - call CCTK_INFO(infoline) - -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,cctk_lsh(3) - do j=1,cctk_lsh(2) - do i=1,cctk_lsh(1) - plus = (kx*x(i,j,k)+ky*y(i,j,k)+kz*z(i,j,k)+w*cctk_time) - minus = (kx*x(i,j,k)+ky*y(i,j,k)+kz*z(i,j,k)-w*cctk_time) - - plusp =(kxp*x(i,j,k)+kyp*y(i,j,k)+kzp*z(i,j,k)+wp*(cctk_time-ra)) - minusp =(kxp*x(i,j,k)+kyp*y(i,j,k)+kzp*z(i,j,k)-wp*(cctk_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. - - kxx(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) - - - kxy(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) - - - kxz(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) - - - kyy(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) - - - kyz(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) - - - kzz(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 - -c initialize the conformal factor - if (use_conformal == 1) then - conformal_state = CONFORMAL_METRIC - psi = 1d0 - psix = 0d0 - psiy = 0d0 - psiz = 0d0 - psixy = 0d0 - psixz = 0d0 - psiyz = 0d0 - psixx = 0d0 - psiyy = 0d0 - psizz = 0d0 - else - conformal_state = NOCONFORMAL_METRIC - end if - - - return - end - - - - - - -- cgit v1.2.3