xtrnl0.f

Go to the documentation of this file.

      SUBROUTINE xtrnl0
*
*
*       External force initialization.
*       ------------------------------
*
      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)
*
*
*       Check option for cluster in circular galactic orbit.
      IF (kz(14).NE.1) go to 20
*
*       Specify Oort's constants (units of km/sec/kpc).
      a = 14.4
      b = -12.0
*       Adopt local density from Gilmore & Kuijken (solar mass/pc**3).
      rho = 0.11
*       Convert rotation constants to units of cm/sec/pc.
      a = 100.0*a
      b = 100.0*b
*
*       Specify the tidal term in star cluster units (solar mass & pc).
      tidal(1) = 4.0*a*(a - b)*(pc/gm)
*
*       Initialize the Y-component to zero.
      tidal(2) = 0.0
*
*       Specify the vertical force gradient.
      tidal(3) = -(2.0*twopi*rho + 2.0*(a - b)*(a + b)*(pc/gm))
*
      fac = 1.0e-10/(pc/gm)
      WRITE (6,5)  zmbar*zmass, fac*tidal(1), fac*tidal(3), pc/gm
    5 FORMAT (/,12x,'TOTAL MASS =',f8.1,'  TIDAL(1&3) =',1p,2e10.2,
     &              '  PC/GM =',e10.2)
*
*       Adopt twice the angular velocity for Coriolis terms.
      tidal(4) = 2.0*(a - b)*sqrt(pc/gm)
*
*       Define time scale in seconds and velocity scale in km/sec.
      tscale = sqrt(pc/gm)*pc
*
*       Convert time scale from units of seconds to million years.
      tscale = tscale/(3.15d+07*1.0d+06)
*
*       Scale to working units of RBAR in pc & ZMBAR in solar masses.
      DO 10 k = 1,3
          tidal(k) = tidal(k)*rbar**3/zmbar
   10 CONTINUE
      tidal(4) = tidal(4)*sqrt(rbar**3/zmbar)
      tscale = tscale*sqrt(rbar**3/(zmass*zmbar))
*
*       Perform an energy scaling iteration to include ETIDE (once is OK).
      e0 = -0.25
      CALL xtrnlv(1,n)
*       Form total energy based on current values and ETIDE (cf. SCALE).
      etot = zkin - pot + etide
      sx = e0/etot
*       Exclude re-scaling for positive energy or binary Plummer cluster.
      IF (etot.GE.0.0) sx = 1.0
      IF (kz(5).EQ.2) sx = 1.0
*
*       Scale coordinates & velocities to yield ETOT = -0.25.
      DO 15 i = 1,n
          DO 12 k = 1,3
              x(k,i) = x(k,i)/sx
              xdot(k,i) = xdot(k,i)*sqrt(sx)
   12     CONTINUE
   15 CONTINUE
*
*       Consider alternatives: circular point-mass orbit or 3D galaxy model.
   20 zmtot = zmass*zmbar
      IF (kz(14).EQ.2) THEN
*
*       Read galaxy mass and central distance (solar units and kpc).
          READ (5,*)  gmg, rg0
*
*       Set circular velocity in km/sec and angular velocity in cgs units.
          vg0 = 1.0d-05*sqrt(gmg/(1000.0*rg0))*sqrt(gm/pc)
          omega = 100.0*vg0/rg0
*
*       Obtain King tidal radius in pc (eq. (9) of Fukushige & Heggie, 1995).
          rt = (zmtot/(3.0*gmg))**0.3333*(1000.0*rg0)
*
          IF (rtide.GT.0.0) THEN
*       Determine RBAR (N-body units) from RT (pc) and King model (see SCALE).
              IF(kz(22).EQ.2) rbar = rt/rtide
          ELSE
              rtide = rt/rbar
          END IF
*
*       Convert from cgs to N-body units.
          omega = omega*sqrt(pc/gm)*sqrt(rbar**3/zmbar)
*
*       Specify the galactic parameters for equations of motion.
          tidal(1) = 3.0*omega**2
          tidal(2) = 0.0d0
          tidal(3) = -omega**2
          tidal(4) = 2.0*omega
          gmg = gmg/zmtot
*
*       Check re-scaling units to current RBAR (i.e. TSTAR & VSTAR).
          IF (kz(22).EQ.2) THEN
              CALL units
          END IF
*
          WRITE (6,35)  gmg, rg0, omega, rtide, rbar
*
*       Perform an energy scaling iteration to include ETIDE (once is OK).
      e0 = -0.25
      CALL xtrnlv(1,n)
*       Form total energy based on current values and ETIDE (cf. SCALE).
      etot = zkin - pot + etide
      sx = e0/etot
      IF (etot.GE.0.0) sx = 1.0
*
*       Scale coordinates & velocities to yield ETOT = -0.25.
      DO 24 i = 1,n
          DO 22 k = 1,3
              x(k,i) = x(k,i)/sx
              xdot(k,i) = xdot(k,i)*sqrt(sx)
   22     CONTINUE
   24 CONTINUE
*
*       Treat the case of 3D orbit for point-mass, bulge, disk and/or halo.
      ELSE IF (kz(14).EQ.3) THEN
*
*       Read all parameters (NB! Do not confuse with Oort's constants A, B).
          READ (5,*)  gmg, disk, a, b, vcirc, rcirc, gmb, ar, gam
*       Gamma/eta model: GMB = mass, AR = scale radius, GAM = exponent.
          READ (5,*)  (rg(k),k=1,3), (vg(k),k=1,3)
*
*       Specify planar motion from SEMI & ECC for no disk & halo if VZ = 0.
          IF (disk + vcirc + gmb.EQ.0.0.AND.vg(3).EQ.0.0d0) THEN
              rap = rg(1)
              ecc = rg(2)
              semi = rap/(1.0 + ecc)
              vg2 = gmg/(1000.0*semi)*(1.0 - ecc)/(1.0 + ecc)
              DO 25 k = 1,3
                  rg(k) = 0.0
                  vg(k) = 0.0
   25         CONTINUE
*       Initialize 2D orbit with given eccentricity at apocentre.
              rg(1) = rap
              vg(2) = 1.0d-05*sqrt(vg2)*sqrt(gm/pc)
          END IF
*
*       Convert from kpc and km/sec to N-body units.
          DO 30 k = 1,3
              rg(k) = 1000.0*rg(k)/rbar
              vg(k) = vg(k)/vstar
   30     CONTINUE
*
*       Define the angular velocity (z-component) and mass in N-body units.
          r02 = rg(1)**2 + rg(2)**2
          omega = (rg(1)*vg(2) - rg(2)*vg(1))/r02
          tidal(4) = 2.0*omega
          gmg = gmg/zmtot
          gmb = gmb/zmtot
          ar = 1000.0*ar/rbar
*
          IF (gmg.GT.0.0) THEN
              WRITE (6,35)  gmg, sqrt(r02), omega, rtide, rbar
   35         FORMAT (/,12x,'POINT-MASS MODEL:    GMG =',1p,e9.1,
     &                      '  RG =',e9.1,'  OMEGA =',e9.1,
     &                      '  RT =',0p,f6.2,'  RBAR =',f6.2)
          END IF
          IF (gmb.GT.0.0d0) THEN
              WRITE (6,36)  gmb, ar, gam
   36         FORMAT (/,12x,'GAMMA/ETA MODEL:    GMB =',1p,e9.1,
     &                      '  AR =',e9.1,'  GAM =',e9.1)
          END IF
*
*       Define disk and/or logarithmic halo parameters in N-body units.
          IF (disk.GT.0.0d0) THEN
              disk = disk/zmtot
              a = 1000.0*a/rbar
              b = 1000.0*b/rbar
              WRITE (6,40)  disk, a, b
   40         FORMAT (/,12x,'DISK MODEL:    MD =',1p,e9.1,
     &                                   '  A =',e9.1,'  B =',e9.1)
          END IF
*
*       Determine local halo velocity from total circular velocity.
          IF (vcirc.GT.0.0d0) THEN
              vcirc = vcirc/vstar
              rcirc = 1000.0*rcirc/rbar
              a2 = rcirc**2 + (a + b)**2
              v02 = vcirc**2 - (gmg/rcirc + disk*rcirc**2/a2**1.5)
*       Include any contribution from bulge potential (V2 = R*F).
              IF (gmb.GT.0.0d0) THEN
                  vb2 = gmb/rcirc*(1.0 + ar/rcirc)**(gam-3.0)
                  v02 = v02 - vb2
              END IF
              IF (v02.LT.0.0d0) THEN
                  WRITE (6,45)  v02, 0.001*rcirc*rbar
   45             FORMAT (' ',' NEGATIVE HALO VELOCITY!    V02 RCIRC ',
     &                                                     1p,2e10.2)
                  stop
              END IF
*       Specify the corresponding scale length of logarithmic halo.
              rl2 = rcirc**2*(vcirc**2 - v02)/v02
*       Define the asymptotic circular velocity due to halo.
              v02 = vcirc**2
*
*       Include table of circular velocity on unit #51 (km/sec & kpc).
              ri = 1000.0/rbar
              dr = 1000.0/rbar
              DO 60 k = 1,30
                  ri2 = ri**2
                  a2 = ri2 + (a + b)**2
                  vb2 = gmb/ri*(1.0 + ar/ri)**(gam-3.0)
                  vcirc2 = gmg/sqrt(ri2) + disk*ri2/a2**1.5 +
     &                     v02*ri2/(rl2 + ri2) + vb2
                  WRITE (52,50)  sqrt(vcirc2)*vstar, ri*rbar/1000.0
   50             FORMAT (' CIRCULAR VELOCITY:    VC R ',f7.1,f7.2)
                  ri = ri + dr
   60         CONTINUE
              CALL flush(52)
*
              a2 = r02 + (a + b)**2
              vb2 = gmb/sqrt(r02)*(1.0 + ar/sqrt(r02))**(gam-3.0)
              vcirc2 = gmg/sqrt(r02) + disk*r02/a2**1.5 +
     &                 v02*r02/(rl2 + r02) + vb2
              vcirc = sqrt(vcirc2)*vstar
              WRITE (6,62)  vcirc, sqrt(r02)/1000.0, sqrt(rl2)/1000.0
   62         FORMAT (/,12x,'CIRCULAR VELOCITY:    VC RG RL',f7.1,2f7.2)
          ELSE
              v02 = 0.0
          END IF
*
*       Initialize F & FDOT of reference frame (point-mass galaxy is OK).
          CALL gcinit
*
*       Form tidal radius from circular angular velocity (assumes apocentre).
          IF (rtide.EQ.0.0d0) rtide = (0.5/omega**2)**0.3333
*
          WRITE (6,65)  (rg(k),k=1,3), (vg(k),k=1,3), sqrt(v02)
   65     FORMAT (/,12x,'SCALED ORBIT:    RG = ',3f7.2,
     &                                 '  VG = ',3f7.1,'  V0 =',f6.1)
      END IF
*
*       Include Plummer potential for 2D and 3D (use MP = 0 if not needed).
      IF (kz(14).EQ.3.OR.kz(14).EQ.4) THEN
*       Check input for Plummer potential.
          READ (5,*)  mp, ap2, mpdot, tdelay
          WRITE (6,70)  mp, ap2, mpdot, tdelay
   70     FORMAT (/,12x,'PLUMMER POTENTIAL:    MP =',f7.3,'  AP =',f6.2,
     &                        '  MPDOT =',f7.3,'  TDELAY =',f5.1)
          mp0 = mp
          ap2 = ap2**2
*       Rescale velocities by including the Plummer & galactic virial energy.
          IF (zkin.GT.0.0d0) THEN
*       Note that QVIR = Q is saved in routine SCALE and VIR < 0 with GPU.
              CALL energy
              vir = pot - vir
              qv = sqrt(qvir*vir/zkin)
              DO 80 i = 1,n
                  DO 78 k = 1,3
                      xdot(k,i) = xdot(k,i)*qv
   78             CONTINUE
   80         CONTINUE
          END IF
          IF (rtide.EQ.0.0d0) rtide = 10.0*rscale
      ELSE
          mp = 0.0
          mpdot = 0.0
      END IF
      rtide0 = rtide
      tscale = tstar
*
*       Define tidal radius in scaled units for linearized field.
      IF (kz(14).LE.2) THEN
          rtide = (zmass/tidal(1))**0.3333
          WRITE (6,90)  (tidal(k),k=1,4), tscale, rtide
   90     FORMAT (/,12x,'TIDAL PARAMETERS:  ',1p,4e10.2,
     &                  '  TSCALE =',e9.2,'  RTIDE =',0p,f6.2,/)
      END IF
*
      RETURN
*
      END