Astro::Coords::Equatorial - Manipulate equatorial coordinates

  $c = new Astro::Coords::Equatorial( name => 'blah',
				      ra   => '05:22:56',
				      dec  => '-26:20:40.4',
				      type => 'B1950'
				      units=> 'sexagesimal');

  $c = new Astro::Coords::Equatorial( name => 'Vega',
                                      ra => ,
                                      dec => ,
                                      type => 'J2000',
                                      units => 'sex',
                                      pm => [ 0.202, 0.286],
                                      parallax => 0.13,
                                      epoch => 2004.529,
                                      );

  $c = new Astro::Coords( ra => '16h24m30.2s',
                          dec => '-00d54m2s',
                          type => 'J2000',
                          rv => 31,
                          vdefn => 'RADIO',
                          vframe => 'LSRK' );

package Astro::Coords::Equatorial;

use 5.006;
use strict;
use warnings;
use warnings::register;
use Carp;

our $VERSION = '0.04';

use Astro::SLA ();
use base qw/ Astro::Coords /;

use overload '""' => "stringify", fallback => 1;

sub new {
  my $proto = shift;
  my $class = ref($proto) || $proto;

  my %args = @_;

  return undef unless exists $args{type};

  # make sure we are upper cased.
  $args{type} = uc($args{type});

  # Convert input args to radians
  $args{ra} = Astro::Coords::Angle::Hour->to_radians($args{ra}, $args{units} )
    if exists $args{ra};
  $args{dec} = Astro::Coords::Angle->to_radians($args{dec}, $args{units} )
    if exists $args{dec};
  $args{long} = Astro::Coords::Angle->to_radians($args{long}, $args{units} )
    if exists $args{long};
  $args{lat} = Astro::Coords::Angle->to_radians($args{lat}, $args{units} )
    if exists $args{lat};

  # Default values for parallax and proper motions
  my( $pm, $parallax );
  if( exists( $args{parallax} ) ) {
    $parallax = $args{parallax};
  } else {
    $parallax = 0;
  }
  if( exists( $args{pm} ) ) {
    $pm = $args{pm};
  } else {
    $pm = [0,0];
  }

  # Try to sort out what we have been given. We need to convert
  # everything to FK5 J2000
  croak "Proper motions are supplied but not as a ref to array"
    unless ref($pm) eq 'ARRAY';

  # Extract the proper motions into convenience variables
  my $pm1 = $pm->[0];
  my $pm2 = $pm->[1];

  my ($ra, $dec, $native);

  if ($args{type} =~ /^j([0-9\.]+)/i) {
    return undef unless exists $args{ra} and exists $args{dec};
    return undef unless defined $args{ra} and defined $args{dec};

    $native = 'radec';

    $ra = $args{ra};
    $dec = $args{dec};

# The equinox is everything after the J.
    my $equinox = $1;

# Wind the RA/Dec to J2000 if the equinox isn't 2000.
    if( $equinox != 2000 ) {
      Astro::SLA::slaPreces( 'FK5', $equinox, '2000.0', $ra, $dec );
    }

# Get the epoch. If it's not given (in $args{epoch}) then it's
# the same as the equinox.
    my $epoch = ( ( exists( $args{epoch} ) && defined( $args{epoch} ) ) ?
                  $args{epoch} :
                  $equinox );

# Wind the RA/Dec to epoch 2000.0 if the epoch isn't 2000.0,
# taking the proper motion and parallax into account.
    if( $epoch != 2000 &&
        ( $pm1 != 0 || $pm2 != 0 || $parallax != 0 ) ) {
      # Assume we are HEL without checking
      my $rv = ( exists $args{rv} && $args{rv} ? $args{rv} : 0);

      my ( $ra0, $dec0 );
      Astro::SLA::slaPm( $ra, $dec,
                         Astro::SLA::DAS2R * $pm1,
                         Astro::SLA::DAS2R * $pm2,
                         Astro::SLA::DAS2R * $parallax,
			 $rv,
                         $epoch, # input epoch
			 2000.0, # output epoch
                         $ra0,
                         $dec0 );
      $ra = $ra0;
      $dec = $dec0;
    }

  } elsif ($args{type} =~ /^b([0-9\.]+)/i) {
    return undef unless exists $args{ra} and exists $args{dec};
    return undef unless defined $args{ra} and defined $args{dec};

    $native = 'radec1950';
    $ra = $args{ra};
    $dec = $args{dec};

# The equinox is everything after the B.
    my $equinox = $1;

# Get the epoch. If it's not given (in $args{epoch}) then it's
# the same as the equinox. Assume supplied epoch is Besselian
    my $epoch = ( ( exists( $args{epoch} ) && defined( $args{epoch} ) ) ?
                  $args{epoch} :
                  $equinox );

    my ( $ra0, $dec0 );

# For the implementation details, see section 4.1 of SUN/67.
    if( $pm1 != 0 || $pm2 != 0 || $parallax != 0 ) {
      # Assume we are HEL without checking
      my $rv = ( exists $args{rv} && $args{rv} ? $args{rv} : 0);

      # We are converting to J2000 but we need to convert that to Besselian epoch
      Astro::SLA::slaPm( $ra, $dec,
                         Astro::SLA::DAS2R * $pm1,
                         Astro::SLA::DAS2R * $pm2,
                         Astro::SLA::DAS2R * $parallax,
                         $rv,
                         $epoch,
                         Astro::SLA::slaEpco('B','J',2000.0), # Besselian epoch
                         $ra0,
                         $dec0 );
      $ra = $ra0;
      $dec = $dec0;
    }

    if( $equinox != 1950 ) {

# Remove the E-terms for the specified Besselian equinox
      my ( $ra0, $dec0 );
      Astro::SLA::slaSubet( $ra, $dec, $equinox, $ra0, $dec0 );
      $ra = $ra0;
      $dec = $dec0;

# Wind the RA/Dec to B1950 if the equinox isn't 1950.
      Astro::SLA::slaPreces( 'FK4', $equinox, 1950.0, $ra, $dec );

# Add the E-terms back in.
      Astro::SLA::slaAddet( $ra, $dec, 1950.0, $ra0, $dec0 );
      $ra = $ra0;
      $dec = $dec0;
    }

# Convert to J2000, no proper motion. We need the epoch at which the
# coordinate was valid
    Astro::SLA::slaFk45z($ra, $dec,
			 $epoch,
                         $ra0, $dec0
                        );
    $ra = $ra0;
    $dec = $dec0;

  } elsif ($args{type} eq "GALACTIC") {
    $native = 'glonglat';
    return undef unless exists $args{long} and exists $args{lat};
    return undef unless defined $args{long} and defined $args{lat};

    Astro::SLA::slaGaleq( $args{long}, $args{lat}, $ra, $dec);

  } elsif ($args{type} eq "SUPERGALACTIC") {
    return undef unless exists $args{long} and exists $args{lat};
    return undef unless defined $args{long} and defined $args{lat};

    $native = 'sglonglat';
    Astro::SLA::slaSupgal( $args{long}, $args{lat}, my $glong, my $glat);
    Astro::SLA::slaGaleq( $glong, $glat, $ra, $dec);

  } else {
    my $type = (defined $args{type} ? $args{type} : "<undef>");
    croak "Supplied coordinate type [$type] not recognized";
  }

  # Now the actual object
  my $c = bless { ra2000 => new Astro::Coords::Angle::Hour($ra, units => 'rad', range => '2PI'),
		  dec2000 => new Astro::Coords::Angle($dec, units => 'rad'),
		  name => $args{name},
		  pm => $args{pm}, parallax => $args{parallax}
		}, $class;

  # Specify the native encoding
  $c->native( $native );

  # Now set the velocity parameters
  if (exists $args{redshift}) {
    $c->_set_redshift( $args{redshift} );
  } else {
    $c->_set_rv( $args{rv} ) if exists $args{rv};
    $c->_set_vdefn( $args{vdefn} ) if exists $args{vdefn};
    $c->_set_vframe( $args{vframe} ) if exists $args{vframe};
  }

  return $c;
}

sub radec {
  my $self = shift;
  my ($sys, $equ) = $self->_parse_equinox( shift || 'J2000' );

  # If we have proper motions we need to take them into account
  # Do this using slaPm rather than via the base class since it
  # must be more efficient than going through apparent
  my @pm = $self->pm;
  my $par = $self->parallax;

  # First convert to J2000 current epoch

  # Fix PM array and parallax if none-defined
  @pm = (0,0) unless @pm;
  $par = 0 unless defined $par;

  # J2000 Epoch 2000.0
  my ($ra,$dec) = $self->radec2000();
  if ($pm[0] != 0 || $pm[1] != 0 || $par != 0) {
    # We have proper motions
    # Radial velocity in HEL frame
    # Note that we need to calculate the RA/Dec to get the HEL frame
    # if the radial velocity is not already in HEL
    # We have to ignore it for now and only use rv if it is 
    # already heliocentric
    my $rv = 0;
    $rv = $self->rv if $self->vframe eq 'HEL';

    # Correct for proper motion
    Astro::SLA::slaPm( $ra, $dec, Astro::SLA::DAS2R * $pm[0], 
		       Astro::SLA::DAS2R * $pm[1], $par, $rv, 2000.0,
		       Astro::SLA::slaEpj($self->_mjd_tt), $ra, $dec );

    # Convert to Angle objects
    $ra = new Astro::Coords::Angle::Hour( $ra, units => 'rad', range => '2PI');
    $dec = new Astro::Coords::Angle( $dec, units => 'rad' );
  }

  # Return it if we have the right answer
  if ($sys eq 'FK5' && $equ == 2000.0) {
    # Already have the right answer
  } elsif ($sys eq 'FK5') {
    # Preces to new equinox
    Astro::SLA::slaPreces( 'FK5', 2000.0, $equ, $ra, $dec );

  } else {
    # Convert to BYYYY
    ($ra, $dec) = $self->_j2000_to_byyyy( $equ, $ra, $dec);

  }

  return (new Astro::Coords::Angle::Hour($ra, units => 'rad', range => '2PI'),
	  new Astro::Coords::Angle($dec, units => 'rad'));

}

sub ra {
  my $self = shift;
  my %opt = @_;
  my ($ra, $dec) = $self->radec;
  my $retval = $ra->in_format( $opt{format} );

  # Tidy up array to remove sign
  shift(@$retval) if ref($retval) eq "ARRAY";
  return $retval;
}

sub dec {
  my $self = shift;
  my %opt = @_;
  my ($ra, $dec) = $self->radec;
  return $dec->in_format( $opt{format} );
}

sub radec2000 {
  my $self = shift;
  return ($self->ra2000, $self->dec2000);
}

sub ra2000 {
  my $self = shift;
  my %opt = @_;
  my $ra = $self->{ra2000};
  my $retval = $ra->in_format( $opt{format} );

  # Tidy up array
  shift(@$retval) if ref($retval) eq "ARRAY";
  return $retval;
}

sub dec2000 {
  my $self = shift;
  my %opt = @_;
  my $dec = $self->{dec2000};
  return $dec->in_format( $opt{format} );
}

sub parallax {
  my $self = shift;
  if (@_) {
    $self->{parallax} = shift;
  }
  return $self->{parallax};
}

sub pm {
  my $self = shift;
  if (@_) {
    my $pm1 = shift;
    my $pm2 = shift;
    if (!defined $pm1) {
      warnings::warnif("Proper motion 1 not defined. Using 0.0 arcsec/year");
      $pm1 = 0.0;
    }
    if (!defined $pm2) {
      warnings::warnif("Proper motion 2 not defined. Using 0.0 arcsec/year");
      $pm2 = 0.0;
    }
    $self->{pm} = [ $pm1, $pm2 ];
  }
  if( !defined( $self->{pm} ) ) { $self->{pm} = []; }
  return @{ $self->{pm} };
}

sub apparent {
  my $self = shift;

  # Assumes that Parallax and proper motions are constants for this object
  my ($ra_app, $dec_app) = $self->_cache_read( "RA_APP", "DEC_APP" );

  if (!defined $ra_app || !defined $dec_app) {

    my $ra = $self->ra2000;
    my $dec = $self->dec2000;
    my $mjd = $self->_mjd_tt;
    my $par = $self->parallax;
    my @pm = $self->pm;

    @pm = (0,0) unless @pm;
    $par = 0.0 unless defined $par;

    # do not attempt to correct for radial velocity unless we are doing parallax or
    # proper motion correction
    my $rv = 0;
    if ($par != 0 || $pm[0] != 0 || $pm[1] != 0 ) {
      # Radial velocity in HEL frame
      # Note that we need to calculate the apparent RA/Dec to get the HEL frame
      # if the radial velocity is not already in HEL
      # We have to ignore it for now and only use rv if it is heliocentric
      $rv = $self->rv if $self->vframe eq 'HEL';
    }

    Astro::SLA::slaMap( $ra, $dec,
			Astro::SLA::DAS2R * $pm[0],
			Astro::SLA::DAS2R * $pm[1], $par, $rv, 2000.0, $mjd,
			$ra_app, $dec_app);

    # Convert from observed to apparent place
    #  Astro::SLA::slaOap("r", $ra_app, $dec_app, $mjd, 0.0, $long, $lat,
    #                     0.0,0.0,0.0,
    #                     0.0,0.0,0.0,0.0,0.0,$ra, $dec);

    $ra_app = new Astro::Coords::Angle::Hour($ra_app, units => 'rad', range => '2PI');
    $dec_app = new Astro::Coords::Angle($dec_app, units => 'rad');

    $self->_cache_write( "RA_APP" => $ra_app, "DEC_APP" => $dec_app );
  }

  return ($ra_app, $dec_app);
}

sub array {
  my $self = shift;
  my ($ra, $dec) = $self->radec2000;
  return ( $self->type, $ra->radians, $dec->radians,
	   undef, undef, undef, undef, undef, undef, undef, undef);
}

sub type {
  return "RADEC";
}

sub stringify {
  my $self = shift;
  my ($ra, $dec) = $self->radec();
  return "$ra $dec";
}

sub summary {
  my $self = shift;
  my $name = $self->name;
  $name = '' unless defined $name;
  my ($ra, $dec) = $self->radec;

  return sprintf("%-16s  %-12s  %-13s  J2000",$name,$ra, $dec);
}

sub set_vel_pars {
  my $self = shift;
  my ($rv, $vdefn, $vframe) = @_;

  $self->_set_rv( $rv ) if defined $rv;
  $self->_set_vdefn( $vdefn ) if defined $vdefn;
  $self->_set_vframe( $vframe ) if defined $vframe;

  return;
}

1;