C @@ C @file Gowdy.F77 C @date December 2002 C @author Denis Pollney C @desc C Cosmological Gowdy metric for a polarized wave in an C expanding universe. See C "Stable 3-level leapfrog integration in numerical relativity" C New, K, Watt, K, Misner, C and Centrella, J, PRD 58, 064022. C @desc C @version $Header$ C @@ #include "cctk.h" #include "cctk_Parameters.h" subroutine Exact__Gowdy_wave( $ x, y, z, t, $ gdtt, gdtx, gdty, gdtz, $ gdxx, gdyy, gdzz, gdxy, gdyz, gdzx, $ gutt, gutx, guty, gutz, $ guxx, guyy, guzz, guxy, guyz, guzx, $ psi, Tmunu_flag) implicit none DECLARE_CCTK_PARAMETERS c input arguments CCTK_REAL z, t CCTK_DECLARE(CCTK_REAL, x,) CCTK_DECLARE(CCTK_REAL, y,) c output arguments CCTK_REAL gdtt, gdtx, gdty, gdtz, $ gdxx, gdyy, gdzz, gdxy, gdyz, gdzx, $ gutt, gutx, guty, gutz, $ guxx, guyy, guzz, guxy, guyz, guzx CCTK_DECLARE(CCTK_REAL, psi,) LOGICAL Tmunu_flag c local static variables logical firstcall CCTK_REAL amp CCTK_REAL PI, twoPI CCTK_REAL Bessel_J0, Bessel_J1 data firstcall /.true./ save firstcall, amp, PI, twoPI, Bessel_J0, Bessel_J1 c$omp threadprivate (firstcall, amp, PI, twoPI, Bessel_J0, Bessel_J1) c local variables CCTK_REAL Bessel_J0_t, Bessel_J1_t CCTK_REAL cosz, eP, lambda CCTK_REAL d1, d2, d3, d4, d5, d6 C This is a vacuum spacetime with no cosmological constant Tmunu_flag = .false. if (firstcall) then amp = Gowdy_wave__amplitude PI = acos(-1.d0) twoPI = 2.d0*PI call jy01a(twoPI, Bessel_J0, d1, Bessel_J1, d2, d3, d4, d5, d6) firstcall = .false. end if if (t.eq.0.d0) then call CCTK_WARN(0, "The Gowdy metric is singular for t=0. You may need to set cactus::cctk_initial_time > 0.") end if call jy01a(twoPI*t, Bessel_J0_t, d1, Bessel_J1_t, d2, d3, d4, d5, d6) cosz = cos(twoPI*z) eP = exp(Bessel_J0_t * cosz) lambda = amp * (-twoPI * t * Bessel_J0_t * Bessel_J1_t * cosz**2 + + twoPI * PI * t**2 * (Bessel_J0_t**2 + Bessel_J1_t**2) + - 0.5d0 * (twoPI**2 * (Bessel_J0**2 + Bessel_J1**2) + - twoPI * Bessel_J0 * Bessel_J1)) c c lower metric components. c gdtt = -exp(0.5d0*lambda)/sqrt(t) gdtx = 0.d0 gdty = 0.d0 gdtz = 0.d0 gdxx = t * eP gdyy = t / eP gdzz = -gdtt gdxy = 0.d0 gdyz = 0.d0 gdzx = 0.d0 c c upper metric components. c gutt = 1.d0 / gdtt gutx = 0.d0 guty = 0.d0 gutz = 0.d0 guxx = 1.d0 / gdxx guyy = 1.d0 / gdyy guzz = 1.d0 / gdzz guxy = 0.d0 guyz = 0.d0 guzx = 0.d0 return end C @@ C @routine Bessel.F77 C @date December 2002 C @author Denis Pollney C @desc C Compute bessel functions of 0th and 1st order. C This routine was downloaded from: C http://iris-lee3.ece.uiuc.edu/~jjin/routines/routines.html C C ======================================================= C Purpose: Compute Bessel functions J0(x), J1(x), Y0(x), C Y1(x), and their derivatives C Input : x --- Argument of Jn(x) & Yn(x) ( x ? 0 ) C Output: BJ0 --- J0(x) C DJ0 --- dJ0(x) C BJ1 --- J1(x) C DJ1 --- dJ1(x) C BY0 --- Y0(x) C DY0 --- dY0(x) C BY1 --- Y1(x) C DY1 --- dY1(x) C ======================================================= C C @desc C @version $Header$ C @@ SUBROUTINE JY01A(X,BJ0,DJ0,BJ1,DJ1,BY0,DY0,BY1,DY1) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION A(12),B(12),A1(12),B1(12) PI=3.141592653589793D0 RP2=0.63661977236758D0 X2=X*X IF (X.EQ.0.0D0) THEN BJ0=1.0D0 BJ1=0.0D0 DJ0=0.0D0 DJ1=0.5D0 BY0=-1.0D+300 BY1=-1.0D+300 DY0=1.0D+300 DY1=1.0D+300 RETURN ENDIF IF (X.LE.12.0D0) THEN BJ0=1.0D0 R=1.0D0 DO 5 K=1,30 R=-0.25D0*R*X2/(K*K) BJ0=BJ0+R IF (DABS(R).LT.DABS(BJ0)*1.0D-15) GO TO 10 5 CONTINUE 10 BJ1=1.0D0 R=1.0D0 DO 15 K=1,30 R=-0.25D0*R*X2/(K*(K+1.0D0)) BJ1=BJ1+R IF (DABS(R).LT.DABS(BJ1)*1.0D-15) GO TO 20 15 CONTINUE 20 BJ1=0.5D0*X*BJ1 EC=DLOG(X/2.0D0)+0.5772156649015329D0 CS0=0.0D0 W0=0.0D0 R0=1.0D0 DO 25 K=1,30 W0=W0+1.0D0/K R0=-0.25D0*R0/(K*K)*X2 R=R0*W0 CS0=CS0+R IF (DABS(R).LT.DABS(CS0)*1.0D-15) GO TO 30 25 CONTINUE 30 BY0=RP2*(EC*BJ0-CS0) CS1=1.0D0 W1=0.0D0 R1=1.0D0 DO 35 K=1,30 W1=W1+1.0D0/K R1=-0.25D0*R1/(K*(K+1))*X2 R=R1*(2.0D0*W1+1.0D0/(K+1.0D0)) CS1=CS1+R IF (DABS(R).LT.DABS(CS1)*1.0D-15) GO TO 40 35 CONTINUE 40 BY1=RP2*(EC*BJ1-1.0D0/X-0.25D0*X*CS1) ELSE DATA A/-.7031250000000000D-01,.1121520996093750D+00, & -.5725014209747314D+00,.6074042001273483D+01, & -.1100171402692467D+03,.3038090510922384D+04, & -.1188384262567832D+06,.6252951493434797D+07, & -.4259392165047669D+09,.3646840080706556D+11, & -.3833534661393944D+13,.4854014686852901D+15/ DATA B/ .7324218750000000D-01,-.2271080017089844D+00, & .1727727502584457D+01,-.2438052969955606D+02, & .5513358961220206D+03,-.1825775547429318D+05, & .8328593040162893D+06,-.5006958953198893D+08, & .3836255180230433D+10,-.3649010818849833D+12, & .4218971570284096D+14,-.5827244631566907D+16/ DATA A1/.1171875000000000D+00,-.1441955566406250D+00, & .6765925884246826D+00,-.6883914268109947D+01, & .1215978918765359D+03,-.3302272294480852D+04, & .1276412726461746D+06,-.6656367718817688D+07, & .4502786003050393D+09,-.3833857520742790D+11, & .4011838599133198D+13,-.5060568503314727D+15/ DATA B1/-.1025390625000000D+00,.2775764465332031D+00, & -.1993531733751297D+01,.2724882731126854D+02, & -.6038440767050702D+03,.1971837591223663D+05, & -.8902978767070678D+06,.5310411010968522D+08, & -.4043620325107754D+10,.3827011346598605D+12, & -.4406481417852278D+14,.6065091351222699D+16/ K0=12 IF (X.GE.35.0) K0=10 IF (X.GE.50.0) K0=8 T1=X-0.25D0*PI P0=1.0D0 Q0=-0.125D0/X DO 45 K=1,K0 P0=P0+A(K)*X**(-2*K) 45 Q0=Q0+B(K)*X**(-2*K-1) CU=DSQRT(RP2/X) BJ0=CU*(P0*DCOS(T1)-Q0*DSIN(T1)) BY0=CU*(P0*DSIN(T1)+Q0*DCOS(T1)) T2=X-0.75D0*PI P1=1.0D0 Q1=0.375D0/X DO 50 K=1,K0 P1=P1+A1(K)*X**(-2*K) 50 Q1=Q1+B1(K)*X**(-2*K-1) CU=DSQRT(RP2/X) BJ1=CU*(P1*DCOS(T2)-Q1*DSIN(T2)) BY1=CU*(P1*DSIN(T2)+Q1*DCOS(T2)) ENDIF DJ0=-BJ1 DJ1=BJ0-BJ1/X DY0=-BY1 DY1=BY0-BY1/X RETURN END