-
Notifications
You must be signed in to change notification settings - Fork 3
/
moon_orbits.f
1125 lines (953 loc) · 35.1 KB
/
moon_orbits.f
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
program moon_orbits
c Ramdom moons and orbits
implicit double precision (a-z)
c================================================
c copyright (c) Selden Ball 2012
character*5 version/'v1.19'/
c ALL RIGHTS RESERVED
c license:
c this program may be freely distributed for educational purposes
c but any commercial use requires written permission from the author
c changes:
c 1.01 cleaned up code so ftnchek doesn't complain so much
c 1.02 replaced calls to ran(dmy) by calls to randf()
c 1.03 explicit output formatting instead of just *
c 1.04 place orbits between existing orbits
c instead of randomly within max radius
c 1.05 prefix moon by name of planet so it'll be unique in this star system
c 1.06 include option name when complaining about it being missing
c 1.07 include TestPlanet specs if appropriate
c 1.08 include program name in catalog file
c 1.09 randomize sma of moon "better"
c 1.10 randomize initial mass available
c 1.11 add atmosphere for large moons
c 1.12 note "Earthlike?" for moons >0.75 Earth masses
c 1.13 more code cleanup
c 1.14 use color
c 1.15 make version variable
c 1.16 limit max planet radius: hill radius calcs fail
c 1.17 display planetary type
c 1.18 add rocky planet radius guesstimate
c 1.19 specify appropriate texture for test planet
c================================================
C this program was developed using g77 v3.4.4
c it seems to work correctly with gfortran v4.5.3 after cleanup
c compilation instructions:
c g77 moon_orbits.f -o moon_orbits
c usage
c ./moon_orbits -?
c Note: all options require a space between the option and its value
c================================================
c someday need to check for
c positive resonances:
c 1:2:4 in orbital periods
c 3:2, 4:3, 3:4, 3:5, 4:7, 1:2, 2:3, 2:4, 2:5
c BUT need adequate separation in orbital phase
c negative resonances:
c 1:1, (except for trojans? -- not likely with many perturbators)
c 3:1, 5:2, 7:3, 2:1, 7:6
integer Nmoons,Nmoon,i,max_moons,Ntest
integer masses(0:11),npc,mmp1,mmp2
DATA MASSES/12*0/
c hill radius = a(1-e)* (m/3M)^(1/3)
c planet & star parameters
parameter (max_moons = 200 )
parameter (mmp1=max_moons+1)
parameter (mmp2=max_moons+2) ! includes minimum and maximum radii
c sma and hill radii of moons
dimension Am (0:mmp1),Hm (0:mmp1)
data Am/mmp2*0.0/,Hm/mmp2*0.0/
character*20 Nameplanet,Namestar,FileName
c---------------------
call get_args
* (FileName,Nameplanet,Mplanet,Rplanet,Aplanet, NameStar,Mstar,
* Nmoons, Mmoons, version)
print *, Nmoons
call str_len(FileName,npc)
c if (Nmoons.gt.200) then
c print *, 'Nmoons (',Nmoons,') is more than 200. Set to 200.'
c Nmoons=200
c endif
open(1,file='Moons_'//FileName(:npc)//'.txt',status='replace')
write (1,*) '# moon catalog for Celestia generated by'
write (1,*) '# moon_orbits '//version
c Scaling parameters
Mjup = 317.82 d0 ! mass of Jupiter in earth masses
Rjup = 71 492.0 d0 ! radius of Jupiter in km (equatorial)
Msol = 332 918.215 d0 ! mass of Sol in earth masses
c Msolkg = 1.98892 d30 ! mass of Sol in kg
kmperau = 149 598 000.0 d0 ! km per au
Mearth = 5.9742 d24 ! mass of Earth in kg
c check mass of sun & Jupiter
c print *, 'Msol, Msolkg, Mearth, Msol*Mearth, mjup/msol'
c print *, Msol,Msolkg, Mearth, Msol*Mearth,
c * (Mjup*Mearth)/(Msol*Mearth)
c=====================================================
c default input parameters
c Mplanet = 300d0 ! planet mass in Earth masses
c Aplanet = 1.0d0 ! au
c Mstar = 1.0d0*Msol ! Star mass
c=====================================================
c hill radius = a(1-e)* (m/3M)^(1/3)
c looks ok for relationship of things within the radius
c like orbits of moons,
c but 1+e is more appropriate for orbital interference, I think.
C Hill radius of planet (max moon distance in au)
Hplanet=
* Aplanet*(Mplanet/(3.0d0*Mstar*Msol))**(1.0d0/3.0d0)
c print *,'Aplanet,Mplanet,Mstar,Msol',Aplanet,Mplanet,Mstar,Msol
c print *,'Planet Hill R (AU) = ',Hplanet
Hplanet = kmperau * Hplanet
c print *,'Planet Hill R (km) = ',Hplanet
c print *,'Hplanet= ',Hplanet
c not stable unless within 1/3 of hill radius per wikipedia
c (also more aesthetic)
Hplanet = Hplanet/3.0d0
Am(1) = Hplanet ! initial maximum radius for placing moons
c call draw_hill(Nameplanet,Hplanet)
c moon parameters
c total mass of all moons = 1.0-3.0 e-4 of planet
c Ref: Canup and Ward (2006)
c A common mass scaling for satellite systems of gaseous planets.
c Nature, 441: 834-839.
c see also:
c http://phl.upr.edu/library/notes/themassandradiusofpotentialexomoons
c
c rm = (1.0d-4) + (1.0d-4)*rgauss(1.0d0)
c rm = 1.5d-4
if (Mplanet.gt.1000.0d0) then
rm = 0.05d0
elseif (Mplanet.gt.10.0d0) then
c rm = (1.1d0 +1.4d0*randf())/10 000.0d0
rm = 1.5d0/10000.0d0
elseif (Mplanet.gt.0.5d0) then
c rm = (1.1d0 +1.4d0*randf())/10 0.0d0
c rm = 1.0d0/100.0d0
if (Rplanet.gt.15000) then
rm = 0.01d0
else
rm = 0.05d0
endif
else
c rm = (1.1d0 +1.4d0*randf())/10 000.0d0
rm = 1.5d0/10000.0d0
endif
c rm = 0.05d0
Mtot = rm *Mplanet ! in units of earth mass
Mtot = Mtot - Mmoons
c Mmin = Mtot/1 000 000.0d0 ! no smaller than 1/1 000 000 of total mass
Mmin = 0.0d0
c Mmin = Mtot/1 000 000 000.0d0;
Mmoonmass = 1.0d0 / 100 000 000 000.0d0
earth_size_moons = 0;
luna_size_moons = 0;
c ------------------------------------------------------------------
c planet radius in km
c radius peaks at 1.65 x Mjup = 1.051 * Rjup
c graph has min at 6Mjup
c Ref: Giant Planet Interior Structure and Thermal Evolution
c Fortney et al. (2009), Fig. 4
c print *, ' '
c if (Mplanet .le. 0.25d0*Mjup) then
c it's rocky
c print *, 'This planet is rocky.'
c Rplanet = ((Mplanet/Mjup)**0.333d0)*Rjup
c elseif (Mplanet .le. 1.65d0*Mjup) then
c assume it's a gas giant
c print *, 'This planet is a gas giant.'
c Rplanet = ((Mplanet/Mjup)**0.10d0)*Rjup
c elseif (Mplanet.le. 6.04d0*Mjup) then
c it's a brown dwarf
c print *, 'This "planet" is a brown dwarf.'
c Rplanet = ((Mplanet/Mjup)**(-0.125d0))*Rjup
c else
c It's a star
c print *, 'This "planet" is a dim star.'
c Rplanet = 0.799d0*Rjup + (Mplanet-(6.04d0*Mjup))*Rjup
c if (Rplanet.gt.10.0d0*Rjup) then
c print *, 'Radius = ',Rplanet,' km = ',Rplanet/Rjup,' x Jup.'
c print *, 'This "planet" is too large and massive: '
c print *, ' Hill radius calculations fail.'
c print *, 'Please use StarGen instead.'
c stop
c endif
c endif
c print *, 'Radius = ',Rplanet,' km = ',Rplanet/Rjup,' x Jup.'
c Density of planet (for Roche Limit calculation)
c mass = density * volume
c density = mass/volume
C Volume = 4/3 pi r^3
c density = mass/ (4.0/3.0* pi *r**3)
c
c but need to get units right......
c density of water = 1 g/cm^3 = 1 kg / m^3
c 1 km^3 = 1E6 m^3
Vplanet = (4.0d0/3.0d0)*3.14159265d0*Rplanet**3 !km^3
Dplanet = (Mearth* Mplanet) / Vplanet ! kg/km^3
c scale:
c 1000 grams/kg = 1d3
c 1000x1000x1000 meters^3/km^3 = 1d9
c 100x100x100 cc / m^3 = 1d6
Dplanet = Dplanet *(1d3/(1d6*1d9)) ! g/cc
c jupiter's radius = ~70 000 km
c metis (innermost moon): 128 000 km (1.8x planet's radius)
c so make innermost orbit 1.5x planet's radius
Amin= 1.5d0* Rplanet
Am(0)=Amin ! initial minimum radius for placing moons
c print *, 'Rplanet,Am(0),Am(1) = ',Rplanet,Am(0),Am(1)
c include dummy planet specs if appropriate
if (Nameplanet.eq.'TestPlanet') then
call write_planet(Nameplanet,Namestar,Rplanet,Aplanet,Mplanet)
endif
c random orbit radii out to hill radius of planet vs its sun
c (rocky moon orbiting gas giant: roche limit is inside planet)
c but put no moons within some multiple of hill radius of each other
c specify masses in terms of earth mass
c convert to kg later
Mremaining = Mtot
Nmoon = 0
Ntest = 0
100 continue
Ntest = Ntest+1
! print *,' Ntest=',Ntest
c if (Ntest.gt.10) stop
c quit if we've run out of mass
c or if the specified number of moons have been created
! print *, Mremaining, ' remaining: ',100*Mremaining/Mtot,'%'
if (Mremaining.lt.Mmin) then
print *, 'Available mass exhausted'
print *, 'after generating ',Nmoon,' moons.'
endif
if (Mremaining.lt.Mmin .or.Nmoon.ge. Nmoons) goto 200
c select mass of this moon:
c a random fraction of remaining mass
r1 = randf()
c r2 = 1.0d0- LOG10(1.0d0+9.0d0*r1)
maxln = exp(1.0d0)
r2 = LOG(1.0d0+maxln*r1)/1.33d0
c print *,'r2= ',r2
Mmoon = r2 *Mremaining
if (Mmoon.lt.Mmoonmass) then
Mmoon = Mmoonmass
endif
c if (Mmoon.gt.2.5d0) then
c if (earth_size_moons.lt.2) then
c Mmoon = 1.4d0 * randf()
c earth_size_moons = earth_size_moons + 1
c elseif (luna_size_moons.lt.6) then
c Mmoon = 0.014d0 * randf()
c luna_size_moons = luna_size_moons + 1
c else
c Mmoon = 0.00014d0 * randf()
c endif
c endif
Mremaining = Mremaining - Mmoon
c print *, 'Mmoon = ', Mmoon,' (mass)'
c print *, 'placing moon'
c insert this moon among the others
call place_moon(Am,Hm,mmp1,
* Nmoon,Amoon,Mmoon,Mtot,Mplanet)
c print *, 'placed moon'
c if (Ntest.gt.5) stop 'Ntest>5'
c add to mass distribution
i = int(Mmoon* (10d0/Mtot))
IF (i.lt.1) i=0
if (i.ge.11) i=11
masses(i)=masses(i)+1
c and write to ssc file (with some other variables)
call write_moon
* (Nmoon,Amoon,Mmoon,Mtot,
* Mplanet,Rplanet,Dplanet,Hplanet,Aplanet,
* Nameplanet, NameStar)
c print *, 'wrote moon'
goto 100
200 continue
print *, ' '
print '(a,f6.3)', ' Total moon mass: ',Mtot
PRINT '(a,12i4)', ' Mass distribution: ',(MASSES(I),I=0,11)
call str_len(FileName,npc)
print *,' Done: ',Nmoon,' moons in ',
* 'Moons_'//FileName(:npc)//'.txt'
print *,'==================================='
print *, ' '
close (1)
end
c-----------------------------------------------------
subroutine draw_hill(Nameplanet,Hplanet)
character*20 nameplanet
double precision Hplanet
call str_len(nameplanet,npl)
write (1,*) '"HillRadius" "Sol/',Nameplanet(:npl),'"'
write (1,*) '{'
write (1,*) ' Class "moon"'
write (1,*) ' Radius 1000'
write (1,*) ' Color [ 1 0 0 ]'
write (1,*) ' EllipticalOrbit '
write (1,*) ' {'
write (1,*) ' Period 1000'
write (1,*) ' SemiMajorAxis ',Hplanet
write (1,*) ' }'
write (1,*) '}'
write (1,*) ' '
RETURN
END
C===========================
c write this moon to the ssc file
c after determining its other parameters
subroutine write_moon
* (Nmoon,Amoon,Mmoon,Mtot,
* Mplanet,Rplanet,Dplanet,Hplanet,Aplanet,
* Nameplanet, NameStar)
implicit double precision (a-z)
integer Nmoon,npc,nsc,nch
character*20 Nameplanet,Namestar,Tmoon
character*13 color
logical lsphere
c derive period
c calculate appropriate gravitational constant
c from the Moon's orbit around the Earth
AM = 384748d0 ! semi-major axis in km
PM = 27.32166155d0 ! sidereal period in days
MM = 0.0123000383d0 !fraction of Earth's mass
K = (PM**2) / ((AM**3)/(1.0d0+MM))
c print *,' K= ',K
c print *, 'Nmoon,Amoon,Mmoon,Mplanet,Aplanet,Dplanet = '
c print *, Nmoon,Amoon,Mmoon,Mplanet,Aplanet,Dplanet
c derive period, ignoring moon's mass
C p^2 = K a^3/m
c moon's period: K (above) SMA (km, provided) Mplanet (earths, provided)
Pmoon = sqrt(K* Amoon**3/Mplanet)
c random composition (density)
c Ramdom density 1-5: ice to metal
c except that the farther out a planet is,
c the lower the density should be (more ices, less metals)
c my guess:
c most likely density = 1.0 for outer planet
c most likely density = 3.3 for inner planet
c fractional range of planetary distances:
c limit planet SMA distance to 40 au
Ap = min(Aplanet,40.0d0)
c and discover its fractional distance
den = (40.0d0-Ap)/15.0d0
c print *, ' den = ',den
c (approximate poisson distribution)
x1 = rgauss (0.54d0)
c print *, x1,x2
IF (x1.lt.0.0d0) X1=-X1
Dmoon = 0.90d0+den*X1
c print *, 'Dmoon = ',Dmoon
c heck with it: limit max density
if (Dmoon.gt.3.3d0) Dmoon = 3.3d0
c Dmoon = 1.0d0+randf()*4.0d0 ! g/cm^3
C calculate roche limit for moon:
c generate rings instead of moon
c if an icy moon is inside its roche limit
c
d = 2.423d0* Rplanet*(2d0*Dplanet/Dmoon)**(1d0/3d0)
c print *,'Dmoon, d, Amoon, diff = ',Dmoon,d,Amoon, Amoon-d
c but not yet
if (Amoon .lt. d .and. Dmoon .lt. 1.2d0) then
print *, 'Rings! (but not yet)'
endif
c Mass as multiple of Earth's mass
c Density in g/cc
c returns Radius in km
call Radius_moon(Mmoon,Dmoon,Rmoon)
c color depends on composition, too
if (Dmoon.lt.2.0d0) then
color = '[0.8 0.8 0.8]'
elseif (Dmoon.lt.3.5d0) then
color = '[0.5 0.4 0.3]'
else
color = '[0.3 0.4 0.5]'
endif
c as does the albedo
Albedo = 0.2d0/Dmoon
Mearth = 5.9742d24 ! earth mass in kg
Kgmoon = Mmoon*Mearth
c and whether it's spherical due to hydrostatic equilibrium
c BUT function isn't entirely obvious
c 250km is close enough for moons of Saturn *shrug*
if (Rmoon.lt.250d0) then
lsphere = .false.
else
lsphere=.true.
endif
c randomized orbital parameters
c probably should have some relationships with each other
c (resonances)
c but not now.
c Inclination
c probably larger with larger distances
c probably inversely proportional to cube of distance
c due to tidal effects
c BUT that expands much too fast
fdist = Amoon/Hplanet
c and larger with smaller masses?
c also expands much too fast
c fmass = Mtot/Mmoon
fmass = (Mtot/Mmoon)**0.2d0
c and larger with smaller oblateness of planet
c but dunno what the oblateness is
ri=rgauss(1.0d0)
Imoon = fdist*fmass*ri
c Eccentricity
c major moons <0.01. more distant minor can be larger.
Emoon = 0.01d0*randf()
c ArgOfPericenter
AOPmoon = 360.0d0*randf()
c AscendingNode
ANmoon = 360.0d0*randf()
c MeanAnomaly
MAmoon = 360.0d0*randf()
c Rotational parameters are omitted
c thus causing all moons to be tidally locked
c atmospheric parameters are omitted for smaller moons
c included if mass & density .ge. titan. see below
call str_len(Nameplanet, npc)
call str_len(Namestar, nsc)
c replace underscores by spaces in star name
call u2sp(Namestar)
call left_justify(Nmoon,Tmoon,nch)
write (1,*) '"'//Nameplanet(:npc)//'_Moon ',Tmoon(:nch),'" "',
* Namestar(:nsc)//'/'//NamePlanet(:npc)//'"'
write (1,*) '{'
write (1,*) ' Class "moon"'
write (1,fmt='(a,e50.32,a,1p,e50.32,a)')
* ' # mass = ',Mmoon,' x Earth = ',Kgmoon,' kg'
write (1,fmt='(a,f6.3)') ' # density = ',Dmoon
if (Mmoon.gt. 0.75d0) then !1.11
print *, ' Earthlike moon?'
write(1,*) '# Earthlike?'
endif
write (1,fmt='(a,f50.32)')' Radius ',Rmoon
IF (.NOT. lsphere) write(1,*) ' Mesh "asteroid.cms"'
if (Dmoon.lt.1.0d0) then
write (1,*) ' Texture "tethys.*"'
elseif (Dmoon.lt.1.5d0) then
write (1,*) ' Texture "enceladus.*"'
elseif (Dmoon.lt.2.5d0) then
write (1,*) ' Texture "callisto.*"'
elseif (Dmoon.lt.3.5d0) then
write (1,*) ' Texture "ganymede.*"'
else
write (1,*) ' Texture "europa.*"'
endif
write (1,*) ' Color '//color ! 1.14
write (1,*) ' BlendTexture true'
write (1,*) ' EllipticalOrbit '
write (1,*) ' {'
write (1,fmt='(a,f50.32)') ' Period ',Pmoon
write (1,fmt='(a,f50.32)') ' SemiMajorAxis ',Amoon
write (1,fmt='(a,f48.32)') ' Eccentricity ',Emoon
write (1,fmt='(a,f48.32)')
* ' Inclination ',Imoon!,' # ',fdist,fmass,ri
write (1,fmt='(a,f48.32)') ' ArgOfPericenter ',AOPmoon
write (1,fmt='(a,f48.32)') ' AscendingNode ',ANmoon
write (1,fmt='(a,f48.32)') ' MeanAnomaly ',MAmoon
write (1,*) ' }'
c all moons are assumed to be tidally locked to the planet
c hence no rotation specificiation
c include atmosphere if mass & density >= titan
if (Kgmoon .gt. 1.3d23 .and. Dmoon .gt. 1.8d0) then
c duplicate Titan's atmosphere for now
write (1,*) ' Atmosphere {'
write (1,*) ' Height 500.0'
write (1,*) ' Lower [ 0.477 0.367 0.211 ]'
write (1,*) ' Upper [ 0.96 0.805 0.461 ]'
write (1,*) ' Sky [ 0.3 0 0 ]'
write (1,*) ' Mie 0.0001'
write (1,*) ' MieAsymmetry -0.55'
write (1,*) ' Rayleigh [ 0.0 0.0 0.00017 ]'
write (1,*) ' Absorption [ 0.000075 0.00030 0.00025 ]'
write (1,*) ' MieScaleHeight 220.0 '
write (1,*) ' }'
endif
write (1,fmt='(a,f6.3)') ' Albedo ',Albedo
write (1,*) '}'
write (1,*) ' '
return
end
c==============================================================
subroutine Radius_moon ( Mmoon,Dmoon, Rmoon)
c calculate radius (depends on composition)
c volume = mass/density = 4/3 pi R^3
c r^3 = (3/4)*(mass/density)/pi
c radius = (0.75d0/3.14159264d0*(mass/density))**(1/3)
c needs to be in km
c
c verified with moons of Jupiter & Saturn
c max deviation = 0.3% (0.003x)
implicit double precision (a-z)
k = 0.75d0/3.14159264d0
Mearth = 5.9742 d27 ! mass of Earth in grams
Mmoon_in_g = Mmoon*Mearth
volume = Mmoon_in_g /Dmoon
c print *, 'Mmmon = ',Mmoon_in_g,' grams, Dmoon = ',Dmoon
c print *, 'volume= ',volume
Rmoon = (k*volume)**0.3333333333333333333d0
c print *, 'Rmoon = ',Rmoon
Rmoon = Rmoon/(100.0d0 * 1000.0d0) ! cm -> m, m -> km
c print *, 'Rmoon = ',Rmoon
return
end
c==============================================================
subroutine get_args
* (FileName,Nameplanet,Mplanet,Rplanet,Aplanet, NameStar,Mstar,
* Nmoons,Mmoons,version)
implicit double precision (a-z)
integer nf,Nmoons,iseed,na,lt,ln
character*20 Nameplanet,Massplanet,SMAplanet,MassStar,NumMoons
character*20 NameStar,RanSeed,tmp,FileName
character*20 Radiusplanet
character*20 MoonMass
character*5 version
character*2 flg
character*1 char
logical Lname,Lmp,Lsmap,Lmstar,Lseed,Lnum,Lnstar,Lfn,Lrad,Lmmoon
data Lname,Lmp,Lsmap,Lmstar,Lseed,Lnum,Lnstar,Lfn,Lrad /9*.false./
c -a planet sma
c -n planet name
c -N number of moons
c -U mass of moons generated by Stargen
c -m planet mass
c -R planet radius
c -M star mass
c -s seed for ran
c -S star name
c default input parameters
Msol = 332 918.215d0 ! earth masses
Mplanet = 300d0 ! planet mass in Earth masses
Aplanet = 1.0d0 ! au
Mstar = 1.0d0*Msol ! Star mass
Rplanet = 1.0d0
c used to calculate default radius
Mjup = 317.82 d0
Rjup = 71 492.0 d0
c=====================================================
print *, ' moon_orbits '//version
do nf = 1,17,2
call getarg(nf,flg)
if (flg.eq.'-a') then
call getarg(nf+1, SMAplanet)
Lsmap=.true.
elseif (flg.eq.'-F') then
call getarg(nf+1, FileName)
Lfn = .true.
elseif (flg.eq.'-m') then
call getarg(nf+1, MassPlanet)
Lmp = .true.
elseif (flg.eq.'-R') then
call getarg(nf+1, RadiusPlanet)
Lrad=.true.
elseif (flg.eq.'-M') then
call getarg(nf+1, MassStar)
Lmstar=.true.
elseif (flg.eq.'-U') then
call getarg(nf+1, MoonMass)
Lmmoon=.true.
elseif (flg.eq.'-n') then
call getarg(nf+1, Nameplanet)
Lname = .true.
elseif (flg.eq.'-N') then
call getarg(nf+1, NumMoons)
Lnum = .true.
elseif (flg.eq.'-s') then
call getarg(nf+1, RanSeed)
Lseed=.true.
elseif (flg.eq.'-S') then
call getarg(nf+1, NameStar)
c assume it has several parts (and is last on line)
do na =2,4
call getarg(nf+na, tmp)
if (tmp.ne.' ' .and. tmp(1:1).ne.'-'
* .and.tmp(1:1).ne.char(13)) then
call str_len(tmp,lt)
if (lt.gt.0) then
call str_len(NameStar,ln)
NameStar(ln+1:)='_'//tmp
endif
endif
enddo
Lnstar=.true.
goto 13
else
if (flg.eq.' ') goto 13
print *, 'Unrecognized option: ',flg
print *, ' '
print *, ' moon_orbits '//version
print *, ' generates SSC for major moons of a gas giant'
print *, ' '
print *, ' put a space between the option flag and its value'
print *, ' -a planet_sma (in AU; default = 5.0)'
print *, ' -F file name (not extension)'
print *, ' -m planet_mass (as fraction of Earth; default=300.0)'
print *, ' -R planet radius (in km)'
print *, ' -M star_mass (as fraction of Sun; default = 1.0)'
print *, ' -n planet_name (max 20 chars; default = TestPlanet )'
print *, ' -N number_of_moons (default = 10)'
print *, ' -s seed (for random number generator; def=14159262)'
print *, ' -S star_name (max 20 chars; default = Sol)'
stop
endif
enddo
13 continue
if (Lname) then
print *, 'Planet name = ',Nameplanet
else
Nameplanet = 'TestPlanet'
print *, 'Planet name (-n) not specified, using ',Nameplanet
print *, ' (Including dummy planet definition in ssc.)'
endif
if (Lmp) then
print *, 'Planet mass = ',MassPlanet
call fix_float(MassPlanet)
else
MassPlanet = '300.0'
print *, 'Planet mass (-m) not specified, using ',MassPlanet
endif
read(MassPlanet,fmt='(d20.10)') Mplanet
if (Lrad) then
print *, 'Planet radius = ',RadiusPlanet
call fix_float(RadiusPlanet)
read(RadiusPlanet,fmt='(d20.10)') Rplanet
else
print *, ' '
if (Mplanet .le. 0.25d0*Mjup) then
print *, 'This planet is rocky.'
Rplanet = ((Mplanet/Mjup)**0.333d0)*Rjup
elseif (Mplanet .le. 1.65d0*Mjup) then
print *, 'This planet is a gas giant.'
Rplanet = ((Mplanet/Mjup)**0.10d0)*Rjup
elseif (Mplanet.le. 6.04d0*Mjup) then
print *, 'This "planet" is a brown dwarf.'
Rplanet = ((Mplanet/Mjup)**(-0.125d0))*Rjup
else
print *, 'This "planet" is a dim star.'
Rplanet = 0.799d0*Rjup + (Mplanet-(6.04d0*Mjup))*Rjup
if (Rplanet.gt.10.0d0*Rjup) then
print *, 'Radius = ',Rplanet,' km = ',Rplanet/Rjup,' x Jup.'
print *, 'This "planet" is too large and massive: '
print *, ' Hill radius calculations fail.'
print *, 'Please use StarGen instead.'
stop
endif
endif
print *, 'Planet radius = ',Rplanet
endif
if (Lsmap) then
print *, 'Planet SMA = ',SMAplanet
call fix_float(SMAPlanet)
else
SMAPlanet = '5.0'
print *, 'Planet SMA (-a) not specified, using ',SMAPlanet
endif
if (Lnstar) then
print *, 'Star name = ',NameStar
else
NameStar = 'Sol'
print *, 'Star name (-S) not specified, using ',NameStar
endif
if (Lmstar) then
print *, 'Star mass = ',MassStar
call fix_float(MassStar)
else
MassStar = '1.0'
print *, 'Star mass (-M) not specified, using ',MassStar
endif
if (Lseed) then
print *, 'Random seed = ',RanSeed
call fix_float(RanSeed)
else
RanSeed = '14159262.0'
print *, 'Seed (-s) not specified, using ',RanSeed
endif
if (Lnum) then
print *, 'Number of Moons = ',NumMoons
call fix_float(NumMoons)
else
NumMoons = '10.0'
print *, 'Number of moons (-N) not specified, using ',NumMoons
endif
if (Lfn) then
print *, 'File name = ',FileName
else
call str_len(NameStar,ln)
FileName = NameStar(:ln)//'_'//Nameplanet
print *, 'File name (-F) not specified, using ',FileName
endif
c translate ascii values into appropriate numbers
read(NumMoons,fmt='(d20.10)') FNmoons
Nmoons = int(fnmoons)
read(MassPlanet,fmt='(d20.10)') Mplanet
read(MassStar,fmt='(d20.10)') Mstar
read(SMAPlanet,fmt='(d20.10)') Aplanet
read(RanSeed,fmt='(d20.10)') Dseed
read(MoonMass,fmt='(d20.10)') Mmoons
iseed = int(Dseed)
call srandf(iseed)
return
END
c-----------------------------------------------------
subroutine fix_float(a0)
c ensure a decimal point in number
character*20 a0
npd = index(a0,'.')
if (npd.ne.0) return
npe = index(a0,'e')
if (npe.eq.0) npe = index(a0,'E')
if (npe.eq.0) npe = index(a0,'d')
if (npe.eq.0) npe = index(a0,'D')
if (npe.eq.0) then
call str_len(a0,np)
a0(np+1:np+2)='.0'
return
else
a0=a0(:npe-1)//'.0'//a0(npe:)
endif
return
end
c-----------------------------------------------------
subroutine str_len(str,ncf)
character*(*) str
lens = len(str)
do nc = lens,1,-1
if (str(nc:nc).ne. ' ') then
ncf = nc
return
endif
enddo
ncf = 1
return
end
c=============================================================
SUBROUTINE SRANDf(ISEED)
C
C This subroutine sets the integer seed to be used with the
C companion RAND function to the value of ISEED. A flag is
C set to indicate that the sequence of pseudo-random numbers
C for the specified seed should start from the beginning.
C
COMMON /SEED/JSEED,IFRST
C
JSEED = ISEED
IFRST = 0
C
RETURN
END
c------------------------
double precision FUNCTION RANDf()
C
C This function returns a pseudo-random number for each invocation.
C It is a FORTRAN 77 adaptation of the "Integer Version 2" minimal
C standard number generator whose Pascal code appears in the article:
C
C Park, Steven K. and Miller, Keith W., "Random Number Generators:
C Good Ones are Hard to Find", Communications of the ACM,
C October, 1988.
C
PARAMETER (MPLIER=16807,MODLUS=2147483647,MOBYMP=127773,
+ MOMDMP=2836)
C
COMMON /SEED/JSEED,IFRST
INTEGER HVLUE, LVLUE, TESTV, NEXTN
SAVE NEXTN
C
IF (IFRST .EQ. 0) THEN
NEXTN = JSEED
IFRST = 1
ENDIF
C
HVLUE = NEXTN / MOBYMP
LVLUE = MOD(NEXTN, MOBYMP)
TESTV = MPLIER*LVLUE - MOMDMP*HVLUE
IF (TESTV .GT. 0) THEN
NEXTN = TESTV
ELSE
NEXTN = TESTV + MODLUS
ENDIF
RANDf = dble(NEXTN)/dble(MODLUS)
C
RETURN
END
BLOCKDATA RANDBD
COMMON /SEED/JSEED,IFRST
C
DATA JSEED,IFRST/123456789,0/
C
END
!****************************************************************************
! RGAUSS
! PURPOSE: generate gaussian random numbers with given sigma
! Uses FUNCTION RANDF() to generate uniform random numbers between 0 and 1
! Uses the Box-Muller algorithm to generate numbers with
! a Gaussian distribution of probability.
!
double precision function rgauss(sigma)
real*8 x1, x2, w, y1, sigma,randf
100 continue
x1 = 2.0d0 * randf() - 1.0d0
x2 = 2.0d0 * randf() - 1.0d0
w = x1 * x1 + x2 * x2
c print *,'x1,x2, w = ',x1,x2,w
if ( (w .ge. 1.0d0).or.(w.eq.0.0d0) ) goto 100
w = sigma*sqrt( (-2.0d0 * log( w ) ) / w )
c print *,' w2 = ',w
y1 = x1 * w
c y2 = x2 * w
rgauss = y1
return
end
c===========================================
subroutine u2sp(Namestar)
c restore spaces in star name
character*20 Namestar
do nc = 1,20
if (Namestar(nc:nc).eq.'_') Namestar(nc:nc)=' '
enddo
return
end
c===========================================
subroutine left_justify(number,string,nch)
c left justifies number in string
integer number,nch,ns,nsp
character*20 string
write(unit=string,fmt=*) number
c eliminate leading spaces
do ns = 1,20
if (string(ns:ns).ne.' ') then
nsp = ns
goto 10
endif
enddo
10 continue
if (nsp.ge.1) string = string(nsp:)
c eliminate trailing spaces
do ns = 1,20
if (string(ns:ns).eq.' ') then
nch = ns-1
goto 20
endif
enddo
20 continue
c print *, string, nsp,nch
return
end
c=================================================
subroutine place_moon(Am,Hm,mmp1,
* Nmoon,Amoon,Mmoon,Mtot,Mplanet)
implicit double precision (a-z)
integer mmp1,Nmoon,nm,MoonThis
dimension Am(0:mmp1),Hm(0:mmp1)
c print *, 'place_moon: Nmoon, am0,am1,hm0,hm1=',
c * Nmoon,am(0),am(1),hm(0),hm(1)
c stop
c Note:
c hill radius = a(1-e)* (m/3M)^(1/3)
c greater eccentricity shrinks the hill radius
c which is OK for stability of orbits within the hill radius
c but the opposite seems more appropriate for
c orbital interference; max e is 1 so use that