c----------------------------------------------------------------------- c linear.f c Mijan F. Huq c Center for Relativity, c University of Texas at Austin c c Program to generate linearized gravitional wave initial data. c References: c Saul A. Teukolsky, Phys Rev D, vol 26, no 4, 1982 c Kenneth Eppley, "Pure Gravitational Waves", in Sources... c c Comments: c Did not fix with L'Hopital's rule terms with x^2 + y^2 in the denominator c instead chose to off-center grid from origin of coordinate system... c----------------------------------------------------------------------- c c x, y, z, r declared as 3d arrays for convenience to convert to f90 c Pass in coordinates and r c temp arrays for storing a, b, c as arrays c subroutine linearwave(nx, ny, nz, amp, width, * x, y, z, r, g11, g12, g13, g22, g23, g33, * k11, k12, k13, k21, k22, k23, k31, k32, k33, * A, B, C) implicit none integer nx, ny, nz real*8 amp, width real*8 x(nx,ny,nz), y(nx,ny,nz), z(nx,ny,nz), r(nx,ny,nz) real*8 g11(nx,ny,nz), g12(nx,ny,nz), g13(nx,ny,nz), * g22(nx,ny,nz), g23(nx,ny,nz), g33(nx,ny,nz) real*8 k11(nx,ny,nz), k12(nx,ny,nz), k13(nx,ny,nz), * k21(nx,ny,nz), k22(nx,ny,nz), k23(nx,ny,nz), * k31(nx,ny,nz), k32(nx,ny,nz), k33(nx,ny,nz) real*8 A(nx,ny,nz), B(nx,ny,nz), C(nx,ny,nz) c Local variables c Set up A, B, C radial functions call setradial(nx, ny, nz, r, a, b, c, amp, width) c Set up metric components call metcomp(nx, ny, nz, a, b, c, x, y, z, r, * g11, g12, g13, g22, g23, g33, * k11, k12, k13, k21, k22, k23, k31, k32, k33) return end c----------------------------------------------------------------------- c Routine to set up Radial functions c c Form time symmetric solution by superposing f(t+r) - f(t-r) c subroutine setradial(nx, ny, nz, r, a, b, c, amp, wid) implicit none integer nx, ny, nz real*8 r(nx, ny, nz), a(nx,ny,nz), b(nx,ny,nz), c(nx,ny,nz), . amp, wid c Local variables integer i, j, k real*8 genf, genf1, genf2, genf3, genf4 real*8 fo, fi, fo1, fi1, fo2, fi2, fo3, fi3, fo4, fi4 do i =1 , nx do j = 1, ny do k = 1, nz c Set up generating function and its derivatives for ingoing and outgoing fi = genf(r(i,j,k), amp, wid) fo = genf(-r(i,j,k), amp, wid) fi1 = genf1(r(i,j,k), amp, wid) fo1 = genf1(-r(i,j,k), amp, wid) fi2 = genf2(r(i,j,k), amp, wid) fo2 = genf2(-r(i,j,k), amp, wid) fi3 = genf3(r(i,j,k), amp, wid) fo3 = genf3(-r(i,j,k), amp, wid) fi4 = genf4(r(i,j,k), amp, wid) fo4 = genf4(-r(i,j,k), amp, wid) c Form a a(i,j,k) = 3.0 * ( . ( fo2 - fi2 )/r(i,j,k)**3 + . ( fo1 + fi1 ) * 3.0 / r(i,j,k)**4 + . ( fo - fi ) * 3.0 / r(i,j,k)**5) c Form b b(i,j,k) = - 1.0 /r(i,j,k)**2 *( . (fo3 + fi3) + 3.0/r(i,j,k)*(fo2-fi2) + . 6.0/r(i,j,k)**2 *(fo1 + fi1) + . 6.0/r(i,j,k)**3 *(fo - fi) ) . c Form c c(i,j,k) = 0.25/r(i,j,k) *( . ( fo4 - fi4) + 2.0/r(i,j,k) *(fo3 + fi3) + . 9.0/r(i,j,k)**2 *(fo2 - fi2) + . 21.0/r(i,j,k)**3 *(fo1 + fi1) + . 21.0/r(i,j,k)**4 *(fo - fi) ) end do end do end do return end c----------------------------------------------------------------------- subroutine metcomp(nx, ny, nz, a, b, c, x, y, z, r, * g11, g12, g13, g22, g23, g33, * k11, k12, k13, k21, k22, k23, k31, k32, k33) implicit none integer nx, ny, nz real*8 g11(nx, ny, nz), g12(nx, ny, nz), g13(nx, ny, nz), . g22(nx, ny, nz), g23(nx, ny, nz), g33(nx, ny, nz) real*8 a(nx, ny, nz), b(nx, ny, nz), c(nx, ny, nz), . x(nx, ny, nz), y(nx, ny, nz), z(nx, ny, nz), . r(nx, ny, nz) real*8 k11(nx, ny, nz), k12(nx, ny, nz), k13(nx, ny, nz), . k21(nx, ny, nz), k22(nx, ny, nz), k23(nx, ny, nz), . k31(nx, ny, nz), k32(nx, ny, nz), k33(nx, ny, nz) c Local variables integer i, j, k, mm real*8 frr, frth, frph, . fthph, d1fthth, d1fphph,d2fthth,d2fphph real*8 grr, grth, grph, gthth, gthph, gphph real*8 Lrx, Lry, Lrz, Ltx, Lty, Ltz, Lpx, Lpy, Lpz c want m=0 waves mm = 0 do i = 1, nx do j = 1, ny do k = 1, nz c Get values of angular functions call angfunc(0,x(i,j,k), y(i,j,k), z(i,j,k), r(i,j,k), . frr, frth, frph, fthph, d1fthth, d1fphph,d2fthth,d2fphph) c Form spherical coordinate metric components grr = 1.0 + a(i,j,k) * frr grth = b(i,j,k) *frth * r(i,j,k) grph = b(i,j,k) *frph * sqrt(x(i,j,k)**2 + y(i,j,k)**2) gthth=( 1.0 + c(i,j,k)*d1fthth + . a(i,j,k) * d2fthth )*r(i,j,k)**2 gthph=((a(i,j,k)-2.0*c(i,j,k))*fthph)* . r(i,j,k)* sqrt(x(i,j,k)**2 + y(i,j,k)**2) gphph=(1.0 + c(i,j,k)*d1fphph + a(i,j,k)*d2fphph)* . (x(i,j,k)**2 + y(i,j,k)**2) c Form transformation matrix Lrx = x(i,j,k) / r(i,j,k) Lry = y(i,j,k) / r(i,j,k) Lrz = z(i,j,k) / r(i,j,k) Ltx = x(i,j,k) * z(i,j,k) / r(i,j,k)**2 / . sqrt(x(i,j,k)**2 + y(i,j,k)**2) Lty = y(i,j,k) * z(i,j,k) / r(i,j,k)**2 / . sqrt(x(i,j,k)**2 + y(i,j,k)**2) Ltz = -sqrt(x(i,j,k)**2 + y(i,j,k)**2) / r(i,j,k)**2 Lpx = -y(i,j,k) / (x(i,j,k)**2 + y(i,j,k)**2) Lpy = x(i,j,k) / (x(i,j,k)**2 + y(i,j,k)**2) Lpz = 0.0 c Form Cartesian metric components g11(i,j,k) = grr * Lrx**2 + 2.0*grth*Lrx*Ltx + . 2.0*grph * Lrx * Lpx + gthth * Ltx**2 + . 2.0 * gthph * Ltx * Lpx + gphph * Lpx**2 g12(i,j,k) = grr * Lrx*Lry + 2.0*grth*Lrx*Lty + . 2.0*grph * Lrx * Lpy + gthth * Ltx*Lty + . 2.0 * gthph * Ltx * Lpy + gphph * Lpx*Lpy g13(i,j,k) = grr * Lrx*Lrz + 2.0*grth*Lrx*Ltz + . 2.0*grph * Lrx * Lpz + gthth * Ltx*Ltz + . 2.0 * gthph * Ltx * Lpz + gphph * Lpx*Lpz g22(i,j,k) = grr * Lry**2 + 2.0*grth*Lry*Lty + . 2.0*grph * Lry * Lpy + gthth * Lty**2 + . 2.0 * gthph * Lty * Lpy + gphph * Lpy**2 g23(i,j,k) = grr * Lry*Lrz + 2.0*grth*Lry*Ltz + . 2.0*grph * Lry * Lpz + gthth * Lty*Ltz + . 2.0 * gthph * Lty * Lpz + gphph * Lpy*Lpz g33(i,j,k) = grr * Lrz**2 + 2.0*grth*Lrz*Ltz + . 2.0*grph * Lrz * Lpz + gthth * Ltz**2 + . 2.0 * gthph * Ltz * Lpz + gphph * Lpz**2 end do end do end do do i = 1, nx do j = 1, ny do k = 1, nz k11(i,j,k) = 0.0 k12(i,j,k) = 0.0 k13(i,j,k) = 0.0 k21(i,j,k) = 0.0 k22(i,j,k) = 0.0 k23(i,j,k) = 0.0 k31(i,j,k) = 0.0 k32(i,j,k) = 0.0 k33(i,j,k) = 0.0 end do end do end do return end c----------------------------------------------------------------------- c Routine to calculate angular functions subroutine angfunc(m,x, y, z, r, . frr, frth, frph, fthph, d1fthth, d1fphph,d2fthth,d2fphph) implicit none integer m real*8 x, y, z, r, . frr, frth, frph, fthph, d1fthth, d1fphph, d2fthth, d2fphph c Local variables real*8 sint, cost, sinp, cosp sint = sqrt((x**2 + y**2)/(x**2 + y**2 + z**2)) cost = z / sqrt(x**2 + y**2 + z**2) sinp = y/sqrt(x**2 + y**2) cosp = x/sqrt(x**2 + y**2) frr = 2.0 - 3.0*sint**2 frth= -3.0 * sint * cost frph = 0.0 fthph= 0.0 d1fthth= 3.0* sint**2 d1fphph= -d1fthth d2fthth=-1.0 d2fphph= 3.0*sint**2 - 1.0 return end c----------------------------------------------------------------------- c Routine to calculate generating function real*8 function genf(x,a,w) implicit none real*8 x, a, w c genf = A*x*exp(-x**2/w**2) genf = -2*A*x/w**2*exp(-x**2/w**2) c genf = A*exp(-x**2/w**2) c write(98,*)x,A,w,genf,A*exp(-x**2/w**2) return end c----------------------------------------------------------------------- c Routine to calculate first derivative of generating function real*8 function genf1(x,a,w) implicit none real*8 x, a, w c genf1 = (w**2 - 2.0 * x**2)*A*exp(-x**2/w**2)/w**2 genf1 = 2*A*exp(-x**2/w**2)*(-w**2+2*x**2)/w**4 c genf1 = -2.0 * A * x/w**2* exp(-x**2/w**2) return end c----------------------------------------------------------------------- c Routine to calculate second derivative of generating function real*8 function genf2(x,a,w) implicit none real*8 x, a, w c genf2 = 2.0 * A* x * exp(-x**2 / w**2)*(2.0*x**2 - 3.0*w**2)/ c . w**4 genf2 = -4*A*x*exp(-x**2/w**2)*(-3*w**2+2*x**2)/w**6 c genf2 = 2.0 * A* exp(-x**2/w**2)* (-w**2+2.0*x**2)/w**4 return end c----------------------------------------------------------------------- c Routine to calculate third derivative of generating function real*8 function genf3(x,a,w) implicit none real*8 x, a, w c genf3 = -2.0 * A * exp(-x**2 / w**2)*(3.0*w**4 - c . 12.0*x**2*w**2 + 4.0*x**4)/w**6 genf3 = 4*A*exp(-x**2/w**2)*(3*w**4-12*x**2*w**2+4*x**4)/w**8 c genf3 = -4.0* A * x* exp(-x**2/w**2)*(-3.0*w**2+2.0*x**2)/w**6 return end c----------------------------------------------------------------------- c Routine to calculate fourth derivative of generating function real*8 function genf4(x,a,w) implicit none real*8 x, a, w c genf4 = 4.0*A*x*exp(-x**2 / w**2)*(15.0*w**4 - c . 20.0*x**2 *w**2 + 4.0*x**4)/w**8 genf4 = -8*A*x*exp(-x**2/w**2)*(15*w**4-20*x**2*w**2+4*x**4) . /w**10 c genf4 = 4.0 * A* exp(-x**2/w**2)* c . (3.0*w**4-12.0*x**2*w**2+4.0*x**4)/w**8 return end c-----------------------------------------------------------------------