# DATA EXAMPLE:
  0  5   .1d0    20000.   0.e-3    0.0   1. 1.e-99 0.e0 1.e4  X  1 0  0.0 .70 0.3  1.e-13 0  10000 /IWR,N,DELTAT,TMAX ,step,soft,cmet,Clight, outfile,Ixc,Nbh ,spin, EPS, ivelocity, KSMX ! AT THE END OF THIS EXAMPLE FILE THERE ARE SOME EXPLANATIONS, please take a look!
  1500000  0.  0.  0.  0.  0.  0. 0    m x y z vx vy vz  soft (and the same for all bodies)
  50000.  39.399994  9.299988  0. -.5800018  .97399998  0.  
  10.  191.  32.5  0. -.55600004  .84599991  0. 
  10.  189. -84.700012  0.  .0033600006  .0017600002  0. 
  0  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0  /STOP ; this would set N=4 even if in first line N=5 !! REMOVE THIS IF U WANT MORE BODIES
  10.  94.099976 -122.  0.  .90300007  .07999992  0. 
  10.  41.799988 -13.  0.  .04699993  0.964799976  0. 
  0  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 /STOP | 

# PLEASE read these expalantions and ALSO look the comments in the code.

## FIRST DATA LINE: 

 UNITS: G=1,  I usually use Solar masses, astronomical units (AU), then the time is such that:  1 Year=2 Pi, i.e. years=time/twopi

 IWR =output index (the larger the more diagnostic output)

 N=number of bodies

 DELTAT = (approx) output interval

 TMAX =maximum time

 step = initial stepsize  (if step=0, the code attempts to estimate a reasonable one), NORMALLY =0

 soft = softening (code works best for soft=0, but softening is not prohibited), NORMALLY = 01

 cmet =three-dimesional vector that defines the used time-transformation  1 0 0 may usually be the best
         (but a small value for cmet(2) may be advisable in case of star-star near collision).

 Clight= speed of light. The code system of units is such that G=1. If AU and M_sun are =1, then CLIGHT \approx 10^4.

 outfile= name of outputfile for time and coords

 IXC= index that tels if exact output times are required, 
      0 => no exact, output when t>tnext, (fast, but sometimes output interval too long) 
      1 => exact time (not always accureate, but faster than the exact iteration)
      2 => exact iteration (time to tnext) (slower!)

 NBH= number of black holes, (which MUST be the bodies # 1,2,.., NBH)

 spin= three-dimensional vector= the (relative) spin of the BH = M1 (i.e. |psin|<1, and only the first body 
       is allowed to be a Kerr-hole, other possible BH's are non-rotating).
 
 EPS = error tolerance for the integration

 ivelocity = index telling that there are v-dependent perturbations (=1) or not (=0)

 KSMX = tells how many (maximum) integration steps the code is at most allowd to take between outputs 
        (this avoidsi almost infinity loops is step has got an unreasonably small valuei, from which it may not recover)

## ---------   END OF FIRST DATA LINE EXPLANATION ----------------

 AFTER the first line of data there is the list of masses and initial conditions
    m x y z  vx vy vz
 for all the bodies from 1 to N.


## ------------------------------------------------------------

 After finishin one simulation the code attempt to read data for the next one.
 It stops when the 1. line is zeros (actually if N<2).

# OUTPUT: (at the moment)

The coordinates goto the OUTPUT file (coords for merged bodies get huge values for moviei/picture makingi i.e. they go out of figure)  
        write(66,234)time,
     & (xwr(k),xwr(k+1),xwr(k+2),k=1,3*n_ini,3)

In addition orbital elements (Keplerian ones) are written to files 71, 72,.. 75 as explaned here:
           write(71,171)time,(ai(k),k=2,N_ini) ! a   write orbital elements (with respect to M1)
           write(72,171)time,(ei(k),k=2,N_ini) ! e
           write(73,171)time,(unci(k),k=2,N_ini) ! i
           write(74,171)time,(Omi(k),k=2,N_ini)  ! \Omega
           write(75,171)time,(ooi(k),k=2,N_ini)  ! \omega
           write(76,*)time,spin,dsp ! spin(k), k=1,3 of M1  (|spin|<1),
                                    ! dsp is error in the length of the spin vector
 

 One may change the main program easily to get other quantities computed from the coords/vels that 'ARC' returns.
 Note that the number of bodies may change if the BH swollows some of them.
  
#  LIST OF SUBROUTINES & FUNCTIONS

         subroutine ARCparams_Dot_CH_is_here ! This subroutine is not used, but the content is simply what is supposed to be in the file  ARCparams.CH.
         Program ARCcode !  This is the main program, not subroutine.
         subroutine Diagnostic output(time,xwr,cmxx)
         subroutine Write Elements(time,iopen) ! with respect to M(1)
         Subroutine MERGE_I1_I2(time)!Merge the colliding body with the BH, ('time' is here only for write(66....)
         SUBROUTINE ARC ! This is the original Alforithmic Regularization Chain code.
         subroutine Iterate2ExactTime(Y0,Nvar,deltaT,f1,d1,f2,d2,x1,x2)
         subroutine LEAPFROG(STEP,Leaps,stime) !Leapfroging steps for the Bulirsch-Stoer extrapolation
         function Wfunction() ! avluate the TTL-function for time tranformation !  not  subroutine
         subroutine omegacoef ! coefficients in \Omega=sum omec(i,j)/r_{ij}
         SUBROUTINE XCMOTION(hs) ! movement of coordinates
         subroutine PUT V 2 W ! these parameters are needed if there are v-ependent forces
         subroutine Velocity Dependent Perturbations ! This name tells what this is!
         SUBROUTINE CHECK SWITCHING CONDITIONS(MUSTSWITCH)! Is it necessary to construct a new CHAIN?
         SUBROUTINE FIND CHAIN INDICES ! Here the new indecies along the new CHAIN are obtained
         SUBROUTINE CHECK CONNECTION(IC,LMIN,LMAX,IJ,LI,SUC,LOOP) ! Loops in CHAIN are not allowed
         SUBROUTINE ARRANGE(N,Array,Indx) ! this is a sorting routine (needed in CHAIN construction)
         SUBROUTINE INITIALIZE XC AND WC  ! find CHAIN coordinates and velocities
         SUBROUTINE UPDATE X AND V ! Normal physical x & v from CHAIN XC and WC
         SUBROUTINE CHAIN TRANSFORMATION ! New CHAIN from the old one
         FUNCTION SQUARE(X,Y) ! |X-Y|^2 !  not  subroutine
         SUBROUTINE DIFSYAB(N,EPS,S,h,t,Y) ! BULIRSCH-STOER EXTRAPOLATOR (i.e. with rationals)
         subroutine SubSteps(Y0,Y,H,Leaps)! substeps for DIFSYAB
         subroutine Initialize increments 2 zero
         subroutine Take Increments 2 Y(Y)
         subroutine Put Y to XC WC (Y,Lmx)
         subroutine Take Y from XC WC (Y,Nvar)
         subroutine Obtain Order of Y(SY) ! an attempt to obtain reasonable comparison values for
         SUBROUTINE EVALUATE X
         SUBROUTINE EVALUATE V(VN,WI)
         SUBROUTINE Relativistic ACCELERATIONS(ACC,ACCGR,Va,spina,dspin)
         subroutine Relativistic terms_not in use ! at the moment this is not used
         subroutine Relativistic terms!_ in use ! this is used for PN-term accelarations
         subroutine Reduce2cm(x,m,nb,cm)
         subroutine cross(a,b,c)
         subroutine gopu_SpinTerms(X,V,r,M1,m2,c,alpha,dv3,dalpha)
         subroutine Q2term(m,r,x,v,c,Q2,e,dvq) ! Quadrupole as C. Will advised
         SUBROUTINE Initial Stepsize(X,V,M,NB,ee,step) ! guesswork for initial stepsize
         subroutine elmnts ! two-body orbital elements (WITHOUT relativistic corrections)
         function atn2(s,c) ! atan 4 interval (0,2Pi) !  not  subroutine
         function oot(alfa,eta,zeta,q,e,sqaf) ! oot=pericentre time !  not  subroutine
         function g3(z) !  not  subroutine 
         SUBROUTINE CONSTANTS OF MOTION(ENE_NB,G,Alag) ! Newtonian constants of motion (i.e. no PN)
         SUBROUTINE FIND BINARIES(time)  ! this is a toy routine for finding binaries
         SUBROUTINE  WCMOTION(hs)  ! This routine evaluates the velocity jumps in leapfrog
         subroutine V_jump(WCj,spinj,cmvj,wttlj,WCi,spini,FCj,acc,dt,
         subroutine V_jACConly(WCj,CMVj,WTTLj,FC,acc,dt,
         subroutine Estimate Stepsize(dtime,step)!Stumpff-Weiss method to estimate the s-step
        subroutine gfunc(xb,al,g) ! G_n=xb^n c_n(al*xb^2)
        SUBROUTINE cfun(z,c)!Stumpff c-functions
        subroutine Xanom(m,r,eta,zet,beta,t,x,rx,g) ! Solving Kepler's equation
        function cdot(a,b) !  not  subroutine 
         subroutine  Non relativistic v_dependent perturbations(df) ! USER DEFINED ( to be programmed by the user)
        subroutine COORDINATE DEPENDENT PERTURBATIONS(ACC) ! USER DEFINED