home *** CD-ROM | disk | FTP | other *** search
/ Fresh Fish 5 / FreshFish_July-August1994.bin / useful / dist / dev / lang / ace / prgs / 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.   '..J