SUBROUTINE INDUCE(IPAIR,JCL,EMAX,EMIN,ICIRC,TC2,ANGLE,TG,EDAV) * * * Induced eccentricity of hierarchical binary. * -------------------------------------------- * INCLUDE 'common6.h' REAL*8 A1(3),A2(3),XREL(3),VREL(3),EI(3),HI(3),HO(3),BHAT(3) * * * Define c.m. & KS indices and evaluate semi-major axis & eccentricity. I = N + IPAIR I1 = 2*IPAIR - 1 I2 = I1 + 1 SEMI = -0.5d0*BODY(I)/H(IPAIR) ECC2 = (1.d0 - R(IPAIR)/SEMI)**2 + TDOT2(IPAIR)**2 & /(BODY(I)*SEMI) ECC = SQRT(ECC2) * * Ensure that KS components are resolved (perturbed not always done). CALL RESOLV(IPAIR,1) * * Obtain inner & outer angular momentum by cyclic notation. RIJ2 = 0.d0 VIJ2 = 0.d0 RDOT = 0.d0 A12 = 0.d0 A22 = 0.d0 A1A2 = 0.d0 RI2 = 0.d0 VI2 = 0.d0 RVI = 0.d0 DO 5 K = 1,3 RIJ2 = RIJ2 + (X(K,I) - X(K,JCL))**2 VIJ2 = VIJ2 + (XDOT(K,I) - XDOT(K,JCL))**2 RDOT = RDOT + (X(K,I) - X(K,JCL))*(XDOT(K,I) -XDOT(K,JCL)) K1 = K + 1 IF (K1.GT.3) K1 = 1 K2 = K1 + 1 IF (K2.GT.3) K2 = 1 A1(K) = (X(K1,I1) - X(K1,I2))*(XDOT(K2,I1) - XDOT(K2,I2)) & - (X(K2,I1) - X(K2,I2))*(XDOT(K1,I1) - XDOT(K1,I2)) A2(K) = (X(K1,JCL) - X(K1,I))*(XDOT(K2,JCL) - XDOT(K2,I)) & - (X(K2,JCL) - X(K2,I))*(XDOT(K1,JCL) - XDOT(K1,I)) A12 = A12 + A1(K)**2 A22 = A22 + A2(K)**2 A1A2 = A1A2 + A1(K)*A2(K) * Form relative vectors and scalars for inner binary. XREL(K) = X(K,I1) - X(K,I2) VREL(K) = XDOT(K,I1) - XDOT(K,I2) RI2 = RI2 + XREL(K)**2 VI2 = VI2 + VREL(K)**2 RVI = RVI + XREL(K)*VREL(K) 5 CONTINUE * * Evaluate orbital parameters for outer orbit. RIJ = SQRT(RIJ2) ZMB = BODY(I) + BODY(JCL) SEMI1 = 2.d0/RIJ - VIJ2/ZMB SEMI1 = 1.d0/SEMI1 ECC1 = SQRT((1.d0 - RIJ/SEMI1)**2 + RDOT**2/(SEMI1*ZMB)) * PMIN1 = SEMI1*(1.d0 - ECC1) * Determine inclination (8 bins of 22.5 degrees). FAC = A1A2/SQRT(A12*A22) FAC = ACOS(FAC) ANGLE = FAC *** ANGLE = FAC*360.d0/TWOPI IN = 1 + FAC*360.0/(TWOPI*22.5) * * Exit on hyperbolic orbit or outer eccentricity near 1. IF (SEMI1.LT.0.0.OR.ECC1.GT.0.9999) THEN EMAX = 0.d0 TC2 = 999.d0 TG = 1.0d+04 EDAV = 1.0d+04 GO TO 30 END IF * * Construct the Runge-Lenz vector (Heggie & Rasio 1995, Eq.(5)). EI2 = 0.d0 DO 10 K = 1,3 EI(K) = (VI2*XREL(K) - RVI*VREL(K))/BODY(I) - & XREL(K)/SQRT(RI2) EI2 = EI2 + EI(K)**2 10 CONTINUE * * Define unit vectors for inner eccentricity and angular momenta. COSJ = 0.d0 SJSG = 0.d0 DO 15 K = 1,3 EI(K) = EI(K)/SQRT(EI2) HI(K) = A1(K)/SQRT(A12) HO(K) = A2(K)/SQRT(A22) COSJ = COSJ + HI(K)*HO(K) SJSG = SJSG + EI(K)*HO(K) 15 CONTINUE * * Form unit vector BHAT and scalars AH & BH (Douglas Heggie, 10/9/96). AH = 0.d0 BH = 0.d0 DO 20 K = 1,3 K1 = K + 1 IF (K1.GT.3) K1 = 1 K2 = K1 + 1 IF (K2.GT.3) K2 = 1 BHAT(K) = HI(K1)*EI(K2) - HI(K2)*EI(K1) AH = AH + EI(K)*HO(K) BH = BH + BHAT(K)*HO(K) 20 CONTINUE * * Evaluate the expressions A & Z. A = COSJ*SQRT(1.d0 - EI2) Z = (1.d0 - EI2)*(2.d0 - COSJ**2) + 5.d0*EI2*SJSG**2 * * Obtain maximum inner eccentricity (Douglas Heggie, Sept. 1995). Z2 = Z**2 + 25.d0 + 16.d0*A**4 - 10.d0*Z - & 20.d0*A**2 - 8.d0*A**2*Z EMAX = ONE6*(Z + 1.d0 - 4.d0*A**2 + SQRT(Z2)) EMAX = MAX(EMAX,0.0001d0) EMAX = SQRT(EMAX) ZI = FAC EM = SQRT(SIN(ZI)**2 + EI2*COS(ZI)**2) * * Form minimum eccentricity (Douglas Heggie, Sept. 1996). AZ2 = A**2 + Z - 2.d0 IF (AZ2.GE.0.0) THEN AZ1 = 1.d0 + Z - 4.d0*A**2 EMIN2 = ONE6*(AZ1 - SQRT(AZ1**2 - 12.d0*AZ2)) ELSE EMIN2 = 1.d0 - 0.5d0*(A**2 + Z) END IF EMIN2 = MAX(EMIN2,0.0001d0) EMIN = SQRT(EMIN2) * * Set impact parameters and estimate eccentricity growth time-scale. PMIN = SEMI*(1.d0 - ECC) PMIN2 = SEMI*(1.d0 - EMAX) TK = TWOPI*SEMI*SQRT(SEMI/BODY(I)) TK1 = TWOPI*SEMI1*SQRT(SEMI1/ZMB) TG = TK1**2*ZMB*(1.d0 - ECC1**2)**1.5/(BODY(JCL)*TK) * CONST = PFAC(A,Z) CONST = CONST*4.d0/(1.5d0*TWOPI*SQRT(6.d0)) * Convert growth time to units of 10**6 yrs. TG = CONST*TG*TSTAR * * Form doubly averaged eccentricity derivative (Douglas Heggie 9/96). YFAC = 15.d0*BODY(JCL)/(4.d0*ZMB)*TWOPI*TK/TK1**2 YFAC = YFAC*ECC*SQRT(1.d0 - ECC**2)/ & (1.d0 - ECC1**2)**(1.5) EDAV = YFAC*AH*BH * * Determine circularization time for current & smallest pericentre. TC2 = 2000.d0 TC = 2000.d0 IF (PMIN2.LT.10.0*MAX(RADIUS(I1),RADIUS(I2)).AND. & KZ(27).EQ.2.AND.ECC.GT.0.002) THEN CALL TCIRC(PMIN,ECC,I1,I2,ICIRC,TC) IF (ICIRC.GE.0) ICIRC = -1 * Note: only use PMIN call because ICIRC may become >0 after first! CALL TCIRC(PMIN2,EMAX,I1,I2,ICIRC,TC2) ELSE IF (ECC.LT.0.002) THEN TC = 0.0 TC2 = 0.0 END IF * * Print diagnostics for high inclinations and TC2 < 10**7 yrs. IF (IN.GE.4.AND.IN.LE.5.AND.TC2.LT.10.0.AND.KSTAR(I).NE.-2.AND. & EI2.GT.4.0D-06) THEN WRITE (44,25) SQRT(EI2), EM, EMAX, KSTAR(I1), KSTAR(I2), & KSTAR(I), SEMI, PMIN2, IN, TG, TC, TC2, TPHYS 25 FORMAT (' INDUCE: E EMX EMAX K* SEMI PM2 IN TG TC TC2 TP ', & 3F8.4,3I3,1P,2E9.1,0P,I3,F7.3,3F7.1) CALL FLUSH(44) END IF * * Use first value for now. TC2 = TC * 30 RETURN * END