c================================================================
c Physics 555 G2 Numerial Relativity Spring 2005
c Code for main_template.f
c
c Integrates the equations for geodesic motion of a massive
c test-particle in 3D (theta0=Pi/2) of
c the Schwarzschild space-time.
c================================================================
c by Ben
c================================================================
program geodesic3d
implicit none
integer iargc
real*8 r8arg
c Number of equations
integer neq
parameter ( neq = 8 )
c================================================================
c Arguments for LSODA
c================================================================
c real*8 cg(4,4,4), gg(4,4), t
external fcn, jac
real*8 y(neq), yprime(neq)
real*8 taubgn, tauend
integer itol
integer itask, istate, iopt
integer lrw
parameter ( lrw = 22 + neq*16 + 100 )
real*8 rwork(lrw)
integer liw
parameter ( liw = 20 + neq +100 )
integer iwork(liw)
integer jt
real*8 rtol
parameter ( rtol = 1.0d-10 )
real*8 atol
parameter ( atol = 1.0d-10 )
c parameter ( atol = 1.0d-6 )
c================================================================
c Other locals
c================================================================
integer i
real*8 tdot0
c================================================================
c Command line arguments.
c
c <L2> is the square of the particle's proyection of
c the angular momentum vector.
c <R0> Initial Radius
c <Rdot0> is the initial time derivative (r. proper time)
c of the particle radius.
c <dtau> is the proper time interval between each output.
c <tau_f> is the final proper time of the integration.
c================================================================
real*8 R0, Rdot0,
& dtau, tau_f
real*8 phidot0, thetadot0, theta0, phi0, t
integer n
real*8 r, rdot, theta, phi, tdot
real*8 phidot,thetadot
real*8 gup(4,4)
c================================================================
c Maximum problem size (50 000)
c================================================================
integer maxnx
parameter ( maxnx = 50000 )
c================================================================
c Argument parsing
c================================================================
if( iargc() .ne. 6 ) go to 900
R0 = r8arg(1,-1.0d0)
if( R0 .lt. 0.0d0 ) go to 900
if( R0 .eq. 1.0d0 ) go to 700
theta0 = r8arg(2,-1.0d0)
phi0 = r8arg(3,-1.0d0)
Rdot0 = r8arg(4,-1.0d0)
thetadot0 = r8arg(5,-1.0d0)
phidot0 = r8arg(6,-1.0d0)
c=====================================================================
c Apply command line information
c=====================================================================
call metric(0.0,R0,theta0,phi0,gup)
tdot0=dsqrt((1.0d0+gup(2,2)*Rdot0**2+
& gup(3,3)*thetadot0**2+gup(4,4)
& *phidot0**2)/(-gup(1,1)))
c=====================================================================
c Set initial conditions using Geometric
c units G = c = 1, with r rescaled so
c that 2M = 1.
c=====================================================================
y(1) = 0.0d0
y(2) = tdot0
y(3) = R0
y(4) = Rdot0
y(5) = phi0
y(6) = phidot0
y(7) = theta0
y(8) = thetadot0
c=====================================================================
c Set LSODA parameters.
c=====================================================================
itol = 1
itask = 1
iopt = 0
jt = 2
c======================================================================
c Integrate to the next output time. (Note that LSODA sets
c taubgn = tauend we need not update taubgn.
c=====================================================================
istate = 1
call lsoda(fcn,neq,y,0.0,1.0,
& itol,rtol,atol,itask,
& istate,iopt,rwork,lrw,iwork,liw,jac,jt)
c=====================================================================
c Check for LSODA error and exit if one occured.
c=====================================================================
if( istate .lt. 0 ) then
write(0,1500) istate, i, n,
& taubgn, tauend
1500 format(/' geodesicp: Error return istate = ',i2,
& ' from integrator LSODA. '/
& ' At time step ',i5,' of ',i5/
& ' Interval ',1pe11.3,' .. ',
& 1pe11.3/)
go to 950
end if
c=====================================================================
c Output time and particle coordinates
c=====================================================================
tdot = y(2)
r = y(3)
rdot =y(4)
phi = y(5)
phidot = y(6)
theta = y(7)
thetadot = y(8)
write(*,*) r*dcos(phi)*dsin(theta),
& r*dsin(phi)*dsin(theta), r*dcos(theta)
write(*,7) r,theta,phi,rdot,thetadot,phidot
7 format(f30.17,1x,f30.17,1x,f30.17,1x,f30.17,
& 1x,f30.17,1x,f30.17)
write(10,*) 'Integration of Geodesic Equations Successfull!'
stop
700 continue
write(0,*) 'Cannot start on the event horizon '
stop
800 continue
write(0,*) 'Ratio of tau_f to dtau is too large.' //
& 'Internal storage space exceeded. '
stop
900 continue
write(0,*) 'usage: orbit <R0> <Rdot0> <theta0> ' //
& '<dtheta0> <phi0> <phidot0>'
write(0,*) R0, theta0, phi0
stop
950 continue
write(0,*) 'geodesicp: Error return from LSODA'
stop
stop
end
c=====================================================================
c Routine that will be called by LSODA. It defines the
c equations of motion for the exterior Schwarzschild
c geometry. Our vector y will contain t,tdot,r,rdot,phi,phidot
c=====================================================================
subroutine fcn(neq,tau,y,yprime)
implicit none
integer neq
real*8 tau, y(neq), yprime(neq)
real*8 cg(4,4,4), gg(4,4), pyprime(neq)
real*8 t, rdot, phidot, u, thetadot, r,
& tdot, v, w,test, theta, phi, z
c=====================================================================
c Below we define the elements of our 1-D array (vector) y(i)
c where i = 1 , neq. In our case, neq (the number of
c equations in our subroutine) is 8.
c====================================================================
t = y(1)
tdot = y(2)
r = y(3)
rdot = y(4)
phi = y(5)
phidot = y(6)
theta = y(7)
thetadot = y(8)
u = tdot
v = rdot
z = phidot
w = thetadot
c=====================================================================
c Below we define the elements of another 1-D array, ydot. These
c are our equations of motion.
c=====================================================================
call cg0 (t,r,theta,phi,cg)
yprime(1)= u
yprime(2)= -(cg(1,1,1)*u**2.0d0 + cg(1,2,2)*v**2.0d0
& +cg(1,3,3)*w**2.0d0
& +cg(1,4,4)*z**2.0d0 + 2*cg(1,1,2)*u*v+2*cg(1,1,3)*u*w
& +2*cg(1,1,4)*u*z + 2*cg(1,2,3)*v*w+2*cg(1,2,4)*v*z
& +2*cg(1,3,4)*w*z)
yprime(3)= v
yprime(4)= -(cg(2,1,1)*u**2.0d0 + cg(2,2,2)*v**2.0d0
& +cg(2,3,3)*w**2.0d0
& +cg(2,4,4)*z**2.0d0 + 2.0d0*cg(2,1,2)*u*v+2.0d0*cg(2,1,3)*u*w
& +2.0d0*cg(2,1,4)*u*z + 2.0d0*cg(2,2,3)*v*w+2.0d0*cg(2,2,4)*v*z
& +2.0d0*cg(2,3,4)*w*z)
yprime(5)= z
yprime(8)=-(cg(3,1,1)*u**2.0d0 + cg(3,2,2)*v**2.0d0
& +cg(3,3,3)*w**2.0d0
& +cg(3,4,4)*z**2.0d0 + 2.0d0*cg(3,1,2)*u*v+2.0d0*cg(3,1,3)*u*w
& +2.0d0*cg(3,1,4)*u*z + 2.0d0*cg(3,2,3)*v*w+2.0d0*cg(3,2,4)*v*z
& +2.0d0*cg(3,3,4)*w*z)
yprime(7)= w
yprime(6)= -(cg(4,1,1)*u**2.0d0 + cg(4,2,2)*v**2.0d0
& +cg(4,3,3)*w**2.0d0
& +cg(4,4,4)*z**2.0d0 + 2.0d0*cg(4,1,2)*u*v+2.0d0*cg(4,1,3)*u*w
& +2.0d0*cg(4,1,4)*u*z + 2.0d0*cg(4,2,3)*v*w+2.0d0*cg(4,2,4)*v*z
& +2.0d0*cg(4,3,4)*w*z)
return
end
c====================================================================
c Dummy Jacobian routine.
c====================================================================
subroutine jac
implicit none
return
end
subroutine metric (t, r, theta, phi, gup)
doubleprecision t
doubleprecision r
doubleprecision theta
doubleprecision phi
doubleprecision gup(4,4)
real*8 m
parameter (m=0.5)
gup(1,1) = -0.1D1 + 0.2D1 * m / r
gup(1,2) = 0.0D0
gup(1,3) = 0.0D0
gup(1,4) = 0.0D0
gup(2,1) = 0.0D0
gup(2,2) = 0.10D1 / (0.1D1 - 0.2D1 * m / r)
gup(2,3) = 0.0D0
gup(2,4) = 0.0D0
gup(3,1) = 0.0D0
gup(3,2) = 0.0D0
gup(3,3) = r ** 2
gup(3,4) = 0.0D0
gup(4,1) = 0.0D0
gup(4,2) = 0.0D0
gup(4,3) = 0.0D0
gup(4,4) = r ** 2 * sin(theta) ** 2
return
end
subroutine cg0 (t, r, theta, phi, cg)
doubleprecision t
doubleprecision r
doubleprecision theta
doubleprecision phi
doubleprecision cg(4,4,4)
real*8 m
parameter (m=0.5)
cg(1,1,1) = 0.0D0
cg(2,1,1) = -(-r + 0.2D1 * m) / r ** 3 * m
cg(3,1,1) = 0.0D0
cg(4,1,1) = 0.0D0
cg(1,1,2) = -0.10D1 / r / (-r + 0.2D1 * m) * m
cg(2,1,2) = 0.0D0
cg(3,1,2) = 0.0D0
cg(4,1,2) = 0.0D0
cg(1,3,1) = 0.0D0
cg(2,1,3) = 0.0D0
cg(3,1,3) = 0.0D0
cg(4,1,3) = 0.0D0
cg(1,4,1) = 0.0D0
cg(2,1,4) = 0.0D0
cg(3,1,4) = 0.0D0
cg(4,1,4) = 0.0D0
cg(1,2,1) = -0.10D1 / r / (-r + 0.2D1 * m) * m
cg(2,1,2) = 0.0D0
cg(3,2,1) = 0.0D0
cg(4,2,1) = 0.0D0
cg(1,2,2) = 0.0D0
cg(2,2,2) = 0.10D1 / r / (-r + 0.2D1 * m) * m
cg(3,2,2) = 0.0D0
cg(4,2,2) = 0.0D0
cg(1,2,3) = 0.0D0
cg(2,3,2) = 0.0D0
cg(3,2,3) = 0.10D1 / r
cg(4,2,3) = 0.0D0
cg(1,2,4) = 0.0D0
cg(2,2,4) = 0.0D0
cg(3,2,4) = 0.0D0
cg(4,2,4) = 0.10D1 / r
cg(1,3,1) = 0.0D0
cg(2,3,1) = 0.0D0
cg(3,1,3) = 0.0D0
cg(4,3,1) = 0.0D0
cg(1,3,2) = 0.0D0
cg(2,3,2) = 0.0D0
cg(3,2,3) = 0.10D1 / r
cg(4,3,2) = 0.0D0
cg(1,3,3) = 0.0D0
cg(2,3,3) = -r + 0.2D1 * m
cg(3,3,3) = 0.0D0
cg(4,3,3) = 0.0D0
cg(1,3,4) = 0.0D0
cg(2,3,4) = 0.0D0
cg(3,4,3) = 0.0D0
cg(4,3,4) = 0.10D1 / sin(theta) * cos(theta)
cg(1,4,1) = 0.0D0
cg(2,4,1) = 0.0D0
cg(3,4,1) = 0.0D0
cg(4,1,4) = 0.0D0
cg(1,4,2) = 0.0D0
cg(2,4,2) = 0.0D0
cg(3,4,2) = 0.0D0
cg(4,2,4) = 0.10D1 / r
cg(1,4,3) = 0.0D0
cg(2,4,3) = 0.0D0
cg(3,4,3) = 0.0D0
cg(1,4,4) = 0.0D0
cg(2,4,4) = (-r + 0.2D1 * m) * sin(theta) ** 2
cg(3,4,4) = -sin(theta) * cos(theta)
cg(4,4,4) = 0.0D0
return
end
c=================================================
c Physics 555 G2 Numerial Relativity
c Spring 2005
c Code for main_template.f
c
c Integrates the equations for geodesic
c motion of a massive
c test-particle in 3D (theta0=Pi/2) of
c the Schwarzschild space-time.
c================================================
c by Ben
c================================================
program geodesic3d
implicit none
integer iargc
real*8 r8arg
c Number of equations
integer neq
parameter ( neq = 8 )
c================================================
c Arguments for LSODA
c================================================
c real*8 cg(4,4,4), gg(4,4), t
external fcn, jac
real*8 y(neq), yprime(neq)
real*8 taubgn, tauend
integer itol
integer itask, istate, iopt
integer lrw
parameter ( lrw = 22 + neq*16 + 100 )
real*8 rwork(lrw)
integer liw
parameter ( liw = 20 + neq +100 )
integer iwork(liw)
integer jt
real*8 rtol
parameter ( rtol = 1.0d-10 )
real*8 atol
parameter ( atol = 1.0d-10 )
c parameter ( atol = 1.0d-6 )
c=========================================================
c Other locals
c========================================================
integer i
real*8 E, Ecurrent, ErrorEmax
real*8 tdot0, L2current, ErrorL2max
c=========================================================
c Command line arguments.
c
c <L2> is the square of the particle's proyection of
c the angular momentum vector.
c <R0> Initial Radius
c <Rdot0> is the initial time derivative (r. proper time)
c of the particle radius.
c <dtau> is the proper time interval between each output.
c <tau_f> is the final proper time of the integration.
c <thetadot0> initial angular vel in theta (wrt tau)
c================================================x=========
real*8 L2, R0, Rdot0,
& dtau, tau_f
real*8 phidot0, thetadot0, theta0
integer n
real*8 r, theta, phi, tdot
real*8 phidot,thetadot
real*8 gup(4,4)
c============================================================
c Maximum problem size (50 000)
c============================================================
integer maxnx
parameter ( maxnx = 50000 )
c============================================================
c Argument parsing
c============================================================
if( iargc() .ne. 6 ) go to 900
L2 = r8arg(1,-1.0d0)
if( L2 .lt. 0.0d0 ) go to 900
R0 = r8arg(2,-1.0d0)
if( R0 .lt. 0.0d0 ) go to 900
if( R0 .eq. 1.0d0 ) go to 700
Rdot0 = r8arg(3,-1.0d0)
dtau = r8arg(4,-1.0d0)
if( dtau .lt. 0.0d0 ) go to 900
tau_f = r8arg(5,-1.0d0)
if( tau_f .lt. 0.0d0 ) go to 900
thetadot0 = r8arg(6,-1.0d0)
if( thetadot0 .lt. 0.0d0 ) go to 900
c==============================================================
c Compute number of uniform time steps requested
c==============================================================
n = tau_f/dtau + 1
if( n .gt. maxnx ) go to 800
c===============================================================
c Dump number of particles, particle masses,
c initial time, and
c initial particle positions to standard output.
c===============================================================
c Output for pp3d or 2d
c
c Number of particles
write(8,*) 1
c
c Size of black hole, size of test particle
c write(8,*) 5
write(8,*) 2
c
c Write initial time, then initial black hole coords, then
c initial test particle coordinates (in 3-D cartesian)
taubgn = 0.0d0
write(8,*) taubgn
c write(8,*) 0, 0
write(8,*) R0/100.0, 0.0, 0.0
c===============================================================
c Apply command line information
c===============================================================
theta0=1.57d0
call metric(taubgn,R0,1.57d0,0.0d0,gup)
write(*,*) "gup",gup(1,1),gup(2,2),gup(3,3),gup(4,4)
c theta0=1.57
phidot0 = (dsqrt(L2))/(R0**2)
c tdot0 = dsqrt((1-1/R0)*(1 + L2/R0**2)
& + Rdot0**2)/(1 - 1/R0)
c tdot0=dsqrt((1.0d0+(1.0d0/(1.0d0-1.0d0/R0))*Rdot0**2+
c & R0**2*thetadot0**2+R0**2*sin(theta0)**2
c & *phidot0**2)/(1.0d0-1.0d0/R0))
tdot0=dsqrt((1.0d0+gup(2,2)*Rdot0**2+
& gup(3,3)*thetadot0**2+gup(4,4)
& *phidot0**2)/(-gup(1,1)))
c tdot0=tdot0-0.82842667641203
c tdot0=2.0
write(*,*) "tdot0",tdot0
E = (1-1/R0)*tdot0
c thetadot0 = 0.0
c=====================================================================
c Set initial conditions using Geometric
c units G = c = 1, with r rescaled so
c that 2M = 1.
c=====================================================================
y(1) = 0.0d0
y(2) = tdot0
y(3) = R0
y(4) = Rdot0
y(5) = 0.0d0
y(6) = phidot0
ccCheck that theta is 1.57 for the first test aware when 3d
y(7) = theta0
y(8) = thetadot0
ErrorEmax = 0
ErrorL2max = 0
c=====================================================================
c Set LSODA parameters.
c=====================================================================
itol = 1
itask = 1
iopt = 0
jt = 2
c======================================================================
c Integration loop.
c======================================================================
do i = 2 , n
c=====================================================================
c Integration
c=====================================================================
tauend = taubgn + dtau
istate = 1
call lsoda(fcn,neq,y,taubgn,tauend,
& itol,rtol,atol,itask,
& istate,iopt,rwork,lrw,iwork,liw,jac,jt)
c=====================================================================
c LSODA errors
c=====================================================================
if( istate .lt. 0 ) then
write(0,1500) istate, i, n,
& taubgn, tauend
1500 format(/' geodesicp: Error return istate = ',i2,
& ' from integrator LSODA. '/
& ' At time step ',i5,' of ',i5/
& ' Interval ',1pe11.3,' .. ',
& 1pe11.3/)
write(0,*)
write(0,*) 'Conserved Quantities of motion:'
write(0,*) 'E = '
write(0,*) E
write(0,*) 'Relative E Error = '
write(0,*) (ErrorEmax/E)
write(0,*) 'L2 = '
write(0,*) L2
write(0,*) 'Relative L2 Error = '
write(0,*) (ErrorL2max/L2)
go to 950
end if
c=====================================================================
c Output time and particle coordinates
c=====================================================================
tdot = y(2)
r = y(3)
phi = y(5)
phidot = y(6)
theta = y(7)
thetadot = y(8)
c Write proper time
write(8,*) tauend
c Write black hole coordinates
c write(8,*) 0.0d0, 0.0d0
c, 0.0d0
c Write test particle coordinates
c write(8,*) r*dcos(phi)/10.0, r*dsin(phi)/10.0, 0.0d0
write(8,*) r*dcos(phi)*dsin(theta)/100.0,
& r*dsin(phi)*dsin(theta)/100.0, r*dcos(theta)/100.0
c===================================================================
c OUTPUT xyz
c===================================================================
c Write test particle coordinates
c write(*,*) "GEODATA:", r*dcos(phi),
c & r*dsin(phi)
c write(*,*) r,phi,theta
write(14,*) r*dcos(phi)*dsin(theta),
& r*dsin(phi)*dsin(theta), r*dcos(theta)
c write event horizon
write(15,*) dcos(phi), dsin(phi)
c Proper time analysis
write(3,*) tauend,dcos(phi)
c=====================================================================
c Monitor constants of motion
c=====================================================================
Ecurrent = (1-1/r)*tdot
L2current = (r**4)*(phidot**2)
if ( abs(Ecurrent - E) .gt. ErrorEmax )
& ErrorEmax = abs(Ecurrent - E)
if ( abs(L2current - L2) .gt. ErrorL2max )
& ErrorL2max = abs(L2current - L2)
c=====================================================================
c End of integration loop.
c=====================================================================
end do
write(10,*) 'Integration of Geodesic Equations Successfull!'
write(10,*)
write(10,*) 'Our Constants of motion result:'
write(10,*) 'E = '
write(10,*) E
write(10,*) 'Relative Error in total energy = '
write(10,*) (ErrorEmax/E)
write(10,*) 'L2 = '
write(10,*) L2
write(10,*) 'Relativw Error in L2 = '
write(10,*) (ErrorL2max/L2)
stop
700 continue
write(0,*) 'Cannot start on the event horizon '
stop
800 continue
write(0,*) 'Ratio of tau_f to dtau is too large.' //
& 'Internal storage space exceeded. '
stop
900 continue
write(0,*) 'usage: orbit <L2> <R0> <Rdot0> ' //
& '<dtau> <tau_f>'
write(0,*) '<L2> is the square of the particle ' //
& 'angular momentum parameter.'
write(0,*) '<R0> is the initial radius of the particle.'
write(0,*) '<Rdot0> is the initial proper time ' //
& 'derivative of the particle radius.'
write(0,*) '<dtau> is the proper time interval ' //
& 'between each output.'
write(0,*) '<tau_f> is the final proper time of the ' //
& 'integration'
write(0,*) 'The orbit is given in geometric units ' //
& 'so m=0.5 '
stop
950 continue
write(0,*) 'geodesicp: Error return from LSODA'
stop
stop
end
c=====================================================================
c Routine that will be called by LSODA. It defines the
c equations of motion for the exterior Schwarzschild
c geometry. Our vector y will contain t,tdot,r,rdot,phi,phidot
c=====================================================================
subroutine fcn(neq,tau,y,yprime)
implicit none
integer neq
real*8 tau, y(neq), yprime(neq)
real*8 cg(4,4,4), gg(4,4), pyprime(neq)
real*8 t, rdot, phidot, u, thetadot, r,
& tdot, v, w,test, theta, phi, z
c=====================================================================
c Below we define the elements of our 1-D array (vector) y(i)
c where i = 1 , neq. In our case, neq (the number of
c equations in our subroutine) is 8.
c====================================================================
t = y(1)
tdot = y(2)
r = y(3)
rdot = y(4)
phi = y(5)
phidot = y(6)
theta = y(7)
thetadot = y(8)
c checar el valor de theta que no sea cero sino 1.57
u = tdot
v = rdot
z = phidot
w = thetadot
c w=0.0
c theta=1.57
c=====================================================================
c Below we define the elements of another 1-D array, ydot. These
c are our equations of motion.
c=====================================================================
c write(*,*) "theta",theta,thetadot
call cg0 (t,r,theta,phi,cg)
yprime(1)= u
yprime(2)= -(cg(1,1,1)*u**2.0d0 + cg(1,2,2)*v**2.0d0
& +cg(1,3,3)*w**2.0d0
& +cg(1,4,4)*z**2.0d0 + 2*cg(1,1,2)*u*v+2*cg(1,1,3)*u*w
& +2*cg(1,1,4)*u*z + 2*cg(1,2,3)*v*w+2*cg(1,2,4)*v*z
& +2*cg(1,3,4)*w*z)
yprime(3)= v
yprime(4)= -(cg(2,1,1)*u**2.0d0 + cg(2,2,2)*v**2.0d0
& +cg(2,3,3)*w**2.0d0
& +cg(2,4,4)*z**2.0d0 + 2.0d0*cg(2,1,2)*u*v+2.0d0*cg(2,1,3)*u*w
& +2.0d0*cg(2,1,4)*u*z + 2.0d0*cg(2,2,3)*v*w+2.0d0*cg(2,2,4)*v*z
& +2.0d0*cg(2,3,4)*w*z)
yprime(5)= z
c yprime(7)= 1.57
c yprime(8)= 0.0
c next eqb for theta, coord 3
yprime(8)=-(cg(3,1,1)*u**2.0d0 + cg(3,2,2)*v**2.0d0
& +cg(3,3,3)*w**2.0d0
& +cg(3,4,4)*z**2.0d0 + 2.0d0*cg(3,1,2)*u*v+2.0d0*cg(3,1,3)*u*w
& +2.0d0*cg(3,1,4)*u*z + 2.0d0*cg(3,2,3)*v*w+2.0d0*cg(3,2,4)*v*z
& +2.0d0*cg(3,3,4)*w*z)
yprime(7)= w
yprime(6)= -(cg(4,1,1)*u**2.0d0 + cg(4,2,2)*v**2.0d0
& +cg(4,3,3)*w**2.0d0
& +cg(4,4,4)*z**2.0d0 + 2.0d0*cg(4,1,2)*u*v+2.0d0*cg(4,1,3)*u*w
& +2.0d0*cg(4,1,4)*u*z + 2.0d0*cg(4,2,3)*v*w+2.0d0*cg(4,2,4)*v*z
& +2.0d0*cg(4,3,4)*w*z)
c=======================================================
c auxiliar code for testing
c======================================================
c yprime(1)= u
c test=-2.0d0*cg(1,1,2)*u*v
c pyprime(2) = tdot*rdot/(r*(1-r))
c pyprime(2) = u*v/(r*(1-r))
c yprime(3)= v
c yprime(4)=(1-r)*(tdot**2.0d0)/(2*r**3)-(rdot**2.0d0)/((2*r)*(1-r))
c & + r*(1 - 1/r)*(phidot**2.0d0)
c yprime(5)= z
cc yprime(7)=w
c yprime(6)= -2*rdot*phidot/r
cc yprime(8) = 0.0
c write(12,*) "PYS" ,yprime(2),pyprime(2),test,theta
return
end
c====================================================================
c Dummy Jacobian routine.
c====================================================================
subroutine jac
implicit none
return
end
subroutine metric (t, r, theta, phi, gup)
doubleprecision t
doubleprecision r
doubleprecision theta
doubleprecision phi
doubleprecision gup(4,4)
real*8 m
parameter (m=0.5)
gup(1,1) = -0.1D1 + 0.2D1 * m / r
gup(1,2) = 0.0D0
gup(1,3) = 0.0D0
gup(1,4) = 0.0D0
gup(2,1) = 0.0D0
gup(2,2) = 0.10D1 / (0.1D1 - 0.2D1 * m / r)
gup(2,3) = 0.0D0
gup(2,4) = 0.0D0
gup(3,1) = 0.0D0
gup(3,2) = 0.0D0
gup(3,3) = r ** 2
gup(3,4) = 0.0D0
gup(4,1) = 0.0D0
gup(4,2) = 0.0D0
gup(4,3) = 0.0D0
gup(4,4) = r ** 2 * sin(theta) ** 2
return
end
subroutine cg0 (t, r, theta, phi, cg)
doubleprecision t
doubleprecision r
doubleprecision theta
doubleprecision phi
doubleprecision cg(4,4,4)
real*8 m
parameter (m=0.5)
cg(1,1,1) = 0.0D0
cg(2,1,1) = -(-r + 0.2D1 * m) / r ** 3 * m
cg(3,1,1) = 0.0D0
cg(4,1,1) = 0.0D0
cg(1,1,2) = -0.10D1 / r / (-r + 0.2D1 * m) * m
cg(2,1,2) = 0.0D0
cg(3,1,2) = 0.0D0
cg(4,1,2) = 0.0D0
cg(1,3,1) = 0.0D0
cg(2,1,3) = 0.0D0
cg(3,1,3) = 0.0D0
cg(4,1,3) = 0.0D0
cg(1,4,1) = 0.0D0
cg(2,1,4) = 0.0D0
cg(3,1,4) = 0.0D0
cg(4,1,4) = 0.0D0
cg(1,2,1) = -0.10D1 / r / (-r + 0.2D1 * m) * m
cg(2,1,2) = 0.0D0
cg(3,2,1) = 0.0D0
cg(4,2,1) = 0.0D0
cg(1,2,2) = 0.0D0
cg(2,2,2) = 0.10D1 / r / (-r + 0.2D1 * m) * m
cg(3,2,2) = 0.0D0
cg(4,2,2) = 0.0D0
cg(1,2,3) = 0.0D0
cg(2,3,2) = 0.0D0
cg(3,2,3) = 0.10D1 / r
cg(4,2,3) = 0.0D0
cg(1,2,4) = 0.0D0
cg(2,2,4) = 0.0D0
cg(3,2,4) = 0.0D0
cg(4,2,4) = 0.10D1 / r
cg(1,3,1) = 0.0D0
cg(2,3,1) = 0.0D0
cg(3,1,3) = 0.0D0
cg(4,3,1) = 0.0D0
cg(1,3,2) = 0.0D0
cg(2,3,2) = 0.0D0
cg(3,2,3) = 0.10D1 / r
cg(4,3,2) = 0.0D0
cg(1,3,3) = 0.0D0
cg(2,3,3) = -r + 0.2D1 * m
cg(3,3,3) = 0.0D0
cg(4,3,3) = 0.0D0
cg(1,3,4) = 0.0D0
cg(2,3,4) = 0.0D0
cg(3,4,3) = 0.0D0
cg(4,3,4) = 0.10D1 / sin(theta) * cos(theta)
cg(1,4,1) = 0.0D0
cg(2,4,1) = 0.0D0
cg(3,4,1) = 0.0D0
cg(4,1,4) = 0.0D0
cg(1,4,2) = 0.0D0
cg(2,4,2) = 0.0D0
cg(3,4,2) = 0.0D0
cg(4,2,4) = 0.10D1 / r
cg(1,4,3) = 0.0D0
cg(2,4,3) = 0.0D0
cg(3,4,3) = 0.0D0
cg(1,4,4) = 0.0D0
cg(2,4,4) = (-r + 0.2D1 * m) * sin(theta) ** 2
cg(3,4,4) = -sin(theta) * cos(theta)
cg(4,4,4) = 0.0D0
return
end
====================================================
#!/bin/sh
#Master driver script Geodesic integrator
# by Ben
#form the GUI this script is triggered
#First we run maple/grii which reads the metricfromguifile
#and outputs the christof_code fortran subroutine
#rm -rf *.ps
maple << END
# VERY IMPORTANT: METRICS ARE LOADED IN THE CURRENT
# WORKING
# GEO4 FOR VER 1.0
#
######################################################
##-- Maple routine to initialize GRTensor
##
## Add to ~/.mapleinit
##
## Intialize via
##
## grtw ();
##
## from the Maple prompt.
##
#grtw := proc ()
# global libname, grOptionMetricPath, grOptionqloadPath:
#
# libname := "/usr/local/maple6/grii/lib", libname:
#
# readlib (grii):
# grtensor():
# grOptionqloadPath := "/d/bh4/home/benjamin/geo4";
# grOptionqloadPath := "/usr/local/maple6/grii/metrics";
# grOptionMetricPath:= "/d/bh4/home/benjamin";
#
#
#end:
###############################################
grtw():
qload(metricfromgui);
grcalc(Chr(dn,dn,up)):
#grdisplay((dn,dn)):
mychristof := proc(gamma::array(1..4,1..4,1..4)
,t,r,theta,phi);
for i from 1 to 4 do
for j from 1 to 4 do
for k from 1 to 4 do
gamma[i,j,k] := grcomponent( Chr(dn,
dn,up),[i,j,k])
end do;
end do;
end do;
return;
end proc;
gammac:=array(1..4,1..4,1..4):
mychristof(gammac,t,r,theta,phi);
body:=[
gammaup[1,1,1]=gammac[1,1,1],
gammaup[2,1,1]=gammac[1,1,2],
gammaup[3,1,1]=gammac[1,1,3],
gammaup[4,1,1]=gammac[1,1,4],
gammaup[1,1,2]=gammac[1,2,1],
gammaup[2,1,2]=gammac[1,2,2],
gammaup[3,1,2]=gammac[1,2,3],
gammaup[4,1,2]=gammac[1,2,4],
gammaup[1,3,1]=gammac[1,3,1],
gammaup[2,1,3]=gammac[1,3,2],
gammaup[3,1,3]=gammac[1,3,3],
gammaup[4,1,3]=gammac[1,3,4],
gammaup[1,4,1]=gammac[1,4,1],
gammaup[2,1,4]=gammac[1,4,2],
gammaup[3,1,4]=gammac[1,4,3],
gammaup[4,1,4]=gammac[1,4,4],
gammaup[1,2,1]=gammac[2,1,1],
gammaup[2,1,2]=gammac[2,1,2],
gammaup[3,2,1]=gammac[2,1,3],
gammaup[4,2,1]=gammac[2,1,4],
gammaup[1,2,2]=gammac[2,2,1],
gammaup[2,2,2]=gammac[2,2,2],
gammaup[3,2,2]=gammac[2,2,3],
gammaup[4,2,2]=gammac[2,2,4],
gammaup[1,2,3]=gammac[2,3,1],
gammaup[2,3,2]=gammac[2,3,2],
gammaup[3,2,3]=gammac[2,3,3],
gammaup[4,2,3]=gammac[2,3,4],
gammaup[1,2,4]=gammac[2,4,1],
gammaup[2,2,4]=gammac[2,4,2],
gammaup[3,2,4]=gammac[2,4,3],
gammaup[4,2,4]=gammac[2,4,4],
gammaup[1,3,1]=gammac[3,1,1],
gammaup[2,3,1]=gammac[3,1,2],
gammaup[3,1,3]=gammac[3,1,3],
gammaup[4,3,1]=gammac[3,1,4],
gammaup[1,3,2]=gammac[3,2,1],
gammaup[2,3,2]=gammac[3,2,2],
gammaup[3,2,3]=gammac[3,2,3],
gammaup[4,3,2]=gammac[3,2,4],
gammaup[1,3,3]=gammac[3,3,1],
gammaup[2,3,3]=gammac[3,3,2],
gammaup[3,3,3]=gammac[3,3,3],
gammaup[4,3,3]=gammac[3,3,4],
gammaup[1,3,4]=gammac[3,4,1],
gammaup[2,3,4]=gammac[3,4,2],
gammaup[3,4,3]=gammac[3,4,3],
gammaup[4,3,4]=gammac[3,4,4],
gammaup[1,4,1]=gammac[4,1,1],
gammaup[2,4,1]=gammac[4,1,2],
gammaup[3,4,1]=gammac[4,1,3],
gammaup[4,1,4]=gammac[4,1,4],
gammaup[1,4,2]=gammac[4,2,1],
gammaup[2,4,2]=gammac[4,2,2],
gammaup[3,4,2]=gammac[4,2,3],
gammaup[4,2,4]=gammac[4,2,4],
gammaup[1,4,3]=gammac[4,3,1],
gammaup[2,4,3]=gammac[4,3,2],
gammaup[3,4,3]=gammac[4,3,3],
#gammaup[4,3,4]=gammac[4,3,4],
gammaup[1,4,4]=gammac[4,4,1],
gammaup[2,4,4]=gammac[4,4,2],
gammaup[3,4,4]=gammac[4,4,3],
gammaup[4,4,4]=gammac[4,4,4]
];
christof:=codegen[makevoid](codegen[makeproc]
(body,[t,r,theta,phi,gammaup::array(1..4,1..4,1..4)])):
CodeGeneration[Fortran](christof,precision=double,
deducetypes=false,output='christof_code'):
#print(gg);
END
maple << END
grtw():
qload(metricfromgui);
grdisplay(g(dn,dn)):
mymetric := proc(gg::array(1..4,1..4),r,theta,phi,t);
gg[1,1]:=grcomponent( g(dn,dn),[1,1]);
gg[1,2]:=grcomponent( g(dn,dn),[1,2]);
gg[1,3]:=grcomponent( g(dn,dn),[1,3]);
gg[1,4]:=grcomponent( g(dn,dn),[1,4]);
gg[2,1]:=grcomponent( g(dn,dn),[2,1]);
gg[2,2]:=grcomponent( g(dn,dn),[2,2]);
gg[2,3]:=grcomponent( g(dn,dn),[2,3]);
gg[2,4]:=grcomponent( g(dn,dn),[2,4]);
gg[3,1]:=grcomponent( g(dn,dn),[3,1]);
gg[3,2]:=grcomponent( g(dn,dn),[3,2]);
gg[3,3]:=grcomponent( g(dn,dn),[3,3]);
gg[3,4]:=grcomponent( g(dn,dn),[3,4]);
gg[4,1]:=grcomponent( g(dn,dn),[4,1]);
gg[4,2]:=grcomponent( g(dn,dn),[4,2]);
gg[4,3]:=grcomponent( g(dn,dn),[4,3]);
gg[4,4]:=grcomponent( g(dn,dn),[4,4]);
return;
end proc;
gg:=array(1..4,1..4):
mymetric(gg,r,theta,phi,t);
body:=[
gup[1,1]=gg[1,1],
gup[1,2]=gg[1,2],
gup[1,3]=gg[1,3],
gup[1,4]=gg[1,4],
gup[2,1]=gg[2,1],
gup[2,2]=gg[2,2],
gup[2,3]=gg[2,3],
gup[2,4]=gg[2,4],
gup[3,1]=gg[3,1],
gup[3,2]=gg[3,2],
gup[3,3]=gg[3,3],
gup[3,4]=gg[3,4],
gup[4,1]=gg[4,1],
gup[4,2]=gg[4,2],
gup[4,3]=gg[4,3],
gup[4,4]=gg[4,4]
];
metric:=codegen[makevoid](codegen[makeproc]
(body,[t,r,theta,phi,gup::array(1..4,1..4)])):
CodeGeneration[Fortran](metric,precision
=double,deducetypes=false,output='metric_code'):
#print(gg);
END
#then we merge the main fortran file with the subroutine
cp main_template.f main.f
sed '6a\ parameter (m=0.5)' christof_code
> christof_code2.f
sed '6a\ real*8 m' christof_code2.f
> christof_code3.f
cat christof_code3.f >> maintemp.f
sed '6a\ parameter (m=0.5)' metric_code
> metric_code2.f
sed '6a\ real*8 m' metric_code2.f
> metric_code3.f
cat metric_code3.f maintemp.f >> main.f
rm -rf christof_code*
rm -rf metric_code*
make >> geodesic.log
echo $1 $2 $3 $4 $5 $6 $#
./main $1 $2 $3 $4 $5 $6
cat fort.10 >> geodesic.log
./3d test
gv test.ps &
#make clean
cp -f main.f main.log.f
rm -rf main.f
rm -rf maintemp.f
mv -f geodesic.log geodesic.0.log
#exit
==============================================