aboutsummaryrefslogtreecommitdiff
path: root/src/metrics/Gowdy_wave.F
diff options
context:
space:
mode:
Diffstat (limited to 'src/metrics/Gowdy_wave.F')
-rw-r--r--src/metrics/Gowdy_wave.F246
1 files changed, 246 insertions, 0 deletions
diff --git a/src/metrics/Gowdy_wave.F b/src/metrics/Gowdy_wave.F
new file mode 100644
index 0000000..580841c
--- /dev/null
+++ b/src/metrics/Gowdy_wave.F
@@ -0,0 +1,246 @@
+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
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-17 22:31:38 +0000