Astro::Coord::ECI::Utils - Utility routines for astronomical calculations

 use Astro::Coord::ECI::Utils qw{:all};
 my $now = time ();
 print "The current Julian day is ", julianday ($now);

 AU = number of kilometers in an astronomical unit
 JD_OF_EPOCH = the Julian Day of Perl epoch 0
 LIGHTYEAR = number of kilometers in a light year
 PARSEC = number of kilometers in a parsec
 PERL2000 = January 1 2000, 12 noon universal, in Perl time
 PI = the circle ratio, computed as atan2 (0, -1)
 PIOVER2 = half the circle ratio
 SECSPERDAY = the number of seconds in a day
 SECS_PER_SIDERIAL_DAY = seconds in a siderial day
 SPEED_OF_LIGHT = speed of light in kilometers per second
 TWOPI = twice the circle ratio

package Astro::Coord::ECI::Utils;

use strict;
use warnings;

our $VERSION = '0.039';
our @ISA = qw{Exporter};

use Carp;
## use Config;
use Data::Dumper;
use POSIX qw{floor strftime};
use Scalar::Util qw{ blessed };

my @time_routines;

BEGIN {

    # NOTE WELL
    #
    # The logic here should be consistent with the optional-module text
    # emitted by inc/Astro/Coord/ECI/Recommend.pm.
    #
    eval {
##	$Config{use64bitint} and return 0;
	require Time::y2038;
	Time::y2038->import();
	@time_routines = ( qw{ gmtime localtime timegm timelocal } );
	1;
    } or do {
	require Time::Local;
	Time::Local->import();
	@time_routines = ( qw{ timegm timelocal } );
    };
}

our @EXPORT;
our @EXPORT_OK = ( qw{
	AU $DATETIMEFORMAT $JD_GREGORIAN JD_OF_EPOCH LIGHTYEAR PARSEC
	PERL2000 PI PIOVER2 SECSPERDAY SECS_PER_SIDERIAL_DAY
	SPEED_OF_LIGHT TWOPI acos asin
	atmospheric_extinction date2epoch date2jd deg2rad distsq
	dynamical_delta embodies epoch2datetime equation_of_time
	find_first_true intensity_to_magnitude jcent2000 jd2date
	jd2datetime jday2000 julianday load_module looks_like_number max
	min mod2pi nutation_in_longitude nutation_in_obliquity obliquity
	omega rad2deg tan theta0 thetag vector_cross_product
	vector_dot_product vector_magnitude vector_unitize
       	_classisa _instance },
	@time_routines );

our %EXPORT_TAGS = (
    all => \@EXPORT_OK,
    params => [ qw{ _classisa _instance } ],
    time => \@time_routines,
);

use constant AU => 149597870;		# 1 astronomical unit, per
					# Meeus, Appendix I pg 407.
use constant LIGHTYEAR => 9.4607e12;	# 1 light-year, per Meeus,
					# Appendix I pg 407.
use constant PARSEC => 30.8568e12;	# 1 parsec, per Meeus,
					# Appendix I pg 407.
use constant PERL2000 => timegm (0, 0, 12, 1, 0, 100);
use constant PI => atan2 (0, -1);
use constant PIOVER2 => PI / 2;
use constant SECSPERDAY => 86400;
use constant SECS_PER_SIDERIAL_DAY => 86164.0905;	# Appendix I, page 408.
use constant SPEED_OF_LIGHT => 299792.458;	# KM/sec, per NIST.
### use constant SOLAR_RADIUS => 1392000 / 2;	# Meeus, Appendix I, page 407.
use constant TWOPI => PI * 2;

sub acos {
    abs ($_[0]) > 1 and confess <<eod;
Programming error - Trying to take the arc cosine of a number greater
        than 1.
eod
    return atan2 (sqrt (1 - $_[0] * $_[0]), $_[0])
}

sub asin {return atan2 ($_[0], sqrt (1 - $_[0] * $_[0]))}

#	Note that the "constant" 0.120 in Aaer (aerosol scattering) is
#	based on a compromise value A0 = 0.050 in Green's equation 3
#	(not exhibited here), which can vary from 0.035 in the winter to
#	0.065 in the summer. This makes a difference of a couple tenths
#	at 20 degrees elevation, but a couple magnitudes at the
#	horizon. Green also remarks that the 1.5 denominator in the
#	same equation (a.k.a. the scale height) can be up to twice
#	that.


sub atmospheric_extinction {
    my ($elevation, $height) = @_;
    my $cosZ = cos (PIOVER2 - $elevation);
    my $X = 1/($cosZ + 0.025 * exp (-11 * $cosZ));	# Green 1
    my $Aray = 0.1451 * exp (-$height / 7.996);	# Green 2
    my $Aaer = 0.120 * exp (-$height / 1.5);	# Green 4
    return ($Aray + $Aaer + 0.016) * $X;	# Green 5, 6
}

our $DATETIMEFORMAT;
our $JD_GREGORIAN;
BEGIN {
    $DATETIMEFORMAT = '%a %b %d %Y %H:%M:%S';
    $JD_GREGORIAN = 2299160.5;
}

sub date2jd {
    my @args = @_;
    unshift @args, 0 while @args < 6;
    my ($sec, $min, $hr, $day, $mon, $yr) = @args;
    $mon++;		# Algorithm expects month 1-12.
    $yr += 1900;	# Algorithm expects zero-based year.
    ($yr < -4712) and croak "Error - Invalid year $yr";
    ($mon < 1 || $mon > 12) and croak "Error - Invalid month $mon";
    if ($mon < 3) {
	--$yr;
	$mon += 12;
    }
    my $A = floor ($yr / 100);
    my $B = 2 - $A + floor ($A / 4);
    my $jd = floor (365.25 * ($yr + 4716)) +
	floor (30.6001 * ($mon + 1)) + $day + $B - 1524.5 +
	((($sec || 0) / 60 + ($min || 0)) / 60 + ($hr || 0)) / 24;
    $jd < $JD_GREGORIAN and
	$jd = floor (365.25 * ($yr + 4716)) +
	floor (30.6001 * ($mon + 1)) + $day - 1524.5 +
	((($sec || 0) / 60 + ($min || 0)) / 60 + ($hr || 0)) / 24;
    return $jd;
}

use constant JD_OF_EPOCH => date2jd (gmtime (0));

sub date2epoch {
    my @args = @_;
    unshift @args, 0 while @args < 6;
    my ($sec, $min, $hr, $day, $mon, $yr) = @args;
    return (date2jd ($day, $mon, $yr) - JD_OF_EPOCH) * SECSPERDAY +
    (($hr || 0) * 60 + ($min || 0)) * 60 + ($sec || 0);
}

sub deg2rad {return $_[0] * PI / 180}

sub distsq {
    my ($a, $b) = @_;
    (ref $a eq 'ARRAY' && ref $b eq 'ARRAY' && @$a == @$b) or
	confess <<eod;
Programming error - Both arguments to distsq must be  references to
        arrays of the same length.
eod

    my $sum = 0;
    my $size = @$a;
    for (my $inx = 0; $inx < $size; $inx++) {
	my $delta = $a->[$inx] - $b->[$inx];
	$sum += $delta * $delta;
    }
    return $sum
}

sub dynamical_delta {
    my ($time) = @_;
    my $year = (gmtime $time)[5] + 1900;
    my $t = ($year - 2000) / 100;
    my $correction = .37 * ($year - 2100);	# Meeus' correction to (10.2)
    return (25.3 * $t + 102) * $t + 102		# Meeus (10.2)
	    + $correction;			# Meeus' correction.
}

sub embodies {
    my ($thingy, $class) = @_;
    my $code = eval {$thingy->can('represents')};
    return $code ? $code->($thingy, $class) : undef;
}

sub epoch2datetime {
    my ($time) = @_;
    my $day = floor ($time / SECSPERDAY);
    my $sec = $time - $day * SECSPERDAY;
    ($day, my $mon, my $yr, my $greg, my $leap) = jd2date (
	my $jd = $day + JD_OF_EPOCH);
    $day = floor ($day + .5);
    my $min = floor ($sec / 60);
    $sec = $sec - $min * 60;
    my $hr = floor ($min / 60);
    $min = $min - $hr * 60;
    my $wday = ($jd + 1.5) % 7;
    my $yd = floor (275 * ($mon + 1) / 9) - (2 - $leap) * floor (($mon +
	10) / 12) + $day - 31;
    wantarray and return ($sec, $min, $hr, $day, $mon, $yr, $wday, $yd,
	0);
    return strftime ($DATETIMEFORMAT, $sec, $min, $hr, $day, $mon, $yr,
	$wday, $yd, 0);
}

sub equation_of_time {

    my $time = shift;

    my $epsilon = obliquity ($time);
    my $y = tan($epsilon / 2);
    $y *= $y;


#	The following algorithm is from Meeus, chapter 25, page, 163 ff.

    my $T = jcent2000($time);				# Meeus (25.1)
    my $L0 = mod2pi(deg2rad((.0003032 * $T + 36000.76983) * $T	# Meeus (25.2)
	    + 280.46646));
    my $M = mod2pi(deg2rad(((-.0001537) * $T + 35999.05029)	# Meeus (25.3)
	    * $T + 357.52911));
    my $e = (-0.0000001267 * $T - 0.000042037) * $T + 0.016708634;# Meeus (25.4)

    my $E = $y * sin (2 * $L0) - 2 * $e * sin ($M) +
	4 * $e * $y * sin ($M) * cos (2 * $L0) -
	$y * $y * .5 * sin (4 * $L0) -
	1.25 * $e * $e * sin (2 * $M);				# Meeus (28.3)

    return $E * SECSPERDAY / TWOPI;	# The formula gives radians.
}

sub find_first_true {
    my ($begin, $end, $test, $limit) = @_;
    $limit ||= 1;
    if ($limit >= 1) {
	if ($begin <= $end) {
	    $begin = floor ($begin);
	    $end = floor ($end) == $end ? $end : floor ($end) + 1;
	} else {
	    $end = floor ($end);
	    $begin = floor ($begin) == $begin ? $begin : floor ($begin) + 1;
	}
    }
    my $iterator = (
	$end > $begin ?
	sub {$end - $begin > $limit} :
	sub {$begin - $end > $limit}
    );
    while ($iterator->()) {
	my $mid = $limit >= 1 ?
	    floor (($begin + $end) / 2) : ($begin + $end) / 2;
	($begin, $end) = ($test->($mid)) ?
	    ($begin, $mid) : ($mid, $end);
    }
    return $end;
}

{	# Begin local symbol block
    my $intensity_to_mag_factor;	# Calculate only if needed.
    sub intensity_to_magnitude {
	return - ($intensity_to_mag_factor ||= 2.5 / log (10)) * log ($_[0]);
    }
}

sub jd2date {
    my ($time) = @_;
    my $mod_jd = $time + 0.5;
    my $Z = floor ($mod_jd);
    my $F = $mod_jd - $Z;
    my $A = (my $julian = $Z < $JD_GREGORIAN) ? $Z : do {
	my $alpha = floor (($Z - 1867216.25)/36524.25);
	$Z + 1 + $alpha - floor ($alpha / 4);
    };
    my $B = $A + 1524;
    my $C = floor (($B - 122.1) / 365.25);
    my $D = floor (365.25 * $C);
    my $E = floor (($B - $D) / 30.6001);
    my $day = $B - $D - floor (30.6001 * $E) + $F;
    my $mon = $E < 14 ? $E - 1 : $E - 13;
    my $yr = $mon > 2 ? $C - 4716 : $C - 4715;
    return ($day, $mon - 1, $yr - 1900, !$julian, ($julian ? !($yr % 4) : (
		$yr % 400 ? $yr % 100 ? !($yr % 4) : 0 : 1)));
##	% 400 ? 1 : $yr % 100 ? 0 : !($yr % 4))));
}

sub jd2datetime {
    my ($time) = @_;
    return epoch2datetime(($time - JD_OF_EPOCH) * SECSPERDAY);
}

sub jcent2000 {return jday2000 ($_[0]) / 36525}

sub jday2000 {return ($_[0] - PERL2000) / SECSPERDAY}

sub julianday {return jday2000($_[0]) + 2_451_545.0}

{	# Local symbol block. Oh, for 5.10 and state variables.
    my %error;
    my %rslt;
    sub load_module {
	my  ($module) = @_;
	exists $error{$module} and croak $error{$module};
	exists $rslt{$module} and return $rslt{$module};
	# I considered Module::Load here, but it appears not to support
	# .pmc files. No, it's not an issue at the moment, but it may be
	# if Perl 6 becomes a reality.
	$rslt{$module} = eval "require $module";
	$@ and croak ($error{$module} = $@);
	return $rslt{$module};
    }
}	# End local symbol block.

unless (eval {require Scalar::Util; Scalar::Util->import
	('looks_like_number'); 1}) {
    no warnings qw{once};
    *looks_like_number = sub {
	local $_ = shift;

	# checks from perlfaq4
	return 0 if !defined($_) || ref($_);
	return 1 if (/^[+-]?\d+$/); # is a +/- integer
	return 1 if (/^([+-]?)(?=\d|\.\d)\d*(\.\d*)?([Ee]([+-]?\d+))?$/); # a C float
	return 1 if ($] >= 5.008 and /^(Inf(inity)?|NaN)$/i)
	    or ($] >= 5.006001 and /^Inf$/i);

	return 0;
    };
}

unless (eval {require List::Util; List::Util->import ('max'); 1}) {
    no warnings qw{once};
    *max = sub {
	my $rslt;
	foreach (@_) {
	    defined $_ or next;
	    if (!defined $rslt || $_ > $rslt) {
		$rslt = $_;
	    }
	}
	$rslt;
    };
}

unless (eval {require List::Util; List::Util->import ('min'); 1}) {
    no warnings qw{once};
    *min = sub {
	my $rslt;
	foreach (@_) {
	    defined $_ or next;
	    if (!defined $rslt || $_ < $rslt) {
		$rslt = $_;
	    }
	}
	$rslt;
    };
}

sub mod2pi {return $_[0] - floor ($_[0] / TWOPI) * TWOPI}

sub nutation_in_longitude {

    my $time = shift;
    my $T = jcent2000 ($time);	# Meeus (22.1)

    my $omega = mod2pi (deg2rad ((($T / 450000 + .0020708) * $T -
	    1934.136261) * $T + 125.04452));

    my $L = mod2pi (deg2rad (36000.7698 * $T + 280.4665));
    my $Lprime = mod2pi (deg2rad (481267.8813 * $T + 218.3165));
    my $delta_psi = deg2rad ((-17.20 * sin ($omega) - 1.32 * sin (2 * $L)
	    - 0.23 * sin (2 * $Lprime) + 0.21 * sin (2 * $omega))/3600);

    return $delta_psi;
}

sub nutation_in_obliquity {

    my $time = shift;
    my $T = jcent2000 ($time);	# Meeus (22.1)

    my $omega = mod2pi (deg2rad ((($T / 450000 + .0020708) * $T -
	    1934.136261) * $T + 125.04452));

    my $L = mod2pi (deg2rad (36000.7698 * $T + 280.4665));
    my $Lprime = mod2pi (deg2rad (481267.8813 * $T + 218.3165));
    my $delta_epsilon = deg2rad ((9.20 * cos ($omega) + 0.57 * cos (2 * $L) +
	    0.10 * cos (2 * $Lprime) - 0.09 * cos (2 * $omega))/3600);

    return $delta_epsilon;
}

use constant E0BASE => (21.446 / 60 + 26) / 60 + 23;

sub obliquity {

    my $time = shift;

    my $T = jcent2000 ($time);	# Meeus (22.1)

    my $delta_epsilon = nutation_in_obliquity ($time);

    my $epsilon0 = deg2rad (((0.001813 * $T - 0.00059) * $T - 46.8150)
	    * $T / 3600 + E0BASE);
    return $epsilon0 + $delta_epsilon;
}

sub omega {
    my $T = jcent2000 (shift);	# Meeus (22.1)
    return mod2pi (deg2rad ((($T / 450000 + .0020708) * $T -
	    1934.136261) * $T + 125.04452));
}

sub rad2deg {return $_[0] / PI * 180}

sub tan {return sin ($_[0]) / cos ($_[0])}

sub theta0 {
    my ($time) = @_;
    return thetag (timegm (0, 0, 0, (gmtime $time)[3 .. 5]));
}

#	Meeus, pg 88, equation 12.4, converted to radians and Perl dates.

sub thetag {
    my ($time) = @_;
    my $T = jcent2000 ($time);
    return mod2pi (4.89496121273579 + 6.30038809898496 *
	    jday2000 ($time))
	    + (6.77070812713916e-06 - 4.5087296615715e-10 * $T) * $T * $T;
}

sub vector_cross_product {
    my ( $b, $c ) = @_;
    @{ $b } == 3 and @{ $c } == 3
	or confess 'Programming error - vector_cross_product arguments',
	    ' must be references to arrays of length 3';
    return [
	$b->[1] * $c->[2] - $b->[2] * $c->[1],
	$b->[2] * $c->[0] - $b->[0] * $c->[2],
	$b->[0] * $c->[1] - $b->[1] * $c->[0],
    ];
}

sub vector_dot_product {
    my ( $b, $c ) = @_;
    @{ $b } == @{ $c }
	or confess 'Programming error - vector_dot_product arguments ',
	    'must be references to arrays of the same length';
    my $prod = 0;
    my $size = @{ $b } - 1;
    foreach my $inx ( 0 .. $size ) {
	$prod += $b->[$inx] * $c->[$inx];
    }
    return $prod;
}

sub vector_magnitude {
    my ( $b ) = @_;
    'ARRAY' eq ref $b
	or confess 'Programming error - vector_magnitude argument ',
    'must be a reference to an array';
    my $mag = 0;
    my $size = @{ $b } - 1;
    foreach my $inx ( 0 .. $size ) {
	$mag += $b->[$inx] * $b->[$inx];
    }
    return sqrt $mag;
}

sub vector_unitize {
    my ( $b ) = @_;
    'ARRAY' eq ref $b
	or confess 'Programming error - vector_unitize argument ',
    'must be a reference to an array';
    my $mag = vector_magnitude( $b );
    return [ map { $_ / $mag } @{ $b } ];
}

#	_classisa( 'Foo', 'Bar' );
#
#	Returns true if Foo->isa( 'Bar' ) is true, and false otherwise.
#	In particular, returns false if the first argument is a
#	reference.

sub _classisa {
    my ( $invocant, $class ) = @_;
    ref $invocant and return;
    defined $invocant or return;
    return $invocant->isa( $class );
}

#	_instance( $foo, 'Bar' );
#
#	Returns true if $foo is an instance of 'Bar', and false
#	otherwise. In particular, returns false if $foo is not a
#	reference, or if it is not blessed.

sub _instance {
    my ( $object, $class ) = @_;
    ref $object or return;
    blessed( $object ) or return;
    return $object->isa( $class );
}

#	$epoch_time = _parse_time_iso_8601
#
#	Parse ISO 8601 date/time. I do not intend to expose this, since
#	it will probably go away when the satpass script is dropped. It
#	would actually be in that script except for the fact that it can
#	be more easily tested here, and because of the possibility that
#	it will be used in App::Satpass2.
{

    my %special_day_offset = (
	yesterday => -SECSPERDAY(),
	today => 0,
	tomorrow => SECSPERDAY(),
    );

    sub _parse_time_iso_8601 {
	my ( $string ) = @_;

	my @zone;
	$string =~ s/ \s* (?: ( Z ) |
				( [+-] ) ( \d{2} ) :? ( \d{2} )? ) \z //smx
	    and @zone = ( $1, $2, $3, $4 );
	my @date;

	# ISO 8601 date
	if ( $string =~ m{ \A
				( \d{4} \D? | \d{2} \D )			# year: $1
				(?: ( \d{1,2} ) \D?				# month: $2
		    		    (?: ( \d{1,2} ) (?: \s* | \D? )		# day: $3
						(?: ( \d{1,2} ) \D?			# hour: $4
			    			    (?: ( \d{1,2} ) \D?			# minute: $5
								(?: ( \d{1,2} ) \D?		# second: $6
				    				    ( \d* )			# fract: $7
								)?
			    			    )?
						)?
		    		    )?
				)?
				\z
	    	    }smx ) {
	    @date = ( $1, $2, $3, $4, $5, $6, $7, undef );

	# special-case 'yesterday', 'today', and 'tomorrow'.
	} elsif ( $string =~ m{ \A
				( yesterday | today | tomorrow )	# day: $1
				(?: \D* ( \d{1,2} ) \D?			# hour: $2
		    		    (?: ( \d{1,2} ) \D?			# minute: $3
						(?: ( \d{1,2} ) \D?		# second: $4
			    			    ( \d* )			# fract: $5
						)?
		    		    )?
				)?
				\z }smx ) {
	    my @today = @zone ? gmtime : localtime;
	    @date = ( $today[5] + 1900, $today[4] + 1, $today[3], $2, $3,
		$4, $5, $special_day_offset{$1} );

	} else {

	    return;

	}


	my $offset = pop @date || 0;
	if ( @zone && !$zone[0] ) {
	    my ( undef, $sign, $hr, $min ) = @zone;
	    $offset -= $sign . ( ( $hr * 60 + ( $min || 0 ) ) * 60 )
	}

	foreach ( @date ) {
	    defined $_ and s/ \D+ //smxg;
	}

	if ( $date[0] < 70 ) {
	    $date[0] += 100;
	} elsif ( $date[0] >= 100 ) {
	    $date[0] -= 1900;
	}
	defined $date[1] and --$date[1];
	defined $date[2] or $date[2] = 1;
	my $frc = pop @date;

	foreach ( @date ) {
	    defined $_ or $_ = 0;
	}

	my $time;
	if ( @zone ) {
	    $time = timegm( reverse @date );
	} else {
	    $time = timelocal( reverse @date );
	}

	if ( defined $frc  && $frc ne '') {
	    my $denom = 1 . ( 0 x length $frc );
	    $time += $frc / $denom;
	}

	return $time + $offset;
    }

}

1;

__END__

# ex: set textwidth=72 :