This program predicts pass times for NOAA polar-orbiting weather satellites (NOAA 6, 8, and 9) from a user-specified ground station location. It reads orbital parameters—period, equatorial crossing time, and longitude—from DATA statements dated December 1, 1985, then iterates orbit by orbit to propagate the satellite’s ground track forward in time. A subroutine at line 1300 computes azimuth and elevation for antenna aiming using spherical trigonometry, accepting the observer’s latitude and longitude as inputs. The program also supports an optional printer for multi-day pass listings and can filter output to daylight-only passes based on whether the satellite is a morning or afternoon orbiter. Orbital longitude is advanced each pass using the ratio of period to minutes-per-day, providing a simple but effective ground-track propagation model without full Keplerian elements.
Program Analysis
Program Structure
The program is organized into a main prediction loop, several display/formatting subroutines, a date-arithmetic package, and a spherical-geometry tracking calculator. Execution begins at line 120 (via the GO TO 120 at line 40), initializes constants, then proceeds through satellite selection, parameter input, and the orbit propagation loop. The overall flow is:
- Lines 120–170: Global initialization (constants, observer position, geometry limits)
- Lines 180–230: Satellite selection menu (NOAA 6, 8, or 9)
- Lines 250–310: Read and display orbital parameters from
DATA - Lines 320–510: Date input, daylight-filter selection, loop setup
- Lines 50–110: Main orbit propagation loop (advancing time and longitude per orbit)
- Lines 520–610: Pass output dispatcher
- Lines 620–740: End-of-day menu (next day, another day, another satellite, end)
- Lines 770–810: Date/header print subroutine
- Lines 820–890: Date increment subroutine (with leap year and month-length logic)
- Lines 900–1040: Day-of-year subroutine (converts month/day to day number)
- Lines 1050–1170: Pass display subroutine with proximity markers
- Lines 1180–1290: Overhead prediction geometry (local pass time and longitude)
- Lines 1300–1430: Azimuth and elevation calculator (spherical trig)
- Lines 1440–1650: Per-minute tracking data generator
- Lines 1670–1690:
DATAfor the three satellites
Orbital Propagation Model
The propagation model is deliberately simple. Each orbit, the equatorial crossing longitude is advanced by inclong, computed at line 490 as (p/md)*360—that is, the fraction of a full day represented by the orbital period, multiplied by 360°. The variable md is fixed at 1440 (minutes per day). This is a first-order approximation that ignores nodal precession and atmospheric drag but is adequate for short-term prediction with fresh reference elements.
The date/time state is stored compactly: d encodes year and day-of-year as YYYY*10000 + DOY, and the custom function FN r(t) at line 30 extracts the day-of-year portion via t - INT(t/1000)*1000. Note that this actually takes t mod 1000, not mod 10000; since day-of-year never exceeds 366, this works correctly in practice for years where the last three digits of YYYY*10000 + DOY are unambiguous.
Spherical Trigonometry (Lines 1300–1430)
The azimuth/elevation subroutine uses standard great-circle formulas. Given the observer’s geodetic latitude/longitude (LLAT, LLONG) and the satellite’s sub-satellite point (SLAT, SLONG), it computes:
COSB— cosine of the central angle between observer and sub-satellite pointANGB— the central angle itself viaACS COSBCOSA— cosine of the azimuth angle via the spherical law of cosinesAZ— azimuth, with quadrant correction usingflag(set when longitude difference ≥ 0)TANE— tangent of elevation using satellite altitudehs=800km and Earth radiusre=6371kmEL— elevation angle viaATN TANE
The subroutine note at line 1300 states it works in the northern hemisphere only, consistent with the azimuth quadrant logic which assumes north-up orientation.
Satellite Ground-Track Computation (Lines 1440–1650)
The per-minute tracking subroutine iterates from a start angle a1 (2 for ascending, 31 for descending passes) to the full orbital period p. At each step f:
- Sub-satellite latitude:
slat = ASN(SIN(ir) * SIN(2*PI*f/p))— sinusoidal latitude variation for a circular orbit at inclinationir - Sub-satellite longitude: derived from
ACS(COS(2*PI*f/p) / COS(slat))with an Earth-rotation correction off/4degrees per minute (approximately correct for a ~101-minute orbit)
The variable minct counts minutes with elevation above 4° as “usable data” time, printed at pass end. The xflg variable at line 1650 is set after the loop completes but the GO TO 1560 it triggers would call the azimuth/elevation subroutine with the last computed position—this appears to be dead code or a vestigial attempt at a wrap-around, since there is no corresponding loop termination guard.
Date Arithmetic
The date increment subroutine (lines 820–890) handles variable month lengths and leap years explicitly with conditional checks rather than a lookup table. The leap year flag ly is set at line 350 using the standard k - INT(k/4)*4 = 0 test (where k is the two-digit year), which correctly identifies years divisible by 4 but does not account for century years—acceptable for the near-term predictions this program targets. The day-of-year subroutine (lines 900–1040) uses a chain of IF l=N THEN LET r=... statements, one per month, with a leap-year correction applied afterward at line 1020.
Proximity Markers and Daylight Filtering
The pass display subroutine (lines 1050–1080) computes ld, the longitude difference between the sub-satellite point and the observer. Passes with ld ≤ 30 are printed; passes with ld ≤ 11 receive a ** marker, and those with ld ≤ 5 receive ***, indicating increasingly close overhead passes.
Daylight filtering is applied in the LPRINT path (lines 1090–1100): NOAA 9 (afternoon satellite, dflg=1) suppresses passes before 11:00 local; NOAA 6 and 8 (morning satellites, dflg=2) suppress passes after 12:00. This reflects the sun-synchronous nature of the NOAA polar orbiters.
Printer Detection
Line 140 uses IN 251 to detect whether a printer interface is connected—a hardware port query returning a value less than 63 is taken as confirmation. When a printer is present (prntr=1), the user is prompted for a number of days (NDAY) to print in batch, and LPRINT statements throughout the program produce a parallel paper copy of all pass data.
Key Variables
| Variable | Purpose |
|---|---|
a | Current orbit time in minutes since midnight |
c | Current equatorial crossing longitude (degrees west) |
d | Reference epoch encoded as YYYY×10000 + DOY |
t | Target date encoded the same way |
p | Orbital period in minutes |
inclong | Longitude advance per orbit (degrees) |
ld | Longitude distance from observer to sub-satellite point |
hd | Hours offset for local time (4 in summer, 5 in winter — Eastern US) |
dflg | Daylight filter flag (0=none, 1=afternoon, 2=morning) |
satinc | Orbital inclination in degrees (from DATA) |
ir | Orbital inclination in radians |
re, hs | Earth radius (km) and satellite altitude (km) |
Notable Anomalies
- Line 30:
FN r(t) = t - INT(t/1000)*1000computest mod 1000, notmod 10000as would be needed for a strict YYYY×10000+DOY encoding. It works only because DOY values (1–366) never exceed 999, so the last three digits are always the day-of-year. - Lines 1640–1650: After the
FOR/NEXTloop at line 1530 completes without an earlyRETURN,xflg=1is set and execution falls toGO TO 1560, which jumps back into the middle of the loop body to call the azimuth subroutine. This creates an infinite loop if no pass is ever detected (elevation never positive), since there is no exit condition for thexflg=1case. - Line 880: The hardcoded Eastern US UTC offset (4 for EDT, 5 for EST) is embedded as
hdand not parameterized, making the program implicitly location-specific to the US Eastern time zone despite accepting arbitrary latitude/longitude at line 150. - Line 1190: The comment “northern hemisphere only” is enforced implicitly by the azimuth quadrant logic but not by any explicit latitude check; negative latitudes would produce incorrect azimuth readings without error.
Content
Image Gallery
Source Code
10 REM Orbiter Charles A. Judge M.D. based on work of Grant Zehr and Ralph Taggart. Math routines based on material from The Satellite Experimenter's Handbook ARRL, Martin Davidoff
20 REM This program will compute arrival times for the NOAA polar orbiting weather satellites. You will need to insert your appropriate latitude and longitude in line 150. The antenna aiming equations work in the nortern hemisphere only. Enter all longitudes as positive west from Greenwich E.G. 170 E = 190 W
30 DEF FN r(t)=t-((INT (t/1000))*1000)
40 GO TO 120
50 LET a=a+p: LET n=n+1: IF a>md THEN LET a=a-md: LET d=d+1: IF FN r(d)>365 THEN LET d=k*10000+1
60 IF w<>1 THEN PRINT AT 15,5;d
70 LET c=c+inclong: IF c>360 THEN LET c=c-360
80 IF d<t THEN GO TO 50
90 IF w<>1 THEN GO SUB 770
100 LET w=1
110 GO TO 520
120 LET w=0: LET md=1440: LET d$="": LET px=1: LET x$="":: LET re=6371: LET hs=800: LET dr=2*PI/360: LET rd=1/dr
130 LET k=85
140 LET prntr=0: IF IN 251<63 THEN LET prntr=1
150 LET lat=38.5: LET long=76.6: REM change for your position!
160 LET lae=long-35: LET law=long+42: LET lde=long+140: LET ldw=long+225
170 CLS : PRINT "Polar Orbit Crossing Program"''"Predicts based on NOAA data as"'"of 12/1/85"
180 PRINT ;''"Which Satellite?"''
190 PRINT TAB 5;"0 for NOAA 6";TAB 5;"1 for NOAA 9";TAB 5;"2 for NOAA 8"
200 LET n$=INKEY$: IF n$="" THEN GO TO 200
210 IF n$>"2" OR n$<"0" THEN GO TO 200
220 LET n3=VAL n$
230 RESTORE 1670+10*n3
250 LET dflg=0
260 IF n3=0 OR n3=2 THEN LET dflg=2
270 IF n3=1 THEN LET dflg=1
280 READ a$,p,l,j,k,n,h,m,c,satinc
290 CLS : PRINT AT 2,0;"Satellite: ";a$;AT 5,0;"Orbital Parameters as of ";l;"/";j;"/";k;AT 8,0;"Period";TAB 13;p;AT 10,0;"Eq. Crossing";AT 11,0;"Time";TAB 13;h;":";INT m;AT 13,0;"Longitude";TAB 13;INT c
300 LPRINT ''"Satellite: ";a$'"Orbital Parameters as of ";l;"/";j;"/";k'"Period";TAB 13;p'"Eq. Crossing"'"Time";TAB 13;h;":";m'"Longitude";TAB 13;c
310 CLS
320 IF dflg=1 THEN PRINT AT 7,0;"Afternoon Satellite - ";AT 10,0;"Do you want daylight passes only?": INPUT q$: IF q$<>"y" AND q$<>"Y" THEN LET dflg=0
330 IF dflg=2 THEN PRINT AT 7,0;"Morning Satellite - ";AT 10,0;"Do you want daylight passes only?": INPUT q$: IF q$<>"y" AND q$<>"Y" THEN LET dflg=0
340 CLS
350 LET ly=0: IF k-(INT (k/4)*4)=0 THEN LET ly=1
360 GO SUB 900: LET d=r+(k*10000)
370 LET a=(h*60)+m
380 PRINT AT 5,0;"Enter Month, Day, and Year ";AT 7,10;"For Predict"
390 INPUT l: PRINT AT 10,15;l;"/";: INPUT j: PRINT j;"/";: INPUT k: PRINT k
400 LET ND=0
410 IF prntr=1 THEN PRINT '"ENTER NUMBER OF DAYS ";: INPUT NDAY: PRINT NDAY
420 PRINT '"Will you want tracking data? (Y/N)": INPUT z$
430 GO SUB 900: LET t=r+(k*10000)
440 IF t>=d THEN GO TO 470
450 CLS : PRINT ''"This program does not compute data prior to reference crossing"
460 GO TO 380
470 CLS : PRINT AT 10,5;"Now Calculating Orbital";AT 12,10;"Crossings"
480 LET p$="": LET ps=0
490 LET inclong=(p/md)*360
500 PRINT AT 15,15;t
510 GO TO 50
520 GO SUB 540
530 GO TO 510
540 LET h=INT (a/60): LET m=INT (a-(h*60)+.5)
550 IF m=60 THEN LET m=0: LET h=h+1
560 LET b=INT (c+.5)
570 GO SUB 1180
580 IF px=0 THEN LET p$=""
590 GO SUB 1050
600 RETURN
610 GO TO 510
620 LET px=0: GO SUB 540: LET px=1: IF prntr=1 THEN LET ND=ND+1: IF ND>=NDAY THEN GO TO 660
630 IF prntr=0 THEN GO TO 670
640 LET t=t+1: GO SUB 820: GO SUB 770
650 GO TO 510
660 IF ps=1 THEN PRINT "@ Indicates previous day": LET ps=0
670 LET w=0: PRINT ''"(N) Next day"'"(O) Another day"'"(S) Another Spacecraft"'"(E) End": PRINT '"Enter Choice ";: FLASH 1: PRINT "?"
680 LET Q$=INKEY$: IF Q$="" THEN GO TO 680
690 FLASH 0
700 IF q$="s" OR q$="S" THEN CLS : GO TO 180
710 IF q$="n" OR q$="N" THEN LET t=t+1: GO SUB 820: GO TO 510
720 IF q$="o" OR q$="O" THEN CLS : GO TO 380
730 IF q$="e" OR q$="E" THEN PRINT "GOODBYE": STOP
740 GO TO 680
750 CLS : GO SUB 900: LET t=r+(k*10000)
760 IF t>d THEN GO SUB 770: GO TO 510
770 CLS
780 PRINT "Date: ";l;"/";j;"/";k;" ";a$: PRINT
790 LPRINT : LPRINT l;"/";j;
800 PRINT "Orbit Local Time Longitude": PRINT
810 RETURN
820 LET j=j+1
830 IF ((l=4) OR (l=6) OR (l=9) OR (l=11)) AND (j>30) THEN LET l=l+1: LET j=1: GO TO 870
840 IF (l=2) AND (j>28) AND ly=0 THEN LET l=l+1: LET j=1: GO TO 870
850 IF (l=2) AND (j>29) AND ly=1 THEN LET l=l+1: LET j=1: GO TO 870
860 IF j>31 THEN LET l=l+1: LET j=1
870 IF l>12 THEN LET l=1: LET k=k+1
880 LET hd=5: IF l>4 AND l<11 THEN LET hd=4
890 RETURN
900 IF l=1 THEN LET r=0
910 IF l=2 THEN LET r=31
920 IF l=3 THEN LET r=59
930 IF l=4 THEN LET r=90
940 IF l=5 THEN LET r=120
950 IF l=6 THEN LET r=151
960 IF l=7 THEN LET r=181
970 IF l=8 THEN LET r=212
980 IF l=9 THEN LET r=243
990 IF l=10 THEN LET r=273
1000 IF l=11 THEN LET r=304
1010 IF l=12 THEN LET r=334
1020 IF l>2 AND ly=1 THEN LET r=r+1
1030 LET hd=5: IF l>4 AND l<11 THEN LET hd=4
1040 LET r=r+j: RETURN
1050 IF ld>30 THEN RETURN
1060 IF ld<=11 THEN LET x$="**"
1070 IF ld<=5 THEN LET x$="***"
1080 PRINT n;TAB 7;p$;TAB 8;sh;":";sm;TAB 19;sl;TAB 23;d$;TAB 28;x$: IF p$="@" THEN LET ps=1
1090 IF dflg=1 AND (sh<11) THEN GO TO 1120
1100 IF dflg=2 AND sh>12 THEN GO TO 1120
1110 LPRINT TAB 7;p$;TAB 8;sh;":";sm;TAB 16;sl;TAB 22;d$;TAB 27;x$: IF p$="@" THEN LET ps=1
1120 LET x$=""
1130 LET b$=INKEY$: IF b$="" THEN GO TO 1160
1140 IF b$="a" THEN LET z$="y": GO TO 1160
1150 LET z$="n"
1160 IF z$="Y" OR z$="y" THEN INPUT "Do you want tracking data",q$: IF q$="y" OR q$="Y" THEN GO SUB 1440
1170 RETURN
1180 IF NOT ((c>lae AND c<law) OR (c>lde AND c<ldw)) THEN LET ld=60: RETURN
1190 IF c>=0 AND c<180 THEN LET td=(lat/360)
1200 IF c>=180 THEN LET td=(180-lat)/360
1210 LET sm=INT (m+(p*td)): IF td<.25 THEN LET sl=INT (c+(inclong/p*td)): LET d$="asc"
1220 IF td>.25 THEN LET sl=INT (c+183.1): LET d$="desc"
1230 IF sl>360 THEN LET sl=sl-360
1240 IF sl<=long THEN LET ld=long-sl
1250 IF sl>long THEN LET ld=sl-long
1260 LET p$="": LET sh=h: LET sm=INT sm
1270 IF sm>=60 THEN LET sm=sm-60: LET sh=h+1
1280 LET sh=sh-hd: IF sh<0 THEN LET sh=sh+24: LET p$="@"
1290 RETURN
1300 REM AZIMUTH AND ELEVATION CALCULATOR
1310 LET LLONG=-LONG*Dr: LET LLAT=LAT*DR
1320 LET Longd=Llong-SLONG: IF LONGD<-PI THEN LET LONGD=LONGD+(2*PI)
1330 IF LONGD>PI THEN LET LONGD=LONGD-(PI*2)
1340 LET flag=0: IF lONGD>=0 THEN LET flag=1
1350 LET COSB=(SIN (LLAT)*SIN (SLAT))+(COS (LLAT)*COS (SLAT)*COS (LONGD))
1360 LET ANGB=ACS COSB
1370 LET COSA=(SIN (SLAT)-(SIN (LLAT)*COSB))/(COS LLAT*SIN ANGB)
1380 LET AZ=ACS COSA
1390 IF flag=1 THEN LET az=(2*PI)-az
1400 LET TANE=(((RE+HS)*COSB)-RE)/((RE+HS)*SIN ANGB)
1410 LET EL=ATN TANE
1420 LET AZ=AZ*RD: LET EL=EL*RD
1430 RETURN
1440 LET ir=satinc*dr
1450 LET pflg=0: LET xflg=0
1460 CLS : PRINT "Date: ";l;"/";j;"/";k;" ";a$: PRINT
1470 LPRINT
1480 PRINT AT 2,0;"TIME";TAB 6;"LAT";TAB 10;"LONG";TAB 17;"AZ";TAB 22;"EL": PRINT
1490 LPRINT : LPRINT "TIME";TAB 6;"LAT";TAB 10;"LONG";TAB 17;"AZ";TAB 22;"EL"
1500 LET minct=0
1510 IF c<180 THEN LET a1=2
1520 IF c>180 THEN LET a1=31
1530 FOR f=a1 TO p
1540 LET slat=ASN (SIN (ir)*SIN (2*PI*f/p)): LET slatd=slat*rd
1550 LET slong=(-c*dr)-(ACS (COS (2*PI*f/p)/(COS slat)))-(f/4*dr): LET slongd=slong*rd
1560 GO SUB 1300
1570 LET lt=a+f: LET h=INT (lt/60): LET m=INT (lt-(h*60)+.5): IF m=60 THEN LET m=0: LET h=h+1
1580 LET h=h-hd: IF h<0 THEN LET h=h+24
1590 IF slongd>360 THEN LET slongd=slongd-360
1600 IF slongd<-360 THEN LET slongd=slongd+360
1610 IF el>0 THEN PRINT h;":";m;TAB 7;INT slatd;TAB 12;INT -slongd;TAB 17;INT az;TAB 22;INT el: LET pflg=1: IF el>4 THEN LET minct=minct+1
1620 IF el>0 THEN LPRINT h;":";m;TAB 7;INT slatd;TAB 12;INT -slongd;TAB 17;INT az;TAB 22;INT el
1630 IF pflg=1 AND el<0 THEN LPRINT : LPRINT minct;" Minutes usable data": LPRINT : PRINT 'minct;" Minutes usable data": PRINT : RETURN
1640 NEXT f
1650 LET xflg=1: GO TO 1560
1660 REM 9000
1670 DATA "NOAA 6 137.5",101.1277,12,1,85,33421,00,31.62,86.22,99
1680 DATA "NOAA 9 137.62",102.0851,12,1,85,4988,0,46.32,152.93,99.0506
1690 DATA "NOAA 8 137.5",101.2979,12,01,85,13909,0,57.76,80.42,99
1700 CLEAR : SAVE "NOAA" LINE 1Note: Type-in program listings on this website use ZMAKEBAS notation for graphics characters.