ksapo.f

Go to the documentation of this file.

      SUBROUTINE ksapo(IPAIR)
*
*
*       Apocentre/pericentre/random KS variables.
*       -----------------------------------------
*
      include 'common6.h'
      REAL*8  ran2
*
*
*       Specify half regularized period (PI/2) or random phase (unperturbed).
      IF (ipair.GT.0) THEN
          theta = 0.25d0*twopi
          ikick = 0
      ELSE
          theta = 0.5*twopi*ran2(idum1)
          ipair = -ipair
*       Note new type not known here but WD case kick decided by option #25.
          ikick = 1
          IF (list(1,2*ipair-1).GT.0) theta = 0.0d0
*       Initialize time because HDOT & TDOT2 not updated for RESOLV.
          t0(2*ipair-1) = time
*       Skip hyperbolic orbit (i.e. kick for second binary component).
          IF (h(ipair).GT.0.0) go to 30
      END IF
*
*       Form transformation coefficients (Stiefel & Scheifele p. 85).
      xc = cos(theta)
      ys = sin(theta)
      ff = sqrt(0.5d0*abs(h(ipair)))
      r(ipair) = 0.0d0
      tdot2(ipair) = 0.0d0
*
*       Generate analytical solutions for U & UDOT using old U0 & UDOT.
      DO 10 k = 1,4
          u(k,ipair) = u0(k,ipair)*xc + udot(k,ipair)*ys/ff
          udot(k,ipair) = udot(k,ipair)*xc - u0(k,ipair)*ys*ff
          u0(k,ipair) = u(k,ipair)
          r(ipair) = r(ipair) + u(k,ipair)**2
          tdot2(ipair) = tdot2(ipair) + 2.0d0*u(k,ipair)*udot(k,ipair)
   10 CONTINUE
*
*       Impose R' < 0 for apocentre procedures (IKICK = 0).
      semi = -0.5d0*body(n+ipair)/h(ipair)
      IF (tdot2(ipair).GT.0.0d0.AND.r(ipair).GT.semi) THEN
          IF (ikick.EQ.0) THEN
              tdot2(ipair) = -1.0e-20
          END IF
      END IF
*
*       Include diagnostic check that correct apocentre has been set.
*     SEMI = -0.5D0*BODY(N+IPAIR)/H(IPAIR)
*     WRITE (6,20)  SEMI, R(IPAIR), H(IPAIR), GAMMA(IPAIR)
*  20 FORMAT (' APOCENTRE:    A RA H G ',1P,2E10.2,2E10.1)
*
*       Save KS parameters for WD or neutron star kick (routine FCORR).
   30 IF (ikick.GT.0) THEN
          CALL kick(ipair,0,0,0.0d0)
      END IF
*
      RETURN
*
      END