#!/usr/bin/perl -w
# NAME: maptest.pl
# AIM: Convert a lat long zoom to a tile number
# 17/02/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; # qw(great_circle_distance great_circle_direction deg2rad rad2deg);
##use POSIX;
use POSIX qw/floor/;

my $in_lat = 45;
my $in_lon = 123;
my $in_zoom = 9;

sub prt($) { print shift; }

sub mlog_10 {
  my $n = shift;
  return log($n)/log(10);
}


# from Theo HTML
#// The math
#// http://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Tile_servers
#	function lon2tile( lon, zoom ) {
#		return Math.floor( ( lon + 180 ) / 360 * pow( 2, zoom ) );
#	}
#	function lat2tile( lat, zoom ) {
#		return Math.floor(( 1 - Math.log( Math.tan( lat * pi / 180) + 1 / cos( lat * pi / 180)) / pi )/2 * pow(2, zoom) );
#	}
#	function tile2lon( x, z ) {
#		return ( x / pow( 2, z ) * 360 - 180 );
#	}
#	function tile2lat( y, z ) {
#		var n = pi - 2 * pi * y / pow( 2, z );
#		return 180 / pi * Math.atan( 0.5 * ( Math.exp( n ) - Math.exp( -n ) ));
#	}
sub lat2tilef($$) {
    my ($lat,$zoom) = @_;
    my $lat_rad = deg2rad($lat);
    my $n = 2 ** $zoom;
    my $ytile = floor((1.0 - mlog_10(tan($lat_rad) + (1 / cos($lat_rad))) / pi) / 2.0 * $n);
    return $ytile;
}

sub lon2tilef($$) {
    my ($lon,$zoom) = @_;
    my $n = 2 ** $zoom;
    my $xtile = floor($n * (($lon + 180) / 360));
    return $xtile;
}


# Python
# Lon./lat. to tile numbers
# import math
# def deg2num(lat_deg, lon_deg, zoom):
#  lat_rad = math.radians(lat_deg)
#  n = 2.0 ** zoom
#  xtile = int((lon_deg + 180.0) / 360.0 * n)
#  ytile = int((1.0 - math.log(math.tan(lat_rad) + (1 / math.cos(lat_rad))) / math.pi) / 2.0 * n)
#  return (xtile, ytile)
sub lat2tile($$) {
    my ($lat,$zoom) = @_;
    my $lat_rad = deg2rad($lat);
    my $n = 2 ** $zoom;
    my $ytile = int((1.0 - mlog_10(tan($lat_rad) + (1 / cos($lat_rad))) / pi) / 2.0 * $n);
    return $ytile;
}

sub lon2tile($$) {
    my ($lon,$zoom) = @_;
    my $n = 2 ** $zoom;
    my $xtile = int($n * (($lon + 180) / 360));
    return $xtile;
}


sub show_mappy() {
    prt("Converting lat $in_lat, lon $in_lon, zoom = $in_zoom to mappy\n");
    my $lat_rad = deg2rad($in_lat);
    my $n = 2 ** $in_zoom;
    my $xtile = int($n * (($in_lon + 180) / 360));
    my $ytile = int((1.0 - mlog_10(tan($lat_rad) + (1 / cos($lat_rad))) / pi) / 2.0 * $n);
    my $x = lon2tile($in_lon,$in_zoom);
    my $y = lat2tile($in_lat,$in_zoom);
    prt("xtile=$xtile, ytile=$ytile ($x/$y)\n");
}

sub test() {
    my $z = 19;
    my ($x,$y,$mx,$my,$tx,$ty,$cx,$cy);
    prt("Running test...\n");
    #for ($x = 0; $x < 180; $x += 0.1) {
    #    for ($y = 0; $y < 90; $y += 0.1) {
    #for ($x = -179; $x < 180; $x += 0.1) {
    #    for ($y = -89; $y < 90; $y += 0.1) {
    # from :
    for ($x = -179; $x < 180; $x += 0.1) {
        for ($y = -89; $y < 90; $y += 0.1) {
            $mx = lon2tile($x,$z);
            $my = lat2tile($y,$z);
            $tx = lon2tilef($x,$z);
            $ty = lat2tilef($y,$z);
            #if (($mx != $tx)||($my != $ty)) {
            #    prt("lon/lat/zoom $x/$y/$z mx=$mx my=$my tx=$tx ty=$ty\n");
            #}
            $cx = sprintf("%.2f",$x);
            $cy = sprintf("%.2f",$y);

            if ($mx != $tx) {
                prt("lon/lat/zoom $cx/$cy/$z DIFF mx/tx=$mx/$tx\n");
            }
            if ($my != $ty) {
                prt("lon/lat/zoom $cx/$cy/$z DIFF my/ty=$my/$ty\n");
            }
            ###prt("lon/lat/zoom $x/$y/$z mx=$mx my=$my\n");
        }
        # prt("Next longitude...\n");
    }
    exit(0);

}

show_mappy();
test();

# eof - maptest.pl