fcorr.f

Go to the documentation of this file.

      SUBROUTINE fcorr(I,DM,KW)
*
*
*       Total force corrections due to masss loss.
*       ------------------------------------------
*
      include 'common6.h'
      REAL*8  a(9)
      LOGICAL ikick
*
*
*       Save the velocity components and square velocity.
      vi2 = 0.0
      DO 1 k = 1,3
          a(k+6) = xdot(k,i)
          vi2 = vi2 + xdot(k,i)**2
    1 CONTINUE
*
*       Include velocity kick in case of new WD, NS, BH or massless SN.
      IF (kw.GE.10.AND.kw.LE.15) THEN
          ikick = .true.
*       Ensure WD kicks are treated consistently (cf. define.f and #25).
          IF ((kw.EQ.10.OR.kw.EQ.11).AND.kz(25).NE.1) ikick = .false.
          IF (kw.EQ.12.AND.kz(25).NE.2) ikick = .false.
*       Distinguish between single star (first time only) and binary.
          IF (i.LE.n.AND.kw.NE.kstar(i).AND.ikick) THEN
              CALL kick(i,1,kw,dm)
          ELSE IF (i.GT.n.AND.ikick) THEN
              ipair = i - n
              CALL kick(ipair,0,kw,dm)
          END IF
      ELSE
          ikick = .false.
      END IF
*
*       Define consistent c.m. variables for KS mass loss (exclude Roche).
      IF (i.GT.n.AND.kstar(i).LE.10) THEN
          i2 = 2*(i - n)
          i1 = i2 - 1
          vf2 = 0.0
          dv2 = 0.0
          bodyi = body(i)
          IF (body(i).EQ.0.0d0) bodyi = body(i1) + body(i2)
          DO 8 k = 1,3
              x(k,i) = (body(i1)*x(k,i1) + body(i2)*x(k,i2))/bodyi
              xdot(k,i) = (body(i1)*xdot(k,i1) + body(i2)*xdot(k,i2))/
     &                                                           bodyi
              x0(k,i) = x(k,i)
              x0dot(k,i) = xdot(k,i)
              vf2 = vf2 + xdot(k,i)**2
              dv2 = dv2 + (xdot(k,i) - a(k+6))**2
    8     CONTINUE
          vfac = sqrt(vf2/vi2)
      END IF
*
*       Correct potential energy, all forces & first derivatives.
      potj = 0.0d0
      DO 40 j = ifirst,ntot
          IF (j.EQ.i) go to 40
          rij2 = 0.0d0
          rijdot = 0.0d0
          rdvdot = 0.0d0
*
          DO 10 k = 1,3
              a(k) = x(k,i) - x(k,j)
              a(k+3) = a(k+6) - xdot(k,j)
              rij2 = rij2 + a(k)**2
              rijdot = rijdot + a(k)*a(k+3)
              rdvdot = rdvdot + a(k)*(xdot(k,i) - a(k+6))
   10     CONTINUE
*
          rij = sqrt(rij2)
          potj = potj + body(j)/rij
          a3 = 1.0/(rij2*rij)
          a4 = body(i)*a3
          a5 = dm*a3
          a6 = 3.0*rijdot/rij2
          a7 = 3.0*rdvdot/rij2
*
          DO 15 k = 1,3
              a(k+3) = (a(k+3) - a6*a(k))*a5
              IF (ikick) THEN
*       Include FDOT corrections due to increased velocity.
                  a(k+3) = a(k+3) + (xdot(k,i) - a(k+6))*a4
                  a(k+3) = a(k+3) - a7*a(k)*a4
              END IF
   15     CONTINUE
*
*       Use neighbour list of #J to distinguish irregular & regular terms.
          nnb1 = list(1,j) + 1
          DO 25 l = 2,nnb1
              IF (list(l,j).EQ.i) THEN
                  DO 20 k = 1,3
                      f(k,j) = f(k,j) - 0.5*a(k)*a5
                      fi(k,j) = fi(k,j) - a(k)*a5
                      fdot(k,j) = fdot(k,j) - one6*a(k+3)
                      d1(k,j) = d1(k,j) - a(k+3)
                      fidot(k,j) = fidot(k,j) - a(k+3)
   20             CONTINUE
                  go to 40
              END IF
   25     CONTINUE      
          DO 30 k = 1,3
              f(k,j) = f(k,j) - 0.5*a(k)*a5
              fr(k,j) = fr(k,j) - a(k)*a5
              fdot(k,j) = fdot(k,j) - one6*a(k+3)
              d1r(k,j) = d1r(k,j) - a(k+3)
              frdot(k,j) = frdot(k,j) - a(k+3)
   30     CONTINUE
   40 CONTINUE
*
*       Update the potential and kinetic energy loss.
      emdot = emdot - dm*potj + 0.5*dm*vi2
*
*       Modify energy loss further for c.m. body (exclude Roche cases).
      IF (i.GT.n.AND.kstar(i).LE.10) THEN
          ecdot = ecdot - 0.5*body(i)*vi2*(vfac**2 - 1.0)
      END IF
*
*       See whether linearized tidal terms should be included.
      IF (kz(14).GT.0.AND.kz(14).LT.3) THEN
          emdot = emdot - 0.5*dm*(tidal(1)*x(1,i)**2 +
     &                            tidal(3)*x(3,i)**2)
      END IF
*
*       Check optional Plummer potential.
      IF (kz(14).EQ.4.OR.kz(14).EQ.3) THEN
          ri2 = ap2
          DO 50 k = 1,3
              ri2 = ri2 + x(k,i)**2
  50      CONTINUE
          emdot = emdot - dm*mp/sqrt(ri2)
      END IF
*
*       Accumulate energy loss for conservation check (not used).
      e(12) = emdot
*
      RETURN
*
      END