Nbody6
xtrnlf.f
Go to the documentation of this file.
1  SUBROUTINE xtrnlf(XI,XIDOT,FIRR,FREG,FD,FDR,KCASE)
2 *
3 *
4 * External force & first derivative.
5 * ----------------------------------
6 *
7  include 'common6.h'
8  common/galaxy/ gmg,rg(3),vg(3),fg(3),fgd(3),tg,
9  & omega,disk,a,b,v02,rl2,gmb,ar,gam,zdum(7)
10  REAL*8 xi(3),xidot(3),firr(3),freg(3),fd(3),fdr(3),
11  & xg(3),xgdot(3),fm(3),fmd(3),fs(3),fsd(3)
12 *
13 *
14 * See whether to include a linearized galactic tidal force (two cases).
15  IF (kz(14).LE.2) THEN
16  firr(1) = firr(1) + tidal(4)*xidot(2)
17  firr(2) = firr(2) - tidal(4)*xidot(1)
18  fd(1) = fd(1) + tidal(4)*(firr(2) + freg(2))
19  fd(2) = fd(2) - tidal(4)*(firr(1) + freg(1))
20 * Add smooth part to the regular components (KCASE = 1 in REGINT).
21  IF (kcase.GT.0) THEN
22  freg(1) = freg(1) + tidal(1)*xi(1)
23  freg(3) = freg(3) + tidal(3)*xi(3)
24  fdr(1) = fdr(1) + tidal(1)*xidot(1)
25  fdr(3) = fdr(3) + tidal(3)*xidot(3)
26  END IF
27  END IF
28 *
29 * Consider point-mass, disk and/or logarithmic halo model.
30  IF (kz(14).EQ.3.AND.kcase.GT.0) THEN
31  DO 5 k = 1,3
32  xg(k) = rg(k) + xi(k)
33  xgdot(k) = vg(k) + xidot(k)
34  5 CONTINUE
35 * Employ differential instead of linearized forms for better accuracy.
36  IF (gmg.GT.0.0d0) THEN
37  CALL fnuc(rg,vg,fs,fsd)
38  CALL fnuc(xg,xgdot,fm,fmd)
39  DO 6 k = 1,3
40  freg(k) = freg(k) + (fm(k) - fs(k))
41  fdr(k) = fdr(k) + (fmd(k) - fsd(k))
42  6 CONTINUE
43  END IF
44 *
45 * Check bulge force.
46  IF (gmb.GT.0.0d0) THEN
47  CALL fbulge(rg,vg,fs,fsd)
48  CALL fbulge(xg,xgdot,fm,fmd)
49  DO 10 k = 1,3
50  freg(k) = freg(k) + (fm(k) - fs(k))
51  fdr(k) = fdr(k) + (fmd(k) - fsd(k))
52  10 CONTINUE
53  END IF
54 *
55 * Include Miyamoto disk for positive disk mass.
56  IF (disk.GT.0.0d0) THEN
57  CALL fdisk(rg,vg,fs,fsd)
58  CALL fdisk(xg,xgdot,fm,fmd)
59  DO 20 k = 1,3
60  freg(k) = freg(k) + (fm(k) - fs(k))
61  fdr(k) = fdr(k) + (fmd(k) - fsd(k))
62  20 CONTINUE
63  END IF
64 *
65 * Check addition of logarithmic halo potential to regular force.
66  IF (v02.GT.0.0d0) THEN
67  CALL fhalo(rg,vg,fs,fsd)
68  CALL fhalo(xg,xgdot,fm,fmd)
69  DO 30 k = 1,3
70  freg(k) = freg(k) + (fm(k) - fs(k))
71  fdr(k) = fdr(k) + (fmd(k) - fsd(k))
72  30 CONTINUE
73  END IF
74  END IF
75 *
76 * Include optional Plummer potential in the regular force.
77  IF (((kz(14).EQ.3.AND.mp.GT.0.0d0).OR.kz(14).EQ.4).
78  & and.kcase.GT.0) THEN
79  ri2 = ap2
80  rrdot = 0.0
81  DO 40 k = 1,3
82  ri2 = ri2 + xi(k)**2
83  rrdot = rrdot + xi(k)*xidot(k)
84  40 CONTINUE
85  IF (time + toff.GT.tdelay) THEN
86  ymdot = -mp0*mpdot/(1.0 + mpdot*(time+toff - tdelay))**2
87 * Form alternative expression (NB! must be consistent with intgrt.f).
88 * DT = TIME + TOFF - TDELAY
89 * YMDOT = -MP0*MPDOT*EXP(-MPDOT*DT)
90  ELSE
91  ymdot = 0.0
92  END IF
93  zf = 1.0/ri2
94  zf2 = zf**1.5
95 * Absorb scaling factor in 1/R3 term as MP*ZF2 (cf. Heggie & Hut p.73).
96  fmp = mp*zf2
97  DO 50 k = 1,3
98  freg(k) = freg(k) - xi(k)*fmp
99  fdr(k) = fdr(k) - (xidot(k) - 3.0*rrdot*zf*xi(k))*fmp
100  fdr(k) = fdr(k) - ymdot*zf2*xi(k)
101  50 CONTINUE
102  END IF
103 *
104  RETURN
105 *
106  END