xtrnlv.f

Go to the documentation of this file.

      SUBROUTINE xtrnlv(I1,I2)
*
*
*       External potential and virial energy.
*       -------------------------------------
*
      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)
      SAVE first
      LOGICAL first
      DATA first /.true./
*
*
      et = 0.0d0
*       Skip external potential during scaling (parameters not defined).
      IF (first.AND.time.EQ.0.0d0) THEN
          first = .false.
          go to 30
*       Treat all cases uniformly (but note Coliolis term in FIRR).
      ELSE IF (kz(14).LE.4.AND.i2.NE.i1) THEN
          vir = 0.0
          DO 5 i = i1,i2
              DO 2 k = 1,3
*       Initialize scalars for force contributions.
                  xi(k) = x(k,i)
                  xidot(k) = xdot(k,i)
                  firr(k) = 0.0
                  freg(k) = 0.0
                  fd(k) = 0.0
                  fdr(k) = 0.0
    2         CONTINUE
*       Evaluate the virial energy using definition sum {m*r*F}.
              CALL xtrnlf(xi,xidot,firr,freg,fd,fdr,1)
*       Note that FIRR(K) only contributes for #14 = 1 or 2.
              DO 4 k = 1,3
                  vir = vir + body(i)*xi(k)*(firr(k) + freg(k))
    4         CONTINUE
    5     CONTINUE
*       Quit for pure 3D tidal field (case of full N summation).
          IF (kz(14).EQ.3) go to 40
*       Note ETIDE is accumulated after each regular step (REGINT/GPUCOR).
      END IF
*
*       See whether to include a linearized galactic tidal force.
      IF (kz(14).LE.2) THEN
          DO 10 i = i1,i2
              et = et - 0.5d0*body(i)*(tidal(1)*x(1,i)**2 +
     &                                 tidal(3)*x(3,i)**2)
   10     CONTINUE
      ELSE IF (kz(14).EQ.3) THEN
*       Retain 3D potential energy expressions for future use.
          i = i1
          rg2 = 0.0
          ri2 = 0.0
          DO 20 k = 1,3
              xi(k) = x(k,i)
              rg2 = rg2 + rg(k)**2
              ri2 = ri2 + (rg(k) + xi(k))**2
   20     CONTINUE
*
*       Include galaxy point mass term for body #I in differential form.
          IF (gmg.GT.0.0d0) THEN
              et = et + gmg*(1.0/sqrt(ri2) - 1.0/sqrt(rg2))
          END IF
*
*       Check contribution from gamma/eta bulge potential.
          IF (gmb.GT.0.0d0) THEN
              IF (gam.NE.2.0d0) THEN
                  rrg = 1.0 - (1.0 + ar/sqrt(rg2))**(gam-2.0)
                  rri = 1.0 - (1.0 + ar/sqrt(ri2))**(gam-2.0)
                  et = et + gmb/(ar*(2.0-gam))*(rri-rrg)
              ELSE
                  rrg = 1.0 + ar/sqrt(rg2)
                  rri = 1.0 + ar/sqrt(ri2)
                  et = et + gmb/ar*log(rrg/rri)
              ENDIF
          END IF
*
*       Add optional Miyamoto disk potential.
          IF (disk.GT.0.0d0) THEN
              r2 = (rg(1) + xi(1))**2 + (rg(2) + xi(2))**2
              bz = sqrt(b**2 + (rg(3) + xi(3))**2)
              az = sqrt(r2 + (a + bz)**2)
              r20 = rg(1)**2 + rg(2)**2
              bz0 = sqrt(b**2 + rg(3)**2)
              az0 = sqrt(r20 + (a + bz0)**2)
              et = et + disk*(1.0/az - 1.0/az0)
          END IF
*
*       Check addition of differential logarithmic potential.
          IF (v02.GT.0.0d0) THEN
              et = et + 0.5*v02*(log(ri2) - log(rg2))
          END IF
*       Form the differential potential energy due to tides.
          et = body(i)*et
      END IF
*
*       Place sum in ETIDE and single particle contribution in HT.
   30 IF (i2.GT.i1) THEN
          etide = et
      ELSE
          ht = et
      END IF
*
   40 RETURN
*
      END