cmd/show/muf.pl

   1 #!/usr/bin/perl
   2 #
   3 # show/muf command
   4 #
   5 # Copyright (c) 1999 Dirk Koopman G1TLH
   6 #
   7 # as fixed by Steve Franke K9AN
   8 #
   9 #
  10 #
  11
  12 use Minimuf;
  13
  14 my ($self, $line) = @_;
  15 my @f = split /\s+/, $line;
  16
  17 my $prefix = uc shift @f;
  18 return (1, $self->msg('e4')) unless $prefix;
  19 my $lp;
  20 my $hr2;
  21
  22 while (@f) {
  23         my $f = shift @f;
  24         $lp++ if $f =~ /^l/;
  25         $hr2 = $f if $f =~ /^\d+$/;
  26 }
  27
  28 $hr2 = 2 if !$hr2 || $hr2 < 2;
  29 $hr2 = 24 if $hr2 > 24;
  30
  31 my @out;
  32
  33 # get prefix data
  34 my ($pre, $a) = Prefix::extract($prefix);
  35
  36 # calc bearings and distance
  37 my ($d, $b1, $b2);                              # distance, bearing from TX and from RX
  38 my ($lat2, $lon2);              # lats and longs in radians
  39 my $lat1 = $self->user->lat;
  40 my $lon1 = $self->user->long;
  41 my $loc1 = $self->user->qth || "unknown";
  42
  43 if (!$lon1 && !$lat1) {
  44         push @out, $self->msg('heade1');
  45         $lat1 = $main::mylatitude;
  46         $lon1 = $main::mylongitude;
  47         $loc1 = $main::myqth;
  48 }
  49 $lat2 = $a->{lat};
  50 $lon2 = $a->{long};
  51 ($b1, $d) = DXBearing::bdist($lat1, $lon1, $lat2, $lon2);
  52 ($b2, undef) = DXBearing::bdist($lat2, $lon2, $lat1, $lon1);
  53
  54 # convert stuff into radians
  55 $lat1 *= $d2r;
  56 $lat2 *= $d2r;
  57 $lon1 *= -$d2r;
  58 $lon2 *= -$d2r;
  59 $b1 *= $d2r;
  60 $b2 *= $d2r;
  61 $d = ($d / $R);
  62
  63 # handle long path
  64 if ($lp) {
  65         $d = $pi2 - $d;
  66         $b1 += $pi;
  67         $b1 -= $pi2 if ($b1 >= $pi2);
  68         $b2 += $pi;
  69         $b2 -= $pi2 if ($b2 >= $pi2);
  70 }
  71
  72
  73 my ($hr1, $day, $month) = (gmtime($main::systime))[2,3,4];
  74 $month++;
  75 my $flux = Geomag::sfi;
  76 my $ssn = Minimuf::spots($flux);
  77
  78 my $theta;                                              # path angle (rad)
  79 $theta=$lon1-$lon2;
  80 $theta=$theta+2.*$pi if( $theta <= -$pi);
  81 $theta=$theta-2.*$pi if( $theta >= $pi);
  82
  83 my ($lats, $lons);                              # subsolar coordinates (rad)
  84 my $dB1 = 26;                                   # transmitter output power (dBW)
  85
  86 my $delay;                                              # path delay (ms)
  87 my $psi;                                                # sun zenith angle (rad)
  88 my ($ftemp, $gtemp);                    # my $temps
  89 my ($i, $j, $h, $n);                    # int temps
  90 my $offset;                                             # offset for local time (hours)
  91 my $fcF;                                                # F-layer critical frequency (MHz)
  92 my $phiF;                                               # F-layer angle of incidence (rad)
  93 my $hop;                                                # number of ray hops
  94 my $beta1;                                              # elevation angle (rad)
  95 my $dhop;                                               # hop great-circle distance (rad)
  96 my $height;                                             # height of F layer (km)
  97 my $time;                                               # time of day (hour)
  98 my $rsens = -128;                               # RX sensitivity
  99
 100
 101 my @freq = qw(1.8 3.5 7.0 10.1 14.0 18.1 21.0 24.9 28.0 50.0); # working frequencies (MHz)
 102 my $nfreq = @freq;                              # number of frequencies
 103 my @mufE;                                               # maximum E-layer MUF (MHz)
 104 my @mufF;                                               # minimum F-layer MUF (MHz)
 105 my @absorp;                                             # ionospheric absorption coefficient
 106 my @dB2;                                                # receive power (dBm)
 107 my @path;                                               # path length (km)
 108 my @beta;                                               # elevation angle (rad)
 109 my @daynight;                                   # path flags
 110
 111 # calculate hops, elevation angle, F-layer incidence, delay.
 112 $hop = int ($d / (2 * acos($R / ($R + $hF))));
 113 $beta1 = 0;
 114 while ($beta1 < $MINBETA) {
 115         $hop++;
 116         $dhop = $d / ($hop * 2);
 117         $beta1 = atan((cos($dhop) - $R / ($R + $hF)) / sin($dhop));
 118 }
 119 $ftemp = $R * cos($beta1) / ($R + $hF);
 120 $phiF = atan($ftemp / sqrt(1 - $ftemp * $ftemp));
 121 $delay = ((2 * $hop * sin($dhop) * ($R + $hF)) / cos($beta1) / $VOFL) * 1e6;
 122
 123 # print summary of data so far
 124 push @out, sprintf("RxSens: $rsens dBM SFI:%4.0f   R:%4.0f   Month: $month   Day: $day", $flux, $ssn);
 125 push @out, sprintf("Power :  %3.0f dBW    Distance:%6.0f km    Delay:%5.1f ms", $dB1, $d * $R, $delay);
 126 push @out, sprintf("Location                       Lat / Long           Azim");
 127 push @out, sprintf("%-30.30s %-18s    %3.0f", $loc1, DXBearing::lltos($lat1*$r2d, -$lon1*$r2d), $b1 * $r2d);
 128 push @out, sprintf("%-30.30s %-18s    %3.0f", $a->name, DXBearing::lltos($lat2*$r2d, -$lon2*$r2d), $b2 * $r2d);
 129 my $head = "UT LT  MUF Zen";
 130 for ($i = 0; $i < $nfreq; $i++) {
 131         $head .= sprintf "%5.1f", $freq[$i];
 132 }
 133 push @out, $head;
 134
 135 my $hour;
 136
 137 # Hour loop: This loop determines the min-hop path and next two
 138 # higher-hop paths. It selects the most likely path for each
 139 # frequency and calculates the receive power. The F-layer
 140 # critical frequency is computed directly from MINIMUF 3.5 and
 141 # the secant law.
 142
 143 $offset = int ($lon2 * 24. / $pi2);
 144 for ($hour = $hr1; $hour < $hr2+$hr1; $hour++) {
 145     my $dh = $hour;
 146         while ($dh >= 24) {
 147                 $dh -= 24;
 148         };
 149         $time = $dh - $offset;
 150         $time += 24 if ($time < 0);
 151         $time -= 24 if ($time >= 24);
 152         my $out = sprintf("%2.0f %2.0f", $dh, $time);
 153         $ftemp = Minimuf::minimuf($flux, $month, $day, $dh, $lat1, $lon1, $lat2, $lon2);
 154         $fcF = $ftemp * cos($phiF);
 155
 156         # Calculate subsolar coordinates.
 157         $ftemp = ($month - 1) * 365.25 / 12. + $day - 80.;
 158         $lats = 23.5 * $d2r * sin($ftemp / 365.25 * $pi2);
 159         $lons = ($dh * 15. - 180.) * $d2r;
 160
 161         # Path loop: This loop determines the geometry of the
 162         # min-hop path and the next two higher-hop paths. It
 163         # calculates the minimum F-layer MUF, maximum E-layer
 164         # MUF and ionospheric absorption factor for each
 165         # geometry.
 166         for ($h = $hop; $h < $hop + 3; $h++) {
 167
 168                 # We assume the F layer height increases during
 169                 # the day and decreases at night, as determined
 170                 # at the midpoint of the path.
 171                 $height = $hF;
 172                 $psi = Minimuf::zenith($d / 2, $lat1, $lon1, $b1, $theta, $lats, $lons);
 173                 if ($psi < 0) {
 174                         $height -= 70.;
 175                 } else {
 176                         $height += 30;
 177                 }
 178                 $dhop = $d / ($h * 2.);
 179                 $beta[$h] = atan((cos($dhop) - $R / ($R + $height)) / sin($dhop));
 180                 $path[$h] = 2 * $h * sin($dhop) * ($R + $height) / cos($beta[$h]);
 181                 Minimuf::ion($h, $d, $fcF, $ssn, $lat1, $lon1, $b1, $theta, $lats, $lons, \@daynight, \@mufE, \@mufF, \@absorp);
 182         }
 183
 184         # Display one line for this hour.
 185         $out .= sprintf("%5.1f%4.0f ", $mufF[$hop], 90 - $psi * $r2d);
 186         $ftemp = $noise;
 187         for ($i = 0; $i < $nfreq; $i++) {
 188                 $n = Minimuf::pathloss($hop, $freq[$i], 20, $rsens, 0,  \@daynight, \@beta, \@path, \@mufF, \@mufE, \@absorp, \@dB2);
 189                 my $s = Minimuf::ds($n, $rsens, \@dB2, \@daynight);
 190                 $out .= " $s";
 191         }
 192         $out =~ s/\s+$//;
 193         push @out, $out;
 194 }
 195
 196 return (1, @out);