libgeo-inverse-perl 0.05-1 / usr / share / perl5 / Geo / Inverse.pm

package Geo::Inverse;

=head1 NAME

Geo::Inverse - Calculate geographic distance from a lat & lon pair.

=head1 SYNOPSIS

  use Geo::Inverse;
  my $obj = Geo::Inverse->new(); # default "WGS84"
  my ($lat1,$lon1,$lat2,$lon2)=(38.87, -77.05, 38.95, -77.23);
  my ($faz, $baz, $dist)=$obj->inverse($lat1,$lon1,$lat2,$lon2); #array context
  my $dist=$obj->inverse($lat1,$lon1,$lat2,$lon2);              #scalar context
  print "Input Lat: $lat1  Lon: $lon1\n";
  print "Input Lat: $lat2 Lon: $lon2\n";
  print "Output Distance: $dist\n";
  print "Output Forward Azimuth: $faz\n";
  print "Output Back Azimuth: $baz\n";

=head1 DESCRIPTION

This module is a pure Perl port of the NGS program in the public domain "inverse" by Robert (Sid) Safford and Stephen J. Frakes.  


=cut

use strict;
use vars qw($VERSION);
use Geo::Constants qw{PI};
use Geo::Functions qw{rad_deg deg_rad};

$VERSION = sprintf("%d.%02d", q{Revision: 0.05} =~ /(\d+)\.(\d+)/);

=head1 CONSTRUCTOR

=head2 new

The new() constructor may be called with any parameter that is appropriate to the ellipsoid method which establishes the ellipsoid.

  my $obj = Geo::Inverse->new(); # default "WGS84"

=cut

sub new {
  my $this = shift();
  my $class = ref($this) || $this;
  my $self = {};
  bless $self, $class;
  $self->initialize(@_);
  return $self;
}

=head1 METHODS

=cut

sub initialize {
  my $self = shift();
  my $param = shift()||undef();
  $self->ellipsoid($param);
}

=head2 ellipsoid

Method to set or retrieve the current ellipsoid object.  The ellipsoid is a Geo::Ellipsoids object.

  my $ellipsoid=$obj->ellipsoid;  #Default is WGS84

  $obj->ellipsoid('Clarke 1866'); #Built in ellipsoids from Geo::Ellipsoids
  $obj->ellipsoid({a=>1});        #Custom Sphere 1 unit radius

=cut

sub ellipsoid {
  my $self = shift();
  if (@_) {
    my $param=shift();
    use Geo::Ellipsoids;
    my $obj=Geo::Ellipsoids->new($param);
    $self->{'ellipsoid'}=$obj;
  }
  return $self->{'ellipsoid'};
}

=head2 inverse

This method is the user frontend to the mathematics. This interface will not change in future versions.

  my ($faz, $baz, $dist)=$obj->inverse($lat1,$lon1,$lat2,$lon2);

=cut

sub inverse {
  my $self=shift();
  my $lat1=shift();      #degrees
  my $lon1=shift();      #degrees
  my $lat2=shift();      #degrees
  my $lon2=shift();      #degrees
  my ($faz, $baz, $dist)=$self->_inverse(rad_deg($lat1), rad_deg($lon1),
                                         rad_deg($lat2), rad_deg($lon2));
  return wantarray ? (deg_rad($faz), deg_rad($baz), $dist) : $dist;
}

########################################################################
#
#   This function was copied from Geo::Ellipsoid
#   Copyright 2005-2006 Jim Gibson, all rights reserved.
#   
#   This program is free software; you can redistribute it and/or modify it
#   under the same terms as Perl itself.
#
#      internal functions
#
#	inverse
#
#	Calculate the displacement from origin to destination.
#	The input to this subroutine is 
#	  ( latitude-1, longitude-1, latitude-2, longitude-2 ) in radians.
#
#	Return the results as the list (range,bearing) with range in meters
#	and bearing in radians.
#
########################################################################

sub _inverse {
  my $self = shift;
  my( $lat1, $lon1, $lat2, $lon2 ) = (@_);

  my $ellipsoid=$self->ellipsoid;
  my $a = $ellipsoid->a;
  my $f = $ellipsoid->f;

  my $eps = 1.0e-23;
  my $max_loop_count = 20;
  my $pi=PI;
  my $twopi = 2 * $pi;

  my $r = 1.0 - $f;
  my $tu1 = $r * sin($lat1) / cos($lat1);
  my $tu2 = $r * sin($lat2) / cos($lat2);
  my $cu1 = 1.0 / ( sqrt(($tu1*$tu1) + 1.0) );
  my $su1 = $cu1 * $tu1;
  my $cu2 = 1.0 / ( sqrt( ($tu2*$tu2) + 1.0 ));
  my $s = $cu1 * $cu2;
  my $baz = $s * $tu2;
  my $faz = $baz * $tu1;
  my $dlon = $lon2 - $lon1;

  my $x = $dlon;
  my $cnt = 0;
  my( $c2a, $c, $cx, $cy, $cz, $d, $del, $e, $sx, $sy, $y );
  do {
    $sx = sin($x);
    $cx = cos($x);
    $tu1 = $cu2*$sx;
    $tu2 = $baz - ($su1*$cu2*$cx);
    $sy = sqrt( $tu1*$tu1 + $tu2*$tu2 );
    $cy = $s*$cx + $faz;
    $y = atan2($sy,$cy);
    my $sa;
    if( $sy == 0.0 ) {
      $sa = 1.0;
    }else{
      $sa = ($s*$sx) / $sy;
    }

    $c2a = 1.0 - ($sa*$sa);
    $cz = $faz + $faz;
    if( $c2a > 0.0 ) {
      $cz = ((-$cz)/$c2a) + $cy;
    }
    $e = ( 2.0 * $cz * $cz ) - 1.0;
    $c = ( ((( (-3.0 * $c2a) + 4.0)*$f) + 4.0) * $c2a * $f )/16.0;
    $d = $x;
    $x = ( ($e * $cy * $c + $cz) * $sy * $c + $y) * $sa;
    $x = ( 1.0 - $c ) * $x * $f + $dlon;
    $del = $d - $x;
 
  } while( (abs($del) > $eps) && ( ++$cnt <= $max_loop_count ) );

  $faz = atan2($tu1,$tu2);
  $baz = atan2($cu1*$sx,($baz*$cx - $su1*$cu2)) + $pi;
  $x = sqrt( ((1.0/($r*$r)) -1.0 ) * $c2a+1.0 ) + 1.0;
  $x = ($x-2.0)/$x;
  $c = 1.0 - $x;
  $c = (($x*$x)/4.0 + 1.0)/$c;
  $d = ((0.375*$x*$x) - 1.0)*$x;
  $x = $e*$cy;

  $s = 1.0 - $e - $e;
  $s = (((((((( $sy * $sy * 4.0 ) - 3.0) * $s * $cz * $d/6.0) - $x) * 
    $d /4.0) + $cz) * $sy * $d) + $y ) * $c * $a * $r;

  # adjust azimuth to (0,360)
  $faz += $twopi if $faz < 0;

  return($faz, $baz, $s);
}

1;

__END__

=head1 TODO

Add more tests.

=head1 BUGS

Please send to the geo-perl email list.

=head1 LIMITS

No guarantees that Perl handles all of the double precision calculations in the same manner as Fortran.

=head1 AUTHOR

Michael R. Davis qw/perl michaelrdavis com/

=head1 LICENSE

Copyright (c) 2006 Michael R. Davis (mrdvt92)

This library is free software; you can redistribute it and/or modify
it under the same terms as Perl itself.

=head1 SEE ALSO

Net::GPSD
Geo::Ellipsoid
GIS::Distance::GeoEllipsoid