#! /usr/bin/env perl
#
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> and
I<2010>.  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
