kspoly.f

Go to the documentation of this file.

      SUBROUTINE kspoly(IPAIR,IMOD)
*
*
*       Initialization of KS polynomials.
*       ---------------------------------
*
      include 'common6.h'
      common/slow0/  range,islow(10)
      REAL*8  a1(3,4),a(8),fp(4),fpdot(3),ui(4),uidot(4)
      REAL*8  b2(3,4),u2(4),q2(4)
*
*
*       Specify indices of c.m. & components and perturber membership (+ 1).
      i = n + ipair
      i2 = 2*ipair
      i1 = i2 - 1
      nnb1 = list(1,i1) + 1
      semi = -0.5*body(i)/h(ipair)
*       Save new slow-down index (note direct CALL KSPOLY).
      kslow(ipair) = imod
*
*       Predict current coordinates & velocities of the perturbers.
      DO 5 l = 2,nnb1
          j = list(l,i1)
          CALL xvpred(j,0)
    5 CONTINUE
*
*       Ensure current values of existing KS solution.
      IF (iphase.EQ.0) THEN
          CALL resolv(ipair,1)
      END IF
*
*       Initialize variables for accumulating contributions.
      DO 10 k = 1,3
          fp(k) = 0.0d0
          fpdot(k) = 0.0d0
   10 CONTINUE
*       Set dummy index for summation of c.m. or resolved components.
      jdum = 0
*
*       Obtain the perturbation & first derivative.
      DO 40 l = 2,nnb1
          j = list(l,i1)
          IF (j.GT.n) THEN
*       See whether c.m. approximation applies.
              rij2 = (x(1,j) - x(1,i))**2 + (x(2,j) - x(2,i))**2 +
     &                                      (x(3,j) - x(3,i))**2
              k = j - n
              IF (rij2.LT.cmsep2*r(k)**2.AND.list(1,2*k-1).GT.0) THEN
                  jdum = 2*k - 1
                  j = jdum
              END IF
          END IF
*
*       Sum over both components (reverse sign for second component).
   20     ii = i1
          DO 30 kcomp = 1,2
              DO 25 k = 1,3
                  a(k) = x(k,j) - x(k,ii)
                  a(k+3) = xdot(k,j) - xdot(k,ii)
   25         CONTINUE
*       Current velocities are predicted in routines INTGRT, KSMOD, etc.
              rij2 = a(1)*a(1) + a(2)*a(2) + a(3)*a(3)
              a8 = body(j)/(rij2*sqrt(rij2))
              a9 = 3.0d0*(a(1)*a(4) + a(2)*a(5) + a(3)*a(6))/rij2
              IF (kcomp.EQ.2) a8 = -a8
*
              DO 28 k = 1,3
                  fp(k) = fp(k) + a(k)*a8
                  fpdot(k) = fpdot(k) + (a(k+3) - a(k)*a9)*a8
   28         CONTINUE
              ii = i2
   30     CONTINUE
*
          IF (j.EQ.jdum) THEN
              j = j + 1
              go to 20
          END IF
   40 CONTINUE
*
*       Check for external perturbations.
      IF (kz(14).GT.0.AND.kz(14).LT.3) THEN
          q1 = x(1,i1) - x(1,i2)
          q3 = x(3,i1) - x(3,i2)
          CALL xtrnlp(q1,q3,fp)
*
*       Use same formalism for the first derivative (omit Coriolis force).
          vx = xdot(1,i1) - xdot(1,i2)
          vz = xdot(3,i1) - xdot(3,i2)
          CALL xtrnlp(vx,vz,fpdot)
      END IF
*
*       Transform to regularized force derivative using T' = R.
      DO 45 k = 1,3
          fpdot(k) = r(ipair)*fpdot(k)
   45 CONTINUE
*
*       Save the relative perturbation.
      fp(4) = sqrt(fp(1)**2 + fp(2)**2 + fp(3)**2)
      gamma(ipair) = fp(4)*r(ipair)**2/body(i)
*
*       Copy U & UDOT to scalars and set current transformation matrix.
      DO 48 k = 1,4
          ui(k) = u(k,ipair)
          uidot(k) = udot(k,ipair)
   48 CONTINUE
      CALL matrix(ui,a1)
*
*       Construct regularized polynomials from explicit derivatives.
      tdot2(ipair) = 0.0d0
      tdot3(ipair) = 0.0d0
      tdot4 = 0.0d0
      tdot5 = 0.0d0
      tdot6 = 0.0d0
      hdot(ipair) = 0.0d0
      hdot2(ipair) = 0.0d0
      hdot3(ipair) = 0.0d0
      hdot4(ipair) = 0.0d0
*
*       Scale perturbing force & first derivative by modification factor.
      IF (imod.GT.1) THEN
          zmod = float(islow(imod))
          DO 50 k = 1,3
              fp(k) = zmod*fp(k)
              fpdot(k) = zmod*fpdot(k)
   50     CONTINUE
      END IF
*
*       Form regularized force & two derivatives of time & binding energy.
      DO 60 k = 1,4
          a(k) = a1(1,k)*fp(1) + a1(2,k)*fp(2) + a1(3,k)*fp(3)
          a(k+4) = a1(1,k)*fpdot(1) + a1(2,k)*fpdot(2) +
     &                                a1(3,k)*fpdot(3)
          fu(k,ipair) = 0.5d0*h(ipair)*u(k,ipair) + 0.5d0*r(ipair)*a(k)
          tdot2(ipair) = tdot2(ipair) + 2.0d0*u(k,ipair)*udot(k,ipair)
          tdot3(ipair) = tdot3(ipair) + 2.0d0*udot(k,ipair)**2 +
     &                                      2.0d0*u(k,ipair)*fu(k,ipair)
          hdot(ipair) = hdot(ipair) + 2.0d0*udot(k,ipair)*a(k)
          hdot2(ipair) = hdot2(ipair) + 2.0d0*fu(k,ipair)*a(k)
          fp0(k,ipair) = 0.5d0*r(ipair)*a(k)
          u2(k) = fu(k,ipair)
   60 CONTINUE
*
*       Set regularized velocity matrix (Levi-Civita matrix not required).
      CALL matrix(uidot,a1)
*
*       Include the whole (L*F)' term in explicit derivatives of FU & H.
      DO 65 k = 1,4
          a(k+4) = a(k+4) + a1(1,k)*fp(1) + a1(2,k)*fp(2) +
     &                                      a1(3,k)*fp(3)
          hdot2(ipair) = hdot2(ipair) + 2.0d0*udot(k,ipair)*a(k+4)
          fudot(k,ipair) = 0.5d0*(h(ipair)*udot(k,ipair) + 
     &                            hdot(ipair)*u(k,ipair) +
     &                            tdot2(ipair)*a(k) + r(ipair)*a(k+4))
          hdot3(ipair) = hdot3(ipair) + 2.0*fudot(k,ipair)*a(k) +
     &                                            4.0*fu(k,ipair)*a(k+4)
          tdot4 = tdot4 + 2.0*fudot(k,ipair)*u(k,ipair) +
     &                                     6.0*fu(k,ipair)*udot(k,ipair)
          fd0(k,ipair) = 0.5d0*(tdot2(ipair)*a(k) + r(ipair)*a(k+4) +
     &                                           hdot(ipair)*u(k,ipair))
   65 CONTINUE
*
*       Form L(U'') and add two of the three terms in (L*F)'' (11/1/98).
      CALL matrix(u2,b2)
      DO 66 k = 1,4
          q2(k) = b2(1,k)*fp(1) + b2(2,k)*fp(2) + b2(3,k)*fp(3) +
     &            2.0*(a1(1,k)*fpdot(1) + a1(2,k)*fpdot(2) +
     &                                    a1(3,k)*fpdot(3))
          hdot3(ipair) = hdot3(ipair) + 2.0*udot(k,ipair)*q2(k)
   66 CONTINUE
*
*       Form higher derivatives by boot-strapping (Q2 added 11/1/98).
      DO 70 k = 1,4
          fudot2(k,ipair) = 0.5*(h(ipair)*fu(k,ipair) +
     &                   hdot2(ipair)*u(k,ipair) + tdot3(ipair)*a(k)) +
     &                   tdot2(ipair)*a(k+4) + hdot(ipair)*udot(k,ipair)
          fudot3(k,ipair) = 0.5*(h(ipair)*fudot(k,ipair) + tdot4*a(k) + 
     &                        hdot3(ipair)*u(k,ipair)) +
     &                   1.5*(hdot(ipair)*fu(k,ipair) +
     &                        hdot2(ipair)*udot(k,ipair) +
     &                        tdot3(ipair)*a(k+4))
          fudot2(k,ipair) = fudot2(k,ipair) + 0.5*r(ipair)*q2(k)
          fudot3(k,ipair) = fudot3(k,ipair) + 1.5*tdot2(ipair)*q2(k)
          hdot3(ipair) = hdot3(ipair) + 2.0*udot(k,ipair)*q2(k)
          hdot4(ipair) = hdot4(ipair) + 2.0*fudot2(k,ipair)*a(k) +
     &                   6.0*(fudot(k,ipair)*a(k+4) + fu(k,ipair)*q2(k))
          tdot5 = tdot5 + 2.0*fudot2(k,ipair)*u(k,ipair) +
     &             8.0*fudot(k,ipair)*udot(k,ipair) + 6.0*fu(k,ipair)**2
          tdot6 = tdot6 + 2.0*fudot3(k,ipair)*u(k,ipair) +
     &                    10.0*fudot2(k,ipair)*udot(k,ipair) +
     &                    20.0*fudot(k,ipair)*fu(k,ipair)
   70 CONTINUE
*
*       Check maximum square step (soft binaries & weak hyperbolic pairs).
      IF (abs(h(ipair)).GT.eclose) THEN
          a2 = 0.5/abs(h(ipair))
      ELSE
          a2 = min(r(ipair)/body(i),0.5/abs(h(ipair)))
      END IF
*
*       Assign a conservative value of the initial step (except IMOD > 1).
      fac = 0.8
      IF (imod.GT.1.OR.h(ipair).GT.0.0) fac = 1.0
      dtu = fac*etau*sqrt(a2)/(1.0 + 1000.0*gamma(ipair))**0.333
*
*       Include convergence criterion DH = H'*DTU + H''*DTU**2/2 = 0.001*|H|.
      IF (imod.EQ.1.AND.hdot2(ipair).NE.0.0d0) THEN
          dh = 1.0e-03*max(abs(h(ipair)),eclose)
          xf = 2.0*dh/abs(hdot2(ipair))
          yf = hdot(ipair)/hdot2(ipair)
          dtu1 = sqrt(xf + yf**2) - abs(yf)
          dtu = min(dtu1,dtu)
      END IF
*
*       Initialize reference energy.
      h0(ipair) = h(ipair)
*
*       Skip Stumpff part and use ZMOD*TK for unperturbed binary.
      IF (list(1,i1).EQ.0.AND.semi.GT.0.0) THEN
          tk = twopi*semi*sqrt(semi/body(i))
          step(i1) = min(tk,step(i))
          go to 80
      END IF
*
*       Obtain Stumpff coefficients.
   75 z = -0.5d0*h(ipair)*dtu**2
      CALL stumpf(ipair,z)
      z5 = sf(6,ipair)
      z6 = sf(7,ipair)
      s6 = z6*tdot6/720.0
*
*       Convert regularized step to physical units and initialize time T0.
      step(i1) = ((((s6 + tdot5*dtu*z5/120.0d0 + tdot4/24.0d0)*dtu +
     &               tdot3(ipair)*one6)*dtu + 0.5d0*tdot2(ipair))*dtu +
     &                                                     r(ipair))*dtu
*
*       Ensure that c.m. step is not exceeded (H > 0 is OK).
      IF (step(i1).GT.step(i).AND.h(ipair).LT.0.0) THEN
          dtu = 0.5*dtu
          go to 75
      END IF
*
*       Initialize step & time and include possible slow-down factor.
   80 dtau(ipair) = dtu
      t0(i1) = time
      IF (imod.GT.1) THEN
          step(i1) = zmod*step(i1)
      END IF
*
*       Include factorials in force and first derivative.
      DO 85 k = 1,4
          fu(k,ipair) = 0.5d0*fu(k,ipair)
          fudot(k,ipair) = one6*fudot(k,ipair)
   85 CONTINUE
*
      RETURN
*
      END