$ 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"
subroutine IDLinearWaves_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
CCTK_REAL plus,minus,plusp,minusp,ain,aout
INTEGER i,j,k
INTEGER CCTK_Equals
pi = 3.14159265358979d0
c Wave characteristics
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 = sqrt(kx*kx+ky*ky+kz*kz)
kxp = 2*pi*sin(the)*cos(phi)/wavep
kyp = 2*pi*sin(the)*sin(phi)/wavep
kzp = 2*pi*cos(the)/wavep
wp = sqrt(kxp*kxp+kyp*kyp+kzp*kzp)
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
c Check if we should create and store conformal factor stuff */
if (CCTK_EQUALS(metric_type, "static conformal")) then
conformal_state = 1
if(CCTK_EQUALS(conformal_storage,"factor+derivs")) then
conformal_state = 2
else if(CCTK_EQUALS(conformal_storage,"factor+derivs+2nd derivs")) then
conformal_state = 3
end if
do k=1,cctk_lsh(3)
do j=1,cctk_lsh(2)
do i=1,cctk_lsh(1)
psi(i,j,k) = 1d0
if(conformal_state .gt. 1) then
psix(i,j,k) = 0d0
psiy(i,j,k) = 0d0
psiz(i,j,k) = 0d0
endif
if(conformal_state .gt. 2) then
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
endif
end do
end do
end do
end if
return
end