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

}
