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