jdo743@uiucuxa.CSO.UIUC.EDU (09/28/86)
To those of you who wanted the MUF source code I talked about in
NET.HAM-RADIO, well here it is.
It is written in Turbo Pascal. I would appreciate any comments concerning
the program.
------------------------------------------------------------------------------
PROGRAM MICROMUF (INPUT,OUTPUT);
{Micromuf - A program to computer the minimum and maximum usable frequencies
of a shortwave path between two specified coordinates. }
CONST PI = 3.14159265;
I = 'INVALID';
VAR a,n,yt,xt,yr,xr,r,x,h,y,u,q,la,ls,hp,sf,fe,se,re,cp,k,l,xz,mf,ff :REAL;
xh,z,fo,yf,ty,tl,yg,zo,yz,yn,yi,xn,sx,xs,wx,lh,lm,ab,rd,d,ex,man :REAL;
transmitter,receiver : string[20];
m,mh,ve,ho,t :integer;
correct :boolean;
key : char;
Function power(man,ex:real):real;
Begin
power:=EXP(ex*LN(man));
End;
Procedure interlat; { Intermediate Latitude & Longitude calculations }
Begin
q:=cos(u*rd)*cos(xt*rd)*sin(k*lm*rd);
x:=q+sin(xt*rd)*cos(k*lm*rd);
xn:=(arctan(x/sqrt(-x*x+1+1E-12)))*d;
q:=(cos(k*lm*rd)-sin(xt*rd)*sin(xn*rd));
yi:=(-arctan(x/sqrt(-x*x+1+1E-12))+(PI/2))*d;
if u < 180.0 then yi:=-yi;
yn:=yt+yi;
if yn > 180.0 then yn:=yn-360.0;
if yn <-180.0 then yn:=yn+360.0;
End;
Procedure mini_f2;
Begin
yz:=yn;
If yn<-160.0 then yz:=yn+360.0;
yg:=(20.0-yz)/50.0;
zo:=20.0*yg/(1+yg+sqr(yg))+5.0*sqr(1-yg/7.0);
z:=xn-zo;
tl:=t-yn/15.0;
if tl > 24.0 then tl:=tl-24.0;
if tl < 0.0 then tl:=tl+24.0;
mh:=m;
If z <= 0.0 then
Begin
z:=-z;
mh:=mh+6;
End;
xh:=cos(30.0*(mh-6.5)*rd); { 1 week delay on equinoxes }
sx:=(abs(xh)+xh)/2.0; { F-layer local summer variance}
wx:=(abs(xh)-xh)/2.0; { F-layer local winter variance}
If z > 77.5 then z:= 77.5;
ty:=tl;
If ty < 5.0 then ty:=tl+24.0;
yf:=(ty-14.0-sx*2.0+wx*2.0-r/175.0)*(7.0-sx*3.0+wx*4.0-r/(150.0-wx*75.0));
If abs(yf) > 60.0 then yf:=60.0;
x:=(1+r/(175.0+sx*175.0));
fo:=6.5*x*cos(yf*rd)*sqrt(cos((z-sx*5.0+wx*5.0)*rd));
ex:=-0.5;
sf:=power((1.0-sqr(cos(a*rd)*6367.0/(6367.0+h))),ex);
ff:=fo*sf;
End;
Procedure e_layer;
Begin
q:=sin(xn*rd)*sin(xs*rd);
x:=q+cos(xn*rd)*cos(xs*rd)*cos((yn-15.0*(t-12.0))*rd);
xz:=(-arctan(x/sqrt(-x*x+1+1E-12))+PI/2)*d;
If xz <= 85.0 then
Begin
ex:=(1.0/3.0);
fe:=3.4*(1.0+0.0016*r)*power(cos(xz*rd),ex);
End
Else
Begin
ex:=-0.5;
fe:=3.4*(1.0+0.0016*r)*power((xz-80.0),ex);
End;
se:=power(1.0-0.965*sqr(cos(a*rd)),ex);
ls:=0.028*sqr(fe)*se;
End;
Begin { Main Program }
rd:=PI/180;
d:=180/PI;
correct:=FALSE;
ClrScr;
Writeln (' *** MICROMUF *** ');
Writeln;
Writeln (' This program calculate the :');
Writeln (' * M. U. F. (Maximum Usable Frequency)');
Writeln;
Writeln (' * L. U. F. (Lowest Usable Frequency)');
Writeln;
Writeln (' of any shortwave sky-wave path.');
Writeln;
Writeln (' Calculations can be done for any month and sunspot number.');
Writeln;
Writeln;
Writeln ('Name transmitter location');
Readln (transmitter);
Writeln;
Repeat;
Writeln ('Transmitter longitude in degrees. (W=+, E=-)');
Readln (yt);
If (yt >=-180.0) and (yt <= 180.0) Then correct:=TRUE
Else Writeln(I);
Until correct = TRUE;
correct:=FALSE;
Writeln;
Repeat;
Writeln ('Transmitter lattitude in degrees. (N=+, S=-)');
Readln(xt);
If (xt >= -90.0) and (xt <= 90.0) Then correct:=TRUE
Else Writeln(I);
Until correct = TRUE;
correct:=FALSE;
Writeln;
Writeln ('Name receiver location.');
Readln(receiver);
Writeln;
Repeat;
Writeln ('Receiver longitude in degrees. (W=+, E=-)');
Readln(yr);
If (yr >= -180.0) and (yr <=180.0) Then correct:=TRUE
Else Writeln(I);
Until correct = TRUE;
correct:=FALSE;
Writeln;
Repeat;
Writeln ('Receiver lattitude in degrees. (N=+, S=-)');
readln(xr);
If (xr >=-90.0) and (xr<=90.0) Then correct:=TRUE
Else Writeln(I);
Until correct = TRUE;
correct:=FALSE;
Writeln;
Repeat;
Writeln ('Sunspot number.');
Readln (r);
If (r >= 1.0) and (r <=180.0) Then correct:=TRUE
Else Writeln(I);
Until correct = TRUE;
correct:=FALSE;
Writeln;
Repeat;
Writeln ('Month.');
Readln (m);
If (m >= 1) and (m <= 12) Then correct:=TRUE
Else Writeln(I);
Until correct = TRUE;
{ Geometry Calculations }
q:=sin(xt*rd)*sin(xr*rd);
x:=q+cos(xt*rd)*cos(xr*rd)*cos(yt*rd-yr*rd);
la:=(-arctan(x/sqrt(-x*x+1+1E-12))+(PI/2))*d;
l:=111.1*la;
q:=(sin(xr*rd)-sin(xt*rd)*cos(la*rd));
x:=q/cos(xt*rd)/sin(la*rd);
u:=(-arctan(x/sqrt(-x*x+1+1E-12))+(PI/2))*d;
If yt-yr <= 0 Then u:=360-u;
h:=275+r/2;
xs:=23.4*cos(30*(m-6.25)*rd);
n:=n+1;
lh:=l/n;
While lh > 4500.0 Do
Begin
n:=n+1;
lh:=l/n;
End; { While }
lm:=la/n;
a:=(arctan((cos(0.5*lm*rd)-6367.0/(h+6367.0))/sin(0.5*lm*rd)))*d;
While a < 1.5 Do
Begin
n:=n+1;
lh:=lh/n;
While lh > 4500.0 Do
Begin
n:=n+1;
lh:=l/n;
End; { While }
lm:=la/n;
a:=(arctan((cos(0.5*lm*rd)-6367.0/(h+6367.0))/sin(0.5*lm*rd)))*d;
End;
{ Plot chart on screen }
ClrScr;
Writeln ('From: ',transmitter,' to: ',receiver);
Write ('Month: ',m);
Writeln (' SSN: ',r:3:0,' Dist: ',Round(l+0.5),' KM');
Writeln ('Azim: ',Round(u+0.5),' degrees. F-hops: ',n:4:2);
ve:=4;
ho:=1;
GotoXY(ho,ve);
q:=34.0;
While q >=2.0 Do
Begin
Writeln ('I I',q:2:0); { 25 spaces }
q:=q-2.0;
End; { While }
Writeln ('---------------------------'); { 27 dashes }
Writeln (' 0 2 4 6 8 10 14 18 22 H (UTC)');
Writeln (' +: MUF -: LUF');
ve:=4;
ho:=32;
GotoXY(ho,ve);
Writeln ('mHz');
For t:=1 to 24 Do
Begin
ab:=0.0;
k:=0.5;
interlat;
mini_f2;
mf:=ff;
k:=n-0.5;
interlat;
mini_f2;
If ff < mf then mf := ff;
ve:=21-Round(mf/2.0+0.5);
ho:=t+1;
if ve < 4 then ve:=4;
GotoXY(ho,ve);
Writeln ('+');
While k <= n-0.25 Do
Begin
interlat;
e_layer;
ab:=ab+ls;
k:=k+0.5;
End;
ve:=20-Round(ab+0.5);
If ve < 4 Then ve:=4;
If ve > 20 Then ve:=20;
GotoXY(ho,ve);
Writeln('-');
End;
Writeln;
Writeln;
Writeln;
Writeln;
Writeln;
Write (' ----- Press a key to continue ----- ');
Readln (key);
End.
{ This program uses 'MINI-F2' devised by R. Fricker (BBC external
services) for FO-F calculations and L.M. Muggleton's formula for
FO-E calculations.
For the L.U.F. a minimum useable fieldstrength of 30 DBUV at the
receiver and 250 KW of transmitter power (aerial gain: 18 DBI) are
assumed. The L.U.F. is derived from absorption calculations based
on the work of Piggot, George, Samuel, and Bradley. In spite of the
program's simplicity it gives a good impression of the ionosphere's
behaviour and can be used for propagation predictions.
Hans Bakhuizen
Propagation Unit; Frequency Bureau
Radio Netherlands
P.O. Box 222
1200 JG Hilversum Holland
(C) Copyright Media Network June 1984
Translation from Basic into TURBO Pascal by Jonathan D Ogden on
September 26, 1986.
Jonathan D Ogden
402 e Daniel
Champaign, Il 61820 USA
}