' PREDICTION 
' First Amiga version by G3WDI
' Added to by A Hewat G8NTH

' Last mod 26/09/88 by ACH

'* this program needs to be up dated regularly 
'* to put in new elements STOP the program 
'* and type LIST keps: 
'* the order of the data items is given at the end 
'* It is important to remember the three ,0 after the day decimal day
'* Also important is the ,0 after the Mean Motion 
 
 
 
 
  DIM T$(20),S$(40),I$(40),C(3,2):P=0
DEFDBL A-Z

C8$=CHR$(10)+CHR$(10)+CHR$(10)+CHR$(10)
50 C9$=CHR$(12)+CHR$(7):PRINT C9$
PRINT "AMSAT ORBITAL PREDICTION PROGRAM"
PRINT  "COPYRIGHT 1980  DR. TOM CLARKE W3IWI"
PRINT  "TRANSLATED FOR AMSAT UK BY G3WDI"
PRINT "COPYRIGHT AMSAT-UK. 1988"
REM SET STARTING DATA
INPUT "start year = ",Y:Y=Y/100:Y=INT(100* (Y-INT(Y)) + 0.1)
IF Y/4=INT(Y/4) THEN F9=1 :ELSE F9=0
INPUT "     month = ",M
INPUT "       day = ",d
 DEF FNT$(d)= CHR$(48+INT(d/10))+CHR$(48+d-10*INT(d/10))
T$=FNT$(d)+"/"+FNT$(M)+"/"+FNT$(Y)+" at "
RESTORE daysofyear:
FOR I= 1 TO M:READ D9:NEXT:D8=d+D9
IF M>2 THEN D8=D8+F9 :PRINT "    day = ",D8
PRINT
INPUT"Start : UTC hours= ",H
INPUT"             Mins= ",M
T7=D8+H/24+M/1440
T$=T$+FNT$(H)+FNT$(M)+":00"
PRINT ''
INPUT"duration hours= ",H1
INPUT"          mins= ",M1:T8=T7+H1/24+M1/1440
 PRINT ''
INPUT "TIMESTEP: Min= ",M2:T9=M2/1440
PRINT '
PRINT "    from ",:PRINT USING"######.####";T7,:PRINT " to ";:PRINT USING"######.####";T8
RESTORE Physicaldata:
READ P1,C,R0,F,G0,G1
P2=2*P1:P0=P1/180:F=1/F
Getyear:
READ Y1,G2:IF Y1=Y THEN statdata:
IF Y1>0 THEN Getyear:
IF Y1=0 THEN PRINT "UNABLE TO FIND YEAR:STOP
statdata:
RESTORE stationdata:
READ L9,W9,H9,C$,E8
GOSUB observerdata: PRINT''
DEFDBL Q,K
GOSUB pobsdata:
PRINT "-------------satellite selection menu -------------------------"
I=0:RESTORE keps:
mnu:
READ S$:IF S$="END" THEN choice: :ELSE READ I$:I=I+1
PRINT "Entry # ";I;" for ";S$:PRINT "   ID= ";I$
READ d,d,d,d,d,d,d,d,d,d,d,d,d,d : GOTO mnu:
choice:
INPUT "Select entry # ",J:RESTORE keps:
IF J<1 OR J>I THEN choice:
FOR I= 1 TO J:READ S$,I$,Y3,D3,H3,M3,S3
READ K0,M0,N0,A0,I0,E0,W0,O0,F1:NEXT
PRINT :PRINT  "Doppler calculated for frequency = ",F1," MHz "
INPUT "Change frequency to (0 for default) ",d
IF d<>0 THEN F1=d
IF Y3=Y THEN 500 :ELSE PRINT "ELEMENTS NOT FROM CURRENT YEAR":STOP
500 IF D3=INT(D3) THEN T0=D3+H3/24+M3/1440+S3/86400 :ELSE T0=D3
PRINT :PRINT "Orbital Elements for ";S$
PRINT "    Reference Id.= ";I$:PRINT 
PRINT "Reference Epoch = ";Y3;" + ";T0
PRINT "Starting  Epoch = ";Y;" + ";T7;" = ";T$
PRINT 
PRINT "Parameter";SPC(9);"Reference";SPC(22);"Starting"
T=T7
GOSUB matrix: : GOSUB meananomoly: 
PRINT "Orbit number ";SPC(10);K0;SPC(10);K
PRINT "Mean Anomaly ";SPC(10);M0;SPC(10);M/P0
PRINT "Inclination  ";SPC(10);I0
PRINT "Eccentricity ";SPC(10);E0
PRINT "Mean Motion  ";SPC(10);N0
PRINT "S.M.A.,Km    ";SPC(10);A0
PRINT "Arg.Perigee  ";SPC(10);W0;SPC(10);W
PRINT " R. A. A. N. ";SPC(10);O0;SPC(10);O
PRINT "Freq. Mhz.   ";SPC(10);F1

K9=5.47102e-315:K8=5.47102e-315:

maincomputationloop:
FOR T=T7 TO T8 STEP T9:K7= INT(T):GOSUB meananomoly
IF K=K9 THEN 720:GOSUB matrix::K8=5.47102e-315:K9=5.47102e-315
720 GOSUB Keplerseq: 
GOSUB parameters: 
IF E9<E8 THEN 830 
IF K7=K8 AND K9=K THEN 770
730 IF K=K9 THEN 760 :ELSE GOSUB pager:
K9=K

740 PRINT " U.T.C.    AZ    EL   DOPPLER  RANGE    HEIGHT   LAT   LONG   PHASE"
750 PRINT "HHMM:SS   deg   deg      Hz      km         km   deg    deg   (256)"
760 PRINT SPC(10);"---DAY# ";K7;"---ORBIT # ";K;"----"
770 K8=K7:T4=T-K7:S4=INT(T4*86400 +0.5):H4=INT(S4/3600+1e-06)
780 M4=INT((S4-H4*3600)/60+1e-06):S4=INT(S4-3600*H4-60*M4)
790 T$=FNT$(H4)+FNT$(M4)+":"+FNT$(S4)
800 F9=-F1*1e+06*R8/C  : DEF FNI(d)=SGN(d)*INT(ABS(d)+0.5)
810 PRINT T$;:PRINT USING "######";FNI(A9);FNI(E9);:PRINT  USING"#######";FNI(F9);
PRINT USING "#########";FNI(R5);FNI(R-R0);:PRINT  USING"######";FNI(L5);FNI(W5);M9
830 NEXT:PRINT : INPUT"PRESS RETURN TO CONTINUE:",d$:IF d$="s" THEN STOP
GOTO 50
pager:
PRINT  C9$;C$;"  lat=";L9;"  W.Long=";W9;"  Ht. = "H9;
P=P+1:PRINT SPC(10);" Pg.="P
PRINT  SPC(10);"---MINIMUM ELEVATION= ";E8;"Deg.--"
PRINT  :RETURN


STOP






observerdata:
L8=L9*P0
S9=SIN(L8)
C9=COS(L8)
S8=SIN(-W9*P0):C8=COS(W9*P0)
R9=R0*(1-(F/2)+(F/2)*COS(2*L8))+H9/1000
L8=ATN((1-F)^2*S9/C9):Z9=R9*SIN(L8)
X9=R9*COS(L8)*C8:Y9=R9*COS(L8)*S8:RETURN
pobsdata:
PRINT C$;"   Lat.="L9;" W.Long.=";W9;"   ht.=";H9
PRINT TAB(5)"---Minimum elevation -";E8;"Deg--"
PRINT 
RETURN
matrix:
IF N0>0.1 THEN A0=((G0/(N0^2))^(1/3))
IF N0<=0.1 THEN N0=SQR(G0/(A0^3))
E2=1-E0^2:E1=SQR(E2):Q0=M0/360+K0
K2=9.95*((R0/A0)^3.5)/(E2^2)
S1=SIN(I0*P0):C1=COS(I0*P0)
O=O0-(T-T0)*K2*C1
S0=SIN(O*P0)
 C0=COS(O*P0)
W=W0+(T-T0)*K2*(2.5 * (C1^2)-0.5)
S2=SIN(W*P0)  
C2=COS(W*P0)
C(1,1)=+(C2*C0)-(S2*S0*C1) 
C(1,2)=-(S2*C0)-(C2*S0*C1)
C(2,1)=+(C2*S0)+(S2*C0*C1) 
C(2,2)=-(S2*S0)+(C2*C0*C1)
C(3,1)=+(S2*S1)           
C(3,2)=+(C2*S1) : RETURN
meananomoly:
Q=N0*(T-T0)+Q0:K=INT(Q):M9=INT((Q-K)*256):M=(Q-K)*P2:RETURN


keps:
DATA "Cosmos 1766"
DATA "LES"
DATA 87,230.74062106,0,0,0
DATA 999
DATA 88.5386
DATA 13.8338111720,0
DATA 82.5567
DATA 0.0011
DATA 271.1364
DATA 29.0651
DATA 137.4

DATA "MIR"
DATA "NASA"
DATA 89,44.56559894,0,0,0
DATA 17185
DATA 89.2202
DATA 15.71251869,0
DATA 51.6221
DATA 0.0011774
DATA 270.7082
DATA 266.8741
DATA 143.625

DATA "METEOR 2/16"
DATA "NASA"
DATA 89,41.77348822,0,0,0
DATA 7503
DATA 174.0910
DATA 13.83420134,0
DATA 82.5593
DATA 0.0011661
DATA 186.0118
DATA 319.7003
DATA 137.850

DATA "METEOR 2/15"
DATA "NASA"
DATA 89,36.98407380,0,0,0
DATA 10547
DATA 88.0889
DATA 13.8366481,0
DATA 82.4678
DATA 0.0012387
DATA 271.8885
DATA 259.0124
DATA 137.85

DATA"NOAA-10"
DATA"NASA"
DATA 89,41.58262838,0,0,0
DATA 12593
DATA 102.3104
DATA 14.2280426,0
DATA 98.651
DATA 0.0012704
DATA 257.6903
DATA 74.4169
DATA 137.5

DATA "NOAA-9"
DATA "NASA"
DATA 88,22.06278960,0,0,0
DATA 16019
DATA 45.5306
DATA 14.11545253,0
DATA 99.0818
DATA 0.0016
DATA 314.4627
DATA 354.8969
DATA 137.62

DATA "METEOR 30"
DATA "NASA"
DATA 87,346.29912200,0,0,0
DATA 40878
DATA 72.8717
DATA 14.98493528,0
DATA 97.7216
DATA 0.0041
DATA 286.7977
DATA 72.4335
DATA 137.010

DATA "SALYUT 7"
DATA "NASA"
DATA 87,350.91797151,0,0,0
DATA 32398
DATA 207.1197
DATA 15.31746328,0
DATA 51.6122
DATA 0.0001
DATA 152.9965
DATA 94.1426
DATA 20.005
DATA "END"                     
'ADD NEW DATA HERE EPOCH ORBIT MA MM INC ECC ARG RAAN FQ




Physicaldata:
DATA 3.1415926535,2.997925E5,6378.160,298.25:'PI/VEL LIGHT/RAD EARTH/FLATTENING
DATA 7.5369793E13,1.0027379093:'GM OF EARTH/SID/SOLAR TIME
siderealdata:
DATA 86 ,.27
DATA 87,0.2753552680
DATA 88,0.2746921630
DATA 89,0.2767669670
DATA 90,0.2767669760
DATA 0,0
stationdata:
' Observers data   lat, long, height, call, min elevation
DATA 51.268333,0.56333,,"G8NTH",1
 
 


daysofyear:
DATA 0,31,59,90,120,151,181,212,243,273,304,334,365

Keplerseq:
E=M+E0*SIN(M)+0.5*(E0^2)*SIN(2*M)
10420 S3=SIN(E):C3=COS(E):R3=1-E0*C3:M1=E-E0*S3
M5=M1-M:IF ABS(M5)<1e-06 THEN 10450 :ELSE E=E-M5/R3:GOTO 10420
10450 X0=A0*(C3-E0):Y0=A0*E1*S3:R=A0*R3
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)
G7=T*G1+G2:G7=(G7-INT(G7))*P2:S7=-SIN(G7):C7=COS(G7)
X=+(X1*C7)-(Y1*S7):Y=+(X1*S7)+(Y1*C7):Z=Z1:RETURN

parameters:
X5=(X-X9):Y5=(Y-Y9):Z5=Z-Z9:R5=SQR(X5^2+Y5^2+Z5^2)
IF T6<>T THEN R8=((R6-R5)/(T6-T))/86400 
IF T6=T THEN R8 = -9e+09
R6=R5:T6=T
Z8=+(X5*C8*C9)+(Y5*S8*C9)+(Z5*C9)
X8=-(X5*C8*S9)-(Y5*S8*S9)+(Z5*C9):Y8=(Y5*C8)-(X5*S8)
S5=Z8/R5:C5=SQR(1-S5*S5):E9=ATN(S5/C5)/P0
A9=ATN(Y8/X8)/P0
W5=(ATN(Y/X))*57.3
IF X<0 THEN W5=180-W5
IF X>0 AND Y<0 THEN W5=-W5
IF  X>0 AND Y>0 THEN W5=360-W5
IF X=0 AND Y>=0 THEN W5=270
IF X=0 AND Y<0 THEN W5=90
IF X8<0 THEN A9=A9+180
IF X8>0 AND Y8<0 THEN A9=360+A9
IF X8=0 AND Y8>=0 THEN A9=90
IF X8=0 AND Y8<0 THEN A9=270
B5=Z/R:L5=ATN(B5/(SQR(1-B5^2)))/P0
RETURN