c 
c relativistic Kerr orbit program using Runge-Kutta integration 
c Stephen C. Bell   12/1998 
c 
c 12/98 - incorporated precession computation 
c 
c The following 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    = no. of solar masses = 10.0 
c  ecc      = beginning eccentricity = 0.5 
c  apans    = begin at apogee? = 'y' 
c  fracm    = apogee distance, in 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 are the values used to generate the example orbit in the  
c paper  "Simulating Relativistic Orbits About a Black Hole." 
c Once you've got the program 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 position (in units of m) to a file (the name of which you specify). 
c These data can then be 3-D plotted. 
c 
c 
      program orbit 
      implicit double precision (a-h,o-z) 
      dimension pos(3),vel(3),begecc(3) 
      dimension use8x1(8),usevec(3),eccvec(3) 
      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,begecc 
      common /write/ iwrite,iwritetot 
 
      pi=3.141592654d0 
      iwritetot=0 
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 
 
c 
cc 
c 
      write (*,*) 
      write (*,*)'Use default variable values? (Entering y <ret>' 
      write (*,*)'will cause 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 an explanation.' 
         write (*,*) 
c 
         facmu = 10.0 
 
c 1.3271544d11 = mu for sun in metric units, from B.M.W. 
c remember, mu has been declared real*8 
         mu=1.3271544d11*facmu 
c 
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 
c 
         goto 7000 
c 
       endif 
 
      write (*,*) 
      write (*,*)'Enter the no. of solar masses', 
     &' of the black hole: ' 
      read (*,*) facmu 
 
c 1.3271544d11 = mu for sun in metric units, from B.M.W. 
c remember, mu has been declared real*8 
      mu=1.3271544d11*facmu 
c 
c compute 1/2 Schw. radius. m is real*8 
      m=mu/c**2 
c 
 
         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 
9476     write (*,*) 
         write (*,*)'Enter beginning inclination (degs):' 
         read (*,*) rinc 
         if (rinc.eq.90.0d0) then 
            write (*,*) 
            write (*,*) 
            write (*,*)'**************************************' 
            write (*,*)'Due to numerical inaccuracies for a 90' 
            write (*,*)'degree inclination, it is not allowed.' 
            write (*,*)'**************************************' 
            write (*,*) 
            write (*,*) 
            goto 9476 
         endif 
         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 
 
c 
c 
      if (apans.eq.'y') then 
         abmw=pos(1)/(1.0d0+ecc) 
      else 
         abmw=pos(1)/(1.0d0-ecc) 
      endif 
      period=2.0d0*pi*abmw**1.5d0/sqrt(mu) 
      write (*,*) 
      write (*,*) 
      write (*,*)'The Newtonian period of one orbit for the data' 
      write (*,*)'you have entered is ',period,' secs.' 
      if (rinc*onerad.lt.22.5d0) then 
         suggest=period/1000.0d0 
      else  
         suggest=period/10000.0d0 
      endif 
 
      if (fracm.le.6.0d0) then 
         suggest=period/50000.0d0 
      endif 
 
 
 
      write (*,*) 
      write (*,*)'The recommended integration step size = ', 
     &           suggest 
      
 
      write (*,*) 
      write (*,*)'Enter total proper time for ', 
     &           'Runge-Kutta integration (secs):' 
      read (*,*) tottprop 
 
 
c 
      write (*,*) 
      write (*,*)'Enter proper time step size: ' 
      read (*,*) deltprop 
       
7000  continue 
 
c Blackburn - the following section of commented out code will 
c run the Blackburn example, found in CIP, Jan/Feb 1994, p. 98. 
c Comment back in if wish to run, and recompile. If commented  
c back in, will override any interactive input. 
 
cc Begin Blackburn 
c         write (*,*) 
c         write (*,*) 
c         write (*,*)'******** RUNNING BLACKBURN EXAMPLE ********' 
c         write (*,*)'* ALL INTERACTIVE INPUT HAS BEEN OVERWRITTEN *' 
c         write (*,*) 
c         write (*,*) 
c         facmu=200000.0d0+3020.9d0 
c         mu=1.3271544d11*facmu 
c         m=mu/c**2 
c         ecc=0.5d0 
c         apans='y' 
c         fracm=40.0d0 
c         pos(1)=fracm*m 
c         pos(2)=0.0d0 
c         pos(3)=0.0d0 
c         vel(1)=0.0d0 
c         vel(2)=0.1d0*m 
c         vel(3)=-0.05d0*m 
c         angmocen=0.25d0 
c         rinc=-26.56505118d0/onerad 
c         tottprop = 10000.0d0 
c         deltprop = 0.125d0 
cc End Blackburn 
 
      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') 
         write (*,*) 
         write (*,*)'By default, the program will write every time' 
         write (*,*)'step when the position of the test particle is' 
         write (*,*)'updated.' 
         write (*,*) 
         write (*,*)'Enter the no. of time steps to skip in between' 
         write (*,*)'when the data is written (entering 1 keeps' 
         write (*,*)'the default, i.e., written at every time step): ' 
         write (*,*) 
         read (*,*) iwrite 
      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 compute unit beginning ecc vec for precession computation 
      call crsprd(pos,vel,usevec) 
      call crsprd(vel,usevec,eccvec) 
      call scavec(eccvec,1.0d0/mu,eccvec) 
      call scavec(pos,1.0d0/rnorm(pos),usevec) 
      begecc(1)=eccvec(1)-usevec(1) 
      begecc(2)=eccvec(2)-usevec(2) 
      begecc(3)=eccvec(3)-usevec(3) 
      call scavec(begecc,1.0d0/rnorm(begecc),begecc) 
 
 
 
 
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,' m/sec' 
      write (*,*)'dtheta/dtau  = ',dtheta_dtau,' rad/sec' 
      write (*,*)'dphi/dtau    = ',dphi_dtau,' rad/sec' 
      write (*,*)'dt/dtau      = ',dt_dtau,' sec/sec' 
 
c 
c final unit ecc vec for precession computation 
       pos(1)=r*sin(theta)*cos(phi) 
       pos(2)=r*sin(theta)*sin(phi) 
       pos(3)=r*cos(theta) 
       vel(1)=dr_dtau*sin(theta)*cos(phi) 
     &             +r*(cos(phi)*cos(theta)*dtheta_dtau 
     &                -sin(theta)*sin(phi)*dphi_dtau) 
       vel(2)=dr_dtau*sin(theta)*sin(phi) 
     &             +r*(sin(theta)*cos(phi)*dphi_dtau 
     &                +sin(phi)*cos(theta)*dtheta_dtau) 
       vel(3)=dr_dtau*cos(theta)-r*dtheta_dtau*sin(theta) 
      call scavec(vel,1.0d0/dt_dtau,vel) 
      call crsprd(pos,vel,usevec) 
      call crsprd(vel,usevec,eccvec) 
      call scavec(eccvec,1.0d0/mu,eccvec) 
      call scavec(pos,1.0d0/rnorm(pos),usevec) 
      eccvec(1)=eccvec(1)-usevec(1) 
      eccvec(2)=eccvec(2)-usevec(2) 
      eccvec(3)=eccvec(3)-usevec(3) 
      call scavec(eccvec,1.0d0/rnorm(eccvec),eccvec) 
      useval=dot(eccvec,begecc) 
      precess=acos(useval)*180.0d0/pi 
      precess=precess*3600.0d0 
      write (*,*)'precession = ',precess,' arcsecs = ', 
     &           precess/3600.0,' degs' 
 
 
      write (*,*) 
      write (*,*)'Enter return to terminate program...' 
      read (*,*) 
      stop 
      end 
c 
c 
c 
      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),begecc(3) 
      dimension usevec(3),eccvec(3) 
      dimension dydx(8),stateout(6),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,begecc 
      common /write/ iwrite,iwritetot 
 
      pi=3.141592654d0 
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 
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 
c initialize variables for keeping track of min and max distances 
c so far computed 
      rmin=r 
      rmax=r 
c 
c write data, if requested.  This writes the very first pos at t = 0 
       if (outans.eq.'y') then  
          write (2,393) 0.0d0,pos(1)/m,pos(2)/m,pos(3)/m,0.0d0 
       endif 
393   format (5(1pe15.7)) 
 
c 
c 
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,' 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 
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're 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) 
       dt_dtau=use8x1(8) 
       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) 
      velout(1)=stateout(4)/dt_dtau 
      velout(2)=stateout(5)/dt_dtau 
      velout(3)=stateout(6)/dt_dtau 
      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) 
c 
      iorb=int(sumang/360.0d0) 
c 
c a mod(float(iorb),1.0) is used to write screen info every new orbit 
c change the 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' 
          write (*,*)'precession = ',precess 
          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 
c compute precession 
      call crsprd(posout,velout,usevec) 
      call crsprd(velout,usevec,eccvec) 
      call scavec(eccvec,1.0d0/mu,eccvec) 
      call scavec(posout,1.0d0/rnorm(posout),usevec) 
      eccvec(1)=eccvec(1)-usevec(1) 
      eccvec(2)=eccvec(2)-usevec(2) 
      eccvec(3)=eccvec(3)-usevec(3) 
      call scavec(eccvec,1.0d0/rnorm(eccvec),eccvec) 
      useval=dot(eccvec,begecc) 
      precess=acos(useval)*180.0d0/pi 
      precess=precess*3600.0d0 
 
c write pos, if requested 
       iwritetot=iwritetot+1 
       if (outans.eq.'y'.and.iwritetot.eq.iwrite) then 
          write (2,393) sumtau,posout(1)/m,posout(2)/m, 
     &                  posout(3)/m,precess/3600.0d0 
       iwritetot=0 
       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 
c 
c 
      SUBROUTINE RK4(Y,DYDX,N,X,H,YOUT) 
c inserted 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       
c 
c general relativistic Kerr acceleration subroutine 
c 
      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 
c following line used to prevent warning during compilation 
      dummy=time 
c 
c 
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 
c 1st, populate the metric matrix (tensor) and others 
c this is Kerr specific 
      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 
c compute the partials 
c the 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 first, compute the partials of sig, delta 
      a=angmocen 
      dsig_dr=2.0d0*r 
      dsig_dtheta=-2.0d0*(a**2)*cos(theta)*sin(theta) 
      ddel_dr=2.0d0*(r-m) 
 
 
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 
      dgr_dr=capsig*ddel_dr/delta 
      dgr_dr=(dgr_dr-dsig_dr)/c**2/delta 
 
 
c numerical 
c      dgphi_dr=(g_phi(r+delr,theta)-g_phi(r,theta))/delr 
c exact 
      use=(r**2+a**2)**2-delta*a**2*(sin(theta))**2 
      use=use*dsig_dr/capsig**2 
      use=use-4.0d0*r*(r**2+a**2)/capsig 
      use=use+a**2*(sin(theta))**2*ddel_dr/capsig 
      dgphi_dr=use*(sin(theta))**2/c**2 
 
 
c numerical 
c      dgtheta_dr=(g_theta(r+delr,theta)-g_theta(r,theta))/delr 
c exact 
      dgtheta_dr=-dsig_dr/c**2 
 
 
 
c numerical 
c      dgt_dr=(g_t(r+delr,theta)-g_t(r,theta))/delr 
c exact 
      use=delta-a**2*(sin(theta))**2 
      use=use*dsig_dr/capsig 
      use=ddel_dr-use 
      dgt_dr=use/capsig 
 
 
       
c numerical 
c      dgphit_dr=(g_phi_t(r+delr,theta)-g_phi_t(r,theta))/delr 
c exact 
      use=-(r**2+a**2-delta)*dsig_dr/capsig 
      use=use+2.0d0*r-ddel_dr 
      dgphit_dr=use*a*(sin(theta))**2/c**2/capsig 
 
 
 
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=-dsig_dtheta/c**2/delta 
 
 
c numerical 
c      dgphi_dtheta=(g_phi(r,theta+deltheta)- 
c     &                g_phi(r,theta))/deltheta 
c exact 
      useval=(r**2+a**2)**2-delta*a**2*(sin(theta))**2 
      use=useval/capsig 
      use=use-delta 
      use=use*sin(theta)*dsig_dtheta 
      use=use-2.0d0*useval*cos(theta) 
      dgphi_dtheta=use*sin(theta)/c**2/capsig 
 
 
 
c numerical 
c      dgtheta_dtheta=(g_theta(r,theta+deltheta)- 
c     &                g_theta(r,theta))/deltheta 
c exact 
      dgtheta_dtheta=-dsig_dtheta/c**2 
 
 
 
c numerical 
c      dgt_dtheta=(g_t(r,theta+deltheta)- 
c     &                g_t(r,theta))/deltheta 
c exact 
      use=a**2*(sin(theta))**2-delta 
      use=use*dsig_dtheta/capsig 
      use=use-2.0d0*a**2*sin(theta)*cos(theta) 
      dgt_dtheta=use/capsig 
 
 
       
c numerical 
c      dgphit_dtheta=(g_phi_t(r,theta+deltheta)- 
c     &                g_phi_t(r,theta))/deltheta 
c exact 
      use=-sin(theta)*(r**2+a**2-delta)*dsig_dtheta/capsig 
      use=use+2.0d0*(r**2+a**2-delta)*cos(theta) 
      dgphit_dtheta=use*a*sin(theta)/c**2/capsig 
 
 
 
c set up dxvec for matrix algebra computations 
       dxvec(1)=dr_dtau 
       dxvec(2)=dtheta_dtau 
       dxvec(3)=dphi_dtau 
       dxvec(4)=statevec(8) 
 
c 1st, compute d2r_dtau2 
c do with Chris symb 
c these two do loops ensure mat4x4 set to zero matrix 
      do i = 1,4 
      do j = 1,4 
         mat4x4(i,j)=0.0d0 
      enddo 
      enddo 
      mat4x4(1,1)=0.5d0*gmetinv(1,1)*dgr_dr 
      mat4x4(1,2)=0.5d0*gmetinv(1,1)*dgr_dtheta 
      mat4x4(1,3)=0.0d0 
      mat4x4(1,4)=0.0d0 
      mat4x4(2,1)=mat4x4(1,2) 
      mat4x4(2,2)=-0.5d0*gmetinv(1,1)*dgtheta_dr 
      mat4x4(2,3)=0.0d0 
      mat4x4(2,4)=0.0d0 
      mat4x4(3,1)=0.0d0 
      mat4x4(3,2)=0.0d0 
      mat4x4(3,3)=-0.5d0*gmetinv(1,1)*dgphi_dr 
      mat4x4(3,4)=-0.5d0*gmetinv(1,1)*dgphit_dr 
      mat4x4(4,1)=0.0d0 
      mat4x4(4,2)=0.0d0 
      mat4x4(4,3)=mat4x4(3,4) 
      mat4x4(4,4)=-0.5d0*gmetinv(1,1)*dgt_dr 
      call matvec(vec4x1,mat4x4,dxvec,4) 
      useval=-dot4(vec4x1,dxvec) 
      d2r_dtau2=useval 
 
 
 
c 2nd, compute d2theta_dtau2 
c do it with Chris symb 
c these two do loops ensure mat4x4 set to zero matrix 
      do i = 1,4 
      do j = 1,4 
         mat4x4(i,j)=0.0d0 
      enddo 
      enddo 
      mat4x4(1,1)=-0.5d0*gmetinv(2,2)*dgr_dtheta 
      mat4x4(1,2)=0.5d0*gmetinv(2,2)*dgtheta_dr 
      mat4x4(1,3)=0.0d0 
      mat4x4(1,4)=0.0d0 
      mat4x4(2,1)=mat4x4(1,2) 
      mat4x4(2,2)=0.5d0*gmetinv(2,2)*dgtheta_dtheta 
      mat4x4(2,3)=0.0d0 
      mat4x4(2,4)=0.0d0 
      mat4x4(3,1)=0.0d0 
      mat4x4(3,2)=0.0d0 
      mat4x4(3,3)=-0.5d0*gmetinv(2,2)*dgphi_dtheta 
      mat4x4(3,4)=-0.5d0*gmetinv(2,2)*dgphit_dtheta 
      mat4x4(4,1)=0.0d0 
      mat4x4(4,2)=0.0d0 
      mat4x4(4,3)=mat4x4(3,4) 
      mat4x4(4,4)=-0.5d0*gmetinv(2,2)*dgt_dtheta 
      call matvec(vec4x1,mat4x4,dxvec,4) 
      useval=-dot4(vec4x1,dxvec) 
      d2theta_dtau2=useval 
 
 
 
 
c 3rd, compute d2phi_dtau2 
c do with Chris symb 
c these two do loops ensure mat4x4 set to zero matrix 
       do i = 1,4 
       do j = 1,4 
          mat4x4(i,j)=0.0d0 
       enddo 
       enddo 
      mat4x4(1,1)=0.0d0 
      mat4x4(1,2)=0.0d0 
      mat4x4(1,3)=0.5d0*gmetinv(3,3)*dgphi_dr 
      mat4x4(1,3)=mat4x4(1,3)+0.5d0*gmetinv(3,4)*dgphit_dr 
      mat4x4(1,4)=0.5d0*gmetinv(3,3)*dgphit_dr 
      mat4x4(1,4)=mat4x4(1,4)+0.5d0*gmetinv(3,4)*dgt_dr 
      mat4x4(2,1)=mat4x4(1,2) 
      mat4x4(2,2)=0.0d0 
      mat4x4(2,3)=0.5d0*gmetinv(3,3)*dgphi_dtheta 
      mat4x4(2,3)=mat4x4(2,3)+0.5d0*gmetinv(3,4)*dgphit_dtheta 
      mat4x4(2,4)=0.5d0*gmetinv(3,3)*dgphit_dtheta 
      mat4x4(2,4)=mat4x4(2,4)+0.5d0*gmetinv(3,4)*dgt_dtheta 
      mat4x4(3,1)=mat4x4(1,3) 
      mat4x4(3,2)=mat4x4(2,3) 
      mat4x4(3,3)=0.0d0 
      mat4x4(3,4)=0.0d0 
      mat4x4(4,1)=mat4x4(1,4) 
      mat4x4(4,2)=mat4x4(2,4) 
      mat4x4(4,3)=0.0d0 
      mat4x4(4,4)=0.0d0 
      call matvec(vec4x1,mat4x4,dxvec,4) 
      useval=-dot4(vec4x1,dxvec) 
      d2phi_dtau2=useval 
 
 
c 4th, compute d2t_dtau2 
c do with Chris symb 
c these two do loops ensure mat4x4 set to zero matrix 
      do i = 1,4 
      do j = 1,4 
         mat4x4(i,j)=0.0d0 
      enddo 
      enddo 
      mat4x4(1,1)=0.0d0 
      mat4x4(1,2)=0.0d0 
      mat4x4(1,3)=0.5d0*gmetinv(4,3)*dgphi_dr 
      mat4x4(1,3)=mat4x4(1,3)+0.5d0*gmetinv(4,4)*dgphit_dr 
      mat4x4(1,4)=0.5d0*gmetinv(4,3)*dgphit_dr 
      mat4x4(1,4)=mat4x4(1,4)+0.5d0*gmetinv(4,4)*dgt_dr 
      mat4x4(2,1)=0.0d0 
      mat4x4(2,2)=0.0d0 
      mat4x4(2,3)=0.5d0*gmetinv(4,3)*dgphi_dtheta 
      mat4x4(2,3)=mat4x4(2,3)+0.5d0*gmetinv(4,4)*dgphit_dtheta 
      mat4x4(2,4)=0.5d0*gmetinv(4,3)*dgphit_dtheta 
      mat4x4(2,4)=mat4x4(2,4)+0.5d0*gmetinv(4,4)*dgt_dtheta 
      mat4x4(3,1)=mat4x4(1,3) 
      mat4x4(3,2)=mat4x4(2,3) 
      mat4x4(3,3)=0.0d0 
      mat4x4(3,4)=0.0d0 
      mat4x4(4,1)=mat4x4(1,4) 
      mat4x4(4,2)=mat4x4(2,4) 
      mat4x4(4,3)=0.0d0 
      mat4x4(4,4)=0.0d0 
      call matvec(vec4x1,mat4x4,dxvec,4) 
      useval=-dot4(vec4x1,dxvec) 
      d2t_dtau2=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 
c 
c 
      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 
c  
      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 
c 
      real*8 function rnorm(x) 
      implicit double precision (a-h,o-z) 
      dimension x(3) 
      rnorm=sqrt(dot(x,x)) 
      return 
      end 
c 
c  
C 
C 
	SUBROUTINE MATINV(M1,A,DET) 
        implicit double precision (a-h,o-z) 
	DIMENSION A(4,4),IPVT(4),PVT(4),IND(4,2) 
C 
C 
	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 
C 
C 
200	RETURN 
	END 
c 
 
c 
        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 
c 
        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 
C 
C 
       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 
       
 
       subroutine crsprd(x,y,z) 
       implicit double precision (a-h,o-z) 
       dimension x(3),y(3),z(3) 
       z(1)=x(2)*y(3)-x(3)*y(2) 
       z(2)=x(3)*y(1)-x(1)*y(3) 
       z(3)=x(1)*y(2)-x(2)*y(1) 
       return 
       end 
 
 
