gcint.f

Go to the documentation of this file.

      SUBROUTINE gcint
*
*
*       Integration of 3D cluster orbit.
*       --------------------------------
*
      include 'common6.h'
      common/galaxy/ gmg,rg(3),vg(3),fg(3),fgd(3),tg,
     &               omega,disk,a,b,v02,rl2,gmb,ar,gam,zdum(7)
      REAL*8  fm(3),fd(3),fs(3),fsd(3)
*
*
*       Predict coordinates and velocities to order FDOT.
      dt = time + toff - tg
*       Note: integration step may exceed STEPX (depends on block-step).
      dt2 = 0.5*dt
      dt3 = one3*dt
      rg2 = 0.0
      rgvg = 0.0
      DO 10 k = 1,3
          rg(k) = ((fgd(k)*dt3 + fg(k))*dt2 + vg(k))*dt + rg(k)
          vg(k) = (fgd(k)*dt2 + fg(k))*dt + vg(k)
          rg2 = rg2 + rg(k)**2
          rgvg = rgvg + rg(k)*vg(k)
   10 CONTINUE
*
*       Obtain force and first derivative of point-mass galaxy.
      IF (kz(14).EQ.3.AND.gmg.GT.0.0d0) THEN
          CALL fnuc(rg,vg,fm,fd)
      ELSE
          DO 15 k = 1,3
              fm(k) = 0.0
              fd(k) = 0.0
   15     CONTINUE
      END IF
*
*       Check bulge force.
      IF (gmb.GT.0.0d0) THEN
          CALL fbulge(rg,vg,fs,fsd)
          DO 20 k = 1,3
              fm(k) = fm(k) + fs(k)
              fd(k) = fd(k) + fsd(k)
   20     CONTINUE
      END IF
*
*       Include optional Miyamoto disk component (Book eq. 8.52).
      IF (disk.GT.0.0d0) THEN
          CALL fdisk(rg,vg,fs,fsd)
          DO 25 k = 1,3
              fm(k) = fm(k) + fs(k)
              fd(k) = fd(k) + fsd(k)
   25     CONTINUE
      END IF
*
*       Check addition of logarithmic galaxy potential (not linearized).
      IF (v02.GT.0.0d0) THEN
          CALL fhalo(rg,vg,fs,fsd)
          DO 30 k = 1,3
              fm(k) = fm(k) + fs(k)
              fd(k) = fd(k) + fsd(k)
   30     CONTINUE
      END IF
*
*       Set time factors for corrector.
      dt13 = one3*dt
      dtsq12 = one12*dt**2  
      tg = tg + dt
*
*       Include the Hermite corrector and update F & FDOT.
      DO 40 k = 1,3
          df = fg(k) - fm(k)
          sum = fgd(k) + fd(k)
          at3 = 2.0*df + sum*dt
          bt2 = -3.0*df - (sum + fgd(k))*dt
          rg(k) = rg(k) + (0.6*at3 + bt2)*dtsq12
          vg(k) = vg(k) + (0.75*at3 + bt2)*dt13
          fg(k) = fm(k)
          fgd(k) = fd(k)
   40 CONTINUE
*
*       Update angular velocity in case of non-circular orbit.
*     OM1 = (RG(2)*VG(3) - RG(3)*VG(2))/RG2
*     OM2 = (RG(3)*VG(1) - RG(1)*VG(3))/RG2
      om3 = (rg(1)*vg(2) - rg(2)*vg(1))/rg2
*     OMEGA2 = OM1**2 + OM2**2 + OM3**2
      omega2 = om3**2
      omega = sqrt(omega2)
      tidal(4) = 2.0*omega
*
      RETURN
*
      END