#!/usr/bin/perl -w
# NAME: solve.pl
# AIM: Given a position, decide the track to take to best arrive at a 'circuit'
# At present only fixed YGIL positions used, but added reading of apt.dat.gz to 
# find an ICAO, and should be able to do it for any airport
# Needs LOTS of work to be fully 'generic'
# Also an example of using ImageMagik 'convert' to generate a GIF image of the circuit.
# Presently ONLY a DEBUG ON version!
# Also note gzip reports a bloken pipe when closed before all lines read, but
# the warning can be ignored.
use strict;
use warnings;
use File::Basename;  # split path ($name,$dir,$ext) = fileparse($file [, qr/\.[^.]*/] )
use Cwd;
use IO::Socket;
use Term::ReadKey;
use Time::HiRes qw( usleep gettimeofday tv_interval );
use Math::Trig;
my $perl_dir = 'C:\GTools\perl';
unshift(@INC, $perl_dir);
require 'lib_utils.pl' or die "Unable to load 'lib_utils.pl'! Check location and \@INC content.\n";
require 'fg_wsg84.pl' or die "Unable to load fg_wsg84.pl ...\n";
# NOT USED require "Bucket2.pm" or die "Unable to load Bucket2.pm ...\n";
# ===============

# log file stuff
our ($LF);
my $pgmname = $0;
if ($pgmname =~ /(\\|\/)/) {
    my @tmpsp = split(/(\\|\/)/,$pgmname);
    $pgmname = $tmpsp[-1];
}
my $outfile = $perl_dir."\\temp.$pgmname.txt";
open_log($outfile);

### constants
my $VER = "0.0.2 2014-11-29";
###my $VER = "0.0.1 2010-09-11";
my $SGD_PI = 3.1415926535;
my $SGD_DEGREES_TO_RADIANS = $SGD_PI / 180.0;
my $SGD_RADIANS_TO_DEGREES = 180.0 / $SGD_PI;
my $DEF_GS = 3;
my $ATAN3 = atan( $DEF_GS * $SGD_DEGREES_TO_RADIANS );
# /** Feet to Meters */
my $SG_FEET_TO_METER = 0.3048;
# /** Meters to Feet */
my $SG_METER_TO_FEET = 3.28083989501312335958;
my $SG_NM_TO_METER = 1852;
my $SG_METER_TO_NM = 0.0005399568034557235;

# user variables
my $load_log = 0;
my $in_file = '';
my $bad_latlon = 200;
my $in_lat = $bad_latlon;
my $in_lon = $bad_latlon;
my $graf_file = "tempgraf.gif";
my $graf_bat = $perl_dir."\\tempgraf.bat";
my $apt_icao = 'YGIL';

my $debug_on = 1;
my $def_latlon = "-31.6,148.6";
my $def_lat2 = "-31.703";
my $def_lon2 = "148.6744";
my $def_lat3 = "-31.8145";
my $def_lon3 = "148.6374";
my $def_lat4 = "-31.756";
my $def_lon4 = "148.559";
my $def_lat5 = "-31.698961";
my $def_lon5 = "148.627853";
my $def_lat6 = "-31.701298";
my $def_lon6 = "148.6045";

my $stand_glide_degs = 3; # degrees
my $stand_patt_alt = 1000; # feet
my $stand_cross_nm = 2.1; # nm, but this will depend on the aircraft

### program variables
my @warnings = ();
my $cwd = cwd();
my $os = $^O;
my @my_airport = ();

my $a_gil_lat = -31.697287500;
my $a_gil_lon = 148.636942500;
my $a_dub_lat = -32.2174865;
my $a_dub_lon = 148.57727;

# rough Gil circuit
my $tl_lat = -31.684063;
my $tl_lon = 148.614120;
my $bl_lat = -31.723495;
my $bl_lon = 148.633003;
my $br_lat = -31.716778;
my $br_lon = 148.666992;
my $tr_lat = -31.672960;
my $tr_lon = 148.649139;
my $use_pattern = 1; # adjust the above values to the computed circuit
my $add_text_count = 1; # add text count
my $try_dash_line = 1;
my $switch_circuit = 0; # try the OTHER circuit 15 (def = 33)

# my $YSDU_metar = "2011/07/22 14:00 YSDU 221400Z AUTO 15010KT 9999 // NCD 06/04 Q1020";
my $WIND_DIR_ABB = 'SSE';
my $WIND_DIR_DEG = 150;
my $WIND_DIR_ENG = 'South/Southeast';
my $WIND_KTS     = 8; # kt 9.2 mph

# RUNWAY ARRAY OFFSETS
my $RW_LEN = 0;
my $RW_HDG = 1;
my $RW_REV = 2;
my $RW_TT1 = 3;
my $RW_TT2 = 4;
my $RW_CLAT = 5;
my $RW_CLON = 6;
my $RW_LLAT = 7;
my $RW_LLON = 8;
my $RW_RLAT = 9;
my $RW_RLON = 10;
my $RW_DONE = 11;
#                 Len    Hdg   Rev  Title  RTit Ctr Lat    Ctr Lon
#                 0      1     2    3     4     5          6           7  8  9  10 11
my @gil_patt = ();
##my @gil_rwys = ( [4204,  162.0, 0, '15', '33', -31.696928, 148.636404, 0, 0, 0, 0, 0 ] );
my @gil_rwys = ( [3984,  162.22, 0, '15', '33', -31.69656323, 148.6363057, 0, 0, 0, 0, 0 ] );

#my @gil_navs = ( ["", 0 ] );
my @gil_navs = ();
#my @gil_rwys = ( [162.0, 4204], [93.0, 1902] );
my @dub_patt = ( [ ] );
my @dub_rwys = ( [5600, 53.61, 0, '05', '23', -32.218265, 148.576145, 0, 0, 0, 0, 0 ] );
my @dub_navs = ( ["VOR", 114.4], ["NDB", 251] );

my $OL_LAT = 0;
my $OL_LON = 1;
my $OL_NAV = 2;
my $OL_RWY = 3;
my $OL_PAT = 4;
my %apt_locations = (
    # key      0           1           2           3           4
    # ICAO => [ Center LAT, LON        NAVAIDS     RUNWAYS     CIRCUIT ]
    'YGIL' => [$a_gil_lat, $a_gil_lon, \@gil_navs, \@gil_rwys, \@gil_patt ],
    'YSDU' => [$a_dub_lat, $a_dub_lon, \@dub_navs, \@dub_rwys, \@dub_patt ]
    );

sub get_locations() { return \%apt_locations; }

# DEBUG
my $dbg_01 = 0; # show pattern values
my $dbg_02 = 0;

sub show_warnings($) {
    my ($val) = @_;
    if (@warnings) {
        prt( "\nGot ".scalar @warnings." WARNINGS...\n" );
        foreach my $itm (@warnings) {
           prt("$itm\n");
        }
        prt("\n");
    } else {
        ###prt( "\nNo warnings issued.\n\n" );
    }
}

sub pgm_exit($$) {
    my ($val,$msg) = @_;
    if (length($msg)) {
        $msg .= "\n" if (!($msg =~ /\n$/));
        prt($msg);
    }
    show_warnings($val);
    close_log($outfile,$load_log);
    exit($val);
}


sub prtw($) {
   my ($tx) = shift;
   $tx =~ s/\n$//;
   prt("$tx\n");
   push(@warnings,$tx);
}


sub in_world_range($$) {
    my ($lat,$lon) = @_;
    if (($lat < -90) ||
        ($lat >  90) ||
        ($lon < -180) ||
        ($lon > 180) ) {
        return 0;
    }
    return 1;
}

# from : http://compsci.ca/v3/viewtopic.php?t=6034
# if a point P is inside triangle ABC, then 
# Area PAB+Area PBC +Area PAC=Area ABC 
# notice that if P is on the edge of AB, BC, or CA, the above hold. 
# But effectively, one of the area PAB, PBC, PAC is 0 (so just make sure you check that). 
# if P is outside, the above equality does NOT hold... 
# for example, if A=(x1,y1) B=(x2,y2), C=(x3,y3) 
# Area= abs(x1*y2 + x2*y3 + x3*y1 - x1*y3 - x3*y2 -x2*y1 ) / 2 
sub Triangle_Area($$$$$$) {
    my ($x1,$y1,$x2,$y2,$x3,$y3) = @_;
    return abs($x1*$y2 + $x2*$y3 + $x3*$y1 - $x1*$y3 - $x3*$y2 - $x2*$y1) / 2;
}
# Lots of stuff at : http://xlinux.nist.gov/dads/
# from : http://www.blackpawn.com/texts/pointinpoly/default.html
# Point in triangle test
# A common way to check if a point is in a triangle is to find 
# the vectors connecting the point to each of the triangle's 
# three vertices and sum the angles between those vectors. 
# If the sum of the angles is 2*pi then the point is inside the triangle, otherwise it is not.
# CrossProduct - A x B = A B sin(&)
# http://tutorial.math.lamar.edu/Classes/CalcII/CrossProduct.aspx
# Two vector: A = (a1,a2,a3), B = (b1,b2,b3)
# A X B = (a2b3 - a3b2, a3b1 - a1b3, a1b2 - a2b1)
# from : http://schools-wikipedia.org/wp/c/Cross_product.htm
# Consider two vectors, a = (1,2,3) and b = (4,5,6). 
# The cross product a × b is 
# a × b = (1,2,3) × (4,5,6) = ((2 × 6 - 3 × 5),-(1 × 6 - 3 × 4),+(1 × 5 - 2 × 4)) = (-3,6,-3).
# Consider two vectors, a = (3,0,0) and b = (0,2,0). The cross product a × b is
# a × b = (3,0,0) × (0,2,0) = ((0 × 0 - 0 × 2), (0 × 0 - 3 × 0), (3 × 2 - 0 × 0)) = (0,0,6).
# from : http://msdn.microsoft.com/en-us/library/system.windows.vector.crossproduct.aspx
# double CrossProduct( Vector & v1, Vector & v2 ) {
#    return (v1.X * v2.Y) - (v1.Y * v2.X); }
# If the value of a1b2 - a2b1 is positive, then it points OUT of the page; 
# if its value is negative, then it points INTO the page.
# from : http://web.mit.edu/wwmath/vectorc/3d/crossp.html
#(v1, v2, v3) × (w1, w2, w3) = ( v2 w3 - v3 w2, v3 w1 - v1 w3, v1 w2 - v2 w1 ) 
# from : http://www.blackpawn.com/texts/pointinpoly/default.html
# function SameSide(p1,p2, a,b)
#    cp1 = CrossProduct(b-a, p1-a)
#    cp2 = CrossProduct(b-a, p2-a)
#    if DotProduct(cp1, cp2) >= 0 then return true
#    else return false
# function PointInTriangle(p, a,b,c)
#    if SameSide(p,a, b,c) and SameSide(p,b, a,c)
#        and SameSide(p,c, a,b) then return true
#    else return false
# Barycentric Technique
# With that in mind we can now describe any point on the plane as 
#    P = A + u * (C - A) + v * (B - A)
#Notice now that if u or v < 0 then we've walked in the wrong direction and must be outside the triangle. Also if u or v > 1 then we've walked too far in a direction and are outside the triangle. Finally if u + v > 1 then we've crossed the edge BC again leaving the triangle. 
#Given u and v we can easily calculate the point P with the above equation, but how can we go in the reverse direction and calculate u and v from a given point P? Time for some math! 
#    P = A + u * (C - A) + v * (B - A)       // Original equation
#    (P - A) = u * (C - A) + v * (B - A)     // Subtract A from both sides
#    v2 = u * v0 + v * v1                    // Substitute v0, v1, v2 for less writing
#    // We have two unknowns (u and v) so we need two equations to solve
#    // for them.  Dot both sides by v0 to get one and dot both sides by
#    // v1 to get a second.
#    (v2) . v0 = (u * v0 + v * v1) . v0
#    (v2) . v1 = (u * v0 + v * v1) . v1
#    // Distribute v0 and v1
#    v2 . v0 = u * (v0 . v0) + v * (v1 . v0)
#    v2 . v1 = u * (v0 . v1) + v * (v1 . v1)
#    // Now we have two equations and two unknowns and can solve one 
#    // equation for one variable and substitute into the other.  Or
#    // if you're lazy like me, fire up Mathematica and save yourself
#    // some handwriting.
#    Solve[v2.v0 == {u(v0.v0) + v(v1.v0), v2.v1 == u(v0.v1) + v(v1.v1)}, {u, v}]
#    u = ((v1.v1)(v2.v0)-(v1.v0)(v2.v1)) / ((v0.v0)(v1.v1) - (v0.v1)(v1.v0))
#    v = ((v0.v0)(v2.v1)-(v0.v1)(v2.v0)) / ((v0.v0)(v1.v1) - (v0.v1)(v1.v0))
# OR
# // Compute vectors        
# v0 = C - A
# v1 = B - A
# v2 = P - A
#// Compute dot products
# dot00 = dot(v0, v0)
# dot01 = dot(v0, v1)
# dot02 = dot(v0, v2)
# dot11 = dot(v1, v1)
# dot12 = dot(v1, v2)
#// Compute barycentric coordinates
#invDenom = 1 / (dot00 * dot11 - dot01 * dot01)
# u = (dot11 * dot02 - dot01 * dot12) * invDenom
# v = (dot00 * dot12 - dot01 * dot02) * invDenom
# // Check if point is in triangle
# return (u > 0) && (v > 0) && (u + v < 1)
# from : http://mathforum.org/library/drmath/view/54735.html
#How can I know whether a point P(x,y) is in a triangle with vertices
#P1(x1,y1), P2(x2,y2), and P3(x3,y3)?
#There are two possibilities:
#     1. P is INSIDE or ON,
#     2. P is OUTSIDE triangle P1P2P3.
#If you draw a sketch and consider triangles,
#     PP1P2, PP2P3, PP3P1
#You will notice the following relationship in their areas.
#If P is ON or INSIDE,
#     A(PP1P2) + A(PP2P3) + A(PP3P1) = A(P1P2P3)
#If P is OUTSIDE then,
#     A(PP1P2) + A(PP2P3) + A(PP3P1) > A(P1P2P3)
#Area of a triangle can be found using the determinant:
#     A = abs(1/2 * |x1  y1  1| )
#                   |x2  y2  1|
#                   |x3  y3  1|
#You can write a function in C++ that accepts the three points
#(as Classes) or their coordinates and returns the area.
#Finally if you have to distinguish between ON and INSIDE, put the point
#in the equations of the three lines that form the sides. If one of them
#is satisfied, P is ON a side. If two are satisfied, P is a vertex of
#the triangle.
sub Point_Inside_Triange($$$$$$$$) {
    my ($px,$py,$x1,$y1,$x2,$y2,$x3,$y3) = @_;
    my $a1 = Triangle_Area($px,$py,$x1,$y1,$x2,$y2);
    my $a2 = Triangle_Area($px,$py,$x2,$y2,$x3,$y3);
    my $a3 = Triangle_Area($px,$py,$x3,$y3,$x1,$y1);
    my $at = Triangle_Area($x1,$y1,$x2,$y2,$x3,$y3);
    my $sum = $a1 + $a2 + $a3;
    my $diff = abs($at - $sum);
#    my $minlat = 200;
#    my $maxlat = -200;
#    my $minlon = 200;
#    my $maxlon = -200;
#    $minlat = $x1 if ($x1 < $minlat);
#    $minlat = $x2 if ($x2 < $minlat);
#    $minlat = $x3 if ($x3 < $minlat);
#    $maxlat = $x1 if ($x1 > $maxlat);
#    $maxlat = $x2 if ($x2 > $maxlat);
#    $maxlat = $x3 if ($x3 > $maxlat);
#    $minlon = $y1 if ($y1 < $minlon);
#    $minlon = $y2 if ($y2 < $minlon);
#    $minlon = $y3 if ($y3 < $minlon);
#    $maxlon = $y1 if ($y1 > $maxlon);
#    $maxlon = $y2 if ($y2 > $maxlon);
#    $maxlon = $y3 if ($y3 > $maxlon);
#    prt("Bound: TL $maxlat,$minlon, BR $minlat,$maxlon\n");
#    if (($px >= $minlat)&&($px <= $maxlat)&&($py >= $minlon)&&($py <= $maxlon)) {
#        prt("Point $px,$py is in bounding box!\n");
#    } else {
#        prt("Point $px,$py is NOT in bounding box!\n");
#    }
#    prt("Areas: if (($a1 + $a2 + $a3) > $at) \n");
#    prt("Total $sum, diff $diff\n");
    #if (($a1 + $a2 + $a3) > $at) { return 0; } # THIS FAILED in some cases
    return 1 if ($diff < 1.0E-010); # take this SMALL value as EQUAL !!! 
    return 0;
}

# ============================================================ #
# SimGear Services, rendered in perl
# ============================================================ #
# dot(const SGVec3<T>& v1, const SGVec3<T>& v2)
# { return v1(0)*v2(0) + v1(1)*v2(1) + v1(2)*v2(2); }
# Given 2 Vectors3, return the dot product
sub scalar_dot_product($$) {
    my ($rv1,$rv2) = @_;
    return ${$rv1}[0] * ${$rv2}[0] + ${$rv1}[1] * ${$rv2}[1] + ${$rv1}[2] * ${$rv2}[2];
}

# In dimension 2, the dot product of vectors [a,b] and [c,d] is ac + bd
sub scalar_dot_product2($$$$) {
    my ($v1x,$v1y,$v2x,$v2y) = @_;
    return ($v1x * $v2x) + ($v1y * $v2y);
}

# The euclidean norm of the vector, that is what most people call length
# norm(const SGVec3<T>& v)
# { return sqrt(dot(v, v)); }
# Given a Vector3, return length
sub norm_vector_length($) {
    my ($rv) = @_;
    return sqrt(scalar_dot_product($rv, $rv));
}

sub norm_vector_length2($$) {
    my ($vx,$vy) = @_;
    return sqrt(scalar_dot_product2($vx, $vy, $vx, $vy));
}


sub process_in_file($) {
    my ($inf) = @_;
    if (! open INF, "<$inf") {
        pgm_exit(1,"ERROR: Unable to open file [$inf]\n"); 
    }
    my @lines = <INF>;
    close INF;
    my $lncnt = scalar @lines;
    prt("Processing $lncnt lines, from [$inf]...\n");
    my ($line,$inc,$lnn);
    $lnn = 0;
    foreach $line (@lines) {
        chomp $line;
        $lnn++;
        if ($line =~ /\s*#\s*include\s+(.+)$/) {
            $inc = $1;
            prt("$lnn: $inc\n");
        }
    }
}

sub set_int_stg($) {
    my $r = shift;
    ${$r} =  int(${$r} + 0.5);
}

sub set_int_dist_stg5($) {
    my $r = shift;
    set_int_stg($r);
    my $r5 = sprintf("%5d",${$r});
    ${$r} = $r5;
}

sub set_hdg_stg($) {
    my ($rh) = @_;
    my $hdg = ${$rh};
    $hdg = 360 if ($hdg == 0); # always replace 000 with 360 ;=))
    $hdg = sprintf("%03d",int($hdg+0.5));
    ${$rh} = $hdg;
}

sub set_lat_stg($) {
    my ($rl) = @_;
    ${$rl} = sprintf("%2.7f",${$rl});
}
sub set_lon_stg($) {
    my ($rl) = @_;
    ${$rl} = sprintf("%3.7f",${$rl});
}

sub set_decimal1_stg($) {
    my $r = shift;
    ${$r} =  int((${$r} + 0.05) * 10) / 10;
    ${$r} = "0.0" if (${$r} == 0);
    ${$r} .= ".0" if !(${$r} =~ /\./);
}

sub get_dist_stg_nm($) {
    my ($dist) = @_;
    my $nm = $dist * $SG_METER_TO_NM;
    set_decimal1_stg(\$nm);
    $nm .= "nm";
    return $nm;
}

sub normalised_hdg($) {
    my $hdg = shift;
    $hdg += 360 if ($hdg < 0);
    $hdg -= 360 if ($hdg >= 360);
    return $hdg;
}

sub set_hdg_stg3($) {
    my $r = shift;
    set_int_stg($r);
    my $r3 = sprintf("%3d",${$r});
    ${$r} = $r3;
}

# ============================================

sub show_distance_heading($$$$) {
    my ($lat1,$lon1,$lat2,$lon2) = @_;
    my ($az1,$az2,$dist);
    fg_geo_inverse_wgs_84 ($lat1,$lon1,$lat2,$lon2,\$az1,\$az2,\$dist);
    $dist = get_dist_stg_nm($dist);
    set_hdg_stg(\$az1);
    prt("Is $dist, on heading $az1\n");
}

sub show_rw_patt($$) {
    my ($key,$rpatts) = @_;
    my $cnt = scalar @{$rpatts};
    prt("Display of $cnt patterns/circuits for $key...\n");
    my ($i,$lat1,$lon1,$lat2,$lon2,$j);
    for ($i = 0; $i < $cnt; $i++) {
        for ($j = 0; $j < 8; $j += 2) {
            $lat1 = ${$rpatts}[$i][$j+0];
            $lon1 = ${$rpatts}[$i][$j+1];
            if ($j == 6) {
                $lat2 = ${$rpatts}[$i][0];
                $lon2 = ${$rpatts}[$i][1];
            } else {
                $lat2 = ${$rpatts}[$i][$j+2];
                $lon2 = ${$rpatts}[$i][$j+3];
            }
            prt("$i:$j: $lat1,$lon1  $lat2,$lon2\n");
            show_distance_heading($lat1,$lon1,$lat2,$lon2);
        }
    }
}

sub set_runway_ends_and_patt($$$$) {
    my ($rrwys,$i,$key,$rpatts) = @_;
    # set ENDS of runway
    my $rlen = ${$rrwys}[$i][$RW_LEN];
    my $rhdg = ${$rrwys}[$i][$RW_HDG];
    my $clat = ${$rrwys}[$i][$RW_CLAT];
    my $clon = ${$rrwys}[$i][$RW_CLON];
    my $rty1 = ${$rrwys}[$i][$RW_TT1];
    my $rty2 = ${$rrwys}[$i][$RW_TT2];
    my $rwlen2 = ($rlen * $SG_FEET_TO_METER) / 2;
    my ($elat1,$elon1,$eaz1,$elat2,$elon2,$eaz2);
    my $hdgr = $rhdg + 180;
    $hdgr -= 360 if ($hdgr >= 360);
    ${$rrwys}[$i][$RW_REV] = $hdgr;

    fg_geo_direct_wgs_84( $clat, $clon, $rhdg, $rwlen2, \$elat1, \$elon1, \$eaz1 );
    fg_geo_direct_wgs_84( $clat, $clon, $hdgr, $rwlen2, \$elat2, \$elon2, \$eaz2 );
    ${$rrwys}[$i][$RW_LLAT] = $elat1;
    ${$rrwys}[$i][$RW_LLON] = $elon1;
    ${$rrwys}[$i][$RW_RLAT] = $elat2;
    ${$rrwys}[$i][$RW_RLON] = $elon2;
    ${$rrwys}[$i][$RW_DONE] = $i + 1;

    my ($az1,$az2,$dist);
    fg_geo_inverse_wgs_84 ($elat1,$elon1,$elat2,$elon2,\$az1,\$az2,\$dist);
    $dist = $dist * $SG_METER_TO_FEET;
    set_int_stg(\$az1);
    set_int_stg(\$az2);
    set_int_stg(\$dist);
    # init: YSDU: 23: -32.2136987804606,148.583432501246 05: -32.2228307960945,148.568856770273 234 5600 54 vs 53.61 5600
    # init: YGIL: 33: -31.7024233216057,148.638492502638 15: -31.6914326394609,148.634315743548 342 4204 162 vs 162 4204
    #prt("init: $key: $rty2: $elat1,$elon1 $az1 $rty1: $elat2,$elon2 $az1 $dist $az2 vs $rhdg $rlen\n");
    prt("init:$i: $key: $rty2:$az1: $elat1,$elon1 $rty1:$az2: $elat2,$elon2\n");
    # We have the RUNWAY ends - now extend out to first turn to crosswind leg, and turn to final
    # but by how MUCH - ok decide from runway end, out to where it is a 3 degree glide from 1000 feet
    $dist = ($stand_patt_alt * $SG_FEET_TO_METER) / tan($stand_glide_degs * $SGD_DEGREES_TO_RADIANS);
    my ($plat11,$plon11,$plat12,$plon12,$plat13,$plon13,$paz1);
    my ($plat21,$plon21,$plat22,$plon22,$plat23,$plon23,$paz2);
    my ($hdg1L,$hdg1R,$crossd);
    fg_geo_direct_wgs_84( $clat, $clon, $rhdg, $rwlen2+$dist, \$plat11, \$plon11, \$paz1 );
    fg_geo_direct_wgs_84( $clat, $clon, $hdgr, $rwlen2+$dist, \$plat21, \$plon21, \$paz2 );
    $hdg1L = normalised_hdg($rhdg - 90);
    $hdg1R = normalised_hdg($rhdg + 90);
    $crossd = $stand_cross_nm * $SG_NM_TO_METER;
    # ON $rhdg to $elat1, $elon1 to ... turn point, go LEFT and to get NEXT points, this end
    fg_geo_direct_wgs_84( $plat11, $plon11, $hdg1L, $crossd, \$plat12, \$plon12, \$paz1 );
    fg_geo_direct_wgs_84( $plat21, $plon21, $hdg1L, $crossd, \$plat13, \$plon13, \$paz1 );

    # from the turn point, go LEFT and RIGHT to get NEXT points, this other end
    fg_geo_direct_wgs_84( $plat21, $plon21, $hdg1R, $crossd, \$plat22, \$plon22, \$paz2 );
    fg_geo_direct_wgs_84( $plat11, $plon11, $hdg1R, $crossd, \$plat23, \$plon23, \$paz2 );

    if ($use_pattern && ($key eq 'YGIL')) {
        if ($switch_circuit) {
            # At YGIL, this is a 15 circuit (the prevailing wind! SSE...
            $tl_lat = $plat12;
            $tl_lon = $plon12;
            $bl_lat = $plat13;
            $bl_lon = $plon13;
            $br_lat = $plat21;
            $br_lon = $plon21;
            $tr_lat = $plat11;
            $tr_lon = $plon11;
        } else {
            # At YGIL, this is a 33 circuit
            $tl_lat = $plat22; #-31.684063;
            $tl_lon = $plon22; #148.614120;
            $bl_lat = $plat23; #-31.723495;
            $bl_lon = $plon23; #148.633003;
            $br_lat = $plat11; #-31.716778;
            $br_lon = $plon11; #148.666992;
            $tr_lat = $plat21; #-31.672960;
            $tr_lon = $plon21; #148.649139;
        }
        prt("Set pattern as the rectangle...\n");
    }

    if ($dbg_01) {
        # now we have 4 points, either side of the runway
        prt("On $rhdg, at $plat11,$plon11 turn $hdg1L to $plat12,$plon12\n");
        show_distance_heading($plat11,$plon11,$plat12,$plon12);
        # This is the LONG downwind side 12 to 13
        prt("On $hdg1L at $plat12,$plon12, turn $hdgr to $plat13,$plon13\n");
        show_distance_heading($plat12,$plon12,$plat13,$plon13);
        prt("On $hdgr at $plat13,$plon13 turn $hdg1R to $plat21,$plon21\n");
        show_distance_heading($plat13,$plon13,$plat21,$plon21);
        prt("On $hdg1R at $plat21,$plon21 turn $rhdg to $elat1,$elon1\n");
        show_distance_heading($plat21,$plon21,$elat1,$elon2);
        prt("\n");
        #E.G. for YGIL - TO 15
        #On 162, at -31.7523059488988,148.65746239832 turn 72 to -31.7414611993009,148.696497359091
        #Is 2.1nm, on heading  72
        #On 72 at -31.7414611993009,148.696497359091, turn 342 to -31.630701184372,148.654359238954
        #Is 7.0nm, on heading 342
        #On 342 at -31.630701184372,148.654359238954 turn 252 to -31.6415460967825,148.615370603846
        #Is 2.1nm, on heading 252
        #On 252 at -31.6415460967825,148.615370603846 turn 162 to -31.7024233216057,148.638492502638
        #Is 3.8nm, on heading 165

        prt("On $hdgr at $plat21,$plon21 turn $hdg1R to $plat22,$plon22\n");
        show_distance_heading($plat21,$plon21,$plat22,$plon22);
        # This is the LONG downwind side 22 to 23
        prt("On $hdg1R at $plat22,$plon22 turn $rhdg to $plat23,$plon23\n");
        show_distance_heading($plat22,$plon22,$plat23,$plon23);
        prt("On $rhdg at $plat23,$plon23 turn $hdg1L to $plat11,$plon11\n");
        show_distance_heading($plat23,$plon23,$plat11,$plon11);
        prt("On $hdg1L at $plat11,$plon11 turn $hdgr to $elat2,$elon2\n");
        show_distance_heading($plat11,$plon11,$elat2,$elon2);
        prt("\n");
        #E.G. for YGIL, TO 33
        #On 342 at -31.6415460967825,148.615370603846 turn 252 to -31.6523790808638,148.576372922039
        #Is 2.1nm, on heading 252
        #On 252 at -31.6523790808638,148.576372922039 turn 162 to -31.7631387187979,148.618418340898
        #Is 7.0nm, on heading 162
        #On 162 at -31.7631387187979,148.618418340898 turn 72 to -31.7523059488988,148.65746239832
        #Is 2.1nm, on heading  72
        #On 72 at -31.7523059488988,148.65746239832 turn 342 to -31.6914326394609,148.634315743548
        #Is 3.8nm, on heading 342
    }
    @{$rpatts} = ();
    # add notional RIGHT side circuit first
    push(@{$rpatts}, [$plat11,$plon11,$plat12,$plon12,$plat13,$plon13,$plat21,$plon21,$clat,$clon,$rlen,$rhdg]);
    # then notional LEFT size circuit
    push(@{$rpatts}, [$plat21,$plon21,$plat22,$plon22,$plat23,$plon23,$plat11,$plon11,$clat,$clon,$rlen,$hdgr]);

}

sub init_runway_array() {
    my $rl = get_locations();
    my ($key,$i,$cnt,$rrwys,$rpatts);
    foreach $key (keys %{$rl}) {
        $rrwys = ${$rl}{$key}[$OL_RWY];
        $rpatts = ${$rl}{$key}[$OL_PAT];
        $cnt = scalar @{$rrwys};
        for ($i = 0; $i < $cnt; $i++) {
            set_runway_ends_and_patt($rrwys,$i,$key,$rpatts);
        }
    }
    if ($dbg_02) {
        foreach $key (keys %{$rl}) {
            $rpatts = ${$rl}{$key}[$OL_PAT];
            show_rw_patt($key,$rpatts);
        }
    }
    # pgm_exit(1,"Temp exit");
}

sub get_runways_and_pattern($$) {
    my ($rh,$key) = @_;
    my $rl = get_locations();
    my ($rrwys,$rpatts);
    if (defined ${$rl}{$key}) {
        $rrwys = ${$rl}{$key}[$OL_RWY];
        $rpatts = ${$rl}{$key}[$OL_PAT];
        ${$rh}{'runways'} = $rrwys;
        ${$rh}{'pattern'} = $rpatts;
        ${$rh}{'airport'} = $key;
    } else {
        pgm_exit(1,"ERROR: Key [$key] NOT in locations!\n");
    }
}

# ============================================

sub set_circuit_values($$) {
    my ($rh,$show) = @_;
    my ($az1,$az2,$dist);
    if ($show) {
        #prt("Set, show circuit...\nTL ".${$rh}{'tl_lat'}.",".${$rh}{'tl_lon'}."\nBL ".
        #    ${$rh}{'bl_lat'}.",".${$rh}{'bl_lon'}."\nBR ".
        #    ${$rh}{'br_lat'}.",".${$rh}{'br_lon'}."\nTR ".
        #    ${$rh}{'tr_lat'}.",".${$rh}{'tr_lon'}."\n");
        my ($tllat,$tllon,$bllat,$bllon,$brlat,$brlon,$trlat,$trlon);
        $tllat = ${$rh}{'tl_lat'};
        $tllon = ${$rh}{'tl_lon'};
        $bllat = ${$rh}{'bl_lat'};
        $bllon = ${$rh}{'bl_lon'};
        $brlat = ${$rh}{'br_lat'};
        $brlon = ${$rh}{'br_lon'};
        $trlat = ${$rh}{'tr_lat'};
        $trlon = ${$rh}{'tr_lon'};
        set_lat_stg(\$tllat);
        set_lat_stg(\$bllat);
        set_lat_stg(\$brlat);
        set_lat_stg(\$trlat);
        set_lon_stg(\$tllon);
        set_lon_stg(\$bllon);
        set_lon_stg(\$brlon);
        set_lon_stg(\$trlon);
        prt("Set, show circuit...\nTL $tllat,$tllon\nBL ".
            "$bllat,$bllon\nBR ".
            "$brlat,$brlon\nTR ".
            "$trlat,$trlon\n");
    }
    # anti-clockwise around circuit
    fg_geo_inverse_wgs_84 (${$rh}{'tl_lat'},${$rh}{'tl_lon'},${$rh}{'bl_lat'},${$rh}{'bl_lon'},\$az1,\$az2,\$dist);
    ${$rh}{'l1_az1'} = $az1;
    ${$rh}{'l1_az2'} = $az2;
    ${$rh}{'l1_dist'} = $dist;
    if ($show) {
        set_int_dist_stg5(\$dist);
        set_hdg_stg3(\$az1);
        prt("l1 $dist m, on $az1 (tl2bl)\n");
    }
    fg_geo_inverse_wgs_84 (${$rh}{'bl_lat'},${$rh}{'bl_lon'},${$rh}{'br_lat'},${$rh}{'br_lon'},\$az1,\$az2,\$dist);
    ${$rh}{'l2_az1'} = $az1;
    ${$rh}{'l2_az2'} = $az2;
    ${$rh}{'l2_dist'} = $dist;
    if ($show) {
        set_int_dist_stg5(\$dist);
        set_hdg_stg3(\$az1);
        prt("l2 $dist m, on $az1 (bl2br)\n");
    }
    fg_geo_inverse_wgs_84 (${$rh}{'br_lat'},${$rh}{'br_lon'},${$rh}{'tr_lat'},${$rh}{'tr_lon'},\$az1,\$az2,\$dist);
    ${$rh}{'l3_az1'} = $az1;
    ${$rh}{'l3_az2'} = $az2;
    ${$rh}{'l3_dist'} = $dist;
    if ($show) {
        set_int_dist_stg5(\$dist);
        set_hdg_stg3(\$az1);
        prt("l3 $dist m, on $az1 (br2tr)\n");
    }
    fg_geo_inverse_wgs_84 (${$rh}{'tr_lat'},${$rh}{'tr_lon'},${$rh}{'tl_lat'},${$rh}{'tl_lon'},\$az1,\$az2,\$dist);
    ${$rh}{'l4_az1'} = $az1;
    ${$rh}{'l4_az2'} = $az2;
    ${$rh}{'l4_dist'} = $dist;
    if ($show) {
        set_int_dist_stg5(\$dist);
        set_hdg_stg3(\$az1);
        prt("l4 $dist m, on $az1 (tr2tl)\n");
    }
}

sub set_most_values_plus($$$$$$) {
    my ($rh,$show,$u_lat,$u_lon,$t_lat,$t_lon) = @_;
    if (!defined ${$rh}{'user_points'}) {
        ${$rh}{'user_points'} = [];
    }
    my $ru = ${$rh}{'user_points'};
    push(@{$ru}, [$u_lat,$u_lon,$t_lat,$t_lon]);
}

sub set_min_max($$$$$$) {
    my ($rmaxlat,$rminlat,$rmaxlon,$rminlon,$lat,$lon) = @_;
    ${$rmaxlat} = $lat if ($lat > ${$rmaxlat});
    ${$rminlat} = $lat if ($lat < ${$rminlat});
    ${$rmaxlon} = $lon if ($lon > ${$rmaxlon});
    ${$rminlon} = $lon if ($lon < ${$rminlon});
}

sub add_img_arrow($$$$) {
    my ($rh,$clat,$clon,$az) = @_;
    my ($w_ind1,$h_ind1);
    my ($maxlat,$minlat,$maxlon,$minlon);
    my ($w_dpp,$h_dpp);
    my ($sqwid,$sqhgt,$rmsg);
    #my $adj = 20; # works for some, not others???
    my $adj = 0;
    $rmsg   = ${$rh}{'rmsg'};
    $maxlat = ${$rh}{'max_lat'};
    $minlat = ${$rh}{'min_lat'};
    $maxlon = ${$rh}{'max_lon'};
    $minlon = ${$rh}{'min_lon'};
    $w_dpp  = ${$rh}{'w_dpp'};
    $h_dpp  = ${$rh}{'h_dpp'};
    $sqwid  = ${$rh}{'sq_wid_adj'};
    $sqhgt  = ${$rh}{'sq_hgt_adj'};

    $w_ind1 = int((($clon - $minlon) * $w_dpp) + 0.5); # get degrees/pixels from left edge
    $h_ind1 = int((($clat - $minlat) * $h_dpp) + 0.5); # get degrees/pixels from bottom edge
    $h_ind1 = $sqhgt - $h_ind1;
    $h_ind1 += 1 if ($h_ind1 == 0);
    $w_ind1 += 1 if ($w_ind1 == 0);
    $h_ind1 -= 1 if ($h_ind1 == $sqhgt);
    $w_ind1 -= 1 if ($w_ind1 == $sqwid);
    # Note that the direction of rotate is clock-wise. 
    # This may seem illogical mathematically, until you realise 
    # that the image coordinate system is relative to the top-left 
    # of the image instead of the mathematical norm of the bottom-left. 
    # The result is the angle of rotation is the oppisite of what you may 
    # logically expect. This is important to keep in mind when dealing 
    # with any form of image rotation, compared to a mathematical rotation.
    my $rot = normalised_hdg(int($az+0.5) - (90+$adj));
    my $arrow_head = "path 'M 0,0 l -15,-5 +5,+5 -5,+5 +15,-5 z'";
    my $msg = "-draw \"push graphic-context stroke blue fill skyblue translate $w_ind1,$h_ind1 rotate $rot ";
    $msg .= "$arrow_head pop graphic-context\" ";
    ${$rmsg} .= $msg;
}

sub add_img_circle($$$) {
    my ($rh,$clat,$clon) = @_;
    my $rmsg = ${$rh}{'rmsg'};
    my ($w_ind1,$h_ind1);
    my ($maxlat,$minlat,$maxlon,$minlon);
    my ($w_dpp,$h_dpp);
    my ($sqwid,$sqhgt);
    $maxlat = ${$rh}{'max_lat'};
    $minlat = ${$rh}{'min_lat'};
    $maxlon = ${$rh}{'max_lon'};
    $minlon = ${$rh}{'min_lon'};
    $w_dpp = ${$rh}{'w_dpp'};
    $h_dpp = ${$rh}{'h_dpp'};
    $sqwid = ${$rh}{'sq_wid_adj'};
    $sqhgt = ${$rh}{'sq_hgt_adj'};

    $w_ind1 = int((($clon - $minlon) * $w_dpp) + 0.5); # get degrees/pixels from left edge
    $h_ind1 = int((($clat - $minlat) * $h_dpp) + 0.5); # get degrees/pixels from bottom edge
    $h_ind1 = $sqhgt - $h_ind1;
    $h_ind1 += 1 if ($h_ind1 == 0);
    $w_ind1 += 1 if ($w_ind1 == 0);
    $h_ind1 -= 1 if ($h_ind1 == $sqhgt);
    $w_ind1 -= 1 if ($w_ind1 == $sqwid);
    ${$rmsg} .= "-draw \"circle ".($w_ind1-1).",$h_ind1 ".($w_ind1+1).",$h_ind1\" ";
}

sub add_arrow_to_center($$$$$) {
    my ($rh,$plat12,$plon12,$plat13,$plon13) = @_;
    my ($clat,$clon,$az1,$az2,$dist);
    # get heading and distance
    fg_geo_inverse_wgs_84($plat12,$plon12,$plat13,$plon13,\$az1,\$az2,\$dist);
    # get center point
    fg_geo_direct_wgs_84($plat12,$plon12,$az1,$dist/2,\$clat,\$clon,\$az2);
    add_img_circle($rh,$clat,$clon);
    add_img_arrow($rh,$clat,$clon,$az1);
}

sub add_img_line($$$$$$) {
    my ($rh,$lat1,$lon1,$lat2,$lon2,$type) = @_;
    my ($w_ind1,$h_ind1,$w_ind2,$h_ind2);
    my ($maxlat,$minlat,$maxlon,$minlon);
    my ($w_dpp,$h_dpp);
    my ($sqwid,$sqhgt,$draw);
    my $rmsg = ${$rh}{'rmsg'};
    $maxlat = ${$rh}{'max_lat'};
    $minlat = ${$rh}{'min_lat'};
    $maxlon = ${$rh}{'max_lon'};
    $minlon = ${$rh}{'min_lon'};
    $w_dpp = ${$rh}{'w_dpp'};
    $h_dpp = ${$rh}{'h_dpp'};
    $sqwid = ${$rh}{'sq_wid_adj'};
    $sqhgt = ${$rh}{'sq_hgt_adj'};

    $w_ind1 = int((($lon1 - $minlon) * $w_dpp) + 0.5); # get degrees/pixels from left edge
    $h_ind1 = int((($lat1 - $minlat) * $h_dpp) + 0.5); # get degrees/pixels from bottom edge
    $h_ind1 = $sqhgt - $h_ind1;
    $h_ind1 += 1 if ($h_ind1 == 0);
    $w_ind1 += 1 if ($w_ind1 == 0);
    $h_ind1 -= 1 if ($h_ind1 == $sqhgt);
    $w_ind1 -= 1 if ($w_ind1 == $sqwid);
    # ${$rmsg} .= "-draw \"circle ".($w_ind1-1).",$h_ind1 ".($w_ind1+1).",$h_ind1\" ";

    $w_ind2 = int((($lon2 - $minlon) * $w_dpp) + 0.5); # get degrees/pixels from left edge
    $h_ind2 = int((($lat2 - $minlat) * $h_dpp) + 0.5); # get degrees/pixels from bottom edge
    $h_ind2 = $sqhgt - $h_ind2;
    #$w_ind2 += 1 if ($w_ind2 == 0);
    $h_ind2 += 1 if ($h_ind2 == 0);
    $w_ind2 += 1 if ($w_ind2 == 0);
    $h_ind2 -= 1 if ($h_ind2 == $sqhgt);
    $w_ind2 -= 1 if ($w_ind2 == $sqwid);
    # ${$rmsg} .= "-draw \"circle ".($w_ind2-1).",$h_ind2 ".($w_ind2+1).",$h_ind2\" ";
    # -draw "stroke-dasharray 5 3       path 'M 10,20 L 90,20'"
    #if ($try_dash_line) {
    if ($type == 1) {
        $draw = "-draw \"stroke-dasharray 5 3 path 'M $w_ind1,$h_ind1 L $w_ind2,$h_ind2'\" ";
    } else {
        $draw = "-draw \"line $w_ind1,$h_ind1 $w_ind2,$h_ind2\" ";
    }
    ${$rmsg} .= $draw;
}

sub draw_text_at_pos($$$$$) {
    my ($rh,$x_txt,$y_txt,$text,$flag) = @_;
    my $rmsg = ${$rh}{'rmsg'};
    # -draw "text 25,65 'Anthony'"
    my $msg = "-draw \"text $x_txt,$y_txt '$text'\" ";
    ${$rh}{'x_txt'} = $x_txt;
    ${$rh}{'y_txt'} = $y_txt;
    ${$rmsg} .= $msg;
}

sub draw_text_at_latlon($$$$$) {
    my ($rh,$lat,$lon,$text,$flag) = @_;
    my ($w_ind1,$h_ind1);
    my ($maxlat,$minlat,$maxlon,$minlon);
    my ($w_dpp,$h_dpp);
    my ($sqwid,$sqhgt,$draw);
    my ($x_txt,$y_txt);
    my $len = length($text);
    return if ($len == 0);
    my $w_adj = $len * 3;
    my $h_adj = 5;

    $maxlat = ${$rh}{'max_lat'};
    $minlat = ${$rh}{'min_lat'};
    $maxlon = ${$rh}{'max_lon'};
    $minlon = ${$rh}{'min_lon'};
    $w_dpp = ${$rh}{'w_dpp'};
    $h_dpp = ${$rh}{'h_dpp'};
    $sqwid = ${$rh}{'sq_wid_adj'};
    $sqhgt = ${$rh}{'sq_hgt_adj'};

    $w_ind1 = int((($lon - $minlon) * $w_dpp) + 0.5); # get degrees/pixels from left edge
    $h_ind1 = int((($lat - $minlat) * $h_dpp) + 0.5); # get degrees/pixels from bottom edge
    $h_ind1 = $sqhgt - $h_ind1;
    
    $x_txt = $w_ind1 + 3;
    $y_txt = $h_ind1 + 4;
    # modify the position to within the page
    $x_txt += 3 if (($x_txt - 3) <= 0);
    $x_txt -= 3 if (($x_txt + 3) >= $sqwid);
    $y_txt -= 5 if (($y_txt + 5) >= $sqhgt);
    $y_txt += 5 if (($y_txt - 5) <= 0);

    draw_text_at_pos($rh,$x_txt,$y_txt,$text,$flag);
}


sub paint_user_points($$) {
    my ($rh,$show) = @_;
    my ($tllat,$tllon,$bllat,$bllon,$brlat,$brlon,$trlat,$trlon);
    $tllat = ${$rh}{'tl_lat'};
    $tllon = ${$rh}{'tl_lon'};
    $bllat = ${$rh}{'bl_lat'};
    $bllon = ${$rh}{'bl_lon'};
    $brlat = ${$rh}{'br_lat'};
    $brlon = ${$rh}{'br_lon'};
    $trlat = ${$rh}{'tr_lat'};
    $trlon = ${$rh}{'tr_lon'};
    my ($minlat,$maxlat,$minlon,$maxlon);
    my ($latdegs,$londegs,$sqwid,$sqhgt);
    $minlat = $bad_latlon;
    $maxlat = -$bad_latlon;
    $minlon = $bad_latlon;
    $maxlon = -$bad_latlon;

    my ($u_lat,$u_lon,$t_lat,$t_lon,$i,$cnt,$ru,$clat,$clon);

    # TEST SET OF 4, PLUS USRE'S POSITION
    $maxlat = $tllat if ($tllat > $maxlat);
    $maxlat = $bllat if ($bllat > $maxlat);
    $maxlat = $brlat if ($brlat > $maxlat);
    $maxlat = $trlat if ($trlat > $maxlat);

    $minlat = $tllat if ($tllat < $minlat);
    $minlat = $bllat if ($bllat < $minlat);
    $minlat = $brlat if ($brlat < $minlat);
    $minlat = $trlat if ($trlat < $minlat);

    $maxlon = $tllon if ($tllon > $maxlon);
    $maxlon = $bllon if ($bllon > $maxlon);
    $maxlon = $brlon if ($brlon > $maxlon);
    $maxlon = $trlon if ($trlon > $maxlon);

    $minlon = $tllon if ($tllon < $minlon);
    $minlon = $bllon if ($bllon < $minlon);
    $minlon = $brlon if ($brlon < $minlon);
    $minlon = $trlon if ($trlon < $minlon);

    # expand by User points
    if (defined ${$rh}{'user_points'}) {
        $ru = ${$rh}{'user_points'};
        $cnt = scalar @{$ru};
        for ($i = 0; $i < $cnt; $i++) {
            $u_lat = ${$ru}[$i][0];
            $u_lon = ${$ru}[$i][1];
            $t_lat = ${$ru}[$i][2];
            $t_lon = ${$ru}[$i][3];
            $maxlat = $u_lat if ($u_lat > $maxlat);
            $minlat = $u_lat if ($u_lat < $minlat);
            $maxlon = $u_lon if ($u_lon > $maxlon);
            $minlon = $u_lon if ($u_lon < $minlon);
        }
    }
    my $key = '';
    my $rcnt = 0;
    my $pcnt = 0;
    my ($rrwys,$rpatts);
    my ($elat1,$elon1,$elat2,$elon2);
    my ($plat11,$plon11,$plat12,$plon12,$plat13,$plon13,$plat21,$plon21);
    #my ($plat21,$plon21,$plat22,$plon22,$plat23,$plon23,$plat11,$plon11);
    if ((defined ${$rh}{'runways'})&&(defined ${$rh}{'pattern'})&&(defined ${$rh}{'airport'})) {
        $rrwys = ${$rh}{'runways'};
        $rpatts = ${$rh}{'pattern'};
        $key = ${$rh}{'airport'};
        $rcnt = scalar @{$rrwys};
        $pcnt = scalar @{$rpatts};
        prt("Adding $rcnt runways, and $pcnt patterns...\n");
        for ($i = 0; $i < $rcnt; $i++) {
            $elat1 = ${$rrwys}[$i][$RW_LLAT];
            $elon1 = ${$rrwys}[$i][$RW_LLON];
            $elat2 = ${$rrwys}[$i][$RW_RLAT];
            $elon2 = ${$rrwys}[$i][$RW_RLON];
            set_min_max(\$maxlat,\$minlat,\$maxlon,\$minlon,$elat1,$elon1);
            set_min_max(\$maxlat,\$minlat,\$maxlon,\$minlon,$elat2,$elon2);
        }
        for ($i = 0; $i < $pcnt; $i++) {
            #                   0       1       2       3       4       5       6       7
            # push(@{$rpatts}, [$plat11,$plon11,$plat12,$plon12,$plat13,$plon13,$plat21,$plon21,$clat,$clon,$rlen,$rhdg]);
            # push(@{$rpatts}, [$plat21,$plon21,$plat22,$plon22,$plat23,$plon23,$plat11,$plon11,$clat,$clon,$rlen,$hdgr]);
            $plat11 = ${$rpatts}[$i][0];
            $plon11 = ${$rpatts}[$i][1];
            $plat12 = ${$rpatts}[$i][2];
            $plon12 = ${$rpatts}[$i][3];
            $plat13 = ${$rpatts}[$i][4];
            $plon13 = ${$rpatts}[$i][5];
            $plat21 = ${$rpatts}[$i][6];
            $plon21 = ${$rpatts}[$i][7];
            set_min_max(\$maxlat,\$minlat,\$maxlon,\$minlon,$plat11,$plon11);
            set_min_max(\$maxlat,\$minlat,\$maxlon,\$minlon,$plat12,$plon12);
            set_min_max(\$maxlat,\$minlat,\$maxlon,\$minlon,$plat13,$plon13);
            set_min_max(\$maxlat,\$minlat,\$maxlon,\$minlon,$plat21,$plon21);
        }
    }

    ${$rh}{'max_lat'} = $maxlat;
    ${$rh}{'min_lat'} = $minlat;
    ${$rh}{'max_lon'} = $maxlon;
    ${$rh}{'min_lon'} = $minlon;
    my $lon_factor = 2;
    $latdegs = $maxlat - $minlat;
    $londegs = ($maxlon - $minlon) * $lon_factor;
    ${$rh}{'lon_degs'} = $londegs;
    ${$rh}{'lat_degs'} = $latdegs;
    $sqhgt = int(($latdegs * 10000) + 0.5);
    $sqwid = int(($londegs * 10000) + 0.5);
    ${$rh}{'sq_wid'} = $sqwid;
    ${$rh}{'sq_hgt'} = $sqhgt;

    # scale the image DOWN/up
    my $targ_wid = 600;
    my $ratio = $sqwid / $sqhgt;
    #if (($sqwid > $targ_wid) || ($sqhgt > $targ_wid)) {
        if ($ratio > 1) {   # width > height
            $sqwid = $targ_wid; # set target width
            $sqhgt = int($targ_wid / $ratio); # and calculate NEW height
	    } else {
            # height > width = set target height, adjust width accordingly
			$sqwid = int($targ_wid * $ratio); # calculate width
			$sqhgt = $targ_wid; # and set target width
        }
    #}
    ${$rh}{'sq_wid_adj'} = $sqwid;
    ${$rh}{'sq_hgt_adj'} = $sqhgt;

    my ($w_dpp,$h_dpp,$w_ind1,$h_ind1,$w_ind2,$h_ind2,$w_ind3,$h_ind3,$w_ind4,$h_ind4,$msg);
    my ($h_indu,$w_indu,$h_indt,$w_indt);
    my ($x_txt,$y_txt,$txt);
    $w_dpp = ($sqwid / $londegs) * $lon_factor;
    $h_dpp = $sqhgt / $latdegs;
    ${$rh}{'w_dpp'} = $w_dpp;
    ${$rh}{'h_dpp'} = $h_dpp;
    ${$rh}{'rmsg'} = \$msg;
    # Imagmagick colors - http://www.imagemagick.org/script/color.php
    #$msg = "convert -size ${sqwid}x${sqhgt} xc:wheat +antialias -fill white ";
    $msg = "convert -size ${sqwid}x${sqhgt} xc:wheat -fill white ";
    #      -stroke Gold        -draw "path 'M 20,70  A 1,1 0 0,1 80,50'" \
    #      -stroke DodgerBlue  -draw "line 30,10  50,80" \
    # a dasshed line -draw "stroke-dasharray 5 3 path 'M 10,20 L 90,20'"
    #      -stroke Red         -draw "circle 80,60  82,60" test.gif
    # -draw 'circle 50,30 50,55' - centered on point 50,30, ???
    #$msg .= "-font Courier -pointsize 12 -strokewidth 0.5 ";
    # As such with a default "-density" of 72dpi (at which 1 point = 1 pixel) 
    # a 12 point font should have 12 pixels separation between the 
    # baselines to two lines of text. 
    $msg .= "-pointsize 12 -strokewidth 0.5 ";
    if (length($key) && $rcnt) {
        if ($pcnt) {
            ${$rh}{'-stroke'} = "SlateGray";
            $msg .= "-stroke SlateGray ";
            for ($i = 0; $i < $pcnt; $i++) {
                #                   0       1       2       3       4       5       6       7       8     9
                # push(@{$rpatts}, [$plat11,$plon11,$plat12,$plon12,$plat13,$plon13,$plat21,$plon21,$clat,$clon]);
                # push(@{$rpatts}, [$plat21,$plon21,$plat22,$plon22,$plat23,$plon23,$plat11,$plon11,$clat,$clon]);
                #      YGIL LAYOUT
                # 22 ----- 21 ----- 13
                #  |                 |
                #  |       15        |
                #  |        |        |
                #  |        33       |
                #  |                 |
                #  23 ----- 11 ----- 12
                $plat11 = ${$rpatts}[$i][0];
                $plon11 = ${$rpatts}[$i][1];
                $plat12 = ${$rpatts}[$i][2];
                $plon12 = ${$rpatts}[$i][3];
                $plat13 = ${$rpatts}[$i][4];
                $plon13 = ${$rpatts}[$i][5];
                $plat21 = ${$rpatts}[$i][6];
                $plon21 = ${$rpatts}[$i][7];
                $clat   = ${$rpatts}[$i][8];
                $clon   = ${$rpatts}[$i][9];
                add_img_line($rh,$plat11,$plon11,$plat12,$plon12,1);
                add_img_line($rh,$plat12,$plon12,$plat13,$plon13,1);
                add_img_line($rh,$plat13,$plon13,$plat21,$plon21,1);
                add_img_line($rh,$plat21,$plon21,$plat11,$plon11,1);
                if ($switch_circuit && ($i == 0)) {
                    add_arrow_to_center($rh,$plat11,$plon11,$plat12,$plon12);
                    add_arrow_to_center($rh,$plat12,$plon12,$plat13,$plon13);
                    add_arrow_to_center($rh,$plat13,$plon13,$plat21,$plon21);
                    add_arrow_to_center($rh,$clat,$clon,$plat11,$plon11);
                    add_arrow_to_center($rh,$plat21,$plon21,$clat,$clon);
                } elsif (!$switch_circuit && ($i == 1)) {
                    add_arrow_to_center($rh,$plat11,$plon11,$plat12,$plon12);
                    add_arrow_to_center($rh,$plat12,$plon12,$plat13,$plon13);
                    add_arrow_to_center($rh,$plat13,$plon13,$plat21,$plon21);
                    add_arrow_to_center($rh,$plat21,$plon21,$clat,$clon);
                    add_arrow_to_center($rh,$clat,$clon,$plat11,$plon11);
                }
            }
        }

        # DRAW THE RUNWAY
        # set color
        $msg .= "-stroke blue ";
        for ($i = 0; $i < $rcnt; $i++) {
            $elat1 = ${$rrwys}[$i][$RW_LLAT];
            $elon1 = ${$rrwys}[$i][$RW_LLON];
            $elat2 = ${$rrwys}[$i][$RW_RLAT];
            $elon2 = ${$rrwys}[$i][$RW_RLON];
            $w_ind1 = int((($elon1 - $minlon) * $w_dpp) + 0.5); # get degrees/pixels from left edge
            $h_ind1 = int((($elat1 - $minlat) * $h_dpp) + 0.5); # get degrees/pixels from bottom edge
            $h_ind1 = $sqhgt - $h_ind1;
            $h_ind1 += 1 if ($h_ind1 == 0);
            $w_ind1 += 1 if ($w_ind1 == 0);
            $h_ind1 -= 1 if ($h_ind1 == $sqhgt);
            $w_ind1 -= 1 if ($w_ind1 == $sqwid);
            $msg .= "-draw \"circle ".($w_ind1-1).",$h_ind1 ".($w_ind1+1).",$h_ind1\" ";

            $w_ind2 = int((($elon2 - $minlon) * $w_dpp) + 0.5); # get degrees/pixels from left edge
            $h_ind2 = int((($elat2 - $minlat) * $h_dpp) + 0.5); # get degrees/pixels from bottom edge
            $h_ind2 = $sqhgt - $h_ind2;
            #$w_ind2 += 1 if ($w_ind2 == 0);
            $h_ind2 += 1 if ($h_ind2 == 0);
            $w_ind2 += 1 if ($w_ind2 == 0);
            $h_ind2 -= 1 if ($h_ind2 == $sqhgt);
            $w_ind2 -= 1 if ($w_ind2 == $sqwid);
            $msg .= "-draw \"circle ".($w_ind2-1).",$h_ind2 ".($w_ind2+1).",$h_ind2\" ";

            $msg .= "-draw \"line $w_ind1,$h_ind1 $w_ind2,$h_ind2\" ";
            my $clat = ${$rrwys}[$i][$RW_CLAT];
            my $clon = ${$rrwys}[$i][$RW_CLON];
            draw_text_at_latlon($rh,$clat,$clon,"YGIL",0);
        }
    }
    # set color
    $msg .= "-stroke black ";

    $x_txt = $sqwid - 100;
    #$y_txt = 12;
    $y_txt = 20;
    draw_text_at_pos($rh,$x_txt,$y_txt,"Wind: SSE 8Kt",0);

    # tllat,$tllon on square
    $w_ind1 = int((($tllon - $minlon) * $w_dpp) + 0.5); # get degrees/pixels from left edge
    $h_ind1 = int((($tllat - $minlat) * $h_dpp) + 0.5); # get degrees/pixels from bottom edge
    $h_ind1 = $sqhgt - $h_ind1;
    $h_ind1 += 1 if ($h_ind1 == 0);
    $w_ind1 += 1 if ($w_ind1 == 0);
    $h_ind1 -= 1 if ($h_ind1 == $sqhgt);
    $w_ind1 -= 1 if ($w_ind1 == $sqwid);
    $msg .= "-draw \"circle ".($w_ind1-1).",$h_ind1 ".($w_ind1+1).",$h_ind1\" ";

    # bllat,$bllon on square
    $w_ind2 = int((($bllon - $minlon) * $w_dpp) + 0.5); # get degrees/pixels from left edge
    $h_ind2 = int((($bllat - $minlat) * $h_dpp) + 0.5); # get degrees/pixels from bottom edge
    $h_ind2 = $sqhgt - $h_ind2;
    #$w_ind2 += 1 if ($w_ind2 == 0);
    $h_ind2 += 1 if ($h_ind2 == 0);
    $w_ind2 += 1 if ($w_ind2 == 0);
    $h_ind2 -= 1 if ($h_ind2 == $sqhgt);
    $w_ind2 -= 1 if ($w_ind2 == $sqwid);
    $msg .= "-draw \"circle ".($w_ind2-1).",$h_ind2 ".($w_ind2+1).",$h_ind2\" ";

    # brlat,$brlon on square
    $w_ind3 = int((($brlon - $minlon) * $w_dpp) + 0.5); # get degrees/pixels from left edge
    $h_ind3 = int((($brlat - $minlat) * $h_dpp) + 0.5); # get degrees/pixels from bottom edge
    $h_ind3 = $sqhgt - $h_ind3;
    #$w_ind3 += 1 if ($w_ind3 == 0);
    $h_ind3 += 1 if ($h_ind3 == 0);
    $w_ind3 += 1 if ($w_ind3 == 0);
    $h_ind3 -= 1 if ($h_ind3 == $sqhgt);
    $w_ind3 -= 1 if ($w_ind3 == $sqwid);
    $msg .= "-draw \"circle ".($w_ind3-1).",$h_ind3 ".($w_ind3+1).",$h_ind3\" ";

    # trlat,$trlon on square
    $w_ind4 = int((($trlon - $minlon) * $w_dpp) + 0.5); # get degrees/pixels from left edge
    $h_ind4 = int((($trlat - $minlat) * $h_dpp) + 0.5); # get degrees/pixels from bottom edge
    $h_ind4 = $sqhgt - $h_ind4;
    #$w_ind4 += 1 if ($w_ind4 == 0);
    $h_ind4 += 1 if ($h_ind4 == 0);
    $w_ind4 += 1 if ($w_ind4 == 0);
    $h_ind4 -= 1 if ($h_ind4 == $sqhgt);
    $w_ind4 -= 1 if ($w_ind4 == $sqwid);
    $msg .= "-draw \"circle ".($w_ind4-1).",$h_ind4 ".($w_ind4+1).",$h_ind4\" ";

    if (!$use_pattern) {
        $msg .= "-draw \"line $w_ind1,$h_ind1 $w_ind2,$h_ind2\" ";
        $msg .= "-draw \"line $w_ind2,$h_ind2 $w_ind3,$h_ind3\" ";
        $msg .= "-draw \"line $w_ind3,$h_ind3 $w_ind4,$h_ind4\" ";
        $msg .= "-draw \"line $w_ind4,$h_ind4 $w_ind1,$h_ind1\" ";
    }
    if (defined ${$rh}{'user_points'}) {
        $ru = ${$rh}{'user_points'};
        $cnt = scalar @{$ru};
        # u_lat,$u_lon on square
        for ($i = 0; $i < $cnt; $i++) {
            $u_lat = ${$ru}[$i][0];
            $u_lon = ${$ru}[$i][1];
            $t_lat = ${$ru}[$i][2];
            $t_lon = ${$ru}[$i][3];
            $w_indu = int((($u_lon - $minlon) * $w_dpp) + 0.5); # get degrees/pixels from left edge
            $h_indu = int((($u_lat - $minlat) * $h_dpp) + 0.5); # get degrees/pixels from bottom edge
            $h_indu = $sqhgt - $h_indu;
            $h_indu += 1 if ($h_indu == 0);
            $w_indu += 1 if ($w_indu == 0);
            $h_indu -= 1 if ($h_indu == $sqhgt);
            $w_indu -= 1 if ($w_indu == $sqwid);
            $msg .= "-stroke blue ";
            $msg .= "-draw \"circle ".($w_indu-1).",$h_indu ".($w_indu+1).",$h_indu\" ";
            if ($add_text_count) {
                # -draw "text 25,65 'Anthony'"
                # this point is the lower left of the text
                $txt = $i + 1;
                $x_txt = $w_indu + 3;
                $y_txt = $h_indu + 4;
                # modify the position to within the page
                $x_txt += 3 if (($x_txt - 3) <= 0);
                $x_txt -= 3 if (($x_txt + 3) >= $sqwid);
                $y_txt -= 5 if (($y_txt + 5) >= $sqhgt);
                $y_txt += 5 if (($y_txt - 5) <= 0);
                $msg .= "-draw \"text $x_txt,$y_txt '$txt'\" ";
            }
            # t_lat,$t_lon on square
            $w_indt = int((($t_lon - $minlon) * $w_dpp) + 0.5); # get degrees/pixels from left edge
            $h_indt = int((($t_lat - $minlat) * $h_dpp) + 0.5); # get degrees/pixels from bottom edge
            $h_indt = $sqhgt - $h_indt;
            $h_indt += 1 if ($h_indt == 0);
            $w_indt += 1 if ($w_indt == 0);
            $h_indt -= 1 if ($h_indt == $sqhgt);
            $w_indt -= 1 if ($w_indt == $sqwid);
            $msg .= "-draw \"circle ".($w_indu-1).",$h_indu ".($w_indu+1).",$h_indu\" ";
            # $msg .= "-stroke red -draw \"line $w_indu,$h_indu $w_indt,$h_indt\" ";
            $msg .= "-stroke red ";
            #add_img_line($rh,$w_indu,$h_indu,$w_indt,$h_indt);
            add_img_line($rh,$u_lat,$u_lon,$t_lat,$t_lon,0);
        }
    }

    $msg .= "$graf_file\n";
    $msg .= "imdisplay $graf_file\n";
    $msg .= "start $graf_file\n";
    write2file($msg,$graf_bat);
    prt("Written $graf_bat\n");
    if ($show) {
        set_lat_stg(\$maxlat);
        set_lat_stg(\$minlat);
        set_lon_stg(\$maxlon);
        set_lon_stg(\$minlon);
        prt("Square wid=$londegs hgt=$latdegs ${sqwid}X${sqhgt}\n");
        prt("TL $maxlat,$minlon\n");
        prt("BL $minlat,$minlon\n");
        prt("BR $minlat,$maxlon\n");
        prt("TR $maxlat,$maxlon\n");
    }
}

sub get_circuit_hash() {
    my %h = ();
    $h{'tl_lat'} = $tl_lat;
    $h{'tl_lon'} = $tl_lon;
    $h{'bl_lat'} = $bl_lat;
    $h{'bl_lon'} = $bl_lon;
    $h{'br_lat'} = $br_lat;
    $h{'br_lon'} = $br_lon;
    $h{'tr_lat'} = $tr_lat;
    $h{'tr_lon'} = $tr_lon;
    set_circuit_values(\%h,1);
    #set_most_values(\%h,1);
    return \%h;
}

sub get_mid_tl2bl($$) {
    my ($rlat,$rlon) = @_;
    my ($az1,$az2,$dist);
    fg_geo_inverse_wgs_84 ($tl_lat,$tl_lon,$bl_lat,$bl_lon,\$az1,\$az2,\$dist);
    my $dist2 = $dist / 2;  # get the center of this line
    my ($clat,$clon);
    #fg_geo_direct_wgs_84( $tl_lat, $tl_lon, $az1, $dist2, $rlat, $rlon, \$az2 );
    fg_geo_direct_wgs_84( $tl_lat, $tl_lon, $az1, $dist2, \$clat, \$clon, \$az2 );
    $az2 = normalised_hdg($az1 + 90); # turn 90 degrees, and get a point
    fg_geo_direct_wgs_84( $clat, $clon, $az2, $dist2, $rlat, $rlon, \$az1 );

}
sub get_mid_bl2br($$) {
    my ($rlat,$rlon) = @_;
    my ($az1,$az2,$dist);
    fg_geo_inverse_wgs_84 ($bl_lat,$bl_lon,$br_lat,$br_lon,\$az1,\$az2,\$dist);
    my $dist2 = $dist / 2;
    fg_geo_direct_wgs_84( $bl_lat, $bl_lon, $az1, $dist2, $rlat, $rlon, \$az2 );
}
sub get_mid_bl2br2($$) {
    my ($rlat,$rlon) = @_;
    my ($az1,$az2,$dist);
    fg_geo_inverse_wgs_84 ($bl_lat,$bl_lon,$br_lat,$br_lon,\$az1,\$az2,\$dist);
    my $dist2 = $dist / 2;
    my ($clat,$clon);
    fg_geo_direct_wgs_84( $bl_lat, $bl_lon, $az1, $dist2, \$clat, \$clon, \$az2 );
    $az2 = normalised_hdg($az1 + 90); # turn 90 degrees, and get a point
    fg_geo_direct_wgs_84( $clat, $clon, $az2, $dist2, $rlat, $rlon, \$az1 );
}

sub set_distances_bearings($$$$) {
    my ($rh,$lat,$lon,$msg) = @_;
    ${$rh}{'usr_lat'} = $lat;
    ${$rh}{'usr_lon'} = $lon;
    ${$rh}{'usr_msg'} = $msg;
    my ($tlat,$tlon);
    my ($az1,$az2,$dist);
    $msg = '';

    $tlat = ${$rh}{'tl_lat'}; # = -31.684063;
    $tlon = ${$rh}{'tl_lon'}; # = 148.614120;
    fg_geo_inverse_wgs_84 ($lat,$lon,$tlat,$tlon,\$az1,\$az2,\$dist);
    ${$rh}{'tl_az1'} = $az1;
    ${$rh}{'tl_az2'} = $az2;
    ${$rh}{'tl_dist'} = $dist;
    set_int_dist_stg5(\$dist);
    $msg .= "TL $dist ";
    $tlat = ${$rh}{'bl_lat'}; # = -31.723495;
    $tlon = ${$rh}{'bl_lon'}; # = 148.633003;
    fg_geo_inverse_wgs_84 ($lat,$lon,$tlat,$tlon,\$az1,\$az2,\$dist);
    ${$rh}{'bl_az1'} = $az1;
    ${$rh}{'bl_az2'} = $az2;
    ${$rh}{'bl_dist'} = $dist;
    set_int_dist_stg5(\$dist);
    $msg .= "BL $dist ";

    $tlat = ${$rh}{'br_lat'}; # = -31.716778;
    $tlon = ${$rh}{'br_lon'}; # = 148.666992;
    fg_geo_inverse_wgs_84 ($lat,$lon,$tlat,$tlon,\$az1,\$az2,\$dist);
    ${$rh}{'br_az1'} = $az1;    # from 'test' to BR point
    ${$rh}{'br_az2'} = $az2;
    ${$rh}{'br_dist'} = $dist;
    set_int_dist_stg5(\$dist);
    $msg .= "BR $dist ";

    $tlat = ${$rh}{'tr_lat'}; # = -31.672960;
    $tlon = ${$rh}{'tr_lon'}; # = 148.649139;
    fg_geo_inverse_wgs_84 ($lat,$lon,$tlat,$tlon,\$az1,\$az2,\$dist);
    ${$rh}{'tr_az1'} = $az1;
    ${$rh}{'tr_az2'} = $az2;
    ${$rh}{'tr_dist'} = $dist;
    set_int_dist_stg5(\$dist);
    $msg .= "TR $dist ";
    prt("Distances: $msg\n");

}

sub dist_less_or_equal5($$) {
    my ($dist,$min_dist) = @_;
    return 1 if ($dist <= $min_dist);
    my $d5 = $dist * 0.95;
    return 2 if ($d5 < $min_dist);
    return 0;
}

sub point_in_circuit($$$$) {
    my ($rh,$lat,$lon,$rres) = @_;
    my ($tl_lat,$tl_lon,$bl_lat,$bl_lon,$br_lat,$br_lon,$tr_lat,$tr_lon);
    my $ret = 0;
    $tl_lat = ${$rh}{'tl_lat'};
    $tl_lon = ${$rh}{'tl_lon'};
    $bl_lat = ${$rh}{'bl_lat'};
    $bl_lon = ${$rh}{'bl_lon'};
    $br_lat = ${$rh}{'br_lat'};
    $br_lon = ${$rh}{'br_lon'};
    $tr_lat = ${$rh}{'tr_lat'};
    $tr_lon = ${$rh}{'tr_lon'};
    my $res = '';
    if (Point_Inside_Triange($lat,$lon,$tl_lat,$tl_lon, $bl_lat,$bl_lon, $br_lat,$br_lon)) {
        $res = "in circuit (1st tri)";
        $ret = 1;
    } elsif (Point_Inside_Triange($lat,$lon,$tl_lat,$tl_lon, $tr_lat,$tr_lon, $br_lat,$br_lon)) {
        $res = "in curcuit (2nd tri)";
        $ret = 1;
    } else {
        $res = "NOT in circuit.";
    }
    ${$rres} = $res;
    return $ret;
}

sub process_lat_lon($$$$) {
    my ($rh,$lat,$lon,$msg) = @_;
    my $res = '';
    my $ptinc = point_in_circuit($rh,$lat,$lon,\$res);
    prt("\nSolving lat=$lat, lon=$lon $msg, $res \n");
    # get distances, and headings...
    my ($tlat,$tlon,$az1,$az2,$dist);
    set_distances_bearings($rh,$lat,$lon,$msg);
    # 1: get closest point, on circuit
    # 1 Mile = 1609.344 Meters.
    my $min_dist = 12000 * 1700;
    my $hdto = 'unsolved';
    my $targ_lon = $bad_latlon;
    my $targ_lat = $bad_latlon;
    my $targ_dist = 1000000000;
    my ($hdg);

    $tlat = ${$rh}{'tl_lat'}; # = -31.684063;
    $tlon = ${$rh}{'tl_lon'}; # = 148.614120;
    $az1  = ${$rh}{'tl_az1'};
    $az2  = ${$rh}{'tl_az2'};
    $dist = ${$rh}{'tl_dist'};
    $res = 0;
    #if ($dist < $min_dist) {
        $min_dist = $dist;
        $hdto = "to top left ($res)";
        $hdg = $az1;
        $targ_lat = $tlat;
        $targ_lon = $tlon;
        $targ_dist = $dist;
    #}

    $tlat = ${$rh}{'bl_lat'}; # = -31.684063;
    $tlon = ${$rh}{'bl_lon'}; # = 148.614120;
    $az1  = ${$rh}{'bl_az1'};
    $az2  = ${$rh}{'bl_az2'};
    $dist = ${$rh}{'bl_dist'};
    #if ($dist <= $min_dist) {
    $res = dist_less_or_equal5($dist,$min_dist);
    if ($res) {
        $min_dist = $dist;
        $hdto = "to bottom left ($res)";
        $hdg = $az1;
        $targ_lat = $tlat;
        $targ_lon = $tlon;
        $targ_dist = $dist;
    }
    $tlat = ${$rh}{'br_lat'}; # = -31.684063;
    $tlon = ${$rh}{'br_lon'}; # = 148.614120;
    $az1  = ${$rh}{'br_az1'};
    $az2  = ${$rh}{'br_az2'};
    $dist = ${$rh}{'br_dist'};
    #if ($dist < $min_dist) {
    $res = dist_less_or_equal5($dist,$min_dist);
    if ($res) {
        $min_dist = $dist;
        $hdto = "to bottom right ($res)";
        $hdg = $az1;
        $targ_lat = $tlat;
        $targ_lon = $tlon;
        $targ_dist = $dist;
    }
    $tlat = ${$rh}{'tr_lat'}; # = -31.684063;
    $tlon = ${$rh}{'tr_lon'}; # = 148.614120;
    $az1  = ${$rh}{'tr_az1'};
    $az2  = ${$rh}{'tr_az2'};
    $dist = ${$rh}{'tr_dist'};
    #if ($dist <= $min_dist) {
    $res = dist_less_or_equal5($dist,$min_dist);
    if ($res) {
        $min_dist = $dist;
        $hdto = "to top right ($res)";
        $hdg = $az1;
        $targ_lat = $tlat;
        $targ_lon = $tlon;
        $targ_dist = $dist;
    }

    $tlat = ${$rh}{'tl_lat'}; # = -31.684063;
    $tlon = ${$rh}{'tl_lon'}; # = 148.614120;
    $az1  = ${$rh}{'tl_az1'};
    $az2  = ${$rh}{'tl_az2'};
    $dist = ${$rh}{'tl_dist'};
    $res = dist_less_or_equal5($dist,$min_dist);
    if ($res) {
        $min_dist = $dist;
        $hdto = "to top left ($res)";
        $hdg = $az1;
        $targ_lat = $tlat;
        $targ_lon = $tlon;
        $targ_dist = $dist;
    }

    if ($ptinc) {
        # if in the circuit, always head for the end of downwind, before turn cross
        $tlat = ${$rh}{'bl_lat'}; # = -31.684063;
        $tlon = ${$rh}{'bl_lon'}; # = 148.614120;
        $az1  = ${$rh}{'bl_az1'};
        $az2  = ${$rh}{'bl_az2'};
        $dist = ${$rh}{'bl_dist'};
        $hdto = "in circuit ($ptinc)";
        $hdg = $az1;
        $targ_lat = $tlat;
        $targ_lon = $tlon;
        $targ_dist = $dist;
    }

    set_most_values_plus($rh,1,$lat,$lon,$targ_lat,$targ_lon);

    # for display
    set_hdg_stg3(\$hdg);
    set_int_dist_stg5(\$targ_dist);
    set_lat_stg(\$targ_lat);
    set_lon_stg(\$targ_lon);
    prt("\nHeading: $hdto, hdg=$hdg, to $targ_lat,$targ_lon, at $targ_dist m.\n");

}

sub is_in_circuit($$$$$) {
    my ($rh,$lat,$lon,$msg,$i) = @_;
    # $tl_lat,$tl_lon, $bl_lat, $bl_lon, $br_lat, $br_lon, $tr_lat, $tr_lon
    # Into two triangles
    # 1 - $tl_lat,$tl_lon, $bl_lat,$bl_lon, $br_lat,$br_lon
    # 2 - $tl_lat,$tl_lon, $tr_lat,$tr_lon, $br_lat,$br_lon
    #sub Point_Inside_Triange($$$$$$$$) {
    # my ($px,$py,$x1,$y1,$x2,$y2,$x3,$y3) = @_;
    if (Point_Inside_Triange($lat,$lon,$tl_lat,$tl_lon, $bl_lat,$bl_lon, $br_lat,$br_lon)) {
        prt("Point $i inside first triangle\n");
    } elsif (Point_Inside_Triange($lat,$lon,$tl_lat,$tl_lon, $tr_lat,$tr_lon, $br_lat,$br_lon)) {
        prt("Point $i inside second triangle\n");
    } else {
        prt("Point $i NOT in circuit\n");
    }
}

# Airport Line. eg '1 5355 1 0 KABQ Albuquerque Intl Sunport'
# 0 1    - this as an airport header line. 16 is a seaplane/floatplane base, 17 a heliport.
# 1 5355 - Airport elevation (in feet above MSL).  
# 2 1    - Airport has a control tower (1=yes, 0=no).
# 3 0   - Display X-Plane’s default airport buildings (1=yes, 0=no).
# 4 KABQ   - Identifying code for the airport (the ICAO code, if one exists).
# 5+Albuquerque Intl Sunport - Airport name.
sub find_apt($) {
    my $icao = shift;
    my $t1 = [gettimeofday];
    my $aptdat = "X:\\fgdata\\Airports\\apt.dat.gz";
    if (! -f $aptdat) {
        prt("Failed to 'stat' file '$aptdat'!\n");
        pgm_exit(1,"");
    }
    prt("Processing file $aptdat... for ICAO $icao... moment...\n");
    if (!open(INF,"gzip -cdq $aptdat|")) {
        prt("Failed to 'open' file '$aptdat'!\n");
        pgm_exit(1,"");
    }
    my ($line,@arr,$type,$len);
    my $ver = 0;
    my $lncnt = 0;
    while ($line = <INF>) {
        $lncnt++;
        chomp $line;
        $line = trim_all($line);
        if ($line =~ /^(\d+)\s+Version/) {
            $ver = $1;
            last;
        }
    }
    my $fnd = 0;
    if ($ver) {
        while ($line = <INF>) {
            $lncnt++;
            chomp $line;
            $line = trim_all($line);
            $len = length($line);
            next if ($len == 0);
            @arr = split(/\s+/,$line);
            $type = $arr[0];
            last if ($type == 99);
            if ((($type == 1)||($type == 16)||($type == 17)) &&
                ($arr[4] eq $icao) ) {
            # if (($type == 1) && ($arr[4] eq $icao))
                my @a = @arr;
                push(@my_airport,\@a);
                while ($line = <INF>) {
                    $lncnt++;
                    chomp $line;
                    $line = trim_all($line);
                    $len = length($line);
                    next if ($len == 0);
                    @arr = split(/\s+/,$line);
                    $type = $arr[0];
                    last if ($type == 99);
                    last if (($type == 1)||($type == 16)||($type == 17));
                    my @a2 = @arr;
                    push(@my_airport,\@a2);
                }
                $fnd = $lncnt;
                last;
            }
        }
    } else {
        pgm_exit(1,"Version NOT found in $aptdat!\n");
    }
    while ($line = <INF>) {
        $lncnt++;
    }
    my $elap = secs_HHMMSS( tv_interval( $t1, [gettimeofday] ) );
    $line = scalar @my_airport;
    if ($fnd) {
        prt("Found apt $icao, $line components, at line $fnd of $lncnt lines... in $elap\n");
    } else {
        pgm_exit(1,"ICAO $icao NOT FOUND in $aptdat!\n");
    }

    close INF;
}

# @runways reference
# 0   1=lat      2=lon      3=s 4=hdg  5=len 6=offsets 7=stopway 8=wid 9=lights 10=surf 11 12 13   14 15
# 10  36.962213  127.031071 14x 131.52  8208 1595.0620 0000.0000 150   321321   1       0  3  0.25 0  0300.0300
# 11=shoulder 12=marks 13=smooth 14=signs 15=GS angles
# 0           3        0.25      0        0300.0300
# ====================================================

# Land Runway
# 0  - 100 - Row code for a land runway (the most common) 100
# 1  - 29.87 - Width of runway in metres Two decimal places recommended. Must be >= 1.00
# 2  - 3 - Code defining the surface type (concrete, asphalt, etc) Integer value for a Surface Type Code
# 3  - 0 - Code defining a runway shoulder surface type 0=no shoulder, 1=asphalt shoulder, 2=concrete shoulder
# 4  - 0.15 - Runway smoothness (not used by X-Plane yet) 0.00 (smooth) to 1.00 (very rough). Default is 0.25
# 5  - 0 - Runway centre-line lights 0=no centerline lights, 1=centre line lights
# 6  - 0 - Runway edge lighting (also implies threshold lights) 0=no edge lights, 2=medium intensity edge lights
# 7  - 1 - Auto-generate distance-remaining signs (turn off if created manually) 0=no auto signs, 1=auto-generate signs
# The following fields are repeated for each end of the runway
# 8  - 13L - Runway number (eg. 31R, 02). Leading zeros are required. Two to three characters. Valid suffixes: L, R or C (or blank)
# 9  - 47.53801700 - Latitude of runway threshold (on runway centerline) in decimal degrees Eight decimal places supported
# 10 - -122.30746100 - Longitude of runway threshold (on runway centerline) in decimal degrees Eight decimal places supported
# 11 - 73.15 - Length of displaced threshold in metres (this is included in implied runway length) Two decimal places (metres). Default is 0.00
# 12 - 0.00 - Length of overrun/blast-pad in metres (not included in implied runway length) Two decimal places (metres). Default is 0.00
# 13 - 2 - Code for runway markings (Visual, non-precision, precision) Integer value for Runway Marking Code
# 14 - 0 - Code for approach lighting for this runway end Integer value for Approach Lighting Code
# 15 - 0 - Flag for runway touchdown zone (TDZ) lighting 0=no TDZ lighting, 1=TDZ lighting
# 16 - 1 - Code for Runway End Identifier Lights (REIL) 0=no REIL, 1=omni-directional REIL, 2=unidirectional REIL
# 17 - 31R
# 18 - 47.52919200
# 19 - -122.30000000
# 20 - 110.95 
# 21 - 0.00 
# 22 - 2
# 23 - 0
# 24 - 0
# 25 - 1

sub find_runway() {
    my ($ra,$type,$clat,$clon,$cnt,$len,$width,$hdg,$tmp,$hdgr);
    my ($rwy1,$lat1,$lon1);
    my ($rwy2,$lat2,$lon2);
    my ($az1,$az2,$dist,$rwlen2);
    my ($isrwy,$dhdg);
    my ($alt,$icao,$name);
    my ($tlat,$tlon,$alat,$alon);
    $cnt = scalar @my_airport;
    prt("Choose runway, from $cnt apt components...\n");
    my @runways = ();
    $tlat = 0;
    $tlon = 0;
    my $gothdr = 0;
    foreach $ra (@my_airport) {
        $cnt = scalar @{$ra};
        $type = ${$ra}[0];
        $isrwy = 0;
        if (($type == 1)||($type == 16)||($type == 17)) {
            # @arr2 = split(/\s+/,$apt);
            $alt = ${$ra}[1]; # Airport (general) ALTITUDE AMSL
            # $actl = $arr2[2]; # control tower
            # $abld = $arr2[3]; # buildings
            $icao = ${$ra}[4]; # ICAO
            $name = join(' ', splice(@{$ra},5)); # Name
            $gothdr = 1;
        } elsif ($type == 10) {
            $rwy1 = ${$ra}[3];  # xxx just a taxiway
            next if ($rwy1 eq 'xxx');
            $clat  = ${$ra}[1];
            $clon  = ${$ra}[2];
            $rwy1 = ${$ra}[3];  # xxx just a taxiway
            $hdg  = ${$ra}[4];  # 4=hdg
            $len  = ${$ra}[5];  # 5=len
            $dist = $len * $SG_FEET_TO_METER;
            $rwlen2 = $dist / 2;
            my $hdgr = $hdg + 180;
            $hdgr -= 360 if ($hdgr >= 360);
            fg_geo_direct_wgs_84( $clat, $clon, $hdg, $rwlen2, \$lat1, \$lon1, \$az1 );
            fg_geo_direct_wgs_84( $clat, $clon, $hdgr, $rwlen2, \$lat2, \$lon2, \$az2 );
            $rwy1 =~ s/x+$//;   # REMOVE any TRAILING 'x', but may have 'L', 'R', 'C', 'S' appended
            if ($rwy1 =~ /^\d+$/) {
                ### $rwy2 = $rwy1 * 10; # 2010-12-15 - get opposite end numbers
                $rwy2 = ($rwy1 > 18) ? $rwy1 - 18 : $rwy1 + 18; 
                $rwy2 = '36' if ($rwy2 == 0);
            } else {
                $rwy2 = $rwy1;
                $rwy2 =~ /(..)(.)/;
                $tmp = ($1 > 18) ? $1 - 18 : $1 + 18;
                $tmp = '0'.$tmp if ($tmp < 10);
                $tmp .= 'L' if ($2 eq 'R');
                $tmp .= 'R' if ($2 eq 'L');
                $tmp .= 'C' if ($2 eq 'C');
                $rwy2 = $tmp; # get opp heading, but may NOT be per numbers
            }

            $isrwy = 1;

        } elsif ($type == 100) {
            $rwy1 = ${$ra}[8]; # 8  - 13L - Runway number (eg. 31R, 02). Leading zeros are required. Two to three characters. Valid suffixes: L, R or C (or blank)
            $lat1 = ${$ra}[9];
            $lon1 = ${$ra}[10];

            $rwy2 = ${$ra}[17];
            $lat2 = ${$ra}[18];
            $lon2 = ${$ra}[19];
            # get LENGTH
            fg_geo_inverse_wgs_84 ($lat1,$lon1,$lat2,$lon2,\$az1,\$az2,\$dist);
            $len = $SG_METER_TO_FEET * $dist;
            $rwlen2 = $len / 2;
            $hdg = $az1;
            $hdgr = $az2;
            # get CENTER
            fg_geo_direct_wgs_84( $lat1, $lon1, $hdg, $rwlen2, \$clat, \$clon, \$az2 );
            $isrwy = 1;
        }
        if ($isrwy) {
            my @a = ();
            $a[$RW_LEN] = $len; # lenght (feet)
            $a[$RW_HDG] = $hdg;
            $a[$RW_REV] = $hdgr;
            $a[$RW_TT1] = $rwy1;
            $a[$RW_TT2] = $rwy2;
            $a[$RW_CLAT] = $clat;
            $a[$RW_CLON] = $clon;
            $a[$RW_LLAT] = $lat1;
            $a[$RW_LLON] = $lon1;
            $a[$RW_RLAT] = $lat2;
            $a[$RW_RLON] = $lat2;
            $a[$RW_DONE] = 0;
            $tlat += $clat;
            $tlon += $clon;
            push(@runways,\@a);
        }
    }
    $cnt = scalar @runways;
    if (!$cnt) {
        pgm_exit(1,"Did NOT find any runways!\n");
    }
    if (!$gothdr) {
        pgm_exit(1,"Did NOT find airport header line!\n");
    }
    $alat = $tlat / $cnt;
    $alon = $tlon / $cnt;
    prt("Found $cnt runways at $icao, $name, ctr $alat,$alon...\n");
    my $off = -1;
    my $max_len = 0;
    $cnt = 0;
    foreach $ra (@runways) {
        $len = ${$ra}[$RW_LEN];
        if ($cnt == 0) {
            $off = $cnt;
            $max_len = $len;
        } elsif ($len > $max_len) {
            $off = $cnt;
            $max_len = $len;
        }
        $cnt++;
    }
    # get chosen ref array
    $ra = $runways[$off];
    $len = ${$ra}[$RW_LEN]; # lenght (feet)
    $hdg = ${$ra}[$RW_HDG];
    $hdgr = ${$ra}[$RW_REV];
    $rwy1 = ${$ra}[$RW_TT1];
    $rwy2 = ${$ra}[$RW_TT2];
    $clat = ${$ra}[$RW_CLAT];
    $clon = ${$ra}[$RW_CLON];
    #        $a[$RW_LLAT] = $lat1;
    #        $a[$RW_LLON] = $lon1;
    #        $a[$RW_RLAT] = $lat2;
    #        $a[$RW_RLON] = $lat2;
    #        $a[$RW_DONE] = 0;
    $dhdg = $hdg;
    set_int_stg(\$len);
    set_hdg_stg(\$hdg);
    set_hdg_stg(\$hdgr);
    set_int_stg(\$alt);
    set_decimal1_stg(\$dhdg);
    prt("At apt $icao, $name, alt $alt, ctr $alat,$alon\n"); 
    prt("Chosen rwy $rwy1, len $len ft, hdg $hdg ($dhdg), ctr $clat,$clon\n");
    prt("Opposite rwy $rwy2, len $len ft, hdg $hdgr\n");
    my @naids = ();
    my @circ  = ();
    my @a2 = ( $alat, $alon, \@naids, $ra, \@circ );
    my $refa = \@a2;    # [ $alat, $alon, \@naids, $ra, \@circ ];
    my $rl = get_locations();
    # key      0           1           2           3           4
    # ICAO => [ Center LAT, LON        NAVAIDS     RUNWAYS     CIRCUIT ]
    if (defined ${$rl}{$icao}) {
        prt("Updating existing entry...\n");
        my $pra = ${$rl}{$icao};
        my $plat = ${$pra}[0];
        my $plon = ${$pra}[1];
        my $prwys = ${$pra}[3];
        my $prcnt = scalar @{$prwys};
        my ($plen,$phdg,$prwy1,$pclat,$pclon,$paz1,$paz2,$pdist);
        $plen = ${$prwys}[0][$RW_LEN]; # lenght (feet)
        $phdg = ${$prwys}[0][$RW_HDG];
        ###$hdgr = ${$prwys}[0][$RW_REV];
        $prwy1 = ${$prwys}[0][$RW_TT1];
        ##$rwy2 = ${$prwys}[0][$RW_TT2];
        $pclat = ${$prwys}[0][$RW_CLAT];
        $pclon = ${$prwys}[0][$RW_CLON];
        fg_geo_inverse_wgs_84 ($clat,$clon,$pclat,$pclon,\$paz1,\$paz2,\$pdist);
        $pdist = get_dist_stg_nm($pdist);
        prt("Was $prwy1, len $plen, hdg $phdg, $pclat,$pclon, cnt $prcnt ($pdist)\n");
        ${$rl}{$icao} = $refa;
    } else {
        prt("Adding a NEW entry...\n");
        ${$rl}{$icao} = $refa;
    }

    pgm_exit(1,"TEMP EXIT\n");
}
#########################################
### MAIN ###
find_apt($apt_icao);    # "YGIL"
find_runway();
parse_args(@ARGV);
init_runway_array();
my $ref_hash = get_circuit_hash();
get_runways_and_pattern($ref_hash,'YGIL');
### is_in_circuit($ref_hash,$in_lat,$in_lon,"User input.",1); # 1
process_lat_lon($ref_hash,$in_lat,$in_lon,"1: User input"); # 1
if ($debug_on) {
    get_mid_tl2bl(\$in_lat,\$in_lon);
    process_lat_lon($ref_hash,$in_lat,$in_lon,"2: Mid tl2bl test point"); # 2
    #get_mid_bl2br(\$in_lat,\$in_lon);
    get_mid_bl2br2(\$in_lat,\$in_lon);
    process_lat_lon($ref_hash,$in_lat,$in_lon,"3: Mid bl2br test point"); # 3
    process_lat_lon($ref_hash,$def_lat2,$def_lon2,"4: 4th test point"); # 4
    process_lat_lon($ref_hash,$def_lat3,$def_lon3,"5: 5th test point"); # 5
    process_lat_lon($ref_hash,$def_lat4,$def_lon4,"6: 6th test point"); # 6
    process_lat_lon($ref_hash,$def_lat5,$def_lon5,"7: 7th test point"); # 7
    process_lat_lon($ref_hash,$def_lat6,$def_lon6,"8: 8th test point"); # 8
}

# create the image of all the points
paint_user_points($ref_hash,1);

pgm_exit(0,"");
########################################
sub give_help {
    prt("$pgmname: version $VER\n");
    prt("Usage: $pgmname [options] lat,lon\n");
    prt("Options:\n");
    prt(" --help (-h or -?) = This help, and exit 0.\n");
}
sub need_arg {
    my ($arg,@av) = @_;
    pgm_exit(1,"ERROR: [$arg] must have following argument!\n") if (!@av);
}

sub parse_args {
    my (@av) = @_;
    my ($arg,$sarg,@arr,$lat,$lon);
    while (@av) {
        $arg = $av[0];
        if ($arg =~ /^-/) {
            $sarg = substr($arg,1);
            $sarg = substr($sarg,1) while ($sarg =~ /^-/);
            if (($sarg =~ /^h/i)||($sarg eq '?')) {
                give_help();
                pgm_exit(0,"Help exit(0)");
            } else {
                pgm_exit(1,"ERROR: Invalid argument [$arg]! Try -?\n");
            }
        } else {
            @arr = split(",",$arg);
            if (scalar @arr == 2) {
                $lat = $arr[0];
                $lon = $arr[1];
            } else {
                pgm_exit(1,"ERROR: Invalid argument [$arg]! Try -?\n");
            }
        }
        shift @av;
    }

    if ((($in_lat == $bad_latlon)||($in_lon == $bad_latlon)) && $debug_on) {
        @arr = split(",",$def_latlon);
        if (scalar @arr == 2) {
            $lat = $arr[0];
            $lon = $arr[1];
            if (in_world_range($lat,$lon)) {
                $in_lat = $lat;
                $in_lon = $lon;
                prt("Set lat=$in_lat, lon=$in_lon\n");
            } else {
                pgm_exit(1,"ERROR: lat=$lat, lon=$lon OUT OF WORLD RANGE!\n");
            }
        } else {
            pgm_exit(1,"ERROR: Default is invlalid\n");
        }
    }
    if (($in_lat == $bad_latlon)||($in_lon == $bad_latlon)) {
        pgm_exit(1,"ERROR: No lat,lon found in command!\n");
    }
}

sub get_v1000_spec() {
    my $txt = <<EOF;
from : http://data.x-plane.com/file_specs/XP%20APT1000%20Spec.pdf
from : http://data.x-plane.com/file_specs/Apt810.htm
Code (apt.dat) Used for 
1  Airport header data. 
16 Seaplane base header data. No airport buildings or boundary fences will be rendered in X-Plane. 
17 Heliport header data.  No airport buildings or boundary fences will be rendered in X-Plane. 
10 Runway or taxiway at an airport. 
14 Tower view location. 
15 Ramp startup position(s) 
18 Airport light beacons (usually "rotating beacons" in the USA).  Different colours may be defined. 
19 Airport windsocks. 
50 to 56 Airport ATC (Air Traffic Control) frequencies. 
Runway line
#   0   1     2 3 4    5 6 7 8   9            10            11    12   13 14 15 16 17  18           19           20     21   22 23 24 25
EG: 100 29.87 3 0 0.00 0 0 0 16  -24.20505300 151.89156100  0.00  0.00 1  0  0  0  34  -24.19732300 151.88585300 0.00   0.00 1  0  0  0
OR: 100 29.87 1 0 0.15 0 2 1 13L 47.53801700  -122.30746100 73.15 0.00 2  0  0  1  31R 47.52919200 -122.30000000 110.95 0.00 2  0  0  1
Land Runway
0  - 100 - Row code for a land runway (the most common) 100
1  - 29.87 - Width of runway in metres Two decimal places recommended. Must be >= 1.00
2  - 3 - Code defining the surface type (concrete, asphalt, etc) Integer value for a Surface Type Code
3  - 0 - Code defining a runway shoulder surface type 0=no shoulder, 1=asphalt shoulder, 2=concrete shoulder
4  - 0.15 - Runway smoothness (not used by X-Plane yet) 0.00 (smooth) to 1.00 (very rough). Default is 0.25
5  - 0 - Runway centre-line lights 0=no centerline lights, 1=centre line lights
6  - 0 - Runway edge lighting (also implies threshold lights) 0=no edge lights, 2=medium intensity edge lights
7  - 1 - Auto-generate distance-remaining signs (turn off if created manually) 0=no auto signs, 1=auto-generate signs

Airport Line. eg '1 5355 1 0 KABQ Albuquerque Intl Sunport'
1	   - this as an airport header line. 16 is a seaplane/floatplane base, 17 a heliport.
5355   - Airport elevation (in feet above MSL).  
1	   - Airport has a control tower (1=yes, 0=no).
0	   - Display X-Plane’s default airport buildings (1=yes, 0=no).
KABQ   - Identifying code for the airport (the ICAO code, if one exists).
Albuquerque Intl Sunport - Airport name.

EOF
    return $txt;
}

# eof - solve.pl
