Astro
view release on metacpan or search on metacpan
Astro/Time.pm view on Meta::CPAN
package Astro::Time;
use strict;
=head1 NAME
Astro::Time - Time based astronomical routines
=head1 SYNOPSIS
use Astro::Time;
$dayno = cal2dayno($day, $month, $year);
print "It's a leap year!\n" if (leap($year));
$lmst = mjd2lst($mjd, $longitude, $dUT1);
$turns = str2turn($string, 'H');
$str = turn2str($turn, 'D', $sig);
=head1 DESCRIPTION
Astro::Time contains an assorted set Perl routines for time based
conversions, such as conversion between calendar dates and Modified
Julian day and conversion of UT to local sidereal time. Include are
routines for conversion between numerical and string representation of
angles.
=head1 AUTHOR
Chris Phillips Chris.Phillips@csiro.au
=head1 FUNCTIONS
=cut
BEGIN {
use Exporter ();
use vars qw($VERSION @ISA @EXPORT @EXPORT_OK @EXPORT_FAIL
$PI $StrSep $StrZero $Quiet );
$VERSION = '1.22';
@ISA = qw(Exporter);
@EXPORT = qw( cal2dayno dayno2cal leap yesterday tomorrow
mjd2cal cal2mjd mjd2dayno dayno2mjd now2mjd mjd2epoch
jd2mjd mjd2jd mjd2time mjd2vextime mjd2weekday mjd2weekdaystr
gst mjd2lst cal2lst dayno2lst rise lst2mjd
turn2str deg2str rad2str str2turn str2deg str2rad
hms2time time2hms month2str str2month
deg2rad rad2deg turn2rad rad2turn turn2deg deg2turn
$PI );
@EXPORT_OK = qw ( daynoOK monthOK dayOK utOK nint $StrSep $StrZero
$Quiet);
@EXPORT_FAIL = qw ( @days @MonthShortStr @MonthStr @WeekShortStr @WeekStr);
use Carp;
use POSIX qw( fmod floor ceil acos );
}
$PI = 3.1415926535897932384626433832795028841971693993751;
$StrZero = 0;
$StrSep = ':';
my $debug = 0; # Used for debugging str2turn
$Quiet = 0;
my @days = (31,28,31,30,31,30,31,31,30,31,30,31);
my @MonthShortStr = ('Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug',
'Sep', 'Oct', 'Nov', 'Dec');
my @MonthStr = ('January', 'February', 'March', 'April', 'May', 'June', 'July',
'August', 'September','October', 'November', 'December');
my @WeekShortStr = ('Mon', 'Tue', 'Wed', 'Thu', 'Fri', 'Sat', 'Sun');
my @WeekStr = ('Monday', 'Tuesday', 'Wednesday', 'Thursday', 'Friday',
'Saturday', 'Sunday');
# Is the dayno valid?
sub daynoOK ($$) {
my ($dayno, $year) = @_;
Astro/Time.pm view on Meta::CPAN
$dayno = 366;
} else {
$dayno = 365;
}
}
if (scalar(@_)==2) {
return ($dayno, $year);
} else {
($day, $month) = dayno2cal($dayno, $year);
return($day, $month, $year);
}
}
=item B<tomorrow>
($dayno, $year) = tomorrow($dayno, $year);
($day, $month, $year) = tomorrow($day, $month, $year);
Advances the day number by one, taking account of year wraps.
$dayno Day number of year
$year Year
$month Month
$day Day of month
=cut
# Verified
sub tomorrow($$;$) {
my ($day, $month, $dayno, $year);
if (scalar(@_)==2) {
($dayno, $year) = @_;
return undef if (!daynoOK($dayno, $year));
} else {
($day, $month, $year) = @_;
return undef if (!monthOK($month));
return undef if (!dayOK($day, $month, $year));
$dayno = cal2dayno($day, $month, $year);
}
$dayno++;
if (($dayno==366 && !leap($year)) || $dayno==367) {
$year++;
$dayno = 1;
}
if (scalar(@_)==2) {
return ($dayno, $year);
} else {
($day, $month) = dayno2cal($dayno, $year);
return($day, $month, $year);
}
}
=item B<mjd2cal>
($day, $month, $year, $ut) = mjd2cal($mjd);
Converts a modified Julian day number into calendar date (universal
time). (based on the slalib routine sla_djcl).
$mjd Modified Julian day (JD-2400000.5)
$day Day of the month.
$month Month of the year.
$year Year
$ut UT day fraction
=cut
# VERIFIED
sub mjd2cal($) {
my $mjd = shift;
my $ut = fmod($mjd,1.0);
if ($ut<0.0) {
$ut += 1.0;
$mjd -= 1;
}
use integer; # Calculations require integer division and modulation
# Get the integral Julian Day number
my $jd = nint($mjd + 2400001);
# Do some rather cryptic calculations
my $temp1 = 4*($jd+((6*(((4*$jd-17918)/146097)))/4+1)/2-37);
my $temp2 = 10*((($temp1-237)%1461)/4)+5;
my $year = $temp1/1461-4712;
my $month =(($temp2/306+2)%12)+1;
my $day = ($temp2%306)/10+1;
return($day, $month, $year, $ut);
}
=item B<cal2mjd>
$mjd = cal2mjd($day, $month, $year, $ut);
Converts a calendar date (universal time) into modified Julian day
number.
$day Day of the month.
$month Month of the year.
$year Year
$ut UT dayfraction
$mjd Modified Julian day (JD-2400000.5)
=cut
# Verified
sub cal2mjd($$$;$) {
my($day, $month, $year, $ut) = @_;
$ut = 0.0 if (!defined $ut);
return undef if (!monthOK($month));
return undef if (!dayOK($day, $month, $year));
return undef if (!utOK($ut));
my ($m, $y);
if ($month <= 2) {
$m = int($month+9);
$y = int($year-1);
} else {
$m = int($month-3);
$y = int($year);
}
my $c = int($y/100);
$y = $y-$c*100;
my $x1 = int(146097.0*$c/4.0);
my $x2 = int(1461.0*$y/4.0);
my $x3 = int((153.0*$m+2.0)/5.0);
return($x1+$x2+$x3+$day-678882+$ut);
}
=item B<mjd2dayno>
($dayno, $year, $ut) = mjd2dayno($mjd);
Converts a modified Julian day number into year and dayno (universal
time).
$mjd Modified Julian day (JD-2400000.5)
$year Year
$dayno Dayno of year
=cut
# NOT Verified
sub mjd2dayno($) {
my $mjd = shift;
my ($day, $month, $year, $ut) = mjd2cal($mjd);
return (cal2dayno($day,$month,$year), $year, $ut);
}
=item B<dayno2mjd>
$mjd = dayno2mjd($dayno, $year, $ut);
Converts a dayno and year to modified Julian day
Astro/Time.pm view on Meta::CPAN
my $d=0.0000062/86400.0;
my $tu = (int($mjd)-($JULIAN_DAY_J2000-2400000.5))/$JULIAN_DAYS_IN_CENTURY;
my $sidtim = $a + $tu*($b + $tu*($e - $tu*$d));
$sidtim -= int($sidtim);
if ($sidtim < 0.0) {$sidtim += 1.0};
my $gmst = $sidtim + ($mjd - int($mjd) + $dUT1/86400.0)*$SOLAR_TO_SIDEREAL;
while ($gmst<0.0) {
$gmst += 1.0;
}
while ($gmst>1.0) {
$gmst -= 1.0;
}
return $gmst;
}
# Not verified - wrapper to gmst
=item B<mjd2lst>
$lst = mjd2lst($mjd, $longitude);
$lmst = mjd2lst($mjd, $longitude, $dUT1);
$ltst = mjd2lst($mjd, $longitude, $dUT1, $eqenx);
Converts a modified Julian day number into local sidereal time (lst),
local mean sidereal time (lmst) or local true sidereal time (ltst).
Unless high precisions is required dUT1 can be omitted (it will always
be in the range -0.5 to 0.5 seconds).
$mjd Modified Julian day (JD-2400000.5)
$longitude Longitude for which the LST is required (turns)
$dUT1 Difference between UTC and UT1 (UT1 = UTC + dUT1)(seconds)
$eqenx Equation of the equinoxes (not yet supported)
$lst Local sidereal time (turns)
$lmst Local mean sidereal time (turns)
$ltst Local true sidereal time (turns)
=cut
sub mjd2lst($$;$$) {
my ($mjd, $longitude, $dUT1, $eqenx) = @_;
my $lst = gst($mjd, $dUT1, $eqenx);
return undef if (!defined $lst);
$lst += $longitude;
while ($lst>1.0) {
$lst-= 1;
}
while ($lst < 0.0) {
$lst += 1;
}
return $lst;
}
=item B<cal2lst>
$lst = cal2lst($day, $month, $year, $ut, $longitude);
$lmst = cal2lst($day, $month, $year, $ut, $longitude, $dUT1);
$ltst = cal2lst($day, $month, $year, $ut, $longitude, $dUT1, $eqenx);
Wrapper to mjd2lst using calendar date rather than mjd
=cut
sub cal2lst($$$$$;$$) {
my ($day, $month, $year, $ut, $longitude, $dUT1, $eqenx) = @_;
my $mjd = cal2mjd($day, $month, $year, $ut);
return undef if (!defined $mjd);
return mjd2lst($mjd, $longitude, $dUT1, $eqenx);
}
=item B<dayno2lst>
$lst = dayno2lst($dayno, $year, $ut, $longitude);
$lmst = dayno2lst($dayno, $year, $ut, $longitude, $dUT1);
$ltst = dayno2lst($dayno, $year, $ut, $longitude, $dUT1, $eqenx);
Wrapper to mjd2lst using calendar date rather than mjd
=cut
sub dayno2lst($$$$;$$) {
my ($dayno, $year, $ut, $longitude, $dUT1, $eqenx) = @_;
my $mjd = dayno2mjd($dayno, $year, $ut);
return undef if (!defined $mjd);
return mjd2lst($mjd, $longitude, $dUT1, $eqenx);
}
# Not verified
=item B<rise>
($lst_rise, $lst_set) = rise($ra, $dec, $obslat, $el_limit);
Return the lst rise and set time of the given source
$lst_rise, $lst_set Rise and set time (turns)
$ra, $dec RA and Dec of source (turns)
$obslat Latitude of observatory (turns)
$el_limit Elevation limit of observatory
(turns, 0 horizontal)
Returns 'Circumpolar' if source circumpolar
Returns undef if source never rises
Uses the formula:
cos $z_limit = sin $obslat * sin $dec + cos $obslat * cos $dec
* cos $HA
where:
$z_limit is the zenith angle limit corresponding to $el_limit
$HA is the Hour Angle of the source
NOTE: For maximum accuracy source coordinated should be precessed to
the current date.
=cut
sub rise ($$$$) {
#print "rise: Got @_\n";
my ($ra, $dec, $obslat, $el_limit) = @_;
$ra = turn2rad($ra);
$dec = turn2rad($dec);
$obslat = turn2rad($obslat);
$el_limit = turn2rad($el_limit);
my $z_limit = $PI/2-$el_limit;
#print "Check it\n";
# Check whether the source ever rises or is circumpolar
my $z = acos(sin($obslat)*sin($dec) + cos($obslat)*cos($dec)); # Highest point
return (undef) if ($z>$z_limit);
#print "Got here\n";
$z = acos(sin($obslat)*sin($dec) - cos($obslat)*cos($dec)); # Lowest point
return ('Circumpolar') if ($z<$z_limit);
my $cos_ha = (cos($z_limit) - sin($obslat)*sin($dec))
( run in 3.456 seconds using v1.01-cache-2.11-cpan-39bf76dae61 )