start.f

Go to the documentation of this file.

      SUBROUTINE start
*
*
*       Initialization of data & polynomials.
*       ------------------------------------
*
      include 'common6.h'
      EXTERNAL scale,merge
      parameter(ns=12)
*
*
*       Initialize global scalars, counters & useful constants.
      CALL zero
*
*       Read input parameters.
      CALL input
*
*       Set initial conditions: BODY(I), X(K,I), XDOT(K,I); I=1,N & K=1,3.
      CALL DATA
*
*       Scale initial conditions to new units.
      CALL scale
*
*       Set total mass in case routines DATA & SCALE are not used.
      zmass = 0.0d0
      DO 10 i = 1,n
          zmass = zmass + body(i)
   10 CONTINUE
*
*       Define mean mass in scaled units and solar mass conversion factor.
      bodym = zmass/float(n)
      IF (kz(5).NE.3) THEN
          zmbar = zmbar/bodym
      END IF
*
*       Introduce scaling factors DAYS, YRS, SU, RAU, SMU, TSTAR & VSTAR.
      CALL units
*
*       Check option for external force.
      IF (kz(14).GT.0) THEN 
          CALL xtrnl0
      END IF 
*
*       Check optional scaling to hot system.
      IF (kz(29).GT.0) THEN
          CALL hotsys
      END IF
*
*       Check option for initial binaries.
      IF (kz(8).EQ.1.OR.kz(8).GE.3) THEN
          CALL binpop
      END IF
*
*       Include stable primordial triples.
      IF (kz(18).GT.1.AND.kz(8).GT.0) THEN
          CALL hipop
      END IF
*
*       Check optional initialization for tidal two-body capture (suppress).
*     IF (KZ(27).GT.0) THEN
*         CALL INTIDE
*     END IF
*
*       Set sequential name, maximum mass & primary velocity.
      body1 = 0.0
      DO 20 i = 1,n
          name(i) = i
          body1 = max(body1,body(i))
          DO 15 k = 1,3
              x0dot(k,i) = xdot(k,i)
   15     CONTINUE
   20 CONTINUE
*
*       Randomize particle indices (dummy routine for standard code).
      CALL swap
*
*       Initialize fixed block steps (40 levels).
      CALL iblock
*
*       Create table of inverse Stumpff coefficients.
      DO 30 i = 1,ns
          scoeff(i) = 1.0d0/((i + 1)*(i + 2))
   30 CONTINUE
*
*       Set optional stellar evolution parameters or define STEPX.
      IF (kz(19).GT.2) THEN
          CALL instar
      ELSE IF (kz(14).GT.1) THEN
          dt = 1.0e-03/tscale
          CALL stepk(dt,dtn)
          stepx = dtn
      END IF
*
*       Initialize optional cloud parameters.
      IF (kz(13).GT.0) THEN
          CALL cloud0
      END IF
*
*       Set initial neighbour list & corresponding radius.
      rs0 = rc
      nnb0 = 0
      DO 40 i = 1,n
          CALL nblist(i,rs0)
          nnb0 = nnb0 + list(1,i)
   40 CONTINUE
*
*       Obtain force & first derivative.
      CALL fpoly1(1,n,0)
*
*       Obtain second & third force derivatives and set time-steps.
      CALL fpoly2(1,n,0)
*
*       Regularize any hard primordial binaries (assume sequential ordering).
      IF (nbin0.GT.0) THEN
          DO 50 ipair = 1,nbin0
              icomp = 2*ipair - 1
              jcomp = 2*ipair
              rij2 = 0.0
*       Include standard distance criterion.
              DO 45 k = 1,3
                  rij2 = rij2 + (x(k,icomp) - x(k,jcomp))**2
   45         CONTINUE
              IF (rij2.LT.rmin**2) THEN
                  CALL ksreg
              END IF
   50     CONTINUE
      END IF
*
*       Include optional regularization of primordial triples.
      IF (kz(18).GT.1.AND.nhi0.GT.0) THEN
          kspair = 1
*       Note that each KS pair will move to the top of the queue.
   60     icomp = 2*kspair - 1
          icm = kspair + n
          rx2 = 1.0
*       Find index of closest outer component without any assumption.
          DO 70 j = ifirst,n
              rij2 = 0.0
              DO 65 k = 1,3
                  rij2 = rij2 + (x(k,icm) - x(k,j))**2
   65         CONTINUE
              IF (rij2.LT.rx2) THEN
                  rx2 = rij2
                  jcomp = j
              END IF
   70     CONTINUE
          IF (sqrt(rx2).GT.rmin) THEN
              kspair = kspair + 1
              IF (kspair.GT.npairs) go to 80
              go to 60
          END IF
*       Evaluate PCRIT for R0(NPAIRS) in MERGE since IMPACT is bypassed.
          CALL histab(kspair,jcomp,pmin,rstab)
*       Initialize the triple (constructed to be stable in HIPOP).
          iphase = 6
          CALL merge
*       Examine the same ICM (moved up after successful MERGE).
          IF (nmerge.LT.nhi0) THEN
              go to 60
          END IF
      END IF
*
*       Check the average neighbour number.
   80 znb = float(nnb0)/float(n)
      IF (znb.LT.0.25*znbmax.OR.znb.LT.0.25*sqrt(float(n))) THEN
          WRITE (6,90)  znb
   90     FORMAT (/,12x,'WARNING!   SMALL NEIGHBOUR NUMBERS   <NNB> =',
     &                                                             f6.1)
      END IF
*
      RETURN
*
      END