-
Notifications
You must be signed in to change notification settings - Fork 24
/
mercury6_2.for
8194 lines (8032 loc) · 259 KB
/
mercury6_2.for
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
c MERCURY6_1.FOR (ErikSoft 3 May 2002)
c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
c Author: John E. Chambers
c
c Mercury is a general-purpose N-body integration package for problems in
c celestial mechanics.
c
c------------------------------------------------------------------------------
c This package contains some subroutines taken from the Swift integration
c package by H.F.Levison and M.J.Duncan (1994) Icarus, vol 108, pp18.
c Routines taken from Swift have names beginning with `drift' or `orbel'.
c
c The standard symplectic (MVS) algorithm is described in J.Widsom and
c M.Holman (1991) Astronomical Journal, vol 102, pp1528.
c
c The hybrid symplectic algorithm is described in J.E.Chambers (1999)
c Monthly Notices of the RAS, vol 304, pp793.
c
c RADAU is described in E.Everhart (1985) in ``The Dynamics of Comets:
c Their Origin and Evolution'' p185-202, eds. A.Carusi & G.B.Valsecchi,
c pub. Reidel.
c
c The Bulirsch-Stoer algorithms are described in W.H.Press et al. (1992)
c ``Numerical Recipes in Fortran'', pub. Cambridge.
c
c texadactyl_20180507.1 Handle warnings in mio_err() calls about character string lengths
c texadactyl_20180507.2 Handle warnings in mio_spl() calls about argument array bounds too small
c texadactyl_20180507.3 nnn open(..., err=nnn) is an infinite loop
c If one cannot open a file for output, let the program fail.
c texadactyl_20180514.1 Once CMAX is exceeded by variable the nclo counter,
c this will cause memory out of bounds errors.
c texadactyl_20180514.2 Return a nonzero error code to the O/S from mio_err().
c Repeat info.out diagnostic on the operator console.
c
c------------------------------------------------------------------------------
c
c Variables:
c ---------
c M = mass (in solar masses)
c XH = coordinates (x,y,z) with respect to the central body (AU)
c VH = velocities (vx,vy,vz) with respect to the central body (AU/day)
c S = spin angular momentum (solar masses AU^2/day)
c RHO = physical density (g/cm^3)
c RCEH = close-encounter limit (Hill radii)
c STAT = status (0 => alive, <>0 => to be removed)
c ID = name of the object (25 characters)
c CE = close encounter status
c NGF = (1-3) cometary non-gravitational (jet) force parameters
c " = (4) beta parameter for radiation pressure and P-R drag
c EPOCH = epoch of orbit (days)
c NBOD = current number of bodies (INCLUDING the central object)
c NBIG = " " " big bodies (ones that perturb everything else)
c TIME = current epoch (days)
c TOUT = time of next output evaluation
c TDUMP = time of next data dump
c TFUN = time of next periodic effect (e.g. next check for ejections)
c H = current integration timestep (days)
c EN(1) = initial energy of the system
c " (2) = current " " " "
c " (3) = energy change due to collisions, ejections etc.
c AM(1,2,3) = as above but for angular momentum
c
c Integration Parameters :
c ----------------------
c ALGOR = 1 -> Mixed-variable symplectic
c 2 -> Bulirsch-Stoer integrator
c 3 -> " " (conservative systems only)
c 4 -> RA15 `radau' integrator
c 10 -> Hybrid MVS/BS (democratic-heliocentric coords)
c 11 -> Close-binary hybrid (close-binary coords)
c 12 -> Wide-binary hybrid (wide-binary coords)
c
c TSTART = epoch of first required output (days)
c TSTOP = " final required output ( " )
c DTOUT = data output interval ( " )
c DTDUMP = data-dump interval ( " )
c DTFUN = interval for other periodic effects (e.g. check for ejections)
c H0 = initial integration timestep (days)
c TOL = Integrator tolerance parameter (approx. error per timestep)
c RMAX = heliocentric distance at which objects are considered ejected (AU)
c RCEN = radius of central body (AU)
c JCEN(1,2,3) = J2,J4,J6 for central body (units of RCEN^i for Ji)
c
c Options:
c OPT(1) = close-encounter option (0=stop after an encounter, 1=continue)
c OPT(2) = collision option (0=no collisions, 1=merge, 2=merge+fragment)
c OPT(3) = time style (0=days 1=Greg.date 2/3=days/years w/respect to start)
c OPT(4) = o/p precision (1,2,3 = 4,9,15 significant figures)
c OPT(5) = < Not used at present >
c OPT(6) = < Not used at present >
c OPT(7) = apply post-Newtonian correction? (0=no, 1=yes)
c OPT(8) = apply user-defined force routine mfo_user? (0=no, 1=yes)
c
c File variables :
c --------------
c OUTFILE (1) = osculating coordinates/velocities and masses
c " (2) = close encounter details
c " (3) = information file
c DUMPFILE (1) = Big-body data
c " (2) = Small-body data
c " (3) = integration parameters
c " (4) = restart file
c
c Flags :
c -----
c NGFLAG = do any bodies experience non-grav. forces?
c ( 0 = no non-grav forces)
c 1 = cometary jets only
c 2 = radiation pressure/P-R drag only
c 3 = both
c OPFLAG = integration mode (-2 = synchronising epochs)
c -1 = integrating towards start epoch
c 0 = main integration, normal output
c 1 = main integration, full output
c
c------------------------------------------------------------------------------
c
implicit none
include 'mercury.inc'
c
integer j,algor,nbod,nbig,opt(8),stat(NMAX),lmem(NMESS)
integer opflag,ngflag,ndump,nfun
real*8 m(NMAX),xh(3,NMAX),vh(3,NMAX),s(3,NMAX),rho(NMAX)
real*8 rceh(NMAX),epoch(NMAX),ngf(4,NMAX),rmax,rcen,jcen(3)
real*8 cefac,time,tstart,tstop,dtout,h0,tol,en(3),am(3)
character*25 id(NMAX)
character*80 outfile(3), dumpfile(4), mem(NMESS)
external mdt_mvs, mdt_bs1, mdt_bs2, mdt_ra15, mdt_hy
external mco_dh2h,mco_h2dh
external mco_b2h,mco_h2b,mco_h2mvs,mco_mvs2h,mco_iden
c
data opt/0,1,1,2,0,1,0,0/
c
c------------------------------------------------------------------------------
c
c Get initial conditions and integration parameters
call mio_in (time,tstart,tstop,dtout,algor,h0,tol,rmax,rcen,jcen,
% en,am,cefac,ndump,nfun,nbod,nbig,m,xh,vh,s,rho,rceh,stat,id,
% epoch,ngf,opt,opflag,ngflag,outfile,dumpfile,lmem,mem)
c
c If this is a new integration, integrate all the objects to a common epoch.
if (opflag.eq.-2) then
c texadactyl_20180507.3
C 20 open (23,file=outfile(3),status='old',access='append',err=20)
open (23,file=outfile(3),status='old',access='append')
write (23,'(/,a)') mem(55)(1:lmem(55))
write (*,'(a)') mem(55)(1:lmem(55))
call mxx_sync (time,tstart,h0,tol,jcen,nbod,nbig,m,xh,vh,s,rho,
% rceh,stat,id,epoch,ngf,opt,ngflag)
write (23,'(/,a,/)') mem(56)(1:lmem(56))
write (*,'(a)') mem(56)(1:lmem(56))
opflag = -1
close (23)
end if
c
c Main integration
if (algor.eq.1) call mal_hcon (time,tstart,tstop,dtout,algor,h0,
% tol,jcen,rcen,rmax,en,am,cefac,ndump,nfun,nbod,nbig,m,xh,vh,s,
% rho,rceh,stat,id,ngf,opt,opflag,ngflag,outfile,dumpfile,mem,
% lmem,mdt_mvs,mco_h2mvs,mco_mvs2h)
c
if (algor.eq.9) call mal_hcon (time,tstart,tstop,dtout,algor,h0,
% tol,jcen,rcen,rmax,en,am,cefac,ndump,nfun,nbod,nbig,m,xh,vh,s,
% rho,rceh,stat,id,ngf,opt,opflag,ngflag,outfile,dumpfile,mem,
% lmem,mdt_mvs,mco_iden,mco_iden)
c
if (algor.eq.2) call mal_hvar (time,tstart,tstop,dtout,algor,h0,
% tol,jcen,rcen,rmax,en,am,cefac,ndump,nfun,nbod,nbig,m,xh,vh,s,
% rho,rceh,stat,id,ngf,opt,opflag,ngflag,outfile,dumpfile,mem,
% lmem,mdt_bs1)
c
if (algor.eq.3) call mal_hvar (time,tstart,tstop,dtout,algor,h0,
% tol,jcen,rcen,rmax,en,am,cefac,ndump,nfun,nbod,nbig,m,xh,vh,s,
% rho,rceh,stat,id,ngf,opt,opflag,ngflag,outfile,dumpfile,mem,
% lmem,mdt_bs2)
c
if (algor.eq.4) call mal_hvar (time,tstart,tstop,dtout,algor,h0,
% tol,jcen,rcen,rmax,en,am,cefac,ndump,nfun,nbod,nbig,m,xh,vh,s,
% rho,rceh,stat,id,ngf,opt,opflag,ngflag,outfile,dumpfile,mem,
% lmem,mdt_ra15)
c
if (algor.eq.10) call mal_hcon (time,tstart,tstop,dtout,algor,h0,
% tol,jcen,rcen,rmax,en,am,cefac,ndump,nfun,nbod,nbig,m,xh,vh,s,
% rho,rceh,stat,id,ngf,opt,opflag,ngflag,outfile,dumpfile,mem,
% lmem,mdt_hy,mco_h2dh,mco_dh2h)
c
c Do a final data dump
do j = 2, nbod
epoch(j) = time
end do
call mio_dump (time,tstart,tstop,dtout,algor,h0,tol,jcen,rcen,
% rmax,en,am,cefac,ndump,nfun,nbod,nbig,m,xh,vh,s,rho,rceh,stat,
% id,ngf,epoch,opt,opflag,dumpfile,mem,lmem)
c
c Calculate and record the overall change in energy and ang. momentum
c texadactyl_20180507.3
c 50 open (23, file=outfile(3), status='old', access='append',
c % err=50)
open (23, file=outfile(3), status='old', access='append')
write (23,'(/,a)') mem(57)(1:lmem(57))
call mxx_en (jcen,nbod,nbig,m,xh,vh,s,en(2),am(2))
c
write (23,231) mem(58)(1:lmem(58)),
% abs((en(2) + en(3) - en(1)) / en(1))
write (23,232) mem(59)(1:lmem(59)),
% abs((am(2) + am(3) - am(1)) / am(1))
write (23,231) mem(60)(1:lmem(60)), abs(en(3) / en(1))
write (23,232) mem(61)(1:lmem(61)), abs(am(3) / am(1))
close (23)
write (*,'(a)') mem(57)(1:lmem(57))
c
c------------------------------------------------------------------------------
c
231 format (/,a,1p1e12.5)
232 format (a,1p1e12.5)
stop
end
c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
c MFO_USER.FOR (ErikSoft 2 March 2001)
c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
c Author: John E. Chambers
c
c Applies an arbitrary force, defined by the user.
c
c If using with the symplectic algorithm MAL_MVS, the force should be
c small compared with the force from the central object.
c If using with the conservative Bulirsch-Stoer algorithm MAL_BS2, the
c force should not be a function of the velocities.
c
c N.B. All coordinates and velocities must be with respect to central body
c ===
c------------------------------------------------------------------------------
c
subroutine mfo_user (time,jcen,nbod,nbig,m,x,v,a)
c
implicit none
include 'mercury.inc'
c
c Input/Output
integer nbod, nbig
real*8 time,jcen(3),m(nbod),x(3,nbod),v(3,nbod),a(3,nbod)
c
c Local
integer j
c
c------------------------------------------------------------------------------
c
do j = 1, nbod
a(1,j) = 0.d0
a(2,j) = 0.d0
a(3,j) = 0.d0
end do
c
c------------------------------------------------------------------------------
c
return
end
c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
c MAL_HVAR.FOR (ErikSoft 4 March 2001)
c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
c Author: John E. Chambers
c
c Does an integration using a variable-timestep integration algorithm. The
c particular integrator routine is ONESTEP and the algorithm must use
c coordinates with respect to the central body.
c
c N.B. This routine is also called by the synchronisation routine mxx_sync,
c === in which case OPFLAG = -2. Beware when making changes involving OPFLAG.
c
c------------------------------------------------------------------------------
c
subroutine mal_hvar (time,tstart,tstop,dtout,algor,h0,tol,jcen,
% rcen,rmax,en,am,cefac,ndump,nfun,nbod,nbig,m,xh,vh,s,rho,rceh,
% stat,id,ngf,opt,opflag,ngflag,outfile,dumpfile,mem,lmem,onestep)
c
implicit none
include 'mercury.inc'
c
c Input/Output
integer algor,nbod,nbig,stat(nbod),opt(8),opflag,ngflag,ndump,nfun
integer lmem(NMESS)
real*8 time,tstart,tstop,dtout,h0,tol,jcen(3),rcen,rmax
real*8 en(3),am(3),cefac,m(nbod),xh(3,nbod),vh(3,nbod)
real*8 s(3,nbod),rho(nbod),rceh(nbod),ngf(4,nbod)
character*25 id(nbod)
character*80 outfile(3),dumpfile(4),mem(NMESS)
c
c Local
integer i,j,k,n,itmp,nhit,ihit(CMAX),jhit(CMAX),chit(CMAX)
integer dtflag,ejflag,nowflag,stopflag,nstored,ce(NMAX)
integer nclo,iclo(CMAX),jclo(CMAX),nce,ice(NMAX),jce(NMAX)
real*8 tmp0,h,hdid,tout,tdump,tfun,tlog,tsmall,dtdump,dtfun
real*8 thit(CMAX),dhit(CMAX),thit1,x0(3,NMAX),v0(3,NMAX)
real*8 rce(NMAX),rphys(NMAX),rcrit(NMAX),a(NMAX)
real*8 dclo(CMAX),tclo(CMAX),epoch(NMAX)
real*8 ixvclo(6,CMAX),jxvclo(6,CMAX)
external mfo_all,onestep
c
c------------------------------------------------------------------------------
c
c Initialize variables. DTFLAG = 0 implies first ever call to ONESTEP
dtout = abs(dtout)
dtdump = abs(h0) * ndump
dtfun = abs(h0) * nfun
dtflag = 0
nstored = 0
tsmall = h0 * 1.d-8
h = h0
do j = 2, nbod
ce(j) = 0.d0
end do
c
c Calculate close-encounter limits and physical radii for massive bodies
call mce_init (tstart,algor,h0,jcen,rcen,rmax,cefac,nbod,nbig,
% m,xh,vh,s,rho,rceh,rphys,rce,rcrit,id,opt,outfile(2),1)
c
c Set up time of next output, times of previous dump, log and periodic effect
if (opflag.eq.-1) then
tout = tstart
else
n = int (abs (time - tstart) / dtout) + 1
tout = tstart + dtout * sign (dble(n), tstop - tstart)
if ((tstop - tstart)*(tout - tstop).gt.0d0) tout = tstop
end if
tdump = time
tfun = time
tlog = time
c
c------------------------------------------------------------------------------
c
c MAIN LOOP STARTS HERE
c
100 continue
c
c Is it time for output ?
if (abs(tout-time).lt.abs(tsmall).and.opflag.ge.-1) then
c
c Beware: the integration may change direction at this point!!!!
if (opflag.eq.-1) dtflag = 0
c
c Output data for all bodies
call mio_out (time,jcen,rcen,rmax,nbod,nbig,m,xh,vh,s,rho,
% stat,id,opt,opflag,algor,outfile(1))
call mio_ce (time,tstart,rcen,rmax,nbod,nbig,m,stat,id,
% 0,iclo,jclo,opt,stopflag,tclo,dclo,ixvclo,jxvclo,mem,lmem,
% outfile,nstored,0)
tmp0 = tstop - tout
tout = tout + sign( min( abs(tmp0), abs(dtout) ), tmp0 )
c
c Update the data dump files
do j = 2, nbod
epoch(j) = time
end do
call mio_dump (time,tstart,tstop,dtout,algor,h,tol,jcen,rcen,
% rmax,en,am,cefac,ndump,nfun,nbod,nbig,m,xh,vh,s,rho,rceh,stat,
% id,ngf,epoch,opt,opflag,dumpfile,mem,lmem)
tdump = time
end if
c
c If integration has finished return to the main part of programme
if (abs(tstop-time).le.abs(tsmall).and.opflag.ne.-1) return
c
c Set the timestep
if (opflag.eq.-1) tmp0 = tstart - time
if (opflag.eq.-2) tmp0 = tstop - time
if (opflag.ge.0) tmp0 = tout - time
h = sign ( max( min( abs(tmp0), abs(h) ), tsmall), tmp0 )
c
c Save the current coordinates and velocities
call mco_iden (time,jcen,nbod,nbig,h,m,xh,vh,x0,v0,ngf,ngflag,opt)
c
c Advance one timestep
call onestep (time,h,hdid,tol,jcen,nbod,nbig,m,xh,vh,s,rphys,
% rcrit,ngf,stat,dtflag,ngflag,opt,nce,ice,jce,mfo_all)
time = time + hdid
c
c Check if close encounters or collisions occurred
nclo = 0
call mce_stat (time,h,rcen,nbod,nbig,m,x0,v0,xh,vh,rce,rphys,
% nclo,iclo,jclo,dclo,tclo,ixvclo,jxvclo,nhit,ihit,jhit,
% chit,dhit,thit,thit1,nowflag,stat,outfile(3),mem,lmem)
c
c------------------------------------------------------------------------------
c
c CLOSE ENCOUNTERS
c
c If encounter minima occurred, output details and decide whether to stop
if (nclo.gt.0.and.opflag.ge.-1) then
itmp = 1
if (nhit.ne.0) itmp = 0
call mio_ce (time,tstart,rcen,rmax,nbod,nbig,m,stat,id,nclo,
% iclo,jclo,opt,stopflag,tclo,dclo,ixvclo,jxvclo,mem,lmem,
% outfile,nstored,itmp)
if (stopflag.eq.1) return
end if
c
c------------------------------------------------------------------------------
c
c COLLISIONS
c
c If a collision occurred, output details and resolve the collision
if (nhit.gt.0.and.opt(2).ne.0) then
do k = 1, nhit
if (chit(k).eq.1) then
i = ihit(k)
j = jhit(k)
call mce_coll (thit(k),tstart,en(3),jcen,i,j,nbod,nbig,m,xh,
% vh,s,rphys,stat,id,opt,mem,lmem,outfile(3))
end if
end do
c
c Remove lost objects, reset flags and recompute Hill and physical radii
call mxx_elim (nbod,nbig,m,xh,vh,s,rho,rceh,rcrit,ngf,stat,
% id,mem,lmem,outfile(3),itmp)
dtflag = 1
if (opflag.ge.0) opflag = 1
call mce_init (tstart,algor,h0,jcen,rcen,rmax,cefac,nbod,nbig,
% m,xh,vh,s,rho,rceh,rphys,rce,rcrit,id,opt,outfile(2),1)
end if
c
c------------------------------------------------------------------------------
c
c COLLISIONS WITH CENTRAL BODY
c
c Check for collisions
call mce_cent (time,hdid,rcen,jcen,2,nbod,nbig,m,x0,v0,xh,vh,nhit,
% jhit,thit,dhit,algor,ngf,ngflag)
c
c Resolve the collisions
if (nhit.gt.0) then
do k = 1, nhit
i = 1
j = jhit(k)
call mce_coll (thit(k),tstart,en(3),jcen,i,j,nbod,nbig,m,xh,
% vh,s,rphys,stat,id,opt,mem,lmem,outfile(3))
end do
c
c Remove lost objects, reset flags and recompute Hill and physical radii
call mxx_elim (nbod,nbig,m,xh,vh,s,rho,rceh,rcrit,ngf,stat,
% id,mem,lmem,outfile(3),itmp)
dtflag = 1
if (opflag.ge.0) opflag = 1
call mce_init (tstart,algor,h0,jcen,rcen,rmax,cefac,nbod,nbig,
% m,xh,vh,s,rho,rceh,rphys,rce,rcrit,id,opt,outfile(2),0)
end if
c
c------------------------------------------------------------------------------
c
c DATA DUMP AND PROGRESS REPORT
c
c Do the data dump
if (abs(time-tdump).ge.abs(dtdump).and.opflag.ge.-1) then
do j = 2, nbod
epoch(j) = time
end do
call mio_ce (time,tstart,rcen,rmax,nbod,nbig,m,stat,id,
% 0,iclo,jclo,opt,stopflag,tclo,dclo,ixvclo,jxvclo,mem,lmem,
% outfile,nstored,0)
call mio_dump (time,tstart,tstop,dtout,algor,h,tol,jcen,rcen,
% rmax,en,am,cefac,ndump,nfun,nbod,nbig,m,xh,vh,s,rho,rceh,stat,
% id,ngf,epoch,opt,opflag,dumpfile,mem,lmem)
tdump = time
end if
c
c Write a progress report to the log file
if (abs(time-tlog).ge.abs(dtdump).and.opflag.ge.0) then
call mxx_en (jcen,nbod,nbig,m,xh,vh,s,en(2),am(2))
call mio_log (time,tstart,en,am,opt,mem,lmem)
tlog = time
end if
c
c------------------------------------------------------------------------------
c
c CHECK FOR EJECTIONS AND DO OTHER PERIODIC EFFECTS
c
if (abs(time-tfun).ge.abs(dtfun).and.opflag.ge.-1) then
c
c Recompute close encounter limits, to allow for changes in Hill radii
call mce_hill (nbod,m,xh,vh,rce,a)
do j = 2, nbod
rce(j) = rce(j) * rceh(j)
end do
c
c Check for ejections
call mxx_ejec (time,tstart,rmax,en,am,jcen,2,nbod,nbig,m,xh,vh,
% s,stat,id,opt,ejflag,outfile(3),mem,lmem)
c
c Remove lost objects, reset flags and recompute Hill and physical radii
if (ejflag.ne.0) then
call mxx_elim (nbod,nbig,m,xh,vh,s,rho,rceh,rcrit,ngf,stat,
% id,mem,lmem,outfile(3),itmp)
dtflag = 1
if (opflag.ge.0) opflag = 1
call mce_init (tstart,algor,h0,jcen,rcen,rmax,cefac,nbod,nbig,
% m,xh,vh,s,rho,rceh,rphys,rce,rcrit,id,opt,outfile(2),0)
end if
tfun = time
end if
c
c Go on to the next time step
goto 100
c
c------------------------------------------------------------------------------
c
end
c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
c MAL_HCON.FOR (ErikSoft 28 March 2001)
c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
c Author: John E. Chambers
c
c Does an integration using an integrator with a constant stepsize H.
c Input and output to this routine use coordinates XH, and velocities VH,
c with respect to the central body, but the integration algorithm uses
c its own internal coordinates X, and velocities V.
c
c The programme uses the transformation routines COORD and BCOORD to change
c to and from the internal coordinates, respectively.
c
c------------------------------------------------------------------------------
c
subroutine mal_hcon (time,tstart,tstop,dtout,algor,h0,tol,jcen,
% rcen,rmax,en,am,cefac,ndump,nfun,nbod,nbig,m,xh,vh,s,rho,rceh,
% stat,id,ngf,opt,opflag,ngflag,outfile,dumpfile,mem,lmem,onestep,
% coord,bcoord)
c
implicit none
include 'mercury.inc'
c
c Input/Output
integer algor,nbod,nbig,stat(nbod),opt(8),opflag,ngflag
integer lmem(NMESS),ndump,nfun
real*8 time,tstart,tstop,dtout,h0,tol,jcen(3),rcen,rmax
real*8 en(3),am(3),cefac,m(nbod),xh(3,nbod),vh(3,nbod)
real*8 s(3,nbod),rho(nbod),rceh(nbod),ngf(4,nbod)
character*25 id(nbod)
character*80 outfile(3),dumpfile(4),mem(NMESS)
c
c Local
integer i,j,k,n,itmp,nclo,nhit,jhit(CMAX),iclo(CMAX),jclo(CMAX)
integer dtflag,ejflag,stopflag,colflag,nstored
real*8 x(3,NMAX),v(3,NMAX),xh0(3,NMAX),vh0(3,NMAX)
real*8 rce(NMAX),rphys(NMAX),rcrit(NMAX),epoch(NMAX)
real*8 hby2,tout,tmp0,tdump,tfun,tlog,dtdump,dtfun
real*8 dclo(CMAX),tclo(CMAX),dhit(CMAX),thit(CMAX)
real*8 ixvclo(6,CMAX),jxvclo(6,CMAX),a(NMAX)
external onestep,coord,bcoord
c
c------------------------------------------------------------------------------
c
c Initialize variables. DTFLAG = 0/2: first call ever/normal
dtout = abs(dtout)
dtdump = abs(h0) * ndump
dtfun = abs(h0) * nfun
dtflag = 0
nstored = 0
hby2 = 0.500001d0 * abs(h0)
c
c Calculate close-encounter limits and physical radii
call mce_init (tstart,algor,h0,jcen,rcen,rmax,cefac,nbod,nbig,
% m,xh,vh,s,rho,rceh,rphys,rce,rcrit,id,opt,outfile(2),1)
c
c Set up time of next output, times of previous dump, log and periodic effect
if (opflag.eq.-1) then
tout = tstart
else
n = int (abs (time-tstart) / dtout) + 1
tout = tstart + dtout * sign (dble(n), tstop - tstart)
if ((tstop-tstart)*(tout-tstop).gt.0) tout = tstop
end if
tdump = time
tfun = time
tlog = time
c
c Convert to internal coordinates and velocities
call coord (time,jcen,nbod,nbig,h0,m,xh,vh,x,v,ngf,ngflag,opt)
c
c------------------------------------------------------------------------------
c
c MAIN LOOP STARTS HERE
c
100 continue
c
c Is it time for output ?
if (abs(tout-time).le.hby2.and.opflag.ge.-1) then
c
c Beware: the integration may change direction at this point!!!!
if (opflag.eq.-1.and.dtflag.ne.0) dtflag = 1
c
c Convert to heliocentric coordinates and output data for all bodies
call bcoord (time,jcen,nbod,nbig,h0,m,x,v,xh,vh,ngf,ngflag,opt)
call mio_out (time,jcen,rcen,rmax,nbod,nbig,m,xh,vh,s,rho,
% stat,id,opt,opflag,algor,outfile(1))
call mio_ce (time,tstart,rcen,rmax,nbod,nbig,m,stat,id,
% 0,iclo,jclo,opt,stopflag,tclo,dclo,ixvclo,jxvclo,mem,lmem,
% outfile,nstored,0)
tmp0 = tstop - tout
tout = tout + sign( min( abs(tmp0), abs(dtout) ), tmp0 )
c
c Update the data dump files
do j = 2, nbod
epoch(j) = time
end do
call mio_dump (time,tstart,tstop,dtout,algor,h0,tol,jcen,rcen,
% rmax,en,am,cefac,ndump,nfun,nbod,nbig,m,xh,vh,s,rho,rceh,stat,
% id,ngf,epoch,opt,opflag,dumpfile,mem,lmem)
tdump = time
end if
c
c If integration has finished, convert to heliocentric coords and return
if (abs(tstop-time).le.hby2.and.opflag.ge.0) then
call bcoord (time,jcen,nbod,nbig,h0,m,x,v,xh,vh,ngf,ngflag,opt)
return
end if
c
c Make sure the integration is heading in the right direction
150 continue
tmp0 = tstop - time
if (opflag.eq.-1) tmp0 = tstart - time
h0 = sign (h0, tmp0)
c
c Save the current heliocentric coordinates and velocities
if (algor.eq.1) then
call mco_iden (time,jcen,nbod,nbig,h0,m,x,v,xh0,vh0,ngf,ngflag,
% opt)
else
call bcoord(time,jcen,nbod,nbig,h0,m,x,v,xh0,vh0,ngf,ngflag,opt)
end if
call onestep (time,tstart,h0,tol,rmax,en,am,jcen,rcen,nbod,nbig,
% m,x,v,s,rphys,rcrit,rce,stat,id,ngf,algor,opt,dtflag,ngflag,
% opflag,colflag,nclo,iclo,jclo,dclo,tclo,ixvclo,jxvclo,outfile,
% mem,lmem)
time = time + h0
c
c------------------------------------------------------------------------------
c
c CLOSE ENCOUNTERS
c
c If encounter minima occurred, output details and decide whether to stop
if (nclo.gt.0.and.opflag.ge.-1) then
itmp = 1
if (colflag.ne.0) itmp = 0
call mio_ce (time,tstart,rcen,rmax,nbod,nbig,m,stat,id,nclo,
% iclo,jclo,opt,stopflag,tclo,dclo,ixvclo,jxvclo,mem,lmem,
% outfile,nstored,itmp)
if (stopflag.eq.1) return
end if
c
c------------------------------------------------------------------------------
c
c COLLISIONS
c
c If collisions occurred, output details and remove lost objects
if (colflag.ne.0) then
c
c Reindex the surviving objects
call bcoord (time,jcen,nbod,nbig,h0,m,x,v,xh,vh,ngf,ngflag,opt)
call mxx_elim (nbod,nbig,m,xh,vh,s,rho,rceh,rcrit,ngf,stat,
% id,mem,lmem,outfile(3),itmp)
c
c Reset flags, and calculate new Hill radii and physical radii
dtflag = 1
if (opflag.ge.0) opflag = 1
call mce_init (tstart,algor,h0,jcen,rcen,rmax,cefac,nbod,nbig,
% m,xh,vh,s,rho,rceh,rphys,rce,rcrit,id,opt,outfile(2),1)
call coord (time,jcen,nbod,nbig,h0,m,xh,vh,x,v,ngf,ngflag,opt)
end if
c
c------------------------------------------------------------------------------
c
c COLLISIONS WITH CENTRAL BODY
c
c Check for collisions with the central body
if (algor.eq.1) then
call mco_iden(time,jcen,nbod,nbig,h0,m,x,v,xh,vh,ngf,ngflag,opt)
else
call bcoord (time,jcen,nbod,nbig,h0,m,x,v,xh,vh,ngf,ngflag,opt)
end if
itmp = 2
if (algor.eq.11.or.algor.eq.12) itmp = 3
call mce_cent (time,h0,rcen,jcen,itmp,nbod,nbig,m,xh0,vh0,xh,vh,
% nhit,jhit,thit,dhit,algor,ngf,ngflag)
c
c If something hit the central body, restore the coords prior to this step
if (nhit.gt.0) then
call mco_iden (time,jcen,nbod,nbig,h0,m,xh0,vh0,xh,vh,ngf,
% ngflag,opt)
time = time - h0
c
c Merge the object(s) with the central body
do k = 1, nhit
i = 1
j = jhit(k)
call mce_coll (thit(k),tstart,en(3),jcen,i,j,nbod,nbig,m,xh,
% vh,s,rphys,stat,id,opt,mem,lmem,outfile(3))
end do
c
c Remove lost objects, reset flags and recompute Hill and physical radii
call mxx_elim (nbod,nbig,m,xh,vh,s,rho,rceh,rcrit,ngf,stat,
% id,mem,lmem,outfile(3),itmp)
if (opflag.ge.0) opflag = 1
dtflag = 1
call mce_init (tstart,algor,h0,jcen,rcen,rmax,cefac,nbod,nbig,
% m,xh,vh,s,rho,rceh,rphys,rce,rcrit,id,opt,outfile(2),0)
if (algor.eq.1) then
call mco_iden (time,jcen,nbod,nbig,h0,m,xh,vh,x,v,ngf,ngflag,
% opt)
else
call coord (time,jcen,nbod,nbig,h0,m,xh,vh,x,v,ngf,ngflag,opt)
end if
c
c Redo that integration time step
goto 150
end if
c
c------------------------------------------------------------------------------
c
c DATA DUMP AND PROGRESS REPORT
c
c Convert to heliocentric coords and do the data dump
if (abs(time-tdump).ge.abs(dtdump).and.opflag.ge.-1) then
call bcoord (time,jcen,nbod,nbig,h0,m,x,v,xh,vh,ngf,ngflag,opt)
do j = 2, nbod
epoch(j) = time
end do
call mio_ce (time,tstart,rcen,rmax,nbod,nbig,m,stat,id,
% 0,iclo,jclo,opt,stopflag,tclo,dclo,ixvclo,jxvclo,mem,lmem,
% outfile,nstored,0)
call mio_dump (time,tstart,tstop,dtout,algor,h0,tol,jcen,rcen,
% rmax,en,am,cefac,ndump,nfun,nbod,nbig,m,xh,vh,s,rho,rceh,stat,
% id,ngf,epoch,opt,opflag,dumpfile,mem,lmem)
tdump = time
end if
c
c Convert to heliocentric coords and write a progress report to the log file
if (abs(time-tlog).ge.abs(dtdump).and.opflag.ge.0) then
call bcoord (time,jcen,nbod,nbig,h0,m,x,v,xh,vh,ngf,ngflag,opt)
call mxx_en (jcen,nbod,nbig,m,xh,vh,s,en(2),am(2))
call mio_log (time,tstart,en,am,opt,mem,lmem)
tlog = time
end if
c
c------------------------------------------------------------------------------
c
c CHECK FOR EJECTIONS AND DO OTHER PERIODIC EFFECTS
c
if (abs(time-tfun).ge.abs(dtfun).and.opflag.ge.-1) then
if (algor.eq.1) then
call mco_iden (time,jcen,nbod,nbig,h0,m,x,v,xh,vh,ngf,ngflag,
% opt)
else
call bcoord(time,jcen,nbod,nbig,h0,m,x,v,xh,vh,ngf,ngflag,opt)
end if
c
c Recompute close encounter limits, to allow for changes in Hill radii
call mce_hill (nbod,m,xh,vh,rce,a)
do j = 2, nbod
rce(j) = rce(j) * rceh(j)
end do
c
c Check for ejections
itmp = 2
if (algor.eq.11.or.algor.eq.12) itmp = 3
call mxx_ejec (time,tstart,rmax,en,am,jcen,itmp,nbod,nbig,m,xh,
% vh,s,stat,id,opt,ejflag,outfile(3),mem,lmem)
c
c Remove ejected objects, reset flags, calculate new Hill and physical radii
if (ejflag.ne.0) then
call mxx_elim (nbod,nbig,m,xh,vh,s,rho,rceh,rcrit,ngf,stat,
% id,mem,lmem,outfile(3),itmp)
if (opflag.ge.0) opflag = 1
dtflag = 1
call mce_init (tstart,algor,h0,jcen,rcen,rmax,cefac,nbod,nbig,
% m,xh,vh,s,rho,rceh,rphys,rce,rcrit,id,opt,outfile(2),0)
if (algor.eq.1) then
call mco_iden (time,jcen,nbod,nbig,h0,m,xh,vh,x,v,ngf,
% ngflag,opt)
else
call coord (time,jcen,nbod,nbig,h0,m,xh,vh,x,v,ngf,ngflag,
% opt)
end if
end if
tfun = time
end if
c
c Go on to the next time step
goto 100
c
c------------------------------------------------------------------------------
c
end
c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
c MCE_BOX.FOR (ErikSoft 30 September 2000)
c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
c Author: John E. Chambers
c
c Given initial and final coordinates and velocities, the routine returns
c the X and Y coordinates of a box bounding the motion in between the
c end points.
c
c If the X or Y velocity changes sign, the routine performs a quadratic
c interpolation to estimate the corresponding extreme value of X or Y.
c
c------------------------------------------------------------------------------
c
subroutine mce_box (nbod,h,x0,v0,x1,v1,xmin,xmax,ymin,ymax)
c
implicit none
include 'mercury.inc'
c
c Input/Output
integer nbod
real*8 h,x0(3,nbod), x1(3,nbod), v0(3,nbod),v1(3,nbod)
real*8 xmin(nbod), xmax(nbod), ymin(nbod),ymax(nbod)
c
c Local
integer j
real*8 temp
c
c------------------------------------------------------------------------------
c
do j = 2, nbod
xmin(j) = min (x0(1,j), x1(1,j))
xmax(j) = max (x0(1,j), x1(1,j))
ymin(j) = min (x0(2,j), x1(2,j))
ymax(j) = max (x0(2,j), x1(2,j))
c
c If velocity changes sign, do an interpolation
if ((v0(1,j).lt.0d0.and.v1(1,j).gt.0d0).or.
% (v0(1,j).gt.0d0.and.v1(1,j).lt.0d0)) then
temp = (v0(1,j)*x1(1,j) - v1(1,j)*x0(1,j)
% - .5d0*h*v0(1,j)*v1(1,j)) / (v0(1,j) - v1(1,j))
xmin(j) = min (xmin(j),temp)
xmax(j) = max (xmax(j),temp)
end if
c
if ((v0(2,j).lt.0d0.and.v1(2,j).gt.0d0).or.
% (v0(2,j).gt.0d0.and.v1(2,j).lt.0d0)) then
temp = (v0(2,j)*x1(2,j) - v1(2,j)*x0(2,j)
% - .5d0*h*v0(2,j)*v1(2,j)) / (v0(2,j) - v1(2,j))
ymin(j) = min (ymin(j),temp)
ymax(j) = max (ymax(j),temp)
end if
end do
c
c------------------------------------------------------------------------------
c
return
end
c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
c MCE_CENT.FOR (ErikSoft 4 March 2001)
c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
c Author: John E. Chambers
c
c Checks all objects with index I >= I0, to see if they have had a collision
c with the central body in a time interval H, when given the initial and
c final coordinates and velocities. The routine uses cubic interpolation
c to estimate the minimum separations.
c
c N.B. All coordinates & velocities must be with respect to the central body!!
c ===
c------------------------------------------------------------------------------
c
subroutine mce_cent (time,h,rcen,jcen,i0,nbod,nbig,m,x0,v0,x1,v1,
% nhit,jhit,thit,dhit,algor,ngf,ngflag)
c
implicit none
include 'mercury.inc'
c
c Input/Output
integer i0, nbod, nbig, nhit, jhit(CMAX), algor, ngflag
real*8 time,h,rcen,jcen(3),m(nbod),x0(3,nbod),v0(3,nbod)
real*8 x1(3,nbod),v1(3,nbod),thit(CMAX),dhit(CMAX),ngf(4,nbod)
c
c Local
integer j
real*8 rcen2,mco_acsh,a,q,u0,uhit,m0,mhit,mm,r0,mcen
real*8 hx,hy,hz,h2,p,rr0,rr1,rv0,rv1,temp,e,v2
real*8 xu0(3,NMAX),xu1(3,NMAX),vu0(3,NMAX),vu1(3,NMAX)
c
c------------------------------------------------------------------------------
c
if (i0.le.0) i0 = 2
nhit = 0
rcen2 = rcen * rcen
mcen = m(1)
c
c If using close-binary code, convert to coords with respect to the binary
c if (algor.eq.11) then
c mcen = m(1) + m(2)
c call mco_h2ub (temp,jcen,nbod,nbig,h,m,x0,v0,xu0,vu0,ngf,ngflag)
c call mco_h2ub (temp,jcen,nbod,nbig,h,m,x1,v1,xu1,vu1,ngf,ngflag)
c end if
c
c Check for collisions with the central body
do j = i0, nbod
if (algor.eq.11) then
rr0 = xu0(1,j)*xu0(1,j) + xu0(2,j)*xu0(2,j) +xu0(3,j)*xu0(3,j)
rr1 = xu1(1,j)*xu1(1,j) + xu1(2,j)*xu1(2,j) +xu1(3,j)*xu1(3,j)
rv0 = vu0(1,j)*xu0(1,j) + vu0(2,j)*xu0(2,j) +vu0(3,j)*xu0(3,j)
rv1 = vu1(1,j)*xu1(1,j) + vu1(2,j)*xu1(2,j) +vu1(3,j)*xu1(3,j)
else
rr0 = x0(1,j)*x0(1,j) + x0(2,j)*x0(2,j) + x0(3,j)*x0(3,j)
rr1 = x1(1,j)*x1(1,j) + x1(2,j)*x1(2,j) + x1(3,j)*x1(3,j)
rv0 = v0(1,j)*x0(1,j) + v0(2,j)*x0(2,j) + v0(3,j)*x0(3,j)
rv1 = v1(1,j)*x1(1,j) + v1(2,j)*x1(2,j) + v1(3,j)*x1(3,j)
end if
c
c If inside the central body, or passing through pericentre, use 2-body approx.
if ((rv0*h.le.0d0.and.rv1*h.ge.0d0)
% .or.min(rr0,rr1).le.rcen2) then
if (algor.eq.11) then
hx = xu0(2,j) * vu0(3,j) - xu0(3,j) * vu0(2,j)
hy = xu0(3,j) * vu0(1,j) - xu0(1,j) * vu0(3,j)
hz = xu0(1,j) * vu0(2,j) - xu0(2,j) * vu0(1,j)
v2 = vu0(1,j)*vu0(1,j) +vu0(2,j)*vu0(2,j) +vu0(3,j)*vu0(3,j)
else
hx = x0(2,j) * v0(3,j) - x0(3,j) * v0(2,j)
hy = x0(3,j) * v0(1,j) - x0(1,j) * v0(3,j)
hz = x0(1,j) * v0(2,j) - x0(2,j) * v0(1,j)
v2 = v0(1,j)*v0(1,j) + v0(2,j)*v0(2,j) + v0(3,j)*v0(3,j)
end if
h2 = hx*hx + hy*hy + hz*hz
p = h2 / (mcen + m(j))
r0 = sqrt(rr0)
temp = 1.d0 + p*(v2/(mcen + m(j)) - 2.d0/r0)
e = sqrt( max(temp,0.d0) )
q = p / (1.d0 + e)
c
c If the object hit the central body
if (q.le.rcen) then
nhit = nhit + 1
jhit(nhit) = j
dhit(nhit) = rcen
c
c Time of impact relative to the end of the timestep
if (e.lt.1d0) then
a = q / (1.d0 - e)
uhit = sign (acos((1.d0 - rcen/a)/e), -h)
u0 = sign (acos((1.d0 - r0/a )/e), rv0)
mhit = mod (uhit - e*sin(uhit) + PI, TWOPI) - PI
m0 = mod (u0 - e*sin(u0) + PI, TWOPI) - PI
else
a = q / (e - 1.d0)
uhit = sign (mco_acsh((1.d0 - rcen/a)/e), -h)
u0 = sign (mco_acsh((1.d0 - r0/a )/e), rv0)
mhit = mod (uhit - e*sinh(uhit) + PI, TWOPI) - PI
m0 = mod (u0 - e*sinh(u0) + PI, TWOPI) - PI
end if
mm = sqrt((mcen + m(j)) / (a*a*a))
thit(nhit) = (mhit - m0) / mm + time
end if
end if
end do
c
c------------------------------------------------------------------------------
c
return
end
c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
c MCE_COLL.FOR (ErikSoft 2 October 2000)
c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
c Author: John E. Chambers
c
c Resolves a collision between two objects, using the collision model chosen
c by the user. Also writes a message to the information file, and updates the
c value of ELOST, the change in energy due to collisions and ejections.
c
c N.B. All coordinates and velocities must be with respect to central body.
c ===
c
c------------------------------------------------------------------------------
c