home *** CD-ROM | disk | FTP | other *** search
/ The Very Best of Atari Inside / The Very Best of Atari Inside 1.iso / sharew / astronom / jupretro / retsourc.lst < prev   
Encoding:
File List  |  1986-07-20  |  12.8 KB  |  488 lines
  1. Dim Dc_year(3)
  2. R_d=180/Pi   ! Constant to change radians to degrees
  3. D_r=Pi/180   ! Constant to change degrees to radians
  4. '
  5. '
  6. '
  7. Rez=Xbios(4) ! Determine the resolution 1 = medium, 2 = High
  8. '
  9. '   Since the plotting and such are set up for Med. and High rez,
  10. '   don't allow the program to be run in low resolution
  11. '
  12. If Xbios(4)=0 Then
  13.   Alert 1,"Please switch to|Medium resolution!",1,"Ok!",Dummy%
  14.   End
  15. Endif
  16. '
  17. '
  18. ' Start the program !!
  19. '
  20. Repeat
  21.   Hidem  ! Hide the mouse cursor or it'll get in the way
  22.   '
  23.   '  Input the Year of curiosity
  24.   '
  25.   Input " Year? ",Y
  26.   '
  27.   ' The next few lines calculate the Julian date for Jan. 1 of the desired year.
  28.   '
  29.   Y1=Y-1
  30.   M=13
  31.   D=1
  32.   Cls
  33.   A=Int(Y1/100)
  34.   B=2-A+Int(A/4)
  35.   Jd=Int(365.25*Y1)+Int(30.6001*(M+1))+D+1720994.5+B
  36.   '
  37.   ' Set up constants to find the date of opposition
  38.   '
  39.   Stepper=32
  40.   Day=-16
  41.   Flag1=True ! Flag to detect when oppostion occurs
  42.   Old_diff=0
  43.   Gosub Bell
  44.   Cls
  45.   Print At(30,12);"  Working..";
  46.   '
  47.   ' Here we enter the loop to find the date of opposition by
  48.   ' quartering the year first by 32 days, then 16, then 8...
  49.   '
  50.   While Stepper>0.25
  51.     Iter=1
  52.     Flag2=False
  53.     Repeat
  54.       Te=Jd+Day-2415020 ! Number of days since 1/1/1900
  55.       Tc=Te/36525       ! Number of centuries since 1/1/1900
  56.       '                 ! Find position of Sun and Jupiter at above date
  57.       Gosub Position(Tc)
  58.       '                  When Diff = 0 then Jupiter is at opposition
  59.       Diff=Abs(Abs(Ra_sun-Ra_jup)-180)
  60.       '   Check flags to see if near or at opposition
  61.       If Flag1=True And Old_diff<Diff Then
  62.         Flag2=False
  63.       Else
  64.         If Flag1=False And Old_diff<Diff Then
  65.           Flag2=True
  66.         Endif
  67.       Endif
  68.       If Old_diff>Diff Then
  69.         Flag1=False
  70.       Else
  71.         Flag1=True
  72.       Endif
  73.       Old_diff=Diff
  74.       Day=Day+Stepper
  75.       Iter=Iter+1
  76.     Until Flag2
  77.     '
  78.     '  Having detected that we're near opposition cut Stepper in two and
  79.     '  start again and go back sufficient time to find actual time of
  80.     '  opposition!
  81.     '
  82.     Print ".";
  83.     Day=Day-3*Stepper
  84.     Stepper=Stepper/2
  85.     '
  86.     ' when stepper meets condition above leave loop
  87.     '
  88.   Wend
  89.   Print
  90.   Print
  91.   ' since there was a redcuton of the value of Day above increment it
  92.   Day=Day+2*Stepper
  93.   Gosub Bell
  94.   Print
  95.   '
  96.   '  calculate calendar date from Julian date
  97.   '
  98.   Cls
  99.   Gosub Opp_date(Day,Jd)
  100.   Print
  101.   '
  102.   '   Print out results
  103.   '
  104.   If Day=0 Then
  105.     Day=1
  106.   Endif
  107.   If Day<=365 And Day>0 Then
  108.     Print "   The opposition of Jupiter occurred on day ";Day;" of the year ";Y;"."
  109.   Endif
  110.   If Day>365 Then
  111.     Day1=Day-365
  112.     Print "   There was no opposition of Jupiter in the year ";Y
  113.     Y=Y+1
  114.     Print "   but occurred on day ";Day1;" of the following year ";Y
  115.   Endif
  116.   Print
  117.   Ra_sun=Ra_sun/15
  118.   Ra_jup=Ra_jup/15
  119.   Ra_s_min=Trunc(60*Frac(Ra_sun))
  120.   Ra_j_min=Trunc(60*Frac(Ra_jup))
  121.   Ra_sun=Trunc(Ra_sun)
  122.   Ra_jup=Trunc(Ra_jup)
  123.   Dc_s_min=Abs(Trunc(60*Frac(Dc_sun)))
  124.   Dc_j_min=Abs(Trunc(60*Frac(Dc_jup)))
  125.   Dc_sun=Trunc(Dc_sun)
  126.   Dc_jup=Trunc(Dc_jup)
  127.   Print
  128.   Print "   R.A. of the Sun at Opposition = ";Ra_sun;"H ";" ";Ra_s_min;"M"
  129.   Print "                 Dec. of the Sun = ";Dc_sun;Chr$(248);" ";Dc_s_min;"'"
  130.   Print
  131.   Print "   R.A. of Jupiter at Opposition = ";Ra_jup;"H ";" ";Ra_j_min;"M"
  132.   Print "                 Dec. of Jupiter = ";Dc_jup;Chr$(248);" ";Dc_j_min;"'"
  133.   Print
  134.   Input "   Press 'Return' to start plot ";A$
  135.   Gosub Bell
  136.   '
  137.   ' Set up the limits for plotting
  138.   '
  139.   Day_opp=Day
  140.   Dc_year(1)=Dc_jup
  141.   Dc_opp=Dc_jup
  142.   Ra_opp=Ra_jup
  143.   '
  144.   ' Calculate the starting position
  145.   '
  146.   Day_start=Day-160
  147.   Te=Jd+Day_start-2415020
  148.   Tc=Te/36525
  149.   Gosub Position(Tc)
  150.   Ra_start=Ra_jup
  151.   Ra_start1=Ra_start/15
  152.   Ra_start$="R.A. Start = "+Str$(Trunc(Ra_start1))+"H "+Str$(Trunc((Ra_start1-Trunc(Ra_start1))*60))+"M"
  153.   Dc_start=Dc_jup
  154.   Dc_year(2)=Dc_start
  155.   '
  156.   ' Calculate the end position
  157.   '
  158.   Day_end=Day+160
  159.   Te=Jd+Day_end-2415020
  160.   Tc=Te/36525
  161.   Gosub Position(Tc)
  162.   Ra_end=Ra_jup
  163.   Dc_end=Dc_jup
  164.   Ra_end1=Ra_end/15
  165.   Ra_end$="R.A. End = "+Str$(Trunc(Ra_end1))+"H "+Str$(Trunc((Ra_end1-Trunc(Ra_end1))*60))+"M"
  166.   Dc_year(3)=Dc_end
  167.   '
  168.   '  Sort Declinations of start, end and opposition positions to determine
  169.   '  the range of degrees of declination traveled to determine Dc_step.
  170.   '
  171.   Gosub Dc_sort
  172.   '
  173.   '   Since R.A.'s go from 0 to 360 degrees[0 to 24 hours] we have to
  174.   '   determine if the plot crosses from 359 to 1 degrees to keep the plot on
  175.   '   scale.
  176.   '
  177.   Fudge=0  ! Fudge factor if R.A. crosses the 0H boundary
  178.   Ra_flag=0
  179.   If Ra_end<Ra_start Then
  180.     Fudge=360
  181.     Ra_flag=1
  182.   Endif
  183.   Ra_start=Ra_start-Fudge
  184.   Ra_diff=Ra_start-Ra_end
  185.   '
  186.   ' Determine the plotting step value for R.A.
  187.   '
  188.   Ra_step=640/Ra_diff
  189.   Diff=0
  190.   '
  191.   ' clear the screen and start drawing
  192.   '
  193.   Cls
  194.   If Rez=1 Then
  195.     Box 0,7,639,199
  196.   Else
  197.     Box 0,10,639,399
  198.   Endif
  199.   Deftext ,0,2700,6
  200.   Text 630,50*Rez,130,Ra_start$
  201.   Deftext ,,900
  202.   Text 8,150*Rez,130,Ra_end$
  203.   Year$=Str$(Y)
  204.   W_title$="Retrograde Motion of Jupiter for the year "+Year$
  205.   If Rez=2 Then
  206.     Deftext ,0,0,6
  207.   Else
  208.     Deftext ,1,0,4
  209.   Endif
  210.   If Rez=1 Then
  211.     Text 150,5,300,W_title$
  212.   Else
  213.     Text 150,8,300,W_title$
  214.   Endif
  215.   Dc_1=(Dc_start-Dc_bottom)*Dc_step
  216.   '
  217.   '  Draw to the starting place and label  the Dec. limit
  218.   '
  219.   Draw 639,200*Rez-Dc_1
  220.   Draw  To 1,200*Rez-Dc_1
  221.   Draw 639,200*Rez-Dc_1
  222.   Dc_bottom1h=Trunc(Dc_bottom)
  223.   Dc_bottom1m=Abs(Trunc((Dc_bottom-Dc_bottom1h)*60))
  224.   Dc_text$=Str$(Dc_bottom1h)+Chr$(248)+" "+Str$(Dc_bottom1m)+"'"
  225.   Dc_text$="Dec. = "+Dc_text$
  226.   If Dc_1<(100*Rez) Then
  227.     Text 20,(Rez*198-Dc_1),100,Dc_text$
  228.   Else
  229.     Text 20,(Rez*206-Dc_1),100,Dc_text$
  230.   Endif
  231.   '
  232.   '  Here be the plotting routine!!
  233.   '
  234.   For Day_pos=Day_start To Day_end
  235.     Te=Jd+Day_pos-2415020
  236.     Tc=Te/36525
  237.     Gosub Position(Tc)
  238.     If Ra_jup>180 And Ra_flag=1 Then
  239.       Ra_jup=Ra_jup-Fudge
  240.     Endif
  241.     Diff_ra_start=Abs(Ra_start-Ra_jup)
  242.     Ra_place=639-Abs(Diff_ra_start*Ra_step)
  243.     Dc_place=Rez*200-Abs(Dc_jup-Dc_bottom)*Dc_step
  244.     Diff=Ra_jup-O_ra_jup
  245.     Draw  To Ra_place,Dc_place
  246.     '
  247.     '  If at the day of opposition draw a tic mark to show it
  248.     '
  249.     If (Day_opp=Day_pos) Then
  250.       Draw  To Ra_place,Dc_place+5
  251.       Draw  To Ra_place,Dc_place-5
  252.       Draw  To Ra_place,Dc_place
  253.     Endif
  254.     O_ra_jup=Ra_jup
  255.   Next Day_pos
  256.   Dc_jup1=Trunc(Dc_jup*100)/100
  257.   '
  258.   '  Draw the end Dec. limit and label it.
  259.   '
  260.   Draw  To 639,Dc_place
  261.   Dc_jup1h=Trunc(Dc_jup)
  262.   Dc_jup1m=Abs(Trunc((Dc_jup-Dc_jup1h)*60))
  263.   Dc_text$=Str$(Dc_jup1h)+Chr$(248)+" "+Str$(Dc_jup1m)+"'"
  264.   Dc_text$="Dec. = "+Dc_text$
  265.   If Dc_place<100 Then
  266.     Text 500,Dc_place+6,100,Dc_text$
  267.   Else
  268.     Text 500,Dc_place-2,100,Dc_text$
  269.   Endif
  270.   Gosub Bell
  271.   What_now:
  272.   Repeat
  273.   Until Mousek=0
  274.   Showm
  275.   Alert 2,"Retrograde Drawing|    Complete| Press any Key|  to Continue",1,"Look|Print|QUIT ",Butn
  276.   Hidem
  277.   If Butn=2 Then
  278.     Hardcopy
  279.     Goto What_now
  280.   Endif
  281.   If Butn=3
  282.     Goto Drop_out
  283.   Else
  284.     Gosub Any_key
  285.     Goto What_now
  286.   Endif
  287.   '
  288.   '  Asks if you want to do it again
  289.   '
  290.   Drop_out:
  291.   Input "Do you wish another year's plot (Y/N) ";Yes$
  292. Until ((Yes$="n") Or (Yes$="N"))
  293. '      Return mouse cursor
  294. Showm
  295. End
  296. '
  297. '   The following Procedure calculates the positions of the Sun and Jupiter
  298. '   for the given Julian Date.
  299. '
  300. Procedure Position(Tx)
  301.   '  The eccentricity of the Earth's orbit's orbit's orbit's orbit's orbitun-E_sun2
  302.     E_sun2=E_sun
  303.   Wend
  304.   ' Radius vector of the Sun
  305.   R_sun=1.0000002*(1-Ecc_earth*Cos(E_sun*D_r))
  306.   ' the true anomaly of the Sun
  307.   V_sun=2*Atn(Sqr((1+Ecc_earth)/(1-Ecc_earth))*Tan(E_sun*D_r/2))*R_d
  308.   ' the true logitude of the Sun
  309.   Sun_true_lng=L_sun+V_sun-M_sun
  310.   '  Correction for nutation and abberation
  311.   Omega_sun=259.18-1934.142*Tx
  312.   ' Obiliquity of the Ecliptic
  313.   Obilquity=Obilquity+0.00256*Cos(Omega_sun*D_r)
  314.   Temp1=Cos(Obilquity*D_r)*Sin(Sun_true_lng*D_r)
  315.   Temp2=Cos(Sun_true_lng*D_r)
  316.   '
  317.   ' Now calculate the right ascension of the Sun and put it in
  318.   ' the right quadrent
  319.   '
  320.   Ra_sun=Atn(Temp1/Temp2)
  321.   If Temp2>0 And Temp1<0 Then
  322.     Ra_sun=Ra_sun+2*Pi
  323.   Endif
  324.   If Temp2<0 Then
  325.     Ra_sun=Ra_sun+Pi
  326.   Endif
  327.   '
  328.   ' Calculate the declination position of the Sun
  329.   '
  330.   Temp1=Sin(Obilquity*D_r)*Sin(Sun_true_lng*D_r)
  331.   Temp2=Sqr(1-Temp1*Temp1)
  332.   Dc_sun=Atn(Temp1/Temp2)*R_d
  333.   Ra_sun=Ra_sun*R_d
  334.   '
  335.   '  Now repeat the above for Jupiter
  336.   '
  337.   L_jup=238.049257+3036.301986*Tx+0.0003347*Tx*Tx-1.65E-06*Tx*Tx*Tx
  338.   A_jup=5.202561
  339.   Ecc_jup=0.04833475+0.00016418*Tx-4.676E-07*Tx*Tx-1.7E-09*Tx*Tx*Tx
  340.   I_jup=1.308736-0.0056961*Tx+3.9E-06*Tx*Tx
  341.   W_jup=273.277558+0.5994317*Tx+0.00070405*Tx*Tx+5.08E-06*Tx*Tx*Tx
  342.   Omega_jup=99.443414+1.01053*Tx+0.00035222*Tx*Tx-8.51E-06*Tx*Tx*Tx
  343.   M_jup=225.322833+3034.69202*Tx-0.000722*Tx*Tx
  344.   E_diff=1
  345.   E_jup2=M_jup
  346.   While Abs(E_diff)>1.0E-06
  347.     E_jup=M_jup+Ecc_jup*R_d*Sin(E_jup2*D_r)
  348.     E_diff=E_jup-E_jup2
  349.     E_jup2=E_jup
  350.   Wend
  351.   Temp1=Sin(E_jup*D_r/2)
  352.   Temp2=Cos(E_jup*D_r/2)
  353.   Temp3=Sqr((1+Ecc_jup)/(1-Ecc_jup))
  354.   V_jup=Atn(Temp3*Temp1/Temp2)
  355.   If (Temp1*Temp3)<0 And Temp2>0 Then
  356.     V_jup=V_jup+2*Pi
  357.   Endif
  358.   If Temp2<0 Then
  359.     V_jup=V_jup+Pi
  360.   Endif
  361.   V_jup=2*V_jup*R_d
  362.   R_jup=A_jup*(1-Ecc_jup*Cos(E_jup*D_r))
  363.   U_jup=L_jup+V_jup-M_jup-Omega_jup
  364.   Temp1=Cos(I_jup*D_r)*Sin(U_jup*D_r)
  365.   Temp2=Cos(U_jup*D_r)
  366.   Long_jup=Atn(Temp1/Temp2)
  367.   If Temp2>0 And Temp1<0 Then
  368.     Long_jup=Long_jup+2*Pi
  369.   Endif
  370.   If Temp2<0 Then
  371.     Long_jup=Long_jup+Pi
  372.   Endif
  373.   Long_jup=Long_jup*R_d
  374.   Long_jup=Long_jup+Omega_jup
  375.   Temp1=Sin(U_jup*D_r)*Sin(I_jup*D_r)
  376.   Temp2=Sqr(1-Temp1*Temp1)
  377.   B_jup=Atn(Temp1/Temp2)*R_d
  378.   Geo_sun_lng=Long_jup-Sun_true_lng
  379.   Temp1=R_jup*Cos(B_jup*D_r)*Sin((Geo_sun_lng)*D_r)
  380.   Temp2=R_jup*Cos(B_jup*D_r)*Cos((Geo_sun_lng)*D_r)+R_sun
  381.   Long_jup_geo=Atn(Temp1/Temp2)
  382.   If Temp2>0 And Temp1<0 Then
  383.     Long_jup_geo=Long_jup_geo+2*Pi
  384.   Endif
  385.   If Temp2<0 Then
  386.     Long_jup_geo=Long_jup_geo+Pi
  387.   Endif
  388.   Long_jup_geo=Long_jup_geo*R_d
  389.   Temp_geo=Long_jup_geo
  390.   Long_jup_geo=Long_jup_geo+Sun_true_lng
  391.   Delta_jup=Sqr(R_sun*R_sun+R_jup*R_jup+2*R_jup*R_sun*Cos(B_jup*D_r)*Cos((Geo_sun_lng)*D_r))
  392.   Temp1=R_jup*Sin(B_jup*D_r)/Delta_jup
  393.   Temp2=Sqr(1-Temp1*Temp1)
  394.   Beta_jup=Atn(Temp1/Temp2)*R_d
  395.   Temp1=Cos(Beta_jup*D_r)*Cos(Temp_geo*D_r)
  396.   Temp2=Sqr(1-Temp1*Temp1)
  397.   Elong_jup=Atn(Temp2/Temp1)
  398.   Temp1=Sin(Long_jup_geo*D_r)*Cos(Obilquity*D_r)-Tan(Beta_jup*D_r)*Sin(Obilquity*D_r)
  399.   Temp2=Cos(Long_jup_geo*D_r)
  400.   Ra_jup=Atn(Temp1/Temp2)
  401.   If Temp2>0 And Temp1<0 Then
  402.     Ra_jup=Ra_jup+2*Pi
  403.   Endif
  404.   If Temp2<0 Then
  405.     Ra_jup=Ra_jup+Pi
  406.   Endif
  407.   Ra_jup=Ra_jup*R_d
  408.   Temp1=Sin(Beta_jup*D_r)*Cos(Obilquity*D_r)+Cos(Beta_jup*D_r)*Sin(Obilquity*D_r)*Sin(Long_jup_geo*D_r)
  409.   Temp2=Sqr(1-Temp1*Temp1)
  410.   Dc_jup=Atn(Temp1/Temp2)*R_d
  411. Return
  412. '
  413. ' Procedure to sort the three important declinations of Jupiter, start,end,
  414. ' and opposition
  415. '
  416. Procedure Dc_sort
  417.   Switch=True
  418.   While Switch=True
  419.     For I%=1 To 2
  420.       If Dc_year(I%)<Dc_year(I%+1) Then
  421.         Dc_temp=Dc_year(I%+1)
  422.         Dc_year(I%+1)=Dc_year(I%)
  423.         Dc_year(I%)=Dc_temp
  424.         Switch=True
  425.       Else
  426.         Switch=False
  427.       Endif
  428.     Next I%
  429.   Wend
  430.   '
  431.   ' Having sorted the Dec.'s calculate the Dec. range and declination
  432.   ' plotting step
  433.   '
  434.   Dc_bottom=Dc_year(3)-0.5
  435.   Dc_range=Abs(Dc_year(1)-Dc_year(3))+1
  436.   Dc_step=200*Rez/Dc_range
  437. Return
  438. '
  439. ' procedure to calculate the calendar date form the Julian Date
  440. '
  441. Procedure Opp_date(Day1,Jd1)
  442.   Jd1=Jd1+Day1+0.5
  443.   Z=Int(Jd1)
  444.   F=Frac(Jd1)
  445.   If Z<2299161 Then
  446.     A=Z
  447.   Else
  448.     Alpha=Int((Z-1867216.25)/36524.25)
  449.     A=Z+1+Alpha-Int(Alpha/4)
  450.   Endif
  451.   B=A+1524
  452.   C=Int((B-122.1)/365.25)
  453.   D=Int(365.25*C)
  454.   E=Int((B-D)/30.6001)
  455.   Day_of_month=B-D-Int(30.6001*E)+F-1
  456.   If E<13.5 Then
  457.     Month=E-1
  458.   Endif
  459.   If E>13.5
  460.     Month=E-13
  461.   Endif
  462.   If Month>2.5 Then
  463.     Year=C-4716
  464.   Endif
  465.   If Month<2.5 Then
  466.     Year=C-4715
  467.   Endif
  468.   Print " >>   Date of Opposition = ";Month;"/";Day_of_month;"/";Year
  469. Return
  470. '
  471. '  Procedure for a bell sound!!!
  472. '
  473. Procedure Bell
  474.   For Y99=1 To 2
  475.     Sound 1,15,1,7
  476.     Wave 1,1,9,15000,100
  477.   Next Y99
  478. Return
  479. '
  480. Procedure Any_key
  481.   Repeat
  482.     Z$=Inkey$
  483.     If Mousek<>0
  484.       Z$=""
  485.     Endif
  486.   Until Z$<>""
  487. Return
  488.