view release on metacpan or  search on metacpan

lib/Astro/Coord/ECI.pm  view on Meta::CPAN

The ellipsoid argument is the name of a L</Reference Ellipsoid> known
to the class, and is optional. If not specified, the most-recently-set
ellipsoid will be used.

The conversion from geocentric to geodetic comes from Kazimierz
Borkowski's "Accurate Algorithms to Transform Geocentric to Geodetic
Coordinates", at F<http://www.astro.uni.torun.pl/~kb/Papers/geod/Geod-BG.htm>.
This is best viewed with Internet Explorer because of its use of Microsoft's
Symbol font.

=cut

sub geodetic {
    my ($self, @args) = @_;

#	Detect and acquire the optional ellipsoid name argument. We do
#	this before the check, since the check expects the extra
#	argument to be a time.

    my $elps = (@args == 1 || @args == 4) ? pop @args : undef;

    $self = $self->_check_coord (geodetic => \@args, $elps ? $elps : ());

#	The following is just a sleazy way to get a consistent
#	error message if the ellipsoid name is unknown.

    $elps && $self->reference_ellipsoid ($elps);

#	If we're fetching the geodetic coordinates
    
    unless (@args) {

#	Return cached coordinates if they exist and we did not
#	override the default ellipsoid.

	return @{$self->{_ECI_cache}{fixed}{geodetic}}
	    if $self->{_ECI_cache}{fixed}{geodetic} && !$elps;
	$self->{debug} and do {
	    require Data::Dumper;
	    local $Data::Dumper::Terse = 1;
	    print "Debug geodetic - explicit ellipsoid ",
		Data::Dumper::Dumper( $elps );
	};

#	Get a reference to the ellipsoid data to use.

	$elps = $elps ? $known_ellipsoid{$elps} : $self;
	$self->{debug} and do {  
	    require Data::Dumper;
	    local $Data::Dumper::Terse = 1;
	    print "Debug geodetic - ellipsoid ", Data::Dumper::Dumper( $elps );
	};

#	Calculate geodetic coordinates.

	my ($phiprime, $lambda, $rho) = $self->geocentric;
	my $r = $rho * cos ($phiprime);
	my $b = $elps->{semimajor} * (1- $elps->{flattening});
	my $a = $elps->{semimajor};
	my $z = $rho * sin ($phiprime);
	# The $b is _not_ a magic variable, since we are not inside
	# any of the specialized blocks that make $b magic. Perl::Critic
	# is simply confused.
	$b = - $b	## no critic (RequireLocalizedPunctuationVars)
	    if $z < 0;	# Per Borkowski, for southern hemisphere.

#	The following algorithm is due to Kazimierz Borkowski's
#	paper "Accurate Algorithms to Transform Geocentric to Geodetic
#	Coordinates", from
#	http://www.astro.uni.torun.pl/~kb/Papers/geod/Geod-BG.htm

	my $bz = $b * $z;
	my $asq_bsq = $a * $a - $b * $b;
	my $ar = $a * $r;
	my $E = ($bz - $asq_bsq) / $ar;		# Borkowski (10)
	my $F = ($bz + $asq_bsq) / $ar;		# Borkowski (11)
	my $Q = ($E * $E - $F * $F) * 2;	# Borkowski (17)
	my $P = ($E * $F + 1) * 4 / 3;		# Borkowski (16)
	my $D = $P * $P * $P + $Q * $Q;		# Borkowski (15)
	my $v = $D >= 0 ? do {
	    my $d = sqrt $D;
	    my $onethird = 1 / 3;
	    my $vp = ($d - $Q) ** $onethird -	# Borkowski (14a)
		($d + $Q) ** $onethird;
	    $vp * $vp >= abs ($P) ? $vp :
		- ($vp * $vp * $vp + 2 * $Q) /	# Borkowski (20)
		    (3 * $P);
	} : do {
	    my $p = - $P;
	    sqrt (cos (acos ($Q /			# Borkowski (14b)
		sqrt ($p * $p * $p)) / 3) * $p) * 2;
	};
	my $G = (sqrt ($E * $E + $v)		# Borkowski (13)
		+ $E) / 2;
	my $t = sqrt ($G * $G + ($F - $v * $G)	# Borkowski (12)
		/ (2 * $G - $E)) - $G;

	# Borkowski (18)
	# equivalent to atan (arg1)
	my $phi = do {
	    my $Y = $a * ( 1 - $t * $t );
	    my $X = 2 * $b * $t;
	    ( $Y, $X ) = ( -$Y, -$X ) if $X < 0;
	    atan2 $Y, $X;
	};

	my $h = ($r - $a * $t) * cos ($phi) +	# Borkowski (19)
	    ($z - $b) * sin ($phi);

#	End of Borkowski's algorthm.

#	Cache the results of the calculation if they were done using
#	the default ellipsoid.

	$self->{_ECI_cache}{fixed}{geodetic} = [$phi, $lambda, $h] unless $elps;

#	Return the results in any event.

	$self->get ('debug') and print <<eod;
Debug geodetic: geocentric to geodetic
    inputs: