Astro::Time::HJD - Perl extension for calculating heliocentric julian date

Astro-Time-HJD documentation Contained in the Astro-Time-HJD distribution.

Index

Code Index:

NAME

Top

Astro::Time::HJD - Perl extension for calculating heliocentric julian date adjustment.

SYNOPSIS

Top

 use Astro::Time::HJD qw( correction );
 ($correction, $jd, $hjd ) = correction( $jd, $ra, $dec );
 $correction = correction( $jd, $ra, $dec );

DESCRIPTION

Top

Given an observation date, right ascention and declination, calculates the correction to be added to the julian date of the observation to adjust it's observation time to a heliocentric julian date.

USAGE

Top

correction

(SCALAR, SCALAR, SCALAR) = correction( SCALAR, SCALAR, SCALAR )
SCALAR = correction( SCALAR, SCALAR, SCALAR )

The first form returns the heliocentric correction (in julian days), the original julian date and the adjusted julian date for the given observation time and RA/DEC location. The second form returns just the heliocentric correction. The correction should be added from the julian date of the observation. The first argument to the function is the date of the observation. The date can be either a julian date or a zero padded date time string of the form 'YYYY-MM-DDTHH:MM:SS'. The second and third arguments are the RA and DEC, respectively, given in decimal degrees.

ACCURACY

Top

The algorithms implemented have a stated accuracy of "better than .1 seconds from 1900 through 2050" per Patrick Wallace. Write me in 2049 and I'll see about updating this.

AUTHOR

Top

Robert Creager <Astro-HJD@LogicalChaos.org>

Patrick Wallace, personal communication Michael Koppelman, personal communication

Astro-Time-HJD documentation Contained in the Astro-Time-HJD distribution. 
package Astro::Time::HJD;

use 5.006;
use strict;
use warnings;

use Math::Trig qw( atan );
use Time::JulianDay;

require Exporter;
use AutoLoader qw(AUTOLOAD);

our @ISA = qw(Exporter);

our @EXPORT_OK = ( qw( correction ) );

our $VERSION = '0.02';

use constant PI  => 3.1415926535897932384626433832795028841971693993751;
use constant PI2 => PI * 2.0;
use constant DEG2RAD => PI / 180.0;
use constant AS2RAD => PI / ( 180.0 * 3600.0 );
use constant SPEED => 1.9913e-7;
use constant REMB => 3.12e-5;
use constant SEMB => 8.31e-11;

sub arctan { return atan( shift ); }

sub mod
   {
   my ( $A, $B ) = @_;
   ( ( $B ) != 0.0 ?
      ( ( $A ) * ( $B ) > 0.0 ? ( $A ) - ( $B ) * int( ( $A ) / ( $B ) ) :
      ( $A ) + ( $B ) * int( -( $A ) / ( $B ) ) ) : ( $A ) );
   }

sub floor
   {
   return int( shift );
   }

sub frac
   {
   my $a = shift;
   return $a - floor( $a );
   }

sub correction
   {
   my $date = shift;
   my $ra = ( shift ) * DEG2RAD;
   my $dec = ( shift ) * DEG2RAD;
   my( @ymd, $jd, $year );

   if ( $date =~ /^[\d.]+$/)
      {
      $jd = $date;
      @ymd = inverse_julian_day( $jd );
      $year = $ymd[0];
      }
   else
      {
      $date =~ /^(\d{4})-(\d{2})-(\d{2})T(\d{2}):(\d{2}):(\d{2})$/i;
      @ymd = ( $1, $2, $3 );
      $year = $1;
      $jd = julian_day( @ymd ) - .5 + $4 / 24.0 + $5 / 1440.0 + $6 / 86400.0;
      }

   my $bjd      = julian_day( $year, 1, 1 ) - .5;
   #
   #  First routine wants the day to start at 1, not zero
   my $day      = floor( $jd - $bjd + 1.0 );
   my $frac_day = frac( $jd - $bjd );
   #
   # The julian year should start at 0, not 1, so we subtract 1 from the
   # day here
   my $jy = $year + (($day - 1.0) + $frac_day) / 365.25;

   my @earth = earth( $year, $day, $frac_day );
   my @precmat = prec( 2000.0, $jy );
   my @v = dcs2c( $ra, $dec );
   my @star = dmxv( \@precmat, \@v );
   my $correction =
      ( $earth[ 0 ] * $star[ 0 ] + $earth[ 1 ] * $star[ 1 ] + $earth[ 2 ] *
      $star[ 2 ] ) * ( 499.004782 / 86400.0 );
   #
   # Return all the info if desired
   if ( wantarray )
      {
      return ( $correction, $jd, $jd + $correction );
      }
   #
   # otherwise, just return the correction
   else
      {
      return $correction;
      }
   }

sub earth
   {
   my ( $iy, $id, $fd ) = @_;
   my @pv;

   my $yi  = $iy - 1900.0;
   my $iy4 = $iy >= 0.0 ? $iy % 4.0 : 3.0 - ( -$iy - 1.0 ) % 4.0;
   my $yf  =
      ( ( 4.0 * ( $id - floor( 1.0 / ( $iy4 + 1.0 ) ) ) - $iy4 - 2.0 ) + 4.0 *
      $fd ) / 1461.0;
   my $t = $yi + $yf;
   my $elm = mod( 4.881627938 + PI2 * $yf + 0.00013420 * $t, PI2 );
   my $gam = 4.90823 + 3.0005e-4 * $t;
   my $em = $elm - $gam;
   my $eps0 = 0.40931975 - 2.27e-6 * $t;
   my $e   = 0.016751 - 4.2e-7 * $t;
   my $esq = $e * $e;
   my $v = $em + 2.0 * $e * sin( $em ) + 1.25 * $esq * sin( 2.0 * $em );
   my $elt = $v + $gam;
   my $r = ( 1.0 - $esq ) / ( 1.0 + $e * cos( $v ) );
   my $elmm = mod( 4.72 + 83.9971 * $t, PI2 );
   my $coselt = cos( $elt );
   my $sineps = sin( $eps0 );
   my $coseps = cos( $eps0 );
   my $w1     = -$r * sin( $elt );
   my $w2     = -( SPEED ) * ( $coselt + $e * cos( $gam ) );
   my $selmm  = sin( $elmm );
   my $celmm  = cos( $elmm );
   $pv[ 0 ] = -$r * $coselt - REMB * $celmm;
   $pv[ 1 ] = ( $w1 - REMB * $selmm ) * $coseps;
   $pv[ 2 ] = $w1 * $sineps;
   $pv[ 3 ] = SPEED * ( sin( $elt ) + $e * sin( $gam ) ) + SEMB * $selmm;
   $pv[ 4 ] = ( $w2 - SEMB * $celmm ) * $coseps;
   $pv[ 5 ] = $w2 * $sineps;
   return @pv;
   }

sub prec
   {
   my ( $ep0, $ep1 ) = @_;
   my $t0 = ( $ep0 - 2000.0 ) / 100.0;
   my $t = ( $ep1 - $ep0 ) / 100.0;
   my $tas2r = $t * AS2RAD;
   my $w     = 2306.2181 + ( ( 1.39656 - ( 0.000139 * $t0 ) ) * $t0 );
   my $zeta  =
      ( $w + ( ( 0.30188 - 0.000344 * $t0 ) + 0.017998 * $t ) * $t ) * $tas2r;
   my $z =
      ( $w + ( ( 1.09468 + 0.000066 * $t0 ) + 0.018203 * $t ) * $t ) * $tas2r;
   my $theta =
      ( ( 2004.3109 + ( -0.85330 - 0.000217 * $t0 ) * $t0 ) +
      ( ( -0.42665 - 0.000217 * $t0 ) - 0.041833 * $t ) * $t ) * $tas2r;
   return deuler( "ZYZ", -$zeta, $theta, -$z, );
   }

sub deuler
   {
   my ( $order, $phi, $theta, $psi ) = @_;
   my @rmat;
   my @result = ( [ 1.0, 0.0, 0.0 ], [ 0.0, 1.0, 0.0 ], [ 0.0, 0.0, 1.0 ] );
   my $l = length $order;
   for ( my $n = 0; $n < 3; ++$n )
      {
      if ( $n <= $l )
         {
         my @rotn = (
                      [ 1.0, 0.0, 0.0 ], [ 0.0, 1.0, 0.0 ],
                      [ 0.0, 0.0, 1.0 ]
         );
         my $angle;
         $_ = $n;
         SWITCH:
            {
            /0/ && do { $angle = $phi;   last SWITCH };
            /1/ && do { $angle = $theta; last SWITCH };
            /2/ && do { $angle = $psi;   last SWITCH };
            die "unknown value: '$_'";
            }
         my $s = sin( $angle );
         my $c = cos( $angle );
         $_ = substr $order, $n, 1;
         SWITCH:
            {
            /[x1]/i && do
               {
               $rotn[ 1 ][ 1 ] = $c;
               $rotn[ 1 ][ 2 ] = $s;
               $rotn[ 2 ][ 1 ] = -$s;
               $rotn[ 2 ][ 2 ] = $c;
               last SWITCH;
            };
            /[y2]/i && do
               {
               $rotn[ 0 ][ 0 ] = $c;
               $rotn[ 0 ][ 2 ] = -$s;
               $rotn[ 2 ][ 0 ] = $s;
               $rotn[ 2 ][ 2 ] = $c;
               last SWITCH;
            };
            /[z3]/i && do
               {
               $rotn[ 0 ][ 0 ] = $c;
               $rotn[ 0 ][ 1 ] = $s;
               $rotn[ 1 ][ 0 ] = -$s;
               $rotn[ 1 ][ 1 ] = $c;
               last SWITCH;
            };
            die "unknown character: '$_'";
            }
         my @wm;
         for ( my $i = 0; $i < 3; ++$i )
            {
            for ( my $j = 0; $j < 3; ++$j )
               {
               my $w = 0.0;
               for ( my $k = 0; $k < 3; ++$k )
                  {
                  $w += $rotn[ $i ][ $k ] * $result[ $k ][ $j ];
                  }
               $wm[ $i ][ $j ] = $w;
               }
            }
         for ( my $j = 0; $j < 3; ++$j )
            {
            for ( my $i = 0; $i < 3; ++$i )
               {
               $result[ $i ][ $j ] = $wm[ $i ][ $j ];
               }
            }
         }
      }
   for ( my $j = 0; $j < 3; ++$j )
      {
      for ( my $i = 0; $i < 3; ++$i )
         {
         $rmat[ $i ][ $j ] = $result[ $i ][ $j ];
         }
      }
   return @rmat;
   }

sub dmxv
   {
   my ( $dm, $va ) = @_;
   my ( @vb, @vw );
   for ( my $j = 0; $j < 3; ++$j )
      {
      my $w = 0.0;
      for ( my $i = 0; $i < 3; ++$i )
         {
         $w += $dm->[ $j ][ $i ] * $va->[ $i ];
         }
      $vw[ $j ] = $w;
      }
   for ( my $j = 0; $j < 3; ++$j )
      {
      $vb[ $j ] = $vw[ $j ];
      }
   return @vb;
   }

sub dcs2c
   {
   my ( $a, $b ) = @_;
   my @v;
   my $cosb = cos( $b );
   $v[ 0 ] = cos( $a ) * $cosb;
   $v[ 1 ] = sin( $a ) * $cosb;
   $v[ 2 ] = sin( $b );
   return @v;
   }

"Should you really be looking here?";
__END__