#!/usr/bin/perl -CSDAL

use strict;
use warnings;
use utf8;

use POSIX qw(tan atan log exp M_PI M_PI_2 M_PI_4);

if( @ARGV != 3 && @ARGV != 1 || $ARGV[0] !~ /^\d+$/ ) {
    print <<EOF;
usage: osmtilecoord.pl <zoom> [ <x> <y> ]
Converts longitude/latitude to and from OpenStreetmap tile index.  Coordinate
arguments are taken to be tile indices if both are integers without decimal
point.  Without coordinate arguments, reads from stdin.
EOF
    exit;
}

my $cmdargs= @ARGV == 3;
my ($z, $x, $y)= splice @ARGV, 0, 3;
my $scale= 1 << $z;
my $line;

while( $cmdargs || defined($line= <>) ) {
    unless( $cmdargs ) {
        chomp $line;
        $line =~ s/^\s+//;
        ($x, $y)= split /\s+/, $line, 2;
    }
    if( $x =~ /\./ || $y =~ /\./ ) {    # lon/lat -> tile ind
        die "Longitude / latitude coordinates have to be in [-180, 180] x [-85, 85]\n"
            unless abs($x) <= 180.0 && abs($y) <= 85.0;
        $x= int( (($x + 180) * $scale) / 360.0 );
        $x= $scale - 1 if $x >= $scale;
        $y= int( ((M_PI - log(tan(0.5 * ($y + 90.0)/180.0*M_PI))) * $scale) / (2*M_PI) );
# see https://en.wikipedia.org/wiki/Web_Mercator and icbm.c (links)
    }
    else {          # tile ind -> lon/lat (centre of tile)
        $x= ($x * 360.0 + 180.0) / $scale - 180.0;
        $y= (2 * atan(exp(M_PI - ($y * 2*M_PI + M_PI) / $scale)) * 180.0 / M_PI) - 90.0;
    }
    print "$x  $y\n";
    last if $cmdargs;
}