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/make.code.defn | 4 +- src/planewaves.F | 256 ----------------------- src/planewaves.F77 | 264 +++++++++++++++++++++++ src/teukwaves.F | 592 ---------------------------------------------------- src/teukwaves.F77 | 601 +++++++++++++++++++++++++++++++++++++++++++++++++++++ 5 files changed, 867 insertions(+), 850 deletions(-) delete mode 100644 src/planewaves.F create mode 100644 src/planewaves.F77 delete mode 100644 src/teukwaves.F create mode 100644 src/teukwaves.F77 diff --git a/src/make.code.defn b/src/make.code.defn index 609248f..a152d03 100644 --- a/src/make.code.defn +++ b/src/make.code.defn @@ -3,8 +3,8 @@ # Source files in this directory SRCS = \ - planewaves.F\ - teukwaves.F + planewaves.F77\ + teukwaves.F77 # Subdirectories containing source files SUBDIRS = 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 - - - - - - diff --git a/src/planewaves.F77 b/src/planewaves.F77 new file mode 100644 index 0000000..79c1d4e --- /dev/null +++ b/src/planewaves.F77 @@ -0,0 +1,264 @@ +/*@@ + @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 + do k=1,cctk_lsh(3) + do j=1,cctk_lsh(2) + do i=1,cctk_lsh(1) + + psi(i,j,k) = 1d0 + psix(i,j,k) = 0d0 + psiy(i,j,k) = 0d0 + psiz(i,j,k) = 0d0 + 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 + + end do + end do + end do + else + conformal_state = NOCONFORMAL_METRIC + end if + + + return + end + + + + + + diff --git a/src/teukwaves.F b/src/teukwaves.F deleted file mode 100644 index 02ff0d7..0000000 --- a/src/teukwaves.F +++ /dev/null @@ -1,592 +0,0 @@ -#include "cctk.h" -#include "cctk_Arguments.h" -#include "cctk_Parameters.h" -c Using macro definitions from Einstein -#include "CactusEinstein/Einstein/src/Einstein.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_ARGUMENTS) - implicit none - - DECLARE_CCTK_ARGUMENTS - DECLARE_CCTK_PARAMETERS - - CCTK_REAL amp,m,ra,pi - CCTK_REAL wave,wave2,wave3,wave4,wave5,wave6,wave7,wave8 - INTEGER ipacket,iparity - INTEGER ingoing,outgoing - CCTK_REAL kappa -c point values of x.y z,r - CCTK_REAL xp,yp,zp,rp - -c spherical metric - CCTK_REAL teuk_grr,teuk_grt,teuk_grp,teuk_gtt,teuk_gtp,teuk_gpp - CCTK_REAL teuk_hrr,teuk_hrt,teuk_hrp,teuk_htt,teuk_htp,teuk_hpp - -c special coefficeints needed for teukolsky waves -c All these were 3d arrays in old versions of Newage,G,etc.. -c uffff.... - CCTK_REAL teuk_tp,teuk_tn - CCTK_REAL tp2,tn2,fp,fn,fpa,fna,fpb,fnb,fpc,fnc,fpd,fnd,fpe,fne - CCTK_REAL ca,cb,cc,ck,cl,frr,frt,frp,ftt1,ftt2,ftp,fpp1,fpp2 - CCTK_REAL cadot,cbdot,ccdot,ckdot,cldot - CCTK_REAL drt,drp,dtt,dtp,dpp,sint,cost,sinp,cosp,tmp - -c from old spheretocart - CCTK_REAL drx,dry,drz,dtx,dty,dtz,dpx,dpy,dpz - - CHARACTER*200 infoline - INTEGER i,j,k - INTEGER CCTK_Equals - -cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - - pi = 3.14159265358979 - kappa = sqrt(143.0d0/pi)/12288.0d0 - - if (wavecenter.ne.0.) then - call 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) - - amp = amplitude - m = mvalue - wave = wavelength - ra = wavecenter - -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') == 1) then - ipacket = 1 - call CCTK_INFO('Teukolsky Packet = eppley') - elseif (CCTK_Equals(packet,'evans') == 1) then - ipacket = 2 - call CCTK_INFO('Teukolsky 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') == 1) then - ipacket = 3 - call CCTK_INFO('Teukolsky Packet = square') - end if - - write(infoline,'(A13,G12.7)') - & ' amplitude = ',amp - call CCTK_INFO(infoline) - write(infoline,'(A11,G12.7)') - & ' m value = ',m - 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) - - if (CCTK_Equals(parity,'even') == 1) then - iparity = 0 - elseif (CCTK_Equals(parity,'odd') == 1) then - iparity = 1 - endif - - if(CCTK_Equals(wavesgoing,'in') == 1) then - ingoing = 1 - endif - - if(CCTK_Equals(wavesgoing,'out') == 1) then - outgoing = 1 - endif - - if(CCTK_Equals(wavesgoing,'both') == 1) then - ingoing = 1 - outgoing = 1 - endif - - if (ipacket.eq.1) then - if (ingoing.eq.0.or.outgoing.eq.0) then - call CCTK_WARN(4,'Epply packet is only non-singular at the origin for ingoing-outgoing combination of waves') - 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,cctk_lsh(3) - do j=1,cctk_lsh(2) - do i=1,cctk_lsh(1) - - xp = x(i,j,k) - yp = y(i,j,k) - zp = z(i,j,k) - rp = r(i,j,k) - - teuk_tp = (cctk_time+rp-ra) - teuk_tn = (cctk_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 = (cctk_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 - - call CCTK_WARN(4,'Need to calculate fpe,fne') - call CCTK_WARN(4,'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 - call CCTK_WARN(0,'m should be one of -2,-1,0,1,2') - endif - -c parity - if (iparity.eq.0) then - teuk_grr = 1 + ca*frr - teuk_grt = cb*frt*rp - teuk_grp = cb*frp*(xp**2+yp**2)**0.5 - teuk_gtt = rp**2*(1 + cc*ftt1 + ca*ftt2) - teuk_gtp = (ca - 2*cc)*ftp*rp*(xp**2+yp**2)**0.5 - teuk_gpp = (xp**2+yp**2)*(1 + cc*fpp1 + ca*fpp2) - else - teuk_grr = 1. - teuk_grt = ck*drt*rp - teuk_grp = ck*drp*rp*sint - teuk_gtt = (1+cl*dtt)*rp**2 - teuk_gtp = cl*dtp*rp**2*sint - teuk_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*teuk_grr+2.*drx*dtx*teuk_grt - & +2.*drx*dpx*teuk_grp+dtx**2*teuk_gtt - & +2.*dtx*dpx*teuk_gtp+dpx**2*teuk_gpp - - gyy(i,j,k) = dry**2*teuk_grr+2.*dry*dty*teuk_grt - & +2.*dry*dpy*teuk_grp+dty**2*teuk_gtt - & +2.*dty*dpy*teuk_gtp+dpy**2*teuk_gpp - - gzz(i,j,k) = drz**2*teuk_grr+2.*drz*dtz*teuk_grt - & +2.*drz*dpz*teuk_grp+dtz**2*teuk_gtt - & +2.*dtz*dpz*teuk_gtp+dpz**2*teuk_gpp - - gxy(i,j,k) = drx*dry*teuk_grr+(drx*dty+dtx*dry)*teuk_grt - & +(drx*dpy+dpx*dry)*teuk_grp - & +dtx*dty*teuk_gtt+(dtx*dpy+dpx*dty)*teuk_gtp+ - & dpx*dpy*teuk_gpp - - gxz(i,j,k) = drx*drz*teuk_grr+(drx*dtz+dtx*drz)*teuk_grt - & +(drx*dpz+dpx*drz)*teuk_grp - & +dtx*dtz*teuk_gtt+(dtx*dpz+dpx*dtz)*teuk_gtp+ - & dpx*dpz*teuk_gpp - - gyz(i,j,k) = dry*drz*teuk_grr+(dry*dtz+dty*drz)*teuk_grt - & +(dry*dpz+dpy*drz)*teuk_grp - & +dty*dtz*teuk_gtt+(dty*dpz+dpy*dtz)*teuk_gtp+ - & dpy*dpz*teuk_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 - teuk_hrr = -0.5*(cadot*frr) - teuk_hrt = -0.5*(cbdot*frt*rp) - teuk_hrp = -0.5*(cbdot*frp*(xp**2+yp**2)**0.5) - teuk_htt = -0.5*(rp**2*(ccdot*ftt1 + cadot*ftt2)) - teuk_htp = - & -0.5*((cadot - 2*ccdot)*ftp*rp*(xp**2+yp**2)**0.5) - teuk_hpp = -0.5*((xp**2+yp**2)*(ccdot*fpp1 + cadot*fpp2)) - else - teuk_hrr = 0. - teuk_hrt = -0.5*(ckdot*drt*rp) - teuk_hrp = -0.5*(ckdot*drp*rp*sint) - teuk_htt = -0.5*((cldot*dtt)*rp**2) - teuk_htp = -0.5*(cldot*dtp*rp**2*sint) - teuk_hpp = -0.5*((cldot*dpp)*rp**2*sint**2) - endif - -c time symmetry - kxx(i,j,k) = drx**2*teuk_hrr+2.*drx*dtx*teuk_hrt - & +2.*drx*dpx*teuk_hrp+dtx**2*teuk_htt - & +2.*dtx*dpx*teuk_htp+dpx**2*teuk_hpp - - kyy(i,j,k) = dry**2*teuk_hrr+2.*dry*dty*teuk_hrt - & +2.*dry*dpy*teuk_hrp+dty**2*teuk_htt - & +2.*dty*dpy*teuk_htp+dpy**2*teuk_hpp - - kzz(i,j,k) = drz**2*teuk_hrr+2.*drz*dtz*teuk_hrt - & +2.*drz*dpz*teuk_hrp+dtz**2*teuk_htt - & +2.*dtz*dpz*teuk_htp+dpz**2*teuk_hpp - - kxy(i,j,k) = drx*dry*teuk_hrr+(drx*dty+dtx*dry)*teuk_hrt - & +(drx*dpy+dpx*dry)*teuk_hrp - & +dtx*dty*teuk_htt+(dtx*dpy+dpx*dty)*teuk_htp+ - & dpx*dpy*teuk_hpp - - kxz(i,j,k) = drx*drz*teuk_hrr+(drx*dtz+dtx*drz)*teuk_hrt - & +(drx*dpz+dpx*drz)*teuk_hrp - & +dtx*dtz*teuk_htt+(dtx*dpz+dpx*dtz)*teuk_htp+ - & dpx*dpz*teuk_hpp - - kyz(i,j,k) = dry*drz*teuk_hrr+(dry*dtz+dty*drz)*teuk_hrt - & +(dry*dpz+dpy*drz)*teuk_hrp - & +dty*dtz*teuk_htt+(dty*dpz+dpy*dtz)*teuk_htp+ - & dpy*dpz*teuk_hpp - - - - 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 - - - diff --git a/src/teukwaves.F77 b/src/teukwaves.F77 new file mode 100644 index 0000000..9e02eb0 --- /dev/null +++ b/src/teukwaves.F77 @@ -0,0 +1,601 @@ +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Parameters.h" +c Using macro definitions from Einstein +#include "CactusEinstein/Einstein/src/Einstein.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_ARGUMENTS) + implicit none + + DECLARE_CCTK_ARGUMENTS + + CCTK_REAL amp,m,ra,pi + CCTK_REAL wave,wave2,wave3,wave4,wave5,wave6,wave7,wave8 + INTEGER ipacket,iparity + INTEGER ingoing,outgoing + CCTK_REAL kappa +c point values of x.y z,r + CCTK_REAL xp,yp,zp,rp + +c spherical metric + CCTK_REAL teuk_grr,teuk_grt,teuk_grp,teuk_gtt,teuk_gtp,teuk_gpp + CCTK_REAL teuk_hrr,teuk_hrt,teuk_hrp,teuk_htt,teuk_htp,teuk_hpp + +c special coefficeints needed for teukolsky waves +c All these were 3d arrays in old versions of Newage,G,etc.. +c uffff.... + CCTK_REAL teuk_tp,teuk_tn + CCTK_REAL tp2,tn2,fp,fn,fpa,fna,fpb,fnb,fpc,fnc,fpd,fnd,fpe,fne + CCTK_REAL ca,cb,cc,ck,cl,frr,frt,frp,ftt1,ftt2,ftp,fpp1,fpp2 + CCTK_REAL cadot,cbdot,ccdot,ckdot,cldot + CCTK_REAL drt,drp,dtt,dtp,dpp,sint,cost,sinp,cosp,tmp + +c from old spheretocart + CCTK_REAL drx,dry,drz,dtx,dty,dtz,dpx,dpy,dpz + + CHARACTER*200 infoline + INTEGER i,j,k + INTEGER CCTK_Equals + + DECLARE_CCTK_PARAMETERS + +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + + pi = 3.14159265358979 + kappa = sqrt(143.0d0/pi)/12288.0d0 + + if (wavecenter.ne.0.) then + call 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) + + amp = amplitude + m = mvalue + wave = wavelength + ra = wavecenter + +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') == 1) then + ipacket = 1 + call CCTK_INFO('Teukolsky Packet = eppley') + elseif (CCTK_Equals(packet,'evans') == 1) then + ipacket = 2 + call CCTK_INFO('Teukolsky 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') == 1) then + ipacket = 3 + call CCTK_INFO('Teukolsky Packet = square') + end if + + write(infoline,'(A13,G12.7)') + & ' amplitude = ',amp + call CCTK_INFO(infoline) + write(infoline,'(A11,G12.7)') + & ' m value = ',m + 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) + + if (CCTK_Equals(parity,'even') == 1) then + iparity = 0 + elseif (CCTK_Equals(parity,'odd') == 1) then + iparity = 1 + endif + + if(CCTK_Equals(wavesgoing,'in') == 1) then + ingoing = 1 + endif + + if(CCTK_Equals(wavesgoing,'out') == 1) then + outgoing = 1 + endif + + if(CCTK_Equals(wavesgoing,'both') == 1) then + ingoing = 1 + outgoing = 1 + endif + + if (ipacket.eq.1) then + if (ingoing.eq.0.or.outgoing.eq.0) then + call CCTK_WARN(4,'Epply packet is only non-singular at the origin for ingoing-outgoing combination of waves') + 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,cctk_lsh(3) + do j=1,cctk_lsh(2) + do i=1,cctk_lsh(1) + + xp = x(i,j,k) + yp = y(i,j,k) + zp = z(i,j,k) + rp = r(i,j,k) + + teuk_tp = (cctk_time+rp-ra) + teuk_tn = (cctk_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 = (cctk_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 + + call CCTK_WARN(4,'Need to calculate fpe,fne') + call CCTK_WARN(4,'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 + call CCTK_WARN(0,'m should be one of -2,-1,0,1,2') + endif + +c parity + if (iparity.eq.0) then + teuk_grr = 1 + ca*frr + teuk_grt = cb*frt*rp + teuk_grp = cb*frp*(xp**2+yp**2)**0.5 + teuk_gtt = rp**2*(1 + cc*ftt1 + ca*ftt2) + teuk_gtp = (ca - 2*cc)*ftp*rp*(xp**2+yp**2)**0.5 + teuk_gpp = (xp**2+yp**2)*(1 + cc*fpp1 + ca*fpp2) + else + teuk_grr = 1. + teuk_grt = ck*drt*rp + teuk_grp = ck*drp*rp*sint + teuk_gtt = (1+cl*dtt)*rp**2 + teuk_gtp = cl*dtp*rp**2*sint + teuk_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*teuk_grr+2.*drx*dtx*teuk_grt + & +2.*drx*dpx*teuk_grp+dtx**2*teuk_gtt + & +2.*dtx*dpx*teuk_gtp+dpx**2*teuk_gpp + + gyy(i,j,k) = dry**2*teuk_grr+2.*dry*dty*teuk_grt + & +2.*dry*dpy*teuk_grp+dty**2*teuk_gtt + & +2.*dty*dpy*teuk_gtp+dpy**2*teuk_gpp + + gzz(i,j,k) = drz**2*teuk_grr+2.*drz*dtz*teuk_grt + & +2.*drz*dpz*teuk_grp+dtz**2*teuk_gtt + & +2.*dtz*dpz*teuk_gtp+dpz**2*teuk_gpp + + gxy(i,j,k) = drx*dry*teuk_grr+(drx*dty+dtx*dry)*teuk_grt + & +(drx*dpy+dpx*dry)*teuk_grp + & +dtx*dty*teuk_gtt+(dtx*dpy+dpx*dty)*teuk_gtp+ + & dpx*dpy*teuk_gpp + + gxz(i,j,k) = drx*drz*teuk_grr+(drx*dtz+dtx*drz)*teuk_grt + & +(drx*dpz+dpx*drz)*teuk_grp + & +dtx*dtz*teuk_gtt+(dtx*dpz+dpx*dtz)*teuk_gtp+ + & dpx*dpz*teuk_gpp + + gyz(i,j,k) = dry*drz*teuk_grr+(dry*dtz+dty*drz)*teuk_grt + & +(dry*dpz+dpy*drz)*teuk_grp + & +dty*dtz*teuk_gtt+(dty*dpz+dpy*dtz)*teuk_gtp+ + & dpy*dpz*teuk_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 + teuk_hrr = -0.5*(cadot*frr) + teuk_hrt = -0.5*(cbdot*frt*rp) + teuk_hrp = -0.5*(cbdot*frp*(xp**2+yp**2)**0.5) + teuk_htt = -0.5*(rp**2*(ccdot*ftt1 + cadot*ftt2)) + teuk_htp = + & -0.5*((cadot - 2*ccdot)*ftp*rp*(xp**2+yp**2)**0.5) + teuk_hpp = -0.5*((xp**2+yp**2)*(ccdot*fpp1 + cadot*fpp2)) + else + teuk_hrr = 0. + teuk_hrt = -0.5*(ckdot*drt*rp) + teuk_hrp = -0.5*(ckdot*drp*rp*sint) + teuk_htt = -0.5*((cldot*dtt)*rp**2) + teuk_htp = -0.5*(cldot*dtp*rp**2*sint) + teuk_hpp = -0.5*((cldot*dpp)*rp**2*sint**2) + endif + +c time symmetry + kxx(i,j,k) = drx**2*teuk_hrr+2.*drx*dtx*teuk_hrt + & +2.*drx*dpx*teuk_hrp+dtx**2*teuk_htt + & +2.*dtx*dpx*teuk_htp+dpx**2*teuk_hpp + + kyy(i,j,k) = dry**2*teuk_hrr+2.*dry*dty*teuk_hrt + & +2.*dry*dpy*teuk_hrp+dty**2*teuk_htt + & +2.*dty*dpy*teuk_htp+dpy**2*teuk_hpp + + kzz(i,j,k) = drz**2*teuk_hrr+2.*drz*dtz*teuk_hrt + & +2.*drz*dpz*teuk_hrp+dtz**2*teuk_htt + & +2.*dtz*dpz*teuk_htp+dpz**2*teuk_hpp + + kxy(i,j,k) = drx*dry*teuk_hrr+(drx*dty+dtx*dry)*teuk_hrt + & +(drx*dpy+dpx*dry)*teuk_hrp + & +dtx*dty*teuk_htt+(dtx*dpy+dpx*dty)*teuk_htp+ + & dpx*dpy*teuk_hpp + + kxz(i,j,k) = drx*drz*teuk_hrr+(drx*dtz+dtx*drz)*teuk_hrt + & +(drx*dpz+dpx*drz)*teuk_hrp + & +dtx*dtz*teuk_htt+(dtx*dpz+dpx*dtz)*teuk_htp+ + & dpx*dpz*teuk_hpp + + kyz(i,j,k) = dry*drz*teuk_hrr+(dry*dtz+dty*drz)*teuk_hrt + & +(dry*dpz+dpy*drz)*teuk_hrp + & +dty*dtz*teuk_htt+(dty*dpz+dpy*dtz)*teuk_htp+ + & dpy*dpz*teuk_hpp + + + + enddo + enddo + enddo + +c initialize the conformal factor + if (use_conformal == 1) then + conformal_state = CONFORMAL_METRIC + do k=1,cctk_lsh(3) + do j=1,cctk_lsh(2) + do i=1,cctk_lsh(1) + + psi(i,j,k) = 1d0 + psix(i,j,k) = 0d0 + psiy(i,j,k) = 0d0 + psiz(i,j,k) = 0d0 + 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 + + end do + end do + end do + else + conformal_state = NOCONFORMAL_METRIC + end if + + + return + end + + + -- cgit v1.2.3