Generated: Sat Oct 24 16:35:15 2020 from fg-navaid.pl 2019/11/20 42.6 KB. text copy
#! /usr/bin/env perl # from : https://github.com/kroki/fg-navaid # 2019-11-20 - Update - added year 2015 # 2013-09-23 - Initial copy use v5.12; use strict; use warnings; use utf8; use open qw(:std :utf8); use Encode qw(decode_utf8); @ARGV = map { decode_utf8($_) } @ARGV; =pod =encoding utf8 =head1 NAME kroki-fg-navaid - navigation helper for FlightGear =head1 SYNOPSIS kroki-fg-navaid [OPTIONS] -- LOCATION =head2 LOCATION Location can be specified either as an I<IDENT> of an object (airport, VOR, fix, etc.), runway as I<ARPT-RW> (like C<KSFO-28L>), or as a pair I<LATITUDE,LONGITUDE>. I<LATITUDE> and I<LONGITUDE> of a location may be given in a number of formats, including N35°23'42.36" -136o 51' 14.04'' s28.567894 w46*36" E 235 * 4.5678 ' I.e. northern hemisphere is specified by prefix B<N>, B<n>, B<+>, or no prefix; southern hemisphere is specified by prefix B<S>, B<s>, or B<->; eastern hemisphere by B<E>, B<e>, B<+> or no prefix; western hemisphere by B<W>, B<w>, or B<->. Coordinates may be given either as a decimal degrees or as a triplet of degrees (denoted by B<°> or B<*> or B<o>), minutes (denoted by B<'>), seconds (denoted by B<"> or B<''>). Each part is optional and the last one may be fractional. Whitespaces are ignored. Coordinates are normalized, for instance B<N91,E096> becomes B<N89,W084>. =head2 OPTIONS =over =item C<-m, --metric> By default distances and heights are given in nautical miles (B<nm>) and feet (B<'>). C<--metric> enables the output in kilometers (B<km>) and meters (B<m>). =item C<-b, --bearings[=NUM]> Report bearings from/to this many beacons (from VORs or to NDBs). For VORs bearing I<from> the station (i.e. current radial) is reported, for NDBs magnetic heading I<to> the beacon is reported. In both cases true heading and the (ground) distance to the station are also reported. When C<NUM> is not provided default is 2. =item C<-r, --range=NM|KM> Maximum range to look for nearest stations. Specified in nautical miles or kilometers depending on the presence of C<--metric>. Default is 200nm (370.4km). =item C<-R, --rsbns[=NUM]> Report Sm, Zm coordinates and UK of RSBN station for NVU correction (think of tu154b aircraft). C<NUM> specifies the number of candidate stations to report and defaults to 1. This option currently works only for runways. =item C<-C, --navdata-cache=PATH> Optional path to I<navdata> cache file. Default is most recent cache file in C<~/.fgfs/>. =item C<-W, --wmm=YEAR> World Magnetic Model coefficients version. Available are I<2005>, I<2010> and I<2015>. Default is I<2005> for compatibility with FlightGear 2.12. =item C<-h, --help> Print this message and exit. =back =head1 BUGS Report bugs at <https://github.com/kroki/fg-navaid/issues> or directly to <tomash.brechko@gmail.com>. =head1 COPYRIGHT & LICENSE Copyright (C) 2013 Tomash Brechko. All rights reserved. This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program. If not, see <http://www.gnu.org/licenses/>. =cut use Getopt::Long qw(:config gnu_getopt); use Pod::Usage; my %option = ( wmm => 2005, ); if (! GetOptions(\%option, qw(metric|m bearings|b:2 range|r=f rsbns|R:1 navdata-cache|C=s wmm|W=i help|h)) || @ARGV != 1) { pod2usage(1); } if (exists $option{help}) { pod2usage(0); exit; } $option{range} //= ($option{metric} ? 370.4 : 200); my $os = $^O; # degree symbol prints badly on windows my $degs = '°'; # $degs = chr(248) if ($os =~ /win/i); # nope, still bad #$degs = 'd' if ($os =~ /win/i); # hmmm, 'd' does not look good $degs = ' ' if ($os =~ /win/i); # maybe space is the best? unless (exists $option{'navdata-cache'}) { my $major = 0; my $minor = 0; if ($os =~ /win/i) { foreach my $cache (<$ENV{APPDATA}\\flightgear.org\\navdata_*_*.cache>) { ### print "Checking [$cache]\n"; if ($cache =~ m|navdata_(\d+)_(\d+)\.cache$| && ($1 > $major || ($1 == $major && $2 > $minor))) { $major = $1; $minor = $2; ### print "Set cache [$cache]\n"; $option{"navdata-cache"} = $cache; } } } else { foreach my $cache (<$ENV{HOME}/.fgfs/navdata_*_*.cache>) { if ($cache =~ m|/navdata_(\d+)_(\d+)\.cache$| && ($1 > $major || ($1 == $major && $2 > $minor))) { $major = $1; $minor = $2; $option{"navdata-cache"} = $cache; } } } } die "No navdata cache file found\n" unless exists $option{'navdata-cache'} && -f $option{'navdata-cache'}; print "Using 'navdata-cache' at ".$option{'navdata-cache'}."\n"; # show which cache used sub NM { sub { ($option{metric} ? ($_[0] * 1.852, "km") : ($_[0], "nm")) } } sub M { sub { ($option{metric} ? ($_[0], "m") : ($_[0] / 0.3048, "'")) } } use DBI; my $dsn = "DBI:SQLite:dbname=$option{'navdata-cache'}"; my $attrs = { ReadOnly => 1, PrintError => 0, RaiseError => 1 }; my $dbh = DBI->connect($dsn, '', '', $attrs); my @id2type = qw( none AIRPORT HELIPORT SEAPORT RUNWAY HELIPAD TAXIWAY PAVEMENT WAYPOINT FIX NDB VOR ILS LOC GS OM MM IM DME TACAN MOBILE_TACAN OBSTACLE TOWER FREQ_GROUND FREQ_TOWER FREQ_ATIS FREQ_AWOS FREQ_APP_DEP FREQ_ENROUTE FREQ_CLEARANCE FREQ_UNICOM PARKING TAXI_NODE COUNTRY CITY TOWN VILLAGE ); my %type2id = map { ($id2type[$_], $_) } (0..$#id2type); my @id2surface = ( "(none)", "asphalt", "concrete", "turf/grass", "dirt", "gravel", "asphalt", "concrete", "turf/grass", "dirt", "asphalt", "concrete", "dry lakebed", "water", "snow/ice", "(unknown)" ); my %info = ( AIRPORT => \&airport_info, HELIPORT => \&airport_info, SEAPORT => \&airport_info, WAYPOINT => \&fix_info, FIX => \&fix_info, NDB => \&nav_info, VOR => \&nav_info, ILS => \&nav_info, LOC => \&nav_info, DME => \&nav_info, TACAN => \&nav_info, MOBILE_TACAN => \&nav_info, ); my $loc_info = $dbh->prepare(qq{ SELECT rowid, type, name, lon, lat, elev_m, cart_x, cart_y, cart_z FROM positioned WHERE ident = ? AND type IN (@{[ join ",", @type2id{(keys %info)} ]}) ORDER BY type }); my $airport_info = $dbh->prepare(qq{ SELECT p.rowid, p.type, p.ident, p.name, r.length_ft AS length_m, r.width_m, r.surface, r.reciprocal, r.ils, c.freq_khz, c.range_nm FROM positioned p LEFT JOIN runway AS r ON p.rowid = r.rowid LEFT JOIN comm AS c ON p.rowid = c.rowid WHERE p.airport = ? AND p.type IN (@{[ join ",", @type2id{qw(RUNWAY HELIPAD FREQ_GROUND FREQ_TOWER FREQ_ATIS FREQ_AWOS FREQ_APP_DEP FREQ_ENROUTE FREQ_CLEARANCE FREQ_UNICOM)} ]}) ORDER BY p.type, p.ident, p.name }); my $fix_info = $dbh->prepare(q{ SELECT DISTINCT a.ident, CASE ae.network WHEN 1 THEN "low" ELSE "high" END FROM airway_edge AS ae JOIN airway AS a ON ae.airway = a.rowid WHERE ae.a = ? ORDER BY a.rowid, a.network }); my $nav_info = $dbh->prepare(q{ SELECT freq, range_nm, multiuse, colocated FROM navaid WHERE rowid = ? }); my $runway_info = $dbh->prepare(qq{ SELECT p.rowid, p.type, p.lon, p.lat, p.elev_m, p.cart_x, p.cart_y, p.cart_z, r.heading, r.length_ft AS length_m, r.width_m, r.surface FROM positioned AS p INDEXED BY pos_apt_type -- force the right index JOIN runway AS r ON p.rowid = r.rowid WHERE p.airport = ? AND p.type = $type2id{RUNWAY} AND p.ident = ? }); my $ils_info = $dbh->prepare(qq{ SELECT p.type, p.ident, p.name, p.lon, p.lat, p.elev_m, n.freq, n.range_nm, n.multiuse, d.multiuse AS dme_bias, dp.lon As dme_lon, dp.lat AS dme_lat FROM positioned AS p JOIN navaid AS n ON p.rowid = n.rowid LEFT JOIN navaid AS d ON n.colocated = d.rowid AND n.freq = d.freq LEFT JOIN positioned AS dp ON d.rowid = dp.rowid WHERE p.airport = ? AND p.type IN (@{[ join ",", @type2id{qw(ILS LOC GS OM MM IM)} ]}) AND n.runway = ? }); my $octree_children = $dbh->prepare(q{ SELECT children FROM octree WHERE rowid = ? }); use Geo::Inverse; my $igeo = Geo::Inverse->new(); my $magdecl = MagneticDeclination->new($igeo->ellipsoid, $option{wmm}); my ($location) = @ARGV; $loc_info->execute($location); my %colocated; while (my $loc = $loc_info->fetchrow_hashref) { next if delete $colocated{$loc->{rowid}}; say "--"; my $typename = $id2type[$loc->{type}]; my $md; if ($typename ne 'TACAN' && $typename ne 'MOBILE_TACAN') { $md = $magdecl->compute($loc->{lat}, $loc->{lon}); say "Coordinates: ", coords2str($loc->{lat}, $loc->{lon}); say "Magnetic declination: ", magdecl2str($md); } say "Ident: $location"; my $info = $info{$typename}; $info->($loc, $typename, \%colocated) if $info; bearings($loc, $md) if defined $md; } my ($apt, $rwy) = split /-/, $location, 2; runway_info($apt, $rwy) if defined $apt and defined $rwy; use POSIX qw(fmod); my ($lat, $lon) = map { str2deg($_) } split /,/, $location, 2; if (defined $lat and defined $lon) { $lat = fmod($lat, 180); if ($lat > 90) { $lat = -$lat + 180; $lon += 180; } elsif ($lat < -90) { $lat = -$lat - 180; $lon += 180; } $lon = fmod($lon, 360); if ($lon > 180) { $lon = $lon - 360; } elsif ($lon < -180) { $lon = $lon + 360; } my $md = $magdecl->compute($lat, $lon); say "Coordinates: ", coords2str($lat, $lon); say "Magnetic declination: ", magdecl2str($md); my $loc = { rowid => 0, lat => $lat, lon => $lon, }; @$loc{qw(cart_x cart_y cart_z)} = geo2xyz($lat, $lon); bearings($loc, $md); } exit; sub str2deg { my ($str) = @_; use Regexp::Grammars; my $deg = qr{ ^ <val> $ <rule: val> <sign>? <deg> <MATCH=(?{ ($MATCH{sign}//1) * ($MATCH{deg}//0) })> <rule: sign> <pos> <MATCH=1> | <neg> <MATCH=-1> <rule: pos> \+ | N | E <rule: neg> \- | S | W <rule: deg> <MATCH=float> [°*o]? | (?: <int> [°*o] )? <min> <MATCH=(?{ ($MATCH{int}//0) + ($MATCH{min}//0) / 60 })> <rule: min> <MATCH=float> \'? | (?: <int> \' )? <sec> <MATCH=(?{ ($MATCH{int}//0) + ($MATCH{sec}//0) / 60 })> <rule: sec> <MATCH=float> (?: \" | \'\' )? | <MATCH=int> (?: \" | \'\' ) | <MATCH=0> <token: float> \d+(?:\.\d*)? | \.\d+ <token: int> \d+ }ix; if ($str =~ $deg) { return $/{val}; } else { return undef; } } sub deg2dms2f { my ($pos, $neg, $deg) = @_; my $sign = $pos; if ($deg < 0) { $sign = $neg; $deg = -$deg; } my $d = int($deg * 360000 + 0.5); my $s = ($d % 6000) / 100; $d = int($d / 6000); my $m = $d % 60; $d = int($d / 60); return ($sign, $d, $m, $s); } sub deg2dm { my ($pos, $neg, $deg) = @_; my $sign = $pos; if ($deg < 0) { $sign = $neg; $deg = -$deg; } my $d = int($deg * 60 + 0.5); my $m = $d % 60; $d = int($d / 60); return ($sign, $d, $m); } sub coords2str { my ($lat, $lon) = @_; my $form ="%s%02d${degs}%02d'%05.2f\", %s%03d${degs}%02d'%05.2f\" (%.7f, %.7f)"; sprintf($form, deg2dms2f('N', 'S', $lat), deg2dms2f('E', 'W', $lon), $lat, $lon); } sub magdecl2str { my ($md) = @_; my $form = "%s%d${degs}%02d' (%.2f)"; sprintf($form, deg2dm('E', 'W', $md), $md); } sub dist2str { my ($d, $units, $decimals) = @_; my ($v, $u) = $units->($d); my $scale = 10**($decimals//0); $v = int($v * $scale) / $scale; sprintf("%.*f%s", ${decimals}//0, $v, $u) } sub mhz2str { my ($mhz) = @_; my $s = sprintf("%.3f", $mhz); $s =~ s/(?<!\.)0+$//; $s . "mhz" } sub khz2str { my ($khz) = @_; $khz . "khz" } sub airport_info { my ($loc, $typename) = @_; say "Object: $loc->{name} ", lc $typename; say "Elevation: ", dist2str($loc->{elev_m}, M); my %reciprocal; my %comm; $airport_info->execute($loc->{rowid}); while (my $info = $airport_info->fetchrow_hashref) { my $tn = $id2type[$info->{type}]; if ($tn eq 'RUNWAY') { my $r = delete $reciprocal{$info->{rowid}}; if ($r) { my $ils = ""; if ($r->{ils} && $info->{ils}) { $ils = ", ILS on both"; } elsif ($r->{ils}) { $ils = ", ILS on $r->{ident}"; } elsif ($info->{ils}) { $ils = ", ILS on $info->{ident}"; } say "Runway $r->{ident}-$info->{ident}: ", dist2str($info->{length_m}, M), " x ", dist2str($info->{width_m}, M), ", $id2surface[$info->{surface}]", $ils; } else { $reciprocal{$info->{reciprocal}} = $info; } } elsif ($tn eq 'HELIPAD') { say "Helipad $info->{ident}: ", dist2str($info->{length_m}, M), " x ", dist2str($info->{width_m}, M), ", $id2surface[$info->{surface}]"; } else { push @{$comm{$info->{name}, $info->{range_nm}}}, $info->{freq_khz}; } } foreach my $c (sort keys %comm) { my ($n, $r) = split "\034", $c, 2; my $f = $comm{$c}; say "$n (", dist2str($r, NM), "): ", join(", ", map { mhz2str($_ / 1000) } @$f); } } sub fix_info { my ($loc, $typename) = @_; say "Object: ", lc $typename; say "Elevation: ", dist2str($loc->{elev_m}, M) if $loc->{elev_m}; my $info = $dbh->selectall_arrayref($fix_info, undef, $loc->{rowid}); if (@$info) { my %info; foreach my $i (@$info) { push @{$info{$i->[0]}}, $i->[1]; } say "Airways: ", join(", ", map { "$_ (" . join(", ", @{$info{$_}}) . ")" } sort keys %info); } } sub nav_info { my ($loc, $typename, $colocated) = @_; say "Object: $loc->{name}"; say "Elevation: ", dist2str($loc->{elev_m}, M); my $info = $dbh->selectrow_hashref($nav_info, undef, $loc->{rowid}); if ($typename eq 'VOR') { say "Orientation: ", magdecl2str($info->{multiuse}); } my $freq = \&mhz2str; if ($typename eq 'NDB' || ($typename eq 'DME' && $info->{freq} % 100 == 0)) { $freq = \&khz2str; } say "Frequency: ", $freq->($info->{freq} / 100), " (", dist2str($info->{range_nm}, NM), ")"; if ($info->{colocated}) { my $ci = $dbh->selectrow_hashref($nav_info, undef, $info->{colocated}); if ($ci->{freq} == $info->{freq} && $ci->{range_nm} == $info->{range_nm}) { $colocated->{$info->{colocated}} = 1; } } } sub deg_round { my ($deg) = @_; $deg = sprintf("%.0f", $deg); $deg += 360 while $deg < 0; $deg -= 360 while $deg >= 360; return $deg; } sub heading2str { my ($heading, $md, $sfx) = @_; $sfx //= 'M'; my $form = "%03d${degs}$sfx (%03d${degs}T)"; sprintf($form, deg_round($heading - $md), deg_round($heading)) } sub runway_info { my ($apt, $rwy) = @_; $loc_info->execute($apt); while (my $loc = $loc_info->fetchrow_hashref) { my $typename = $id2type[$loc->{type}]; next if $typename ne 'AIRPORT' and $typename ne 'HELIPORT' and $typename ne 'SEAPORT'; $runway_info->execute($loc->{rowid}, $rwy); while (my $r = $runway_info->fetchrow_hashref) { $r->{elev_m} ||= $loc->{elev_m}; say "--"; my $md = $magdecl->compute($r->{lat}, $r->{lon}); say "Coordinates: ", coords2str($r->{lat}, $r->{lon}); say "Magnetic declination: ", magdecl2str($md); say "Ident: $apt $rwy"; say "Object: $loc->{name} ", lc $typename, ", ", lc $id2type[$r->{type}], " $rwy"; say "Surface: $id2surface[$r->{surface}]"; say "Dimensions: ", dist2str($r->{length_m}, M), " x ", dist2str($r->{width_m}, M); say "Elevation: ", dist2str($r->{elev_m}, M); say "Heading: ", heading2str($r->{heading}, $md); my %ils; $ils_info->execute($loc->{rowid}, $r->{rowid}); while (my $i = $ils_info->fetchrow_hashref) { $ils{$id2type[$i->{type}]} = $i; } my $dme_bias; my $dme_dist; foreach my $n (qw(LOC ILS)) { my $i = $ils{$n} or next; $dme_bias = $i->{dme_bias}; if (defined $dme_bias) { my (undef, undef, $dist) = $igeo->inverse($i->{dme_lat}, $i->{dme_lon}, $r->{lat}, $r->{lon}); $dme_dist = $dist / 1852 - $dme_bias; } my $name; if ($n eq 'ILS') { $name = $i->{name}; $name =~ s/.+ //; } else { $name = "Localizer"; } say "$name (", dist2str($i->{range_nm}, NM), "): $i->{ident} ", mhz2str($i->{freq} / 100), " ", heading2str($i->{multiuse}, $md); } my $gs = $ils{GS}; my $touchdown; my $tan_slope; if ($gs) { $gs->{elev_m} = $r->{elev_m}; my (undef, $baz, $dist) = $igeo->inverse($gs->{lat}, $gs->{lon}, $r->{lat}, $r->{lon}); $touchdown = $dist * cos(Math::Trig::deg2rad($baz - $r->{heading})); my $deg = $gs->{multiuse} / 100000; $tan_slope = Math::Trig::tan(Math::Trig::deg2rad($deg)); my $slope = sprintf "%.1f", $deg; $slope =~ s/\.?0+$//; say "Glideslope: $slope${degs}"; } foreach my $n (["OM", "Outer"], ["MM", "Middle"], ["IM", "Inner"]) { my $i = $ils{$n->[0]} or next; $i->{elev_m} ||= $r->{elev_m}; my (undef, undef, $dist) = $igeo->inverse($i->{lat}, $i->{lon}, $r->{lat}, $r->{lon}); print "$n->[1] marker: "; if (defined $dme_dist) { print "DME ", dist2str($dist / 1852 + $dme_dist, NM, 1), ", "; } if ($gs) { my $h = ($dist + $touchdown) * $tan_slope; $h -= $i->{elev_m} - $gs->{elev_m}; print "RA ", dist2str($h, M), ", "; } print dist2str($dist / 1852, NM, 1), " to runway edge\n"; } if (defined $dme_dist or $gs) { print "Runway edge: "; if (defined $dme_dist) { print "DME ", dist2str($dme_dist, NM, 1); } if ($gs) { my $h = $touchdown * $tan_slope; $h -= $r->{elev_m} - $gs->{elev_m}; print ", " if defined $dme_dist; print "RA ", dist2str($h, M); } print "\n"; } if ($gs) { say "Touchdown: ", dist2str($touchdown / 1852, NM, 1), " past runway edge"; } bearings($r, $md); rsbns($r, $md, $r->{heading}); } } } sub octree_nodes { my ($p, $r2, $id, $b, $e, $d2, $nodes) = @_; $octree_children->execute($id); my ($children) = $octree_children->fetchrow_array; if ($children) { for my $o (0..7) { next unless $children & (1 << $o); my @b = @$b; for my $c (0..2) { $b[$c] -= $e if $o & (1 << $c); } my $d2 = 0; for my $c (0..2) { my $d = $p->[$c] - $b[$c]; $d = -$d - $e if $d < 0; $d2 += $d**2 if $d > 0; } if ($d2 <= $r2) { octree_nodes($p, $r2, ($id << 3) | $o, \@b, $e / 2, $d2, $nodes); } } } else { push @$nodes, [$id, $d2]; } } sub nav_in_range { my ($loc, $count, $types, $has_colocated) = @_; return unless $count; my $range = $option{range} * ($option{metric} ? 1000 : 1852); octree_nodes([@$loc{qw(cart_x cart_y cart_z)}], $range**2, 1, [7000000, 7000000, 7000000], 7000000, 0, \my @nodes); my $nav_in_node = $dbh->prepare_cached(qq{ SELECT p.rowid, p.type, p.ident, p.name, p.lon, p.lat, p.elev_m, n.freq, n.range_nm, n.multiuse, n.colocated FROM positioned AS p JOIN navaid AS n ON p.rowid = n.rowid WHERE p.octree_node = ? AND p.type IN (@{[ join ",", @type2id{@$types} ]}) }); my @nav; foreach my $n (sort { $a->[1] <=> $b->[1] } @nodes) { last if $n->[1] > $range**2; $nav_in_node->execute($n->[0]); while (my $ni = $nav_in_node->fetchrow_hashref) { next if $ni->{rowid} == $loc->{rowid}; next if $has_colocated and not $ni->{colocated}; my @geo = $igeo->inverse($loc->{lat}, $loc->{lon}, $ni->{lat}, $ni->{lon}); next if $geo[2] > $range; next if 1 / ($geo[2]/1852 / $ni->{range_nm})**2 <= 0.2; push @geo, $ni; my $i = 0; ++$i while $i < @nav and $nav[$i][2] <= $geo[2]; splice @nav, $i, 0, \@geo; if (@nav > $count) { pop @nav; $range = $nav[-1][2]; } } } return @nav; } sub bearings { my ($loc, $md) = @_; my @nav = nav_in_range($loc, $option{bearings}, ['VOR', 'NDB']); foreach my $n (@nav) { my ($faz, $baz, $dist, $i) = @$n; my $tn = $id2type[$i->{type}]; my ($freq, $bearing); my $dir; if ($tn eq 'VOR') { $dir = "from"; $freq = mhz2str($i->{freq} / 100); $bearing = heading2str($baz, $i->{multiuse}, 'R'); } else { $dir = "to"; $freq = khz2str($i->{freq} / 100); $bearing = heading2str($faz, $md); } printf("Bearing $dir $i->{name} ($i->{ident} $freq %s): $bearing %s\n", dist2str($i->{range_nm}, NM), dist2str($dist / 1852, NM, 1)); } } sub geo2xyz { my ($lat, $lon) = @_; $lat = Math::Trig::deg2rad($lat); $lon = Math::Trig::deg2rad($lon); my $cos_lat = cos($lat); my $sin_lat = sin($lat); my $cos_lon = cos($lon); my $sin_lon = sin($lon); my $e = $igeo->ellipsoid; my $v = $e->a / sqrt(1 - $e->e2 * $sin_lat**2); my $x = $v * $cos_lat * $cos_lon; my $y = $v * $cos_lat * $sin_lon; my $z = $v * (1 - $e->e2) * $sin_lat; return ($x, $y, $z); } sub rsbns { my ($loc, $md, $heading) = @_; my @vor = nav_in_range($loc, $option{rsbns}, ['VOR'], 1); if (@vor) { my $zpu = $heading - $md; $zpu += 360 while $zpu < 0; $zpu -= 360 while $zpu >= 360; printf("NVU heading: ZPU:%03d${degs}%02d'M UK:%03d${degs}%02d'T\n", (deg2dm(undef, undef, $zpu))[1, 2], (deg2dm(undef, undef, $heading))[1, 2]); } foreach my $v (@vor) { my ($faz, $baz, $dist, $i) = @$v; my $Sm = $dist * cos(Math::Trig::deg2rad($faz - $heading)) / 1000; my $Zm = $dist * sin(Math::Trig::deg2rad($faz - $heading)) / 1000; printf("$i->{name} ($i->{ident} %s %s): Sm:%+.2fkm Zm:%+.2fkm (%s)\n", mhz2str($i->{freq} / 100), dist2str($i->{range_nm}, NM), $Sm, $Zm, dist2str($dist / 1852, NM, 1)); } } ###################################################################### # # MagneticDeclination computes magnetic declination of a location. # # Computations are copied from simgear/magvar/coremag.cxx. # package MagneticDeclination; use DateTime; use Math::Trig qw(:pi); use constant NMAX => 12; sub new { my $class = shift; my ($ellipsoid, $wmm) = @_; my $self = { r_0 => 6371.2, a => $ellipsoid->a / 1000, b => $ellipsoid->b / 1000, }; foreach my $n (1..NMAX) { foreach my $m (0..NMAX) { $self->{gnm}[$n][$m] = 0; $self->{hnm}[$n][$m] = 0; } } my $yearfrac = (DateTime->now->mjd - DateTime->new(year => $wmm)->mjd) / 365.25; while (<DATA>) { if (/\sWMM-$wmm\s/ .. /^9{10}/) { my @c = split; next unless @c == 6; $self->{gnm}[$c[0]][$c[1]] = $c[2] + $c[4] * $yearfrac; $self->{hnm}[$c[0]][$c[1]] = $c[3] + $c[5] * $yearfrac; } } foreach my $n (2..NMAX) { $self->{root}[$n] = sqrt((2*$n - 1) / (2*$n)); } foreach my $m (0..NMAX) { foreach my $n (($m < 2 ? 2 : $m + 1)..NMAX) { $self->{roots}[$m][$n] = [ sqrt(($n - 1)**2 - $m**2), 1 / sqrt($n**2 - $m**2) ]; } } return bless($self, $class); } sub compute { my $self = shift; my ($lat, $lon) = @_; $lat = Math::Trig::deg2rad($lat); $lon = Math::Trig::deg2rad($lon); my ($a, $b, $root, $roots, $gnm, $hnm, $r_0) = @$self{qw(a b root roots gnm hnm r_0)}; my $sinlat = sin($lat); my $coslat = cos($lat); my $theta = atan2($coslat * $a**2, $sinlat * $b**2); my $r = sqrt(($a**4 - ($a**4 - $b**4) * $sinlat**2) / ($a**2 - ($a**2 - $b**2) * $sinlat**2)); my $c = cos($theta); my $s = sin($theta); my $inv_s = 1 / ($s or 1e-8); my (@P, @DP); foreach my $n (0..NMAX) { foreach my $m (0..NMAX) { $P[$n][$m] = 0; $DP[$n][$m] = 0; } } $P[0][0] = 1; $P[1][0] = $c; $P[1][1] = $s; $DP[0][0] = 0; $DP[1][0] = -$s; $DP[1][1] = $c; foreach my $n (2..NMAX) { $P[$n][$n] = $P[$n-1][$n-1] * $s * $root->[$n]; $DP[$n][$n] = ($DP[$n-1][$n-1] * $s + $P[$n-1][$n-1] * $c) * $root->[$n]; } foreach my $m (0..NMAX) { foreach my $n (($m < 2 ? 2 : $m + 1)..NMAX) { $P[$n][$m] = ($P[$n-1][$m] * $c * (2*$n - 1) - $P[$n-2][$m] * $roots->[$m][$n][0]) * $roots->[$m][$n][1]; $DP[$n][$m] = (($DP[$n-1][$m] * $c - $P[$n-1][$m] * $s) * (2*$n - 1) - $DP[$n-2][$m] * $roots->[$m][$n][0]) * $roots->[$m][$n][1]; } } my (@sm, @cm); foreach my $m (0..NMAX) { $sm[$m] = sin($m * $lon); $cm[$m] = cos($m * $lon); } my ($B_r, $B_theta, $B_phi) = (0, 0, 0); my $fn_0 = $r_0 / $r; my $fn = $fn_0**2; foreach my $n (1..NMAX) { my ($c1_n, $c2_n, $c3_n) = (0, 0, 0); foreach my $m (0..NMAX) { my $t = $gnm->[$n][$m] * $cm[$m] + $hnm->[$n][$m] * $sm[$m]; $c1_n += $t * $P[$n][$m]; $c2_n += $t * $DP[$n][$m]; $c3_n += $m * ($gnm->[$n][$m] * $sm[$m] - $hnm->[$n][$m] * $cm[$m]) * $P[$n][$m]; } $fn *= $fn_0; $B_r += ($n + 1) * $c1_n * $fn; $B_theta -= $c2_n * $fn; $B_phi += $c3_n * $fn * $inv_s; } my $psi = $theta - (pip2 - $lat); my $X = -$B_theta * cos($psi) - $B_r * sin($psi); my $Y = $B_phi; return ($X != 0 || $Y != 0 ? Math::Trig::rad2deg(atan2($Y, $X)) : 0); } 1; __DATA__ # WMM2005COF: 2005.0 WMM-2005 12/8/2004 1 0 -29556.8 0.0 8.0 0.0 1 1 -1671.7 5079.8 10.6 -20.9 2 0 -2340.6 0.0 -15.1 0.0 2 1 3046.9 -2594.7 -7.8 -23.2 2 2 1657.0 -516.7 -0.8 -14.6 3 0 1335.4 0.0 0.4 0.0 3 1 -2305.1 -199.9 -2.6 5.0 3 2 1246.7 269.3 -1.2 -7.0 3 3 674.0 -524.2 -6.5 -0.6 4 0 919.8 0.0 -2.5 0.0 4 1 798.1 281.5 2.8 2.2 4 2 211.3 -226.0 -7.0 1.6 4 3 -379.4 145.8 6.2 5.8 4 4 100.0 -304.7 -3.8 0.1 5 0 -227.4 0.0 -2.8 0.0 5 1 354.6 42.4 0.7 0.0 5 2 208.7 179.8 -3.2 1.7 5 3 -136.5 -123.0 -1.1 2.1 5 4 -168.3 -19.5 0.1 4.8 5 5 -14.1 103.6 -0.8 -1.1 6 0 73.2 0.0 -0.7 0.0 6 1 69.7 -20.3 0.4 -0.6 6 2 76.7 54.7 -0.3 -1.9 6 3 -151.2 63.6 2.3 -0.4 6 4 -14.9 -63.4 -2.1 -0.5 6 5 14.6 -0.1 -0.6 -0.3 6 6 -86.3 50.4 1.4 0.7 7 0 80.1 0.0 0.2 0.0 7 1 -74.5 -61.5 -0.1 0.6 7 2 -1.4 -22.4 -0.3 0.4 7 3 38.5 7.2 1.1 0.2 7 4 12.4 25.4 0.6 0.3 7 5 9.5 11.0 0.5 -0.8 7 6 5.7 -26.4 -0.4 -0.2 7 7 1.8 -5.1 0.6 0.1 8 0 24.9 0.0 0.1 0.0 8 1 7.7 11.2 0.3 -0.2 8 2 -11.6 -21.0 -0.4 0.1 8 3 -6.9 9.6 0.3 0.3 8 4 -18.2 -19.8 -0.3 0.4 8 5 10.0 16.1 0.2 0.1 8 6 9.2 7.7 0.4 -0.2 8 7 -11.6 -12.9 -0.7 0.4 8 8 -5.2 -0.2 0.4 0.4 9 0 5.6 0.0 0.0 0.0 9 1 9.9 -20.1 0.0 0.0 9 2 3.5 12.9 0.0 0.0 9 3 -7.0 12.6 0.0 0.0 9 4 5.1 -6.7 0.0 0.0 9 5 -10.8 -8.1 0.0 0.0 9 6 -1.3 8.0 0.0 0.0 9 7 8.8 2.9 0.0 0.0 9 8 -6.7 -7.9 0.0 0.0 9 9 -9.1 6.0 0.0 0.0 10 0 -2.3 0.0 0.0 0.0 10 1 -6.3 2.4 0.0 0.0 10 2 1.6 0.2 0.0 0.0 10 3 -2.6 4.4 0.0 0.0 10 4 0.0 4.8 0.0 0.0 10 5 3.1 -6.5 0.0 0.0 10 6 0.4 -1.1 0.0 0.0 10 7 2.1 -3.4 0.0 0.0 10 8 3.9 -0.8 0.0 0.0 10 9 -0.1 -2.3 0.0 0.0 10 10 -2.3 -7.9 0.0 0.0 11 0 2.8 0.0 0.0 0.0 11 1 -1.6 0.3 0.0 0.0 11 2 -1.7 1.2 0.0 0.0 11 3 1.7 -0.8 0.0 0.0 11 4 -0.1 -2.5 0.0 0.0 11 5 0.1 0.9 0.0 0.0 11 6 -0.7 -0.6 0.0 0.0 11 7 0.7 -2.7 0.0 0.0 11 8 1.8 -0.9 0.0 0.0 11 9 0.0 -1.3 0.0 0.0 11 10 1.1 -2.0 0.0 0.0 11 11 4.1 -1.2 0.0 0.0 12 0 -2.4 0.0 0.0 0.0 12 1 -0.4 -0.4 0.0 0.0 12 2 0.2 0.3 0.0 0.0 12 3 0.8 2.4 0.0 0.0 12 4 -0.3 -2.6 0.0 0.0 12 5 1.1 0.6 0.0 0.0 12 6 -0.5 0.3 0.0 0.0 12 7 0.4 0.0 0.0 0.0 12 8 -0.3 0.0 0.0 0.0 12 9 -0.3 0.3 0.0 0.0 12 10 -0.1 -0.9 0.0 0.0 12 11 -0.3 -0.4 0.0 0.0 12 12 -0.1 0.8 0.0 0.0 999999999999999999999999999999999999999999999999 999999999999999999999999999999999999999999999999 # WMM2010COF: 2010.0 WMM-2010 11/20/2009 1 0 -29496.6 0.0 11.6 0.0 1 1 -1586.3 4944.4 16.5 -25.9 2 0 -2396.6 0.0 -12.1 0.0 2 1 3026.1 -2707.7 -4.4 -22.5 2 2 1668.6 -576.1 1.9 -11.8 3 0 1340.1 0.0 0.4 0.0 3 1 -2326.2 -160.2 -4.1 7.3 3 2 1231.9 251.9 -2.9 -3.9 3 3 634.0 -536.6 -7.7 -2.6 4 0 912.6 0.0 -1.8 0.0 4 1 808.9 286.4 2.3 1.1 4 2 166.7 -211.2 -8.7 2.7 4 3 -357.1 164.3 4.6 3.9 4 4 89.4 -309.1 -2.1 -0.8 5 0 -230.9 0.0 -1.0 0.0 5 1 357.2 44.6 0.6 0.4 5 2 200.3 188.9 -1.8 1.8 5 3 -141.1 -118.2 -1.0 1.2 5 4 -163.0 0.0 0.9 4.0 5 5 -7.8 100.9 1.0 -0.6 6 0 72.8 0.0 -0.2 0.0 6 1 68.6 -20.8 -0.2 -0.2 6 2 76.0 44.1 -0.1 -2.1 6 3 -141.4 61.5 2.0 -0.4 6 4 -22.8 -66.3 -1.7 -0.6 6 5 13.2 3.1 -0.3 0.5 6 6 -77.9 55.0 1.7 0.9 7 0 80.5 0.0 0.1 0.0 7 1 -75.1 -57.9 -0.1 0.7 7 2 -4.7 -21.1 -0.6 0.3 7 3 45.3 6.5 1.3 -0.1 7 4 13.9 24.9 0.4 -0.1 7 5 10.4 7.0 0.3 -0.8 7 6 1.7 -27.7 -0.7 -0.3 7 7 4.9 -3.3 0.6 0.3 8 0 24.4 0.0 -0.1 0.0 8 1 8.1 11.0 0.1 -0.1 8 2 -14.5 -20.0 -0.6 0.2 8 3 -5.6 11.9 0.2 0.4 8 4 -19.3 -17.4 -0.2 0.4 8 5 11.5 16.7 0.3 0.1 8 6 10.9 7.0 0.3 -0.1 8 7 -14.1 -10.8 -0.6 0.4 8 8 -3.7 1.7 0.2 0.3 9 0 5.4 0.0 0.0 0.0 9 1 9.4 -20.5 -0.1 0.0 9 2 3.4 11.5 0.0 -0.2 9 3 -5.2 12.8 0.3 0.0 9 4 3.1 -7.2 -0.4 -0.1 9 5 -12.4 -7.4 -0.3 0.1 9 6 -0.7 8.0 0.1 0.0 9 7 8.4 2.1 -0.1 -0.2 9 8 -8.5 -6.1 -0.4 0.3 9 9 -10.1 7.0 -0.2 0.2 10 0 -2.0 0.0 0.0 0.0 10 1 -6.3 2.8 0.0 0.1 10 2 0.9 -0.1 -0.1 -0.1 10 3 -1.1 4.7 0.2 0.0 10 4 -0.2 4.4 0.0 -0.1 10 5 2.5 -7.2 -0.1 -0.1 10 6 -0.3 -1.0 -0.2 0.0 10 7 2.2 -3.9 0.0 -0.1 10 8 3.1 -2.0 -0.1 -0.2 10 9 -1.0 -2.0 -0.2 0.0 10 10 -2.8 -8.3 -0.2 -0.1 11 0 3.0 0.0 0.0 0.0 11 1 -1.5 0.2 0.0 0.0 11 2 -2.1 1.7 0.0 0.1 11 3 1.7 -0.6 0.1 0.0 11 4 -0.5 -1.8 0.0 0.1 11 5 0.5 0.9 0.0 0.0 11 6 -0.8 -0.4 0.0 0.1 11 7 0.4 -2.5 0.0 0.0 11 8 1.8 -1.3 0.0 -0.1 11 9 0.1 -2.1 0.0 -0.1 11 10 0.7 -1.9 -0.1 0.0 11 11 3.8 -1.8 0.0 -0.1 12 0 -2.2 0.0 0.0 0.0 12 1 -0.2 -0.9 0.0 0.0 12 2 0.3 0.3 0.1 0.0 12 3 1.0 2.1 0.1 0.0 12 4 -0.6 -2.5 -0.1 0.0 12 5 0.9 0.5 0.0 0.0 12 6 -0.1 0.6 0.0 0.1 12 7 0.5 0.0 0.0 0.0 12 8 -0.4 0.1 0.0 0.0 12 9 -0.4 0.3 0.0 0.0 12 10 0.2 -0.9 0.0 0.0 12 11 -0.8 -0.2 -0.1 0.0 12 12 0.0 0.9 0.1 0.0 999999999999999999999999999999999999999999999999 999999999999999999999999999999999999999999999999 # WMM2015COF: 2015.0 WMM-2015 12/15/2014 1 0 -29438.5 0.0 10.7 0.0 1 1 -1501.1 4796.2 17.9 -26.8 2 0 -2445.3 0.0 -8.6 0.0 2 1 3012.5 -2845.6 -3.3 -27.1 2 2 1676.6 -642.0 2.4 -13.3 3 0 1351.1 0.0 3.1 0.0 3 1 -2352.3 -115.3 -6.2 8.4 3 2 1225.6 245.0 -0.4 -0.4 3 3 581.9 -538.3 -10.4 2.3 4 0 907.2 0.0 -0.4 0.0 4 1 813.7 283.4 0.8 -0.6 4 2 120.3 -188.6 -9.2 5.3 4 3 -335.0 180.9 4.0 3.0 4 4 70.3 -329.5 -4.2 -5.3 5 0 -232.6 0.0 -0.2 0.0 5 1 360.1 47.4 0.1 0.4 5 2 192.4 196.9 -1.4 1.6 5 3 -141.0 -119.4 0.0 -1.1 5 4 -157.4 16.1 1.3 3.3 5 5 4.3 100.1 3.8 0.1 6 0 69.5 0.0 -0.5 0.0 6 1 67.4 -20.7 -0.2 0.0 6 2 72.8 33.2 -0.6 -2.2 6 3 -129.8 58.8 2.4 -0.7 6 4 -29.0 -66.5 -1.1 0.1 6 5 13.2 7.3 0.3 1.0 6 6 -70.9 62.5 1.5 1.3 7 0 81.6 0.0 0.2 0.0 7 1 -76.1 -54.1 -0.2 0.7 7 2 -6.8 -19.4 -0.4 0.5 7 3 51.9 5.6 1.3 -0.2 7 4 15.0 24.4 0.2 -0.1 7 5 9.3 3.3 -0.4 -0.7 7 6 -2.8 -27.5 -0.9 0.1 7 7 6.7 -2.3 0.3 0.1 8 0 24.0 0.0 0.0 0.0 8 1 8.6 10.2 0.1 -0.3 8 2 -16.9 -18.1 -0.5 0.3 8 3 -3.2 13.2 0.5 0.3 8 4 -20.6 -14.6 -0.2 0.6 8 5 13.3 16.2 0.4 -0.1 8 6 11.7 5.7 0.2 -0.2 8 7 -16.0 -9.1 -0.4 0.3 8 8 -2.0 2.2 0.3 0.0 9 0 5.4 0.0 0.0 0.0 9 1 8.8 -21.6 -0.1 -0.2 9 2 3.1 10.8 -0.1 -0.1 9 3 -3.1 11.7 0.4 -0.2 9 4 0.6 -6.8 -0.5 0.1 9 5 -13.3 -6.9 -0.2 0.1 9 6 -0.1 7.8 0.1 0.0 9 7 8.7 1.0 0.0 -0.2 9 8 -9.1 -3.9 -0.2 0.4 9 9 -10.5 8.5 -0.1 0.3 10 0 -1.9 0.0 0.0 0.0 10 1 -6.5 3.3 0.0 0.1 10 2 0.2 -0.3 -0.1 -0.1 10 3 0.6 4.6 0.3 0.0 10 4 -0.6 4.4 -0.1 0.0 10 5 1.7 -7.9 -0.1 -0.2 10 6 -0.7 -0.6 -0.1 0.1 10 7 2.1 -4.1 0.0 -0.1 10 8 2.3 -2.8 -0.2 -0.2 10 9 -1.8 -1.1 -0.1 0.1 10 10 -3.6 -8.7 -0.2 -0.1 11 0 3.1 0.0 0.0 0.0 11 1 -1.5 -0.1 0.0 0.0 11 2 -2.3 2.1 -0.1 0.1 11 3 2.1 -0.7 0.1 0.0 11 4 -0.9 -1.1 0.0 0.1 11 5 0.6 0.7 0.0 0.0 11 6 -0.7 -0.2 0.0 0.0 11 7 0.2 -2.1 0.0 0.1 11 8 1.7 -1.5 0.0 0.0 11 9 -0.2 -2.5 0.0 -0.1 11 10 0.4 -2.0 -0.1 0.0 11 11 3.5 -2.3 -0.1 -0.1 12 0 -2.0 0.0 0.1 0.0 12 1 -0.3 -1.0 0.0 0.0 12 2 0.4 0.5 0.0 0.0 12 3 1.3 1.8 0.1 -0.1 12 4 -0.9 -2.2 -0.1 0.0 12 5 0.9 0.3 0.0 0.0 12 6 0.1 0.7 0.1 0.0 12 7 0.5 -0.1 0.0 0.0 12 8 -0.4 0.3 0.0 0.0 12 9 -0.4 0.2 0.0 0.0 12 10 0.2 -0.9 0.0 0.0 12 11 -0.9 -0.2 0.0 0.0 12 12 0.0 0.7 0.0 0.0 999999999999999999999999999999999999999999999999 999999999999999999999999999999999999999999999999