]> dxcluster.org Git - spider.git/commitdiff
started the process of adding minimuf, translated the minimuf subroutines
authordjk <djk>
Tue, 19 Oct 1999 23:48:19 +0000 (23:48 +0000)
committerdjk <djk>
Tue, 19 Oct 1999 23:48:19 +0000 (23:48 +0000)
and created the Minimuf.pm file

Changes
cmd/show/muf.pl [new file with mode: 0644]
perl/Minimuf.pm [new file with mode: 0644]
perl/proto.html

diff --git a/Changes b/Changes
index 93142c226cfc7b166f7cbcd6f22fb366ddc94e31..06b5c67ffba69e30a72996a5f0b9e53203b59d88 100644 (file)
--- a/Changes
+++ b/Changes
@@ -1,3 +1,5 @@
+20Oct99=======================================================================
+1. Translated all the subroutines of minimuf into perl as Minimuf.pm
 18Oct99=======================================================================
 1. changed help command so that it works correctly with multiple title lines.
 2. added to address to the list of things a message checks to see whether it
diff --git a/cmd/show/muf.pl b/cmd/show/muf.pl
new file mode 100644 (file)
index 0000000..e69de29
diff --git a/perl/Minimuf.pm b/perl/Minimuf.pm
new file mode 100644 (file)
index 0000000..e493e5b
--- /dev/null
@@ -0,0 +1,407 @@
+#!/usr/bin/perl -w
+# A perl Minimuf calculator, nicked from the minimuf program written in
+# C.
+#
+# Translated and modified for my own purposes by Dirk Koopman G1TLH
+#
+# Copyright (c) 1999 Dirk Koopman G1TLH
+#
+# The original copyright:-
+#/***********************************************************************
+# *                                                                     *
+# * Copyright (c) David L. Mills 1994-1998                              *
+# *                                                                     *
+# * Permission to use, copy, modify, and distribute this software and   *
+# * its documentation for any purpose and without fee is hereby         *
+# * granted, provided that the above copyright notice appears in all    *
+# * copies and that both the copyright notice and this permission       *
+# * notice appear in supporting documentation, and that the name        *
+# * University of Delaware not be used in advertising or publicity      *
+# * pertaining to distribution of the software without specific,        *
+# * written prior permission.  The University of Delaware makes no      *
+# * representations about the suitability this software for any         *
+# * purpose. It is provided "as is" without express or implied          *
+# * warranty.                                                           *
+# *                                                                     *
+# ***********************************************************************
+#
+# MINIMUF 3.5 from QST December 1982
+# (originally in BASIC)
+#
+# $Id$
+#
+#
+
+package Minimuf;
+
+use strict;
+use POSIX;
+use vars qw($pi $d2r $r2d $halfpi $pi2 $VOFL $R $hE $hF $GAMMA $LN10
+                   $MINBETA $BOLTZ $NTEMP $DELTAF $MPATH $GLOSS $SLOSS
+            $noise);
+
+$pi = 3.141592653589;
+$d2r = ($pi/180);
+$r2d = (180/$pi);
+$halfpi = $pi/2;
+$pi2 = $pi*2;
+$VOFL = 2.9979250e8;                   # velocity of light
+$R = 6371.2;                                   # radius of the Earth (km)  
+$hE = 110;                                             # mean height of E layer (km) 
+$hF = 320;                                             # mean height of F layer (km) 
+$GAMMA = 1.42;                                 # geomagnetic constant 
+$LN10 = 2.302585;                              # natural logarithm of 10 
+$MINBETA = (10 * $d2r);                        # min elevation angle (rad) 
+$BOLTZ = 1.380622e-23;                 # Boltzmann's constant 
+$NTEMP = 290;                                  # receiver noise temperature (K) 
+$DELTAF = 2500;                                        # communication bandwidth (Hz) 
+$MPATH = 3;                                            # multipath threshold (dB) 
+$GLOSS = 3;                                            # ground-reflection loss (dB) 
+$SLOSS = 10;                                   # excess system loss 
+$noise = 10 * log10($BOLTZ * $NTEMP * $DELTAF) + 30;
+
+# basic SGN function
+sub SGN
+{
+       my $x = shift;
+       return 0 if $x == 0;
+       return ($x > 0) ? 1 : -1;
+}
+
+#
+# MINIMUF 3.5 (From QST December 1982, originally in BASIC)
+#
+
+sub minimuf
+{
+       my $flux = shift;               # 10-cm solar flux 
+       my $month = shift;              # month of year (1 - 12) 
+       my $day = shift;                # day of month (1 - 31) 
+       my $hour = shift;               # hour of day (utc) (0 - 23) 
+       my $lat1 = shift;               # transmitter latitude (deg n) 
+       my $lon1 = shift;               # transmitter longitude (deg w) 
+       my $lat2 = shift;               # receiver latitude (deg n) 
+       my $lon2 = shift;               # receiver longitude (deg w) 
+       
+       my $ssn;                # sunspot number dervived from flux 
+       my $muf;                # maximum usable frequency 
+       my $dist;               # path angle (rad) 
+       my ($a, $p, $q);                # unfathomable local variables 
+       my ($y1, $y2, $y3);
+       my ($t, $t4, $t9);
+       my ($g0, $g8);
+       my ($k1, $k6, $k8, $k9);
+       my ($m9, $c0);
+       my ($ftemp, $gtemp);    # volatile temps 
+       
+       # Determine geometry and invariant coefficients
+       $ssn = spots($flux);
+       $ftemp = sin($lat1) * sin($lat2) + cos($lat1) * cos($lat2) *
+           cos($lon2 - $lon1);
+       $ftemp = -1 if ($ftemp < -1);
+       $ftemp = 1 if ($ftemp > 1);
+       $dist = acos($ftemp);
+       $k6 = 1.59 * $dist;
+       $k6 = 1 if ($k6 < 1);
+       $p = sin($lat2);
+       $q = cos($lat2);
+       $a = (sin($lat1) - $p * cos($dist)) / ($q * sin($dist));
+       $y1 = 0.0172 * (10 + ($month - 1) * 30.4 + $day);
+       $y2 = 0.409 * cos($y1);
+       $ftemp = 2.5 * $dist / $k6;
+       $ftemp = $halfpi if ($ftemp > $halfpi);
+       $ftemp = sin($ftemp);
+       $m9 = 1 + 2.5 * $ftemp * sqrt($ftemp);
+       $muf = 100;
+
+       # Loop along path
+       for ($k1 = 1 / (2 * $k6); $k1 <= 1 - 1 / (2 * $k6); $k1 += abs(0.9999 - 1 / $k6)) {
+               $gtemp = $dist * $k1;
+               $ftemp = $p * cos($gtemp) + $q * sin($gtemp) * $a;
+               $ftemp = -1 if ($ftemp < -1);
+               $ftemp = 1 if ($ftemp > 1);
+               $y3 = $halfpi - acos($ftemp);
+               $ftemp = (cos($gtemp) - $ftemp * $p) / ($q * sqrt(1 - $ftemp * $ftemp));
+               $ftemp = -1 if ($ftemp < -1);
+               $ftemp = 1 if ($ftemp > 1);
+               $ftemp = $lon2 + SGN(sin($lon1 - $lon2)) * acos($ftemp);
+               $ftemp += $pi2 if ($ftemp < 0);
+               $ftemp -= $pi2 if ($ftemp >= $pi2);
+               $ftemp = 3.82 * $ftemp + 12 + 0.13 * (sin($y1) + 1.2 * sin(2 * $y1));
+               $k8 = $ftemp - 12 * (1 + SGN($ftemp - 24)) * SGN(abs($ftemp - 24));
+               if (cos($y3 + $y2) <= -0.26) {
+                       $k9 = 0;
+                       $g0 = 0;
+               } else {
+                       $ftemp = (-0.26 + sin($y2) * sin($y3)) / (cos($y2) * cos($y3) + 0.001);
+                       $k9 = 12 - atan($ftemp / sqrt(abs(1 - $ftemp * $ftemp))) * 7.639437;
+                       $t = $k8 - $k9 / 2 + 12 * (1 - SGN($k8 - $k9 / 2)) * SGN(abs($k8 - $k9 / 2));
+                       $t4 = $k8 + $k9 / 2 - 12 * (1 + SGN($k8 + $k9 / 2 - 24)) * SGN(abs($k8 + $k9 / 2 - 24));
+                       $c0 = abs(cos($y3 + $y2));
+                       $t9 = 9.7 * pow($c0, 9.6);
+                       $t9 = 0.1 if ($t9 < 0.1);
+                       
+                       $g8 = $pi * $t9 / $k9;
+                       if (($t4 < $t && ($hour - $t4) * ($t - $hour) > 0.) || ($t4 >= $t && ($hour - $t) * ($t4 - $hour) <= 0)) {
+                               $ftemp = $hour + 12 * (1 + SGN($t4 - $hour)) * SGN(abs($t4 - $hour));
+                               $ftemp = ($t4 - $ftemp) / 2;
+                               $g0 = $c0 * ($g8 * (exp(-$k9 / $t9) + 1)) * exp($ftemp) / (1 + $g8 * $g8);
+                       } else {
+                               $ftemp = $hour + 12 * (1 + SGN($t - $hour)) * SGN(abs($t - $hour));
+                               $gtemp = $pi * ($ftemp - $t) / $k9;
+                               $ftemp = ($t - $ftemp) / $t9;
+                               $g0 = $c0 * (sin($gtemp) + $g8 * (exp($ftemp) - cos($gtemp))) / (1 + $g8 * $g8);
+                               $ftemp = $c0 * ($g8 * (exp(-$k9 / $t9) + 1)) * exp(($k9 - 24) / 2) / (1 + $g8 * $g8);
+                               $g0 = $ftemp if ($g0 < $ftemp);
+                       }
+               }
+               $ftemp = (1 + $ssn / 250) * $m9 * sqrt(6 + 58 * sqrt($g0));
+               $ftemp *= 1 - 0.1 * exp(($k9 - 24) / 3);
+               $ftemp *= 1 + 0.1 * (1 - SGN($lat1) * SGN($lat2));
+               $ftemp *= 1 - 0.1 * (1 + SGN(abs(sin($y3)) - cos($y3)));
+               $muf = $ftemp if ($ftemp < $muf);
+       }
+       return $muf;
+}
+
+#
+# spots(flux) - Routine to map solar flux to sunspot number.
+#
+# THis routine was done by eyeball and graph on p. 22-6 of the 1991
+# ARRL Handbook. The nice curve fitting was done using Mathematica.
+# 
+sub spots
+{
+       my $flux = shift; # 10-cm solar flux 
+       my $ftemp;                      # double temp 
+
+       return 0 if ($flux < 65);
+       if ($flux < 110) {
+               $ftemp = $flux - 200.6;
+               $ftemp = 108.36 - .005896 * $ftemp * $ftemp;
+       } elsif ($flux < 213) {
+               $ftemp = 60 + 1.0680 * ($flux - 110);
+       } else {
+               $ftemp = $flux - 652.9;
+               $ftemp = 384.0 - 0.0011059 * $ftemp * $ftemp;
+       }
+       return $ftemp;
+}
+
+# ion - determine paratmeters for hop h
+#
+# This routine determines the reflection zones for each hop along the
+# path and computes the minimum F-layer MUF, maximum E-layer MUF,
+# ionospheric absorption factor and day/night flags for the entire
+# path.
+
+sub ion
+{
+       my $h = shift;                          # hop index
+       my $d = shift;                          # path angle (rad)
+       my $fcF = shift;                        # F-layer critical frequency 
+       my $ssn = shift;            # current sunspot number
+    my $daynight = shift;              # ref to daynight array one per hop
+       
+       # various refs to arrays
+       my $mufE = shift;
+       my $mufF = shift;
+       my $absorp = shift;
+       
+       my $beta;               # elevation angle (rad) 
+       my $psi;                # sun zenith angle (rad) 
+       my $dhop;               # hop angle / 2 (rad) 
+       my $dist;               # path angle (rad) 
+       my $phiF;               # F-layer angle of incidence (rad) 
+       my $phiE;               # E-layer angle of incidence (rad) 
+       my $fcE;                # E-layer critical frequency (MHz) 
+       my $ftemp;              # double temp 
+    
+
+       # Determine the path geometry, E-layer angle of incidence and
+       # minimum F-layer MUF. The F-layer MUF is determined from the
+       # F-layer critical frequency previously calculated by MINIMUF
+       # 3.5 and the secant law and so depends only on the F-layer
+       # angle of incidence. This is somewhat of a crock; however,
+       # doing it with MINIMUF 3.5 on a hop-by-hop basis results in
+       # rather serious errors.
+        
+
+       $dhop = $d / ($h * 2);
+       $beta = atan((cos($dhop) - $R / ($R + $hF)) / sin($dhop));
+       $ftemp = $R * cos($beta) / ($R + $hE);
+       $phiE = atan($ftemp / sqrt(1 - $ftemp * $ftemp));
+       $ftemp = $R * cos($beta) / ($R + $hF);
+       $phiF = atan($ftemp / sqrt(1 - $ftemp * $ftemp));
+       $$mufF->[$h] = $fcF / cos($phiF);;
+       for ($dist = $dhop; $dist < $d; $dist += $dhop * 2) {
+
+               # Calculate the E-layer critical frequency and MUF.
+                
+               $fcE = 0;
+               $psi = zenith($dist);
+               $ftemp = cos($psi);
+               $fcE = .9 * pow((180. + 1.44 * $ssn) * $ftemp, .25) if ($ftemp > 0);
+               $fcE = .005 * $ssn if ($fcE < .005 * $ssn);
+               $ftemp = $fcE / cos($phiE);
+               $mufE->[$h] = $ftemp if ($ftemp > $mufE->[$h]);
+
+               # Calculate ionospheric absorption coefficient and
+               # day/night indicators. Note that some hops along a
+               # path can be in daytime and others in nighttime.
+
+               $ftemp = $psi;
+               if ($ftemp > 100.8 * $d2r) {
+                       $ftemp = 100.8 * $d2r;
+                       $daynight->[$h] |= 2;
+               } else {
+                       $daynight->[$h] |= 1;
+               }
+               $ftemp = cos(90. / 100.8 * $ftemp);
+               $ftemp = 0. if ($ftemp < 0.);
+               $ftemp = (1. + .0037 * $ssn) * pow($ftemp, 1.3);
+               $ftemp = .1 if ($ftemp < .1);
+               $absorp->[$h] += $ftemp;
+       }
+}
+
+
+#
+# pathloss(freq, hop) - Compute receive power for given path.
+#
+# This routine determines which of the three ray paths determined
+# previously are usable. It returns the hop index of the best of these
+# or zero if none are found.
+
+sub pathloss
+{
+       my $hop = shift;                        # minimum hops 
+       my $freq = shift;                       # frequency
+    my $txpower = shift || 20; # transmit power 
+    my $rsens = shift || -123; # receiver sensitivity
+    my $antgain = shift || 0;   # antenna gain
+               
+    my $daynight = shift;              # ref to daynight array one per hop
+    my $beta = shift;
+       my $path = shift;
+       my $mufF = shift;
+       my $mufE = shift;
+    my $absorp = shift;
+    my $dB2 = shift;
+                       
+       my $h;                                          # hop number 
+       my $level;                                      # max signal (dBm) 
+       my $signal;                                     # receive signal (dBm) 
+       my $ftemp;                                      # double temp 
+       my $j;                                          # index temp 
+
+       #
+       # Calculate signal and noise for all hops. The noise level is
+       # -140 dBm for a receiver bandwidth of 2500 Hz and noise
+       # temperature 290 K. The receiver sensitivity is assumed -123
+       # dBm (0.15 V at 50 Ohm for 10 dB S/N). Paths where the signal
+       # is less than the noise or when the frequency exceeds the F-
+       # layer MUF are considered unusable.
+        
+       $level = $noise;
+       $j = 0;
+       for ($h = $hop; $h < $hop + 3; $h++) {
+#              $daynight->[$h] &= ~(P_E | P_S | P_M);
+               if ($freq < 0.85 * $mufF->[$h]) {
+
+                       # Transmit power (dBm)
+                        
+                       $signal = $txpower + $antgain + 30;
+
+                       # Path loss
+                        
+                       $signal -= 32.44 + 20 * log10($path->[$h] * $freq) + $SLOSS;
+
+                       # Ionospheric loss
+                        
+                       $ftemp = $R * cos($beta->[$h]) / ($R + $hE);
+                       $ftemp = atan($ftemp / sqrt(1 - $ftemp * $ftemp));
+                       $signal -= 677.2 * $absorp->[$h] / cos($ftemp) / (pow(($freq + $GAMMA), 1.98) + 10.2);
+
+                       # Ground reflection loss
+                        
+                       $signal -= $h * $GLOSS;
+                       $dB2->[$h] = $signal;
+
+                       # Paths where the signal is greater than the
+                       # noise, but less than the receiver sensitivity
+                       # are marked 's'. Paths below the E-layer MUF
+                       # are marked 'e'. When comparing for maximum
+                       # signal, The signal for these paths is reduced
+                       # by 3 dB so they will be used only as a last
+                       # resort.
+                        
+                       
+                       $daynight->[$h] |= 4 if ($signal < $rsens);
+                       if ($freq < $mufE->[$h]) {
+                               $daynight->[$h] |= 8;
+                               $signal -= $MPATH;
+                       }
+                       if ($signal > $level) {
+                               $level = $signal;
+                               $j = $h;
+                       }
+               }
+       }
+
+       # We have found the best path. If this path is less than 3 dB
+       # above the RMS sum of the other paths, the path is marked 'm'.
+        
+       return 0 if ($j == 0);
+
+       $ftemp = 0;
+       for ($h = $hop; $h < $hop + 3; $h++) {
+               $ftemp += exp(2 / 10 * $dB2->[$h] * $LN10) if ($h != $j);
+       }
+       $ftemp = 10 / 2 * log10($ftemp);
+       $daynight->[$j] |= 16 if ($level < $ftemp + $MPATH);
+       
+       return $j;
+}
+
+# zenith(dist) - Determine sun zenith angle at reflection zone.
+
+sub zenith
+{
+       my $dist = shift;                       # path angle 
+       my $txlat = shift;          # tx latitude (rad)
+       my $txlong = shift;         # tx longitude (rad)
+       my $txbearing = shift;      # tx bearing
+       my $pathangle = shift;      # 'b1'
+       my $lats = shift;           # subsolar latitude
+       my $lons = shift;           # subsolar longitude
+                               
+       my ($latr, $lonr);                      # reflection zone coordinates (rad) 
+       my $thetar;                                     # reflection zone angle (rad) 
+       my $psi;                                        # sun zenith angle (rad) 
+
+       # Calculate reflection zone coordinates.
+        
+       $latr = acos(cos($dist) * sin($txlat) + sin($dist) * cos($txlat) * cos($txbearing));
+       $latr += $pi if ($latr < 0);
+       $latr = $halfpi - $latr;
+       $lonr = acos((cos($dist) - sin($latr) * sin($txlat)) / (cos($latr) * cos($txlat)));
+       $lonr += $pi if ($lonr < 0);
+       $lonr = - $lonr if ($pathangle < 0);
+       $lonr = $txlong - $lonr;
+       $lonr -= $pi2 if ($lonr >= $pi);
+       $lonr += $pi2 if ($lonr <= -$pi);
+       $thetar = $lons - $lonr;
+       $thetar = $pi2 - $thetar if ($thetar > $pi);
+       $thetar -= $pi2 if ($thetar < - $pi);
+       
+       # Calculate sun zenith angle.
+        
+       $psi = acos(sin($latr) * sin($lats) + cos($latr) * cos($lats) * cos($thetar));
+       $psi += $pi if ($psi < 0);
+       return($psi);
+}
+               
+               
+
+1;
index 0c3933ce23cd7857ba28e8f4623a4103c99103d8..d3f32a36b3230192422268947973affcd8ff6047 100644 (file)
@@ -11,7 +11,7 @@
     <address><a href="mailto:djk@tobit.co.uk"></a>Dirk Koopman G1TLH</address><br>
 <!-- Created: Sat Sep  4 21:11:44 BST 1999 -->
 <!-- hhmts start -->
-Last modified: Sat Sep  4 22:50:53 BST 1999
+Last modified: Sat Sep  4 23:50:50 BST 1999
 <!-- hhmts end -->
        <hr>
        <h4>Some thoughts</h4>
@@ -41,15 +41,29 @@ Last modified: Sat Sep  4 22:50:53 BST 1999
                <tr>
                  <td width="20%"><b>OP Code</b></td>
                  <td width="80%">
+                       
                        This will be 'B' for now, but will be extended when multicasting
-                       comes into being
+                       comes into being, the possible values are:-<br>
+                       <table>
+                         <tr><td>A</td><td>ACK</td></tr>
+                         <tr><td>B</td><td>Broadcast</td></tr>
+                         <tr><td>C</td><td>Catch up</td></tr>
+                         <tr><td>J</td><td>Join</td></tr>
+                         <tr><td>L</td><td>Leave</td></tr>
+                         <tr><td>M</td><td>Multicast</td></tr>
+                         <tr><td>N</td><td>NAK</td></tr>
+                         <tr><td>P</td><td>Point to Point</td></tr>
+                         <tr><td>T</td><td>Time (ping with time)</td></tr>
+                       </table>
                  </td>
                </tr>
                <tr>
                  <td width="20%"><b>From Address</b></td>
                  <td width="80%">
-                       This will always be the originating cluster callsign and shall 
-                       <emph>never</emph> be changed
+                       This will always be the originating (cluster) callsign and shall 
+                       <emph>never</emph> be changed. The word cluster is in brackets because
+                       it is envisaged that 'connected' callsigns will eventually speak the
+                       the same protocol as the clusters.
                  </td>
                </tr>
                <tr>
@@ -107,11 +121,11 @@ Last modified: Sat Sep  4 22:50:53 BST 1999
 
        <p>
          A typical packet might look like:- <br><br>
-         <quote>
+         <code>
                BGB7DJK||8BCF65DE|DX^G1TLH^M0BAA^144123^Humungous Signal!|A8<BR>
                BGB7DJK|SYSOP|8BCF65FC|AN^G1TLH^What @$%7C%5E!** condxs?|5C<BR>
                BGB7DJK|GB7BAA|8BCF670102|TA^G1TLH^G8TIC^Baaaaaaaaaaa|FD<BR>
-      </quote>
+      </code>
     <p>
 
     <p>