xtrnlf.f

Go to the documentation of this file.

      SUBROUTINE xtrnlf(XI,XIDOT,FIRR,FREG,FD,FDR,KCASE)
*
*
*       External force & first derivative.
*       ----------------------------------
*
      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  xi(3),xidot(3),firr(3),freg(3),fd(3),fdr(3),
     &        xg(3),xgdot(3),fm(3),fmd(3),fs(3),fsd(3)
*
*
*       See whether to include a linearized galactic tidal force (two cases).
      IF (kz(14).LE.2) THEN
          firr(1) = firr(1) + tidal(4)*xidot(2)
          firr(2) = firr(2) - tidal(4)*xidot(1)
          fd(1) = fd(1) + tidal(4)*(firr(2) + freg(2))
          fd(2) = fd(2) - tidal(4)*(firr(1) + freg(1))
*       Add smooth part to the regular components (KCASE = 1 in REGINT).
          IF (kcase.GT.0) THEN
              freg(1) = freg(1) + tidal(1)*xi(1)
              freg(3) = freg(3) + tidal(3)*xi(3)
              fdr(1) = fdr(1) + tidal(1)*xidot(1)
              fdr(3) = fdr(3) + tidal(3)*xidot(3)
          END IF
      END IF
*
*       Consider point-mass, disk and/or logarithmic halo model.
      IF (kz(14).EQ.3.AND.kcase.GT.0) THEN
          DO 5 k = 1,3
              xg(k) = rg(k) + xi(k)
              xgdot(k) = vg(k) + xidot(k)
    5     CONTINUE
*       Employ differential instead of linearized forms for better accuracy.
          IF (gmg.GT.0.0d0) THEN
              CALL fnuc(rg,vg,fs,fsd)
              CALL fnuc(xg,xgdot,fm,fmd)
              DO 6 k = 1,3
                  freg(k) = freg(k) + (fm(k) - fs(k))
                  fdr(k) = fdr(k) + (fmd(k) - fsd(k))
    6         CONTINUE
          END IF
*
*       Check bulge force.
          IF (gmb.GT.0.0d0) THEN
              CALL fbulge(rg,vg,fs,fsd)
              CALL fbulge(xg,xgdot,fm,fmd)
              DO 10 k = 1,3
                  freg(k) = freg(k) + (fm(k) - fs(k))
                  fdr(k) = fdr(k) + (fmd(k) - fsd(k))
   10         CONTINUE
          END IF
*
*       Include Miyamoto disk for positive disk mass.
          IF (disk.GT.0.0d0) THEN
              CALL fdisk(rg,vg,fs,fsd)
              CALL fdisk(xg,xgdot,fm,fmd)
              DO 20 k = 1,3
                  freg(k) = freg(k) + (fm(k) - fs(k))
                  fdr(k) = fdr(k) + (fmd(k) - fsd(k))
   20         CONTINUE
          END IF
*
*       Check addition of logarithmic halo potential to regular force.
          IF (v02.GT.0.0d0) THEN
              CALL fhalo(rg,vg,fs,fsd)
              CALL fhalo(xg,xgdot,fm,fmd)
              DO 30 k = 1,3
                  freg(k) = freg(k) + (fm(k) - fs(k))
                  fdr(k) = fdr(k) + (fmd(k) - fsd(k))
   30         CONTINUE
          END IF
      END IF
*
*       Include optional Plummer potential in the regular force.
      IF (((kz(14).EQ.3.AND.mp.GT.0.0d0).OR.kz(14).EQ.4).
     &    and.kcase.GT.0) THEN
          ri2 = ap2
          rrdot = 0.0
          DO 40 k = 1,3
              ri2 = ri2 + xi(k)**2
              rrdot = rrdot + xi(k)*xidot(k)
   40     CONTINUE
          IF (time + toff.GT.tdelay) THEN
              ymdot = -mp0*mpdot/(1.0 + mpdot*(time+toff - tdelay))**2
*       Form alternative expression (NB! must be consistent with intgrt.f).
*             DT = TIME + TOFF - TDELAY
*             YMDOT = -MP0*MPDOT*EXP(-MPDOT*DT)
          ELSE
              ymdot = 0.0
          END IF
          zf = 1.0/ri2
          zf2 = zf**1.5
*       Absorb scaling factor in 1/R3 term as MP*ZF2 (cf. Heggie & Hut p.73).
          fmp = mp*zf2
          DO 50 k = 1,3
              freg(k) = freg(k) - xi(k)*fmp
              fdr(k) = fdr(k) - (xidot(k) - 3.0*rrdot*zf*xi(k))*fmp
              fdr(k) = fdr(k) - ymdot*zf2*xi(k)
   50     CONTINUE
      END IF
*
      RETURN
*
      END