mirror of
git://scm.dxcluster.org/scm/spider
synced 2024-09-21 15:57:12 +00:00
489 lines
14 KiB
Perl
489 lines
14 KiB
Perl
#!/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
|
|
#
|
|
# as fixed by Steve Franke K9AN
|
|
#
|
|
# 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;
|
|
|
|
|
|
require Exporter;
|
|
@ISA = qw(Exporter);
|
|
@EXPORT = qw($pi $d2r $r2d $halfpi $pi2 $VOFL $R $hE $hF $GAMMA $LN10
|
|
$MINBETA $BOLTZ $NTEMP $DELTAF $MPATH $GLOSS $SLOSS
|
|
$noise);
|
|
|
|
use strict;
|
|
|
|
use DXDebug;
|
|
use POSIX qw(:math_h);
|
|
|
|
use vars qw($VERSION $BRANCH);
|
|
$VERSION = sprintf( "%d.%03d", q$Revision$ =~ /(\d+)\.(\d+)/ );
|
|
$BRANCH = sprintf( "%d.%03d", q$Revision$ =~ /\d+\.\d+\.(\d+)\.(\d+)/ || (0,0));
|
|
$main::build += $VERSION;
|
|
$main::branch += $BRANCH;
|
|
|
|
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 $lat1 = shift;
|
|
my $lon1 = shift;
|
|
my $b1 = shift;
|
|
my $b2 = shift;
|
|
my $lats = shift;
|
|
my $lons = shift;
|
|
|
|
# various refs to arrays
|
|
my $daynight = shift; # ref to daynight array one per hop
|
|
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));
|
|
$absorp->[$h] = $mufE->[$h] = $daynight->[$h] = 0;
|
|
$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, $lat1, $lon1, $b1, $b2, $lats, $lons);
|
|
$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] &= ~(4 | 8 | 16);
|
|
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);
|
|
}
|
|
|
|
# official minimuf version of display
|
|
sub dsx
|
|
{
|
|
my $h = shift;
|
|
my $rsens = shift;
|
|
my $dB2 = shift;
|
|
my $daynight = shift;
|
|
|
|
my $c1;
|
|
my $c2;
|
|
|
|
return " " unless $h;
|
|
|
|
if (($daynight->[$h] & 3) == 3) {
|
|
$c1 = 'x';
|
|
} elsif ($daynight->[$h] & 1) {
|
|
$c1 = 'j';
|
|
} elsif ($daynight->[$h] & 2) {
|
|
$c1 = 'n';
|
|
}
|
|
if ($daynight->[$h] & 4) {
|
|
$c2 = 's';
|
|
} elsif ($daynight->[$h] & 16) {
|
|
$c2 = 'm';
|
|
} else {
|
|
$c2 = ' ';
|
|
}
|
|
return sprintf("%4.0f%s%1d%s", $dB2->[$h] - $rsens, $c1, $h, $c2)
|
|
}
|
|
|
|
# my version
|
|
sub ds
|
|
{
|
|
my $h = shift;
|
|
my $rsens = shift;
|
|
my $dB2 = shift;
|
|
my $daynight = shift;
|
|
|
|
my $c2;
|
|
|
|
return " " unless $h;
|
|
|
|
if ($daynight->[$h] & 4) {
|
|
$c2 = 's';
|
|
} elsif ($daynight->[$h] & 16) {
|
|
$c2 = 'm';
|
|
} else {
|
|
$c2 = ' ';
|
|
}
|
|
my $l = $dB2->[$h] - $rsens;
|
|
my $s = int $l / 6;
|
|
$s = 9 if $s > 9;
|
|
$s = 0 if $s < 0;
|
|
my $plus = (($l / 6) >= $s + 0.5) ? '+' : ' ';
|
|
|
|
return "$c2". "S$s$plus";
|
|
}
|
|
|
|
1;
|