view release on metacpan or  search on metacpan

Astro/Time.pm  view on Meta::CPAN

sub mjd2vextime($;$) {
  my ($dayno, $year, $ut) = mjd2dayno(shift);
  my $np = shift;
  $np = 0 if (! defined $np);
  return sprintf("%dy%03dd%s", $year, $dayno, turn2str($ut, 'H', $np, 'hms'));
}

=item B<mjd2epoch>

  $time = mjd2epoch($mjd);

 Converts a Modified Julian day to unix Epoch (seconds sinve 1 Jan 1970)
 Rounded to the nearest second
    $mjd     Modified Julian day
    $tie     Seconds since 1 Jan 1970

=cut

sub mjd2epoch($) {
  my $mjd = shift;
  my $epoch = ($mjd - 40587)*24*60*60;
  return int($epoch + $epoch/abs($epoch*2)); # Work even if epoch is negative
}

=item B<gst>

  $gst  = gst($mjd);
  $gmst = gst($mjd, $dUT1);
  $gtst = gst($mjd, $dUT1, $eqenx);

 Converts a modified Julian day number to Greenwich sidereal time
   $mjd     modified Julian day (JD-2400000.5)
   $dUT1    difference between UTC and UT1 (UT1 = UTC + dUT1) (seconds)
   $eqenx   Equation of the equinoxes (not yet supported)
   $gst     Greenwich sidereal time (turns)
   $gmst    Greenwich mean sidereal time (turns)
   $gtst    Greenwich true sidereal time (turns)

=cut

# Verified
sub gst($;$$) {
  my ($mjd, $dUT1, $eqenx) = @_;
  $dUT1 = 0.0 if (! defined $dUT1);
  if ($dUT1 > 0.5 || $dUT1 < -0.5) {
    carp '$dUT1 out of range';
    return undef;
  }
  if (defined $eqenx) {
    croak '$eqenx is not supported yet';
  }

  my $JULIAN_DAY_J2000          = 2451545.0;
  my $JULIAN_DAYS_IN_CENTURY    = 36525.0;
  my $SOLAR_TO_SIDEREAL         = 1.002737909350795;

  my $a=101.0 + 24110.54841/86400.0;
  my $b=8640184.812866/86400.0;
  my $e=0.093104/86400.0;
  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