| 0,0 → 1,950 |
|---|
| # Geo::Ellipsoid |
| # |
| # This package implements an Ellipsoid class to perform latitude |
| # and longitude calculations on the surface of an ellipsoid. |
| # |
| # This is a Perl conversion of existing Fortran code (see |
| # ACKNOWLEDGEMENTS) and the author of this class makes no |
| # claims of originality. Nor can he even vouch for the |
| # results of the calculations, although they do seem to |
| # work for him and have been tested against other methods. |
| |
| package Geo::Ellipsoid; |
| |
| use warnings; |
| use strict; |
| use 5.006_00; |
| |
| use Scalar::Util 'looks_like_number'; |
| use Math::Trig; |
| use Carp; |
| |
| =head1 NAME |
| |
| Geo::Ellipsoid - Longitude and latitude calculations using an ellipsoid model. |
| |
| =head1 VERSION |
| |
| Version 1.12, released July 4, 2008. |
| |
| =cut |
| |
| our $VERSION = '1.12'; |
| our $DEBUG = 0; |
| |
| =head1 SYNOPSIS |
| |
| use Geo::Ellipsoid; |
| $geo = Geo::Ellipsoid->new(ellipsoid=>'NAD27', units=>'degrees'); |
| @origin = ( 37.619002, -122.374843 ); # SFO |
| @dest = ( 33.942536, -118.408074 ); # LAX |
| ( $range, $bearing ) = $geo->to( @origin, @dest ); |
| ($lat,$lon) = $geo->at( @origin, 2000, 45.0 ); |
| ( $x, $y ) = $geo->displacement( @origin, $lat, $lon ); |
| @pos = $geo->location( $lat, $lon, $x, $y ); |
| |
| =head1 DESCRIPTION |
| |
| Geo::Ellipsoid performs geometrical calculations on the surface of |
| an ellipsoid. An ellipsoid is a three-dimension object formed from |
| the rotation of an ellipse about one of its axes. The approximate |
| shape of the earth is an ellipsoid, so Geo::Ellipsoid can accurately |
| calculate distance and bearing between two widely-separated locations |
| on the earth's surface. |
| |
| The shape of an ellipsoid is defined by the lengths of its |
| semi-major and semi-minor axes. The shape may also be specifed by |
| the flattening ratio C<f> as: |
| |
| f = ( semi-major - semi-minor ) / semi-major |
| |
| which, since f is a small number, is normally given as the reciprocal |
| of the flattening C<1/f>. |
| |
| The shape of the earth has been surveyed and estimated differently |
| at different times over the years. The two most common sets of values |
| used to describe the size and shape of the earth in the United States |
| are 'NAD27', dating from 1927, and 'WGS84', from 1984. United States |
| Geological Survey topographical maps, for example, use one or the |
| other of these values, and commonly-available Global Positioning |
| System (GPS) units can be set to use one or the other. |
| See L<"DEFINED ELLIPSOIDS"> below for the ellipsoid survey values |
| that may be selected for use by Geo::Ellipsoid. |
| |
| =cut |
| |
| # class data and constants |
| our $degrees_per_radian = 180/pi; |
| our $eps = 1.0e-23; |
| our $max_loop_count = 20; |
| our $twopi = 2 * pi; |
| our $halfpi = pi/2; |
| our %defaults = ( |
| ellipsoid => 'WGS84', |
| units => 'radians', |
| distance_units => 'meter', |
| longitude => 0, |
| latitude => 1, # allows use of _normalize_output |
| bearing => 0, |
| ); |
| our %distance = ( |
| 'foot' => 0.3048, |
| 'kilometer' => 1_000, |
| 'meter' => 1.0, |
| 'mile' => 1_609.344, |
| 'nm' => 1_852, |
| ); |
| |
| # set of ellipsoids that can be used. |
| # values are |
| # 1) a = semi-major (equatorial) radius of Ellipsoid |
| # 2) 1/f = reciprocal of flattening (f), the ratio of the semi-minor |
| # (polar) radius to the semi-major (equatorial) axis, or |
| # polar radius = equatorial radius * ( 1 - f ) |
| |
| our %ellipsoids = ( |
| 'AIRY' => [ 6377563.396, 299.3249646 ], |
| 'AIRY-MODIFIED' => [ 6377340.189, 299.3249646 ], |
| 'AUSTRALIAN' => [ 6378160.0, 298.25 ], |
| 'BESSEL-1841' => [ 6377397.155, 299.1528128 ], |
| 'CLARKE-1880' => [ 6378249.145, 293.465 ], |
| 'EVEREST-1830' => [ 6377276.345, 300.8017 ], |
| 'EVEREST-MODIFIED' => [ 6377304.063, 300.8017 ], |
| 'FISHER-1960' => [ 6378166.0, 298.3 ], |
| 'FISHER-1968' => [ 6378150.0, 298.3 ], |
| 'GRS80' => [ 6378137.0, 298.25722210088 ], |
| 'HOUGH-1956' => [ 6378270.0, 297.0 ], |
| 'HAYFORD' => [ 6378388.0, 297.0 ], |
| 'IAU76' => [ 6378140.0, 298.257 ], |
| 'KRASSOVSKY-1938' => [ 6378245.0, 298.3 ], |
| 'NAD27' => [ 6378206.4, 294.9786982138 ], |
| 'NWL-9D' => [ 6378145.0, 298.25 ], |
| 'SOUTHAMERICAN-1969' => [ 6378160.0, 298.25 ], |
| 'SOVIET-1985' => [ 6378136.0, 298.257 ], |
| 'WGS72' => [ 6378135.0, 298.26 ], |
| 'WGS84' => [ 6378137.0, 298.257223563 ], |
| ); |
| |
| =head1 CONSTRUCTOR |
| |
| =head2 new |
| |
| The new() constructor may be called with a hash list to set the value of the |
| ellipsoid to be used, the value of the units to be used for angles and |
| distances, and whether or not the output range of longitudes and bearing |
| angles should be symmetric around zero or always greater than zero. |
| The initial default constructor is equivalent to the following: |
| |
| my $geo = Geo::Ellipsoid->new( |
| ellipsoid => 'WGS84', |
| units => 'radians' , |
| distance_units => 'meter', |
| longitude => 0, |
| bearing => 0, |
| ); |
| |
| The constructor arguments may be of any case and, with the exception of |
| the ellipsoid value, abbreviated to their first three characters. |
| Thus, ( UNI => 'DEG', DIS => 'FEE', Lon => 1, ell => 'NAD27', bEA => 0 ) |
| is valid. |
| |
| =cut |
| |
| sub new |
| { |
| my( $class, %args ) = @_; |
| my $self = {%defaults}; |
| print "new: @_\n" if $DEBUG; |
| foreach my $key ( keys %args ) { |
| my $val = $args{$key}; |
| if( $key =~ /^ell/i ) { |
| $self->{ellipsoid} = uc $args{$key}; |
| }elsif( $key =~ /^uni/i ) { |
| $self->{units} = $args{$key}; |
| }elsif( $key =~ /^dis/i ) { |
| $self->{distance_units} = $args{$key}; |
| }elsif( $key =~ /^lon/i ) { |
| $self->{longitude} = $args{$key}; |
| }elsif( $key =~ /^bea/i ) { |
| $self->{bearing} = $args{$key}; |
| }else{ |
| carp("Unknown argument: $key => $args{$key}"); |
| } |
| } |
| set_units($self,$self->{units}); |
| set_ellipsoid($self,$self->{ellipsoid}); |
| set_distance_unit($self,$self->{distance_units}); |
| set_longitude_symmetric($self,$self->{longitude}); |
| set_bearing_symmetric($self,$self->{bearing}); |
| print |
| "Ellipsoid(units=>$self->{units},distance_units=>" . |
| "$self->{distance_units},ellipsoid=>$self->{ellipsoid}," . |
| "longitude=>$self->{longitude},bearing=>$self->{bearing})\n" if $DEBUG; |
| bless $self,$class; |
| return $self; |
| } |
| |
| =head1 METHODS |
| |
| =head2 set_units |
| |
| Set the angle units used by the Geo::Ellipsoid object. The units may |
| also be set in the constructor of the object. The allowable values are |
| 'degrees' or 'radians'. The default is 'radians'. The units value is |
| not case sensitive and may be abbreviated to 3 letters. The units of |
| angle apply to both input and output latitude, longitude, and bearing |
| values. |
| |
| $geo->set_units('degrees'); |
| |
| =cut |
| |
| sub set_units |
| { |
| my $self = shift; |
| my $units = shift; |
| if( $units =~ /deg/i ) { |
| $units = 'degrees'; |
| }elsif( $units =~ /rad/i ) { |
| $units = 'radians'; |
| }else{ |
| croak("Invalid units specifier '$units' - please use either " . |
| "degrees or radians (the default)") unless $units =~ /rad/i; |
| } |
| $self->{units} = $units; |
| } |
| |
| =head2 set_distance_unit |
| |
| Set the distance unit used by the Geo::Ellipsoid object. The unit of |
| distance may also be set in the constructor of the object. The recognized |
| values are 'meter', 'kilometer', 'mile', 'nm' (nautical mile), or 'foot'. |
| The default is 'meter'. The value is not case sensitive and may be |
| abbreviated to 3 letters. |
| |
| $geo->set_distance_unit('kilometer'); |
| |
| For any other unit of distance not recogized by this method, pass a |
| numerical argument representing the length of the distance unit in |
| meters. For example, to use units of furlongs, call |
| |
| $geo->set_distance_unit(201.168); |
| |
| The distance conversion factors used by this module are as follows: |
| |
| Unit Units per meter |
| -------- --------------- |
| foot 0.3048 |
| kilometer 1000.0 |
| mile 1609.344 |
| nm 1852.0 |
| |
| =cut |
| |
| sub set_distance_unit |
| { |
| my $self = shift; |
| my $unit = shift; |
| print "distance unit = $unit\n" if $DEBUG; |
| |
| my $conversion = 0; |
| |
| if( defined $unit ) { |
| |
| my( $key, $val ); |
| while( ($key,$val) = each %distance ) { |
| my $re = substr($key,0,3); |
| print "trying ($key,$re,$val)\n" if $DEBUG; |
| if( $unit =~ /^$re/i ) { |
| $self->{distance_units} = $unit; |
| $conversion = $val; |
| |
| # finish iterating to reset 'each' function call |
| while( each %distance ) {} |
| last; |
| } |
| } |
| |
| if( $conversion == 0 ) { |
| if( looks_like_number($unit) ) { |
| $conversion = $unit; |
| }else{ |
| carp("Unknown argument to set_distance_unit: $unit\nAssuming meters"); |
| $conversion = 1.0; |
| } |
| } |
| }else{ |
| carp("Missing or undefined argument to set_distance_unit: ". |
| "$unit\nAssuming meters"); |
| $conversion = 1.0; |
| } |
| $self->{conversion} = $conversion; |
| } |
| |
| =head2 set_ellipsoid |
| |
| Set the ellipsoid to be used by the Geo::Ellipsoid object. See |
| L<"DEFINED ELLIPSOIDS"> below for the allowable values. The value |
| may also be set by the constructor. The default value is 'WGS84'. |
| |
| $geo->set_ellipsoid('NAD27'); |
| |
| =cut |
| |
| sub set_ellipsoid |
| { |
| my $self = shift; |
| my $ellipsoid = uc shift || $defaults{ellipsoid}; |
| print " set ellipsoid to $ellipsoid\n" if $DEBUG; |
| unless( exists $ellipsoids{$ellipsoid} ) { |
| croak("Ellipsoid $ellipsoid does not exist - please use " . |
| "set_custom_ellipsoid to use an ellipsoid not in valid set"); |
| } |
| $self->{ellipsoid} = $ellipsoid; |
| my( $major, $recip ) = @{$ellipsoids{$ellipsoid}}; |
| $self->{equatorial} = $major; |
| if( $recip == 0 ) { |
| carp("Infinite flattening specified by ellipsoid -- assuming a sphere"); |
| $self->{polar} = $self->{equatorial}; |
| $self->{flattening} = 0; |
| $self->{eccentricity} = 0; |
| }else{ |
| $self->{flattening} = ( 1.0 / $ellipsoids{$ellipsoid}[1]); |
| $self->{polar} = $self->{equatorial} * ( 1.0 - $self->{flattening} ); |
| $self->{eccentricity} = sqrt( 2.0 * $self->{flattening} - |
| ( $self->{flattening} * $self->{flattening} ) ); |
| } |
| } |
| |
| =head2 set_custom_ellipsoid |
| |
| Sets the ellipsoid parameters to the specified ( major semiaxis and |
| reciprocal flattening. A zero value for the reciprocal flattening |
| will result in a sphere for the ellipsoid, and a warning message |
| will be issued. |
| |
| $geo->set_custom_ellipsoid( 'sphere', 6378137, 0 ); |
| |
| =cut |
| |
| sub set_custom_ellipsoid |
| { |
| my $self = shift; |
| my( $name, $major, $recip ) = @_; |
| $name = uc $name; |
| $recip = 0 unless defined $recip; |
| if( $major ) { |
| $ellipsoids{$name} = [ $major, $recip ]; |
| }else{ |
| croak("set_custom_ellipsoid called without semi-major radius parameter"); |
| } |
| set_ellipsoid($self,$name); |
| } |
| |
| =head2 set_longitude_symmetric |
| |
| If called with no argument or a true argument, sets the range of output |
| values for longitude to be in the range [-pi,+pi) radians. If called with |
| a false or undefined argument, sets the output angle range to be |
| [0,2*pi) radians. |
| |
| $geo->set_longitude_symmetric(1); |
| |
| =cut |
| |
| sub set_longitude_symmetric |
| { |
| my( $self, $sym ) = @_; |
| # see if argument passed |
| if( $#_ > 0 ) { |
| # yes -- use value passed |
| $self->{longitude} = $sym; |
| }else{ |
| # no -- set to true |
| $self->{longitude} = 1; |
| } |
| } |
| |
| =head2 set_bearing_symmetric |
| |
| If called with no argument or a true argument, sets the range of output |
| values for bearing to be in the range [-pi,+pi) radians. If called with |
| a false or undefined argument, sets the output angle range to be |
| [0,2*pi) radians. |
| |
| $geo->set_bearing_symmetric(1); |
| |
| =cut |
| |
| sub set_bearing_symmetric |
| { |
| my( $self, $sym ) = @_; |
| # see if argument passed |
| if( $#_ > 0 ) { |
| # yes -- use value passed |
| $self->{bearing} = $sym; |
| }else{ |
| # no -- set to true |
| $self->{bearing} = 1; |
| } |
| } |
| |
| =head2 set_defaults |
| |
| Sets the defaults for the new method. Call with key, value pairs similar to |
| new. |
| |
| $Geo::Ellipsoid->set_defaults( |
| units => 'degrees', |
| ellipsoid => 'GRS80', |
| distance_units => 'kilometer', |
| longitude => 1, |
| bearing => 0 |
| ); |
| |
| Keys and string values (except for the ellipsoid identifier) may be shortened |
| to their first three letters and are case-insensitive: |
| |
| $Geo::Ellipsoid->set_defaults( |
| uni => 'deg', |
| ell => 'GRS80', |
| dis => 'kil', |
| lon => 1, |
| bea => 0 |
| ); |
| |
| =cut |
| |
| sub set_defaults |
| { |
| my $self = shift; |
| my %args = @_; |
| foreach my $key ( keys %args ) { |
| if( $key =~ /^ell/i ) { |
| $defaults{ellipsoid} = uc $args{$key}; |
| }elsif( $key =~ /^uni/i ) { |
| $defaults{units} = $args{$key}; |
| }elsif( $key =~ /^dis/i ) { |
| $defaults{distance_units} = $args{$key}; |
| }elsif( $key =~ /^lon/i ) { |
| $defaults{longitude} = $args{$key}; |
| }elsif( $key =~ /^bea/i ) { |
| $defaults{bearing} = $args{$key}; |
| }else{ |
| croak("Geo::Ellipsoid::set_defaults called with invalid key: $key"); |
| } |
| } |
| print "Defaults set to ($defaults{ellipsoid},$defaults{units}\n" |
| if $DEBUG; |
| } |
| |
| =head2 scales |
| |
| Returns a list consisting of the distance unit per angle of latitude |
| and longitude (degrees or radians) at the specified latitude. |
| These values may be used for fast approximations of distance |
| calculations in the vicinity of some location. |
| |
| ( $lat_scale, $lon_scale ) = $geo->scales($lat0); |
| $x = $lon_scale * ($lon - $lon0); |
| $y = $lat_scale * ($lat - $lat0); |
| |
| =cut |
| |
| sub scales |
| { |
| my $self = shift; |
| my $units = $self->{units}; |
| my $lat = $_[0]; |
| if( defined $lat ) { |
| $lat /= $degrees_per_radian if( $units eq 'degrees' ); |
| }else{ |
| carp("scales() method requires latitude argument; assuming 0"); |
| $lat = 0; |
| } |
| |
| my $aa = $self->{equatorial}; |
| my $bb = $self->{polar}; |
| my $a2 = $aa*$aa; |
| my $b2 = $bb*$bb; |
| my $d1 = $aa * cos($lat); |
| my $d2 = $bb * sin($lat); |
| my $d3 = $d1*$d1 + $d2*$d2; |
| my $d4 = sqrt($d3); |
| my $n1 = $aa * $bb; |
| my $latscl = ( $n1 * $n1 ) / ( $d3 * $d4 * $self->{conversion} ); |
| my $lonscl = ( $aa * $d1 ) / ( $d4 * $self->{conversion} ); |
| |
| if( $DEBUG ) { |
| print "lat=$lat, aa=$aa, bb=$bb\nd1=$d1, d2=$d2, d3=$d3, d4=$d4\n"; |
| print "latscl=$latscl, lonscl=$lonscl\n"; |
| } |
| |
| if( $self->{units} eq 'degrees' ) { |
| $latscl /= $degrees_per_radian; |
| $lonscl /= $degrees_per_radian; |
| } |
| return ( $latscl, $lonscl ); |
| } |
| |
| =head2 range |
| |
| Returns the range in distance units between two specified locations given |
| as latitude, longitude pairs. |
| |
| my $dist = $geo->range( $lat1, $lon1, $lat2, $lon2 ); |
| my $dist = $geo->range( @origin, @destination ); |
| |
| =cut |
| |
| sub range |
| { |
| my $self = shift; |
| my @args = _normalize_input($self->{units},@_); |
| my($range,$bearing) = _inverse($self,@args); |
| print "inverse(@_[1..4]) returns($range,$bearing)\n" if $DEBUG; |
| return $range; |
| } |
| |
| =head2 bearing |
| |
| Returns the bearing in degrees or radians from the first location to |
| the second. Zero bearing is true north. |
| |
| my $bearing = $geo->bearing( $lat1, $lon1, $lat2, $lon2 ); |
| |
| =cut |
| |
| sub bearing |
| { |
| my $self = shift; |
| my $units = $self->{units}; |
| my @args = _normalize_input($units,@_); |
| my($range,$bearing) = _inverse($self,@args); |
| print "inverse(@args) returns($range,$bearing)\n" if $DEBUG; |
| my $t = $bearing; |
| $self->_normalize_output('bearing',$bearing); |
| print "_normalize_output($t) returns($bearing)\n" if $DEBUG; |
| return $bearing; |
| } |
| |
| |
| =head2 at |
| |
| Returns the list (latitude,longitude) in degrees or radians that is a |
| specified range and bearing from a given location. |
| |
| my( $lat2, $lon2 ) = $geo->at( $lat1, $lon1, $range, $bearing ); |
| |
| =cut |
| |
| sub at |
| { |
| my $self = shift; |
| my $units = $self->{units}; |
| my( $lat, $lon, $az ) = _normalize_input($units,@_[0,1,3]); |
| my $r = $_[2]; |
| print "at($lat,$lon,$r,$az)\n" if $DEBUG; |
| my( $lat2, $lon2 ) = _forward($self,$lat,$lon,$r,$az); |
| print "_forward returns ($lat2,$lon2)\n" if $DEBUG; |
| $self->_normalize_output('longitude',$lon2); |
| $self->_normalize_output('latitude',$lat2); |
| return ( $lat2, $lon2 ); |
| } |
| |
| =head2 to |
| |
| In list context, returns (range, bearing) between two specified locations. |
| In scalar context, returns just the range. |
| |
| my( $dist, $theta ) = $geo->to( $lat1, $lon1, $lat2, $lon2 ); |
| my $dist = $geo->to( $lat1, $lon1, $lat2, $lon2 ); |
| |
| =cut |
| |
| sub to |
| { |
| my $self = shift; |
| my $units = $self->{units}; |
| my @args = _normalize_input($units,@_); |
| print "to($units,@args)\n" if $DEBUG; |
| my($range,$bearing) = _inverse($self,@args); |
| print "to: inverse(@args) returns($range,$bearing)\n" if $DEBUG; |
| #$bearing *= $degrees_per_radian if $units eq 'degrees'; |
| $self->_normalize_output('bearing',$bearing); |
| if( wantarray() ) { |
| return ( $range, $bearing ); |
| }else{ |
| return $range; |
| } |
| } |
| |
| =head2 displacement |
| |
| Returns the (x,y) displacement in distance units between the two specified |
| locations. |
| |
| my( $x, $y ) = $geo->displacement( $lat1, $lon1, $lat2, $lon2 ); |
| |
| NOTE: The x and y displacements are only approximations and only valid |
| between two locations that are fairly near to each other. Beyond 10 kilometers |
| or more, the concept of X and Y on a curved surface loses its meaning. |
| |
| =cut |
| |
| sub displacement |
| { |
| my $self = shift; |
| print "displacement(",join(',',@_),"\n" if $DEBUG; |
| my @args = _normalize_input($self->{units},@_); |
| print "call _inverse(@args)\n" if $DEBUG; |
| my( $range, $bearing ) = _inverse($self,@args); |
| print "disp: _inverse(@args) returns ($range,$bearing)\n" if $DEBUG; |
| my $x = $range * sin($bearing); |
| my $y = $range * cos($bearing); |
| return ($x,$y); |
| } |
| |
| =head2 location |
| |
| Returns the list (latitude,longitude) of a location at a given (x,y) |
| displacement from a given location. |
| |
| my @loc = $geo->location( $lat, $lon, $x, $y ); |
| |
| =cut |
| |
| sub location |
| { |
| my $self = shift; |
| my $units = $self->{units}; |
| my($lat,$lon,$x,$y) = @_; |
| my $range = sqrt( $x*$x+ $y*$y ); |
| my $bearing = atan2($x,$y); |
| $bearing *= $degrees_per_radian if $units eq 'degrees'; |
| print "location($lat,$lon,$x,$y,$range,$bearing)\n" if $DEBUG; |
| return $self->at($lat,$lon,$range,$bearing); |
| } |
| |
| ######################################################################## |
| # |
| # 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 the |
| # current specified distance unit and bearing in radians. |
| |
| sub _inverse() |
| { |
| my $self = shift; |
| my( $lat1, $lon1, $lat2, $lon2 ) = (@_); |
| print "_inverse($lat1,$lon1,$lat2,$lon2)\n" if $DEBUG; |
| |
| my $a = $self->{equatorial}; |
| my $f = $self->{flattening}; |
| |
| 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; |
| |
| if( $DEBUG ) { |
| printf "lat1=%.8f, lon1=%.8f\n", $lat1, $lon1; |
| printf "lat2=%.8f, lon2=%.8f\n", $lat2, $lon2; |
| printf "r=%.8f, tu1=%.8f, tu2=%.8f\n", $r, $tu1, $tu2; |
| printf "faz=%.8f, dlon=%.8f\n", $faz, $dlon; |
| } |
| |
| my $x = $dlon; |
| my $cnt = 0; |
| print "enter loop:\n" if $DEBUG; |
| my( $c2a, $c, $cx, $cy, $cz, $d, $del, $e, $sx, $sy, $y ); |
| do { |
| printf " x=%.8f\n", $x if $DEBUG; |
| $sx = sin($x); |
| $cx = cos($x); |
| $tu1 = $cu2*$sx; |
| $tu2 = $baz - ($su1*$cu2*$cx); |
| |
| printf " sx=%.8f, cx=%.8f, tu1=%.8f, tu2=%.8f\n", |
| $sx, $cx, $tu1, $tu2 if $DEBUG; |
| |
| $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; |
| } |
| |
| printf " sy=%.8f, cy=%.8f, y=%.8f, sa=%.8f\n", $sy, $cy, $y, $sa |
| if $DEBUG; |
| |
| $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; |
| |
| if( $DEBUG ) { |
| printf " c2a=%.8f, cz=%.8f\n", $c2a, $cz; |
| printf " e=%.8f, d=%.8f\n", $e, $d; |
| printf " (d-x)=%.8g\n", $del; |
| } |
| |
| }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; |
| |
| if( $DEBUG ) { |
| printf "e=%.8f, cy=%.8f, x=%.8f\n", $e, $cy, $x; |
| printf "sy=%.8f, c=%.8f, d=%.8f\n", $sy, $c, $d; |
| printf "cz=%.8f, a=%.8f, r=%.8f\n", $cz, $a, $r; |
| } |
| |
| $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; |
| |
| printf "s=%.8f\n", $s if $DEBUG; |
| |
| # adjust azimuth to (0,360) or (-180,180) as specified |
| if( $self->{symmetric} ) { |
| $faz += $twopi if $faz < -(pi); |
| $faz -= $twopi if $faz >= pi; |
| }else{ |
| $faz += $twopi if $faz < 0; |
| $faz -= $twopi if $faz >= $twopi; |
| } |
| |
| # return result |
| my @disp = ( ($s/$self->{conversion}), $faz ); |
| print "disp = (@disp)\n" if $DEBUG; |
| return @disp; |
| } |
| |
| # _forward |
| # |
| # Calculate the location (latitue,longitude) of a point |
| # given a starting point and a displacement from that |
| # point as (range,bearing) |
| # |
| sub _forward |
| { |
| my $self = shift; |
| my( $lat1, $lon1, $range, $bearing ) = @_; |
| |
| if( $DEBUG ) { |
| printf "_forward(lat1=%.8f,lon1=%.8f,range=%.8f,bearing=%.8f)\n", |
| $lat1, $lon1, $range, $bearing; |
| } |
| |
| my $eps = 0.5e-13; |
| |
| my $a = $self->{equatorial}; |
| my $f = $self->{flattening}; |
| my $r = 1.0 - $f; |
| |
| my $tu = $r * sin($lat1) / cos($lat1); |
| my $faz = $bearing; |
| my $s = $self->{conversion} * $range; |
| my $sf = sin($faz); |
| my $cf = cos($faz); |
| |
| my $baz = 0.0; |
| $baz = 2.0 * atan2($tu,$cf) if( $cf != 0.0 ); |
| |
| my $cu = 1.0 / sqrt(1.0 + $tu*$tu); |
| my $su = $tu * $cu; |
| my $sa = $cu * $sf; |
| my $c2a = 1.0 - ($sa*$sa); |
| my $x = 1.0 + sqrt( (((1.0/($r*$r)) - 1.0 )*$c2a) +1.0); |
| $x = ($x-2.0)/$x; |
| my $c = 1.0 - $x; |
| $c = ((($x*$x)/4.0) + 1.0)/$c; |
| my $d = $x * ((0.375*$x*$x)-1.0); |
| $tu = (($s/$r)/$a)/$c; |
| my $y = $tu; |
| |
| if( $DEBUG ) { |
| printf "r=%.8f, tu=%.8f, faz=%.8f\n", $r, $tu, $faz; |
| printf "baz=%.8f, sf=%.8f, cf=%.8f\n", $baz, $sf, $cf; |
| printf "cu=%.8f, su=%.8f, sa=%.8f\n", $cu, $su, $sa; |
| printf "x=%.8f, c=%.8f, y=%.8f\n", $x, $c, $y; |
| } |
| |
| my( $cy, $cz, $e, $sy ); |
| do { |
| $sy = sin($y); |
| $cy = cos($y); |
| $cz = cos($baz+$y); |
| $e = (2.0*$cz*$cz)-1.0; |
| $c = $y; |
| $x = $e * $cy; |
| $y = (2.0 * $e) - 1.0; |
| $y = ((((((((($sy*$sy*4.0)-3.0)*$y*$cz*$d)/6.0)+$x)*$d)/4.0)-$cz)*$sy*$d) + |
| $tu; |
| } while( abs($y-$c) > $eps ); |
| |
| $baz = ($cu*$cy*$cf) - ($su*$sy); |
| $c = $r*sqrt(($sa*$sa) + ($baz*$baz)); |
| $d = $su*$cy + $cu*$sy*$cf; |
| my $lat2 = atan2($d,$c); |
| $c = $cu*$cy - $su*$sy*$cf; |
| $x = atan2($sy*$sf,$c); |
| $c = (((((-3.0*$c2a)+4.0)*$f)+4.0)*$c2a*$f)/16.0; |
| $d = (((($e*$cy*$c) + $cz)*$sy*$c)+$y)*$sa; |
| my $lon2 = $lon1 + $x - (1.0-$c)*$d*$f; |
| #$baz = atan2($sa,$baz) + pi; |
| |
| # return result |
| return ($lat2,$lon2); |
| |
| } |
| |
| # _normalize_input |
| # |
| # Normalize a set of input angle values by converting to |
| # radians if given in degrees and by converting to the |
| # range [0,2pi), i.e. greater than or equal to zero and |
| # less than two pi. |
| # |
| sub _normalize_input |
| { |
| my $units = shift; |
| my @args = @_; |
| return map { |
| $_ = deg2rad($_) if $units eq 'degrees'; |
| while( $_ < 0 ) { $_ += $twopi } |
| while( $_ >= $twopi ) { $_ -= $twopi } |
| $_ |
| } @args; |
| } |
| |
| # _normalize_output |
| # |
| # Normalize a set of output angle values by converting to |
| # degrees if needed and by converting to the range [-pi,+pi) or |
| # [0,2pi) as needed. |
| # |
| sub _normalize_output |
| { |
| my $self = shift; |
| my $elem = shift; # 'bearing' or 'longitude' |
| # adjust remaining input values by reference |
| for ( @_ ) { |
| if( $self->{$elem} ) { |
| # normalize to range [-pi,pi) |
| while( $_ < -(pi) ) { $_ += $twopi } |
| while( $_ >= pi ) { $_ -= $twopi } |
| }else{ |
| # normalize to range [0,2*pi) |
| while( $_ < 0 ) { $_ += $twopi } |
| while( $_ >= $twopi ) { $_ -= $twopi } |
| } |
| $_ = rad2deg($_) if $self->{units} eq 'degrees'; |
| } |
| } |
| |
| =head1 DEFINED ELLIPSOIDS |
| |
| The following ellipsoids are defined in Geo::Ellipsoid, with the |
| semi-major axis in meters and the reciprocal flattening as shown. |
| The default ellipsoid is WGS84. |
| |
| Ellipsoid Semi-Major Axis (m.) 1/Flattening |
| --------- ------------------- --------------- |
| AIRY 6377563.396 299.3249646 |
| AIRY-MODIFIED 6377340.189 299.3249646 |
| AUSTRALIAN 6378160.0 298.25 |
| BESSEL-1841 6377397.155 299.1528128 |
| CLARKE-1880 6378249.145 293.465 |
| EVEREST-1830 6377276.345 290.8017 |
| EVEREST-MODIFIED 6377304.063 290.8017 |
| FISHER-1960 6378166.0 298.3 |
| FISHER-1968 6378150.0 298.3 |
| GRS80 6378137.0 298.25722210088 |
| HOUGH-1956 6378270.0 297.0 |
| HAYFORD 6378388.0 297.0 |
| IAU76 6378140.0 298.257 |
| KRASSOVSKY-1938 6378245.0 298.3 |
| NAD27 6378206.4 294.9786982138 |
| NWL-9D 6378145.0 298.25 |
| SOUTHAMERICAN-1969 6378160.0 298.25 |
| SOVIET-1985 6378136.0 298.257 |
| WGS72 6378135.0 298.26 |
| WGS84 6378137.0 298.257223563 |
| |
| =head1 LIMITATIONS |
| |
| The methods should not be used on points which are too near the poles |
| (above or below 89 degrees), and should not be used on points which |
| are antipodal, i.e., exactly on opposite sides of the ellipsoid. The |
| methods will not return valid results in these cases. |
| |
| =head1 ACKNOWLEDGEMENTS |
| |
| The conversion algorithms used here are Perl translations of Fortran |
| routines written by LCDR S<L. Pfeifer> NGS Rockville MD that implement |
| S<T. Vincenty's> Modified Rainsford's method with Helmert's elliptical |
| terms as published in "Direct and Inverse Solutions of Ellipsoid on |
| the Ellipsoid with Application of Nested Equations", S<T. Vincenty,> |
| Survey Review, April 1975. |
| |
| The Fortran source code files inverse.for and forward.for |
| may be obtained from |
| |
| ftp://ftp.ngs.noaa.gov/pub/pcsoft/for_inv.3d/source/ |
| |
| =head1 AUTHOR |
| |
| Jim Gibson, C<< <Jim@Gibson.org> >> |
| |
| =head1 BUGS |
| |
| See LIMITATIONS, above. |
| |
| Please report any bugs or feature requests to |
| C<bug-geo-ellipsoid@rt.cpan.org>, or through the web interface at |
| L<http://rt.cpan.org/NoAuth/ReportBug.html?Queue=Geo-Ellipsoid>. |
| |
| =head1 COPYRIGHT & LICENSE |
| |
| Copyright 2005-2008 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. |
| |
| =head1 SEE ALSO |
| |
| Geo::Distance, Geo::Ellipsoids |
| |
| =cut |
| |
| 1; # End of Geo::Ellipsoid |