Astro-satpass
view release on metacpan or search on metacpan
lib/Astro/Coord/ECI/Sun.pm view on Meta::CPAN
use Astro::Coord::ECI::Utils qw{ @CARP_NOT :mainstream };
use Carp;
use POSIX qw{ ceil };
use constant MEAN_MAGNITUDE => -26.8;
my %attrib = (
iterate_for_quarters => 1,
);
my %static = (
id => 'Sun',
name => 'Sun',
diameter => 1392000,
iterate_for_quarters => undef,
);
my $weaken = eval {
require Scalar::Util;
Scalar::Util->can('weaken');
};
our $Singleton = $weaken;
my %object; # By class
=item $sun = Astro::Coord::ECI::Sun->new();
This method instantiates an object to represent the coordinates of the
Sun. This is a subclass of L<Astro::Coord::ECI|Astro::Coord::ECI>, with
the id and name attributes set to 'Sun', and the diameter attribute set
to 1392000 km per Jean Meeus' "Astronomical Algorithms", 2nd Edition,
Appendix I, page 407.
Any arguments are passed to the set() method once the object has been
instantiated. Yes, you can override the "hard-wired" id, name, and so
forth in this way.
If C<$Astro::Coord::ECI::Sun::Singleton> is true, you get a singleton
object; that is, only one object is instantiated and subsequent calls to
C<new()> just return that object. If higher-accuracy subclasses are ever
implemented, there will be one singleton for each class.
The singleton logic only works if L<Scalar::Util|Scalar::Util> exports
C<weaken()>. If it does not, the setting of
C<$Astro::Coord::ECI::Sun::Singleton> is silently ignored. The default
is true if L<Scalar::Util|Scalar::Util> can be loaded and exports
C<weaken()>, and false otherwise.
=cut
sub new {
my ($class, @args) = @_;
ref $class and $class = ref $class;
if ( $Singleton && $weaken && __classisa( $class, __PACKAGE__ ) ) {
my $self;
if ( $self = $object{$class} ) {
$self->set( @args ) if @args;
return $self;
} else {
$self = $object{$class} = $class->SUPER::new (%static, @args);
$weaken->( $object{$class} );
return $self;
}
} else {
return $class->SUPER::new (%static, @args);
}
}
=item @almanac = $sun->almanac( $station, $start, $end );
This method produces almanac data for the Sun for the given observing
station, between the given start and end times. The station is assumed
to be Earth-Fixed - that is, you can't do this for something in orbit.
The C<$station> argument may be omitted if the C<station> attribute has
been set. That is, this method can also be called as
@almanac = $sun->almanac( $start, $end )
The start time defaults to the current time setting of the $sun
object, and the end time defaults to a day after the start time.
The almanac data consists of a list of list references. Each list
reference points to a list containing the following elements:
[0] => time
[1] => event (string)
[2] => detail (integer)
[3] => description (string)
The @almanac list is returned sorted by time.
The following events, details, and descriptions are at least
potentially returned:
horizon: 0 = Sunset, 1 = Sunrise;
transit: 0 = local midnight, 1 = local noon;
twilight: 0 = end twilight, 1 = begin twilight;
quarter: 0 = spring equinox, 1 = summer solstice,
2 = fall equinox, 3 = winter solstice.
Twilight is calculated based on the current value of the twilight
attribute of the $sun object. This attribute is inherited from
L<Astro::Coord::ECI|Astro::Coord::ECI>, and documented there.
=cut
sub __almanac_event_type_iterator {
my ( $self, $station ) = @_;
my $inx = 0;
my $horizon = $station->__get_almanac_horizon();
my @events = (
[ $station, next_elevation => [ $self, $horizon, 1 ],
horizon => '__horizon_name' ],
[ $station, next_meridian => [ $self ],
transit => '__transit_name' ],
[ $station, next_elevation =>
[ $self, $self->get( 'twilight' ) + $horizon, 0 ],
twilight => '__twilight_name' ],
[ $self, next_quarter => [], quarter => '__quarter_name', ],
);
return sub {
$inx < @events
and return @{ $events[$inx++] };
return;
};
}
use Astro::Coord::ECI::Mixin qw{ almanac };
=item @almanac = $sun->almanac_hash( $station, $start, $end );
This convenience method wraps $sun->almanac(), but returns a list of
hash references, sort of like Astro::Coord::ECI::TLE->pass()
does. The hashes contain the following keys:
{almanac} => {
{event} => the event type;
{detail} => the event detail (typically 0 or 1);
{description} => the event description;
}
{body} => the original object ($sun);
{station} => the observing station;
{time} => the time the quarter occurred.
The {time}, {event}, {detail}, and {description} keys correspond to
elements 0 through 3 of the list returned by almanac().
=cut
use Astro::Coord::ECI::Mixin qw{ almanac_hash };
=item $coord2 = $coord->clone ();
If singleton objects are enabled, this override of the superclass'
method simply returns the invocant. Otherwise it does a deep clone of an
object, producing a different but identical object.
Prior to version 0.099_01 it always returned a clone. Yes,
this is a change in long-standing functionality, but a long-standing bug
is still a bug.
=cut
sub clone {
my ( $self ) = @_;
$Singleton
and $weaken
and return $self;
return $self->SUPER::clone();
}
=item $elevation = $tle->correct_for_refraction( $elevation )
This override of the superclass' method simply returns the elevation
passed to it. I have no algorithm for refraction at the surface of the
photosphere or anywhere else in the environs of the Sun, and explaining
why I make no correction at all seemed easier than explaining why I make
an incorrect correction.
See the L<Astro::Coord::ECI|Astro::Coord::ECI> C<azel()> and
C<azel_offset()> documentation for whether this class'
C<correct_for_refraction()> method is actually called by those methods.
=cut
sub correct_for_refraction {
my ( undef, $elevation ) = @_; # Invocant unused
return $elevation;
}
=item $long = $sun->geometric_longitude ()
This method returns the geometric longitude of the Sun in radians at
the last time set.
=cut
sub geometric_longitude {
my $self = shift;
croak <<eod unless defined $self->{_sun_geometric_longitude};
Error - You must set the time of the Sun object before the geometric
longitude can be returned.
eod
return $self->{_sun_geometric_longitude};
}
=item $sun->get( ... )
This method has been overridden to return the invocant as the C<'sun'>
attribute.
=cut
sub get {
my ( $self, @args ) = @_;
my @rslt;
foreach my $name ( @args ) {
push @rslt, 'sun' eq $name ? $self :
$attrib{$name} ?
ref $self ? $self->{$name} : $static{$name} :
$self->SUPER::get( $name );
}
return wantarray ? @rslt : $rslt[0];
}
=item ($point, $intens, $central) = $sun->magnitude ($theta, $omega);
This method returns the magnitude of the Sun at a point $theta radians
from the center of its disk, given that the disk's angular radius
(B<not> diameter) is $omega radians. The returned $point is the
magnitude at the given point (undef if $theta > $omega), $intens is the
ratio of the intensity at the given point to the central intensity (0
if $theta > $omega), and $central is the central magnitude.
If this method is called in scalar context, it returns $point, the point
magnitude.
If the $omega argument is omitted or undefined, it is calculated based
on the geocentric range to the Sun at the current time setting of the
object.
If the $theta argument is omitted or undefined, the method returns
the average magnitude of the Sun, which is taken to be -26.8.
The limb-darkening algorithm and the associated constants come from
L<https://en.wikipedia.org/wiki/Limb_darkening>.
For consistency's sake, an observing station can optionally be passed as
the first argument (i.e. before C<$theta>). This is currently ignored.
=cut
{ # Begin local symbol block
my $central_mag;
my @limb_darkening = (.3, .93, -.23);
sub magnitude {
my ( $self, @args ) = @_;
# We want to accept a station as the second argument for
# consistency's sake, though we do not (at this point) use it.
embodies( $args[0], 'Astro::Coord::ECI' )
and shift @args;
my ( $theta, $omega ) = @args;
return MEAN_MAGNITUDE unless defined $theta;
unless (defined $omega) {
my @eci = $self->eci ();
$omega = $self->get ('diameter') / 2 /
sqrt (distsq (\@eci[0 .. 2], [0, 0, 0]));
}
unless (defined $central_mag) {
my $sum = 0;
my $quotient = 2;
foreach my $a (@limb_darkening) {
$sum += $a / $quotient++;
}
$central_mag = MEAN_MAGNITUDE - intensity_to_magnitude (2 * $sum);
}
my $intens = 0;
my $point;
if ($theta < $omega) {
lib/Astro/Coord/ECI/Sun.pm view on Meta::CPAN
# The above calculation gives the wrong $season between the
# December equinox and the end of the year. We fix it up here:
if ( $time - $season > SECSPERDAY * 180 ) {
$year++;
$season = $self->season( $year, $quarter );
}
# If we're a quarter too early, add one and repeat the
# calculation. We shouldn't have to do this more than once,
# since our maximum error even with the fudge factor is a day.
if ( $season < $time ) {
$quarter++;
if ( $quarter > 3 ) {
$quarter -= 4;
$year++;
}
$season = $self->season( $year, $quarter );
}
}
$season = ceil( $season ); # Make sure we're AFTER.
$self->universal( $season );
return wantarray ? ( $season, $quarter, $self->__quarter_name(
$quarter ) ) : $season;
}
=item $hash_reference = $sun->next_quarter_hash($want);
This convenience method wraps $sun->next_quarter(), but returns the
data in a hash reference, sort of like Astro::Coord::ECI::TLE->pass()
does. The hash contains the following keys:
{body} => the original object ($sun);
{almanac} => {
{event} => 'quarter',
{detail} => the quarter number (0 through 3);
{description} => the quarter description;
}
{time} => the time the quarter occurred.
The {time}, {detail}, and {description} keys correspond to elements 0
through 2 of the list returned by next_quarter().
=cut
use Astro::Coord::ECI::Mixin qw{ next_quarter_hash };
=item $period = $sun->period ()
Although this method is attached to an object that represents the
Sun, what it actually returns is the sidereal period of the Earth,
per Appendix I (pg 408) of Jean Meeus' "Astronomical Algorithms,"
2nd edition.
=cut
sub period {return 31558149.7632} # 365.256363 * 86400
sub __horizon_name_tplt {
my ( $self ) = @_;
return $self->__object_is_self_named() ?
[ '%sset', '%srise' ] :
$self->SUPER::__horizon_name_tplt();
}
sub __quarter_name {
my ( $self, $event, $tplt ) = @_;
$tplt ||= $self->__object_is_self_named() ?
[
'Spring equinox', 'Summer solstice',
'Fall equinox', 'Winter solstice',
] : [
'%s Spring equinox', '%s Summer solstice',
'%s Fall equinox', '%s Winter solstice',
];
my $station;
$station = $self->get( 'station' )
and ( $station->geodetic() )[0] < 0
and $event = ( $event + @{ $tplt } / 2 ) % @{ $tplt };
return $self->__event_name( $event, $tplt );
}
sub __transit_name_tplt {
my ( $self ) = @_;
return $self->__object_is_self_named() ?
[ 'local midnight', 'local noon' ] :
[ '%s transits nadir', '%s transits meridian' ];
}
sub __twilight_name {
my ( $self, $event, $tplt ) = @_;
$tplt ||= $self->__object_is_self_named() ?
[ 'end twilight', 'begin twilight' ] :
[ 'end %s twilight', 'begin %s twilight' ];
return $self->__event_name( $event, $tplt );
}
=begin comment
=item $time = $self->season( $year, $season );
This method calculates the time of the given season of the given
Gregorian year. The $season is an integer from 0 to 3, with 0 being the
astronomical Spring equinox (first point of Aries), and so on.
The algorithm comes from Jean Meeus' "Astronomical Algorithms", 2nd
Edition, Chapter 27 ("Equinoxes and Solstices), pages 278ff.
THIS METHOD IS UNSUPPORTED. I understand the temptation to call it if
all you want are the seasons, but if possible I would like to be able to
remove it if its use in next_quarter() turns out to be a bad idea. I am
not unwilling to support it; if you want me to, please contact me.
Because it is unsupported, its name may change without warning if
L<Test::Pod::Coverage|Test::Pod::Coverage> becomes smart enough to
realize that the =begin/end comment markers mean that this method is not
documented after all.
=end comment
=cut
sub season {
lib/Astro/Coord/ECI/Sun.pm view on Meta::CPAN
[ 203, 337.23, 32964.467 ],
[ 199, 342.08, 20.186 ],
[ 182, 27.85, 445267.112 ],
[ 156, 73.14, 45036.886 ],
[ 136, 171.52, 22518.443 ],
[ 77, 222.54, 65928.934 ],
[ 74, 296.72, 3034.906 ],
[ 70, 243.58, 9037.513 ],
[ 58, 119.81, 33718.147 ],
[ 52, 297.17, 150.678 ],
[ 50, 21.02, 2281.226 ],
[ 45, 247.54, 29929.562 ],
[ 44, 325.15, 31555.956 ],
[ 29, 60.93, 4443.417 ],
[ 18, 155.12, 67555.328 ],
[ 17, 288.79, 4562.452 ],
[ 16, 198.04, 62894.029 ],
[ 14, 199.76, 31436.921 ],
[ 12, 95.39, 14577.848 ],
[ 12, 287.11, 31931.756 ],
[ 12, 320.81, 34777.259 ],
[ 9, 227.73, 1222.114 ],
[ 8, 15.45, 16859.074 ],
) {
$S += $term->[0] * cos( mod2pi( deg2rad( $term->[1] + $term->[2]
* $T ) ) );
}
$self->{debug}
and print "Debug - S = $S\n";
my $JDE = 0.00001 * $S / $delta_lambda + $JDE0;
$self->{debug}
and print "Debug - JDE = $JDE\n";
my $time = ( $JDE - JD_OF_EPOCH ) * SECSPERDAY;
# Note that gmtime() in the following needs the parens because it
# might have come from Time::y2038, which appears to take more than
# one argument -- even though, as I read it, its prototype is (;$).
$self->{debug}
and print "Debug - dynamical date is ", scalar gmtime( $time ), "\n";
return $time - dynamical_delta( $time ); # Not quite right.
}
=item $sun->set( ... )
This method has been overridden to silently ignore any attempt to set
the C<'sun'> attribute.
=cut
sub set {
my ( $self, @args ) = @_;
while ( @args ) {
my ( $name, $val ) = splice @args, 0, 2;
if ( 'sun' eq $name ) {
} elsif ( $attrib{$name} ) {
if ( ref $self ) {
$self->{$name} = $val;
} else {
$static{$name} = $val;
}
} else {
$self->SUPER::set( $name, $val );
}
}
return $self;
}
=item $sun->time_set ()
This method sets coordinates of the object to the coordinates of the
Sun at the object's currently-set universal time. The velocity
components are arbitrarily set to 0. The 'equinox_dynamical' attribute
is set to the object's currently-set dynamical time.
Although there's no reason this method can't be called directly, it
exists to take advantage of the hook in the B<Astro::Coord::ECI>
object, to allow the position of the Sun to be computed when the
object's time is set.
The algorithm comes from Jean Meeus' "Astronomical Algorithms", 2nd
Edition, Chapter 25, pages 163ff.
=cut
# The following constants are used in the time_set() method,
# because Meeus' equations are in degrees, I was too lazy
# to hand-convert them to radians, but I didn't want to
# penalize the user for the conversion every time.
use constant SUN_C1_0 => deg2rad (1.914602);
use constant SUN_C1_1 => deg2rad (-0.004817);
use constant SUN_C1_2 => deg2rad (-0.000014);
use constant SUN_C2_0 => deg2rad (0.019993);
use constant SUN_C2_1 => deg2rad (0.000101);
use constant SUN_C3_0 => deg2rad (0.000289);
use constant SUN_LON_2000 => deg2rad (- 0.01397);
sub time_set {
my $self = shift;
my $time = $self->dynamical;
# The following algorithm is from Meeus, chapter 25, page, 163 ff.
my $T = jcent2000 ($time); # Meeus (25.1)
my $L0 = mod2pi(deg2rad((.0003032 * $T + 36000.76983) * $T # Meeus (25.2)
+ 280.46646));
my $M = mod2pi(deg2rad(((-.0001537) * $T + 35999.05029) # Meeus (25.3)
* $T + 357.52911));
my $e = (-0.0000001267 * $T - 0.000042037) * $T + 0.016708634;# Meeus (25.4)
my $C = ((SUN_C1_2 * $T + SUN_C1_1) * $T + SUN_C1_0) * sin ($M)
+ (SUN_C2_1 * $T + SUN_C2_0) * sin (2 * $M)
+ SUN_C3_0 * sin (3 * $M);
my $O = $self->{_sun_geometric_longitude} = $L0 + $C;
my $omega = mod2pi (deg2rad (125.04 - 1934.136 * $T));
my $lambda = mod2pi ($O - deg2rad (0.00569 + 0.00478 * sin ($omega)));
my $nu = $M + $C;
my $R = (1.000_001_018 * (1 - $e * $e)) / (1 + $e * cos ($nu))
* AU;
$self->{debug} and print <<eod;
Debug sun - @{[ gm_strftime '%d-%b-%Y %H:%M:%S', $time ]}
T = $T
L0 = @{[rad2deg ($L0)]} degrees
( run in 3.916 seconds using v1.01-cache-2.11-cpan-140bd7fdf52 )