Astro::Coords::Elements - Specify astronomical coordinates using orbital elements

  $c = new Astro::Coords::Elements( elements => \%elements );

package Astro::Coords::Elements;

use 5.006;
use strict;
use warnings;
use Carp;

our $VERSION = '0.04';

# Need working slaPlante
use Astro::SLA 0.95 ();
use Astro::Coords::Angle;
use Time::Piece qw/ :override /;

use base qw/ Astro::Coords /;

use overload '""' => "stringify";

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

  my %opts = @_;

  # "elements" key must be defined and point to a reference
  return undef unless (exists $opts{elements} && ref($opts{elements}));

  # We allow use of an array reference or a hash
  # The array ref is converted to a hash
  my %elements;
  if (ref($opts{elements}) eq 'HASH') {
    %elements = %{ $opts{elements} };
  } elsif (ref($opts{elements}) eq 'ARRAY' &&
	   $opts{elements}->[0] eq 'ELEMENTS') {

    my $i = 3;
    for my $key (qw/ EPOCH ORBINC ANODE PERIH AORQ E AORL /) {
      $elements{$key} = $opts{elements}->[$i];
      $i++;
    }
    # assign EPOCHPERIH
    if (!defined $elements{AORL} && defined $opts{elements}->[$i]) {
      $elements{EPOCHPERIH} = $opts{elements}->[$i];
    } else {
      $elements{DM} = $opts{elements}->[$i];
    }

  } else {
    return undef;
  }

  # Sanity check
  for (qw/ ORBINC ANODE PERIH AORQ E/) {
    #return undef unless exists $elements{$_};
    croak "Must supply element $_ to constructor"
      unless exists $elements{$_};
  }

  # Fix up EPOCHs if it has been specified as a string
  for my $key (qw/ EPOCH EPOCHPERIH / ) {
    next unless exists $elements{$key};
    my $epoch = $elements{$key};

    # if we are missing one of the EPOCHs that is okay
    # so just skip
    if (! defined $epoch) {
      # and delete it from the hash as if it was never supplied
      # this avoids complications later
      delete $elements{$key};
      next;
    }

    if ($epoch =~ /^\d+\.\d+$/ || $epoch =~ /^\d+$/) {
      # an MJD so do not modify
    } elsif ($epoch =~ /\d\d\d\d \w\w\w \d+\.\d+/) {
      # has letters in it so try to parse
      # Split on decimal point
      my ($date, $frac) = split(/\./,$epoch,2);
      $frac = "0.". $frac; # preserve as decimal fraction
      my $format = '%Y %B %d';
      #print "EPOCH : $epoch and $date and $frac\n";
      my $obj = Time::Piece->strptime($date, $format);
      my $tzoffset = $obj->tzoffset;
      $tzoffset = $tzoffset->seconds if defined $tzoffset;
      $obj = gmtime($obj->epoch() + $tzoffset);

      # get the MJD and add on the fraction
      my $mjd = $obj->mjd() + $frac;
      $elements{$key} = $mjd;
      #print "MJD: $mjd\n";

    } else {
      # do not understand the format so return undef
      warn "Unable to recognize format for elements $key [$epoch]";
      return undef;
    }

    # Convert JD to MJD
    if ($elements{$key} > 2400000.5) {
      $elements{$key} -= 2400000.5;
    }

  }

  # but complain if we do not have one of them
  croak "Must supply one of EPOCH or EPOCHPERIH - both were undefined"
    if (!exists $elements{EPOCH} &&
	!exists $elements{EPOCHPERIH});

  # create the object
  bless { elements => \%elements, name => $opts{name} }, $class;

}

sub elements {
  my $self = shift;
  return %{ $self->{elements}};
}

sub array {
  my $self = shift;
  my %el = $self->elements;

  # use EPOCHPERIH if EPOCH is not defined
  my $epoch = $el{EPOCH};
  $epoch = $el{EPOCHPERIH} unless $epoch;

  # the 8th element can be the EPOCHPERIH or DM field
  # dependent on the element type.
  my $lastel;
  if (defined $el{EPOCHPERIH} && !defined $el{DM}
      && !defined $el{AORL}) {
    $lastel = $el{EPOCHPERIH};
  } else {
    $lastel = $el{DM};
  }

  return ( $self->type, undef, undef,
	   $epoch, $el{ORBINC}, $el{ANODE}, $el{PERIH}, 
	   $el{AORQ}, $el{E}, $el{AORL}, $lastel);
}

sub type {
  return "ELEMENTS";
}

sub stringify {
  my $self = shift;
  my %el = $self->elements;
  my $str = join(",", map { (defined $_ ? $_ : 'UNDEF' ) } values %el);
  return $str;
}

sub summary {
  my $self = shift;
  my $name = $self->name;
  $name = '' unless defined $name;
  return sprintf("%-16s  %-12s  %-11s ELEMENTS",$name,'','');
}

sub apparent {
  my $self = shift;

  my ($ra_app,$dec_app) = $self->_cache_read( "RA_APP", "DEC_APP" );

  # not in cache so must calculate it
  if (!defined $ra_app || !defined $dec_app) {

    my $tel = $self->telescope;
    my $long = (defined $tel ? $tel->long : 0.0 );
    my $lat = (defined $tel ? $tel->lat : 0.0 );
    my %el = $self->elements;
    my $jform;
    if (exists $el{DM} and defined $el{DM}) {
      # major planets
      $jform = 1;
    } elsif (exists $el{AORL} and defined $el{AORL}) {
      # minor planets
      $jform = 2;
      $el{DM} = 0;
    } else {
      # comets
      $jform = 3;
      $el{DM} = 0;
      $el{AORL} = 0;
    }

    # synch epoch if need be
    if (!exists $el{EPOCH} || !exists $el{EPOCHPERIH}) {
      if (exists $el{EPOCH}) {
	$el{EPOCHPERIH} = $el{EPOCH};
      } else {
	$el{EPOCH} = $el{EPOCHPERIH};
      }
    }

    # First have to perturb the elements to the current epoch
    # if we have a minor planet or comet
    if ( $jform == 2 || $jform == 3) {
      # for now we do not have enough information for jform=3
      # so just assume the EPOCH is the same
#    use Data::Dumper;
#    print "Before perturbing: ". Dumper(\%el);
#    print "MJD ref : " . $self->_mjd_tt . " and jform = $jform\n";
      Astro::SLA::slaPertel($jform,$el{EPOCH},$self->_mjd_tt,
			    $el{EPOCHPERIH},$el{ORBINC},$el{ANODE},
			    $el{PERIH},$el{AORQ},$el{E},$el{AORL},
			    $el{EPOCH},$el{ORBINC}, $el{ANODE},
			    $el{PERIH},$el{AORQ},$el{E},$el{AORL},
			    my $jstat);
#    print "After perturbing: " .Dumper(\%el);
      croak "Error perturbing elements for target ".
	(defined $self->name ? $self->name : '' )
	  ." [status=$jstat]" 
	    if $jstat != 0;
    }


  # Print out the values
  #print "EPOCH:  $el{EPOCH}\n";
  #print "ORBINC: ". ($el{ORBINC}*Astro::SLA::DR2D) . "\n";
  #print "ANODE:  ". ($el{ANODE}*Astro::SLA::DR2D) . "\n";
  #print "PERIH : ". ($el{PERIH}*Astro::SLA::DR2D) . "\n";
  #print "AORQ:   $el{AORQ}\n";
  #print "E:      $el{E}\n";

    Astro::SLA::slaPlante($self->_mjd_tt, $long, $lat, $jform,
			  $el{EPOCH}, $el{ORBINC}, $el{ANODE}, $el{PERIH}, 
			  $el{AORQ}, $el{E}, $el{AORL}, $el{DM}, 
			  my $ra, my $dec, my $dist, my $j);

    croak "Error determining apparent RA/Dec for target ".
      (defined $self->name ? $self->name : '' )
	."[status=$j]" if $j != 0;

    # Convert from observed to apparent place
    Astro::SLA::slaOap("r", $ra, $dec, $self->_mjd_tt, 0.0, $long, $lat, 
		       0.0,0.0,0.0,
		       0.0,0.0,0.0,0.0,0.0,$ra, $dec);

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

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

  return( $ra_app, $dec_app );
}

sub rv {
  croak "Not yet implemented element radial velocities";
}

sub vdefn {
  return 'RADIO';
}

sub vframe {
  return 'GEO';
}

1;