https://github.com/florentrenaud/nbody6tt
Raw File
Tip revision: 8a4716382ead3ece116c48a4ae5c65f8c9534437 authored by Florent on 29 January 2015, 12:19:28 UTC
Nbody6 - 29 January 2015 (added GPU2/Build/.gitkeep)
Tip revision: 8a47163
ksapo.f
      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 IF (IPAIR.LT.0) THEN
          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
      ELSE
*       Include small angle for moving away from pericentre.
          THETA = 1.0
          IKICK = -1
          IPAIR = KSPAIR
*       Increase angle near small pericentre.
          SEMI = -0.5D0*BODY(N+IPAIR)/H(IPAIR)
          ECC2 = (1.0 - R(IPAIR)/SEMI)**2 +
     &                          TDOT2(IPAIR)**2/(BODY(N+IPAIR)*SEMI)
          ECC = SQRT(ECC2)
          IF (ECC.LT.1.0) THEN
              EFAC = SQRT((1.0 + ECC)/(1.0 - ECC))
              THETA = 2.0*ATAN(EFAC)
          END IF
          WRITE (6,66)  ECC, THETA, SEMI, SEMI*(1.0 - ECC)
   66     FORMAT (' KSAPO    E THETA A PM ',F10.6,F7.3,1P,2E10.2)
      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
back to top