package Astrometry;
#
# implemented functions:
#
#  * RA2Deg ...... convert RA from hms to deg
#  * RA2hms ...... convert RA from string to hms
#  * Dec2Deg ..... convert Dec from dms to deg
#  * Dec2hms ..... convert Dec from string to hms
#  * HA2hms ...... convert HA from string to hms
#  * MJD ......... calculate the modified Julian date
#  * JD .......... calculate the Julian date
#
# still to come:
#  - Bessel2Julian ... convert Besselian epoch to Julian epoch
#  - Julian2Bessel ... convert Julian epoch to Besselian epoch
#  - Precess ......... precess coordinates (in degrees) between epochs
#  - Deg2RA
#  - Deg2Dec

require Exporter;

@EXPORT = qw( Dec2Deg Dec2hms HA2hms MJD Precess RA2Deg RA2hms );

######################################################################
# Split RA string into hour, minutes and seconds
######################################################################
sub RA2hms
{
if ( $#_ < 0 ) {return (0,0,0);}	# return 0 if called without parameters
my $ra = $_[0];				# get input value
my @a = ( $ra =~ /[^0-9]*([0-9]+)[^0-9]+([0-9]+)[^0-9]+([0-9]+.?[0-9]*)/ );
					# split into hms fields
if ( $#a == 0 ) {return ($a[0],0,0);}	# only hour
if ( $#a == 1 ) {return ($a[0],$a[1],0);}# hour and minutes
return ($a[0],$a[1],$a[2]);		# hour, minutes and seconds
}

######################################################################
# Calculation of RA in degrees from hour, minutes and seconds
######################################################################
sub RA2Deg
{
if ($#_ < 0) {return 0;}
if ($#_ == 0) {@_ = RA2hms( $_[0] );}
#if ($#_ == 0) {return $_[0]*15;}
if ($#_ == 1) {return 15*($_[0] + $_[1]/60);}
if ($#_ == 2) {return 15*($_[0] + $_[1]/60 + $_[2]/3600);}
if ($#_ > 2) {die "Too many arguments for RA2Deg";}
}

######################################################################
# Split Dec string into hour, minutes and seconds
######################################################################
sub Dec2hms
{
if ( $#_ < 0 ) {return (0,0,0);}	# return 0 if called without parameters
my $dec = $_[0];			# get input value
my @a = ( $dec =~ /[^\+\-0-9]*([\+\-]?[0-9]+)[^0-9]+([0-9]+)[^0-9]+([0-9]+.?[0-9]*)/ );
					# split into hms fields
if ( $#a == 0 ) {return ($a[0],0,0);}	# only hour
if ( $#a == 1 ) {return ($a[0],$a[1],0);}# hour and minutes
return ($a[0],$a[1],$a[2]);		# hour, minutes and seconds
}

######################################################################
# Calculation of Declination in degrees from degree, minutes and seconds
######################################################################
sub Dec2Deg
{
if ($#_ < 0) {return 0;}
if ($#_ == 0) {@_ = Dec2hms( $_[0] );}
my $deg=$_[0];
$deg=~s/ //g;
if ($#_ == 0) {return $deg;}
my $neg=1;
if ($deg < 0) {$deg=-$deg;$neg=-1;};
if ($#_ == 1) {return $neg*($deg + $_[1]/60);}
if ($#_ == 2) {return $neg*($deg + $_[1]/60 + $_[2]/3600);}
if ($#_ > 2) {die "Too many arguments for Dec2Deg";}
}

######################################################################
# Split HA string into hour, minutes and seconds
######################################################################
sub HA2hms
{
if ( $#_ < 0 ) {return (0,0,0);}	# return 0 if called without parameters
my $ha = $_[0];				# get input value
my @a = ( $ha =~ /[^\+\-0-9]*([\+-]?[0-9]+)[^0-9]+([0-9]+)[^0-9]+([0-9]+.?[0-9]*)/ );
					# split into hms fields
if ( $#a == 0 ) {return ($a[0],0,0);}	# only hour
if ( $#a == 1 ) {return ($a[0],$a[1],0);}# hour and minutes
return ($a[0],$a[1],$a[2]);		# hour, minutes and seconds
}

######################################################################
# calculation of the modified Julian date
# MJD = JD - 2400000.5
# ie. time since 17. November 1858 0h.
#
# * dates after 1. January 4713 BC
# * julian  calendar   until  4. October 1582
# * gregorian calendar since 15. October 1582
######################################################################
sub MJD
  {
  ($year, $month, $day, $hour, $minute, $second) = (0,0,0,0,0,0);
  ($year) = @_ if ($#_ == 0);
  ($year, $month) = @_ if ($#_ == 1);
  ($year, $month, $day) = @_ if ($#_ == 2);
  ($year, $month, $day, $hour) = @_ if ($#_ == 3);
  ($year, $month, $day, $hour, $minute) = @_ if ($#_ == 4);
  ($year, $month, $day, $hour, $minute, $second) = @_ if ($#_ == 5);
  $a = 10000 * $year + 100 * $month + $day;
  if ($month <= 2)
    {
    $month += 12;
    $year--;
    }
  if ($a <= 15821004.1)
    {
    $b = -2 + (int ($year + 4716) / 4) - 1179;
    }
  else
    {
    $b = (int $year/400) - (int $year/100) + (int $year/4);
    }
  $a = 365 * $year - 679004;
  $time = ( $hour + ($minute + ($second/60))/60 )/24;
  $mjd = $a + $b + (int 30.6001*($month + 1)) +
         $day + $time;
  return ($mjd);
  }

######################################################################
# calculation of the Julian date
######################################################################
sub JD
  {
  return (MJD(@_) + 2400000.5);
  }

######################################################################
# calculation change in coordinates due to precession
# Parameters: RA, Declination, epoch from, epoch to
# RA and Dec have to be in degrees
# if no "to" epoch is specified J2000.0 is assumed
######################################################################
sub Precess
  {
  use POSIX;
  $pi =3.141592654;
  $pi2=6.283185307;

  ($rai, $deci) = @_;
  $rai=$rai*$pi/180;
  $deci=$deci*$pi/180;
  if ($#_ == 2)
    {
    $from=$_[2];
    $to=2000;
    }
  if ($#_ == 3)
    {
    $from=$_[2];
    $to=$_[3];
    }

##########
#convert coordinates to J2000.0
##########
  if ($from != 2000)
    {
#calculte T=(t-2000.0)/100
    $tdt=($from-2000)/100;
    $tdt2=$tdt * $tdt;
    $tdt3=$tdt2 * $tdt;
#calculate angles Xi,z and Theta
    $tzeta  = 0.6406161*$tdt + 0.0000839*$tdt2 + 0.0000050*$tdt3;
    $tz     = 0.6406161*$tdt + 0.0003041*$tdt2 + 0.0000051*$tdt3;
    $ttheta = 0.5567530*$tdt + 0.0001185*$tdt2 + 0.0000116*$tdt3;
#convert to radians
    $fzeta  = $fzeta*$pi/180;
    $fz     = $fz*$pi/180;
    $ftheta = $ftheta*$pi/180;

    $cosdec = POSIX::cos $deci;
    $sindec = POSIX::sin $deci;
    $costtheta = POSIX::cos $ttheta;
    $sinttheta = POSIX::sin $ttheta;

    $dec0 = POSIX::asin ($costtheta*$sindec - (POSIX::cos $rai-$tz)*$sinttheta*$cosdec);
    $sinra0zeta = ((POSIX::sin $rai-$tz)*$cosdec) / (POSIX::cos $dec0);
    $cosra0zeta = ((POSIX::cos $rai-$tz)*$costtheta*$cosdec + $sinttheta*$sindec) / (POSIX::cos $dec0);

    $asin = POSIX::asin $sinra0zeta;
    $acos = POSIX::acos $cosra0zeta;
    $ra0zeta = $acos;
    if ($asin < 0) {$ra0zeta = $pi2 - $ra0zeta};
    $ra0=$ra0zeta - $tzeta;
    }
  else
    {
    $ra0=$rai;
    $dec0=$deci;
    }

##########
#convert coordinates from J2000.0
##########
  if ($to != 2000)
    {
#calculte T=(t-2000.0)/100
    $tdt=($to-2000)/100;
    $tdt2=$tdt * $tdt;
    $tdt3=$tdt2 * $tdt;
#calculate angles Xi,z and Theta
    $fzeta  = 0.6406161*$tdt + 0.0000839*$tdt2 + 0.0000050*$tdt3;
    $fz     = 0.6406161*$tdt + 0.0003041*$tdt2 + 0.0000051*$tdt3;
    $ftheta = 0.5567530*$tdt + 0.0001185*$tdt2 + 0.0000116*$tdt3;
#convert to radians
    $fzeta  = $fzeta*$pi/180;
    $fz     = $fz*$pi/180;
    $ftheta = $ftheta*$pi/180;

    $cosdec0 = POSIX::cos $dec0;
    $sindec0 = POSIX::sin $dec0;
    $cosftheta = POSIX::cos $ftheta;
    $sinftheta = POSIX::sin $ftheta;

#
#use sin(d) = cos(a0+Zeta)*sin(Theta)*cos(d0) + cos(Theta)*sin(d0)
#
    $dec = POSIX::asin ($cosftheta*$sindec0 + (POSIX::cos $ra0+$fzeta)*$sinftheta*$cosdec0);
    $sinraz = ((POSIX::sin $ra0+$fzeta)*$cosdec0) / (POSIX::cos $dec);
    $cosraz = ((POSIX::cos $ra0+$fzeta)*$cosftheta*$cosdec0 - $sinftheta*$sindec0) / (POSIX::cos $dec);

    $asin = POSIX::asin $sinraz;
    $acos = POSIX::acos $cosraz;
    $raz = $acos;
    if ($asin < 0) {$raz = $pi2 - $raz}
    $ra=$raz + $fz;
    }
  else
    {
    $ra=$ra0;
    $dec=$dec0;
    }

  if ($ra < 0) {$ra += $pi2}
  return ($ra*180/$pi, $dec*180/$pi);
  }