home *** CD-ROM | disk | FTP | other *** search
/ Phoenix CD 2.0 / Phoenix_CD.cdr / 02a / orbsrc14.zip / ORBP.BAS < prev    next >
BASIC Source File  |  1987-11-16  |  24KB  |  724 lines
  1.     '= = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
  2.     ' ORBP.BAS
  3.     ' Compile as .TBC file
  4.     '
  5.     '                              O R B P
  6.     '             ORBITAL PREDICTION PROGRAM WITH TIME WINDOWS
  7.     '                    for the IBM PC & Turbo BASIC
  8.     '                             11/1/1987
  9.     '871116-1
  10.     'Derived from an original publication in AMSAT and ORBIT magazine
  11.     ' by Dr. Thomas A. Clark (W3IWI).  See ORBIT magazine issue #6, 4/81.
  12.     '
  13.     'Converted the program to Turbo BASIC.  Elements are read in from the
  14.     ' file KEPLER.ORB, which is generated directly from AMSAT's Orbital
  15.     ' Element Bulletin contained in a file named ELEMENTS (as read from
  16.     ' a BBS).  The routine which converts ELEMENTS to KEPLER.ORB is
  17.     ' ORBUPDAT.EXE.
  18.     '
  19.     'Station data is now contained in a file named ORBS.CFG.  GMST is
  20.     ' now calculated.  Sidereal time tables no longer needed.
  21.     '      W0PN, 8/14/87
  22.     '= = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
  23.     '
  24.     'Customization of Station data is accomplished by modifying the
  25.     ' file named ORBS.CFG...
  26.     ' Station DATA items are: <call letters>, <Lat>, <W. Lon (in degrees)>,
  27.     ' <Height (meters)>, <min elevation to print in degrees>, <CR/LF>.
  28.  
  29.     ' Satellite records in KEPLER.ORB contain:  <Satellite name>, <Catalog #>,
  30.     ' <epoch [year + Julian day # + decimal portion of day]>,
  31.     ' <Element set number (NASA)>,
  32.     ' <epoch [inclination],  [Right Ascension of Ascending Node (RAAN)]>,
  33.     ' <epoch [eccentricity], [Argument of Perigee],>
  34.     ' <epoch [Mean Anomaly], [Mean Motion],  [Decay Rate]>,
  35.     ' <epoch [Orbit number]>
  36.     ' <F1 Beacon freq in MHz>(optional),
  37.     ' <F2 Beacon freq in MHz>(optional), <CR/LF>.
  38.     '
  39.     '= = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
  40.  
  41.      DEFDBL A-Z
  42.      DIM T$(20), S$(40), I$(40), C(3,2)
  43.  
  44.      CLS
  45.      P% = 1                                       'Page counter zeroed
  46.      LN0% = 60        'LNO% = Lines printed per page, =24 for terminal
  47.      ln%  = 0         'ln% = number of lines printed this page
  48.      dimflag = 0      'Set to 0 to force initialization of DOWTBL
  49.  
  50.      DEF FNT$(D)=chr$(48+fix(D/10))+chr$(48+D-10*fix(D/10))  '2 char string
  51.      DEF FNI(D)=D
  52.  
  53.      PRINT
  54.      color 0,7
  55.      PRINT"= = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
  56.      PRINT"=                                                                   =
  57.      PRINT"=                             O R B S                               =
  58.      PRINT"=                                                                   =
  59.      PRINT"=          ORBITAL PREDICTION PROGRAM WITH TIME WINDOWS             =
  60.      PRINT"=                            11/1/1987                              =
  61.      PRINT"=                                                  Ron Dunbar W0PN  =
  62.      PRINT"= = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
  63.      color 7,0
  64.      PRINT
  65.      PRINT
  66.      PRINT "   This program prints tracking data for earth-orbiting satellites"
  67.      PRINT
  68.      PRINT
  69.      PRINT tab(23);
  70.      color 0,7
  71.      print "SATELLITE SELECTION MENU"
  72.      color 7,0
  73.      PRINT
  74.      I%=0
  75.  
  76.      '- - - - - Select parameters for satellite of interest - - - - -
  77.  
  78.      open "KEPLER.ORB" for INPUT as #1        'read Keplerians from disk
  79.      WHILE not eof(1)
  80.     input# 1, S$,ident,d3#,S0,D,D,D,D,D,D,D,D,D,D  'Ds are dummies
  81.     I$ = str$(ident)            'convert for compatibility
  82.         julian# = fix(d3#)
  83.         gosub julian                'convert the date
  84.         I%=I%+1
  85.         print tab(10);STR$(I%)+".";tab(16);S$;tab(34);"ID ="I$;_
  86.                  "   Set";S0;tab(58);dt$    'display satellite selections
  87.      WEND
  88.     close #1
  89.  
  90.      '- - - now choose specific satellite - - -
  91.  
  92. selsat:
  93.      SP$="SELECT SATELLITE"
  94.      PRINT
  95.      PRINT tab(35-(LEN(SP$)/2));
  96.      color 0,7
  97.      print SP$;
  98.      INPUT; J%
  99.      color 7,0
  100.      print
  101.      PRINT
  102.      IF J%<1 OR J%>I% THEN selsat               'Hiccup if out of bounds
  103.  
  104.      OPEN "KEPLER.ORB" for INPUT as #1
  105.      FOR I%=1 TO J%
  106.     input# 1,S$,ident,d3#,S0,incl,raan#,E0,argp,M0,mmot,drag#,K0,F1,F2
  107.     I$ = str$(ident)            'convert for compatibility
  108.      NEXT I%
  109.  
  110.      S$ = ucase$(S$)                'make sat name upper case
  111.      CLOSE #1
  112.  
  113.      'Convert d3# to Y3,d3#
  114.      y3 = fix(d3#/1000)
  115.      d3# = d3# - (y3*1000)
  116.         julian# = fix(d3#)
  117.         gosub julian                'convert the date
  118.  
  119.      cls
  120.      locate 5,1
  121.      print tab(35 - (len(S$) / 2));
  122.      color 0,7
  123.      print " " ; S$ ; " ";
  124.      color 7,0
  125.      print
  126.      print
  127.      print "      Beacon frequencies (MHz) are [F1]";
  128.      color 0,7
  129.      print using "####.###";F1;
  130.      color 7,0
  131.      print " and [F2]";
  132.      color 0,7
  133.      print using "####.###";F2;
  134.      color 7,0
  135.      print
  136.      print
  137.      print "(Respond to this query with <ENTER> if F2 is NOT to be changed..)"
  138.      INPUT "Change beacon frequency F2 to ";D9
  139.      IF D9 > 0 THEN F2 = D9
  140.     color 7,0
  141.      print
  142.  
  143.      'Set starting date/time, duration, stepsize and initialize
  144.  
  145.      INPUT "Start: year = ";yr
  146.      yr=yr/100
  147.      yr=fix(100*(yr-fix(yr))+.1)        'allows either YYYY or YY input
  148.      IF yr/4=fix(yr/4) THEN F9%=1 ELSE F9%=0    'LEAP-YEAR FLAG
  149.      IF (yr+1900)/400 = fix((yr+1900)/400) then f9% = 0  'Every 400th NOT leap!
  150.  
  151.      INPUT "Month(1-12) = ";M
  152.      INPUT "        Day = ";D
  153.  
  154.      mon  = m
  155.      day  = d
  156.      year = yr
  157.      gosub datecalc        'get all the dates you will ever need!
  158.  
  159.      daynr = doy
  160.      K79 = daynr                                      'k79= old day#
  161.      D8 = doy                        'Day Of Year
  162.      DW$ = dow$
  163.  
  164.      PRINT  "Julian date";julian#
  165.  
  166.      color 0,7
  167.      IF Y3 = yr THEN getcfg ELSE PRINT "ELEMENTS NOT FROM CURRENT YEAR"
  168.      print " Cannot continue."
  169.      color 7,0
  170.      STOP
  171.  
  172. getcfg:
  173.      color 7,0
  174.      open "ORBS.CFG" for input as #1
  175.      input# 1,C$,L9,wlong,obsrvralt,minelev,sflag    'Read station data
  176.      close #1
  177.  
  178.      print
  179. 2280 INPUT"  Start: UTC Hours= ";H
  180.      IF H >= 24 THEN 2280
  181.  
  182.      INPUT"          Minutes = ";M
  183.      Stime# = D8 + H / 24 + M / 1440        'Stime# = Start Epoch
  184.      T$ = T$ + FNT$(H) + FNT$(M) + ":00"        'For printing
  185.  
  186.      INPUT"  Duration: Hours = ";H1
  187.      INPUT"           Minutes= ";M1
  188.      Etime# = Stime# + H1/24 + M1/1440          'Etime# = End Epoch
  189.  
  190.      INPUT"Time step: Minutes= ";M2
  191.      stepamt# = M2/1440                            'stepamt# = Time step, minutes
  192.      PRINT
  193.  
  194. 2450 INPUT "Print only specific time window (y/N)";N$
  195.      n$ = ucase$(n$)                'convert to uppercase
  196.      IF N$ <> "Y" THEN rdphys                 'not specific window
  197.  
  198. 2460 INPUT "Start of window:   Hours = ";SWH
  199.      IF SWH >= 24 THEN 2460
  200.  
  201. 2470 INPUT"                  Minutes= ";SWM
  202.      if swm > 59 then 2470            'Do it right!
  203.      sw# = SWH/24 + SWM/1440
  204.      SW$ = FNT$(SWH) + FNT$(SWM) + ":00"        'start window for printing
  205.  
  206. 2520 INPUT"  End of window:   Hours = ";EWH
  207.      IF EWH >= 24 THEN 2520
  208.  
  209. 2540 INPUT"                  Minutes= ";EWM
  210.      if ewm > 59 then 2540            'Do it right!
  211.      ew# = EWH/24 + EWM/1440
  212.      EW$ = FNT$(EWH) + FNT$(EWM) + ":00"        'end window for printing
  213.      PRINT
  214.      if ew# < sw# then ew# = ew# + 1.0        'stop time is in next day
  215.      swd# = sw# + fix(Stime#)            'set 1st start window time
  216.      ewd# = ew# + fix(Stime#)            'set 1st window stop time
  217.  
  218. rdphys:
  219.      RESTORE 6460
  220.  
  221.      READ P1,C,R0,F,G0,G1                'Get physical constants
  222.      P2 = 2 * P1
  223.      P0 = P1/180
  224.      F  = 1/F
  225.      g2 = gmst#                    'got gmst from datesubs
  226.  
  227.      GOSUB initobsv                            'Init Observer's vector
  228.  
  229.      orbs$ = "*** O R B S ***"
  230.      lprint tab(40-(LEN(orbs$)/2));orbs$    'center the print line
  231.      lprint
  232.      ln% = ln% + 2
  233.  
  234.      GOSUB inithdr                           'Print initial header line
  235.  
  236.      T0# = d3#                            'Epoch start
  237.  
  238.      SET$ = "Element set" + STR$(S0)
  239.      lprint tab(31-(LEN(SET$)/2));SET$;" is";
  240.      lprint using "###.#";(Stime#-T0#);
  241.      lprint " days old"
  242.  
  243.      PRINT "Listing covers the period ";
  244.      PRINT using "###.####";Stime#;
  245.      PRINT " to ";
  246.      PRINT using "###.####";Etime#;
  247.      PRINT " or";
  248.      PRINT using "####.##";(Etime# - Stime#) * 24;
  249.      PRINT " hrs."
  250.      locate 1,1
  251.      print "                                        ";
  252.      print "                                        "
  253.      locate 1,49
  254.      print "Julian search time:"
  255.  
  256.      lprint tab(15) "Listing covers the period ";
  257.      lprint using "###.####";Stime#;
  258.      lprint " to ";
  259.      lprint using "###.####";Etime#;
  260.      lprint " or";
  261.      lprint using "####.##";(Etime# - Stime#) * 24;
  262.      lprint " hrs."
  263.      lprint tab(15) "This listing run at ";TIME$;" on ";DATE$+"."
  264.      lprint
  265.  
  266.      lprint
  267.      lprint tab(15) "PARAMETER";tab(35);"REFERENCE";tab(55);"STARTING"
  268.      lprint tab(15) "---------";tab(35);"---------";tab(55);"--------"
  269.      I = I + 1
  270.  
  271.      time# = Stime#
  272.      GOSUB ncsub        'Update elements to current epoch
  273.      GOSUB nmsub            'Evaluate Mean Anomaly, calc orbit number
  274.  
  275.      lprint tab(15) "Catalog number";tab(35);I$
  276.      lprint tab(15) "Epoch        ";
  277.      lprint tab(35);using "#####.########";(Y3 * 1000) + T0#;_
  278.                         tab(55);(yr * 1000) + Stime#
  279.      lprint tab(15) "Inclination ";
  280.      lprint tab(35);using "###.####";incl ;
  281.      lprint tab(55);" [ No change ]"
  282.      lprint tab(15) "R. A. A. N. ";
  283.      lprint tab(35);using "###.####";raan# ;tab(55);traan#
  284.      lprint tab(15) "Eccentricity ";
  285.      lprint tab(35);using "#.########";E0;
  286.      lprint tab(55);" [ No change ]"
  287.      lprint tab(15) "Arg. Perigee";
  288.      lprint tab(35);using "###.####";argp ;tab(55);targp
  289.      lprint tab(15) "Mean Anomaly ";
  290.      lprint tab(35);using "###.####";M0;tab(55);manom/P0
  291.      lprint tab(15) "Mean Motion ";
  292.      lprint tab(35);using "##.########";mmot ;tab(55);tmmot
  293.      lprint tab(15) "Drag Correction";tab(35);drag#;tab(55);" [ No change ]"
  294.      lprint tab(15) "Orbit number";tab(35);K0;tab(55);orbnr
  295.      lprint tab(15) "S.M.A.,(km)";
  296.      lprint tab(35);using "#####.####";sma;tab(55);tsma
  297.      lprint tab(15) "Perigee Ht, km";
  298.      lprint tab(35);using "#####.####";sma*(1-E0)-R0;tab(55);tsma*(1-E0)-R0
  299.      lprint tab(15) "Apogee  Ht, km";
  300.      lprint tab(35);using "#####.####";sma*(1+E0)-R0;tab(55);tsma*(1+E0)-R0
  301.      if f1 > 0 or f2 > 0 then lprint tab(15);"Freq (MHz) ";
  302.      if f1 > 0 then lprint tab(35);f1;" (F1)";
  303.      if f2 > 0 then lprint tab(55);f2;" (F2)";
  304.      lprint
  305.  
  306.      IF N$ <> "Y" THEN skpprt        'not specific window
  307.      lprint tab(13)"*************************************************************"
  308.      lprint tab(15)"Printing only the window between ";SW$;" and ";EW$;" UTC."
  309.      LN% = LN% + 3     'update line count
  310.      lprint tab(13)"*************************************************************"
  311.  
  312. skpprt:
  313.      K9 = 8.999999E+09
  314.      K8 = 8.999999E+09
  315.      GOSUB utchdr                'print UTC, etc header
  316.      LN% = LN% + 20
  317.      dwc# = fix(Stime#)
  318.  
  319.     '= = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
  320.      '- - - Here follows the actual computation loop - - -
  321.  
  322.      time# = Stime#        'Can't use FOR/NEXT with double precision variables.
  323.      d0 = 3           'disallow treset on first pass
  324.  
  325. comploop:
  326.      daynr = fix(time#)            'RRD.. 9/22/87
  327.      IF N$ <> "Y" THEN inwindow        'Not in window mode.. prt all visible
  328.      if time# > ewd# then        'if so, time to update both swd, ewd
  329.         swd# = swd# + 1.0        'bump to next day
  330.         ewd# = ewd# + 1.0        'bump to next day
  331.      end if
  332.  
  333.      IF NOT (time# >= swd# AND time# <= ewd#) THEN stepadj 'not in window..skip
  334.  
  335. inwindow:
  336.      GOSUB nmsub                        'Evaluate Mean Anomaly, calc orbnr
  337.      IF orbnr = k9 THEN sameorb        'orbnr = Orbit number
  338.      GOSUB ncsub                   'Update elements to current epoch
  339.      k8 = 8.999999E+09            'set day number to force new day prt
  340.      k9 = 8.999999E+09                  'set orb number to force new orb prt
  341.  
  342. sameorb:
  343.      GOSUB nkmsub            'Solve Kepler's equation given MA
  344.      GOSUB nxtsub            'Calc range, velocity, az, el, SSP
  345.  
  346.      d = satelev - minelev                'check elevation
  347.      if d < 0 then quickstep        'below horizon; accelerate stepping
  348.      if d0 = 0 then treset  else d0 = 2        'if quickstepping, reset time
  349.                                                 ' because satellite has risen
  350.      IF k9 = orbnr AND daynr = k8 THEN nochg    'If no chg in orb or day, skip
  351.      IF k9 = orbnr THEN chkpage ELSE k9 = orbnr        'If just new day, print
  352.      GOSUB neworbit                 'new orbit heading
  353.      GOTO  chkpage                'skip around next two routines
  354.  
  355. '- - - - - - - - - - - - - - - - - - - - -
  356. treset:                    'This code resets time to be in sync
  357.      time# = time# - dsave
  358.      steps = fix((time# - fix(time#)) / stepamt#)
  359.      time# = fix(time#) + (stepamt# * steps)
  360.      if time# < Stime# then time# = Stime#    'Don't back up too far
  361.      d0 = 1
  362.      locate 1,70
  363.     color 0,7
  364.      print using "####.###";time#    'show julian time upper right screen
  365.     color 7,0
  366.      goto comploop            'recalc & test new step
  367. '- - - - -
  368. quickstep:                'This algorithm speeds the stepping
  369.      if d0 = 1 then stepadj  else d = range * d^2 * 1e-09
  370.      d0 = 0
  371.      if d > 0.2/mmot  then time# = time# + 0.2/mmot  else time# = time# + d
  372.      dsave = d            'save quickstep step amount for Treset
  373.      goto stepadj
  374.                 'D0 = 0 :below horizon, in Quickstep mode
  375.                 'D0 = 1 :time has just been re-synced
  376.                 'D0 = 2 :above horizon, time in step
  377. '- - - - - - - - - - - - - - - - - - - - -
  378. chkpage:
  379.      IF (LN0% - LN%) < 5 THEN GOSUB newpage ELSE lprint;    'need new page?
  380.      IF LN% < 5 THEN 5740 ELSE lprint;
  381.      IF fix(time#) = dwc# THEN 5120
  382.  
  383.      dwc# = fix(time#)
  384.      julian# = fix((yr * 1000) + dwc#)
  385.  
  386.      gosub julian            'get all the dates again
  387.      daynr = fix(time#)            'daynr = Day number (integer)
  388.  
  389. 5120 lprint tab(5); "- - - Orbit #";orbnr;"- - - Day #";julian#;_
  390.                     "- - - ";dow$;",";dt$;" - - -"
  391.      LN% = LN%+1
  392.  
  393. nochg:
  394.      K8 = daynr
  395.      T4 = time# - daynr
  396.      S4 = fix(T4 * 86400! + .5)
  397.      H4 = fix(S4/3600 + .000001)
  398.      M4 = fix((S4 - H4 * 3600) / 60 + .000001)
  399.      S4 = S4 - 3600 * H4 - 60 * M4
  400.      T$ = FNT$(H4) + FNT$(M4) + ":" + FNT$(S4)      'Printable time string
  401.      F9 = -F1 * 1000000! * velocity / C             '=Doppler F1 * VELOCITY/C
  402.      F10 = -F2 * 1000000! * velocity / C            '=Doppler F2 * VELOCITY/C
  403.      IF LN%> = LN0% THEN GOSUB 5740
  404.  
  405.      LN% = LN% + 1
  406.      lprint tab(4)T$;               'print the detail line
  407.      lprint using "#####";FNI(az),
  408.      lprint using "####";FNI(satelev),
  409.      lprint using "+########";FNI(F9),
  410.      lprint using "+########";FNI(F10),
  411.      lprint using "#######";FNI(range),
  412.      lprint using "######";FNI(R-R0),
  413.      lprint using "#####.#";FNI(L5),
  414.      lprint using "####.#";FNI(W5),
  415.      lprint using "#####";M9;
  416.      lprint "   ";
  417.      IF F2 = 0 THEN f2skip
  418.  
  419.      FTK = (F2 * velocity / C) + F2            'print tracking freq
  420.      lprint using "####.###";FTK;
  421.  
  422. f2skip:
  423.      lprint
  424.  
  425. stepadj:
  426.      time# = time# + stepamt#           'Update time by [stepamt#]
  427.      locate  1,70
  428.     color 0,7
  429.      print using "####.###";time#;            'show julian day/time
  430.     color 7,0
  431.      if satelev < minelev then visible = 0 else visible = 1
  432.      print using "##";visible            '=1 if above min elevation
  433.  
  434.  IF time# < Etime# THEN GOTO comploop               'This replaces NEXT
  435.  
  436.     '= = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
  437.  
  438.      PRINT "END OF ";S$
  439.      lprint
  440.      lprint "End of ";S$;" listing."
  441.      lprint chr$(12);
  442.      PRINT chr$(7)
  443.      cls
  444.      locate 12,5
  445.      print "Thanks for using the   O R B S   tracking program..  Ron, W0PN"
  446.      print
  447.      print
  448.      RUN "orbs.exe"            'EXIT to menu program
  449. '- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  450. newpage:
  451.      lprint                    'start new page
  452.      LN% = LN% + 1
  453.      IF (LN0% - LN%) <= 0 THEN 5740 ELSE newpage
  454.  
  455. 5740 LN% = 1
  456.      GOSUB prthdr                          'Print header on new page
  457.  
  458. neworbit:
  459.      IF (LN0% - LN%) <= 0 THEN 5740
  460.      IF (LN0% - LN%) <  8 THEN newpage
  461.      IF LN% > 6 THEN 6160        'RETURN
  462.  
  463.      '- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  464. utchdr:
  465.      lprint
  466.      lprint tab(4)" U T C    AZ  EL ";
  467.      lprint using "####.###";F1;
  468.      lprint " ";
  469.      lprint using "####.###";F2;
  470.      lprint " ";
  471.      lprint " RANGE HEIGHT   LAT  LONG PHASE  ";
  472.      IF F2 = 0 THEN 5960
  473.  
  474.      lprint "F2 TRACK";
  475. 5960 lprint
  476.      LN% = LN% + 2
  477.      lprint tab(4)"hhmm:ss  deg deg  ";
  478.      lprint "DPLR Hz  DPLR Hz";
  479.      lprint "    km    km     deg  deg  <256>  ";
  480.      IF F2=0 THEN 6100 else lprint using "####.###";F2;
  481. 6100 lprint
  482.      lprint
  483.      LN% = LN% + 2
  484. 6160 RETURN
  485.  
  486.      '- - - - - Page header subroutine - - - - -
  487.  
  488. prthdr:
  489.      IF (LN0% - LN%) = LN0% - 1  THEN 6220  ELSE newpage
  490. 6220 P% = P% + 1
  491.      lprint chr$(12)            'printer TOF (Top Of Forms)
  492.  
  493. inithdr:
  494.      SHDR$=S$+" Orbital Predictions"
  495.      lprint tab(40-(LEN(SHDR$)/2)); SHDR$; tab(70); "Page #"; P%
  496.      lprint
  497.      lprint tab(22);"for ";C$;" at Lat:";
  498.      lprint using "###.###"; L9;
  499.      lprint "  Lon:";
  500.      lprint using "###.###"; wlong
  501.      LN% = LN% + 4
  502.      lprint
  503.      RETURN
  504.      '- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  505.  
  506.      '- - - Physical and numerical constants - - -
  507.  
  508.      'Pi, Velocity of Light, Earth Radius, L/Earth flattening coefficient
  509. 6460 DATA   3.1415926535,   2.997925E5,   6378.160,   298.25
  510.  
  511.      'GM of Earth in (orbits/day)^2/KM^3 and Siderial/Solar time rate ratio
  512.      DATA   7.5369793d13,   1.0027379093
  513.  
  514.      '- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  515.  
  516.      '- - - - - Orbit determination routines - - - - -
  517.      'FNA
  518.      'Calculates inverse tangent in proper quadrant ala Fortran ATNZ
  519.      'cases lying in quadrants 2 and 3
  520.  
  521. nasub:
  522.      IF DX < 0 THEN Q2orQ3
  523.      IF DX > 0 THEN Q1test
  524.  
  525.      'The two cases for DX = 0
  526.      IF DY >= 0 THEN deg90
  527.      D = 3 * P1/2
  528.      RETURN
  529.  
  530.      'Cases lying in quadrants 1 and 4
  531. Q1test:
  532.      IF DY >= 0 then Q1 ELSE Q4
  533. Q1:
  534.      D = ATN(DY/DX)
  535.      RETURN          'Q1
  536. Q2orQ3:
  537.      D = P1 + ATN(DY/DX)
  538.      RETURN          'Q2 OR Q3
  539. deg90:
  540.      D = P1/2
  541.      RETURN          '90 DEG
  542. Q4:
  543.      D = P2 + ATN(DY/DX)
  544.      RETURN          'Q4
  545.  
  546.      '- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  547.      'FNC(T)
  548.      'Routine to initialize the C(J,K) coordinate rotation matrix
  549.      'and other parameters associated with the orbital elements.
  550.  
  551. ncsub:
  552.      IF mmot > .1 THEN  sma = ((G0/(mmot^2))^(1/3))        'Given mmot=Mean Motion
  553.      IF mmot <=.1 THEN mmot = SQR(G0/(sma^3))        'Given sma=Semi-major axis
  554.      tmmot = mmot + 2 * (time#-T0#) * drag#          'Mean Motion at epoch time#
  555.      tsma = ((G0/(tmmot^2))^(1/3))                   'SMA at epoch time#
  556.      E2 = 1 - E0^2
  557.      E1 = SQR(E2)
  558.      Q0 = M0/360 + K0                    'Q0 = initial orbit phase
  559.  
  560.      'Account for nodal effects due to lumpy gravity field due to
  561.      'the flattened, oblate spheroid shape of Earth.
  562.  
  563.      K2 = 9.95 * ((R0 / tsma)^3.5) / (E2^2)
  564.  
  565.      'Update elements to current epoch and evaluate their SIN/COSs
  566.  
  567.      S1 = SIN(incl * P0)
  568.      C1 = COS(incl * P0)                   'incl = inclination
  569.      traan# = raan# - (time# - T0#) * K2 * C1
  570.      S0 = SIN(traan# * P0)
  571.      C0 = COS(traan# * P0)           'traan# = R.A.A.N.(deg) at epoch T
  572.      targp = argp + (time# - T0#) * K2 * (2.5 * (C1^2) - 0.5)
  573.      S2 = SIN(targp * P0)
  574.      C2 = COS(targp * P0)            'targp  = Arg of Perigee
  575.  
  576.      'Set up coordinate rotation matrix for the current orbit
  577.  
  578.      C(1,1) = +(C2 * C0) - (S2 * S0 * C1)
  579.      C(1,2) = -(S2 * C0) - (C2 * S0 * C1)
  580.      C(2,1) = +(C2 * S0) + (S2 * C0 * C1)
  581.      C(2,2) = -(S2 * S0) + (C2 * C0 * C1)
  582.      C(3,1) = +(S2 * S1)
  583.      C(3,2) = +(C2 * S1)
  584.      RETURN
  585.  
  586.      '- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  587.      'DEF FNM(T)
  588.      'Function to evaluate manom = Mean Anomaly in 0-2PI range
  589.      'Revised to incl drag (drag#)
  590. nmsub:
  591.      Q  = Q0 + mmot * (time# - T0#) + drag# * ((time# - T0#)^2)
  592.      orbnr  = fix(Q + .000001)
  593.      M9 = fix((Q - orbnr + .000001) * 256)
  594.      manom  = (Q - orbnr) * P2
  595.      RETURN
  596.      '- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  597.  
  598.      'DEF FNK(M)
  599.      'Routine to solve Kepler's equation given manom
  600.  
  601.      'Returns Satellite's GeoCentric coordinates
  602.  
  603. nkmsub:                            'Initial trial value
  604.      E = manom + E0 * SIN(manom) + 0.5 * (E0^2) * SIN(2 * manom)
  605.  
  606. nkmsub1:    'Iteration loop solves Kepler's transcendental equation
  607.      S3 = SIN(E)
  608.      C3 = COS(E)
  609.      R3 = 1 - E0 * C3
  610.      M1 = E - E0 * S3
  611.      M5 = M1 - manom
  612.      IF ABS(M5) < 9.999999E-06 THEN nkmsub2 ELSE E = E - M5/R3
  613.      GOTO nkmsub1
  614.  
  615. nkmsub2:
  616.      X0 = tsma * (C3 - E0)
  617.      Y0 = tsma * E1 * S3
  618.      R  = tsma * R3           'In the plane of the orbit
  619.      X1 = X0 * C(1,1) + Y0 * C(1,2)
  620.      Y1 = X0 * C(2,1) + Y0 * C(2,2)
  621.      Z1 = X0 * C(3,1) + Y0 * C(3,2)
  622.  
  623.      'Rotate through current GHA of Aries, convert to Geocentric coordinates
  624.  
  625.      G7 = time# * G1 + G2
  626.      G7 = (G7 - fix(G7)) * P2
  627.      S7 = -SIN(G7)
  628.      C7 = COS(G7)
  629.      X  = +(X1 * C7) - (Y1 * S7)
  630.      Y  = +(X1 * S7) + (Y1 * C7)
  631.      Z  = Z1
  632.      RETURN
  633.      '- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  634.  
  635.      'DEF FNO(D)
  636.      'Routine to evaluate observer's geocentric coordinates, where
  637.      'X-Axis=Greenwich, Y-Axis goes thru India, Z-Axis=North Pole
  638.  
  639. initobsv:
  640.      L8 = L9 * P0            'obsrvralt=Height in meters
  641.      S9 = SIN(L8)            'wlong= West Longitude
  642.      C9 = COS(L8)                      'Initial geocentric coordinates
  643.      S8 = SIN(-wlong * P0)
  644.      C8 = COS(wlong * P0)                     
  645.      R9 = R0 * (1 - (F/2) + (F/2) * COS(2 * L8)) + obsrvralt/1000
  646.  
  647.      'Now to make L8 the geocentric latitude
  648.  
  649.      L8 = ATN((1 - F)^2 * S9/C9)
  650.      Z9 = R9 * SIN(L8)
  651.      X9 = R9 * COS(L8) * C8
  652.      Y9 = R9 * COS(L8) * S8
  653.      RETURN
  654.      '- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  655.  
  656.      'DEF FNX(T)
  657.      'Routine to extract all the parameters you might never need
  658.  
  659.      'First get vector from observer to satellite
  660. nxtsub:
  661.      X5 = (X - X9)
  662.      Y5 = (Y - Y9)
  663.      Z5 = (Z - Z9)
  664.      range = SQR(X5^2 + Y5^2 + Z5^2)
  665.  
  666.      'Finite difference the range (range) to get the velocity (velocity)
  667.  
  668.      IF oldtime <> time# THEN
  669.                        velocity = ((oldrange - range)/(oldtime - time#))/86400!
  670.                   ELSE
  671.                        velocity = -8.999999E+09
  672.      END IF
  673.  
  674.      oldrange = range        'Save current time and range for next time
  675.      oldtime = time#
  676.  
  677.      'Now rotate into observer's local coordinates
  678.      'where X8 = North, Y8 = East, Z8 = Up (Left handed system)
  679.  
  680.      Z8 = +(X5 * C8 * C9) + (Y5 * S8 * C9) + (Z5 * S9)
  681.      X8 = -(X5 * C8 * S9) - (Y5 * S8 * S9) + (Z5 * C9)
  682.      Y8 = +(Y5 * C8) - (X5 * S8)
  683.      S5 = Z8/range
  684.      C5 = SQR(1 - S5 * S5)
  685.      satelev = (ATN(S5/C5))/P0            'satelev=Elevation
  686.      DX = X8
  687.      DY = Y8
  688.      GOSUB nasub
  689.  
  690.      az = D/P0              'FNA resolves quadrant az = Azimuth
  691.      DX = X
  692.      DY = Y
  693.      GOSUB nasub
  694.  
  695.      W5 = 360 - D/P0              'W5=Subsatellite point W.Long
  696.      B5 = Z/R
  697.      L5 = ATN(B5/(SQR(1 - B5^2)))/P0     'L5=SSP LAT
  698.      RETURN                     '---NOTE R-R0=Sat. altitude above mean spheroid
  699.  
  700.      '- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  701.  
  702.      'FNT and FNI moved to front of program so they'll be initialized.
  703.  
  704.     '* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
  705.  
  706. $include "DATESUBS.BAS"
  707.  
  708.  
  709. '= = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
  710.  
  711. '            LABEL CROSS REFERENCE
  712.  
  713. '    OLD     NEW           OLD       NEW
  714. '-----------------------------------------------------------
  715. '    A0    tsma        O    traan
  716. '    E8    minelev        O0    raan
  717. '    E9    satelev        R5    range
  718. '    K    orbnr        R6    oldrange
  719. '    K8    (dayNRsav)    R8    velocity
  720. '    K9    (orbNRsav)    T    time#
  721. '    M    manom        T6    oldtime
  722. '    N    tmmot        W    targp
  723. '    N0    mmot        W0    argp
  724. '