perl/Minimuf.pm

   1 #!/usr/bin/perl -w
   2 # A perl Minimuf calculator, nicked from the minimuf program written in
   3 # C.
   4 #
   5 # Translated and modified for my own purposes by Dirk Koopman G1TLH
   6 #
   7 # Copyright (c) 1999 Dirk Koopman G1TLH
   8 #
   9 # The original copyright:-
  10 #/***********************************************************************
  11 # *                                                                     *
  12 # * Copyright (c) David L. Mills 1994-1998                              *
  13 # *                                                                     *
  14 # * Permission to use, copy, modify, and distribute this software and   *
  15 # * its documentation for any purpose and without fee is hereby         *
  16 # * granted, provided that the above copyright notice appears in all    *
  17 # * copies and that both the copyright notice and this permission       *
  18 # * notice appear in supporting documentation, and that the name        *
  19 # * University of Delaware not be used in advertising or publicity      *
  20 # * pertaining to distribution of the software without specific,        *
  21 # * written prior permission.  The University of Delaware makes no      *
  22 # * representations about the suitability this software for any         *
  23 # * purpose. It is provided "as is" without express or implied          *
  24 # * warranty.                                                           *
  25 # *                                                                     *
  26 # ***********************************************************************
  27 #
  28 # MINIMUF 3.5 from QST December 1982
  29 # (originally in BASIC)
  30 #
  31 # $Id$
  32 #
  33 #
  34
  35 package Minimuf;
  36
  37 use strict;
  38 use POSIX;
  39 use vars qw($pi $d2r $r2d $halfpi $pi2 $VOFL $R $hE $hF $GAMMA $LN10
  40                     $MINBETA $BOLTZ $NTEMP $DELTAF $MPATH $GLOSS $SLOSS
  41             $noise);
  42
  43 $pi = 3.141592653589;
  44 $d2r = ($pi/180);
  45 $r2d = (180/$pi);
  46 $halfpi = $pi/2;
  47 $pi2 = $pi*2;
  48 $VOFL = 2.9979250e8;                    # velocity of light
  49 $R = 6371.2;                                    # radius of the Earth (km)
  50 $hE = 110;                                              # mean height of E layer (km)
  51 $hF = 320;                                              # mean height of F layer (km)
  52 $GAMMA = 1.42;                                  # geomagnetic constant
  53 $LN10 = 2.302585;                               # natural logarithm of 10
  54 $MINBETA = (10 * $d2r);                 # min elevation angle (rad)
  55 $BOLTZ = 1.380622e-23;                  # Boltzmann's constant
  56 $NTEMP = 290;                                   # receiver noise temperature (K)
  57 $DELTAF = 2500;                                 # communication bandwidth (Hz)
  58 $MPATH = 3;                                             # multipath threshold (dB)
  59 $GLOSS = 3;                                             # ground-reflection loss (dB)
  60 $SLOSS = 10;                                    # excess system loss
  61 $noise = 10 * log10($BOLTZ * $NTEMP * $DELTAF) + 30;
  62
  63 # basic SGN function
  64 sub SGN
  65 {
  66         my $x = shift;
  67         return 0 if $x == 0;
  68         return ($x > 0) ? 1 : -1;
  69 }
  70
  71 #
  72 # MINIMUF 3.5 (From QST December 1982, originally in BASIC)
  73 #
  74
  75 sub minimuf
  76 {
  77         my $flux = shift;               # 10-cm solar flux
  78         my $month = shift;              # month of year (1 - 12)
  79         my $day = shift;                # day of month (1 - 31)
  80         my $hour = shift;               # hour of day (utc) (0 - 23)
  81         my $lat1 = shift;               # transmitter latitude (deg n)
  82         my $lon1 = shift;               # transmitter longitude (deg w)
  83         my $lat2 = shift;               # receiver latitude (deg n)
  84         my $lon2 = shift;               # receiver longitude (deg w)
  85
  86         my $ssn;                # sunspot number dervived from flux
  87         my $muf;                # maximum usable frequency
  88         my $dist;               # path angle (rad)
  89         my ($a, $p, $q);                # unfathomable local variables
  90         my ($y1, $y2, $y3);
  91         my ($t, $t4, $t9);
  92         my ($g0, $g8);
  93         my ($k1, $k6, $k8, $k9);
  94         my ($m9, $c0);
  95         my ($ftemp, $gtemp);    # volatile temps
  96
  97         # Determine geometry and invariant coefficients
  98         $ssn = spots($flux);
  99         $ftemp = sin($lat1) * sin($lat2) + cos($lat1) * cos($lat2) *
 100             cos($lon2 - $lon1);
 101         $ftemp = -1 if ($ftemp < -1);
 102         $ftemp = 1 if ($ftemp > 1);
 103         $dist = acos($ftemp);
 104         $k6 = 1.59 * $dist;
 105         $k6 = 1 if ($k6 < 1);
 106         $p = sin($lat2);
 107         $q = cos($lat2);
 108         $a = (sin($lat1) - $p * cos($dist)) / ($q * sin($dist));
 109         $y1 = 0.0172 * (10 + ($month - 1) * 30.4 + $day);
 110         $y2 = 0.409 * cos($y1);
 111         $ftemp = 2.5 * $dist / $k6;
 112         $ftemp = $halfpi if ($ftemp > $halfpi);
 113         $ftemp = sin($ftemp);
 114         $m9 = 1 + 2.5 * $ftemp * sqrt($ftemp);
 115         $muf = 100;
 116
 117         # Loop along path
 118         for ($k1 = 1 / (2 * $k6); $k1 <= 1 - 1 / (2 * $k6); $k1 += abs(0.9999 - 1 / $k6)) {
 119                 $gtemp = $dist * $k1;
 120                 $ftemp = $p * cos($gtemp) + $q * sin($gtemp) * $a;
 121                 $ftemp = -1 if ($ftemp < -1);
 122                 $ftemp = 1 if ($ftemp > 1);
 123                 $y3 = $halfpi - acos($ftemp);
 124                 $ftemp = (cos($gtemp) - $ftemp * $p) / ($q * sqrt(1 - $ftemp * $ftemp));
 125                 $ftemp = -1 if ($ftemp < -1);
 126                 $ftemp = 1 if ($ftemp > 1);
 127                 $ftemp = $lon2 + SGN(sin($lon1 - $lon2)) * acos($ftemp);
 128                 $ftemp += $pi2 if ($ftemp < 0);
 129                 $ftemp -= $pi2 if ($ftemp >= $pi2);
 130                 $ftemp = 3.82 * $ftemp + 12 + 0.13 * (sin($y1) + 1.2 * sin(2 * $y1));
 131                 $k8 = $ftemp - 12 * (1 + SGN($ftemp - 24)) * SGN(abs($ftemp - 24));
 132                 if (cos($y3 + $y2) <= -0.26) {
 133                         $k9 = 0;
 134                         $g0 = 0;
 135                 } else {
 136                         $ftemp = (-0.26 + sin($y2) * sin($y3)) / (cos($y2) * cos($y3) + 0.001);
 137                         $k9 = 12 - atan($ftemp / sqrt(abs(1 - $ftemp * $ftemp))) * 7.639437;
 138                         $t = $k8 - $k9 / 2 + 12 * (1 - SGN($k8 - $k9 / 2)) * SGN(abs($k8 - $k9 / 2));
 139                         $t4 = $k8 + $k9 / 2 - 12 * (1 + SGN($k8 + $k9 / 2 - 24)) * SGN(abs($k8 + $k9 / 2 - 24));
 140                         $c0 = abs(cos($y3 + $y2));
 141                         $t9 = 9.7 * pow($c0, 9.6);
 142                         $t9 = 0.1 if ($t9 < 0.1);
 143
 144                         $g8 = $pi * $t9 / $k9;
 145                         if (($t4 < $t && ($hour - $t4) * ($t - $hour) > 0.) || ($t4 >= $t && ($hour - $t) * ($t4 - $hour) <= 0)) {
 146                                 $ftemp = $hour + 12 * (1 + SGN($t4 - $hour)) * SGN(abs($t4 - $hour));
 147                                 $ftemp = ($t4 - $ftemp) / 2;
 148                                 $g0 = $c0 * ($g8 * (exp(-$k9 / $t9) + 1)) * exp($ftemp) / (1 + $g8 * $g8);
 149                         } else {
 150                                 $ftemp = $hour + 12 * (1 + SGN($t - $hour)) * SGN(abs($t - $hour));
 151                                 $gtemp = $pi * ($ftemp - $t) / $k9;
 152                                 $ftemp = ($t - $ftemp) / $t9;
 153                                 $g0 = $c0 * (sin($gtemp) + $g8 * (exp($ftemp) - cos($gtemp))) / (1 + $g8 * $g8);
 154                                 $ftemp = $c0 * ($g8 * (exp(-$k9 / $t9) + 1)) * exp(($k9 - 24) / 2) / (1 + $g8 * $g8);
 155                                 $g0 = $ftemp if ($g0 < $ftemp);
 156                         }
 157                 }
 158                 $ftemp = (1 + $ssn / 250) * $m9 * sqrt(6 + 58 * sqrt($g0));
 159                 $ftemp *= 1 - 0.1 * exp(($k9 - 24) / 3);
 160                 $ftemp *= 1 + 0.1 * (1 - SGN($lat1) * SGN($lat2));
 161                 $ftemp *= 1 - 0.1 * (1 + SGN(abs(sin($y3)) - cos($y3)));
 162                 $muf = $ftemp if ($ftemp < $muf);
 163         }
 164         return $muf;
 165 }
 166
 167 #
 168 # spots(flux) - Routine to map solar flux to sunspot number.
 169 #
 170 # THis routine was done by eyeball and graph on p. 22-6 of the 1991
 171 # ARRL Handbook. The nice curve fitting was done using Mathematica.
 172 #
 173 sub spots
 174 {
 175         my $flux = shift; # 10-cm solar flux
 176         my $ftemp;                      # double temp
 177
 178         return 0 if ($flux < 65);
 179         if ($flux < 110) {
 180                 $ftemp = $flux - 200.6;
 181                 $ftemp = 108.36 - .005896 * $ftemp * $ftemp;
 182         } elsif ($flux < 213) {
 183                 $ftemp = 60 + 1.0680 * ($flux - 110);
 184         } else {
 185                 $ftemp = $flux - 652.9;
 186                 $ftemp = 384.0 - 0.0011059 * $ftemp * $ftemp;
 187         }
 188         return $ftemp;
 189 }
 190
 191 # ion - determine paratmeters for hop h
 192 #
 193 # This routine determines the reflection zones for each hop along the
 194 # path and computes the minimum F-layer MUF, maximum E-layer MUF,
 195 # ionospheric absorption factor and day/night flags for the entire
 196 # path.
 197
 198 sub ion
 199 {
 200         my $h = shift;                          # hop index
 201         my $d = shift;                          # path angle (rad)
 202         my $fcF = shift;                        # F-layer critical frequency
 203         my $ssn = shift;            # current sunspot number
 204     my $daynight = shift;               # ref to daynight array one per hop
 205
 206         # various refs to arrays
 207         my $mufE = shift;
 208         my $mufF = shift;
 209         my $absorp = shift;
 210
 211         my $beta;               # elevation angle (rad)
 212         my $psi;                # sun zenith angle (rad)
 213         my $dhop;               # hop angle / 2 (rad)
 214         my $dist;               # path angle (rad)
 215         my $phiF;               # F-layer angle of incidence (rad)
 216         my $phiE;               # E-layer angle of incidence (rad)
 217         my $fcE;                # E-layer critical frequency (MHz)
 218         my $ftemp;              # double temp
 219
 220
 221         # Determine the path geometry, E-layer angle of incidence and
 222         # minimum F-layer MUF. The F-layer MUF is determined from the
 223         # F-layer critical frequency previously calculated by MINIMUF
 224         # 3.5 and the secant law and so depends only on the F-layer
 225         # angle of incidence. This is somewhat of a crock; however,
 226         # doing it with MINIMUF 3.5 on a hop-by-hop basis results in
 227         # rather serious errors.
 228
 229
 230         $dhop = $d / ($h * 2);
 231         $beta = atan((cos($dhop) - $R / ($R + $hF)) / sin($dhop));
 232         $ftemp = $R * cos($beta) / ($R + $hE);
 233         $phiE = atan($ftemp / sqrt(1 - $ftemp * $ftemp));
 234         $ftemp = $R * cos($beta) / ($R + $hF);
 235         $phiF = atan($ftemp / sqrt(1 - $ftemp * $ftemp));
 236         $$mufF->[$h] = $fcF / cos($phiF);;
 237         for ($dist = $dhop; $dist < $d; $dist += $dhop * 2) {
 238
 239                 # Calculate the E-layer critical frequency and MUF.
 240
 241                 $fcE = 0;
 242                 $psi = zenith($dist);
 243                 $ftemp = cos($psi);
 244                 $fcE = .9 * pow((180. + 1.44 * $ssn) * $ftemp, .25) if ($ftemp > 0);
 245                 $fcE = .005 * $ssn if ($fcE < .005 * $ssn);
 246                 $ftemp = $fcE / cos($phiE);
 247                 $mufE->[$h] = $ftemp if ($ftemp > $mufE->[$h]);
 248
 249                 # Calculate ionospheric absorption coefficient and
 250                 # day/night indicators. Note that some hops along a
 251                 # path can be in daytime and others in nighttime.
 252
 253                 $ftemp = $psi;
 254                 if ($ftemp > 100.8 * $d2r) {
 255                         $ftemp = 100.8 * $d2r;
 256                         $daynight->[$h] |= 2;
 257                 } else {
 258                         $daynight->[$h] |= 1;
 259                 }
 260                 $ftemp = cos(90. / 100.8 * $ftemp);
 261                 $ftemp = 0. if ($ftemp < 0.);
 262                 $ftemp = (1. + .0037 * $ssn) * pow($ftemp, 1.3);
 263                 $ftemp = .1 if ($ftemp < .1);
 264                 $absorp->[$h] += $ftemp;
 265         }
 266 }
 267
 268
 269 #
 270 # pathloss(freq, hop) - Compute receive power for given path.
 271 #
 272 # This routine determines which of the three ray paths determined
 273 # previously are usable. It returns the hop index of the best of these
 274 # or zero if none are found.
 275
 276 sub pathloss
 277 {
 278         my $hop = shift;                        # minimum hops
 279         my $freq = shift;                       # frequency
 280     my $txpower = shift || 20;  # transmit power
 281     my $rsens = shift || -123;  # receiver sensitivity
 282     my $antgain = shift || 0;   # antenna gain
 283
 284     my $daynight = shift;               # ref to daynight array one per hop
 285     my $beta = shift;
 286         my $path = shift;
 287         my $mufF = shift;
 288         my $mufE = shift;
 289     my $absorp = shift;
 290     my $dB2 = shift;
 291
 292         my $h;                                          # hop number
 293         my $level;                                      # max signal (dBm)
 294         my $signal;                                     # receive signal (dBm)
 295         my $ftemp;                                      # double temp
 296         my $j;                                          # index temp
 297
 298         #
 299         # Calculate signal and noise for all hops. The noise level is
 300         # -140 dBm for a receiver bandwidth of 2500 Hz and noise
 301         # temperature 290 K. The receiver sensitivity is assumed -123
 302         # dBm (0.15 V at 50 Ohm for 10 dB S/N). Paths where the signal
 303         # is less than the noise or when the frequency exceeds the F-
 304         # layer MUF are considered unusable.
 305
 306         $level = $noise;
 307         $j = 0;
 308         for ($h = $hop; $h < $hop + 3; $h++) {
 309 #               $daynight->[$h] &= ~(P_E | P_S | P_M);
 310                 if ($freq < 0.85 * $mufF->[$h]) {
 311
 312                         # Transmit power (dBm)
 313
 314                         $signal = $txpower + $antgain + 30;
 315
 316                         # Path loss
 317
 318                         $signal -= 32.44 + 20 * log10($path->[$h] * $freq) + $SLOSS;
 319
 320                         # Ionospheric loss
 321
 322                         $ftemp = $R * cos($beta->[$h]) / ($R + $hE);
 323                         $ftemp = atan($ftemp / sqrt(1 - $ftemp * $ftemp));
 324                         $signal -= 677.2 * $absorp->[$h] / cos($ftemp) / (pow(($freq + $GAMMA), 1.98) + 10.2);
 325
 326                         # Ground reflection loss
 327
 328                         $signal -= $h * $GLOSS;
 329                         $dB2->[$h] = $signal;
 330
 331                         # Paths where the signal is greater than the
 332                         # noise, but less than the receiver sensitivity
 333                         # are marked 's'. Paths below the E-layer MUF
 334                         # are marked 'e'. When comparing for maximum
 335                         # signal, The signal for these paths is reduced
 336                         # by 3 dB so they will be used only as a last
 337                         # resort.
 338
 339
 340                         $daynight->[$h] |= 4 if ($signal < $rsens);
 341                         if ($freq < $mufE->[$h]) {
 342                                 $daynight->[$h] |= 8;
 343                                 $signal -= $MPATH;
 344                         }
 345                         if ($signal > $level) {
 346                                 $level = $signal;
 347                                 $j = $h;
 348                         }
 349                 }
 350         }
 351
 352         # We have found the best path. If this path is less than 3 dB
 353         # above the RMS sum of the other paths, the path is marked 'm'.
 354
 355         return 0 if ($j == 0);
 356
 357         $ftemp = 0;
 358         for ($h = $hop; $h < $hop + 3; $h++) {
 359                 $ftemp += exp(2 / 10 * $dB2->[$h] * $LN10) if ($h != $j);
 360         }
 361         $ftemp = 10 / 2 * log10($ftemp);
 362         $daynight->[$j] |= 16 if ($level < $ftemp + $MPATH);
 363
 364         return $j;
 365 }
 366
 367 # zenith(dist) - Determine sun zenith angle at reflection zone.
 368
 369 sub zenith
 370 {
 371         my $dist = shift;                       # path angle
 372         my $txlat = shift;          # tx latitude (rad)
 373         my $txlong = shift;         # tx longitude (rad)
 374         my $txbearing = shift;      # tx bearing
 375         my $pathangle = shift;      # 'b1'
 376         my $lats = shift;           # subsolar latitude
 377         my $lons = shift;           # subsolar longitude
 378
 379         my ($latr, $lonr);                      # reflection zone coordinates (rad)
 380         my $thetar;                                     # reflection zone angle (rad)
 381         my $psi;                                        # sun zenith angle (rad)
 382
 383         # Calculate reflection zone coordinates.
 384
 385         $latr = acos(cos($dist) * sin($txlat) + sin($dist) * cos($txlat) * cos($txbearing));
 386         $latr += $pi if ($latr < 0);
 387         $latr = $halfpi - $latr;
 388         $lonr = acos((cos($dist) - sin($latr) * sin($txlat)) / (cos($latr) * cos($txlat)));
 389         $lonr += $pi if ($lonr < 0);
 390         $lonr = - $lonr if ($pathangle < 0);
 391         $lonr = $txlong - $lonr;
 392         $lonr -= $pi2 if ($lonr >= $pi);
 393         $lonr += $pi2 if ($lonr <= -$pi);
 394         $thetar = $lons - $lonr;
 395         $thetar = $pi2 - $thetar if ($thetar > $pi);
 396         $thetar -= $pi2 if ($thetar < - $pi);
 397
 398         # Calculate sun zenith angle.
 399
 400         $psi = acos(sin($latr) * sin($lats) + cos($latr) * cos($lats) * cos($thetar));
 401         $psi += $pi if ($psi < 0);
 402         return($psi);
 403 }
 404
 405
 406
 407 1;