*** EOOH *** From: James Miller Date: Tue, 18 May 93 22:38:33 BST To: amsat-bb@amsat.org Subject: Tracking with a Calculator Cc: james.miller@mrc-apu.cam.ac.uk A recent contributor asked for satellite tracking algorithms, in order to transfer them to a calculator. Perhaps it is timely to confirm that it is indeed not necessary to use a 256 Kbyte program running on a 486-50DX PC to compute a satellite's position, albeit 300 times a second. It can be done swiftly (and probably more accurately) on a modest HP15C in two dozen lines of code. Below is a program originally written by Karl Meinzer DJ4ZC a few years ago. It was in BASIC for the Sharp PC1246, and is extremely popular. It computes range, azimuth, elevation and MA for Oscar-13. It is virtually as Karl presented it, but I have removed the GOSUBs for clarity (but not the GOTO). The encoding is for AO-13, though obviously it can be adapted for any other satellite, and I/O bells and whistles added. If you are at all familiar with the entrails of computing a satellite's position, you will know that apart from the cosmetics, most of it concerns coordinate transformations. Good scientific pocket calculators feature Rectangular-to-Polar and Polar-to-Rectangular transformations, and exploiting this contributes to Karl's tight coding (also in the AO-13 Flight Software). The two equivalent calls in the BASIC below are PROCrtp and PROCptr. Calculators also have H.M.S. time parsers, and I have used DEF FNdms(H) for this purpose. Now follows the code, and a sample run. After that are some comments on the numeric constants embedded in the program. --------------------------------------------------------------------- 1REM Based on Karl Meinzer's (DJ4ZC) Algorithms for Sharp PC1246 2REM ----------------------------------------------------------- 3REM 4REM Minor changes added by James Miller G3RUH 1989-1993 5REM 10Y=93:G=126.522217:I=57.76:O=318.04:E=0.7243: REM AO-13 Smoothed 11W=314.56:M=5.35:N=2.097208 13G=INT(365.25*(Y-1))-28125 + G - M/N/360:N=N*2*PI 14A=4.0424:Q=-0.1846:V=0.0774 16: 100B=RAD 52.21:L=RAD 0.06:PRINT "G3RUH" 320INPUT "DAY=";D,"MTH=";M,"YR=";Y 330M=M+1: IF M<4 THEN Y=Y-1:M=M+12 340D=D+ INT(Y*365.25)+INT(M*30.6)-28553 350INPUT "UTC(HH.MM)= ";U: T=INT U: T= (T+(U-T)/.6)/24 + D 360U=24*(T-INT T) 370PRINT "DY="; INT T;" @ "; FNdms(U); 380: 410D=T-G: K=RAD(D*Q+O-100.29-T*360.985647-L) 420M=D*N: R=E: P=M 430REPEAT PROCptr:H=(M-P+Y)/(1-X):P=P+H 440UNTIL ABS H<0.001 450R=1:PROCptr:Y=Y*SQR(1-E*E):X=X-E 460PROCrtp:R=A*R:P=P+RAD(W+D*V) 470PROCptr:H=X:R=Y:P=RAD(I) 480PROCptr:S=Y:Y=X:X=H 490PROCrtp:P=P+K:K=P+L:U=R 500PROCptr:Z=S-SIN B:X=X-COS B 510H=Y:Y=Z: PROCrtp:P=P-B+PI/2 520PROCptr:J=Y:Y=H:PROCrtp 530C=PI-P:X=R:Y=J:PROCrtp 540M=M/(2*PI):M=256*(M-INT M) 550H=T:PRINT " RAEM = "; INT(R*6378);" ";INT(DEG C+0.5);" ";INT(DEG P+.5);" ";INT(M+.5) 580K=K/2/PI:K=(K-INT K)*2*PI 590X=U:Y=S: PROCrtp 600PRINT"LA= ";INT(DEG P+.5);" LO=";INT(DEG K+.5) 610PRINT 630T=T+1/96: GOTO 360 640: 669DEF PROCptr 670Y=R*SIN P: X=R*COS P: ENDPROC 678: 679DEF PROCrtp 680R=SQR(X*X+Y*Y): IF X=0 THEN P=PI/2*SGN Y:ENDPROC 690P=ATN(Y/X): IF X>0 THEN ENDPROC 700IF Y>0 THEN P=P+PI: ENDPROC 710P=P-PI: ENDPROC 800 10000DEF FNdms(H) 10005H=H+.5/3600 10010h=INT(H):m=INT(60*(H-h)) 10020H$=STR$(h): IF LEN(H$)<2 H$="0"+H$ 10030M$=STR$(m): IF LEN(M$)<2 M$="0"+M$ 10050=H$+M$ --------------------------------------------------------------------- RUN PC1246 ( for 1993 May 18 @ 2100 utc ) G3RUH DAY=?18 MTH=?5 YR=?93 UTC(HH.MM)= ?21.00 DY=5616 @ 2100 RAEM = 18127 225 -9 236 << 5616 is Amsat Day No. LA= -15 LO=314 << Sub sat point DY=5616 @ 2115 RAEM = 15784 216 -20 241 LA= -26 LO=319 DY=5616 @ 2130 RAEM = 13589 200 -36 247 LA= -43 LO=332 etc Numerical Constants, comments etc: ---------------------------------- Line 10,11,14 keplerian elements: these are Oscar-13 Smoothed. Valid for the rest of 1993 I should think: Y=93 Epoch year years G=126.522217 Epoch Day day I=57.76 Inclination deg O=318.04 R.A.A.N. deg E=0.7243 Eccentricity - W=314.56 Argument of Perigee deg M=5.35 Mean Anomaly deg N=2.097208 Mean Motion rev/day Line 14, auxiliary data. Hardly changes, if ever. A= 4.0424 Semi-major axis Earth radii Q=-0.1846 Rate of change of RAAN deg/day V= 0.0774 Rate of change of Arg Peri deg/day Line 100, Observer's location deg. +N, +E degrees B=RAD 52.21:L=RAD 0.06:PRINT "G3RUH" Other interesting constants: Line 410 100.29 GHAA for 1978 Jan 01 @ 0000 utc, the "siderial whatsit" 360.985647 Earth's rotation rate, deg/day Line 550 6378 Earth's radius, km RAEM Range, Azimuth, Elevation, MA ( km, deg, deg, /256) Line 630 T=T+1/96: GOTO 360 Time increment, days. The constant PI = 3.141592654 RAD(X) is degrees to radians DEG(X) is radians to degrees INT(X) is the integer part, the largest integer smaller than X. HP calculators get this wrong for negative arguments. Plunder, enjoy, improve! ------------------------------------------------------------------------------ If you are interested my "Pamphlet" mentioned in recent amsat-bb correspondence, I will air-mail you a copy in return for $2 US or 5 IRC (for postage), and of course your address. The 12 page booklet documents all the calculations needed to track a satellite, plus solar illumination, squint, visibility, eclipses, footprint, direct doppler and much more. The routines form the basis of a good many programs in common use. Not many people know that. +-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-+ | 73 de James G3RUH | | james.miller@mrc-apu.cam.ac.uk | | StarDate: 1993 May 18 [Tue] 2049 utc | +-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-+ In-Reply-To: <9304261923.AA03699@silver.lcs.mit.edu> Date: Tue, 27 Apr 93 23:05:50 -0400 To: "David Moisan" From: "Seth M. Dworken" Organization: KC3IL, Fairless Hills PA, USA Reply-To: p00123@psilink.com Subject: Re: Shouldn't ask this, but I need W3IWI's Orbit code... X-Mailer: PSILink-DOS (3.4) Here you go -- this is the North Star Basic version. Good luck! -Seth ===== 10 REM --------------------------------------------------------- 20 REM ******* ******* 30 REM ******* AMSAT ORBITAL PREDICTION PROGRAM de W3IWI ******* 40 REM ******* ******* 50 REM ******* F I L E = O R B I T ******* 60 REM ******* ******* 70 REM ******* d e W 3 I W I 1 0 J a n . 8 2 ******* 80 REM ******* ******* 90 REM --------------------------------------------------------- 100REM 110 PRINT "ERROR MEANS YOU FORGOT TO APPEND ORBITAL ELEMENTS !!",CHR$(7) 120 RESTORE 20000 \ PRINT" THAT'S GOOD!",CHR$(7) 130 DIM T$(20),S$(40),I$(40),C(3,2) \ P=0 \ REM--P=Page # counter 140 DIM C1$(50),X(4,9) \ REM PLACE-HOLDERS FOR STATION NAME/COORD 150 DIM S(10),T1$(20) \ REM--PLACE-HOLDERS FOR MULTI-SATELLITE LOOP 160REM 170 REM ------------------------------ 180 GOSUB 2660 \ REM*** (RE)Initialization for multiple passes *** 190REM 200 REM ------------------------------ 210 REM ***** SIGN-ON HEADER AND DISCLAIMER ***** 220 PRINT" AMSAT ORBITAL PREDICTION PROGRAM de W3IWI -- May, 1980" 230 PRINT" This Version Prepared for Distribution 10 Jan. 1982" 240 PRINT 250 PRINT" COPYRIGHT 1980 by Dr. Thomas A. Clark, W3IWI" 260 PRINT" 6388 Guilford Road" 270 PRINT" Clarksville, MD 21029" 280 PRINT" Permission granted for non-commercial use providing" 290 PRINT" credit is given to the author, AMSAT and ORBIT Magazine." 300 PRINT CHR$(7) 310REM 320 REM ------------------------------ 330 REM--***** REVISIONS FROM THE VERSION IN ORBIT #6 ***** 340 REM 350 REM (1) ORBITAL ELEMENTS MOVED FROM THE LINES RANGING FROM 360 REM LINES 1100 THRU 9999 TO LINES BEGINNING 20000. 370 REM THIS ALLOWS THE ORBITAL ELEMENTS TO BE SEPARATE 380 REM FILES WHICH CAN BE MERGED WITH "APPEND". W3IWI 15 July 81 390 REM 400 REM (2) A PROCEEDURE HAS BEEN ADDED TO ALLOW THE PROGRAM 410 REM TO STEP MORE RAPIDLY WHILE THE SATELLITE IS BELOW 420 REM THE ELEVATION LIMIT (E8). THE PROCEEDURE INVOLVES 430 REM FORCIBLY MANIPULATING THE TIME VARIABLE (T) 440 REM WHILE THE SATELLITE IS BELOW THE LIMIT BY THE AD HOC 450 REM ALGORITHM IN LINE 2120. ONCE THE SATELLITE RISES, THE 460 REM TIME IS RE-SYNCHRONIZED IN LINE 2100. THE 3-POSITION 470 REM SWITCH (D0) IS USED TO SELECT THE ACTION NEXT TIME 480 REM THRU THE LOOP. 490 REM W3IWI 15 July 81 500 REM 510 REM (3) INCLUSION OF A PERIOD DERIVATIVE TERM IN MEAN MOTION 520 REM TO BETTER MOm[G+ LOW-ALTITUDE SATELLITES LIKE UoSAT. 530 REM THIS INVOLVES CHANGES IN THE FNx "SUBROUTINES" IN LINES 540 REM 10140-10156, 10190, 10345-10350, 10450, PLUS THE NEW 550 REM VARIABLES N & A (MEAN MOTION & SMA AT T) & N1 (DRAG TERM). 560 REM NOTE THAT YOU MUST ADD THE N1 DRAG TERM (ORBITS/DAY^2) 570 REM TO THE ORBITAL ELEMENTS DATA STATEMENTS. THE 1st PAGE 580 REM SUMMARY OF THE ELEMENTS HAS BEEN MADE MORE INFORMATIVE & 590 REM THE S.S.P. LONG & LAT NOW HAVE ONE MORE DIGIT THAN 600 REM IN THE PRINTOUTS IN ORBIT#6 ARTICLE. 610 REM W3IWI 15 NOV 81 620 REM 630 REM (4) IT IS NOW POSSIBLE TO SPECIFY MORE THAN ONE STATION. 640 REM AFTER CALCULATIONS FOR THE "MASTER" STATION HAVE BEEN 650 REM MADE, THE ADDITIONAL STATIONS ARE TESTED TO SEE IF THE 660 REM SATELLITE IS ALSO VISIBLE AT THE OTHER STATIONS. 670 REM 680 REM (5) SINCE RS-3 THRU RS-8 WERE LAUNCHED, WE FELT A NEED 690 REM TO BE ABLE TO SET UP A "BATCH" PROCESSING MODE TO 700 REM HANDLE MULTIPLE SATELLITES. YOU ARE NOW ABLE TO 710 REM SPECIFY A m[G+ OF UP TO 11 SATELLITES FOR PROCESSING. 720 REM 730 REM (6) WITH THE REVISIONS CITED ABOVE, THE m[G+ING WAS GETTING 740 REM RATHER CLUTTERED, SO THE LINE NUMBERS HAVE BEEN CLEANED 750 REM UP BY m[G+MUBERING. THE SUBROUTINES AT 10000 ARE STILL 760 REM AS THEY APPEARED IN ORBIT, WITH SLIGHT REVISIONS MADE 770 REM NECESSARY TO ACCOMMODATE THESE CHANGES. 780 REM W3IWI 10 JAN 82 790 REM 800 REM 73 es gud luck de W3IWI 810 REM ------------------------------ 820REM 830 REM ***** Set starting Day/Time, duration & step size & initialize ***** 840 INPUT1 "Start: Year = ",Y \ Y=Y/100 \ Y2=INT(100* (Y-INT(Y)) + .1) 850 IF Y2/4=INT(Y2/4) THEN F9=1 ELSE F9=0 \ REM--F9=Leap-year flag 860 INPUT1 " Month (1-12) = ",M \ INPUT1 " Day = ",D 870 T1$=FNT$(Y2) + "/" + FNT$(M) + "/" + FNT$(D) + " at " 880 REM--Calculate D8=Day# (Note that month #13 = January of following year) 890 RESTORE 2480 \ FOR I=1 TO M \ READ D9 \ NEXT \ D8=D+D9 900 IF M>2 THEN D8=D8+F9 \ PRINT " Day #",D8 \ PRINT CHR$(7) 910 INPUT1"Start: UTC Hours = ",H 920 INPUT " Min. = ",M \ T7=D8 + H/24 + M/1440 930 T1$=T1$ + FNT$(H) + FNT$(M) + ":00" \ REM--For printing 940 INPUT1"Duration: Hours = ",H1 950 INPUT " Min. = ",M1\ T8=T7 + H1/24 + M1/1440 960 INPUT1"Time Step: Min. = ",M2 \ T9=M2/1440 970 PRINT %10F4," From ",T7," to ",T8 \ PRINT CHR$(7) 980 INPUT "Output Unit # ",P9 \ IF P9<>0 THEN C9$=C8$ 990 INPUT"Print Summary of Elements (Y/N)? ",Q2$ 1000 REM ------------------------------ 1010 REM ***** GET VARIOUS PHYSICAL CONSTANTS ***** 1020 RESTORE 2510 \ READ P1,C,R0,F,G0,G1 1030 P2=2*P1 \ P0=P1/180 \ F=1/F 1040 READ Y1,G2 \ IF Y1=Y2 THEN 1090 ELSE IF Y1>0 THEN 1040 1050 PRINT "Unable to find Year ",Y2," in Sidereal Time Table" \ STOP 1060REM 1070 REM ------------------------------ 1080 REM ***** SELECT OBSERVING STATIONS ***** 1090 RESTORE 15000 \ PRINT C9$, "STATION SELECTION MENU" \ I=1 \ PRINT 1100 PRINT "Station # 0 Will Be The Master Station" \ PRINT 1110 READ C$ \ IF C$(1,3)="END" THEN 1170 ELSE READ W9,L9,H9,E8 1120 C$=C$+" " \ C$=C$(1,10) \ REM--Pad Length 1130 REM--Now pick up up to 4 additional stations 1140 PRINT "STN # ",I," = ",C$,"at ",%5F1,W9," WLong", 1150 PRINT %5F1," ",L9," Lat ",E8," El" 1160 I=I+1 \ GOTO 1110 1170 PRINT 1180 FOR K9=0 TO 4\ REM--FILL IN STATION m[G+ 1190 PRINT CHR$(7),"Select Station # ",K9," <0 to end> ", \ INPUT " ",K 1200 IF K=0 THEN EXIT 1270 \ IF K>I THEN 1190 \ RESTORE 15000 \ K5=K9 1210 FOR J=1 TO K \ READ C$,W9,L9,H9,E8 \ NEXT J 1220 D=FNO(D) \ REM--Calculate each observer's coordinates 1230 X(K5,0)=0 \ X(K5,1)=X9 \ X(K5,2)=Y9 \ X(K5,3)=Z9 \ X(K5,4)=0 1240 X(K5,5)=E8 \ X(K5,6)=S8 \ X(K5,7)=C8 \ X(K5,8)=S9 \ X(K5,9)=C9 1250 C$=C$+" " \ C$=C$(1,10) \ REM--Pad Length 1260 C1$(K5*10+1,K5*10+10)=C$ \ NEXT K9 1270 PRINT "THE FOLLOWING STATIONS HAVE BEEN SELECTED " 1280 PRINT " ",C1$ \ C$=C1$(1,10) \ IF K5<0 THEN 180 1290REM 1300 REM ------------------------------ 1310 REM ****** Select parameters for satellite(s) of interest ****** 1320 PRINT "SATELLITE SELECTION MENU" \ I=0 \ RESTORE 20000 \ PRINT 1330 READ S$ \ IF S$(1,3)="END" THEN 1360 ELSE READ I$ \ I=I+1 1340 PRINT "Entry #",I," for ",S$ \ PRINT " ID = ",I$ 1350 READ D,D,D,D,D,D,D,D,D,D,D,D,D,D,D \ GOTO 1330 \ REM--D's are dummy 1360 IF I>1 THEN 1370 ELSE S(0)=1 \ S(1)=1 \ GOTO 1470 1370 PRINT 1380 PRINT CHR$(7), "Select 1st Satellite # <0 for all> ", \ INPUT " ",J 1390 IF J>0 THEN 1410 ELSE S(0)=I 1400 FOR J=1 TO I \ S(J)=J \ NEXT J \ GOTO 1470 1410 S(0)=1 \ S(1)=J \ PRINT "Select Additional Satellites <0 ends list>" 1420 FOR J=2 TO I \ PRINT CHR$(7), "Satellite # ",J,TAB(37), \ INPUT "",K 1430 IF K<=0 THEN EXIT 1470 ELSE S(0)=J \ S(J)=K \ NEXT J 1440REM 1450 REM ------------------------------ 1460 REM ***** OUTER LOOP THROUGH ALL REQUESTED SATELLITES ***** 1470 FOR J4=1 TO S(0) \ RESTORE 20000 1480 FOR I=1 TO S(J4) \ READ S$,I$,Y3,D3,H3,M3,S3 1490 READ K0,M0,N0,A0,N1,I0,E0,W0,O0,F1 \ NEXT I 1500 IF J4>1 THEN 1540 \ REM--SUPRESS DOPPLER INQUIRY IN "BATCH" MODE 1510 PRINT \ PRINT "Doppler calculated for freq = ",F1," MHz" 1520 REM \ INPUT " Change frequency to (0 for default) ",D 1530 IF D<>0 THEN F1=D 1540 IF Y3=Y2 THEN 1590 ELSE PRINT "ELEMENTS NOT FROM CURRENT YEAR" 1550 GOTO 2360 1560 REM--IF D3 is an Integer then epoch is assumed to be in HH,MM,SS 1570 REM as a time-of-day. 1580 REM IF D3 is not Integer, then D3 is epoch in days + fractions. 1590 IF D3=INT(D3) THEN T0=D3 + H3/24 + M3/1440 + S3/86400 ELSE T0=D3 1600REM 1605 IF Q2$<>"Y" THEN 1800 \ REM--Bypass printout summary 1610 REM ------------------------------ 1620 REM ***** PRINTOUT SUMMARY OF ELEMENTS, ETC ***** 1630 PRINT#P9,C7$ \ GOSUB 2390 \ PRINT#P9 \ PRINT#P9 1640 PRINT#P9,"North-Star BASIC Version of W3IWI Orbit Program" 1650 PRINT#P9,"Revisions as of 10 January 1982 Included by W3IWI" 1660 PRINT#P9 \ PRINT#P9,"Orbital Elements for ",S$ 1670 PRINT#P9," Reference ID= ",I$ \ PRINT#P9 1680 PRINT#P9,"Reference Epoch = ",Y3," +",T0 1690 PRINT#P9,"Starting Epoch = ",Y2," +",T7," = ",T1$ 1700 PRINT#P9,"Elements are",%7F2,T7-T0," Days Old at Starting Epoch." 1710 PRINT#P9 1720 PRINT#P9,"Duration =",%9F1,(T8-T7)*1440," Minutes = ", 1730 PRINT#P9,(T8-T7)*24," Hours" 1740 PRINT#P9,"Step Size =",%8F2,T9*1440," Minutes." 1750 PRINT#P9 1760 PRINT#P9,"Parameter",TAB(20),"Reference Epoch", 1770 PRINT#P9,TAB(40),"Starting Epoch" 1780 PRINT#P9,"---------",TAB(20),"---------------", 1790 PRINT#P9,TAB(40),"--------------" 1800 REM--Initialize to epoch T by calling FNC & FNM-- 1810 T=T7 \ D=FNC(T) \ D=FNM(T) 1815 IF Q2$<>"Y" THEN 1960 1820 PRINT#P9,"Orbit Number ",TAB(20),K0,TAB(40),K 1830 PRINT#P9,"Mean Anomaly ",TAB(20),M0,TAB(40),M/P0 1840 PRINT#P9,"Inclination ",TAB(20),I0,TAB(40)," [ No Change ]" 1850 PRINT#P9,"Eccentricity ",TAB(20),E0,TAB(40)," [ No Change ]" 1860 PRINT#P9,"Mean Motion ",TAB(20),N0,TAB(40),N 1870 PRINT#P9,"Drag Correction",TAB(20),N1," <+ For normal drag speed-up>" 1880 PRINT#P9,"S.M.A.,km ",TAB(20),A0,TAB(40),A 1890 PRINT#P9,"Perigee Ht,km",TAB(20),A0*(1-E0)-R0,TAB(40),A*(1-E0)-R0 1900 PRINT#P9,"Apogee Ht,km",TAB(20),A0*(1+E0)-R0,TAB(40),A*(1+E0)-R0 1910 PRINT#P9,"Arg. Perigee ",TAB(20),W0,TAB(40),W 1920 PRINT#P9,"R. A. A. N. ",TAB(20),O0,TAB(40),O 1930 PRINT#P9,"Freq.,MHz ",TAB(20),F1 \ K9=9E9 \ K8=9E9 \ D0=2 1940REM 1950 REM ------------------------------ 1960 REM ****** Here follows the actual computation loop ****** 1970 FOR T=T7 TO T8 STEP T9 \ IF NOT K5 THEN 2000 1980 R6=X(0,0) \ X9=X(0,1) \ Y9=X(0,2) \ Z9=X(0,3) \ T6=X(0,4) 1990 E8=X(0,5) \ S8=X(0,6) \ C8=X(0,7) \ S9=X(0,8) \ C9=X(0,9) 2000 K7=INT(T) \ D=FNM(T) \ IF D0 THEN 2040 2010 REM--New orbit? If so, update the values for W & O for new orbit 2020 IF K=K9 THEN 2040 \ D=FNC(T) \ K8=9E9 \ K9=9E9 2030 REM--Above Elevation limit? New date? 2040 D=FNK(M) \ D=FNX(T) \ D=E9-E8 \ IF D<0 THEN 2120 2050 REM--IF IT'S BEEN BELOW HORIZON,THEN RESET TIME IN #726 2060 IF D0=0 THEN 2100 ELSE D0=2 2070 REM--ABOVE HORIZON--TEST ORBIT ## AND DATE 2080 IF K7=K8 AND K9=K THEN 2200 ELSE 2140 2090 REM--RESET TIME TO BE IN STEP 2100 T=T7 + T9*INT(((T-T7)/T9)-1) \ D0=1 \ GOTO 2000 2110 REM--THIS AD-HOC ALGORITHM SPEEDS UP THE STEPPING 2120 IF D0=1 THEN 2350 ELSE D= R5*D*D *1E-9 \ D0=0 2130 IF D > 0.2/N0 THEN T=T+ 0.2/N0 ELSE T=T+D \ GOTO 2350 2140 IF K=K9 THEN 2190 ELSE GOSUB 2390 \ K9=K 2150 PRINT#P9," U.T.C. AZ EL DOPPLER RANGE HEIGHT", 2160 PRINT#P9," LAT LONG PHASE" 2170 PRINT#P9,"HHMM:SS deg deg Hz km km ", 2180 PRINT#P9," deg deg <256>" 2190 PRINT#P9,TAB(10),"- - - DAY #",K7," - - - ORBIT #",K," - - -" 2200 K8=K7 \ T4=T-K7 \ S4=INT(T4*86400 + .5) \ H4=INT(S4/3600 + 1E-6) 2210 M4=INT((S4 - H4*3600) / 60 + 1E-6) \ S4=S4 - 3600*H4 - 60*M4 2220 T$=FNT$(H4) + FNT$(M4) + ":" + FNT$(S4) \ REM--Printable time string 2230 F9=-F1*1E6 * R8/C \ REM--F9=Doppler(Hz)=Frequency * Velocity/C 2240 PRINT#P9,T$, %6I,FNI(A9),FNI(E9), %7I,FNI(F9), 2250 PRINT#P9, %9I,FNI(R5),FNI(R-R0), %6F1,L5,W5, %6I,M9 2260 IF K5=0 THEN 2350 ELSE X(0,0)=R6 \ X(0,4)=T6 2270 FOR I=1 TO K5 \ REM--LOOP THRU OTHER STATIONS TO SEE IF VISIBLE 2280 R6=X(I,0) \ X9=X(I,1) \ Y9=X(I,2) \ Z9=X(I,3) \ T6=X(I,4) 2290 E8=X(I,5) \ S8=X(I,6) \ C8=X(I,7) \ S9=X(I,8) \ C9=X(I,9) 2300 D=FNX(T) \ IF E90 THEN 10090 10060 REM--The two cases for X=0-- 10070 IF Y>=0 THEN RETURN P1/2 ELSE RETURN 3*P1/2 10080 REM--Cases lying in Quadrants 1&4-- 10090 IF Y>=0 THEN RETURN ATN(Y/X) ELSE RETURN P2 + ATN(Y/X) 10100 FNEND ---------------------------- 10110 DEF FNC(T) 10120 REM--Routine to initialize the C(J,K) coordinate rotation matrix 10130 REM and other parameters associated with the orbital elements. 10139 REM--Given N0=Mean Motion, execute Lines 1014x 10140 REM--Input Elements give either A0= Semi-Major Axis in km OR 10142 REM N0=Mean Motion in Orbits/Day. Decide which was given and 10144 REM calculate the other, including drag corrections -- 10146 IF N0 > 0.1 THEN 10150 10148 N0=SQRT(G0 / (A0^3)) \ GOTO 10154 \ Calculate Mean Motion 10150 A0=((G0 / (N0*N0)) ^ (1/3)) \ REM--Calculate SMA 10152 REM--The following correct for drag in low-altitude satellites-- 10154 N=N0 + 2*(T-T0)*N1 \ REM--Mean Motion at T 10156 A =((G0 / (N *N )) ^ (1/3)) \ REM--SMA at Epoch T 10160 E2=1-E0^2 \ E1=SQRT(E2) \ Q0=M0/360 + K0 \ REM--Q0=Initial orbit phase 10170 REM--Account for nodal effects due to lumpy gravity field due to the 10180 REM flattened, oblate spheroidal, figure of the earth 10190 K2=9.95 * ((R0/A)^3.5) / (E2^2) 10200 REM--Update elements to current epoch & evaluate their SIN/COS's 10210 S1=SIN(I0*P0) \ C1=COS(I0*P0) \ REM--I0=Inclination (deg) 10220 O=O0 - (T-T0) * K2 * C1 10230 S0=SIN(O*P0) \ C0=COS(O*P0) \ REM--O=R.A.A.N. (degrees) 10240 W=W0 + (T-T0) * K2 * (2.5 * (C1^2) - .5) 10250 S2=SIN(W*P0) \ C2=COS(W*P0) \ REM--W=Argument of Perigee (degrees) 10260 REM--Set up coordinate rotation matrix for the current orbit 10270 C(1,1)= + (C2*C0) - (S2*S0*C1) \ C(1,2)= - (S2*C0) - (C2*S0*C1) 10280 C(2,1)= + (C2*S0) + (S2*C0*C1) \ C(2,2)= - (S2*S0) + (C2*C0*C1) 10290 C(3,1)= + (S2*S1) \ C(3,2)= + (C2*S1) \ RETURN 0 10300 FNEND ---------------------------- 10310 DEF FNM(T) 10320 REM--Function to evaluate M=MEAN ANOMALY in (0-2*PI) range, 10330 REM K=Perigee Passage Kounter (=Orbit ##) & M9=modulo 256 10340 REM orbital phase compatible with General Beacon telemetry. 10345 Q=Q0 + N0*(T-T0) + N1*((T-T0)^2) \ REM--REVISED TO INCLUDE DRAG 10350 K=INT(Q+1E-6) \ M9=INT((Q-K+1E-6)*256) \ M=(Q-K)*P2 \ RETURN 0 10360 FNEND ---------------------------- 10370 DEF FNK(M) 10380 REM--Routine to solve KEPLER's equation given M & return 10390 REM the Satellite's GEOCENTRIC coordinates. 10400 E=M + E0*SIN(M) + .5 * (E0^2) * SIN(2*M) \ REM--Initial trial value 10410 REM--Iteration loop to solve trancendental Kepler's equation-- 10420 S3=SIN(E) \ C3=COS(E) \ R3=1 - E0*C3 \ M1=E - E0*S3 10430 M5=M1-M \ IF ABS(M5) < 1E-6 THEN 10450 ELSE E=E - M5/R3 \ GOTO 10420 10440 REM--Now get satellite's XYZ coordinates-- 10450 X0=A*(C3-E0) \ Y0=A*E1*S3 \ R=A*R3 \ REM--In the plane of the orbit 10460 REM--Rotate from Orbit plane to INERTIAL CELESTIAL Coordinates. 10470 X1=X0*C(1,1)+Y0*C(1,2) \ Y1=X0*C(2,1)+Y0*C(2,2) \ Z1=X0*C(3,1)+Y0*C(3,2) 10480 REM--Rotate thru current GHA of ARIES, convert to GEOCENTRIC coordinates 10490 G7=T*G1 + G2 \ G7=(G7-INT(G7)) * P2 \ S7=-SIN(G7) \ C7=COS(G7) 10500 X= + (X1*C7) - (Y1*S7) \ Y= + (X1*S7) + (Y1*C7) \ Z=Z1 \ RETURN 0 10510 FNEND ---------------------------- 10520 DEF FNO(D) 10530 REM--Routine to evaluate OBSERVER's GEOCENTRIC coordinates, 10540 REM where X-axis=GREENWICH, Y-axis=VU-land & Z-axis=North pole. 10550 L8=L9*P0 \ S9=SIN(L8) \ C9=COS(L8) \ REM--Initial GEODETIC coordinates 10560 S8=SIN(-W9*P0) \ C8=COS(W9*P0) \ REM--W9=West Longitude 10570 REM--Following accounts for flattened oblate spheroidal earth-- 10580 R9=R0*(1 - (F/2) + (F/2) *COS(2*L8)) + H9/1000 \ REM--H9=height(meters) 10590 REM--Following makes L8 be the GEOCENTRIC latitude-- 10600 L8=ATN( (1-F)^2 * S9/C9 ) \ Z9=R9 * SIN(L8) 10610 X9=R9 * COS(L8) * C8 \ Y9=R9 * COS(L8) * S8 \ RETURN 0 10620 FNEND ---------------------------- 10630 DEF FNX(T) 10640 REM--Routine to extract all the parameters you might ever need. 10650 REM First get vector from observer to satellite-- 10660 X5= (X-X9) \ Y5= (Y-Y9) \ Z5= (Z-Z9) \ R5=SQRT(X5^2 + Y5^2 + Z5^2) 10670 REM--Finite difference the range (R5) to get the velocity (R8)-- 10680 IF T6<>T THEN R8=((R6-R5)/(T6-T)) / 86400 ELSE R8=-9E9 10690 R6=R5 \ T6=T \ REM--& Save current range & time for next time thru 10700 REM--Now rotate into OBSERVER's LOCAL coordinates 10710 REM where X8=North, Y8=East & Z8=Up (Left-Handed system!) 10720 Z8= + (X5*C8*C9) + (Y5*S8*C9) + (Z5*S9) 10730 X8= - (X5*C8*S9) - (Y5*S8*S9) + (Z5*C9) \ Y8= + (Y5*C8) - (X5*S8) 10740 S5=Z8/R5 \ C5=SQRT(1-S5*S5) \ E9=ATN(S5/C5) / P0 \ REM--E9=Elevation 10750 A9=FNA(X8,Y8) / P0 \ REM--FNA resolves correct quadrant for A9=Azimuth 10760 W5=360-FNA(X,Y) / P0 \ REM--W5=Sub-Satellite Point (SSP) W.Longitude 10770 B5=Z/R \ L5=ATN(B5 / (SQRT(1 - B5^2))) / P0 \ REM--L5=SSP Latitude 10780 RETURN 0 \ REM--Note R-R0=Satellite's altitude above mean spheroid. 10790 FNEND ---------------------------- 10800 REM--FNT$ Returns two character printable string corresponding to D-- 10810 DEF FNT$(D) = CHR$(48 + INT(D/10)) + CHR$(48 + D - 10 * INT(D/10)) 10820 REM ------------------------------ 10830 REM--FNI Rounds to nearest integer, accounting for sign-- 10840 DEF FNI(D)=SGN(D) * INT (ABS(D) + .5) 10850 END 15000 REM---ALL STATIONS 15001 DATA "KC3IL",74.84,40.1725,31,0 15002 REM CALL,LONG,LAT,ANT HT METERS,MIN EL. 15003 DATA "WB6DEB",122.53084,37.975,0,0 15004 DATA "K1HTV",76.79,38.97,100,0 15010 DATA "VK2SYD", 209,-34,0,0 15014 DATA "ZL1AOX",185,-37,0,0 15016 DATA "H44PT",200,-10,0,0 15022 DATA "RA0LFI",228,43.2,0,0 15026 DATA "JA1ANG",220,35.7,0,0 15040 DATA "UL7GBD",282,43,0,0 15044 DATA "SV1AB",326.2,38,0,0 15054 DATA "ZS1FE",341.5,-34,0,0 15056 DATA "GB2RS",0,51.5,0,0 15062 DATA "I1ROM",347.5,42,0,0 15070 DATA "PY2RIO",43,-23,0,0 15072 DATA "CE3XK",70,-34,0,0 16000 DATA "END" \ REM--END OF STATION LIST ============================================================================= Internet: p00123@psilink.com -or- kc3il@amsat.org AX.25net: kc3il@wa3dsp.#epa.pa.usa.na AMPR net: kc3il@wa3dsp.ampr.org Packet to Internet gateway: kc3il@w2xo.#wpa.pa.usa.na |-Internet to Packet gateway: bbs@w2xo.pgh.pa.us |- w/ first line of message: sp kc3il@wa3dsp.#epa.pa.usa.na ============================================================================= ARRL Life Member / AMSAT Life Member # 784 / QCWA Life Member # 23542 ============================================================================= "Now is the test of the boomerang, thought jewels polished and gleaming..." ============================================================================= "It is conceived of by him by whom It is not conceived of. He by whom It is conceived of, knows It not. It is not understood by those who [say they] understand It. It is understood by those who [say they] understand It not." (-KENA UPANISHAD 11.) ============================================================================= Date: Tue, 27 Apr 1993 18:08:08 -0700 (PDT) From: Dave Reeves Sender: Dave Reeves Reply-To: Dave Reeves Subject: Re: Shouldn't ask this, but I need W3IWI's Orbit code... To: David Moisan In-Reply-To: <9304261923.AA03699@silver.lcs.mit.edu> Mime-Version: 1.0 Content-Type: TEXT/PLAIN; CHARSET=US-ASCII On Mon, 26 Apr 1993, David Moisan wrote: > I shouldn't be a "me too" poster, but does anyone know where I can > find Tom Clark's Orbit source? I'd like to port it to my friend's > new programmable calc. > > Thanks and 73's, Dave N1KGH > > | David Moisan, N1KGH /^\_/^\ moisan@silver.lcs.mit.edu | > | 86 Essex St. Apt #204 ( o ^ o ) n1kgh@amsat.org | > | Salem. MA 01970-5225 | | | > \-----------------------------------------------------------------/ > > > > > Dave, You'll probably get twenty copies but here goes: 73, Dave KF6PJ kf6pj@amsat.org 100 ' Orbit Predictor from W3IWI program. 102 REM INSTRUCTIONS FOR DELAY OF LAUNCH: 104 REM ADD THE DIFFERENCE IN TIME TO THE EPOCH 106 REM INCREASE R.A.A.N. BY 15.04107 DEG/HOUR (RESOLVE TO 0-360) 108 REM IF LAUNCH AT 18:50 UTC, EPOCH:85/212.74725211;RAAN:112.4790 110 '========================================== 120 CLS: KEY OFF: CLOSE: CLEAR: OPTION BASE 1 130 DEFDBL P,G: DEFINT J 140 DIM MD%(13), MONTH$(12), MENU$(11,3), DAY$(7), S$(15) 150 DEF FNI!(X)=SGN(X)*INT(ABS(X)+.5) 160 DEF FNT$(D%)=CHR$(48+INT(D%/10))+CHR$(48+D%-10*INT(D%/10)) 165 GOSUB 2000 'Choose Satellite 170 GOSUB 10000 'Initialize 300 ' 310 FOR TINC=TSTART TO TEND STEP TSTEP 320 K7=INT(TINC): GOSUB 7010: IF K<>K9 THEN GOSUB 9040 330 GOSUB 7040: GOSUB 8010 340 T4=TINC-K7: S4=INT(T4*86400!+59): H4=INT(S4/3600!+.000001) 350 M4=INT((S4-H4*3600)/60+.000001): S4=S4-3600*H4-60*M4 360 KTZ=K7: JTZ=H4+MYTZ: IF JTZ<0 THEN JTZ=JTZ+24: KTZ=KTZ-1 370 IF E9KTZ THEN GOSUB 5020: K9=K: GOTO 410 ' Titles for new day 390 IF K=K9 THEN 450 ELSE K9=K: GOTO 450 400 ' 410 GOSUB 16020 ' Get Gregorian date 420 PRINT #1," ",TAB(12),"- - - ";G$;" - - -";TAB(1);" " 430 PRINT #1," PDT Az El Doppler Range Height Lat Long Phase Orbit" 440 PRINT #1,"HH:MM deg deg Hz km km deg deg <256> Number Mode" 450 K8=KTZ 460 T$=FNT$(JTZ)+":"+FNT$(M4): F9=-FREQ#*1000000!*R8/299792.5# 470 PRINT #1,T$; 480 PRINT #1,USING"######";FNI!(A9!);FNI!(E9-.5); 490 PRINT #1,USING"########";FNI!(F9);: PRINT #1," "; 500 PRINT #1,USING "########";FNI!(R5);FNI!(R-R0); 510 MODE$="SSTV" 550 PRINT #1,USING"######";FNI!(L5!);FNI!(W5);M9;K; 560 PRINT #1," ";MODE$ 570 NEXT 580 ' 590 PRINT #1,LF$;LF$ 600 CLOSE: END 2000 REM ** Choose Satellite ** 2400 S$(1)="Oscar-9" 2410 S$(2)="Oscar-10" 2420 S$(3)="Oscar-11" 2430 S$(4)="Oscar-12" 2440 S$(5)="Russian Sputnik-5" 2450 S$(6)="Russian Sputnik-7" 2460 S$(7)="Russian Sputnik-8" 2470 S$(8)="NOAA-6" 2480 S$(9)="NOAA-9" 2490 S$(10)="MIR" 2500 S$(11)="Salut-7" 2510 S$(12)="JAS-1" 2520 S$(13)="Challenger" 2530 PRINT"":PRINT:PRINT:PRINT"Satellite Selection Menu":PRINT 2540 FOR I=1 TO 9:PRINT" ";I;" - ";S$(I):NEXT 2550 FOR I=10 TO 13:PRINT I;" - ";S$(I):NEXT:PRINT 2560 INPUT"Select Entry # ";J 2570 IF J<1 OR J>13 GOTO 2560 2575 S$=S$(J) 2580 ON J GOSUB 2610,2730,2850,2970,3090,3210,3330,3450,3570,3690,3810,3930,4050 2585 CLS 2590 RETURN 2600 REM***** keplerian elements ******* 2610 REM********* oscar-9 ************** 2620 Y3=87:D3!=51.04303340#:H3=0:M3=0:S3=0:REM *yr,day,hms* 2630 REFORB#=29883#:REM *orbit number* 2640 MA#=54.1769#:REM *mean anomaly, deg.* 2650 MM#=15.29230010#:REM *mean motion, rev/day* 2660 INC#=97.6532#:REM *inclination, deg.* 2670 ECC#=.0001124#:REM *eccentricity* 2680 PERG#=305.9335#:REM *arg. of perigee, deg.* 2690 RAAN#=64.9846#:REM *r.a.a.n., deg.* 2700 FREQ#=145.825#:REM *beacon/downlink freq.* 2710 N1=1.388E-05:REM *orbit decay, rev/day/day* 2720 RETURN 2730 REM********* oscar-10 ************* 2740 Y3=87:D3!=23.13103915#:H3=0:M3=0:S3=0 2750 REFORB#=3047# 2760 MA#=183.0539# 2770 MM#=2.05877501# 2780 INC#=27.1040# 2790 ECC#=.6022546# 2800 PERG#=178.8941# 2810 RAAN#=38.6101# 2820 FREQ#=145.81# 2830 N1=-8.70000E-07 2840 RETURN 2850 REM********* oscar-11 ************* 2860 Y3=87:D3!=43.25738475#:H3=0:M3=0:S3=0 2870 REFORB#=15744# 2880 MA#=232 6872# 2890 MM#=14.62103402# 2900 INC#=98.1158# 2910 ECC#=.0014214# 2920 PERG#=127.5624# 2930 RAAN#=111.4578# 2940 FREQ#=145.826# 2950 N1=2.5E-07 2960 RETURN 2970 REM********* oscar-12 ************* 2980 Y3=87:D3!=51.94676829#:H3=0:M3=0:S3=0 2990 REFORB#=2391# 3000 MA#=11.2854# 3010 MM#=12.44393206# 3020 INC#=50.0167# 3030 ECC#=.0010878# 3040 PERG#=348.7725# 3050 RAAN#=22.6436# 3060 FREQ#=435.795# 3070 N1=-2.5E-07 3080 RETURN 3090 REM********* rs-5 ***************** 3100 Y3=87:D3!=48.14511068#:H3=0:M3=0:S3=0 3110 REFORB#=22737# 3120 MA#=79.2958# 3130 MM#=12.05052555# 3140 INC#=82.9558# 3150 ECC#=.0008619# 3160 PERG#=280.7135# 3170 RAAN#=333.5271# 3180 FREQ#=29.45#:REM also 29.330 3190 N1=1.2E-07 3200 RETURN 3210 REM********* rs-7 ***************** 3220 Y3=87:D3!=49.1767659#:H3=0:M3=0:S3=0 3230 REFORB#=22818# 3240 MA#=188.0657# 3250 MM#=12.08701241# 3260 INC#=82.9609# 3270 ECC#=.0021217# 3280 PERG#=172.0742# 3290 RAAN#=325.6488# 3300 FREQ#=29.5#:REM also 29.340 3310 N1=1.3E-07 3320 RETURN 3330 REM********* rs-8 ***************** 3340 Y3=87:D3!=51.36260285#:H3=0:M3=0:S3=0 3350 REFORB#=22736# 3360 MA#=43.4113# 3370 MM#=12.02960515# 3380 INC#=82.9647# 3390 ECC#=.0018991# 3400 PERG#=316.5469# 3410 RAAN#=335.8218# 3420 FREQ#=29.5# 3430 N1=1.2E-07 3440 RETURN 3450 REM********* noaa-6 *************** 3460 Y3=86:D3!=307.07480828#:H3=0:M3=0:S3=0 3470 REFORB#=38157# 3480 MA#=326.7897# 3490 MM#=14.24959779# 3500 INC#=98.5# 3510 ECC#=.0012787# 3520 PERG#=33.4085# 3530 RAAN#=318.1542# 3540 FREQ#=137.5# 3550 N1=.0000011 3560 RETURN 3570 REM********* noaa-9 *************** 3580 Y3=86:D3!=307.89222401#:H3=0:M3=0:S3=0 3590 REFORB#=9753# 3600 MA#=244.3005# 3610 MM#=14.11457336# 3620 INC#=99.0189# 3630 ECC#=.0015717# 3640 PERG#=115.9787# 3650 RAAN#=267.0251# 3660 FREQ#=137.62#: REM also 137.500 3670 N1=1.62E-06 3680 RETURN 3690 REM*********** mir **************** 3700 Y3=87:D3!=054.97294408#:H3=0:M3=0:S3=0 3710 REFORB#=5824# 3720 MA#=341.5803# 3730 MM#=15.70570004# 3740 INC#=51.613# 3750 ECC#=.0017911# 3760 PERG#=18.5479# 3770 RAAN#=16.9270# 3780 FREQ#=143.625# 3790 N1=4.6651E-04 3800 RETURN 3810 REM********* salyut-7 ************* 3820 Y3=87:D3!=054.94462430#:H3=0:M3=0:S3=0 3830 REFORB#=27863# 3840 MA#=202.3983# 3850 MM#=15.30894852# 3860 INC#=51.6133# 3870 ECC#=.0001314# 3880 PERG#=157.7087# 3890 RAAN#=80.2189# 3900 FREQ#=0#:REM no beacon listed 3910 N1=1.412E-05 3920 RETURN 3930 REM********* jas-1 **************** 3940 Y3=86:D3!=224.9077199075#:H3=0:M3=0:S3=0 3950 REFORB#=1# 3960 MA#=187.2703# 3970 MM#=2#:REM ?? mean motion not listed 3980 INC#=49.97123# 3990 ECC#=1.260921E-03 4000 PERG#=144.959# 4010 RAAN#=249.567# 4020 FREQ#=0#:REM no beacon listed 4030 N1=0 4040 RETURN 4050 REM********* shuttle ************** 4060 Y3=85:D3!=213.62777778#:H3=0:M3=0:S3=0 4070 REFORB#=44# 4080 MA#=317.7854# 4090 MM#=15.853206# 4100 INC#=49.6015# 4110 ECC#=.0003231# 4120 PERG#=328.8386# 4130 RAAN#=146.516# 4140 FREQ#=145.55# 4150 N1=0# 4160 RETURN 5000 ' 5010 ' Titles for new page 5020 IF PT$="S" THEN PRINT #1,CHR$(10);CHR$(10);CHR$(10);:RETURN 5030 PRINT #1,CTL$; 5040 CTL$=CHR$(30)+CHR$(27)+CHR$(54) '+CHR$(12) 10CPI,6LPI,LF 5050 PRINT #1," ",TAB(7);S$+" - Orbit Predictions"; 5060 PRINT #1,TAB(36);"For ";C$; 5070 P%=P%+1: PRINT #1," Pg. #";P%: PRINT #1," ": PRINT #1," " 5080 RETURN 7000 ' 7010 Q=MM#*(TINC-T0)+Q0+N1*((TINC-T0)^2): K=INT(Q): M9=INT((Q-K)*256): M=(Q-K)*P2 7020 RETURN 7030 ' 7040 E#=M+ECC#*SIN(M)+.5*(ECC#^2)*SIN(2*M) 7050 S3=SIN(E#): C3#=COS(E#): R3=1-ECC#*C3#: M1=E#-ECC#*S3: M5=M1-M 7060 IF ABS(M5)<.000001 THEN 7070 ELSE E#=E#-M5/R3: GOTO 7050 7070 X0=A!*(C3#-ECC#): Y0=SMA!*E1*S3: R=SMA!*R3: X1=X0*C11#+Y0*C12# 7080 Y1=X0*C21#+Y0*C22#: Z1=X0*C31#+Y0*C32#: G7=TINC*G1+G2 7090 G7=(G7-INT(G7))*P2: S7=-SIN(G7): C7#=COS(G7) 7100 X=+(X1*C7#)-(Y1*S7): Y=+(X1*S7)+(Y1*C7#): Z=Z1 7110 RETURN 8000 ' 8010 X5=(X-X9): Y5=(Y-Y9): Z5=(Z-Z9): R5=SQR(X5^2+Y5^2+Z5^2) 8020 IF T6<>TINC THEN R8=((R6-R5)/(T6-TINC))/86400! ELSE R8=-9.000001E+09 8030 R6=R5: T6=TINC 8040 Z8=+(X5*C8#*C9#)+(Y5*S8*C9#)+(Z5*S9):X8=-(X5*C8#*S9)-(Y5*S8*S9)+(Z5*C9#) 8050 Y8=+(Y5*C8#)-(X5*S8): S5=Z8/R5: C5#=SQR(1-S5*S5): E9=ATN(S5/C5#)/P0 8060 V1=X8: V2=Y8: GOSUB 9010: A9!=VR/P0: V1=X: V2=Y: GOSUB 9010 8070 W5=360-VR/P0: B5!=Z/R: L5!=ATN(B5!/(SQR(1-B5!^2)))/P0 8080 RETURN 9000 ' 9010 IF V1<0 THEN VR=P1+ATN(V2/V1):RETURN ELSE IF V1>0 THEN GOTO 9030 9020 IF V2>=0 THEN VR=P1/2:RETURN ELSE VR=3*P1/2: RETURN 9030 IF V2>=0 THEN VR=ATN(V2/V1):RETURN ELSE VR=P2+ATN(V2/V1): RETURN 9040 IF MM#>.1 THEN SMA!=((G0/(MM#^2))^(1/3)) 9050 IF MM#>.1 THEN GOTO 9090 9060 MM#=SQR(G0/(SMA!^3)) 9070 GOTO 9100 9080 ' 9090 SMA!=((G0/(MM#^2))^(1/3)) 9100 N=MM#+2*(TINC-T0)*N1: A!=((G0/(N^2))^(1/3)): E2=1-ECC#^2: E1=SQR(E2) 9110 Q0=MA#/360+REFORB#: K2!=9.95*((R0/A!)^3.5)/(E2^2): S1=SIN(INC#*P0) 9120 C1#=COS(INC#*P0): O=RAAN#-(TINC-T0)*K2!*C1#: S0=SIN(O*P0): C0#=COS(O*P0) 9130 W=PERG#+(TINC-T0)*K2!*(2.5*(C1#^2)-.5): S2=SIN(W*P0): C2#=COS(W*P0) 9140 C11#=+(C2#*C0#)-(S2*S0*C1#): C12#=-(S2*C0#)-(C2#*S0*C1#) 9150 C21#=+(C2#*S0)+(S2*C0#*C1#): C22#=-(S2*S0)+(C2#*C0#*C1#) 9160 C31#=+(S2*S1): C32#=+(C2#*S1) 9170 RETURN 10000 ' 10010 ' Initialization 10020 ' 10030 CTL$=CHR$(27)+CHR$(40)+CHR$(30)+CHR$(27)+CHR$(54) '10CPI, 6LPI for Oki ML92 10040 P%=0: P1=3.1415926535#: R0=6378.16: LF$=CHR$(10) 10050 F=298.25: G0=75369793000000# 10060 G1=1.0027379093#: P2=2*P1: P0=P1/180: F=1/F 10070 '*************************************** 10080 FOR J=1 TO 13: READ MD%(J): NEXT ' Days in months table 10090 FOR J=1 TO 12: READ MONTH$(J): NEXT ' Names of months 10100 FOR J=1 TO 7 : READ DAY$(J): NEXT ' Names of days 10110 FOR J=1 TO 11: READ MENU$(J,1): NEXT ' Parameter titles 10120 '*************************************** 10130 MENU$(6,2)="34.209" ' Latitude 10140 MENU$(7,2)="118.663" ' Longitude 10150 MENU$(8,2)="975" ' Antenna height above MSL 10160 MENU$(9,2)="-2" ' Minimum elevation angle 10170 MENU$(10,2)="WA6BYE" 10180 MENU$(11,2)="-8" ' Time zone offset, PST 10190 '*************************************** 10200 GOSUB 11020 ' Process menu options 10210 GOSUB 12020 ' Calculate station location 10220 GOSUB 5020 ' Print headings 10230 GOSUB 14130 ' Set Keplerian elements 10240 GOSUB 15010 ' Print orbital parameters 10250 RETURN 11000 ' 11010 ' Process menu options 11020 D$=DATE$: MID$(D$,7,4)=MID$(D$,9,2)+" ": MENU$(1,2)=D$ 11030 MENU$(2,2)="0000": MENU$(3,2)="24": MENU$(4,2)="30": MENU$(5,2)="S" 11040 LOCATE 1,25,0: PRINT S$+" - Orbit Predictor" 11050 FOR J=1 TO 11 11060 LOCATE J+5,12,0: PRINT MENU$(J,1): LOCATE J+5,52: PRINT MENU$(J,2) 11070 NEXT J 11080 LOCATE 20,25,0: PRINT "Enter ESCape to begin" 11090 ' 11100 C%=52: R%=6 11110 IF R%>16 THEN R%=6 11120 IF R%<6 THEN R%=16 11130 IF C%<52 THEN C%=65: R%=R%-1: GOTO 11110 11140 IF C%>65 THEN C%=52: R%=R%+1: GOTO 11110 11150 LOCATE R%,C%,1 11160 A$=INKEY$: IF A$="" THEN 11160 11170 IF A$=CHR$(13) THEN R%=R%+1: C%=52: GOTO 11110 11180 IF A$=CHR$(27) THEN GOTO 11340 ' ESCape to exit 11190 IF A$=CHR$(8) THEN C%=C%-1: GOTO 11150 11200 IF LEN(A$)=1 THEN 11270 11210 J=ASC(MID$(A$,2,1)) 11220 IF J=77 THEN C%=C%+1: GOTO 11110 ' cursor right 11230 IF J=75 THEN C%=C%-1: GOTO 11110 ' cursor left 11240 IF J=72 THEN R%=R%-1: GOTO 11110 ' cursor up 11250 IF J=80 THEN R%=R%+1: GOTO 11110 ' cursor down 11260 BEEP: LOCATE ,42: GOTO 11160 11270 C%=POS(0): R%=CSRLIN 11280 IF C%=52 THEN MENU$(R%-5,2)=STRING$(14," "): PRINT STRING$(14," ") 11290 LOCATE R%,C%,0: PRINT A$ 11300 MID$(MENU$(R%-5,2),C%-51,1)=A$: C%=C%+1 11310 GOTO 11110 11320 ' 11330 ' Convert input parameters 11340 LAT!=VAL(MENU$(6,2)): LON!=VAL(MENU$(7,2)): ANTHT#=VAL(MENU$(8,2)) 11350 MINEL#=VAL(MENU$(9,2)): C$=MENU$(10,2) : MYTZ=VAL(MENU$(11,2)) 11360 D$=MENU$(1,2): MM$=MID$(D$,1,2): DD$=MID$(D$,4,2): YY$=MID$(D$,7,2) 11370 MN=VAL(MM$): D%=VAL(DD$): Y=VAL(YY$): Y=Y+1900!: Y=Y/100 11380 Y=INT(100*(Y-INT(Y))+.1) 11390 IF Y/4=INT(Y/4) THEN F9%=1 ELSE F9%=0 11400 T$=FNT$(Y)+"/"+FNT$(MN)+"/"+FNT$(D%)+" at " 11410 D8!=D%+MD%(MN) 11420 IF MN>2 THEN D8!=D8!+F9%: FOR J=2 TO 13: MD%(J)=MD%(J)+F9%:NEXT J 11430 D$=MENU$(2,2): IF LEN(D$)=0 THEN D$="0000" 11440 HR%=VAL(MID$(D$,1,2)): M=VAL(MID$(D$,3,2)) 11450 TSTART=D8!+(HR%/24#)+(M#/1440) : T$=T$+FNT$(HR%)+FNT$(M)+":00" 11460 D$=MENU$(3,2): IF LEN(D$)=0 THEN D$="24" 11470 H1=VAL(D$): M1=0#: TEND=TSTART+(H1/24#)+(M1/1440#) 11480 M2=VAL(MENU$(4,2)): IF M2=0 THEN M2=30 11490 TSTEP=M2: TSTEP=TSTEP/1440# 11500 PT$=MENU$(5,2):PT$=CHR$(ASC(PT$) AND &HDF) 11510 IF PT$="P" THEN FIL$="LPT1:" ELSE FIL$="SCRN:" 11520 OPEN FIL$ FOR OUTPUT AS #1 11530 RETURN 12000 ' 12010 ' Calculate station location 12020 ON (Y-78) GOSUB 12100,12110,12120,12130,12140,12150,12160,12165,12167,12170 12030 L8!=LAT!*P0:S9=SIN(L8!):C9#=COS(L8!): S8=SIN(-LON!*P0): C8#=COS(LON!*P0) 12040 R9=R0*(1-(F/2)*COS(S*L8!))+ANTHT#/1000: L8!=ATN((1-F)^2*S9/C9#) 12050 Z9=R9*SIN(L8!): X9=R9*COS(L8!)*C8#: Y9=R9*COS(L8!)*S8 12060 RETURN 12100 G2=.2751843198# :RETURN 12110 G2=.2745212008# :RETURN 12120 G2=.2765959911# :RETURN 12130 G2=.2759328721# :RETURN 12140 G2=.2752697531# :RETURN 12150 G2=.2746066342# :RETURN 12160 G2=.2766814244# :RETURN 12165 G2=.2755417421# :RETURN 12167 G2=.2755417421# :RETURN 12170 PRINT"Unable to find year ";Y;" in sidereal time table.": STOP 14000 ' 14010 ' Current orbit parameters 14130 SMA!=0! 'Semi Major Axis 14140 IF Y3<>Y THEN PRINT TAB(1);"Elements not from current year.":STOP 14150 IF D3!=INT(D3!) THEN T0=D3!+H3/24+M3/1440+S3/86400! ELSE T0=D3! 14160 RETURN 15000 ' 15010 CLS ' PRINT ORBITAL PARAMETERS 15020 PRINT #1," ": PRINT #1,"Orbital Elements For "+S$ 15030 PRINT #1,TAB(10);'Reference ID = MH 2-12-84": PRINT #1," " 15040 PRINT #1,TAB(1);"Reference epoch ";Y3;" + ";T0 15050 PRINT #1,TAB(1);"Starting epoch ";Y;" + ";TSTART;" = ";T$ 15060 PRINT #1," " 15070 PRINT #1,"Parameter Reference Starting" 15080 TINC=TSTART: GOSUB 9040: GOSUB 7010 15090 PRINT #1,"Orbit Number ";REFORB#;TAB(40)K 15100 PRINT #1,"Mean Anomoly ";MA#;TAB(40)M/P0 15110 PRINT #1,"Inclination ";INC# 15120 PRINT #1,"Eccentricity ";ECC# 15130 PRINT #1,"Mean Motion ";MM# 15140 PRINT #1,"S.M.A., Km ";SMA! 15150 PRINT #1,"Arg. Perigee ";PERG#;TAB(40) W 15160 PRINT #1,"R.A.A.N. ";RAAN#;TAB(40) O 15170 PRINT #1,"Ref Freq. MHz ";FREQ# 15180 PRINT #1," " : PRINT #1,"Station parameters:" 15190 PRINT #1," Latitude ";LAT! 15200 PRINT #1," Longtitude ";LON! 15210 PRINT #1," Min. elevation ";MINEL# 15220 PRINT #1," Ant. height ";ANTHT# 15230 PRINT #1," Time zone ";MYTZ 15240 K9=9.000001E+09: K8=9.000001E+09 15250 RETURN 16000 ' 16010 ' Given KTZ, return G$ Gregorian date 16020 FOR J=2 TO 13: IF KTZ>MD%(J) THEN NEXT J 16030 ' 16040 ' Compute Day of Week using Zeller's Congruence 16050 MM=J-1: DD=KTZ-MD%(J-1): YY=VAL(YY$) 16060 IF YY<100 THEN YY=YY+1900 ' assume 20th century 16070 M=MM:D=DD:Y=YY ' use separate work area 16080 IF M<3 THEN M=M+12:Y=Y-1 ' Jan & Feb = 13th & 14th month of prior year 16090 J=INT(D+2*M+INT(3*(M+1)/5)+Y+INT(Y/4)-INT(Y/100)+INT(Y/400)+1) 16100 J=J MOD 7 16110 IF MM>12 THEN MM=MM-12: YY=YY+1 16120 G$=DAY$(J+1)+", "+MONTH$(MM)+STR$(DD)+","+STR$(YY) 16130 RETURN 20000 ' 20010 DATA 0,31,59,90,120,151,181,212,243,273,304,334,365 20020 DATA "January","February","March","April","May","June" 20030 DATA "July","August","September","October","November","December" 20040 DATA "Sunday","Monday","Tuesday","Wednesday","Thursday" 20050 DATA "Friday","Saturday" 20060 DATA "Starting date, MM/DD/YY","Starting time, HHMM" 20070 DATA "Duration, hours","Interval, mins" 20080 DATA "Screen or Printer, S/P","Latitude, degrees" 20090 DATA "Longitude, degrees","Antenna height, feet" 20100 DATA "Minimum elevation, degrees","Call sign(s)" 20110 DATA "Time zone offset"