n2ic@druxm.UUCP (LondonSM) (02/13/84)
Here is a program for MUF estimation which I have found to give
considerably better results than MINIMUF.
It is reputedly from one of the RSGB publications.
The version listed below runs on Commodore 64 BASIC.
I find the "Short Path" results to be rather accurate, but the
"Long Path" results somewhat questionable.
Any comments, questions, etc. would be welcome.
Steve London, N2IC/0
90 P1=3.141592654:J9=P1/180
100 INPUT "ROUTE";L$
110 REM North Latitude is positive, south latitude is negative
130 INPUT "LATITUDE";L6
135 REM West Longitude is negative, East Longitude is positive
140 INPUT "LONGITUDE";W6
160 W6=-W6
165 REM Input month as a number from 1 to 12
170 INPUT "MONTH";M0
180 INPUT "SOLAR FLUX";SF
185 S0=625*(SQR((0.73)^2-.0032*(65-SF))-.73)
190 REM L5 is hard-coded for the originating station latitude
200 L5=40
205 REM W5 is hard coded for the originating station longitude
210 W5=-105
220 P$="S"
230 IF ABS(W5-W6)<120 THEN 260
240 INPUT "LONG OR SHORT PATH";P$
260 D0=15
270 W1=-(W5-180*(1+SGN(W5-.001)))*J9
280 W2=-(W6-180*(1+SGN(W6-.001)))*J9
290 L1=L5*J9:L2=L6*J9
300 PRINT CHR$(147):PRINT "ROUTE: ";L$;" MONTH:";M0;" S.S. NO: ";S0
330 W3=W2-W1+.001
340 W3=P1*(1-SGN(W3))+W3
360 H1=SIN(L1)*SIN(L2)+COS(L1)*COS(L2)*COS(W3)
370 G1=ATN(SQR(1-H1*H1)/H1)+P1/2*(1-SGN(H1))
380 IF P$="L" THEN G1=P1+P1-G1
400 H0=INT(1.59*G1)+1
420 H9=(SIN(L2)-H1*SIN(L1))/SIN(G1)/COS(L1)
430 H9=ATN(SQR(1-H9*H9)/H9)+P1/2*(1-SGN(H9))
440 H9=H9*SGN(W3-P1)*SGN(P1-G1)
450 H9=H9+P1*(1-SGN(H9))
460 PRINT "LONG/SHORT: ";P$;" BEARING: ";INT(H9/J9+.5):PRINT
480 PRINT"GMT HPF MUF LUF":PRINT
500 Y6=ATN(1/TAN(G1/(H0+1))-.952/SIN(G1/(H0+1)))
510 IF Y6<.314 THEN Y6=.314
520 Y6=1/SQR(1-.965*COS(Y6)^2)
530 Y1=.0172*(10+(M0-1)*30.4+D0)
540 Y2=.409*COS(Y1)
550 Y1=.13*SIN(Y1)+.156*SIN(Y1+Y1)
570 H9=(SIN(L1)-COS(G1)*SIN(L2))/SIN(G1)/COS(L2)
590 M9=SIN(2.5*G1/H0)
600 M9=1+2.5*M9*SQR(M9)
610 K1=1-.5/H0
620 FOR N=0 TO 1
640 L9=COS(G1*K1)*SIN(L2)+SIN(G1*K1)*COS(L2)*H9
650 L0=P1/2-(ATN(SQR(1-L9*L9)/L9)+P1/2*(1-SGN(L9)))
660 W0=(COS(G1*K1)-SIN(L2)*L9)/COS(L2)/COS(L0)
670 W0=ATN(SQR(1-W0*W0)/W0)+P1/2*(1-SGN(W0))
680 W0=P1-SGN(P1-G1*K1)*(P1-W0)
690 W0=W3+W0*SGN(W3-P1)*SGN(P1-G1)+W1-.001
700 W0=W0-P1*(1-SGN(P1+P1-W0))
720 T0=3.82*W0+12+Y1
730 T0=T0-12*(1+SGN(T0-24))*SGN(ABS(T0-24))
740 IF COS(L0+Y2)>-.26 THEN 780
750 T1(N)=0
760 GOTO 900
780 T1(N)=(-.26+SIN(Y2)*L9)/(COS(Y2)*COS(L0)+.001)
790 T1(N)=12-ATN(T1(N)/SQR(ABS(1-T1(N))))*24/P1
810 T7=T0-T1(N)/2
820 T3(N)=T7+12*(1-SGN(T7))*SGN(ABS(T7))
840 T7=T0+T1(N)/2
850 T4(N)=T7-12*(1+SGN(T7-24))*SGN(ABS(T7-24))
860 C0(N)=ABS(COS(L0+Y2))
880 T9(N)=9.7*(C0(N)^8)
890 IF T9(N)<.1 THEN T9(N)=.1
900 K1=1-K1
920 M(N)=M9*.75*((12/T1(N)-1)*SGN(INT(12/T1(N)))+1)
930 M(N)=M(N)*(1+S0/100*(1-(T1(N)/12-1)*SGN(INT(T1(N)/12))))
940 L9=ABS(L0+.21*SIN(W0+.35))
950 G2=.5
960 IF L9<P1/4 THEN 990
970 M(N)=M(N)*(1-.1*(1+COS(L9*4)))
980 G2=.2
990 L(N)=SIN(L9*4)*G2
1010 G8(N)=P1*T9(N)/T1(N)
1020 T7=T1(N)/T9(N)
1030 IF T7>85 THEN T7=85
1040 G7(N)=C0(N)*G8(N)*(EXP(-T7)+1)
1050 G6(N)=G7(N)*EXP((T1(N)-24)/2)
1060 NEXT N
1070 FOR T5=0 TO 23
1080 J9=100
1090 K9=0
1100 FOR N=0 TO 1
1110 G0=0
1120 G3=P1/2
1130 IF T1(N)=0 THEN 1320
1140 IF T4(N)<T3(N) THEN 1190
1160 IF (T5-T3(N))*(T4(N)-T5)>0 THEN 1210
1170 GOTO 1370
1190 IF (T5-T4(N))*(T3(N)-T5)>0 THEN 1370
1210 T6=T5+12*(1+SGN(T3(N)-T5))*SGN(ABS(T3(N)-T5))
1220 G4=P1*(T6-T3(N))/T1(N)
1230 T8=(T3(N)-T6)/T9(N)
1240 IF ABS(T8)>85 THEN T8=85*SGN(T8)
1250 G0=C0(N)*(SIN(G4)+G8(N)*(EXP(T8)-COS(G4)))
1260 G3=P1/2
1270 IF T6-T3(N)>T1(N)/2+3 THEN 1290
1280 G3=(T6-T3(N))/(T1(N)/2+3)*G3
1290 G3=G3*(1+SGN(L(N)))
1300 IF G0<G6(N) THEN G0=G6(N)
1320 G2=SQR(7+45*SQR(G0/(1+G8(N)*G8(N))))
1340 G2=G2*M(N)*1.27*(1+SIN(G3)*L(N))
1350 GOTO 1430
1370 T6=T5+12*(1+SGN(T4(N)-T5))*SGN(ABS(T4(N)-T5))
1380 G4=P1*(T6-T4(N))/(24-T1(N))
1390 G0=G7(N)*EXP((T4(N)-T6)/2)
1400 G3=G4+(P1-G4)/4*(1+SGN(L(N)))
1410 G4=0
1420 GOTO 1320
1430 IF G2<J9 THEN J9=G2
1450 Y8=.2:IF T1(N)=0 THEN 1520
1460 IF T1(N)*G4=0 THEN 1520
1470 Y9=C0(N)*SIN(P1*(T6-T3(N))/T1(N))
1480 IF Y9>.174 THEN 1510
1490 Y8=(ATN(SQR(1-Y9*Y9)/Y9)*180/P1-76)^(-.4)
1500 GOTO 1520
1510 Y8=Y9^(.3)
1520 Y9=(3.4+.00544*S0)*Y8*Y6
1530 IF Y9>7 THEN 1560
1540 Y9=.91*Y9-.37
1550 GOTO 1570
1560 Y9=(1.33*Y9-3.31)^(2)/7
1570 IF K9<Y9 THEN K9=Y9
1580 NEXT N
1590 PRINT T5,INT(J9+.5),INT(J9/1.27+.5),INT(K9+.5)
1620 NEXT T5
1640 INPUT "ENTER TO CONTINUE";II$:GOTO 90