# GeoCoord.pm
# Geo coordinates class
# 30/09/2011 geoff mclane http://geoffair.net/mperl
package GeoCoord;

use strict;
use warnings;
use Carp;
use Math::Trig;
#use constant PI => 4 * atan2(1, 1);

require 'fg_wsg84.pl' or die "ERROR: Unable to load 'fg_wsg84.pl'!\n";

#my $PI = PI;
#my $PI = 4 * atan2(1, 1);
my $PI = 3.14159265358979323846264338327950288419716939937510;
my $D2R = $PI / 180;
my $R2D = 180 / $PI;
my $EPSILON = 2.22507e-308;
my $ERAD = 6378137.0;

our ($VERSION);
$VERSION = "0.1";

sub new {
    my $class = shift;
    my $self = {};

    $self->{'LAT'} = undef;   # degrees
    $self->{'LON'} = undef;   # degrees
    $self->{'ALT'} = undef;   # feet
    $self->{'X'} = undef;
    $self->{'Y'} = undef;
    $self->{'Z'} = undef;
    $self->{'PD'} = 1;    # polar - deg,deg,feet
    $self->{'CD'} = 1;    # cartesian

    bless ($self, $class);

    return $self;
}

sub lon {
    my $self = shift;
    if (@_) { 
        $self->{'LON'} = shift;
        $self->{'CD'} = 1;
        #print "Set LON = ".$self->{'LON'}." ".$self->{'CD'}."\n";
    } else {
        pupdate($self);
    }
    return $self->{'LON'};
}

sub lat {
    my $self = shift;
    if (@_) { 
        $self->{'LAT'} = shift;
        $self->{'CD'} = 1;
        #print "Set LAT = ".$self->{'LAT'}." ".$self->{'CD'}."\n";
    } else {
        pupdate($self);
    }
    return $self->{'LAT'};
}

sub alt {
    my $self = shift;
    if (@_) { 
        $self->{'ALT'} = shift;
        $self->{'CD'} = 1;
        #print "Set ALT = ".$self->{'ALT'}." ".$self->{'CD'}."\n";
    } else {
        pupdate($self);
    }
    return $self->{'ALT'};
}

sub cupdate {
    my $self = shift;
    if ( $self->{'CD'} ) {
        my @geod = ( $self->{'LAT'}, $self->{'LON'}, $self->{'ALT'} );
        my @cart = ( 0, 0, 0 );
        myGeod_DEG_FT_ToCart(\@geod,\@cart);
        $self->{'X'} = $cart[0];
        $self->{'Y'} = $cart[1];
        $self->{'Z'} = $cart[2];
        $self->{'CD'} = 0;
    #} else {
    #    print "NOT CD ". $self->{'CD'} . "\n";
    }
}

sub pupdate {
    my $self = shift;
    if ( $self->{'PD'} ) {
        my @geod = ( 0, 0, 0 );
        my @cart = ( $self->{'X'}, $self->{'Y'}, $self->{'Z'} );
        myCart_To_Geod(\@cart,\@geod);
        $self->{'LAT'} = $geod[0];
        $self->{'LON'} = $geod[1];
        $self->{'ALT'} = $geod[2];
        $self->{'PD'} = 0;
        #print "\@cart = ( ".$self->{'X'}.", ".$self->{'Y'}.", ".$self->{'Z'}."\n";
        #print "\@geod = ( ".$self->{'LAT'}.", ".$self->{'LON'}.", ".$self->{'ALT'}."\n";
    #} else {
        #print "NOT PD ". $self->{'PD'} . "\n";
    }
}


sub get_x {
    my $self = shift;
    $self->cupdate();
    return $self->{'X'};
}

sub get_y {
    my $self = shift;
    $self->cupdate();
    return $self->{'Y'};
}

sub get_z {
    my $self = shift;
    $self->cupdate();
    return $self->{'Z'};
}

sub set_x {
    my $self = shift;
    my $val = shift;
    $self->{'X'} = $val;
    $self->{'PD'} = 1;
    return $self->{'X'};
}

sub set_y {
    my $self = shift;
    my $val = shift;
    $self->{'Y'} = $val;
    $self->{'PD'} = 1;
    return $self->{'Y'};
}

sub set_z {
    my $self = shift;
    my $val = shift;
    $self->{'Z'} = $val;
    $self->{'PD'} = 1;
    return $self->{'Z'};
}

sub apply_course_distance {
    my $self = shift;
    my ($course,$dist) = @_; 
    my $lat = $self->{'LAT'};
    my $lon = $self->{'LON'};
    my ($nlat,$nlon,$naz);
    fg_geo_direct_wgs_84( $lat, $lon, $course, $dist, \$nlat, \$nlon, \$naz );
    # print "fg_geo_direct_wgs_84( $lat, $lon, $course, $dist, $nlat, $nlon, $naz );\n";
    $self->{'LAT'} = $nlat;
    $self->{'LON'} = $nlon;
    $self->{'CD'} = 1;
    return $self;
}

# this does not seem to give the right results???
sub apply_course_distance2 {
    my $self = shift;
    my ($course,$dist) = @_; 
#	apply_course_distance: func(course, dist) {
#		me._pupdate();
    $self->pupdate();
#		course *= D2R;
    $course *= $D2R;
#		dist /= ERAD;
    $dist /= $ERAD;
    my $latR = $self->{'LAT'} * $D2R;
#		me._lat = math.asin(math.sin(me._lat) * math.cos(dist)
#				+ math.cos(me._lat) * math.sin(dist) * math.cos(course));
    my $lat = asin(sin($latR) * cos($dist) +
        cos($latR) * sin($dist) * cos($course));
#		if (math.cos(me._lat) > EPSILON)
    my $lon = $self->{'LON'} * $D2R;
    if (cos($lat) > $EPSILON) {
#			me._lon = math.pi - math.mod(math.pi - me._lon
#					- math.asin(math.sin(course) * math.sin(dist)
#					/ math.cos(me._lat)), 2 * math.pi);
#        $lon = $PI - mod($PI - $lon 
#					- asin(sin($course) * sin($dist)
#					/ cos($lat)), 2 * $PI);
        #$lon = $PI - ($PI - $lon 
		#			- asin(sin($course) * sin($dist)
		#			/ cos($lat))) & (2 * $PI);
        my $term1 = ($PI - $lon - asin(sin($course) * sin($dist) / cos($lat)));
        my $term2 = (2 * $PI);
        my $mod = (($term1 / $term2) - int($term1 / $term2));
        $lon = $PI - $mod; 
    }
#
#		me._cdirty = 1;
#		me;
    $self->{'LAT'} = $lat * $R2D;
    $self->{'LON'} = $lon * $R2D;
    $self->{'CD'} = 1;
    return $self;
}

# normalize degree to 0 <= angle < 360
#
sub normdeg {
    my $self = shift;
    my $angle = shift;
	while ($angle < 0) {
		$angle += 360;
    }
	while ($angle >= 360) {
		$angle -= 360;
    }
	return $angle;
}

sub bucket_span {
    my $self = shift;
    my ($mlat) = @_;
    my $span = 360;
    if ($mlat >= 89.9) {
		$span = 360.0;
	} elsif ($mlat >= 88.0) {
		$span = 8.0;
    } elsif ($mlat >= 86.0) {
        $span = 4.0;
    } elsif ($mlat >= 83.0) {
		$span = 2.0;
    } elsif ($mlat >= 76.0) {
        $span = 1.0;
    } elsif ($mlat >= 62.0) {
		$span = 0.5;
    } elsif ($mlat >= 22.0) {
        $span = 0.25;
    } elsif ($mlat >= -22.0) {
		$span = 0.125;
    } elsif ($mlat >= -62.0) {
		$span = 0.25;
    } elsif ($mlat >= -76.0) {
        $span = 0.5;
    } elsif ($mlat >= -83.0) {
		$span = 1.0;
    } elsif ($mlat >= -86.0) {
        $span = 2.0;
    } elsif ($mlat >= -88.0) {
		$span = 4.0;
    } elsif ($mlat >= -89.0) {
		$span = 8.0;
	} else {
		$span = 360.0;
    }
    return $span;
}


sub tile_index {
    my $self = shift;
    my ($lat, $lon) = @_;
    ### print "Tile index for [$lat] [$lon]\n";
	### my $lat_floor = int($lat); # floor($lat);
	### my $lon_floor = int($lon); # floor($lon);
	my $lat_floor = ($lat < 0.0) ? -int(-$lat) - 1 : int($lat); # floor($lat);
	my $lon_floor = ($lon < 0.0) ? -int(-$lat) - 1 : int($lon); # floor($lon);
	my $span = $self->bucket_span($lat);
	my $xx = 0;
	if ($span < 0.0000001) {
		$lon = 0;
	} elsif ($span <= 1.0) {
		$xx = int(($lon - $lon_floor) / $span);
	} else {
		if ($lon >= 0) {
			$lon = int(int($lon / $span) * $span);
		} else {
			$lon = int(int(($lon + 1) / $span) * $span - $span);
			if ($lon < -180) {
				$lon = -180;
            }
		}
	}
	my $yy = int(($lat - $lat_floor) * 8);
    return ($lon_floor + 180) * 16384 + ($lat_floor + 90) * 64 + $yy * 8 + $xx;
}

1; # for require or use to succeed
# eof - GeoCoord.pm
