program orbit c relativistic Kerr orbit program using Runge-Kutta integration c Stephen C. Bell 06/04/96 c The program uses interactive input from the keyboard. c Answering y to the question: "Use default variable values?" c will cause the program to use the following default values: c c facmu = number of solar masses = 10.0 c ecc = beginning eccentricity = 0.5 c apans = begin at apogee? = 'y' c fracm = apogee distance (units of m) = 25.0 c rinc = beginning inclination (degs) = 45.0 c angmocen = angular momentum constant = 0.5 c tottprop = total proper time for integration (secs) = 0.04 c deltprop = integration step size (secs) = 1.0d-5 c c These values are used to generate the example orbit in the c paper, "Simulating Relativistic Orbits About a Black Hole." c Once the program is running correctly with the default c values, answer y to the "Output (t,x,y,z) data?" question, and the c program will write the proper time and the orbiting bodies (x,y,z) c coordinates (in units of m) to a file. The data can be plotted separately. c implicit double precision (a-h,o-z) dimension pos(3),vel(3) dimension use8x1(8) character fn*40,outans*1,defans*1 character apans*1 real*8 mu,m common /relstuff/ c,mu,facmu,m common /newans/ outans common /nu/ sumang,onerad common /time/ sumtau,timint common /hval/ deltprop common /am/ angmocen common /results/ use8x1 c initialize summing variables sumang=0.0d0 sumtau=0.0d0 c speed of light, km/sec from Fund. of Phys, p. 733 c=2.99792459d5 c one rad in degrees onerad=57.29577951d0 write (*,*) write (*,*) 'Use default variable values? (Entering y ' write (*,*) 'causes program to generate example orbit):' read (*,1) defans 1 format (a) if (defans .eq. 'y' .or. defans .eq. 'Y') then write (*,*) write (*,*) 'Using example defaults. See comments at' write (*,*) 'beginning of source code for explanation.' write (*,*) facmu = 10.0 c 1.3271544d11 = mu for sun in metric units, from B.M.W. c mu has been declared real*8 mu=1.3271544d11*facmu c compute 1/2 Schw. radius. m is real*8 m=mu/c**2 ecc = 0.5 apans = 'y' fracm = 25.0 pos(1)=fracm*m pos(2)=0.0d0 pos(3)=0.0d0 rinc = 45.0 rinc=rinc/onerad angmocen = 0.5 tottprop = 0.04 deltprop = 1.0d-5 goto 7000 endif write (*,*) write (*,*) 'Enter no. of solar masses', &' of the black hole: ' read (*,*) facmu c 1.3271544d11 = mu for sun in metric units, from B.M.W. mu=1.3271544d11*facmu c compute 1/2 Schw. radius. m is real*8 m=mu/c**2 write (*,*) write (*,*) 'Enter eccentricity (0 - 1): ' read (*,*) ecc write (*,*) write (*,*) 'Begin at apogee? ' read (*,1) apans write (*,*) if (apans .eq. 'y' .or. apans .eq. 'Y') then apans='y' write (*,*) 'Enter apogee distance as an m-multiple ' write (*,*) '(m = ',m,' km):' else apans='n' write (*,*) 'Enter perigee distance as an m-multiple ' write (*,*) '(m = ',m,' km):' endif read (*,*) fracm pos(1)=fracm*m write (*,*) write (*,*) 'Enter beginning inclination (degs):' read (*,*) rinc write (*,*) rinc=rinc/onerad write (*,*) 'Enter p, the angular momentum constant ' write (*,*) 'of the rotating black hole (0 - 1):' read (*,*) angmocen pos(2)=0.0d0 pos(3)=0.0d0 write (*,*) write (*,*) 'Enter total proper time for ', & 'Runge-Kutta integration (secs):' read (*,*) tottprop write (*,*) write (*,*) 'Enter proper time step size: ' read (*,*) deltprop 7000 write (*,*) write (*,*)'Output (t,x,y,z) data? ' read (*,1) outans if (outans.eq.'y'.or.outans.eq.'Y') then outans='y' write (*,*) write (*,*)'Enter output file name: ' read (*,1) fn open (2,file=fn,status='unknown') endif if (apans.eq.'y') then rectum=pos(1)*(1.0d0-ecc) else rectum=pos(1)*(1.0d0+ecc) endif h=sqrt(rectum*mu) vmag=h/pos(1) c write out rectangulars write (*,*) write (*,*) 'Beginning rectangular (x,y,z) position:' write (*,*) ' x = ',pos(1),' km = ',pos(1)/m,' m' write (*,*) ' y = ',pos(2),' km = ',pos(2)/m,' m' write (*,*) ' z = ',pos(3),' km = ',pos(3)/m,' m' vel(1)=0.0d0 vel(2)=vmag*cos(rinc) vel(3)=vmag*sin(rinc) write (*,*) write (*,*) 'Beginning rectangular (vx,vy,vz) velocity:' write (*,*) ' vx = ',vel(1),' km/sec = ', & vel(1)/m,' m/sec' write (*,*)' vy = ',vel(2),' km/sec = ', & vel(2)/m,' m/sec' write (*,*)' vz = ',vel(3),' km/sec = ', & vel(3)/m,' m/sec' write (*,*)' velocity magnitude = ',vmag,' km/sec' write (*,*)' = ',vmag/m,' m/sec' write (*,5970) 100*vmag/c 5970 format( ' = ',f6.2,' %c') c call relativistic state propagator c rectangular state passed in and returned via arguments c polar state returned via common for use8x1 call relprop(pos,vel,tottprop,pos,vel) write (*,*) write (*,*)'Runge-Kutta integration completed' write (*,*) write (*,*)'Hit return to see the results' read (*,*) write (*,*) write (*,*)'Integrated coordinate time = ',timint,' sec' write (*,*)'Total proper time = ',sumtau,' sec' write (*,*)' Difference = ',timint-sumtau,' sec' write (*,*) write (*,*)'Final proper time polar state:' r=use8x1(1) theta=use8x1(2) phi=use8x1(3) tcoor=use8x1(4) dr_dtau=use8x1(5) dtheta_dtau=use8x1(6) dphi_dtau=use8x1(7) dt_dtau=use8x1(8) write (*,*)'r = ',r/m,' m' write (*,*)'theta = ',theta,' rad' write (*,*)'phi = ',phi,' rad' write (*,*)'coordinate t = ',tcoor,' sec' write (*,*)'dr/dtau = ',dr_dtau,' m/sec' write (*,*)'dtheta/dtau = ',dtheta_dtau,' rad/sec' write (*,*)'dphi/dtau = ',dphi_dtau,' rad/sec' write (*,*)'dt/dtau = ',dt_dtau,' sec/sec' stop end subroutine relprop(pos,vel,tottprop,posout,velout) c state (pos,vel) coming in and going out is rectangular, BUT c integration is performed in polar coordinates implicit double precision (a-h,o-z) dimension posin(3),posout(3),velout(3) dimension dydx(8),stateout(8),pos(3),vel(3) real*8 mu,m,use8x1(8),vec4x1(4),gmet(4,4),dxvec(4) character outans*1 common /relstuff/ c,mu,facmu,m common /nu/ sumang,onerad common /time/ sumtau,timint c dtau = proper time step size common /hval/ dtau common /newans/ outans common /am/ angmocen common /results/ use8x1 c initialize variables for keeping track of orbit no. stoiorb=1 iorb=0 c convert to geometrized units call scavec(pos,1.0d0/c,pos) call scavec(vel,1.0d0/c,vel) m=m/c mu=m angmocen=angmocen*m c=1.0d0 c compute initial Kerr dt/dtau & coordinate time polars r=rnorm(pos) dr_dt=dot(pos,vel)/r phi=atan2(pos(2),pos(1)) dphi_dt=vel(2)/pos(1) dphi_dt=dphi_dt-pos(2)*vel(1)/pos(1)**2 dphi_dt=dphi_dt*(pos(1)**2/(pos(1)**2+pos(2)**2)) theta=acos(pos(3)/r) dtheta_dt=pos(3)/r dtheta_dt=dtheta_dt*dr_dt dtheta_dt=dtheta_dt-vel(3) dtheta_dt=dtheta_dt/sqrt(r**2-pos(3)**2) c set up Kerr metric tensor do 4953 i=1,4 do 4953 j=1,4 4953 gmet(i,j)=0.0d0 gmet(1,1)=g_r(r,theta) gmet(2,2)=g_theta(r,theta) gmet(3,3)=g_phi(r,theta) gmet(4,4)=g_t(r,theta) gmet(3,4)=g_phi_t(r,theta) gmet(4,3)=gmet(3,4) dxvec(1)=dr_dt dxvec(2)=dtheta_dt dxvec(3)=dphi_dt dxvec(4)=1.0d0 call matvec(vec4x1,gmet,dxvec,4) c compute Kerr dt/dtau dt_dtau=1.0d0/sqrt(dot4(dxvec,vec4x1)) c use dt/dtau to initialize proper time 1st derivatives dr_dtau=dr_dt*dt_dtau dphi_dtau=dphi_dt*dt_dtau dtheta_dtau=dtheta_dt*dt_dtau c initialize time summer flttime=0.0d0 iexit=0 c initialize variables for keeping track of min and max distances c so far computed rmin=r rmax=r c write data, if requested. This writes the very first pos at t = 0 if (outans .eq. 'y') & write (2,393) 0.0d0,pos(1)/m,pos(2)/m,pos(3)/m 393 format(5(1pe15.7)) c put polar state vector into 8x1 array for integration use8x1(1)=r use8x1(2)=theta use8x1(3)=phi use8x1(4)=0.0d0 use8x1(5)=dr_dtau use8x1(6)=dtheta_dtau use8x1(7)=dphi_dtau use8x1(8)=dt_dtau write (*,*) write (*,*)'Hit return for a listing of corresponding ' write (*,*)'beginning proper time polar state ' read (*,*) write (*,*) write (*,*) 'r = ',r/m,' m' write (*,*) 'theta = ',theta,' rad' write (*,*) 'phi = ',phi,' rad' write (*,*) 'coordinate t = ',0.0d0,' sec' write (*,*) 'dr/dtau = ',dr_dtau,' m/sec' write (*,*) 'dtheta/dtau = ',dtheta_dtau,' rad/sec' write (*,*) 'dphi/dtau = ',dphi_dtau,' rad/sec' write (*,*) 'dt/dtau = ',dt_dtau,' sec/sec' write (*,*) write (*,*) 'Event horizon radius = ', & (m+sqrt(m**2-angmocen**2))/m write (*,*) 'Ergosphere radius = ', & (m+sqrt(m**2-angmocen**2*cos(theta)**2))/m c write (*,*) write (*,*)'Hit return to begin the simulation' read (*,*) c initialize phi info for orbit no. tracking oldphi=phi*onerad if (oldphi.lt.0.0d0) oldphi=360.0d0+oldphi c begin integration 100 flttime=flttime+dtau if (abs(flttime).gt.abs(tottprop)) then iexit=1 dtau=abs(dtau)-(abs(flttime)-abs(tottprop)) dtau=sign(dtau,tottprop) endif c get derivatives before Runge-Kutta call c proper time polar state use8x1 is passed in c proper time polar accelerations dydx are returned c dummy time is passed in to maintain argument structure of Num. Recip. call derivs(time,use8x1,dydx) c store rectangular position for keeping track of which orbit c we are on call scavec(pos,1.0d0,posin) c call rk4 with rel state c dummy time is used to maintain argument structure of Num. Recip. call rk4(use8x1,dydx,8,time,dtau,use8x1) timint=use8x1(4) c compute rectangulars for writing output & orbit no. tracking theta=use8x1(2) phi=use8x1(3) dr_dtau=use8x1(5) dtheta_dtau=use8x1(6) dphi_dtau=use8x1(7) stateout(1)=r*sin(theta)*cos(phi) stateout(2)=r*sin(theta)*sin(phi) stateout(3)=r*cos(theta) stateout(4)=dr_dtau*sin(theta)*cos(phi) & +r*(cos(phi)*cos(theta)*dtheta_dtau & -sin(theta)*sin(phi)*dphi_dtau) stateout(5)=dr_dtau*sin(theta)*sin(phi) & +r*(sin(theta)*cos(phi)*dphi_dtau & +sin(phi)*cos(theta)*dtheta_dtau) stateout(6)=dr_dtau*cos(theta)-r*dtheta_dtau*sin(theta) c sum up proper time sumtau=sumtau+dtau c keep track of what orbit we're on posout(1)=stateout(1) posout(2)=stateout(2) posout(3)=stateout(3) phi=phi*onerad if (phi.lt.0.0d0) phi=360.0d0+phi delphi=phi-oldphi if (phi.lt.oldphi) delphi=phi+360.0d0-oldphi c compute delphi using posin and posout useval=dot(posin,posout)/rnorm(posin)/rnorm(posout) delphi=acos(useval)*onerad sumang=sumang+delphi oldphi=phi call scavec(posout,1.0d0,posin) iorb=int(sumang/360.0d0) c a mod(float(iorb),1.0) is used to write screen info every new orbit c change 1.0 to a 5.0 (or a 10.0, etc) to write c every 5th orbit (or 10th orbit, etc) if (mod(float(iorb),1.0).eq.0.0.and.stoiorb.ne.iorb) then write (*,*) if (iorb+1.eq.1) then write (*,*)'on orbit',iorb+1 write (*,*)'coordinate time = ',0.0d0,' sec' write (*,*)'proper time = ',0.0d0,' sec' write (*,*)'minimum distance so far computed = ',rmax/m,' m' write (*,*)'maximum distance so far computed = ',rmax/m,' m' else write (*,*)'on orbit',iorb+1 write (*,*)'coordinate time = ',use8x1(4),' sec' write (*,*)'proper time = ',flttime,' sec' write (*,*)'minimum distance so far computed = ',rmin/m,' m' write (*,*)'maximum distance so far computed = ',rmax/m,' m' endif stoiorb=iorb endif c keep track of max and min distances so far computed r=use8x1(1) if (r.le.rmin) rmin=r if (r.ge.rmax) rmax=r c write pos, if requested if (outans.eq.'y') then write (2,393) sumtau,posout(1)/m,posout(2)/m, & posout(3)/m endif c exit if conditions are correct if (iexit.eq.1.or.flttime.eq.tottprop) goto 200 goto 100 200 velout(1)=stateout(4) velout(2)=stateout(5) velout(3)=stateout(6) return end SUBROUTINE RK4(Y,DYDX,N,X,H,YOUT) c insert implicit statement to make RK4 double precision implicit double precision (a-h,o-z) PARAMETER (NMAX=10) DIMENSION Y(N),DYDX(N),YOUT(N),YT(NMAX),DYT(NMAX),DYM(NMAX) HH=H*0.5 H6=H/6. XH=X+HH DO 11 I=1,N YT(I)=Y(I)+HH*DYDX(I) 11 CONTINUE CALL DERIVS(XH,YT,DYT) DO 12 I=1,N YT(I)=Y(I)+HH*DYT(I) 12 CONTINUE CALL DERIVS(XH,YT,DYM) DO 13 I=1,N YT(I)=Y(I)+H*DYM(I) DYM(I)=DYT(I)+DYM(I) 13 CONTINUE CALL DERIVS(X+H,YT,DYT) DO 14 I=1,N YOUT(I)=Y(I)+H6*(DYDX(I)+DYT(I)+2.*DYM(I)) 14 CONTINUE RETURN END c general relativistic Kerr acceleration subroutine subroutine derivs(time,statevec,dydx) implicit double precision (a-h,o-z) double precision mu,m,mat4x4(4,4) dimension statevec(8),dydx(8) dimension gmet(4,4),vec4x1(4),dxvec(4),gmetinv(4,4) character outans*1 common /relstuff/ c,mu,facmu,m common /newans/ outans common /nu/ sumang,onerad common /am/ angmocen c following line used to prevent warning during compilation dummy=time c assumes polar coordinates passed in r=statevec(1) theta=statevec(2) phi=statevec(3) dr_dtau=statevec(5) dtheta_dtau=statevec(6) dphi_dtau=statevec(7) dt_dtau=statevec(8) c 1st, populate the metric matrix (tensor) and others npass=4 do 3209 i=1,npass do 3209 j=1,npass gmet(i,j)=0.0d0 mat4x4(i,j)=0.0d0 3209 gmetinv(i,j)=gmet(i,j) gmet(1,1)=g_r(r,theta) gmetinv(1,1)=gmet(1,1) gmet(2,2)=g_theta(r,theta) gmetinv(2,2)=gmet(2,2) gmet(3,3)=g_phi(r,theta) gmetinv(3,3)=gmet(3,3) gmet(4,4)=g_t(r,theta) gmetinv(4,4)=gmet(4,4) gmet(4,3)=g_phi_t(r,theta) gmetinv(4,3)=gmet(4,3) gmet(3,4)=gmet(4,3) gmetinv(3,4)=gmet(3,4) c invert the metric matrix call matinv(npass,gmetinv,det) c compute the partials c following deltas are good to use if c numerical differentiation is used delr=1.0d-5 deltheta=1.0d-5 c NOTE: In the following, the numerical c derivative computations are denoted with the word "numerical" c and the exact derivative computations with the word "exact". c The numerical code as been commented out, and the c exact derivative computations are performed. c derivatives w.r.t. r c numerical c dgr_dr=(g_r(r+delr,theta)-g_r(r,theta))/delr c exact capsig=r**2+a**2*cos(theta)**2 delta=r**2-2.0d0*m*r+a**2 a=angmocen dgr_dr=capsig*(r-m)/delta**2 dgr_dr=dgr_dr-(r/delta) dgr_dr=dgr_dr*2.0d0/c**2 c numerical c dgphi_dr=(g_phi(r+delr,theta)-g_phi(r,theta))/delr c exact dgphi_dr=m*a**2*sin(theta)**2/capsig dgphi_dr=dgphi_dr+r-2.0d0*m*r**2*a**2*sin(theta)**2/ & capsig**2 dgphi_dr=-2.0d0*sin(theta)**2*dgphi_dr/c**2 c numerical c dgtheta_dr=(g_theta(r+delr,theta)-g_theta(r,theta))/delr c exact dgtheta_dr=-2.0d0*r/c**2 c numerical c dgt_dr=(g_t(r+delr,theta)-g_t(r,theta))/delr c exact dgt_dr=(r**2-a**2*cos(theta)**2)/capsig**2 dgt_dr=dgt_dr*2.0d0*m c numerical c dgphit_dr=(g_phi_t(r+delr,theta)-g_phi_t(r,theta))/delr c exact use1=-4.0d0*a*m*r**2*sin(theta)**2/capsig**2 use2=2.0d0*a*m*sin(theta)**2/capsig dgphit_dr=(use1+use2)/c**2 c derivatives w.r.t. theta c numerical c dgr_dtheta=(g_r(r,theta+deltheta)- c & g_r(r,theta))/deltheta c exact dgr_dtheta=2.0d0*a**2*sin(theta)*cos(theta)/c**2/delta c numerical c dgphi_dtheta=(g_phi(r,theta+deltheta)- c & g_phi(r,theta))/deltheta c exact dgphi_dtheta=1.0d0+a**2*sin(theta)**2/capsig dgphi_dtheta=dgphi_dtheta*4.0d0*m*r*a**2*sin(theta)**3* & cos(theta)/capsig useval=r**2+a**2+(2.0d0*m*r*a**2*sin(theta)**2)/capsig useval=useval*2.0d0*sin(theta)*cos(theta) dgphi_dtheta=dgphi_dtheta+useval dgphi_dtheta=-dgphi_dtheta/c**2 c numerical c dgtheta_dtheta=(g_theta(r,theta+deltheta)- c & g_theta(r,theta))/deltheta c exact dgtheta_dtheta=2.0d0*a**2*cos(theta)*sin(theta)/c**2 c numerical c dgt_dtheta=(g_t(r,theta+deltheta)- c & g_t(r,theta))/deltheta c exact dgt_dtheta=-4.0d0*m*r*a**2*cos(theta)*sin(theta)/capsig**2 c numerical c dgphit_dtheta=(g_phi_t(r,theta+deltheta)- c & g_phi_t(r,theta))/deltheta c exact use1=1.0d0/capsig use2=a**2*sin(theta)**2/capsig**2 dgphit_dtheta=4.0d0*a*m*r*sin(theta)*cos(theta)*(use1+use2) dgphit_dtheta=dgphit_dtheta/c**2 c set up dxvec for lin alg computations dxvec(1)=dr_dtau dxvec(2)=dtheta_dtau dxvec(3)=dphi_dtau dxvec(4)=statevec(8) c 1st compute d2r_dtau2 vec4x1(1)=dgr_dr vec4x1(2)=dgr_dtheta c next 2 for Kerr metric vec4x1(3)=0.0d0 vec4x1(4)=0.0d0 useval=dot4(vec4x1,dxvec) useval=-useval*2.0d0*dr_dtau mat4x4(1,1)=dgr_dr mat4x4(2,2)=dgtheta_dr mat4x4(3,3)=dgphi_dr mat4x4(4,4)=dgt_dr mat4x4(3,4)=dgphit_dr mat4x4(4,3)=mat4x4(3,4) call matvec(vec4x1,mat4x4,dxvec,4) useval2=dot4(vec4x1,dxvec) useval=useval+useval2 d2r_dtau2=0.5d0*gmetinv(1,1)*useval c 2nd, compute d2theta_dtau2 vec4x1(1)=dgtheta_dr vec4x1(2)=dgtheta_dtheta c next 2 for Kerr metric vec4x1(3)=0.0d0 vec4x1(4)=0.0d0 useval=dot4(vec4x1,dxvec) useval=-useval*2.0d0*dtheta_dtau mat4x4(1,1)=dgr_dtheta mat4x4(2,2)=dgtheta_dtheta mat4x4(3,3)=dgphi_dtheta mat4x4(4,4)=dgt_dtheta mat4x4(3,4)=dgphit_dtheta mat4x4(4,3)=mat4x4(3,4) call matvec(vec4x1,mat4x4,dxvec,4) useval2=dot4(vec4x1,dxvec) useval=useval+useval2 d2theta_dtau2=0.5d0*gmetinv(2,2)*useval c 3rd, compute d2phi_dtau2 vec4x1(1)=dgphi_dr vec4x1(2)=dgphi_dtheta c next 2 for Kerr metric vec4x1(3)=0.0d0 vec4x1(4)=0.0d0 useval=dot4(vec4x1,dxvec) useval=-useval*2.0d0*dphi_dtau mat4x4(1,1)=0.0d0 mat4x4(2,2)=0.0d0 mat4x4(3,3)=0.0d0 mat4x4(4,4)=0.0d0 mat4x4(3,4)=0.0d0 mat4x4(4,3)=mat4x4(3,4) call matvec(vec4x1,mat4x4,dxvec,4) useval2=dot4(vec4x1,dxvec) useval=useval+useval2 d2phi_dtau2=0.5d0*gmetinv(3,3)*useval c 4th, compute d2t_dtau2 vec4x1(1)=dgt_dr vec4x1(2)=dgt_dtheta c next 2 for Kerr metric vec4x1(3)=0.0d0 vec4x1(4)=0.0d0 useval=dot4(vec4x1,dxvec) useval=-useval*2.0d0*dt_dtau mat4x4(1,1)=0.0d0 mat4x4(2,2)=0.0d0 mat4x4(3,3)=0.0d0 mat4x4(4,4)=0.0d0 mat4x4(3,4)=0.0d0 mat4x4(4,3)=mat4x4(3,4) call matvec(vec4x1,mat4x4,dxvec,4) useval2=dot4(vec4x1,dxvec) useval=useval+useval2 d2t_dtau2=0.5d0*gmetinv(4,4)*useval c set proper time polar coordinate accels and return dydx(1)=dr_dtau dydx(2)=dtheta_dtau dydx(3)=dphi_dtau dydx(4)=dt_dtau dydx(5)=d2r_dtau2 dydx(6)=d2theta_dtau2 dydx(7)=d2phi_dtau2 dydx(8)=d2t_dtau2 return end real*8 function dot(x,y) implicit double precision (a-h,o-z) dimension x(3),y(3) dot=0.0d0 do 1 i=1,3 1 dot=dot+x(i)*y(i) return end real*8 function dot4(x,y) implicit double precision (a-h,o-z) dimension x(4),y(4) dot4=0.0d0 do 1 i=1,4 1 dot4=dot4+x(i)*y(i) return end subroutine scavec(x,sca,y) implicit double precision (a-h,o-z) dimension x(3),y(3) do 1 i=1,3 1 y(i)=sca*x(i) return end real*8 function rnorm(x) implicit double precision (a-h,o-z) dimension x(3) rnorm=sqrt(dot(x,x)) return end SUBROUTINE MATINV(M1,A,DET) implicit double precision (a-h,o-z) DIMENSION A(4,4),IPVT(4),PVT(4),IND(4,2) DET=1. DO 1 J=1,M1 1 IPVT(J)=0 DO 10 I=1,M1 AMAX=0. DO 5 J=1,M1 IF (IPVT(J)-1) 2,5,2 2 DO 55 K=1,M1 IF (IPVT(K)-1) 3,55,200 3 IF (ABS(AMAX)-ABS(A(J,K))) 4,55,55 4 IROW=J ICOL=K AMAX=A(J,K) 55 CONTINUE 5 CONTINUE IPVT(ICOL)=IPVT(ICOL)+1 IF (IROW-ICOL) 6,8,6 6 DET=-DET DO 7 L=1,M1 SWAP=A(IROW,L) A(IROW,L)=A(ICOL,L) 7 A(ICOL,L)=SWAP 8 IND(I,1)=IROW IND(I,2)=ICOL PVT(I)=A(ICOL,ICOL) DET=DET*PVT(I) A(ICOL,ICOL)=1. DO 9 L=1,M1 9 A(ICOL,L)=A(ICOL,L)/PVT(I) DO 100 L1=1,M1 IF (L1-ICOL) 11,100,11 11 SWAP=A(L1,ICOL) A(L1,ICOL)=0. DO 12 L=1,M1 12 A(L1,L)=A(L1,L)-A(ICOL,L)*SWAP 100 CONTINUE 10 CONTINUE DO 201 I=1,M1 L=M1+1-I IF (IND(L,1)-IND(L,2)) 13,201,13 13 IROW=IND(L,1) ICOL=IND(L,2) DO 20 K=1,M1 SWAP=A(K,IROW) A(K,IROW)=A(K,ICOL) A(K,ICOL)=SWAP 20 CONTINUE 201 CONTINUE 200 RETURN END real*8 function g_r(r,theta) implicit double precision (a-h,o-z) real*8 m,mu common /am/ angmocen common /relstuff/ c,mu,facmu,m a=angmocen capsig=r**2+a**2*cos(theta)**2 delta=r**2-2.0d0*m*r+a**2 g_r=-capsig/c**2/delta return end real*8 function g_theta(r,theta) implicit double precision (a-h,o-z) real*8 m,mu common /am/ angmocen common /relstuff/ c,mu,facmu,m a=angmocen capsig=r**2+a**2*cos(theta)**2 g_theta=-capsig/c**2 return end real*8 function g_phi(r,theta) implicit double precision (a-h,o-z) real*8 m,mu common /am/ angmocen common /relstuff/ c,mu,facmu,m a=angmocen capsig=r**2+a**2*cos(theta)**2 delta=r**2-2.0d0*m*r+a**2 c like Wald g_phi=(r**2+a**2)**2 g_phi=g_phi-delta*a**2*sin(theta)**2 g_phi=g_phi/capsig g_phi=-g_phi*sin(theta)**2/c**2 return end real*8 function g_t(r,theta) implicit double precision (a-h,o-z) real*8 m,mu common /am/ angmocen common /relstuff/ c,mu,facmu,m a=angmocen capsig=r**2+a**2*cos(theta)**2 delta=r**2-2.0d0*m*r+a**2 c like Wald g_t=delta-a**2*sin(theta)**2 g_t=g_t/capsig return end real*8 function g_phi_t(r,theta) implicit double precision (a-h,o-z) real*8 m,mu common /am/ angmocen common /relstuff/ c,mu,facmu,m a=angmocen capsig=r**2+a**2*cos(theta)**2 delta=r**2-2.0d0*m*r+a**2 c like Wald g_phi_t=r**2+a**2-delta g_phi_t=g_phi_t*a*sin(theta)**2 g_phi_t=g_phi_t/capsig/c**2 return end SUBROUTINE MATVEC(RES,A,X,M) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION X(M),RES(M),A(M,M) DO 1 I=1,M RES(I)=0.D0 DO 2 J=1,M RES(I)=RES(I)+A(I,J)*X(J) 2 CONTINUE 1 CONTINUE RETURN END