freeze.f

Go to the documentation of this file.

      SUBROUTINE freeze(IPAIR)
*
*
*       Partial reflection of KS binary.
*       --------------------------------
*
      include 'common6.h'
      SAVE nksref
      DATA nksref /0/
*
*
*       Set c.m. & component indices and form semi-major axis.
      icm = n + ipair
      i2 = 2*ipair
      i1 = i2 - 1
      semi = -0.5d0*body(icm)/h(ipair)
*
*       Define auxiliary quantities and obtain the eccentricity.
      zeta = 1.0 - r(ipair)/semi
      psi = tdot2(ipair)/sqrt(body(icm))
      ecc = sqrt(zeta**2 + psi**2/semi)
*
*       Skip reflection of hyperbolic orbit (just in case).
      IF (ecc.GE.1.0) go to 100
*
*       Determine the eccentric anomaly with respect to pericentre (0,PI).
      theta = atan2(abs(psi)/sqrt(semi),zeta)
*
*       Obtain total reflection time from Kepler's equation.
      dt = 2.0d0*semi*sqrt(semi/body(icm))*(theta - abs(psi)/sqrt(semi))
*
*       Specify regularized time (based on Baumgarte & Stielel, 1974).
      dtu = -2.0d0*(h(ipair)*dt + tdot2(ipair))/body(icm)
*
*       Skip reflection near pericentre.
      IF (dtu.LT.4.0*dtau(ipair)) go to 100
*
*       Predict c.m. coordinates and resolve unreflected KS components.
      CALL xvpred(icm,0)
*
*       Specify transformation coefficients (Mikkola's proceure).
      xc = zeta/ecc
      ys = 2.0d0*psi/(ecc*sqrt(body(icm)))
      zz = body(icm)/(4.0*semi)
      r(ipair) = 0.0d0
*
*       Generate analytical solutions for U & UDOT using old U & UDOT.
      DO 10 k = 1,4
          u(k,ipair) = u0(k,ipair)*xc - udot(k,ipair)*ys
          udot(k,ipair) = u0(k,ipair)*ys*zz + udot(k,ipair)*xc
          u0(k,ipair) = u(k,ipair)
          r(ipair) = r(ipair) + u(k,ipair)**2
*       Consistent R = U*U for stabilization procedure.
   10 CONTINUE
*
*       Copy significant perturbers for tidal correction (R-dependent).
      corr = 2.0 + 2.0*(semi - r(ipair))/(semi*(1.0 + ecc))
      rcrit2 = cmsep2*(corr*r(ipair))**2
      nnb1 = list(1,i1) + 1
      np = 0
      DO 20 l = 2,nnb1
          j = list(l,i1)
          rij2 = (x(1,icm) - x(1,j))**2 + (x(2,icm) - x(2,j))**2 +
     &                                    (x(3,icm) - x(3,j))**2
          IF (rij2.LT.rcrit2) THEN
              np = np + 1
              jpert(np) = j
          END IF
   20 CONTINUE
*
*       Ensure that at least one perturber is retained.
      IF (np.EQ.0) THEN
          np = 1
          jpert(np) = list(2,i1)
      END IF
*
*       Reverse second time derivative and specify unperturbed motion.
      tdot2(ipair) = -tdot2(ipair)
      list(1,i1) = 0
*
*       Set new step in physical units and increase counter.
      step(i1) = dt
      nksref = nksref + 1
*
*       Check minimum two-body distance.
      dmin2 = min(dmin2,semi*(1.0d0 - ecc))
*
*       Obtain potential energy with respect to the components.
      jlist(1) = i1
      jlist(2) = i2
      CALL nbpot(2,np,pot1)
*
*       Transform to reflected coordinates and repeat the summation.
      CALL ksres(ipair,j1,j2,rij2)
      CALL nbpot(2,np,pot2)
*
*       Form correction factor due to the tidal potential.
      vi2 = x0dot(1,icm)**2 + x0dot(2,icm)**2 + x0dot(3,icm)**2
      corr = 1.0 + 2.0*(pot2 - pot1)/(body(icm)*vi2)
      etcorr = etcorr + (pot2 - pot1)
      IF (corr.GT.0.0d0) THEN
          corr = sqrt(corr)
      ELSE
          corr = 0.0d0
      END IF
*
*       Modify the c.m. velocity.
      DO 30 k = 1,3
          x0dot(k,icm) = corr*x0dot(k,icm)
   30 CONTINUE
*
  100 RETURN
*
      END