Astro::Coords - Class for handling astronomical coordinates

  use Astro::Coords;

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

  $c = new Astro::Coords( long => '05:22:56',
                          lat  => '-26:20:40.4',
                          type => 'galactic');

  $c = new Astro::Coords( planet => 'mars' );

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

  $c = new Astro::Coords( az => 345, el => 45 );

  # Associate with an observer location
  $c->telescope( new Astro::Telescope( 'JCMT' ));

  # ...and a reference epoch for all calculations
  $date = Time::Piece->strptime($string, $format);
  $c->datetime( $date );

  # or use DateTime
  $date = DateTime->from_epoch( epoch => $epoch, time_zone => 'UTC' );
  $c->datetime( $date );

  # Return coordinates J2000, for the epoch stored in the datetime
  # object. This will work for all variants.
  ($ra, $dec) = $c->radec();
  $radians = $ra->radians;

  # or individually
  $ra = $c->ra();  # returns Astro::Coords::Angle::Hour object
  $dec = $c->dec( format => 'deg' );

  # Return coordinates J2000, epoch 2000.0
  $ra = $c->ra2000();
  $dec = $c->dec2000();

  # Return coordinats apparent, reference epoch, from location
  # In sexagesimal format.
  ($ra_app, $dec_app) = $c->apparent;
  $ra_app = $c->ra_app( format => 's');
  $dec_app = $c->dec_app( format => 's' );

  # Azimuth and elevation for reference epoch from observer location
  ($az, $el) = $c->azel;
  my $az = $c->az;
  my $el = $c->el;

  # obtain summary string of object
  $summary = "$c";

  # Obtain full summary as an array
  @summary = $c->array;

  # See if the target is observable for the current time
  # and telescope
  $obs = 1 if $c->isObservable;

  # Calculate distance to another coordinate (in radians)
  $distance = $c->distance( $c2 );

  # Calculate the rise and set time of the source
  $tr = $c->rise_time;
  $ts = $c->set_time;

  # transit elevation
  $trans = $c->transit_el;

  # transit time
  $mtime = $c->meridian_time();

  + Equatorial RA/Dec, galactic (including proper motions and parallax)
  + Planets
  + Comets/Asteroids
  + Fixed locations in azimuth and elevations
  + interpolated apparent coordinates

  $radians = $c->ra( format => 'rad' );
  $string  = $c->ra( format => 'sex' );
  $deg     = $c->ra( format => 'deg' );
  $object  = $c->ra();

  SUN_RISE_SET: Position of Sun for sunrise or sunset (-50 arcminutes)
  CIVIL_TWILIGHT: Civil twilight (-6 degrees)
  NAUT_TWILIGHT: Nautical twilight (-12 degrees)
  AST_TWILIGHT: Astronomical twilight (-18 degrees)

  $set = $c->set_time( horizon => Astro::Coords::AST_TWILIGHT );

package Astro::Coords;

use 5.006;
use strict;
use warnings;
use warnings::register;
use Carp;
use vars qw/ $DEBUG /;
$DEBUG = 0;

our $VERSION = '0.12';

use Math::Trig qw/ acos /;
use Astro::SLA ();
use Astro::Coords::Angle;
use Astro::Coords::Angle::Hour;
use Astro::Coords::Equatorial;
use Astro::Coords::Elements;
use Astro::Coords::Planet;
use Astro::Coords::Interpolated;
use Astro::Coords::Fixed;
use Astro::Coords::Calibration;

use Scalar::Util qw/ blessed /;
use DateTime;
use Time::Piece;

# Constants for Sun rise/set and twilight definitions
# Elevation in radians
# See http://aa.usno.navy.mil/faq/docs/RST_defs.html
use constant SUN_RISE_SET => ( - (50 * 60) * Astro::SLA::DAS2R); # 50 arcmin
use constant CIVIL_TWILIGHT => ( - (6 * 3600) * Astro::SLA::DAS2R); # 6 deg
use constant NAUT_TWILIGHT => ( - (12 * 3600) * Astro::SLA::DAS2R); # 12 deg
use constant AST_TWILIGHT => ( - (18 * 3600) * Astro::SLA::DAS2R); # 18 deg

# This is a fudge. Not accurate
use constant MOON_RISE_SET => ( 5 * 60 * Astro::SLA::DAS2R);

# Number of km in one Astronomical Unit
use constant AU2KM => 149.59787066e6;

# Speed of light ( km/s )
use constant CLIGHT => 2.99792458e5;

# UTC TimeZone
use constant UTC => new DateTime::TimeZone( name => 'UTC' );

sub new {
  my $class = shift;

  my %args = @_;

  my $obj;

  # Always try for a planet object first if $args{planet} is used
  # (it might be that ra/dec are being specified and planet is a target
  # name - this allows all the keys to be specified at once and the
  # object can decide the most likely coordinate object to use
  # This has the distinct disadvantage that planet is always tried
  # even though it is rare. We want to be able to throw anything
  # at this without knowing what we are.
  if (exists $args{planet} and defined $args{planet}) {
    $obj = new Astro::Coords::Planet( $args{planet} );
  }

  # planet did not work. Try something else.
  unless (defined $obj) {

    # For elements we must not only check for the elements key
    # but also make sure that that key points to a hash containing
    # at least the EPOCH or EPOCHPERIH key
    if (exists $args{elements} and defined $args{elements}
       && (UNIVERSAL::isa($args{elements},"HASH") # A hash
       &&  ((exists $args{elements}{EPOCH}
       and defined $args{elements}{EPOCH})
       ||  (exists $args{elements}{EPOCHPERIH}
       and defined $args{elements}{EPOCHPERIH}))) ||
	( ref($args{elements}) eq 'ARRAY' && # ->array() ref
	  $args{elements}->[0] eq 'ELEMENTS')
     ) {

      $obj = new Astro::Coords::Elements( %args );

    } elsif (exists $args{mjd1}) {

      $obj = new Astro::Coords::Interpolated( %args );

    } elsif (exists $args{type} and defined $args{type}) {

      $obj = new Astro::Coords::Equatorial( %args );

    } elsif (exists $args{az} or exists $args{el} or exists $args{ha}) {

      $obj = new Astro::Coords::Fixed( %args );

    } elsif ( scalar keys %args == 0 ) {

      $obj = new Astro::Coords::Calibration();

    } else {
    # unable to work out what you are asking for
      return undef;

    }
  }

  return $obj;
}

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

sub telescope {
  my $self = shift;
  if (@_) { 
    my $tel = shift;
    return undef unless UNIVERSAL::isa($tel, "Astro::Telescope");
    $self->{Telescope} = $tel;
    # Telescope is part of the caching scheme
    $self->_calc_cache_key();
  }
  return $self->{Telescope};
}

sub datetime {
  my $self = shift;
  if (@_) {
    my $time = shift;

    # undef is okay
    croak "datetime: Argument does not have an mjd() method [class="
      . ( ref($time) ? ref($time) : $time) ."]"
      if (defined $time && !UNIVERSAL::can($time, "mjd"));
    $self->{DateTime} = _clone_time( $time );
    $self->usenow(0);
    # Changing usenow will have forced a recalculation of the cache key
    # so no update is required here
    # $self->_calc_cache_key;
  }
  if (defined $self->{DateTime} && ! $self->usenow) {
    return $self->{DateTime};
  } else {
    return DateTime->now( time_zone => UTC );
  }
}

sub has_datetime {
  my $self = shift;
  return (defined $self->{DateTime});
}

sub usenow {
  my $self = shift;
  if (@_) {
    $self->{UseNow} = shift;
    # Since this affects caching we need to force a key check if this value
    # is updated
    $self->_calc_cache_key();
  }
  return $self->{UseNow};
}

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

sub comment {
  my $self = shift;
  if (@_) {
    $self->{Comment} = shift;
  }
  my $com = $self->{Comment};
  $com = '' unless defined $com;
  return $com;
}

sub native {
  my $self = shift;
  if (@_) {
    $self->{NativeMethod} = shift;
  }
  return (defined $self->{NativeMethod} ? $self->{NativeMethod} : 'apparent' );
}

sub azel {
  my $self = shift;

  my ($az,$el) = $self->_cache_read( "AZ", "EL" );

  if (!defined $az || !defined $el) {

    my $ha = $self->ha;
    my $dec = $self->dec_app;
    my $tel = $self->telescope;
    my $lat = ( defined $tel ? $tel->lat : 0.0);
    Astro::SLA::slaDe2h( $ha, $dec, $lat, $az, $el );
    $az = new Astro::Coords::Angle( $az, units => 'rad', range => '2PI' );
    $el = new Astro::Coords::Angle( $el, units => 'rad' );

    $self->_cache_write( "AZ" => $az, "EL" => $el );
  }
  return ($az, $el);
}

sub ra_app {
  my $self = shift;
  my %opt = @_;

  my $ra = $self->_cache_read( "RA_APP" );
  if (!defined $ra) {
    $ra = ($self->apparent)[0];
  }
  return $ra->in_format( $opt{format} );
}

sub dec_app {
  my $self = shift;
  my %opt = @_;
  my $dec = $self->_cache_read( "DEC_APP" );
  if (!defined $dec) {
    $dec = ($self->apparent)[1];
  }
  return $dec->in_format( $opt{format} );
}

# normalize key was supported but should its absence imply no normalization?

sub ha {
  my $self = shift;
  my %opt = @_;

  my $ha = $self->_cache_read( "HA" );
  if (!defined $ha) {
    $ha = $self->_lst - $self->ra_app;
    # Always normalize?
    $ha = new Astro::Coords::Angle::Hour( $ha, units => 'rad', range => 'PI' );
  }
  return $ha->in_format( $opt{format} );
}

sub az {
  my $self = shift;
  my %opt = @_;
  my $az = $self->_cache_read( "AZ" );
  if (!defined $az) {
    ($az, my $el) = $self->azel();
  }
  return $az->in_format( $opt{format} );
}

sub el {
  my $self = shift;
  my %opt = @_;
  my $el = $self->_cache_read( "EL" );

  if (!defined $el) {
    (my $az, $el) = $self->azel();
  }
  return $el->in_format( $opt{format} );
}

sub airmass {
  my $self = shift;
  my $el = $self->el;
  my $zd = Astro::SLA::DPIBY2 - $el;
  return Astro::SLA::slaAirmas( $zd );
}

sub radec {
  my $self = shift;
  my ($sys, $equ) = $self->_parse_equinox( shift || 'J2000' );
  my ($ra_app, $dec_app) = $self->apparent;
  my $mjd = $self->_mjd_tt;
  my ($rm, $dm);
  if ($sys eq 'FK5') {
    # Julian epoch
    Astro::SLA::slaAmp($ra_app, $dec_app, $mjd, $equ, $rm, $dm);
  } elsif ($sys eq 'FK4') {
    # Convert to J2000 and then convert to Besselian epoch
    Astro::SLA::slaAmp($ra_app, $dec_app, $mjd, 2000.0, $rm, $dm);

    ($rm, $dm) = $self->_j2000_to_byyyy( $equ, $rm, $dm);
  }

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

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

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

sub glong {
  my $self = shift;
  my %opt = @_;
  my ($glong,$glat) = $self->glonglat();
  return $glong->in_format( $opt{format} );
}

sub glat {
  my $self = shift;
  my %opt = @_;
  my ($glong,$glat) = $self->glonglat();
  return $glat->in_format( $opt{format} );
}

sub sglong {
  my $self = shift;
  my %opt = @_;
  my ($sglong,$sglat) = $self->sglonglat();
  return $sglong->in_format( $opt{format} );
}

sub sglat {
  my $self = shift;
  my %opt = @_;
  my ($sglong,$sglat) = $self->sglonglat();
  return $sglat->in_format( $opt{format} );
}

sub ecllong {
  my $self = shift;
  my %opt = @_;
  my ($eclong,$eclat) = $self->ecllonglat();
  return $eclong->in_format( $opt{format} );
}

sub ecllat {
  my $self = shift;
  my %opt = @_;
  my ($eclong,$eclat) = $self->ecllonglat();
  return $eclat->in_format( $opt{format} );
}

sub glonglat {
  my $self = shift;
  my ($ra,$dec) = $self->radec;
  Astro::SLA::slaEqgal( $ra, $dec, my $long, my $lat );
  return (new Astro::Coords::Angle($long, units => 'rad', range => '2PI'),
	  new Astro::Coords::Angle($lat, units => 'rad'));
}

sub sglonglat {
  my $self = shift;
  my ($glong, $glat) = $self->glonglat();
  Astro::SLA::slaGalsup( $glong, $glat, my $sglong, my $sglat);
  return (new Astro::Coords::Angle($sglong, units => 'rad', range => '2PI'),
	  new Astro::Coords::Angle($sglat, units => 'rad'));
}

sub ecllonglat {
  my $self = shift;
  my ($ra, $dec) = $self->radec;
  Astro::SLA::slaEqecl( $ra, $dec, $self->_mjd_tt, my $long, my $lat );
  return (new Astro::Coords::Angle($long, units => 'rad', range => '2PI'),
	  new Astro::Coords::Angle($lat, units => 'rad'));
}

sub radec2000 {
  my $self = shift;

  # store current configuration
  my $reftime = ( $self->has_datetime ? $self->datetime : undef);

  # Create new time
  $self->datetime( DateTime->new( year => 2000, month => 1,
				  day => 1, hour => 12) );

  # Ask for the answer
  my ($ra, $dec) = $self->radec( 'J2000' );

  # restore the date state
  $self->datetime( $reftime );

  return ($ra, $dec);
}

sub radec1950 {
  my $self = shift;
  my ($ra, $dec) = $self->radec;

  # No E-terms or precession since we are going to B1950 epoch 1950
  Astro::SLA::slaFk54z($ra,$dec,1950.0,my $r1950, 
		       my $d1950, my $dr1950, my $dd1950);

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

sub pa {
  my $self = shift;
  my %opt = @_;
  my $ha = $self->ha;
  my $dec = $self->dec_app;
  my $tel = $self->telescope;
  my $lat = ( defined $tel ? $tel->lat : 0.0);
  my $pa = Astro::SLA::slaPa($ha, $dec, $lat);
  $pa = new Astro::Coords::Angle( $pa, units => 'rad' );
  return $pa->in_format( $opt{format} );
}

sub isObservable {
  my $self = shift;

  # Get the telescope
  my $tel = $self->telescope;
  return 0 unless defined $tel;

  # Get the limits hash
  my %limits = $tel->limits;

  if (exists $limits{type}) {

    if ($limits{type} eq 'AZEL') {

      # Get the current elevation of the source
      my $el = $self->el;

      if ($el > $limits{el}{min} and $el < $limits{el}{max}) {
	return 1;
      } else {
	return 0;
      }

    } elsif ($limits{type} eq 'HADEC') {

      # Get the current HA
      my $ha = $self->ha( normalize => 1 );

      if ( $ha > $limits{ha}{min} and $ha < $limits{ha}{max}) {
	my $dec= $self->dec_app;

	if ($dec > $limits{dec}{min} and $dec < $limits{dec}{max}) {
	  return 1;
	} else {
	  return 0;
	}

      } else {
	return 0;
      }

    } else {
      # have no idea
      return 0;
    }

  } else {
    return 0;
  }

}

sub array {
  my $self = shift;
  croak "The method array() must be subclassed\n";
}

sub distance {
  my $self = shift;
  my $offset = shift;

  my( $ra, $dec ) = $self->radec;
  my( $ra_off, $dec_off ) = $offset->radec;

  Astro::SLA::slaDs2tp($ra_off, $dec_off,
		       $ra, $dec,
		       my $xi, my $eta, my $j);

  return () unless $j == 0;

  if (wantarray) {
    return (new Astro::Coords::Angle($xi, units => 'rad'),
	    new Astro::Coords::Angle($eta, units => 'rad'));
  } else {
    my $dist = ($xi**2 + $eta**2)**0.5;
    return new Astro::Coords::Angle( $dist, units => 'rad' );
  }
}

sub status {
  my $self = shift;
  my $string;

  $string .= "Target name:    " . $self->name . "\n"
    if $self->name;

  $string .= "Coordinate type:" . $self->type ."\n";

  if ($self->type ne 'CAL') {

    my ($az,$el) = $self->azel;
    $string .= "Elevation:      " . $el->degrees ." deg\n";
    $string .= "Azimuth  :      " . $az->degrees ." deg\n";
    my $ha = $self->ha->hours;
    $string .= "Hour angle:     " . $ha ." hrs\n";
    my ($ra_app, $dec_app) = $self->apparent;
    $string .= "Apparent RA :   " . $ra_app->string . "\n";
    $string .= "Apparent dec:   " . $dec_app->string ."\n";

    # Transit time
    my $mt = $self->meridian_time;
    $string .= "Time of next transit: " . 
      (defined $mt ? $mt->datetime : "<never>") ."\n";

    my $t_el = $self->transit_el(format=>'d');
    $string .= "Transit El:     " . (defined $t_el ? "$t_el deg"
				     : "<undefined>") ."\n";
    my $ha_set = $self->ha_set( format => 'hour');
    $string .= "Hour Ang. (set):" . (defined $ha_set ? $ha_set : '??')." hrs\n";

    my $t = $self->rise_time;
    $string .= "Next Rise time:      " . $t->datetime . "\n" if defined $t;
    $t = $self->set_time;
    $string .= "Next Set time:       " . $t->datetime . "\n" if defined $t;

    # This check was here before we added a RA/Dec to the
    # base class.
    if ($self->can('radec')) {
      my ($ra, $dec) = $self->radec;
      $string .= "RA (J2000):     " . $ra->string . "\n";
      $string .= "Dec(J2000):     " . $dec->string . "\n";
    }
  }

  if (defined $self->telescope) {
    my $name = (defined $self->telescope->fullname ?
		$self->telescope->fullname : $self->telescope->name );
    $string .= "Telescope:      $name\n";
    if ($self->isObservable) {
      $string .= "The target is currently observable\n";
    } else {
      $string .= "The target is not currently observable\n";
    }
  }

  $string .= "For time ". $self->datetime->datetime ."\n";
  my $fmt = 's';
  $string .= "LST: ". $self->_lst->hours ."\n";

  return $string;
}

sub calculate {
  my $self = shift;

  my %opt = @_;

  croak "No start time specified" unless exists $opt{start};
  croak "No end time specified" unless exists $opt{end};
  croak "No time increment specified" unless exists $opt{inc};

  # Get the increment as an integer (DateTime::Duration or Time::Seconds)
  my $inc = $opt{inc};
  if (UNIVERSAL::can($inc, "seconds")) {
    $inc = $inc->seconds;
  }
  croak "Increment must be greater than zero" unless $inc > 0;

  # Determine date class to use for calculations
  my $dateclass = blessed( $opt{start} );
  croak "Start time must be either Time::Piece or DateTime object"
    if (!$dateclass || 
	($dateclass ne "Time::Piece" && $dateclass ne 'DateTime' ));

  my @data;

  # Get the current datetime
  my $original;
  my $usenow = $self->usenow;
  $original = $self->datetime() unless $usenow;

  # Get a private copy of the date object for calculations
  # (copy constructor)
  my $current = _clone_time( $opt{start} );

  while ( $current->epoch <= $opt{end}->epoch ) {

    # Hash for storing the data
    my %timestep;

    # store a copy of the time
    $timestep{time} = _clone_time( $current );

    # Set the time in the object
    # [standard problem with knowing whether we are overriding
    # another setting]
    $self->datetime( $current );

    # Now calculate the positions
    ($timestep{azimuth}, $timestep{elevation}) = $self->azel();
    $timestep{parang} = $self->pa;

    if (defined $opt{units} ) {
      $timestep{azimuth} = $timestep{azimuth}->in_format( $opt{units} );
      $timestep{elevation} = $timestep{elevation}->in_format( $opt{units} );
      $timestep{parang} = $timestep{parang}->in_format( $opt{units} );
    }

    $timestep{lst} = $self->_lst();

    # store the timestep
    push(@data, \%timestep);

    # increment the time
    $current = _inc_time( $current, $inc );

  }

  # Restore the time
  if ($usenow) {
    $self->datetime( undef );
  } else {
    $self->datetime( $original );
  }
  $self->usenow( $usenow );

  return @data;

}

sub rise_time {
  my $self = shift;
  return $self->_rise_set_time( 1, @_);
}

sub set_time {
  my $self = shift;
  return $self->_rise_set_time( 0, @_);
}

sub ha_set {
  my $self = shift;

  # Get the reference horizon elevation
  my %opt = @_;

  my $horizon = (defined $opt{horizon} ? $opt{horizon} :
		 $self->_default_horizon );

  # Get the telescope position
  my $tel = $self->telescope;

  # Get the longitude (in radians)
  my $lat = (defined $tel ? $tel->lat : 0.0 );

  # Declination
  my $dec = $self->dec_app;

  # Calculate the hour angle for this elevation
  # See http://www.faqs.org/faqs/astronomy/faq/part3/section-5.html
  my $cos_ha0 = ( sin($horizon) - sin($lat)*sin( $dec ) ) /
    ( cos($lat) * cos($dec) );

  # Make sure we have a valid number for the cosine
  if (defined $self->name && lc($self->name) eq 'moon' && abs($cos_ha0) > 1) {
    # for the moon this routine can incorrectly determine
    # cos HA near transit [in fact it always will be inaccurate
    # but near transit it won't return any value at all]
    # Calculate transit elevation and if it is greater than the
    # requested "horizon" use an iterative technique to find the
    # set time.
    if ($self->transit_el > $horizon) {
      my $oldtime = ( $self->has_datetime ? $self->datetime : undef);
      my $mt = $self->meridian_time();
      $self->datetime( $mt );
      print "# Calculating iterative set time for HA determination\n"
         if $DEBUG;
      my $convok = $self->_iterative_el( $horizon, -1 );
      return undef unless $convok;
      my $seconds = $self->datetime->epoch - $mt->epoch;
      $cos_ha0 = cos( $seconds * Astro::SLA::DS2R );
      $self->datetime( $oldtime );
    }
  }

  return undef if abs($cos_ha0) > 1;

  # Work out the hour angle for this elevation
  my $ha0 = acos( $cos_ha0 );

  # If we are the Sun we need to convert this to solar time
  # time from sidereal time
  $ha0 *= 365.2422/366.2422
    unless (defined $self->name && 
            lc($self->name) eq 'sun' && $self->isa("Astro::Coords::Planet"));


#  print "HA 0 is $ha0\n";
#  print "#### in hours: ". ( $ha0 * Astro::SLA::DR2S / 3600)."\n";

  # return the result (converting if necessary)
  return Astro::Coords::Angle::Hour->new( $ha0, units => 'rad',
					  range => 'PI')->in_format($opt{format});

}

sub meridian_time {
  my $self = shift;
  my %opt = @_;

  # extract the true value of "event" given all the backwards compatibility
  # problems
  my $event = $self->_extract_event( %opt );

  # Get the current time (do not modify it since we need to put it back)
  my $oldtime = ( $self->has_datetime ? $self->datetime : undef);

  # Need a reference time for the calculation
  my $reftime = $self->datetime;

  # For fast moving objects such as planets, we need to calculate
  # the transit time iteratively since the apparent RA/Dec will change
  # slightly during the night so we need to adjust the internal clock
  # to get it close to the actual transit time. We also need to make sure
  # that we are starting at the correct reference time so start at the
  # current time and look forward until we get a transit time > than
  # our start time

  # calculate the required times
  my $mtime;
  if ($event == 1 || $event == -1) {
    # previous or next
    $mtime = $self->_calc_mtime( $reftime, $event );

    # Reset the clock
    $self->datetime( $oldtime );

  } elsif ($event == 0) {
    # nearest requires both
    my $prev = $self->_calc_mtime( $reftime, -1 );
    my $next = $self->_calc_mtime( $reftime,  1 );

    # Reset the clock (before we possibly croak)
    $self->datetime( $oldtime );

    # calculate the diff
    if (defined $prev && defined $next) {
      my $prev_diff = abs( $reftime->epoch - $prev->epoch );
      my $next_diff = abs( $reftime->epoch - $next->epoch );

      if ($prev_diff < $next_diff) {
	$mtime = $prev;
      } else {
	$mtime = $next;
      }

    } elsif (defined $prev) {
      $mtime = $prev;
    } elsif (defined $next) {
      $mtime = $next;
    } else {
      croak "Should not occur in meridian_time\n";
    }
  } else {
    croak "Unrecognized value for event: $event\n";
  }

  return $mtime;
}

sub _calc_mtime {
  my $self = shift;
  my ($reftime, $event ) = @_;

  # event must be 1 or -1
  if (!defined $event || ($event != 1 && $event != -1)) {
    croak "Event must be either +1 or -1";
  }

  # do we have DateTime objects
  my $dtime = $self->_isdt();

  # Somewhere to store the previous time so we can make sure things
  # are iterating nicely
  my $prevtime;

  # The current best guess of the meridian time
  my $mtime;

  # Number of times we want to loop before aborting
  my $max = 10;

  # Tolerance for good convergence
  my $tol = 1;

  # Increment (in hours) to jump forward each loop
  # Need to make sure we lock onto the correct transit so I'm
  # wary of jumping forward by exactly 24 hours
  my $inc = 12 * $event;
  $inc /= 2 if (defined $self->name && lc($self->name) eq 'moon');

  # Loop until mtime is greater than the reftime
  # and (mtime - prevtime) is smaller than a second
  # and we have not looped more than $max times
  # There is probably an analytical solution. The problem is that
  # the apparent RA depends on the current time yet the apparent RA
  # varies with time
  my $count = 0;
  print "Looping..............".$reftime->datetime."\n" if $DEBUG;
  while ( $count <= $max ) {
    $count++;

    if (defined $mtime) {
      $prevtime = _clone_time( $mtime );
      $self->datetime( $mtime );
    }
    $mtime = $self->_local_mtcalc();
    print "New meridian time: ".$mtime->datetime ."\n" if $DEBUG;

    # make sure we are going the correct way
    # since we are enforced to only find a transit in the direction
    # governed by "event"

    # Calculate the difference in epoch seconds before the current
    # object reference time and the calculate transit time.
    # Use ->epoch rather than overload since I'm having problems
    # with Duration objects
    my $diff = $reftime->epoch - $mtime->epoch;
    $diff *= $event; # make the comparison correct sense
    if ($diff > 0) {
      print "Need to offset....[diff = $diff sec]\n" if $DEBUG;
      # this is an earlier transit time
      # Need to keep jumping forward until we lock on to a meridian
      # time that ismore recent than the ref time
      if ($dtime) {
	$mtime->add( hours => ($count * $inc));
      } else {
	$mtime = $mtime + ($count * $inc * Time::Seconds::ONE_HOUR);
      }
    }

    # End loop if the difference between meridian time and calculated
    # previous time is less than the acceptable tolerance
    if (defined $prevtime && defined $mtime) {
      last if (abs($mtime->epoch - $prevtime->epoch) <= $tol);
    }
  }

  # warn if we did not converge
  warnings::warnif "Meridian time calculation failed to converge"
     if $count > $max;

  # return the time
  return $mtime;
}

# Returns true if 
#    time - reftime is negative

# Returns RA-LST added on to reference time
sub _local_mtcalc {
  my $self = shift;

  # Now calculate the offset from the RA of the source.
  # Note that RA should be apparent RA and so the time should
  # match the actual time stored in the object.
  # Make sure the LST and Apparent RA are -PI to +PI
  # so that we do not jump whole days
  my $lst = Astro::SLA::slaDrange($self->_lst);
  my $ra_app = Astro::SLA::slaDrange( $self->ra_app );
  my $offset = $ra_app - $lst;

  # This is in radians. Need to convert it to seconds
  my $offset_sec = $offset * Astro::SLA::DR2S;

#  print "LST:            $lst\n";
#  print "RA App:         ". $self->ra_app ."\n";
#  print "Offset radians: $offset\n";
#  print "Offset seconds: $offset_sec\n";

  # If we are not the Sun we need to convert this to sidereal
  # time from solar time
  $offset_sec *= 365.2422/366.2422
    unless (defined $self->name &&
            lc($self->name) eq 'sun' && $self->isa("Astro::Coords::Planet"));

  my $datetime = $self->datetime;
  if ($self->_isdt()) {
    return $datetime->clone->add( seconds => $offset_sec );
  } else {
    return ($datetime + $offset_sec);
  }

#  return $mtime;
}

sub transit_el {
  my $self = shift;

  # Get meridian time
  my $mtime = $self->meridian_time();

  # Cache the current time if required
  # Note that we can leave $cache as undef if there is no
  # real time.
  my $cache;
  $cache = $self->datetime if $self->has_datetime;

  # set the new time
  $self->datetime( $mtime );

  # calculate the elevation
  my $el = $self->el( @_ );

  # fix the time back to what it was (including an undef value
  # if we did not read the cache).
  $self->datetime( $cache );

  return $el;
}

sub rv {
  my $self = shift;
  return (defined $self->{RadialVelocity} ? $self->{RadialVelocity} : 0 );
}

# internal set routine
sub _set_rv {
  my $self = shift;
  $self->{RadialVelocity} = shift;
}

sub redshift {
  my $self = shift;
  my $vd = $self->vdefn;
  if ($vd eq 'REDSHIFT' || $vd eq 'OPTICAL') {
    return ( $self->rv / CLIGHT );
  } elsif ($vd eq 'RELATIVISTIC') {
    # need to add
    return undef;
  } else {
    my $rv = $self->rv;
    # 1% of light speed
    if ( $rv > ( 0.01 * CLIGHT) ) {
      my $vopt = $rv / ( 1 - ( $rv / CLIGHT ) );
      return ( $vopt / CLIGHT );
    } else {
      return undef;
    }
  }
}

# internal set routine
sub _set_redshift {
  my $self = shift;
  my $z = shift;
  $z = 0 unless defined $z;
  $self->_set_rv( CLIGHT * $z );
  $self->_set_vdefn( 'REDSHIFT' );
  $self->_set_vframe( 'HEL' );
}

sub vdefn {
  my $self = shift;
  return (defined $self->{VelocityDefinition} ? $self->{VelocityDefinition} : 'OPTICAL' );
}

# internal set routine
sub _set_vdefn {
  my $self = shift;
  my $defn = shift;
  # undef resets to default
  if (defined $defn) {
    $defn = $self->_normalise_vdefn( $defn );
  }
  $self->{VelocityDefinition} = $defn;
}

sub vframe {
  my $self = shift;
  return (defined $self->{VelocityFrame} ? $self->{VelocityFrame} : 'HEL' );
}

# internal set routine
sub _set_vframe {
  my $self = shift;
  my $frame = shift;
  if (defined $frame) {
    # undef resets to default
    $frame = $self->_normalise_vframe( $frame );
  }
  $self->{VelocityFrame} = $frame;
}

sub obsvel {
  my $self = shift;
  my $vdefn = $self->vdefn;
  my $vframe = $self->vframe;
  my $rv = $self->rv;

  # Now we need to calculate the observer velocity in the
  # target frame
  my $vobs = $self->vdiff( '', 'TOP' );

  # Total velocity between observer and target
  my $vtotal = $vobs + $rv;

  return $vtotal;
}

sub doppler {
  my $self = shift;
  my $vdefn = $self->vdefn;
  my $obsvel = $self->obsvel;

  # Doppler correction depends on definition
  my $doppler;
  if ( $vdefn eq 'RADIO' ) {
    $doppler = 1 - ( $obsvel / CLIGHT );
  } elsif ( $vdefn eq 'OPTICAL' || $vdefn eq 'REDSHIFT' ) {
    if ( $obsvel > (0.01 * CLIGHT)) {
      # Relativistic velocity
      # First calculate the redshift correction
      my $zcorr = 1 / ( 1 + $self->redshift );

      # Now the observer doppler correction to Heliocentric frame
      my $vhel = $self->vhelio;
      my $obscorr = 1 - ( $vhel / ( CLIGHT * $vhel) );

      $doppler = $zcorr * $obscorr;

    } else {
      # small radial velocity, use standard doppler formula
      $doppler = 1 - ( $obsvel / ( CLIGHT + $obsvel ) );
    }
  } elsif ( $vdefn eq 'RELATIVISTIC' ) {
    # do we need to use the same correction as for OPTICAL and REDSHIFT?
    # presumably
    $doppler = sqrt( ( CLIGHT - $obsvel ) / ( CLIGHT + $obsvel ) );
  } else {
    croak "Can not calculate doppler correction for unsupported definition $vdefn\n";
  }
  return $doppler;
}

sub vdiff {
  my $self = shift;
  my $f1 = ( shift || $self->vframe );
  my $f2 = ( shift || $self->vframe );

  # convert the arguments to standardised frames
  $f1 = $self->_normalise_vframe( $f1 );
  $f2 = $self->_normalise_vframe( $f2 );

  return 0 if $f1 eq $f2;

  # put all the supported answers in a hash relative to TOP
  my %vel;
  $vel{TOP} = 0;
  $vel{GEO} = $self->verot();
  $vel{HEL} = $self->vhelio;
  $vel{BAR} = $self->vbary;
  $vel{LSRK} = $self->vlsrk;
  $vel{LSRD} = $self->vlsrd;
  $vel{GAL}  = $self->vgalc;
  $vel{LG}   = $self->vlg;

  # now the difference is easy
  return ( $vel{$f1} - $vel{$f2} );
}

sub verot {
  my $self = shift;

  # Local Sidereal Time
  my $lst = $self->_lst;

  # Observer location
  my $tel = $self->telescope;
  my $lat = (defined $tel ? $tel->lat : 0 );

  # apparent ra dec
  my ($ra, $dec) = $self->apparent();

  return Astro::SLA::slaRverot( $lat, $ra, $dec, $lst );
}

sub vorb {
  my $self = shift;
  my $usebary = shift;

  # Earth velocity (and position)
  my @vb = (0,0,0);
  my @pb = (0,0,0);
  my @vh = (0,0,0);
  my @ph = (0,0,0);
  Astro::SLA::slaEvp($self->_mjd_tt(), 2000.0,@vb,@pb,@vh,@ph);

  # Barycentric or heliocentric velocity?
  my @v = ( $usebary ? @vb : @vh );

  # Convert spherical source coords to cartesian
  my ($ra, $dec) = $self->radec;
  my @cart = (0,0,0);
  Astro::SLA::slaDcs2c($ra,$dec,@cart);

  # Velocity due to Earth's orbit is scalar product of the star position
  # with the Earth's velocity
  my $vorb = - Astro::SLA::slaDvdv(@cart,@v)* AU2KM;
  return $vorb;
}

sub vhelio {
  my $self = shift;
  return ($self->verot + $self->vorb);
}

sub vbary {
  my $self = shift;
  return ($self->verot + $self->vorb(1));
}

sub vlsrk {
  my $self = shift;
  my ($ra, $dec) = $self->radec;
  return (Astro::SLA::slaRvlsrk( $ra, $dec ) + $self->vhelio);
}

sub vlsrd {
  my $self = shift;
  my ($ra, $dec) = $self->radec;
  return (Astro::SLA::slaRvlsrd( $ra, $dec ) + $self->vhelio);
}

sub vgalc {
  my $self = shift;
  my ($ra, $dec) = $self->radec;
  return (Astro::SLA::slaRvgalc( $ra, $dec ) + $self->vlsrd);
}

sub vlg {
  my $self = shift;
  my ($ra, $dec) = $self->radec;
  return (Astro::SLA::slaRvlg( $ra, $dec ) + $self->vhelio);
}

sub _lst {
  my $self = shift;
  my $time = $self->datetime;
  my $tel = $self->telescope;

  # Get the longitude (in radians)
  my $long = (defined $tel ? $tel->long : 0.0 );

  # Seconds can be floating point.
  my $sec = $time->sec;
  if ($time->can("nanosecond")) {
    $sec += 1E-9 * $time->nanosecond;
  }

  # Return the first arg
  # Note that we guarantee a UT time representation
  my $lst = (Astro::SLA::ut2lst( $time->year, $time->mon,
				 $time->mday, $time->hour,
				 $time->min, $sec, $long))[0];

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

sub _mjd_tt {
  my $self = shift;
  my $mjd = $self->datetime->mjd;
  my $offset = Astro::SLA::slaDtt( $mjd );
  $mjd += ($offset / (86_400));
  return $mjd;
}

sub _clone_time {
  my $input = shift;
  return unless defined $input;

  if (UNIVERSAL::isa($input, "Time::Piece")) {
    return Time::Piece::gmtime( $input->epoch );
  } elsif (UNIVERSAL::isa($input, "DateTime")) {
    # Use proper clone method
    return $input->clone;
  }
  return;
}

sub _inc_time {
  my $date = shift;
  my $delta = shift;

  # convert to integer seconds
  $delta = $delta->seconds() if UNIVERSAL::can( $delta, "seconds" );

  # Check for Time::Piece, else assume DateTime
  if (UNIVERSAL::isa( $date, "Time::Piece" ) ) {
    # can not add time in place
    $date += $delta;
  } else {
    # increment the time (this happens in place)
    if (abs($delta) > 1) {
      $date->add( seconds => "$delta" );
    } else {
      $date->add( nanoseconds => ( $delta * 1E9 ) );
    }
  }

  return $date;
}

sub _cmp_time {
  my $t1 = shift;
  my $t2 = shift;
  my $e1 = (defined $t1 ? $t1->epoch : 0);
  my $e2 = (defined $t2 ? $t2->epoch : 0);
  return $e1 <=> $e2;
}

sub _default_horizon {
  return 0;
}

sub _rise_set_time {
  my $self = shift;
  my $rise = shift;
  my %opt = @_;

  # Calculate the HA required for setting
  my $ha_set = $self->ha_set( %opt, format=> 'radians' );
  return if ! defined $ha_set;

  # extract the event information
  my $event = $self->_extract_event( %opt );

  # and convert to seconds
  $ha_set *= Astro::SLA::DR2S;

  # and thence to a duration if required
  if ($self->_isdt()) {
    $ha_set = new DateTime::Duration( seconds => $ha_set );
  }

  # Get the current time (do not modify it since we need to put it back)
  my $reftime = $self->datetime;

  # Determine whether we have to remember the cache
  my $havetime = $self->has_datetime;

  # Need the requested horizon
  my $refel = (defined $opt{horizon} ? $opt{horizon} :
		 $self->_default_horizon );

  # somewhere to store the times
  my @times;

  # Calculate the set or rise time for both the previous and next meridian time
  # so that we can work out which set time to return. There is probably
  # an analytic short cut that can be used to decide whether the
  # correct set time will be related to the previous or next meridian
  # time simply by calculating the meridian time once and shifting
  # it by 24 hours

  for my $n (-1, 1) {
    # Calculate the transit time
    print "Calculating " .($n == -1 ? "previous" : "next") . " meridian time\n"
      if $DEBUG;
    my $mt = $self->meridian_time( event => $n );

    # First guestimate for the event
    my $event_time;
    if ($rise) {
      $event_time = $mt - $ha_set;
    } else {
      $event_time = $mt + $ha_set;
    }

    # Now converge about this value
    $self->datetime( $event_time );
    my $convok = $self->_iterative_el( $refel, ($rise ? 1 : -1) );
    $event_time = ( $convok ? $self->datetime : undef );
    print "**** Failed to converge\n" if ($DEBUG && !$convok);

    # and restore the reference date to reset the clock
    $self->datetime( ( $havetime ? $reftime : undef ) );

    # store the event time
    push(@times, $event_time) if defined $event_time;

    print "Determined ".($rise ? "rise" : "set" ) . " time of ".
       (defined $event_time ? $event_time : 'undef')."\n"
       if $DEBUG;
  }

  # choose a value
  my $event_time;
  if ($event == -1) {
    # We want the nearest time that is earlier than reference time
    # Start from the end and jump out when 
    for my $t (reverse @times) {
      if ($t->epoch <= $reftime->epoch) {
	$event_time = $t;
	last;
      }
    }
  } elsif ($event == 1) {
    # start from the oldest until we hit a new time
    for my $t (@times) {
      if ($t->epoch >= $reftime->epoch) {
	$event_time = $t;
	last;
      }
    }
  } elsif ($event == 0) {
    # calculate the diff
    my %diffs = map { abs( $reftime->epoch - $_->epoch), $_ } @times;
    my @sort = sort { $a <=> $b } keys %diffs;
    $event_time = $diffs{$sort[0]};

  } else {
    croak "Unrecognized event specifier in set_time";
  }


  return $event_time;
}

sub _iterative_el {
  my $self = shift;
  my $refel = shift;
  my $grad = shift;

  # See what type of date object we are dealing with
  my $use_dt = $self->_isdt();

  # We are tweaking DateTime objects without the cache knowing about
  # it (because it is faster than lots of object cloning) so we must
  # turn off caching
  my $old_unsafe = $self->datetime_is_unsafe;
  $self->datetime_is_unsafe( 1 );

  # Calculate current elevation
  my $el = $self->el;

  # Tolerance (in arcseconds)
  my $tol = 5 * Astro::SLA::DAS2R;

  # smallest support increment
  my $smallinc = 0.2;

  # Get the estimated time for this elevation
  my $time = $self->datetime;

  # now compare the requested elevation with the actual elevation for the
  # previously calculated rise time
  if (abs($el - $refel) > $tol ) {
    if ($DEBUG) {
      print "# ================================ -> ".$self->name."\n";
      print "# Requested elevation: " . (Astro::SLA::DR2D * $refel) ."\n";
      print "# Elevation out of range: ". $self->el(format => 'deg')."\n";
      print "# For " . ($grad > 0 ? "rise" : "set")." time: ". $time->datetime ."\n";
    }

    # use 1 minute for all except the moon
    my $inc = 60; # seconds
    $inc *= 10 if (defined $self->name && lc($self->name) eq 'moon');

    my $sign = ($el < $refel ? 1 : -1); # incrementing or decrementing time
    my $prevel; # previous elevation

    # Indicate whether we have straddled the correct answer
    my $has_been_high;
    my $has_been_low;

    # This is a very simple convergence algorithm.
    # Newton-Raphson would be much faster given that the function
    # is almost linear for most elevations.
    while (abs($el-$refel) > $tol) {
      if (defined $prevel) {

	# should check sign of gradient to make sure we are not
	# running away to an incorrect gradient. Note that this function
	# can be highly curved if we are near ante transit.

	# There are two constraints that have to be used to control
	# the convergence
	#  - make sure we are not diverging with each iteration
	#  - make sure we are not bouncing around missing the final
	#    value due to inaccuracies in the time resolution (this is
	#    important for fast moving objects like the sun). Currently
	#    we are not really doing anything useful by tweaking by less
	#    than a second

	# The problem is that we can not stop on bracketing the correct
	# value if we are never going to reach the value itself because
	# it never rises or never sets

	#          \   /
	#           \ /
	#            -   <--  minimum elevation
	#                <--  required elevation

	# in the above case each step will diverge but will never
	# straddle

	# The solution is to first iterate by using a divergence
	# criterion. Then, if we determine that we have straddled the
	# correct result, we switch to a stopping criterion that
	# uses the time resolution if we never fully converge (the
	# important thing is to return a time not an elevation)

	my $diff_to_prev = $prevel - $refel;
	my $diff_to_curr = $el - $refel;

	# set the straddle parameters
	if ($diff_to_curr < 0) {
	  $has_been_low = 1;
	} else {
	  $has_been_high = 1;
	}

	# if we have straddled, stop on increment being small

	# now determine whether we should reverse
	# if we know we are bouncing around the correct value
	# we can reverse as soon as the previous el and the current el
	# are either side of the correct answer
	my $reverse;
	if ($has_been_low && $has_been_high) {
	  # we can abort if we are below the time resolution
	  last if $inc < $smallinc;

	  # reverse 
	  $reverse = 1 if ($diff_to_prev / $diff_to_curr < 0);
	} else {
	  # abort if the inc is too small
	  if ($inc < $smallinc) {
	    $self->datetime_is_unsafe( $old_unsafe );
	    return undef;
	  }

	  # if we have not straddled, we reverse if the previous el
	  # is closer than this el
	  $reverse = 1 if abs($diff_to_prev) < abs($diff_to_curr);
	}

	# if the increment is small, bounce


	# reverse
	if ( $reverse ) {

	  # the gap between the previous measurement and the reference
 	  # is smaller than the current gap. We seem to be diverging.
	  # Change direction
	  $sign *= -1;
	  # and use half the step size
	  $inc /= 2;

	  # in the linear approximation
	  # we know the gradient

	}
      }

      # Now calculate a new time
      my $delta = $sign * $inc;
      if (!$use_dt) {
	$time = $time + $delta;
	# we have created a new object so need to store it for next time
	# round
	$self->datetime( $time );
      } else {
	# increment the time (this happens in place so we do not need to
	# register the change with the datetime method
	if (abs($delta) > 1) {
	  $time->add( seconds => "$delta" );
	} else {
	  $time->add( nanoseconds => ( $delta * 1E9 ) );
	}

      }
      # recalculate the elevation, storing the previous as reference
      $prevel = $el;
      $el = $self->el;
      print "# New elevation: ". $self->el(format=>'deg')." \t@ ".$time->datetime." with delta $delta sec\n"
	if $DEBUG;
    }

  }

  # reset the unsafeness level
  $self->datetime_is_unsafe( $old_unsafe );

  return 1;
}

sub _isdt {
  my $self = shift;
  my $test = shift;

  $test = $self->datetime unless defined $test;

  # Check Time::Piece first since there is a possibility that 
  # this is really a subclass of DateTime
  my $dtime;
  if ($test->isa( "Time::Piece")) {
    $dtime = 0;
  } elsif ($test->isa("DateTime")) {
    $dtime = 1;
  } else {
    croak "Unknown DateTime object class ". blessed($test);
  }

  return $dtime;
}

sub _extract_event {
  my $self = shift;

  # make sure the $opt{event} always exists
  my $default = 1;
  my %opt = (event => $default, @_);

  if (exists $opt{nearest}) {
    if ($opt{nearest}) {
      # request for nearest
      return 0;
    } else {
      # request for not nearest
      if ($opt{event} == 0) {
	# this is incompatible with nearest=0 so return default
	return $default;
      } else {
	return $opt{event};
      }
    }

  } else {
    return $opt{event};
  }
}

sub _normalise_vframe {
  my $self = shift;
  my $in = shift;

  croak "Velocity frame not defined. Can not normalise" unless defined $in;

  # upper case
  $in = uc( $in );

  # LSRK or LSRD need no normalisation
  return $in if ($in eq 'LSRK' || $in eq 'LSRD' || $in eq 'LG');

  # Truncate
  my $trunc = substr( $in, 0, 3 );

  # Verify
  croak "Unrecognized velocity frame '$trunc'"
    unless $trunc =~ /^(GEO|TOP|HEL|LSR|GAL|BAR)/;

  # special case
  $trunc = 'LSRK' if $trunc eq 'LSR';

  # okay
  return $trunc;
}

sub _normalise_vdefn {
  my $self = shift;
  my $in = shift;

  croak "Velocity definition not defined. Can not normalise" unless defined $in;

  # upper case
  $in = uc( $in );

  # Truncate
  my $trunc = substr( $in, 0, 3 );

  # Verify
  if ($trunc eq 'RAD') {
    return 'RADIO';
  } elsif ($trunc eq 'OPT') {
    return 'OPTICAL';
  } elsif ($trunc eq 'RED') {
    return 'REDSHIFT';
  } elsif ($trunc eq 'REL') {
    return 'RELATIVISTIC';
  } else {
    croak "Unrecognized velocity definition '$trunc'";
  }
}

sub _parse_equinox {
  my $self = shift;
  my $str = shift;
  my ($sys, $epoch) = ('FK5', 2000.0);
  if ($str =~ /^([BJ]?)(\d+(\.\d+)?)$/i) {
    my $typ = $1;
    $sys = ($typ eq 'B' ? 'FK4' : 'FK5' );
    $epoch = $2;
  } else {
    warnings::warnif( "Supplied equinox '$str' does not look like an equinox");
  }
  return ($sys, $epoch);
}

sub _j2000_to_byyyy {
  my $self = shift;
  my ($equ, $ra2000, $dec2000) = @_;

  # First to 1950
  Astro::SLA::slaFk54z($ra2000, $dec2000, 
		       Astro::SLA::slaEpb( $self->_mjd_tt ),
		       my $rb, my $db, my $drb, my $drd);

  # Then preces to reference epoch frame
  # I do not know whether fictitious proper motions should be included
  # here with slaPm or whether it is enough to use non-1950 epoch
  # in slaFk54z and then preces 1950 to the same epoch. Not enough test
  # data for this rare case.
  if ($equ != 1950) {
    # Add E-terms
    Astro::SLA::slaSubet( $rb, $db, 1950.0, my $rnoe, my $dnoe);

    # preces
    Astro::SLA::slaPreces( 'FK4', 1950, $equ, $rnoe, $dnoe);

    # Add E-terms
    Astro::SLA::slaAddet( $rnoe, $dnoe, $equ, $rb, $db);

  }
  return ($rb, $db);
}

sub _cache_write {
  my $self = shift;
  my %add = @_;

  my $primary = $self->_cache_key;
  return () unless defined $primary;

  my $C = $self->_cache_ref;

  if (!exists $C->{$primary}) {
    $C->{$primary} = {};
  }

  my $local = $C->{$primary};

  for my $key (keys %add) {
    $local->{$key} = $add{$key};
  }
  return;
}

sub _cache_read {
  my $self = shift;
  my @keys = @_;

  my $primary = $self->_cache_key;
  return () unless defined $primary;

  my $C = $self->_cache_ref;

  return () unless exists $C->{$primary};
  my $local = $C->{$primary};
  return () unless defined $local;

  my @answer =  map { $local->{$_} } @keys;
#  print "Caller: ". join(":", caller() ) ."\n";
#  print "Getting cache values for ". join(":",@keys) ."\n";
#  print "Getting cache values for ". join(":",@answer) ."\n";
  return (wantarray() ? @answer : $answer[0] );
}

sub _calc_cache_key {
  my $self = shift;
#  return; # This will disable caching

  # if we have a floating clock we do not want to cache
  if (!$self->has_datetime || $self->usenow || $self->datetime_is_unsafe) {
    $self->_cache_key( undef );
    return;
  }

  # The cache key currently uses the time and the telescope name
  my $dt = $self->datetime;
  my $tel = $self->telescope;
  my $telName = (defined $tel ? $tel->name : "NONE" );

  print "# Calculating cache key using $telName and ". $dt->epoch ."\n"
    if $DEBUG;

  # usually epoch is quicker to generate than ISO date but we have to
  # be careful about fractional seconds in DateTime (Time::Piece can
  # not do it)

  my $addendum = "";
  $addendum = $dt->nanosecond if $dt->can("nanosecond");

  # Use date + telescope name as key
  $self->_cache_key($telName . "_" . $dt->epoch. $addendum); 
}

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

{
  my $KEY = "__AC_CACHE";
  sub _cache_ref {
    my $self = shift;
    if (!exists $self->{$KEY} ) {
      $self->{$KEY} = {};
    }
    return $self->{$KEY};
  }
}