#!/usr/bin/perl -w
# NAME: slippydirs.pl
# AIM: Just some 'tests' on slippy directories, and other slippy map things
# 18/07/2020 - review
# 11/07/2016 - Review mainly of output
# 15/05/2014 - Added /5/23/13/png as another way to enter 'lat lon zoom'
# 15/04/2014 geoff mclane http://geoffair.net/mperl
use strict;
use warnings;
use File::Basename;  # split path ($name,$dir,$ext) = fileparse($file [, qr/\.[^.]*/] )
use Math::Trig;
use Cwd;
my $os = $^O;
my $perl_dir = '/home/geoff/bin';
my $PATH_SEP = '/';
my $temp_dir = '/tmp';
if ($os =~ /win/i) {
    $perl_dir = 'C:\GTools\perl';
    $temp_dir = $perl_dir;
    $PATH_SEP = "\\";
}
unshift(@INC, $perl_dir);
require 'lib_utils.pl' or die "Unable to load 'lib_utils.pl' Check paths in \@INC...\n";
require 'fg_wsg84.pl' or die "Unable to load fg_wsg84.pl ...\n";
# log file stuff
our ($LF);
my $pgmname = $0;
if ($pgmname =~ /(\\|\/)/) {
    my @tmpsp = split(/(\\|\/)/,$pgmname);
    $pgmname = $tmpsp[-1];
}
my $outfile = $temp_dir.$PATH_SEP."temp.$pgmname.txt";
open_log($outfile);

# user variables
my $VERS = "0.0.5 2020-07-18";
###my $VERS = "0.0.4 2016-07-11";
###my $VERS = "0.0.3 2014-05-15";
###my $VERS = "0.0.2 2014-04-15";
my $load_log = 0;
my $in_file = '';
my $verbosity = 0;
my $out_file = '';
my $slippy_dir = '';
my $list_touching_tiles = 0;

my ($in_lat,$in_lon,$in_zoom);

my $MAX_SLIPPY_LAT = 85.0511;
my $MIN_SLIPPY_LAT = -85.0511;
my $MAX_SLIPPY_LON = 180.0;
my $MIN_SLIPPY_LON = -180.0;

my $srtm1_dir = "F:\\data\\dem1\\srtm1";
my $ferr3_dir  = "F:\\data\\dem3\\hgt";
my $ferr30_dir = "F:\\data\\srtm30_plus"; # Bathymetry.srtm

my $root_dir = 'F:\FGx\fgx-terrain';

my ($g_min_lat,$g_max_lat,$g_min_lon,$g_max_lon);


# ### DEBUG ###
my $debug_on = 0;
my $def_file = 'def_file';

### program variables
my @warnings = ();
my $cwd = cwd();
my $SG_EPSILON = 0.000001;

sub VERB1() { return $verbosity >= 1; }
sub VERB2() { return $verbosity >= 2; }
sub VERB5() { return $verbosity >= 5; }
sub VERB9() { return $verbosity >= 9; }

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" ) if (VERB9());
    }
}

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);
}

#           prefix, min.lat, max.lat, min.lon max.lon
#           0       1        2        3       4
my @a1 = ( "w180n90", 40, 90, -180, -140, -22, 6098, 448, 482 );
my @a2 = ( "w140n90", 40, 90, -140, -100, -108, 4635, 731, 596 );
my @a3 = ( "w100n90", 40, 90, -100, -60, -35, 2416, 337, 280 );
my @a4 = ( "w060n90", 40, 90, -60, -20, -13, 3940, 1626, 932 );
my @a5 = ( "w020n90", 40, 90, -20, 20, -179, 4536, 402, 426 );
my @a6 = ( "e020n90", 40, 90, 20, 60, -188, 5472, 213, 312 );
my @a7 = ( "e060n90", 40, 90, 60, 100, -156, 7169, 509, 697 );
my @a8 = ( "e100n90", 40, 90, 100, 140, -110, 3901, 596, 455 );
my @a9 = ( "e140n90", 40, 90, 140, 180, -26, 4578, 415, 401 );
my @a10 = ( "w180n40", -10, 40, -180, -140, -3, 4120, 832, 860 );
my @a11 = ( "w140n40", -10, 40, -140, -100, -174, 4228, 1322, 745 );
my @a12 = ( "w100n40", -10, 40, -100, -60, -171, 6543, 367, 609 );
my @a13 = ( "w060n40", -10, 40, -60, -20, -22, 2504, 217, 160 );
my @a14 = ( "w020n40", -10, 40, -20, 20, -138, 3958, 438, 298 );
my @a15 = ( "e020n40", -10, 40, 20, 60, -422, 5778, 724, 557 );
my @a16 = ( "e060n40", -10, 40, 60, 100, -46, 8685, 1807, 1889 );
my @a17 = ( "e100n40", -10, 40, 100, 140, -147, 7213, 690, 911 );
my @a18 = ( "e140n40", -10, 40, 140, 180, -42, 4650, 530, 728 );
my @a19 = ( "w180s10", -60, -10, -180, -140, -41, 1784, 191, 294 );
my @a20 = ( "w140s10", -60, -10, -140, -100, -5, 910, 79, 133 );
my @a21 = ( "w100s10", -60, -10, -100, -60, -752, 6813, 1080, 1359 );
my @a22 = ( "w060s10", -60, -10, -60, -20, -127, 2823, 411, 294 );
my @a23 = ( "w020s10", -60, -10, -20, 20, -24, 2498, 1088, 404 );
my @a24 = ( "e020s10", -60, -10, 20, 60, -26, 3408, 889, 453 );
my @a25 = ( "e060s10", -60, -10, 60, 100, -3, 2557, 251, 262 );
my @a26 = ( "e100s10", -60, -10, 100, 140, -33, 1360, 290, 172 );
my @a27 = ( "e140s10", -60, -10, 140, 180, -43, 3119, 278, 265 );
my @a28 = ( "w180s60", -90, -60, -180, -120 );
my @a29 = ( "w120s60", -90, -60, -120, -60 );
my @a30 = ( "e060s60", -90, -60, -60, 0 );
my @a31 = ( "w000s60", -90, -60, 0, 60 );
my @a32 = ( "w060s60", -90, -60, 60, 120 );
my @a33 = ( "e120s60", -90, -60, 120, 180 );

my @srtm30_tile_table = ( \@a1, \@a2, \@a3, \@a4, \@a5, \@a6, \@a7, \@a8, \@a9, \@a10,
    \@a11, \@a12, \@a13, \@a14, \@a15, \@a16, \@a17, \@a18, \@a19, \@a20,
    \@a21, \@a22, \@a23, \@a24, \@a25, \@a26, \@a27, \@a28, \@a29, \@a30,
    \@a31, \@a32, \@a33 );

sub getTileNumber($$$) {
    my ($lat,$lon,$zoom) = @_;
    my $z2 = 2 ** $zoom;
    my $xtile = int( ($lon+180)/360 * $z2 ) ;
    my $ytile = int( (1 - log(tan(deg2rad($lat)) + sec(deg2rad($lat)))/pi)/2 * $z2 ) ;
    return ($xtile, $ytile);
}

sub getLonLat {
	my ($xtile, $ytile, $zoom) = @_;
	my $n = 2 ** $zoom;
	my $lon_deg = $xtile / $n * 360.0 - 180.0;
	my $lat_deg = rad2deg(atan(sinh(pi * (1 - 2 * $ytile / $n))));
	return ($lon_deg, $lat_deg);
}

# https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames
# convert from permalink OSM format like:
# http://www.openstreetmap.org/?lat=43.731049999999996&lon=15.79375&zoom=13&layers=M
# to OSM "Export" iframe embedded bbox format like:
# http://www.openstreetmap.org/export/embed.html?bbox=15.7444,43.708,15.8431,43.7541&layer=mapnik
#
# BBOX left lon, bottom lat, right lon, top lat 
# Zoom 0 - Span x 360 y 180 = whole world
# Zoom 1 - 2 x 2 y - Span x 180 y 90 - 0/0.png -180,0,0,90  - 1/0.png 0,0,180,90
#                                    - 0/1.png -180,-90,0,0 - 1/1.png 0,-90,180,0
# Zoom 2 - 4 x 4 y 
# 2 0 3 -225,-40.9798980696201,-135,-79.1713346408195 - span x 90 y 38.1914 lon -180 lat -85
# 2 0 2 -225,40.9798980696201,-135,-40.9798980696201 - span x 90 y 81.9598 lon -180 lat -66.5
# 2 0 1 -225,79.1713346408195,-135,40.9798980696201 - span x 90 y 38.1914 lon -180 lat 0.5
# 2 0 0 -225,87.7425091586603,-135,79.1713346408195 - span x 90 y 8.5712 lon -180 lat 67
# ...
# Zoom  6: got -
#   36 xtiles 0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,1...,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35
#   64 ytiles 0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,1...,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63
# 6 0 63 -182.8125,-84.2671724043167,-177.1875,-84.8024737243345 - span x 5.625 y 0.5353 lon -180 lat -85
# 6 0 62 -182.8125,-83.6769430484155,-177.1875,-84.2671724043167 - span x 5.625 y 0.5902 lon -180 lat -84.5
#     ---
# 6 0 28 -182.8125,24.5271348225978,-177.1875,19.3111433550646 - span x 5.625 y 5.216 lon -180 lat 17
#     ---
# 6 0 1 -182.8125,84.8024737243345,-177.1875,84.2671724043167 - span x 5.625 y 0.5353 lon -180 lat 84
# 6 0 0 -182.8125,85.287916121237,-177.1875,84.8024737243345 - span x 5.625 y 0.4854 lon -180 lat 85
# ...
# Zoom 19: got -
#  360 xtiles 0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,1...8,349,350,351,352,353,354,355,356,357,358,359
#  341 ytiles 858,8822,16095,22787,28985,34758,40160,45236,...127,489529,495302,501500,508192,515465,523429
# 19 0 523429 -180.000343322754,-84.9999544617033,-179.999656677246,-85.0000143069889 - span x 0.0007 y 0.0001 lon -180 lat -85
#       ---
# 19 0 262872 -180.000343322754,-0.499528278561968,-179.999656677246,-0.50021489793785 - span x 0.0007 y 0.0007 lon -180 lat -0.5
#       ---
# 19 0 8822 -180.000343322754,84.5000882101567,-179.999656677246,84.500022398761 - span x 0.0007 y 0.0001 lon -180 lat 84.5
# 19 0 858 -180.000343322754,85.0000741515601,-179.999656677246,85.0000143069889 - span x 0.0007 y 0.0001 lon -180 lat 85
#

sub LatLon_to_bbox {
	my ($lat, $lon, $zoom) = @_;
 
	#my $width = 425; my $height = 350;	# note: must modify this to match your embed map width/height in pixels
	my $width = 256; my $height = 256;	# note: must modify this to match your embed map width/height in pixels
	my $tile_size = 256;
 
	my ($xtile, $ytile) = getTileNumber ($lat, $lon, $zoom);
    # what is the operator pecedence here??????????
	my $xtile_s = ($xtile * $tile_size - $width/2) / $tile_size;
	my $ytile_s = ($ytile * $tile_size - $height/2) / $tile_size;
	my $xtile_e = ($xtile * $tile_size + $width/2) / $tile_size;
	my $ytile_e = ($ytile * $tile_size + $height/2) / $tile_size;
 
	my ($lon_s, $lat_s) = getLonLat($xtile_s, $ytile_s, $zoom);
	my ($lon_e, $lat_e) = getLonLat($xtile_e, $ytile_e, $zoom);
 
	my $bbox = "$lon_s,$lat_s,$lon_e,$lat_e";
	return $bbox;
}

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 num_sort {
    return -1 if ($a < $b);
    return 1 if ($a > $b);
    return 0;
}

sub max_line($$) {
    my ($ln,$max) = @_;
    my $len = length($ln);
    if ($len < $max) {
        return $ln;
    }
    my $h = $max / 2;
    my $nl = substr($ln,0,$h);
    $nl .= '...';
    $nl .= substr($ln,($len - $h));
    return $nl;
}

# the /zoom/x/y.png tile name
sub do_test() {
    my ($lat,$lon,$num,$z,$cnt,@arr,$pnum,$err);
    my ($cz,$ccnt,$line,$dir,$ok,$chked,$fnd,$chkedz,$fndz);
    my ($ok2,$x,$y,$yline,$bbox,$cycnt,$ra);
    my ($xs,$ys);
    $z = 0;
    my %zooms = ();
    $chked = 0;
    $fnd = 0;
    $chkedz = 0;
    $fndz = 0;
    my @bboxes = ();
    for ($z = 0; $z < 20; $z++) {
    #for ($z = 0; $z < 8; $z++) {
        my %nums = ();
        $dir = $root_dir.$PATH_SEP.$z;
        $ok = 'NF';
        $chkedz++;
        if (-d $dir) {
            $ok = 'ok';
            $fndz++;
        }
        my %ynums = ();
        for ($lon = -180; $lon < 180; $lon++) {
            $num = int(($z**2)*(($lon + 180) / 360));
            if (!defined $nums{$num}) {
                if ($ok eq 'ok') {
                    $chked++;
                    $ok2 = 'NF';
                    $dir = $root_dir.$PATH_SEP.$z.$PATH_SEP.$num;
                    if (-d $dir) {
                        $ok2 = 'ok';
                        $fnd++;
                    }
                }
                $nums{$num} = $ok;
                ##prt("z=$z: from lon $lon, got number $num\n");
            }
            for ($lat = -85; $lat <= 85; $lat += 0.5) {
                ($x,$y) = getTileNumber($lat,$lon,$z);
                if (!defined $ynums{$y}) {
                    ##prt("For zoom $z, lat $lat, lon $lon, xtile $x, ytile $y\n");
                    $ynums{$y} = 1;
                    $bbox = LatLon_to_bbox( $lat, $lon, $z );
                    push(@bboxes, [$lat, $lon, $z, $bbox, $x, $y]);

                }
            }
        }
        ##############################################
        # for this ZOOM
        ##############################################
        @arr = sort num_sort keys %nums;
        $cnt = scalar @arr;
        $cz = sprintf("%2d",$z);
        $ccnt = sprintf("%4d",$cnt);
        $line = max_line(join(",",@arr),90);
        @arr = sort num_sort keys %ynums;
        $cnt = scalar @arr;
        $cycnt = sprintf("%4d",$cnt);
        $yline = max_line(join(",",@arr),90);

        prt("Zoom $cz: got -\n $ccnt xtiles $line\n $cycnt ytiles $yline\n");
        $cnt = 0;
        $pnum = 0;
        $err = 0;
        @arr = sort num_sort keys %nums;
        foreach $num (@arr) {
            if ($num != $cnt) {
                $err++;
                prt("$err: Not sequention $pnum $num $cnt\n");
            }
            $cnt++;
            $pnum = $num;
        }
        $zooms{$z} = \%nums;
        ##############################################
    }
    prt("Of $chkedz zooms checked, found $fndz. Of $chked checked, found $fnd\n");
    $cnt = scalar @bboxes;
    prt("Got $cnt bounding boxes...\n");
    #               0     1     2   3
    #push(@bboxes, [$lat, $lon, $z, $bbox]);
    $cnt = 0;
    foreach $ra (@bboxes) {
        $lat = ${$ra}[0];
        $lon = ${$ra}[1];
        $z   = ${$ra}[2];
        $bbox = ${$ra}[3];
        $x   = ${$ra}[4];
        $y   = ${$ra}[5];
        if ($z != $cnt) {
            prt("\n");
            $cnt = $z;
        }
        @arr = split(",",$bbox);
        $xs = ($arr[2] < $arr[0]) ? $arr[0] - $arr[2] : $arr[2] - $arr[0];
        $ys = ($arr[3] < $arr[1]) ? $arr[1] - $arr[3] : $arr[3] - $arr[1];
        $xs = int(($xs + 0.000005) * 100000) / 100000;
        $ys = int(($ys + 0.000005) * 100000) / 100000;
        #prt("$z $lat $lon $bbox\n"); 
        prt("$z $x $y $bbox - span x $xs y $ys lon $lon lat $lat\n"); 
    }
    $load_log = 1;
}

sub do_test2() {
    my ($z,$lat,$lon,$tx,$ty,$clat,$tc,$ctx,$cty,$cz,$ctc);
    my ($lonsp,$latsp,$clonsp,$clatsp,$bbox,$i);
    my $minlat = -85.0511;
    my $maxlon = 180;
    my $latmax = 85.0511 * 2;
    my $lonmax = 360;
    # Mt Everest
    my $tlat = 27.988;
    my $tlon = 86.9253;
    my @bboxes = ();
    for ($z = 0; $z < 20; $z++) {
        $cz = sprintf("%2d",$z);
        $tc = (2 ** $z) * (2 ** $z);
        $ctc = sprintf("%12d",$tc);
        ($lon,$lat) = getLonLat(0,0,$z);
        $clat = sprintf("%.4f",$lat);
        ($tx,$ty) = getTileNumber($minlat,$maxlon,$z);
        $ty++;
        $lonsp = $lonmax / $tx;
        $latsp = $latmax / $ty;
        $ctx = sprintf("%6d",$tx);
        $cty = sprintf("%6d",$ty);
        $clonsp = sprintf("%10.6f",$lonsp);
        $clatsp = sprintf("%10.6f",$latsp);
        prt("Z=$cz x,y 0,0 lon,lat $lon,$clat -- $ctx,$cty $maxlon,$minlat $ctc $clonsp,$clatsp\n");
        $bbox = LatLon_to_bbox( $tlat, $tlon, $z );
        push(@bboxes,$bbox);
    }
    # without format sucks
    ##prt(join("\n",@bboxes)."\n");
    my (@arr);
    foreach $bbox (@bboxes) {
        @arr = split(",",$bbox);
        for ($i = 0; $i < 4; $i++) {
            $arr[$i] = sprintf("%14.9f",$arr[$i]);
        }
        prt(join(",",@arr)."\n");
    }

    $load_log = 1;
}

sub get_double($) {
    my $d = shift;
    my $cd = sprintf("%4.10f",$d);
    return $cd;
}

my %hgt_shown = ();
# HGT file format :
#  Its south western corner can be deduced from its file name: 
# for example, n51e002.hgt covers the area between N 51 E 2 and N 52 E 3, 
# and s14w077.hgt covers S 14 W 77 to S 13 W 76
#
sub show_hgt($$) {
    my ($inlat,$inlon) = @_;
    my $lat = $inlat;
    my $lon = $inlon;
    my ($ilat,$ilon);
    my ($nlat,$nlon);
    #  Its south western corner can be deduced from its file name: 
    # for example, n51e002.hgt covers the area between N 51 E 2 and N 52 E 3, 
    # and s14w077.hgt covers S 14 W 77 to S 13 W 76 -14,-77 to -13,-76
    my $ns = 'N';
    my $ew = 'E';
    if ($lat < 0.0) {
        $ns = 'S';
        $lat *= -1;
        $ilat = int($lat);
        if ($lat - $ilat) {
            $lat += 1;
        }
    }
    if ($lon < 0.0) {
        $ew = 'W';
        $lon *= -1;
        $ilon = int($lon);
        if ($lon - $ilon) {
            $lon += 1;
        }
    }
    $ilon = int($lon);
    $ilat = int($lat);
    my $file = $ns.sprintf("%02d",$ilat).$ew.sprintf("%03d",$ilon).".hgt";
    if ($inlat < 0.0) {
        $nlat = $ilat - 1;
    } else {
        $nlat = $ilat + 1;
    }
    if ($inlon < 0.0) {
        $nlon = $ilon - 1;
    } else {
        $nlon = $ilon + 1;
    }
    if (!defined $hgt_shown{$file}) {
        $hgt_shown{$file} = 1;
        prt("HGT file $file Range BL $ns $ilat $ew $ilon to TR $ns $nlat $ew $nlon\n");
    }
}

sub get_hgt_file($$) {
    my ($inlat,$inlon) = @_;
    my $lat = $inlat;
    my $lon = $inlon;
    my ($ilat,$ilon);
    my ($nlat,$nlon);
    #  Its south western corner can be deduced from its file name: 
    # for example, n51e002.hgt covers the area between N 51 E 2 and N 52 E 3, 
    # and s14w077.hgt covers S 14 W 77 to S 13 W 76 -14,-77 to -13,-76
    my $ns = 'N';
    my $ew = 'E';
    if ($lat < 0.0) {
        $ns = 'S';
        $lat *= -1;
        $ilat = int($lat);
        if ($lat - $ilat) {
            $lat += 1;
        }
    }
    if ($lon < 0.0) {
        $ew = 'W';
        $lon *= -1;
        $ilon = int($lon);
        if ($lon - $ilon) {
            $lon += 1;
        }
    }
    $ilon = int($lon);
    $ilat = int($lat);
    my $file = $ns.sprintf("%02d",$ilat).$ew.sprintf("%03d",$ilon).".hgt";
    if ($inlat < 0.0) {
        $nlat = $ilat - 1;
    } else {
        $nlat = $ilat + 1;
    }
    if ($inlon < 0.0) {
        $nlon = $ilon - 1;
    } else {
        $nlon = $ilon + 1;
    }
    $g_min_lat = $ilat;
    $g_max_lat = $nlat;
    $g_min_lon = $ilon;
    $g_max_lon = $nlon;
    if ($inlat < 0.0) {
        $g_min_lat *= -1;
        $g_max_lat *= -1;
    }
    if ($inlon < 0.0) {
        $g_min_lon *= -1;
        $g_max_lon *= -1;
    }
    return $file;
}


#typedef struct tagSRTM30 {
#    const char *tile; 0
my $lat_min = 1;
my $lat_max = 2;
my $lon_min = 3;
my $lon_max = 4;
#    int elev_min;
#    int elev_max;
#    int mean;
#    int std_dev;
#} SRTM30, *PSRTM30;

sub get_srtm30_tile($$) {
    my ($lat,$lon) = @_;
    my ($ra);
    foreach $ra (@srtm30_tile_table) {
        if (($lat <= ${$ra}[$lat_max]) &&
            ($lat >  ${$ra}[$lat_min]) &&
            ($lon <  ${$ra}[$lon_max]) &&
            ($lon >= ${$ra}[$lon_min])) {
            return ${$ra}[0];
        }
    }
    return "not found";
}
sub get_srtm30_ra($$$) {
    my ($lat,$lon,$rra) = @_;
    my ($ra);
    foreach $ra (@srtm30_tile_table) {
        if (($lat <= ${$ra}[$lat_max]) &&
            ($lat >  ${$ra}[$lat_min]) &&
            ($lon <  ${$ra}[$lon_max]) &&
            ($lon >= ${$ra}[$lon_min])) {
            ${$rra} = $ra;
            return 1;
        }
    }
    return 0;   # "not found";
}

###############################################################
### from say cfcsvlogs.pl
sub get_decimal_stg($$$) {
    my ($dec,$il,$dl) = @_;
    my (@arr);
    if ($dec =~ /\./) {
        @arr = split(/\./,$dec);
        if (scalar @arr == 2) {
            $arr[0] = " ".$arr[0] while (length($arr[0]) < $il);
            $dec = $arr[0];
            if ($dl > 0) {
                $dec .= ".";
                $arr[1] = substr($arr[1],0,$dl) if (length($arr[1]) > $dl);
                $dec .= $arr[1];
            }
        }
    } else {
        $dec = " $dec" while (length($dec) < $il);
        if ($dl) {
            $dec .= ".";
            while ($dl--) {
                $dec .= "0";
            }
        }
    }
    return $dec;
}

sub get_sg_dist_stg($) {
    my ($sg_dist) = @_;
    my $sg_km = $sg_dist / 1000;
    my $sg_im = int($sg_dist);
    my $sg_ikm = int($sg_km + 0.5);
    my $dlen = 5;
    # if (abs($sg_pdist) < $CP_EPSILON)
    my $sg_dist_stg = "";
    if (abs($sg_km) > $SG_EPSILON) { # = 0.0000001; # EQUALS SG_EPSILON 20101121
        if ($sg_ikm && ($sg_km >= 1)) {
            $sg_km = int(($sg_km * 10) + 0.05) / 10;
            #$sg_dist_stg .= get_decimal_stg($sg_km,5,1)." km";
            $sg_dist_stg .= get_decimal_stg($sg_km,($dlen - 2),1)." km";
        } else {
            #$sg_dist_stg .= "$sg_im m, <1km";
            #$sg_dist_stg .= get_decimal_stg($sg_im,7,0)." m.";
            $sg_dist_stg .= get_decimal_stg($sg_im,$dlen,0)." m.";
        }
    } else {
        #$sg_dist_stg .= "0 m";
        #$sg_dist_stg .= get_decimal_stg('0',7,0)." m.";
        $sg_dist_stg .= get_decimal_stg('0',$dlen,0)." m.";
    }
    return $sg_dist_stg;
}


###############################################################

sub print_slippy_dir($$$) {
    my ($lat,$lon,$z) = @_;
    my ($x, $y) = getTileNumber($lat,$lon,$z);
    my $dir = "/$z/$x/$y.png";
    prt("For lat,lon,zoom ".get_double($lat).",".get_double($lon).",$z, dir is [$dir], y,x is $y,$x\n");
}

# // { "Mont Blanc", 45.832704, 6.864797, 4810 }, // N45E006.hgt
#   370 KKIC MESA DEL REY             36.228056000,-121.121945000 83 km on 71.9 d.
sub show_values() {
    my $lat = $in_lat;
    my $lon = $in_lon;
    my $z   = $in_zoom;

    my ($maxx,$maxy) = getTileNumber($MIN_SLIPPY_LAT,$MAX_SLIPPY_LON,$z);
    my ($minx,$miny) = getTileNumber($MAX_SLIPPY_LAT,$MIN_SLIPPY_LON,$z);

    my ($x, $y) = getTileNumber($lat,$lon,$z);
    my $dir = "/$z/$x/$y.png";
    my ($tl_lon, $tl_lat) = getLonLat($x+0,$y+0,$z);
    my ($bl_lon, $bl_lat) = getLonLat($x+0,$y+1,$z);
    my ($tr_lon, $tr_lat) = getLonLat($x+1,$y+0,$z);
    my ($br_lon, $br_lat) = getLonLat($x+1,$y+1,$z);
    my $ctr_lat = ($tl_lat + $br_lat) / 2;
    my $ctr_lon = ($tl_lon+$br_lon) / 2;
    my ($cx, $cy) = getTileNumber($ctr_lat,$ctr_lon,$z);
    prt("\n");
    prt("For lat,lon,zoom $lat,$lon,$z, dir is [$dir], y,x is $y,$x ($miny,$minx,$maxy,$maxx)\n");
    prt("CENTER:  lat,lon ".get_double($ctr_lat).",".get_double($ctr_lon));
    ##prt(" (".get_double(($bl_lat + $tr_lat) / 2).",".get_double(($bl_lon+$tr_lon) / 2));
    prt(" /z/cx/cy.png /$z/$cx/$cy.png\n");
    prt("\n");

    #    TL 27.6835280838,86.6601562500 y+0,x+0 3034,1720
    prt("   latitude,     longitude,            x, y\n");
    prt("TL ".get_double($tl_lat).",".get_double($tl_lon)." y+0,x+0 ".($x+0).",".($y+0)."\n");
    prt("BL ".get_double($bl_lat).",".get_double($bl_lon)." y+1,x+0 ".($x+0).",".($y+1)."\n");
    prt("TR ".get_double($tr_lat).",".get_double($tr_lon)." y+0,x+1 ".($x+1).",".($y+0)."\n");
    prt("BR ".get_double($br_lat).",".get_double($br_lon)." y+1,x+1 ".($x+1).",".($y+1)."\n");

    prt("\n");
    my $xspan = ($tr_lon - $tl_lon);    # width at top
    my $yspan = ($tl_lat - $bl_lat);    # height at top
    # get WIDTH and HEIGHT
    my ($res,$az1,$az2,$width,$height,$cwid,$chgt);
    $res = fg_geo_inverse_wgs_84( $tl_lat, $tl_lon, $tl_lat, $tr_lon, \$az1, \$az2, \$width );
    $res = fg_geo_inverse_wgs_84( $tl_lat, $tl_lon, $bl_lat, $tl_lon, \$az1, \$az2, \$height );
    $cwid = trim_all(get_sg_dist_stg($width));
    $chgt = trim_all(get_sg_dist_stg($height));

    #my $xas = $xspan * 3600;
    #my $yas = $yspan * 3600;
    #my $xas = $xspan * 1201 / 3600;
    #my $yas = $yspan * 1201 / 3600;
    my $xas = 1201 / 3600;
    my $yas = 3601 / 3600;
    prt("SPAN x=".get_double($xspan)." y=".get_double($yspan)." degs".
        " width=$cwid height=$chgt".
        " (3=".get_double($xas).",1=".get_double($yas).")\n");
    #prt("SPAN x = $xspan - y = $yspan degrees\n");
    prt("\n");
    prt("bbox(southLatitude,westLongitude,northLatitude,eastLongitude)\n");
    prt("BBOX=$bl_lat,$bl_lon,$tr_lat,$tr_lon (BLTR)\n");
    prt("\n");
    prt("osm bbox=westLongitude,southLatitude,eastLongitude,northLatitude)\n");
    prt("bbox=$bl_lon,$bl_lat,$tr_lon,$tr_lat (LBRT)\n");

    # HGT file name
    show_hgt($lat,$lon);
    show_hgt($bl_lat,$bl_lon);
    show_hgt($tr_lat,$tr_lon);

    my ($ra,$tile,$msg,$file);
    my ($min_lat,$max_lat,$min_lon,$max_lon);
    prt("\nGiven lat,lon $lat,$lon the DEM is :-\n");
    if (get_srtm30_ra($lat,$lon,\$ra)) {
        #           prefix, min.lat, max.lat, min.lon max.lon
        #           0       1        2        3       4
        #my  = ( "w180n90", 40,     90,     -180,   -140, -22, 6098, 448, 482 );
        $tile = ${$ra}[0];
        $min_lat = ${$ra}[1];
        $max_lat = ${$ra}[2];
        $min_lon = ${$ra}[3];
        $max_lon = ${$ra}[4];

        $file = $tile.".Bathymetry.srtm";
        $msg = "N/A";
        if (-d $ferr30_dir) {
            $msg = "NF";
            $file = $ferr30_dir.'\\'.$file;
            if (-f $file) {
                $msg = 'ok';
            }
        }
        $msg .= " lat,lon min $min_lat,$min_lon max $max_lat,$max_lon";
        prt("DEM30: $file $msg\n");
        $file = get_hgt_file($lat,$lon);
        $msg = "N/A";
        if (-d $ferr3_dir) {
            $msg = "NF";
            $file = $ferr3_dir.'\\'.$file;
            if (-f $file) {
                $msg = 'ok';
            }
        }
        $msg .= " lat,lon min $g_min_lat,$g_min_lon max $g_max_lat,$g_max_lon";
        prt("DEM 3: $file $msg\n");
        $file = get_hgt_file($lat,$lon);
        $msg = "N/A";
        if (-d $srtm1_dir) {
            $msg = "NF";
            $file = $srtm1_dir.'\\'.$file;
            if (-f $file) {
                $msg = 'ok';
            }
        }
        $msg .= " lat,lon min $g_min_lat,$g_min_lon max $g_max_lat,$g_max_lon";
        prt("DEM 1: $file $msg\n");

    } else {
        $tile = get_srtm30_tile($lat,$lon);
        prt("$tile".".bathymetric.srtm\n");
    }

    if ($list_touching_tiles) {
        prt("\nListing of 9 tiles, forming a block...\n");
        my $TL_lat = $lat - $yspan;
        my $TL_lon = $lon - $xspan;
        prt("Print top 3 tiles...\n");
        print_slippy_dir($TL_lat,$TL_lon,$z);
        $TL_lon += $xspan;
        print_slippy_dir($TL_lat,$TL_lon,$z);
        $TL_lon += $xspan;
        print_slippy_dir($TL_lat,$TL_lon,$z);

        prt("Print center 3 tiles...\n");
        # reset lon, adjust lat
        $TL_lon = $lon - $xspan;
        $TL_lat += $yspan;
        print_slippy_dir($TL_lat,$TL_lon,$z);
        $TL_lon += $xspan;
        print_slippy_dir($TL_lat,$TL_lon,$z);
        $TL_lon += $xspan;
        print_slippy_dir($TL_lat,$TL_lon,$z);
        prt("Print bottom 3 tiles...\n");
        # reset lon, adjust lat
        $TL_lon = $lon - $xspan;
        $TL_lat += $yspan;
        print_slippy_dir($TL_lat,$TL_lon,$z);
        $TL_lon += $xspan;
        print_slippy_dir($TL_lat,$TL_lon,$z);
        $TL_lon += $xspan;
        print_slippy_dir($TL_lat,$TL_lon,$z);

        prt("\n");
    }
}

sub show_values_OK() {
    my $lat = $in_lat;
    my $lon = $in_lon;
    my $z   = $in_zoom;
    my ($x, $y) = getTileNumber($lat,$lon,$z);
    my $dir = "$z/$x/$y.png";
    prt("For lat,lon,zoom $lat,$lon,$z, dir is [$dir], y,x is $y,$x\n");
    prt("\n");
    #                                  x    y    z
    my ($tl_lon, $tl_lat) = getLonLat($x+0,$y+0,$z);
    prt("TL y+0,x+0 ".($x+0).",".($y+0)." $tl_lat,$tl_lon\n");
    my ($bl_lon, $bl_lat) = getLonLat($x+0,$y+1,$z);
    prt("BL y+1,x+0 ".($x+0).",".($y+1)." $bl_lat,$bl_lon\n");
    my ($tr_lon, $tr_lat) = getLonLat($x+1,$y+0,$z);
    prt("TR y+0,x+1 ".($x+1).",".($y+0)." $tr_lat,$tr_lon\n");
    my ($br_lon, $br_lat) = getLonLat($x+1,$y+1,$z);
    prt("BR y+1,x+1 ".($x+1).",".($y+1)." $br_lat,$br_lon\n");
    prt("\n");
    prt("SPAN x = ".($tr_lon - $tl_lon)." y = ".($tl_lat - $bl_lat)."\n");
    prt("\n");
    prt("bbox(southLatitude,westLongitude,northLatitude,eastLongitude)\n");
    prt("BBOX=$bl_lat,$bl_lon,$tr_lat,$tr_lon (BLTR)\n");
    prt("\n");
    prt("osm bbox=westLongitude,southLatitude,eastLongitude,northLatitude)\n");
    prt("bbox=$bl_lon,$bl_lat,$tr_lon,$tr_lat (LBRT)\n");

}


sub do_test3() {
    # center
    my $lat = 51.562;
    my $lon = -1.747;
    my $z = 7;
    my ($xtile, $ytile) = getTileNumber($lat,$lon,$z);
    my $dir = "$z/$xtile/$ytile.png";
    prt("lat,lon,z $lat,$lon,$z = $dir ie y,x $ytile,$xtile\n");
    prt("How to get the extent of that tile?\n");
    my ($lon_deg, $lat_deg);
    my ($x,$y);
    for ($y = $ytile - 1; $y <= $ytile + 1; $y++) {
        for ($x = $xtile - 1; $x <= $xtile + 1; $x++) {
            ($lon_deg, $lat_deg) = getLonLat($x,$y,$z);
            prt("y,x $y,$x lat,lon $lat_deg,$lon_deg\n");
        }
    }

    prt("Oxford (appr) lat,lon,z $lat,$lon,$z = $dir is y,x $ytile,$xtile\n");
    my ($tl_lon, $tl_lat) = getLonLat(63,42,$z);
    prt("TL y+0,x+0 42,63 $tl_lat,$tl_lon\n");
    my ($bl_lon, $bl_lat) = getLonLat(63,43,$z);
    prt("BL y+1,x+0 43,63 $bl_lat,$bl_lon\n");
    my ($tr_lon, $tr_lat) = getLonLat(64,42,$z);
    prt("TR y+0,x+1 42,64 $tr_lat,$tr_lon\n");
    my ($br_lon, $br_lat) = getLonLat(64,43,$z);
    prt("BR y+1,x+1 43,64 $br_lat,$br_lon\n");
    prt("x span ".($tr_lon - $tl_lon)." y span ".($tl_lat - $bl_lat)."\n");
    prt("x span ".($br_lon - $bl_lon)." y span ".($tr_lat - $br_lat)."\n");
    prt("\n");
    $in_lat = $lat;
    $in_lon = $lon;
    $in_zoom = $z;
    show_values();
}

sub show_hgt_ok($$) {
    my ($inlat,$inlon) = @_;
    my $lat = $inlat;
    my $lon = $inlon;
    my ($ilat,$ilon);
    my ($nlat,$nlon);
    #  Its south western corner can be deduced from its file name: 
    # for example, n51e002.hgt covers the area between N 51 E 2 and N 52 E 3, 
    # and s14w077.hgt covers S 14 W 77 to S 13 W 76 -14,-77 to -13,-76
    my $ns = 'N';
    my $ew = 'E';
    if ($lat < 0.0) {
        $ns = 'S';
        $lat *= -1;
        $ilat = int($lat);
        if ($lat - $ilat) {
            $lat += 1;
        }
    }
    if ($lon < 0.0) {
        $ew = 'W';
        $lon *= -1;
        $ilon = int($lon);
        if ($lon - $ilon) {
            $lon += 1;
        }
    }
    $ilon = int($lon);
    $ilat = int($lat);
    my $file = $ns.sprintf("%02d",$ilat).$ew.sprintf("%03d",$ilon).".hgt";
    if ($inlat < 0.0) {
        $nlat = $ilat - 1;
    } else {
        $nlat = $ilat + 1;
    }
    if ($inlon < 0.0) {
        $nlon = $ilon - 1;
    } else {
        $nlon = $ilon + 1;
    }
    prt("lat,lon $inlat,$inlon, ilat,ilon $ilat,$ilon, file $file to \n");
    prt("Range BL $ilat $ns $ilon $ew to TR $nlat $ns $nlon $ew\n");
}

sub do_test4() {
    # center
    my $inlat = 51.562;
    my $inlon = -1.747;
    show_hgt($inlat,$inlon);

    $inlon = -1;
    show_hgt($inlat,$inlon);
    $inlat = -14;
    $inlon = -77;
    show_hgt($inlat,$inlon);
    prt("s14w077.hgt covers S 14 W 77 to S 13 W 76 -14,-77 to -13,-76\n");

}

sub do_test5() {
    my ($as1,$xdim,$dist,$circ,$rad,$degs,$cols);
    $rad = 6378100;
    $circ = 2 * $rad * pi;
    $as1 = 30.87; # meters at the equator
    $as1 = $circ / (3500*360);
    my $HGT1_XDIM = 0.000277778;
    prt("from : http://www.webgis.com/srtm30.html SRTM30 - Shuttle Radar Topography Mission\n");
    prt("Global Coverage (~900m) 30 arc-sec\n");
    prt("The radius of the earth at the equator is 6378100 meters (approx.)\n");
    prt("The circumference is 40074784 (2*radius*pi) meters ($circ)\n");
    prt("The number of meters to an ark second is (40074784/(3600*360)) or 30.922 m ($as1)\n"); 
    prt("1 Nautical Mile = 1852 Meters\n");
    $cols = 10800;
    $degs = 90;
    $xdim = ($degs / $cols);   # GTOPO_XDEM
    $dist = (int(($as1 * $xdim * $cols) * 10) / 10);
    prt("For 90 degrees in 10800 columns xdim = $xdim, $dist meters\n"); 
    #define SRTM30_COLS 4800    // each 40 degrees
    #define SRTM30_ROWS 6000    // each 50 degrees
    #define SRTM30_COL2 7200    // each 60 degrees
    #define SRTM30_ROW2 3600    // each 30 degrees
    $degs = 40;
    $cols = 4800;
    $xdim = ($degs / $cols);
    #$dist = $as1 * $xdim * $cols;
    $dist = (int(($as1 * $xdim * $cols) * 10) / 10);
    prt("For 40 degrees in  4800 columns xdim = $xdim, $dist meters\n"); 
    $degs = 1;
    $cols = 1201;
    $xdim = ($degs / $cols);
    $dist = (int(($as1 * $xdim * $cols) * 10) / 10);
    prt("For  1 degrees in  1201 columns xdim = $xdim, $dist meters\n"); 
    $degs = 1;
    $cols = 3601;
    $xdim = ($degs / $cols);
    $dist = (int(($as1 * $xdim * $cols) * 10) / 10);
    prt("For  1 degrees in  3601 columns xdim = $xdim, $dist meters ($HGT1_XDIM)\n"); 

    pgm_exit(1,"");
}


sub do_test6() {
    my $cnt = scalar @srtm30_tile_table;
    my $lat = 27.988;
    my $lon = 86.9253;
    my $tile = get_srtm30_tile($lat,$lon);
    prt("Got $cnt tile sets. lat,lon $lat,$lon is tile $tile\n");
    pgm_exit(1,"");
}

sub show_path($) {
    my $path = shift;
    my @arr = split(/\//,$path);
    my $cnt = scalar @arr;
    my ($i);
    if ($cnt != 4) {
        prt("Expected split to yield 4 items. Got $cnt\n");
        for ($i = 0; $i < $cnt; $i++) {
            prt("$i [".$arr[$i]."]\n");
        }
    } else {
        $in_zoom = $arr[1];
        my $xlon = $arr[2];
        my $ylat = $arr[3];
        $ylat =~ s/\..*$//;
        ($in_lon,$in_lat) = getLonLat($xlon,$ylat,$in_zoom);
        $in_lat -= 0.00005; # hmmm, need to REDUCE lat
        show_values();
        my ($z,$mz);
        if ($in_zoom > 4) {
            $z = $in_zoom - 3;
        } else {
            $z = 0;
        }
        $mz = $z + 6;
        prt("\nUsing lat,lon ".get_double($in_lat).",".get_double($in_lon)." for zooms $z to $mz\nPaths ");
        for (; $z <= $mz; $z++) {
            my ($x, $y) = getTileNumber($in_lat,$in_lon,$z);
            my $dir = "/$z/$x/$y.png";
            prt("$dir ");
        }
        prt("\n");
    }
}

#########################################
### MAIN ###
#do_test();
#do_test2();
#do_test3();
#do_test4();
#do_test5();
#do_test6();
parse_args(@ARGV);

if (length($slippy_dir)) {
    show_path($slippy_dir);
} else {
    if ($in_zoom =~ /:/) {
        my @arr = split(":",$in_zoom);
        my ($z);
        for ($z = $arr[0]; $z <= $arr[1]; $z++) {
            $in_zoom = $z;
            prt("\nProcessing zoom $z\n");
            show_values();
        }
    } else {
        show_values();
    }
}
pgm_exit(0,"\n");
########################################

sub need_arg {
    my ($arg,@av) = @_;
    pgm_exit(1,"ERROR: [$arg] must have a following argument!\n") if (!@av);
}

sub parse_args {
    my (@av) = @_;
    my ($arg,$sarg);
    my $verb = VERB2();
    my $cnt = 0;
    while (@av) {
        $arg = $av[0];
        if ( ($arg =~ /^-/) && !($arg =~ /^-\d+/) ) {
        #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)");
            } elsif ($sarg =~ /^v/) {
                if ($sarg =~ /^v.*(\d+)$/) {
                    $verbosity = $1;
                } else {
                    while ($sarg =~ /^v/) {
                        $verbosity++;
                        $sarg = substr($sarg,1);
                    }
                }
                $verb = VERB2();
                prt("Verbosity = $verbosity\n") if ($verb);
            } elsif ($sarg =~ /^l/) {
                if ($sarg =~ /^ll/) {
                    $load_log = 2;
                } else {
                    $load_log = 1;
                }
                prt("Set to load log at end. ($load_log)\n") if ($verb);
            } elsif ($sarg =~ /^o/) {
                need_arg(@av);
                shift @av;
                $sarg = $av[0];
                $out_file = $sarg;
                prt("Set out file to [$out_file].\n") if ($verb);
            } elsif ($sarg =~ /^b/) {
                $list_touching_tiles = 1;
                prt("Set to show 8 touching tiles\n") if ($verb);
            } else {
                pgm_exit(1,"ERROR: Invalid argument [$arg]! Try -?\n");
            }
        } elsif ($arg =~ /^\//) {
            $slippy_dir = $arg;
        } else {
            if ($cnt == 0) {
                $in_lat = $arg;
                prt("Set input lat to [$in_lat]\n") if ($verb);
            } elsif ($cnt == 1) {
                $in_lon = $arg;
                prt("Set input lon to [$in_lon]\n") if ($verb);
            } elsif ($cnt == 2) {
                $in_zoom = $arg;
                prt("Set input zoom to [$in_zoom]\n") if ($verb);
            } else {
                pgm_exit(1,"What is this $arg! Already have lat,lon,zoom $in_lat,$in_lon,$in_zoom\n");
            }
            $cnt++;
        }
        shift @av;
    }

    if (($cnt != 3) && (length($slippy_dir) == 0)) {
        prt("ERROR: No slippy directory, must begin with '/', and\n");
        pgm_exit(1,"       No 'lat lon zoom' found in command!\n");
    }

}

sub give_help {
    prt("$pgmname: version $VERS\n");
    prt("Usage: $pgmname [options] lat lon zoom, OR\n");
    prt("Usage: $pgmname [options] /7/23/17.png (ie a slippy path - zoom/xlon/ylat.png)\n");
    prt("Options:\n");
    prt(" --help  (-h or -?) = This help, and exit 0.\n");
    prt(" --verb[n]     (-v) = Bump [or set] verbosity. def=$verbosity\n");
    prt(" --load        (-l) = Load LOG at end. ($outfile)\n");
    prt(" --out <file>  (-o) = Write output to this file.\n");
    prt(" --block       (-b) = Show paths to block of 9 tiles. (def=$list_touching_tiles)\n");
    prt("\n");
    prt("  Given lat lon zoom values, show the bounding box of that tile.\n");
    prt("  A zoom range can be given by 1:5. Will process zooms 1 to 5 inclusive.\n");

}

# eof - slippydirs.pl