view release on metacpan or  search on metacpan

lib/Astro/Coord/ECI/VSOP87D.pm  view on Meta::CPAN

package Astro::Coord::ECI::VSOP87D;

use 5.008;

use strict;
use warnings;

use utf8;

use Astro::Coord::ECI::Utils qw{
    AU PI SECSPERDAY TWOPI
    asin deg2rad jcent2000 julianday
    load_module looks_like_number mod2pi
    rad2deg rad2dms rad2hms tan
};
use Exporter qw{ import };
use Carp;
use POSIX qw{ floor };
use Storable qw{ dclone };

use constant DAYS_PER_JULIAN_MILENNIUM	=> 365250;
use constant CODE_REF	=> ref sub {};
use constant HASH_REF	=> ref {};
use constant SUN_CLASS	=> __PACKAGE__ . '::Sun';

BEGIN {
    __PACKAGE__->can( 'DEBUG' )
	or constant->import( DEBUG => 0 );
}

our $VERSION = '0.008';

my @basic_export = qw{
    SUN_CLASS
    model_cutoff_definition
    nutation obliquity order period time_set year
    __default
    __get_attr
    __mutate_model_cutoff __mutate_nutation_cutoff
};

our @EXPORT_OK = (
    @basic_export,
    qw{
	geometric_longitude
	synodic_period
	__angle_subtended_from_earth
	__longitude_from_sun
	__model
    },
);
our %EXPORT_TAGS = (
    mixin	=> [ @basic_export, qw{
	synodic_period
	__angle_subtended_from_earth
	__longitude_from_sun
	__model
    } ],
    sun		=> [ @basic_export, qw{ geometric_longitude } ],
);

# We want to ensure SUN_CLASS is loaded because it is our default Sun
# implementation, but we do it at run time to avoid a circular
# dependency.
load_module( SUN_CLASS );

sub time_set {
    my ( $self ) = @_;
    my $attr = $self->__get_attr();
    my $time = $self->dynamical();
    my $cutoff = $self->get( 'model_cutoff' );
    my $cutoff_def = $self->model_cutoff_definition();
    my $sun = $self->get( 'sun' );

lib/Astro/Coord/ECI/VSOP87D.pm  view on Meta::CPAN

	}

	# The model computes nutations in milli arc seconds, but we
	# return them in radians
	return ( @memoize_result =
	    map { deg2rad( $_ / 3600_0000 ) } $delta_psi, $delta_eps );
    }

}

sub __get_attr {
    my ( $self ) = @_;
    my $ref = ref $self
	or confess 'Can not call as static method';
    return $self->{$ref} ||= {
	model_cutoff_definition	=> dclone( $self->__model_definition(
		'default_model_cutoff' ) ),
    };
}

sub __default {
    my ( $class, $arg ) = @_;

    my $name = $class->__model_definition( 'body' );
    defined $arg->{id}
	or $arg->{id} = $name;
    defined $arg->{name}
	or $arg->{name} = $name;

    defined $arg->{diameter}
	or $arg->{diameter} = $class->__model_definition( 'diameter' );

    defined $arg->{model_cutoff}
	or $arg->{model_cutoff} = 'Meeus';

    defined $arg->{nutation_cutoff}
	or $arg->{nutation_cutoff} = 3;

    return;
}

sub __mutate_model_cutoff {
    my ( $self, $name, $val ) = @_;
    defined $val
	or croak "model cutoff must be defined";
    $self->model_cutoff_definition( $val )
	or croak "model cutoff '$val' is unknown";
    $self->__get_attr()->{$name} = $val;
    return $self;
}

sub model_cutoff_definition {
    my ( $self, $name, @arg ) = @_;
    defined $name
	or $name = $self->get( 'model_cutoff' );
    my $attr = $self->__get_attr();
    if ( @arg ) {
	if ( defined( my $val = $arg[0] ) ) {
	    unless ( ref $val ) {
		looks_like_number( $val )
		    and $val !~ m/ \A Inf (?: inity )? | NaN \z /smx
		    or croak 'Scalar model cutoff definition must be a number';
		my $num = $val;
		$val = sub {
		    my ( @model ) = @_;
		    my %cutoff;
		    foreach my $series ( @model ) {
			my $count = 0;
			foreach my $term ( @{ $series->{terms} } ) {
			    last if $term->[0] < $num;
			    $count++;
			}
			$count
			    and $cutoff{$series->{series}} = $count;
		    }
		    return \%cutoff;
		};
	    }
	    if ( CODE_REF eq ref $val ) {
		$val = $val->(
		    map { @{ $_ } } @{
			$self->__model_definition( 'model' ) }
		);
		$val->{name} = $name;
	    }
	    HASH_REF eq ref $val
		or croak 'The model cutoff definition value must be a hash ref';
	    my $terms = $self->__model_definition(
		'default_model_cutoff')->{none};
	    foreach my $name ( keys %{ $val } ) {
		'name' eq $name
		    and next;
		exists $terms->{$name}
		    or croak "Series '$name' not in this model";
		$val->{$name} > $terms->{$name}
		    and croak "Series '$name' has only $terms->{$name} terms";
	    }
	    $attr->{model_cutoff_definition}{$name} = $val;
	} else {
	    $self->__model_definition( 'default_model_cutoff' )->{$name}
		and croak "You may not delete model cutoff definition '$name'";
	    delete $attr->{model_cutoff_definition}{$name};
	}
	return $self;
    } else {
	return $attr->{model_cutoff_definition}{$name};
    }
}

# Static method
# Given dynamical time in seconds, return Ecliptic position in Ecliptic
# longitude (radians), Ecliptic latitude (radians), Range (AU) and 
# AU and velocity in AU/day
sub __model {
    my ( $self, $time, %arg ) = @_;

    DEBUG
	and printf <<'EOD',

__model:
    invocant: %s
    model_cutoff: %s
EOD
    $self,
    ( $arg{model_cutoff_definition} ?
	( $arg{model_cutoff_definition}{name} || '<anonymous>' ) :
	'<unspecified>' ),
    ;

    my $jm = jcent2000( $time ) / 10;	# Meeus 32.1
    my @p_vec;
    my @v_vec;
    foreach my $coord ( @{ $self->__model_definition( 'model' ) } ) {
	my $dT = my $exponent = 0;
	my $T = 1;
	my $pos = my $vel = 0;
	foreach my $series ( @{ $coord } ) {
	    my $limit = $arg{model_cutoff_definition} ?
		( $arg{model_cutoff_definition}{$series->{series}} || 0 ) :
		@{ $series->{terms} }
		or next;
	    --$limit;
	    foreach my $inx ( 0 .. $limit ) {
		my $term = $series->{terms}[$inx];
		my $u = $term->[1] + $term->[2] * $jm;
		my $cos_u = cos $u;
		$pos += $term->[0] * $cos_u * $T;
		my $sin_u = sin $u;
		$vel += $dT * $exponent * $term->[0] * $cos_u -
		    $T * $term->[0] * $term->[2] * $sin_u;
	    }
	    $dT = $T;
	    $exponent++;
	    $T *= $jm;
	}
	$vel /= DAYS_PER_JULIAN_MILENNIUM;	# units/millennium -> units/day
	push @p_vec, $pos;
	push @v_vec, $vel;
    }
    $p_vec[0] = mod2pi( $p_vec[0] );
    return ( @p_vec, @v_vec );
}

sub __mutate_nutation_cutoff {
    my ( $self, undef, $val ) = @_;
    defined $val
	or croak 'Nutation cutoff must be defined';
    looks_like_number( $val )
	and $val >= 0
	or croak 'Nutation cutoff must be a non-negative number';
    $self->__get_attr()->{nutation_cutoff} = $val;
    return $self;
}

sub obliquity {
    my ( $self, $time ) = @_;

    defined $time
	or $time = $self->dynamical();
    # Obliquity per Meeus 22.3
    my $U = jcent2000( $time ) / 100;
    my $epsilon_0 = 0;
    $epsilon_0 = $epsilon_0 * $U + $_ for qw{ 2.45 5.79 27.87 7.12
	-39.05 -249.67 -51.38 1999.25 -1.55 -4680.93 84381.448 };
    $epsilon_0 = deg2rad( $epsilon_0 / 3600 );
    my ( undef, $delta_eps ) = $self->nutation( $time );
    my $epsilon = $epsilon_0 + $delta_eps;

    DEBUG
	and printf <<"EOD",

Obliquity:
𝛆₀ = %.7f
𝛆 = %.7f
  = %s
EOD
	rad2deg( $epsilon_0 ), rad2deg( $epsilon ), rad2dms( $epsilon );

    return $epsilon;
}

sub order {
    my ( $self ) = @_;
    return $self->__model_definition( 'order' );
}

# Calculate the period of the orbit. The orbital velocity in radians per
# Julian centutu is just the inverse of the coefficient of the constant
# L1 term -- that is, the A value of the one with B and C both zero,
# since that computes as A * cos( B + C * tau ) = A * cos( 0 ) = A
#
# So the tropical year is just TWOPI divided by this. But we want the
# Siderial year, so we have to add in the precession. I found no
# reference for this, but to a first approximation the additional angle
# to travel is the solar year times the rate of precession, and the
# siderial year is the solar year plus the additional distance divided
# by the orbital velocity.
#
# But this calculation is done when the model parameters are built, and
# period() just returns the resultant value.
#
# The rate of pecession is the IAU 2003 value, taken from fapa03.c.

sub period {
    my ( $self ) = @_;
    return $self->__model_definition( 'sidereal_period' );
}