From 22970505cb6d1babd4ddb5af18c15c408f3cf7ac Mon Sep 17 00:00:00 2001 From: djk Date: Tue, 19 Oct 1999 23:48:19 +0000 Subject: [PATCH] started the process of adding minimuf, translated the minimuf subroutines and created the Minimuf.pm file --- Changes | 2 + cmd/show/muf.pl | 0 perl/Minimuf.pm | 407 ++++++++++++++++++++++++++++++++++++++++++++++++ perl/proto.html | 26 +++- 4 files changed, 429 insertions(+), 6 deletions(-) create mode 100644 cmd/show/muf.pl create mode 100644 perl/Minimuf.pm diff --git a/Changes b/Changes index 93142c22..06b5c67f 100644 --- 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 index 00000000..e69de29b diff --git a/perl/Minimuf.pm b/perl/Minimuf.pm new file mode 100644 index 00000000..e493e5b9 --- /dev/null +++ b/perl/Minimuf.pm @@ -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; diff --git a/perl/proto.html b/perl/proto.html index 0c3933ce..d3f32a36 100644 --- a/perl/proto.html +++ b/perl/proto.html @@ -11,7 +11,7 @@
Dirk Koopman G1TLH

-Last modified: Sat Sep 4 22:50:53 BST 1999 +Last modified: Sat Sep 4 23:50:50 BST 1999

Some thoughts

@@ -41,15 +41,29 @@ Last modified: Sat Sep 4 22:50:53 BST 1999 OP Code + This will be 'B' for now, but will be extended when multicasting - comes into being + comes into being, the possible values are:-
+ + + + + + + + + + +
AACK
BBroadcast
CCatch up
JJoin
LLeave
MMulticast
NNAK
PPoint to Point
TTime (ping with time)
From Address - This will always be the originating cluster callsign and shall - never be changed + This will always be the originating (cluster) callsign and shall + never 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. @@ -107,11 +121,11 @@ Last modified: Sat Sep 4 22:50:53 BST 1999

A typical packet might look like:-

- + BGB7DJK||8BCF65DE|DX^G1TLH^M0BAA^144123^Humungous Signal!|A8
BGB7DJK|SYSOP|8BCF65FC|AN^G1TLH^What @$%7C%5E!** condxs?|5C
BGB7DJK|GB7BAA|8BCF670102|TA^G1TLH^G8TIC^Baaaaaaaaaaa|FD
- +

-- 2.43.0