home *** CD-ROM | disk | FTP | other *** search
/ Frozen Fish 1: Amiga / FrozenFish-Apr94.iso / bbs / alib / d9xx / d969 / ace.lha / ACE / ACE-2.0.lha / PRGS.lha / Astronomy / Jovian.b < prev    next >
Text File  |  1994-01-10  |  12KB  |  522 lines
  1. { Displays the positions of the Galilean
  2.   satellites relative to Jupiter for a 
  3.   given date and period of time.
  4.  
  5.   The view is as it would appear through
  6.   an inverting telescope in the Southern
  7.   Hemisphere.
  8.  
  9.   All dates and times must be entered in 
  10.   Universal (Greenwich Mean) Time.
  11.  
  12.   Method taken from Jean Meeus' 
  13.   "Astronomical Algorithms", ch 43.
  14.  
  15.   Author: David J Benn
  16.     Date: 11th,12th,13th,17th,26th July 1993,
  17.       18th December 1993 }
  18.  
  19. CONST radconv=57.295779
  20. CONST xorigin=320!, yorigin=116!
  21. CONST radius=5
  22. CONST true=-1&, false=0&
  23. CONST JovianColor=2
  24. CONST black=0
  25. CONST xscale=10, yscale=7.5
  26.  
  27. SINGLE   JDE
  28. SINGLE   d,V,M,N,J,A,B,K,R,rr,delta,psi,DE
  29. SINGLE   first_u1,first_u2,first_u3,first_u4
  30. SINGLE   u1,u2,u3,u4
  31. SINGLE   G,H
  32. SINGLE   r1,r2,r3,r4
  33. LONGINT  first
  34. SHORTINT moon
  35.  
  36. DIM x(4),y(4),lastx(4)
  37.  
  38. SUB SINGLE decimal_hours(hrs,mins,secs)
  39. '..return decimal hours
  40.  
  41.   decimal_hours = hrs + mins/60 + secs/3600
  42. END SUB
  43.  
  44. SUB SINGLE JulianDay(dy,mn,yr) 
  45. SINGLE m1,y1,a,b,c,d,dj
  46. '..This routine calculates the number of days elapsed 
  47. '..since the epoch 1900 January 0.5 (ie: 1200 GMT, 31st Dec 1899). 
  48.  
  49.  if yr = 0 then 
  50.    dj=-1  '..error 
  51.  else 
  52.    m1=mn : y1=yr : b=0 
  53.  
  54.   if y1 < 1 then ++y1 
  55.   if mn < 3 then m1=mn+12 : --y1 
  56.  
  57.   if y1 > 1582 or mn > 10 or dy >= 15 then  
  58.     a=int(y1/100) : b=2-a+int(a/4) 
  59.     c=int(365.25*y1)-694025 
  60.     if y1 < 0 then c=fix((365.25*y1)-0.75)-694025 
  61.     d=int(30.6001*(m1+1))
  62.     dj=b+c+d+dy-0.5 
  63.   else 
  64.     if (y1<1582 or (y1=1582 and mn<10) or (y1=1582 and mn=10 and dy<5)) then
  65.       c=int(365.25*y1)-694025 
  66.       if y1 < 0 then c=fix((365.25*y1)-0.75)-694025 
  67.       d=int(30.6001*(m1+1)); dj=b+c+d+dy-0.5 
  68.     else 
  69.       dj=-1  '..error
  70.     end if
  71.   end if 
  72.  end if
  73.  
  74.  JulianDay = dj  '..Return Julian Date (error = -1)
  75. END SUB
  76.  
  77. SUB STRING time_from_day_fraction(SINGLE fd)
  78. '..return time string from day fraction.
  79. SINGLE hrs,mins
  80.   
  81.   hrs = 24*fd
  82.   mins = 60*(hrs-fix(hrs))
  83.   hr$ = str$(fix(hrs))
  84.   min$ = str$(fix(mins))
  85.   hr$ = right$(hr$,len(hr$)-1)
  86.   min$ = right$(min$,len(min$)-1)
  87.   if len(hr$)=1 then hr$="0"+hr$
  88.   if len(min$)=1 then min$="0"+min$
  89.  
  90.   time_from_day_fraction = hr$+":"+min$
  91. END SUB
  92.  
  93. SUB STRING date_and_time(dj!)
  94. '..This routine converts the number of (Julian) days since 
  95. '..1900 January 0.5 into the calendar date and time.
  96. SINGLE a,b,c,d,g,i,fd
  97.  
  98.  d=dj!+0.5 : i=int(d) : fd=d-i 
  99.  
  100.  if fd = 1 then fd=0 : ++i 
  101.  
  102.  if i > -115860 then 
  103.    a=int((i/36524.25)+9.9835726e-1)+14 
  104.    i=i+1+a-int(a/4) 
  105.  end if 
  106.  
  107.  b=int((i/365.25)+8.02601e-1) 
  108.  c=i-int((365.25*b)+7.50001e-1)+416 
  109.  g=int(c/30.6001) : mn=g-1 
  110.  dy=c-int(30.6001*g)+fd : yr=b+1899 
  111.  if g > 13.5 then mn=g-13 
  112.  if mn < 2.5 then yr=b+1900 
  113.  if yr < 1 then --yr 
  114.  
  115.  '..return a date string (whole days only) 
  116.  dy$=str$(int(dy)) : if dy < 10 then dy$="0"+right$(dy$,1)
  117.  dy$=right$(dy$,2)
  118.  mn$=str$(int(mn)) : if mn < 10 then mn$="0"+right$(mn$,1)
  119.  mn$=right$(mn$,2)
  120.  yr$=str$(int(yr))
  121.  yr$=right$(yr$,4)
  122.  
  123.  date_and_time = dy$+"-"+mn$+"-"+yr$+"   "+time_from_day_fraction(fd)
  124. END SUB 
  125.  
  126. SUB SINGLE in_range(n)
  127. '..ensure n is in the range 0..360 degrees
  128.   if n<0! then 
  129.     in_range = 360! + (n mod 360!)
  130.   else
  131.     if n>360! then in_range = n mod 360!
  132.   end if
  133. END SUB
  134.  
  135. SUB SINGLE sinh(x)
  136. '..return hyperbolic sine of x
  137.   sinh = (exp(x)-exp(-x))/2!
  138. END SUB
  139.  
  140. SUB JovianEphemeris(SINGLE JDE)
  141. SHARED d,V,M,N,J,A,B,K,R,rr,delta,psi,DE
  142. '..calculate circumstances of Jupiter at JDE
  143.  
  144.   '..days since 2000 January 1, 12h
  145.   d = JDE - 36525.0
  146.  
  147.   '..argument for long-period term in motion of Jupiter    
  148.   V = in_range(172.74 + 0.00111588*d)    
  149.  
  150.   '..mean anomalies of Earth and Jupiter
  151.   M = in_range(357.529 + 0.9856003*d)
  152.   N = 20.02 + 0.0830853*d + 0.329*sin(V/radconv)
  153.  
  154.   '..difference between mean heliocentric 
  155.   '..longitudes of Earth and Jupiter
  156.   J = in_range(66.115 + 0.9025179*d - 0.329*sin(V/radconv))
  157.  
  158.   '..equations of the center of Earth and Jupiter
  159.   A = 1.915*sin(M/radconv) + 0.02*sin((2*M)/radconv)
  160.   B = 5.555*sin(N/radconv) + 0.168*sin((2*N)/radconv)
  161.   K=J+A-B
  162.  
  163.   '..radius vector of Earth
  164.   R = 1.00014 - 0.01671*cos(M/radconv) - 0.00014*cos((2*M)/radconv)
  165.  
  166.   '..radius vector of Jupiter
  167.   rr = 5.20872 - 0.25208*cos(N/radconv) - 0.00611*cos((2*N)/radconv)  
  168.  
  169.   '..Earth-Jupiter distance
  170.   delta = ABS(SQR(rr*rr + R*R - 2*rr*R*cos(K/radconv)))
  171.  
  172.   '..phase angle (Earth-Jupiter-Sun)
  173.   sin_of_psi = (R/delta)*sin(K/radconv)
  174.   psi = sinh(sin_of_psi)*radconv
  175.  
  176.   '..longitudes of central meridian in systems I and II
  177.   w1 = in_range(210.98 + 877.8169088*(d-(delta/173!)) + psi - B)
  178.   w2 = in_range(187.23 + 870.1869088*(d-(delta/173!)) + psi - B)
  179.  
  180.   '..heliocentric longitude
  181.   lambda = 34.35 + 0.083091*d + 0.329*sin((V+B)/radconv)
  182.  
  183.   '..planetocentric declination
  184.   DS = 3.12*sin((lambda+42.8)/radconv)
  185.   DE = DS - 2.22*sin(psi/radconv)*cos((lambda+22!)/radconv)
  186.   DE = DE - 1.3*((rr-delta)/delta)*sin((lambda-100.5)/radconv)  
  187. END SUB
  188.  
  189. SUB AngleFromInfConj
  190. '..calculate angle from inferior conjunction
  191. SHARED d,delta,psi,B
  192. SHARED first_u1,first_u2,first_u3,first_u4
  193. SHARED u1,u2,u3,u4
  194. SHARED G,H
  195.  
  196.   deltaterm = d - (delta/173)
  197.  
  198.   first_u1 = in_range(163.8067 + 203.4058643*deltaterm + psi - B)
  199.   first_u2 = in_range(358.4108 + 101.2916334*deltaterm + psi - B)
  200.   first_u3 = in_range(5.7129 + 50.2345179*deltaterm + psi - B)
  201.   first_u4 = in_range(224.8151 + 21.4879801*deltaterm + psi - B)
  202.  
  203.   '..correct for mutual perturbations
  204.   G = 331.18 + 50.310482*deltaterm
  205.   H = 87.4 + 21.569231*deltaterm
  206.  
  207.   u1 = first_u1 + 0.473*sin((2*(first_u1-first_u2))/radconv)
  208.   u2 = first_u2 + 1.065*sin((2*(first_u2-first_u3))/radconv)
  209.   u3 = first_u3 + 0.165*sin(G/radconv)
  210.   u4 = first_u4 + 0.841*sin(H/radconv)
  211. END SUB
  212.  
  213. SUB DistFromCenter
  214. '..calculate distance of each satellite from
  215. '..center of Jupiter
  216. SHARED first_u1,first_u2,first_u3,first_u4
  217. SHARED r1,r2,r3,r4
  218. SHARED G,H
  219.   
  220.   r1 = 5.9073 - 0.0244*cos((2*(first_u1-first_u2))/radconv)
  221.   r2 = 9.3991 - 0.0882*cos((2*(first_u2-first_u3))/radconv) 
  222.   r3 = 14.9924 - 0.0216*cos(G/radconv)
  223.   r4 = 26.3699 - 0.1935*cos(H/radconv)
  224. END SUB
  225.  
  226. SUB calc_x_y(SHORTINT n)
  227. '..calculate rectangular coordinates 
  228. '..of four satellites
  229. SHARED x,y
  230. SHARED r1,r2,r3,r4
  231. SHARED u1,u2,u3,u4
  232. SHARED DE
  233.  
  234.   case
  235.     n=1 : r=r1:u=u1
  236.     n=2 : r=r2:u=u2
  237.     n=3 : r=r3:u=u3
  238.     n=4 : r=r4:u=u4
  239.   end case 
  240.  
  241.   x(n) = r*sin(u/radconv)
  242.   y(n) = -r*cos(u/radconv)*sin(DE/radconv)   
  243. END SUB
  244.  
  245. SUB LONGINT moving_east(SHORTINT moon)
  246. '..return true if moon is moving east.
  247. SHARED lastx,x  
  248.  
  249.   if lastx(moon) > x(moon) then 
  250.     moving_east = true
  251.   else
  252.     moving_east = false
  253.   end if
  254. END SUB
  255.  
  256. SUB LONGINT in_disk_region(SHORTINT x,SHORTINT y)
  257. '..return true or false according to whether
  258. '..a satellite is in the region of the disk
  259. '..of Jupiter.
  260.  
  261.   if point(x-1,y)=JovianColor and point(x+1,y)=JovianColor then
  262.     in_disk_region = true
  263.   else
  264.     in_disk_region = false
  265.   end if    
  266. END SUB
  267.  
  268. SUB LONGINT behind_disk(SHORTINT xcoord,SHORTINT ycoord,SHORTINT moon)
  269. '..return true or false according to whether
  270. '..a satellite is behind the disk of Jupiter.
  271.  
  272.   if in_disk_region(xcoord,ycoord) and moving_east(moon) then
  273.     behind_disk = true
  274.   else
  275.     behind_disk = false
  276.   end if   
  277. END SUB
  278.  
  279. SUB RemoveMoons
  280. SHARED x,y
  281. SHORTINT moon,xx,yy,colr
  282.  
  283.   '..clear moons
  284.   for moon=1 to 4
  285.     xx = xorigin+x(moon)*xscale
  286.     yy = yorigin-y(moon)*yscale
  287.     
  288.     if in_disk_region(xx,yy) then
  289.       '..in transit across disk or behind it
  290.       colr = JovianColor
  291.     else
  292.       colr = black
  293.     end if
  294.     pset (xx,yy),colr
  295.   next
  296. END SUB
  297.  
  298. SUB ShowJovianSpace
  299. '..display Jupiter and the Galilean satellites
  300. SHARED x,y
  301. SHORTINT xx,yy,moon
  302.  
  303.   '..plot moons
  304.   for moon=1 to 4
  305.     xx = xorigin+x(moon)*xscale
  306.     yy = yorigin-y(moon)*yscale
  307.     if not behind_disk(xx,yy,moon) then pset (xx,yy),moon
  308.   next
  309.  
  310.   '..draw Jupiter
  311.   circle (xorigin,yorigin),radius,JovianColor
  312.   paint (xorigin,yorigin),JovianColor
  313. END SUB
  314.  
  315. SUB DisplayData
  316. '..display Jupiter/satellite data 
  317. CONST startcol=15
  318. SHARED JDE,x,y,rr,delta
  319. SHORTINT moon  
  320.  
  321.   '..Date & Time
  322.   locate 2,10
  323.   color 1,0
  324.   print date_and_time(JDE);" UT"
  325.   
  326.   '..Earth-Jupiter distance  
  327.   locate 4,10
  328.   color 4,0
  329.   print "Earth-Jupiter Distance (AU): ";
  330.   color 5,0
  331.   print delta;"    " 
  332.  
  333.   '..Jupiter's Radius Vector
  334.   locate 5,10
  335.   color 6,0
  336.   print "  Sun-Jupiter Distance (AU): ";
  337.   color 5,0
  338.   print rr;"    "
  339.  
  340.   '..headings
  341.   locate 7,startcol
  342.   color 7,0
  343.   print " Io"
  344.   locate 7,startcol+15
  345.   color 6,0
  346.   print " Europa"
  347.   locate 7,startcol+30
  348.   color 7,0
  349.   print " Ganymede"
  350.   locate 7,startcol+45
  351.   color 6,0
  352.   print " Callisto"
  353.   print
  354.  
  355.   '..satellite's X coordinate
  356.   locate csrlin,10:print "X: ";
  357.   col=startcol
  358.   for moon=1 to 4
  359.     locate csrlin,col
  360.     color moon
  361.     print x(moon);"    ";
  362.     col=col+15 
  363.   next
  364.  
  365.   '..satellite's Y coordinate
  366.   print
  367.   locate csrlin,10:print "Y: ";
  368.   col=startcol
  369.   for moon=1 to 4
  370.     locate csrlin,col
  371.     color moon
  372.     print y(moon);"    ";
  373.     col=col+15  
  374.   next
  375. END SUB
  376.  
  377. '..main
  378. if arg$(1)="?" then 
  379.   print
  380.   print "usage: Jovian date start duration interval"
  381.   print
  382.   print "       where: date is of the form dd mm yyyy"
  383.   print
  384.   print "              start, duration and interval"
  385.   print "              consist of hh mm ss "
  386.   print
  387.   print "              Date and Time are taken to be UT.
  388.   print        
  389.   stop
  390. end if
  391.  
  392. if argcount=12 then  
  393.   dd = val(arg$(1))
  394.   mm = val(arg$(2))
  395.   yy = val(arg$(3))
  396.  
  397.   hrs = val(arg$(4))
  398.   mins = val(arg$(5))
  399.   secs = val(arg$(6))
  400.  
  401.   hr_dur = val(arg$(7))
  402.   min_dur = val(arg$(8))
  403.   sec_dur = val(arg$(9))
  404.  
  405.   hr_int = val(arg$(10))
  406.   min_int = val(arg$(11))
  407.   sec_int = val(arg$(12))
  408. else
  409.   window 1,"Jovian Satellite Simulation",(0,0)-(640,70)
  410.   print
  411.   print "Start Date (UT):"
  412.   print
  413.   input "Enter day   ",dd
  414.   input "Enter month ",mm
  415.   input "Enter year  ",yy
  416.  
  417.   cls
  418.   print
  419.   print "Start Time (UT):"
  420.   print
  421.   input "Enter hours   ",hrs
  422.   input "Enter minutes ",mins
  423.   input "Enter seconds ",secs
  424.  
  425.   cls
  426.   print
  427.   print "Duration:"
  428.   print
  429.   input "Enter hours   ",hr_dur
  430.   input "Enter minutes ",min_dur
  431.   input "Enter seconds ",sec_dur
  432.  
  433.   cls
  434.   print
  435.   print "Interval:"
  436.   print
  437.   input "Enter hours   ",hr_int
  438.   input "Enter minutes ",min_int
  439.   input "Enter seconds ",sec_int
  440.   window close 1
  441.  
  442.   '..An interval of 6 minutes is the finest resolution 
  443.   '..possible with single-precision floating-point math!
  444.   if hr_int=0 and min_int<6 then min_int=6
  445. end if
  446.  
  447. day_fraction = decimal_hours(hrs,mins,secs) / 24!
  448. JDE = JulianDay(dd,mm,yy) + day_fraction
  449.  
  450. duration = decimal_hours(hr_dur,min_dur,sec_dur) / 24!
  451. end_point = JDE+duration
  452.  
  453. interval = decimal_hours(hr_int,min_int,sec_int) / 24!
  454.  
  455. screen 1,640,225,4,2
  456.  
  457. palette 0,0,0,0        '..black
  458. palette 1,1,.73,0    '..orange (Io)
  459. palette 2,1,.87,.73    '..tan (Europa, Jupiter)
  460. palette 3,.93,.2,0    '..fire engine red (Ganymede) 
  461. palette 4,.8,.6,.53    '..brown (Callisto) 
  462. palette 5,1,1,1        '..white
  463. palette 6,0,.93,.87    '..aqua
  464. palette 7,.33,.87,0    '..green
  465. palette 8,0,0,1        '..blue
  466. palette 9,1,1,.13    '..yellow
  467.  
  468. '..border
  469. line (5,5)-(635,195),8,b
  470.  
  471. '..N,S,E,W markers
  472. color 9,0
  473. locate 15,5
  474. print "E"
  475. locate 15,76
  476. print "W"
  477. locate 12,41
  478. print "N"
  479. locate 18,41
  480. print "S"
  481.  
  482. locate 23,3
  483. color 6
  484. print "Hit key/left mouse button."
  485.  
  486. '..initialise lastX array
  487. for moon=1 to 4
  488.   lastx(moon) = -99
  489. next
  490.  
  491. first=true
  492.  
  493. repeat
  494.   JovianEphemeris(JDE)
  495.   AngleFromInfConj
  496.   DistFromCenter
  497.  
  498.   if not first then 
  499.     RemoveMoons
  500.   else 
  501.     first=false
  502.   end if
  503.  
  504.   for moon=1 to 4
  505.     if not first then lastx(moon) = x(moon)
  506.     calc_x_y(moon)
  507.   next
  508.  
  509.   ShowJovianSpace
  510.  
  511.   DisplayData
  512.  
  513.   JDE = JDE + interval
  514. until mouse(0) or inkey$<>"" or JDE > end_point
  515.  
  516. locate 23,3
  517. color 7
  518. print "Press Q key.               "
  519. while ucase$(inkey$)<>"Q":wend
  520.  
  521. screen close 1
  522.