line |
stmt |
bran |
cond |
sub |
pod |
time |
code |
1
|
|
|
|
|
|
|
=head1 NAME |
2
|
|
|
|
|
|
|
|
3
|
|
|
|
|
|
|
Astro::Coord::ECI::Sun - Compute the position of the Sun. |
4
|
|
|
|
|
|
|
|
5
|
|
|
|
|
|
|
=head1 SYNOPSIS |
6
|
|
|
|
|
|
|
|
7
|
|
|
|
|
|
|
use Astro::Coord::ECI; |
8
|
|
|
|
|
|
|
use Astro::Coord::ECI::Sun; |
9
|
|
|
|
|
|
|
use Astro::Coord::ECI::Utils qw{deg2rad}; |
10
|
|
|
|
|
|
|
|
11
|
|
|
|
|
|
|
# 1600 Pennsylvania Ave, Washington DC USA |
12
|
|
|
|
|
|
|
# latitude 38.899 N, longitude 77.038 W, |
13
|
|
|
|
|
|
|
# altitude 16.68 meters above sea level |
14
|
|
|
|
|
|
|
my $lat = deg2rad (38.899); # Radians |
15
|
|
|
|
|
|
|
my $long = deg2rad (-77.038); # Radians |
16
|
|
|
|
|
|
|
my $alt = 16.68 / 1000; # Kilometers |
17
|
|
|
|
|
|
|
my $sun = Astro::Coord::ECI::Sun->new (); |
18
|
|
|
|
|
|
|
my $sta = Astro::Coord::ECI-> |
19
|
|
|
|
|
|
|
universal (time ())-> |
20
|
|
|
|
|
|
|
geodetic ($lat, $long, $alt); |
21
|
|
|
|
|
|
|
my ($time, $rise) = $sta->next_elevation ($sun); |
22
|
|
|
|
|
|
|
print "Sun @{[$rise ? 'rise' : 'set']} is ", |
23
|
|
|
|
|
|
|
scalar gmtime $time, " UT\n"; |
24
|
|
|
|
|
|
|
|
25
|
|
|
|
|
|
|
Although this example computes the Sun rise or set in Washington D.C. |
26
|
|
|
|
|
|
|
USA, the time is displayed in Universal Time. This is because I did not |
27
|
|
|
|
|
|
|
want to complicate the example by adding machinery to convert the time |
28
|
|
|
|
|
|
|
to the correct zone for Washington D.C. (which is UT - 5 except when |
29
|
|
|
|
|
|
|
Summer Time is in effect, when it is UT - 4). |
30
|
|
|
|
|
|
|
|
31
|
|
|
|
|
|
|
=head1 DESCRIPTION |
32
|
|
|
|
|
|
|
|
33
|
|
|
|
|
|
|
This module implements the position of the Sun as a function of time, as |
34
|
|
|
|
|
|
|
described in Jean Meeus' "Astronomical Algorithms," second edition. It |
35
|
|
|
|
|
|
|
is a subclass of B, with the id, name, and diameter |
36
|
|
|
|
|
|
|
attributes initialized appropriately, and the time_set() method |
37
|
|
|
|
|
|
|
overridden to compute the position of the Sun at the given time. |
38
|
|
|
|
|
|
|
|
39
|
|
|
|
|
|
|
=head2 Methods |
40
|
|
|
|
|
|
|
|
41
|
|
|
|
|
|
|
The following methods should be considered public: |
42
|
|
|
|
|
|
|
|
43
|
|
|
|
|
|
|
=over |
44
|
|
|
|
|
|
|
|
45
|
|
|
|
|
|
|
=cut |
46
|
|
|
|
|
|
|
|
47
|
|
|
|
|
|
|
package Astro::Coord::ECI::Sun; |
48
|
|
|
|
|
|
|
|
49
|
16
|
|
|
16
|
|
1925
|
use strict; |
|
16
|
|
|
|
|
35
|
|
|
16
|
|
|
|
|
527
|
|
50
|
16
|
|
|
16
|
|
90
|
use warnings; |
|
16
|
|
|
|
|
34
|
|
|
16
|
|
|
|
|
825
|
|
51
|
|
|
|
|
|
|
|
52
|
|
|
|
|
|
|
our $VERSION = '0.129'; |
53
|
|
|
|
|
|
|
|
54
|
16
|
|
|
16
|
|
105
|
use base qw{Astro::Coord::ECI}; |
|
16
|
|
|
|
|
35
|
|
|
16
|
|
|
|
|
1964
|
|
55
|
|
|
|
|
|
|
|
56
|
16
|
|
|
16
|
|
108
|
use Astro::Coord::ECI::Utils qw{ @CARP_NOT :mainstream }; |
|
16
|
|
|
|
|
45
|
|
|
16
|
|
|
|
|
7542
|
|
57
|
16
|
|
|
16
|
|
127
|
use Carp; |
|
16
|
|
|
|
|
34
|
|
|
16
|
|
|
|
|
1102
|
|
58
|
16
|
|
|
16
|
|
103
|
use POSIX qw{ ceil floor strftime }; |
|
16
|
|
|
|
|
33
|
|
|
16
|
|
|
|
|
136
|
|
59
|
|
|
|
|
|
|
|
60
|
16
|
|
|
16
|
|
1489
|
use constant MEAN_MAGNITUDE => -26.8; |
|
16
|
|
|
|
|
35
|
|
|
16
|
|
|
|
|
6640
|
|
61
|
|
|
|
|
|
|
|
62
|
|
|
|
|
|
|
my %attrib = ( |
63
|
|
|
|
|
|
|
iterate_for_quarters => 1, |
64
|
|
|
|
|
|
|
); |
65
|
|
|
|
|
|
|
|
66
|
|
|
|
|
|
|
my %static = ( |
67
|
|
|
|
|
|
|
id => 'Sun', |
68
|
|
|
|
|
|
|
name => 'Sun', |
69
|
|
|
|
|
|
|
diameter => 1392000, |
70
|
|
|
|
|
|
|
iterate_for_quarters => undef, |
71
|
|
|
|
|
|
|
); |
72
|
|
|
|
|
|
|
|
73
|
|
|
|
|
|
|
my $weaken = eval { |
74
|
|
|
|
|
|
|
require Scalar::Util; |
75
|
|
|
|
|
|
|
Scalar::Util->can('weaken'); |
76
|
|
|
|
|
|
|
}; |
77
|
|
|
|
|
|
|
|
78
|
|
|
|
|
|
|
our $Singleton = $weaken; |
79
|
|
|
|
|
|
|
|
80
|
|
|
|
|
|
|
my %object; # By class |
81
|
|
|
|
|
|
|
|
82
|
|
|
|
|
|
|
=item $sun = Astro::Coord::ECI::Sun->new(); |
83
|
|
|
|
|
|
|
|
84
|
|
|
|
|
|
|
This method instantiates an object to represent the coordinates of the |
85
|
|
|
|
|
|
|
Sun. This is a subclass of L, with |
86
|
|
|
|
|
|
|
the id and name attributes set to 'Sun', and the diameter attribute set |
87
|
|
|
|
|
|
|
to 1392000 km per Jean Meeus' "Astronomical Algorithms", 2nd Edition, |
88
|
|
|
|
|
|
|
Appendix I, page 407. |
89
|
|
|
|
|
|
|
|
90
|
|
|
|
|
|
|
Any arguments are passed to the set() method once the object has been |
91
|
|
|
|
|
|
|
instantiated. Yes, you can override the "hard-wired" id, name, and so |
92
|
|
|
|
|
|
|
forth in this way. |
93
|
|
|
|
|
|
|
|
94
|
|
|
|
|
|
|
If C<$Astro::Coord::ECI::Sun::Singleton> is true, you get a singleton |
95
|
|
|
|
|
|
|
object; that is, only one object is instantiated and subsequent calls to |
96
|
|
|
|
|
|
|
C just return that object. If higher-accuracy subclasses are ever |
97
|
|
|
|
|
|
|
implemented, there will be one singleton for each class. |
98
|
|
|
|
|
|
|
|
99
|
|
|
|
|
|
|
The singleton logic only works if L exports |
100
|
|
|
|
|
|
|
C. If it does not, the setting of |
101
|
|
|
|
|
|
|
C<$Astro::Coord::ECI::Sun::Singleton> is silently ignored. The default |
102
|
|
|
|
|
|
|
is true if L can be loaded and exports |
103
|
|
|
|
|
|
|
C, and false otherwise. |
104
|
|
|
|
|
|
|
|
105
|
|
|
|
|
|
|
=cut |
106
|
|
|
|
|
|
|
|
107
|
|
|
|
|
|
|
sub new { |
108
|
134
|
|
|
134
|
1
|
1949
|
my ($class, @args) = @_; |
109
|
134
|
50
|
|
|
|
311
|
ref $class and $class = ref $class; |
110
|
134
|
100
|
66
|
|
|
688
|
if ( $Singleton && $weaken && __classisa( $class, __PACKAGE__ ) ) { |
|
|
|
66
|
|
|
|
|
111
|
132
|
|
|
|
|
214
|
my $self; |
112
|
132
|
100
|
|
|
|
362
|
if ( $self = $object{$class} ) { |
113
|
97
|
100
|
|
|
|
225
|
$self->set( @args ) if @args; |
114
|
97
|
|
|
|
|
293
|
return $self; |
115
|
|
|
|
|
|
|
} else { |
116
|
35
|
|
|
|
|
380
|
$self = $object{$class} = $class->SUPER::new (%static, @args); |
117
|
35
|
|
|
|
|
240
|
$weaken->( $object{$class} ); |
118
|
35
|
|
|
|
|
156
|
return $self; |
119
|
|
|
|
|
|
|
} |
120
|
|
|
|
|
|
|
} else { |
121
|
2
|
|
|
|
|
9
|
return $class->SUPER::new (%static, @args); |
122
|
|
|
|
|
|
|
} |
123
|
|
|
|
|
|
|
} |
124
|
|
|
|
|
|
|
|
125
|
|
|
|
|
|
|
=item @almanac = $sun->almanac( $station, $start, $end ); |
126
|
|
|
|
|
|
|
|
127
|
|
|
|
|
|
|
This method produces almanac data for the Sun for the given observing |
128
|
|
|
|
|
|
|
station, between the given start and end times. The station is assumed |
129
|
|
|
|
|
|
|
to be Earth-Fixed - that is, you can't do this for something in orbit. |
130
|
|
|
|
|
|
|
|
131
|
|
|
|
|
|
|
The C<$station> argument may be omitted if the C attribute has |
132
|
|
|
|
|
|
|
been set. That is, this method can also be called as |
133
|
|
|
|
|
|
|
|
134
|
|
|
|
|
|
|
@almanac = $sun->almanac( $start, $end ) |
135
|
|
|
|
|
|
|
|
136
|
|
|
|
|
|
|
The start time defaults to the current time setting of the $sun |
137
|
|
|
|
|
|
|
object, and the end time defaults to a day after the start time. |
138
|
|
|
|
|
|
|
|
139
|
|
|
|
|
|
|
The almanac data consists of a list of list references. Each list |
140
|
|
|
|
|
|
|
reference points to a list containing the following elements: |
141
|
|
|
|
|
|
|
|
142
|
|
|
|
|
|
|
[0] => time |
143
|
|
|
|
|
|
|
[1] => event (string) |
144
|
|
|
|
|
|
|
[2] => detail (integer) |
145
|
|
|
|
|
|
|
[3] => description (string) |
146
|
|
|
|
|
|
|
|
147
|
|
|
|
|
|
|
The @almanac list is returned sorted by time. |
148
|
|
|
|
|
|
|
|
149
|
|
|
|
|
|
|
The following events, details, and descriptions are at least |
150
|
|
|
|
|
|
|
potentially returned: |
151
|
|
|
|
|
|
|
|
152
|
|
|
|
|
|
|
horizon: 0 = Sunset, 1 = Sunrise; |
153
|
|
|
|
|
|
|
transit: 0 = local midnight, 1 = local noon; |
154
|
|
|
|
|
|
|
twilight: 0 = end twilight, 1 = begin twilight; |
155
|
|
|
|
|
|
|
quarter: 0 = spring equinox, 1 = summer solstice, |
156
|
|
|
|
|
|
|
2 = fall equinox, 3 = winter solstice. |
157
|
|
|
|
|
|
|
|
158
|
|
|
|
|
|
|
Twilight is calculated based on the current value of the twilight |
159
|
|
|
|
|
|
|
attribute of the $sun object. This attribute is inherited from |
160
|
|
|
|
|
|
|
L, and documented there. |
161
|
|
|
|
|
|
|
|
162
|
|
|
|
|
|
|
=cut |
163
|
|
|
|
|
|
|
|
164
|
|
|
|
|
|
|
sub __almanac_event_type_iterator { |
165
|
2
|
|
|
2
|
|
6
|
my ( $self, $station ) = @_; |
166
|
|
|
|
|
|
|
|
167
|
2
|
|
|
|
|
4
|
my $inx = 0; |
168
|
|
|
|
|
|
|
|
169
|
2
|
|
|
|
|
8
|
my $horizon = $station->__get_almanac_horizon(); |
170
|
|
|
|
|
|
|
|
171
|
2
|
|
|
|
|
13
|
my @events = ( |
172
|
|
|
|
|
|
|
[ $station, next_elevation => [ $self, $horizon, 1 ], |
173
|
|
|
|
|
|
|
horizon => '__horizon_name' ], |
174
|
|
|
|
|
|
|
[ $station, next_meridian => [ $self ], |
175
|
|
|
|
|
|
|
transit => '__transit_name' ], |
176
|
|
|
|
|
|
|
[ $station, next_elevation => |
177
|
|
|
|
|
|
|
[ $self, $self->get( 'twilight' ) + $horizon, 0 ], |
178
|
|
|
|
|
|
|
twilight => '__twilight_name' ], |
179
|
|
|
|
|
|
|
[ $self, next_quarter => [], quarter => '__quarter_name', ], |
180
|
|
|
|
|
|
|
); |
181
|
|
|
|
|
|
|
|
182
|
|
|
|
|
|
|
return sub { |
183
|
|
|
|
|
|
|
$inx < @events |
184
|
10
|
100
|
|
10
|
|
38
|
and return @{ $events[$inx++] }; |
|
8
|
|
|
|
|
53
|
|
185
|
2
|
|
|
|
|
9
|
return; |
186
|
2
|
|
|
|
|
13
|
}; |
187
|
|
|
|
|
|
|
} |
188
|
|
|
|
|
|
|
|
189
|
16
|
|
|
16
|
|
5961
|
use Astro::Coord::ECI::Mixin qw{ almanac }; |
|
16
|
|
|
|
|
47
|
|
|
16
|
|
|
|
|
973
|
|
190
|
|
|
|
|
|
|
|
191
|
|
|
|
|
|
|
=item @almanac = $sun->almanac_hash( $station, $start, $end ); |
192
|
|
|
|
|
|
|
|
193
|
|
|
|
|
|
|
This convenience method wraps $sun->almanac(), but returns a list of |
194
|
|
|
|
|
|
|
hash references, sort of like Astro::Coord::ECI::TLE->pass() |
195
|
|
|
|
|
|
|
does. The hashes contain the following keys: |
196
|
|
|
|
|
|
|
|
197
|
|
|
|
|
|
|
{almanac} => { |
198
|
|
|
|
|
|
|
{event} => the event type; |
199
|
|
|
|
|
|
|
{detail} => the event detail (typically 0 or 1); |
200
|
|
|
|
|
|
|
{description} => the event description; |
201
|
|
|
|
|
|
|
} |
202
|
|
|
|
|
|
|
{body} => the original object ($sun); |
203
|
|
|
|
|
|
|
{station} => the observing station; |
204
|
|
|
|
|
|
|
{time} => the time the quarter occurred. |
205
|
|
|
|
|
|
|
|
206
|
|
|
|
|
|
|
The {time}, {event}, {detail}, and {description} keys correspond to |
207
|
|
|
|
|
|
|
elements 0 through 3 of the list returned by almanac(). |
208
|
|
|
|
|
|
|
|
209
|
|
|
|
|
|
|
=cut |
210
|
|
|
|
|
|
|
|
211
|
16
|
|
|
16
|
|
119
|
use Astro::Coord::ECI::Mixin qw{ almanac_hash }; |
|
16
|
|
|
|
|
34
|
|
|
16
|
|
|
|
|
9255
|
|
212
|
|
|
|
|
|
|
|
213
|
|
|
|
|
|
|
=item $coord2 = $coord->clone (); |
214
|
|
|
|
|
|
|
|
215
|
|
|
|
|
|
|
If singleton objects are enabled, this override of the superclass' |
216
|
|
|
|
|
|
|
method simply returns the invocant. Otherwise it does a deep clone of an |
217
|
|
|
|
|
|
|
object, producing a different but identical object. |
218
|
|
|
|
|
|
|
|
219
|
|
|
|
|
|
|
Prior to version 0.099_01 it always returned a clone. Yes, |
220
|
|
|
|
|
|
|
this is a change in long-standing functionality, but a long-standing bug |
221
|
|
|
|
|
|
|
is still a bug. |
222
|
|
|
|
|
|
|
|
223
|
|
|
|
|
|
|
=cut |
224
|
|
|
|
|
|
|
|
225
|
|
|
|
|
|
|
sub clone { |
226
|
2
|
|
|
2
|
1
|
612
|
my ( $self ) = @_; |
227
|
2
|
100
|
66
|
|
|
11
|
$Singleton |
228
|
|
|
|
|
|
|
and $weaken |
229
|
|
|
|
|
|
|
and return $self; |
230
|
1
|
|
|
|
|
10
|
return $self->SUPER::clone(); |
231
|
|
|
|
|
|
|
} |
232
|
|
|
|
|
|
|
|
233
|
|
|
|
|
|
|
=item $elevation = $tle->correct_for_refraction( $elevation ) |
234
|
|
|
|
|
|
|
|
235
|
|
|
|
|
|
|
This override of the superclass' method simply returns the elevation |
236
|
|
|
|
|
|
|
passed to it. I have no algorithm for refraction at the surface of the |
237
|
|
|
|
|
|
|
photosphere or anywhere else in the environs of the Sun, and explaining |
238
|
|
|
|
|
|
|
why I make no correction at all seemed easier than explaining why I make |
239
|
|
|
|
|
|
|
an incorrect correction. |
240
|
|
|
|
|
|
|
|
241
|
|
|
|
|
|
|
See the L C and |
242
|
|
|
|
|
|
|
C documentation for whether this class' |
243
|
|
|
|
|
|
|
C method is actually called by those methods. |
244
|
|
|
|
|
|
|
|
245
|
|
|
|
|
|
|
=cut |
246
|
|
|
|
|
|
|
|
247
|
|
|
|
|
|
|
sub correct_for_refraction { |
248
|
0
|
|
|
0
|
1
|
0
|
my ( undef, $elevation ) = @_; # Invocant unused |
249
|
0
|
|
|
|
|
0
|
return $elevation; |
250
|
|
|
|
|
|
|
} |
251
|
|
|
|
|
|
|
|
252
|
|
|
|
|
|
|
=item $long = $sun->geometric_longitude () |
253
|
|
|
|
|
|
|
|
254
|
|
|
|
|
|
|
This method returns the geometric longitude of the Sun in radians at |
255
|
|
|
|
|
|
|
the last time set. |
256
|
|
|
|
|
|
|
|
257
|
|
|
|
|
|
|
=cut |
258
|
|
|
|
|
|
|
|
259
|
|
|
|
|
|
|
sub geometric_longitude { |
260
|
1125
|
|
|
1125
|
1
|
1748
|
my $self = shift; |
261
|
1125
|
50
|
|
|
|
2387
|
croak <{_sun_geometric_longitude}; |
262
|
|
|
|
|
|
|
Error - You must set the time of the Sun object before the geometric |
263
|
|
|
|
|
|
|
longitude can be returned. |
264
|
|
|
|
|
|
|
eod |
265
|
|
|
|
|
|
|
|
266
|
1125
|
|
|
|
|
2416
|
return $self->{_sun_geometric_longitude}; |
267
|
|
|
|
|
|
|
} |
268
|
|
|
|
|
|
|
|
269
|
|
|
|
|
|
|
=item $sun->get( ... ) |
270
|
|
|
|
|
|
|
|
271
|
|
|
|
|
|
|
This method has been overridden to return the invocant as the C<'sun'> |
272
|
|
|
|
|
|
|
attribute. |
273
|
|
|
|
|
|
|
|
274
|
|
|
|
|
|
|
=cut |
275
|
|
|
|
|
|
|
|
276
|
|
|
|
|
|
|
sub get { |
277
|
6742
|
|
|
6742
|
1
|
13362
|
my ( $self, @args ) = @_; |
278
|
6742
|
|
|
|
|
9649
|
my @rslt; |
279
|
6742
|
|
|
|
|
11651
|
foreach my $name ( @args ) { |
280
|
|
|
|
|
|
|
push @rslt, 'sun' eq $name ? $self : |
281
|
|
|
|
|
|
|
$attrib{$name} ? |
282
|
6742
|
0
|
|
|
|
25575
|
ref $self ? $self->{$name} : $static{$name} : |
|
|
50
|
|
|
|
|
|
|
|
50
|
|
|
|
|
|
283
|
|
|
|
|
|
|
$self->SUPER::get( $name ); |
284
|
|
|
|
|
|
|
} |
285
|
6742
|
100
|
|
|
|
21342
|
return wantarray ? @rslt : $rslt[0]; |
286
|
|
|
|
|
|
|
} |
287
|
|
|
|
|
|
|
|
288
|
|
|
|
|
|
|
=item ($point, $intens, $central) = $sun->magnitude ($theta, $omega); |
289
|
|
|
|
|
|
|
|
290
|
|
|
|
|
|
|
This method returns the magnitude of the Sun at a point $theta radians |
291
|
|
|
|
|
|
|
from the center of its disk, given that the disk's angular radius |
292
|
|
|
|
|
|
|
(B diameter) is $omega radians. The returned $point is the |
293
|
|
|
|
|
|
|
magnitude at the given point (undef if $theta > $omega), $intens is the |
294
|
|
|
|
|
|
|
ratio of the intensity at the given point to the central intensity (0 |
295
|
|
|
|
|
|
|
if $theta > $omega), and $central is the central magnitude. |
296
|
|
|
|
|
|
|
|
297
|
|
|
|
|
|
|
If this method is called in scalar context, it returns $point, the point |
298
|
|
|
|
|
|
|
magnitude. |
299
|
|
|
|
|
|
|
|
300
|
|
|
|
|
|
|
If the $omega argument is omitted or undefined, it is calculated based |
301
|
|
|
|
|
|
|
on the geocentric range to the Sun at the current time setting of the |
302
|
|
|
|
|
|
|
object. |
303
|
|
|
|
|
|
|
|
304
|
|
|
|
|
|
|
If the $theta argument is omitted or undefined, the method returns |
305
|
|
|
|
|
|
|
the average magnitude of the Sun, which is taken to be -26.8. |
306
|
|
|
|
|
|
|
|
307
|
|
|
|
|
|
|
The limb-darkening algorithm and the associated constants come from |
308
|
|
|
|
|
|
|
L. |
309
|
|
|
|
|
|
|
|
310
|
|
|
|
|
|
|
For consistency's sake, an observing station can optionally be passed as |
311
|
|
|
|
|
|
|
the first argument (i.e. before C<$theta>). This is currently ignored. |
312
|
|
|
|
|
|
|
|
313
|
|
|
|
|
|
|
=cut |
314
|
|
|
|
|
|
|
|
315
|
|
|
|
|
|
|
{ # Begin local symbol block |
316
|
|
|
|
|
|
|
|
317
|
|
|
|
|
|
|
my $central_mag; |
318
|
|
|
|
|
|
|
my @limb_darkening = (.3, .93, -.23); |
319
|
|
|
|
|
|
|
|
320
|
|
|
|
|
|
|
sub magnitude { |
321
|
0
|
|
|
0
|
1
|
0
|
my ( $self, @args ) = @_; |
322
|
|
|
|
|
|
|
# We want to accept a station as the second argument for |
323
|
|
|
|
|
|
|
# consistency's sake, though we do not (at this point) use it. |
324
|
0
|
0
|
|
|
|
0
|
embodies( $args[0], 'Astro::Coord::ECI' ) |
325
|
|
|
|
|
|
|
and shift @args; |
326
|
0
|
|
|
|
|
0
|
my ( $theta, $omega ) = @args; |
327
|
0
|
0
|
|
|
|
0
|
return MEAN_MAGNITUDE unless defined $theta; |
328
|
0
|
0
|
|
|
|
0
|
unless (defined $omega) { |
329
|
0
|
|
|
|
|
0
|
my @eci = $self->eci (); |
330
|
0
|
|
|
|
|
0
|
$omega = $self->get ('diameter') / 2 / |
331
|
|
|
|
|
|
|
sqrt (distsq (\@eci[0 .. 2], [0, 0, 0])); |
332
|
|
|
|
|
|
|
} |
333
|
0
|
0
|
|
|
|
0
|
unless (defined $central_mag) { |
334
|
0
|
|
|
|
|
0
|
my $sum = 0; |
335
|
0
|
|
|
|
|
0
|
my $quotient = 2; |
336
|
0
|
|
|
|
|
0
|
foreach my $a (@limb_darkening) { |
337
|
0
|
|
|
|
|
0
|
$sum += $a / $quotient++; |
338
|
|
|
|
|
|
|
} |
339
|
0
|
|
|
|
|
0
|
$central_mag = MEAN_MAGNITUDE - intensity_to_magnitude (2 * $sum); |
340
|
|
|
|
|
|
|
} |
341
|
0
|
|
|
|
|
0
|
my $intens = 0; |
342
|
0
|
|
|
|
|
0
|
my $point; |
343
|
0
|
0
|
|
|
|
0
|
if ($theta < $omega) { |
344
|
0
|
|
|
|
|
0
|
my $costheta = cos ($theta); |
345
|
0
|
|
|
|
|
0
|
my $cosomega = cos ($omega); |
346
|
0
|
|
|
|
|
0
|
my $sinomega = sin ($omega); |
347
|
0
|
|
|
|
|
0
|
my $cospsi = sqrt ($costheta * $costheta - |
348
|
|
|
|
|
|
|
$cosomega * $cosomega) / $sinomega; |
349
|
0
|
|
|
|
|
0
|
my $psiterm = 1; |
350
|
0
|
|
|
|
|
0
|
foreach my $a (@limb_darkening) { |
351
|
0
|
|
|
|
|
0
|
$intens += $a * $psiterm; |
352
|
0
|
|
|
|
|
0
|
$psiterm *= $cospsi; |
353
|
|
|
|
|
|
|
} |
354
|
0
|
|
|
|
|
0
|
$point = $central_mag + intensity_to_magnitude ($intens); |
355
|
|
|
|
|
|
|
} |
356
|
0
|
0
|
|
|
|
0
|
return wantarray ? ($point, $intens, $central_mag) : $point; |
357
|
|
|
|
|
|
|
} |
358
|
|
|
|
|
|
|
} # End local symbol block. |
359
|
|
|
|
|
|
|
|
360
|
|
|
|
|
|
|
=item ($time, $quarter, $desc) = $sun->next_quarter($want); |
361
|
|
|
|
|
|
|
|
362
|
|
|
|
|
|
|
This method calculates the time of the next equinox or solstice after |
363
|
|
|
|
|
|
|
the current time setting of the $sun object. The return is the time, |
364
|
|
|
|
|
|
|
which equinox or solstice it is as a number from 0 (March equinox) to 3 |
365
|
|
|
|
|
|
|
(December solstice), and a string describing the equinox or solstice. If |
366
|
|
|
|
|
|
|
called in scalar context, you just get the time. |
367
|
|
|
|
|
|
|
|
368
|
|
|
|
|
|
|
If the C attribute is not set or set to a location on or north |
369
|
|
|
|
|
|
|
of the Equator, the descriptor strings are |
370
|
|
|
|
|
|
|
|
371
|
|
|
|
|
|
|
0 - Spring equinox |
372
|
|
|
|
|
|
|
1 - Summer solstice |
373
|
|
|
|
|
|
|
2 - Fall equinox |
374
|
|
|
|
|
|
|
3 - Winter solstice |
375
|
|
|
|
|
|
|
|
376
|
|
|
|
|
|
|
If the C attribute is set to a location south of the Equator, |
377
|
|
|
|
|
|
|
the descriptor strings are |
378
|
|
|
|
|
|
|
|
379
|
|
|
|
|
|
|
0 - Fall equinox |
380
|
|
|
|
|
|
|
1 - Winter solstice |
381
|
|
|
|
|
|
|
2 - Spring equinox |
382
|
|
|
|
|
|
|
3 - Summer solstice |
383
|
|
|
|
|
|
|
|
384
|
|
|
|
|
|
|
The optional $want argument says which equinox or solstice you want, as |
385
|
|
|
|
|
|
|
a number from 0 through 3. |
386
|
|
|
|
|
|
|
|
387
|
|
|
|
|
|
|
As a side effect, the time of the $sun object ends up set to the |
388
|
|
|
|
|
|
|
returned time. |
389
|
|
|
|
|
|
|
|
390
|
|
|
|
|
|
|
As of version 0.088_01, the algorithm given in Jean Meeus' "Astronomical |
391
|
|
|
|
|
|
|
Algorithms", 2nd Edition, Chapter 27 ("Equinoxes and Solstices"), pages |
392
|
|
|
|
|
|
|
278ff is used. This should be good for the range -1000 to 3000 |
393
|
|
|
|
|
|
|
Gregorian, and good to within a minute or so within the range 1951 to |
394
|
|
|
|
|
|
|
2050 Gregorian, but the longitude of the Sun at the calculated time may |
395
|
|
|
|
|
|
|
be as much as 0.01 degree off the exact time for the event. |
396
|
|
|
|
|
|
|
|
397
|
|
|
|
|
|
|
If you take the United States Naval Observatory's times (given to the |
398
|
|
|
|
|
|
|
nearest minute) as the standard, the maximum deviation from that |
399
|
|
|
|
|
|
|
standard in the range 1700 to 2100 is 226 seconds. I have no information |
400
|
|
|
|
|
|
|
on this algorithm's accuracy outside that range. I. |
401
|
|
|
|
|
|
|
|
402
|
|
|
|
|
|
|
In version 0.088 and before, this calculation was done by successive |
403
|
|
|
|
|
|
|
approximation based on the position of the Sun, and was good to about 15 |
404
|
|
|
|
|
|
|
minutes. |
405
|
|
|
|
|
|
|
|
406
|
|
|
|
|
|
|
If you want the old iterative version back, set attribute |
407
|
|
|
|
|
|
|
C to a true value. |
408
|
|
|
|
|
|
|
|
409
|
|
|
|
|
|
|
=cut |
410
|
|
|
|
|
|
|
|
411
|
16
|
|
|
16
|
|
127
|
use constant NEXT_QUARTER_INCREMENT => 86400 * 85; # 85 days. |
|
16
|
|
|
|
|
34
|
|
|
16
|
|
|
|
|
4959
|
|
412
|
|
|
|
|
|
|
|
413
|
|
|
|
|
|
|
*__next_quarter_coordinate = __PACKAGE__->can( |
414
|
|
|
|
|
|
|
'ecliptic_longitude' ); |
415
|
|
|
|
|
|
|
|
416
|
|
|
|
|
|
|
# use Astro::Coord::ECI::Mixin qw{ next_quarter }; |
417
|
|
|
|
|
|
|
|
418
|
|
|
|
|
|
|
sub next_quarter { |
419
|
7
|
|
|
7
|
1
|
1308
|
my ( $self, $quarter ) = @_; |
420
|
|
|
|
|
|
|
$self->{iterate_for_quarters} |
421
|
7
|
50
|
|
|
|
21
|
and goto &Astro::Coord::ECI::Mixin::next_quarter; |
422
|
7
|
|
|
|
|
21
|
my $time = $self->universal(); |
423
|
7
|
|
|
|
|
31
|
my $year = ( gmtime( $time ) )[5] + 1900; |
424
|
7
|
|
|
|
|
15
|
my $season; |
425
|
7
|
50
|
|
|
|
18
|
if ( defined $quarter ) { |
426
|
|
|
|
|
|
|
# I can't think of an edge case that makes the first calculation |
427
|
|
|
|
|
|
|
# give a quarter that is too late. Wish that made me feel better |
428
|
|
|
|
|
|
|
# than it in fact does. |
429
|
0
|
|
|
|
|
0
|
$season = $self->season( $year, $quarter ); |
430
|
0
|
0
|
|
|
|
0
|
$season < $time |
431
|
|
|
|
|
|
|
and $season = $self->season( $year + 1, $quarter ); |
432
|
0
|
|
|
|
|
0
|
$self->universal( $season ); |
433
|
|
|
|
|
|
|
} else { |
434
|
7
|
|
|
|
|
19
|
my ( undef, $lon ) = $self->ecliptic(); |
435
|
|
|
|
|
|
|
# The fudged-in 359 is because I am worried about the limited |
436
|
|
|
|
|
|
|
# accuracy of the longitude calculation causing me to pick the |
437
|
|
|
|
|
|
|
# wrong quarter. So I essentially back the Sun up by about a |
438
|
|
|
|
|
|
|
# day. That may indeed put me a quarter too early, but I have to |
439
|
|
|
|
|
|
|
# check for that anyway because of the accuracy (or lack |
440
|
|
|
|
|
|
|
# thereof) of the calculated longitude. |
441
|
7
|
|
|
|
|
20
|
$quarter = ceil( ( rad2deg( $lon ) + 359 ) / 90 ) % 4; |
442
|
7
|
|
|
|
|
26
|
$season = $self->season( $year, $quarter ); |
443
|
|
|
|
|
|
|
# The above calculation gives the wrong $season between the |
444
|
|
|
|
|
|
|
# December equinox and the end of the year. We fix it up here: |
445
|
7
|
100
|
|
|
|
20
|
if ( $time - $season > SECSPERDAY * 180 ) { |
446
|
1
|
|
|
|
|
4
|
$year++; |
447
|
1
|
|
|
|
|
5
|
$season = $self->season( $year, $quarter ); |
448
|
|
|
|
|
|
|
} |
449
|
|
|
|
|
|
|
# If we're a quarter too early, add one and repeat the |
450
|
|
|
|
|
|
|
# calculation. We shouldn't have to do this more than once, |
451
|
|
|
|
|
|
|
# since our maximum error even with the fudge factor is a day. |
452
|
7
|
100
|
|
|
|
22
|
if ( $season < $time ) { |
453
|
3
|
|
|
|
|
5
|
$quarter++; |
454
|
3
|
50
|
|
|
|
8
|
if ( $quarter > 3 ) { |
455
|
0
|
|
|
|
|
0
|
$quarter -= 4; |
456
|
0
|
|
|
|
|
0
|
$year++; |
457
|
|
|
|
|
|
|
} |
458
|
3
|
|
|
|
|
8
|
$season = $self->season( $year, $quarter ); |
459
|
|
|
|
|
|
|
} |
460
|
|
|
|
|
|
|
} |
461
|
7
|
|
|
|
|
19
|
$season = ceil( $season ); # Make sure we're AFTER. |
462
|
7
|
|
|
|
|
24
|
$self->universal( $season ); |
463
|
7
|
100
|
|
|
|
35
|
return wantarray ? ( $season, $quarter, $self->__quarter_name( |
464
|
|
|
|
|
|
|
$quarter ) ) : $season; |
465
|
|
|
|
|
|
|
} |
466
|
|
|
|
|
|
|
|
467
|
|
|
|
|
|
|
=item $hash_reference = $sun->next_quarter_hash($want); |
468
|
|
|
|
|
|
|
|
469
|
|
|
|
|
|
|
This convenience method wraps $sun->next_quarter(), but returns the |
470
|
|
|
|
|
|
|
data in a hash reference, sort of like Astro::Coord::ECI::TLE->pass() |
471
|
|
|
|
|
|
|
does. The hash contains the following keys: |
472
|
|
|
|
|
|
|
|
473
|
|
|
|
|
|
|
{body} => the original object ($sun); |
474
|
|
|
|
|
|
|
{almanac} => { |
475
|
|
|
|
|
|
|
{event} => 'quarter', |
476
|
|
|
|
|
|
|
{detail} => the quarter number (0 through 3); |
477
|
|
|
|
|
|
|
{description} => the quarter description; |
478
|
|
|
|
|
|
|
} |
479
|
|
|
|
|
|
|
{time} => the time the quarter occurred. |
480
|
|
|
|
|
|
|
|
481
|
|
|
|
|
|
|
The {time}, {detail}, and {description} keys correspond to elements 0 |
482
|
|
|
|
|
|
|
through 2 of the list returned by next_quarter(). |
483
|
|
|
|
|
|
|
|
484
|
|
|
|
|
|
|
=cut |
485
|
|
|
|
|
|
|
|
486
|
16
|
|
|
16
|
|
133
|
use Astro::Coord::ECI::Mixin qw{ next_quarter_hash }; |
|
16
|
|
|
|
|
35
|
|
|
16
|
|
|
|
|
13487
|
|
487
|
|
|
|
|
|
|
|
488
|
|
|
|
|
|
|
=item $period = $sun->period () |
489
|
|
|
|
|
|
|
|
490
|
|
|
|
|
|
|
Although this method is attached to an object that represents the |
491
|
|
|
|
|
|
|
Sun, what it actually returns is the sidereal period of the Earth, |
492
|
|
|
|
|
|
|
per Appendix I (pg 408) of Jean Meeus' "Astronomical Algorithms," |
493
|
|
|
|
|
|
|
2nd edition. |
494
|
|
|
|
|
|
|
|
495
|
|
|
|
|
|
|
=cut |
496
|
|
|
|
|
|
|
|
497
|
246
|
|
|
246
|
1
|
963
|
sub period {return 31558149.7632} # 365.256363 * 86400 |
498
|
|
|
|
|
|
|
|
499
|
|
|
|
|
|
|
sub __horizon_name_tplt { |
500
|
4
|
|
|
4
|
|
11
|
my ( $self ) = @_; |
501
|
4
|
50
|
|
|
|
25
|
return $self->__object_is_self_named() ? |
502
|
|
|
|
|
|
|
[ '%sset', '%srise' ] : |
503
|
|
|
|
|
|
|
$self->SUPER::__horizon_name_tplt(); |
504
|
|
|
|
|
|
|
} |
505
|
|
|
|
|
|
|
|
506
|
|
|
|
|
|
|
sub __quarter_name { |
507
|
2
|
|
|
2
|
|
7
|
my ( $self, $event, $tplt ) = @_; |
508
|
2
|
50
|
33
|
|
|
15
|
$tplt ||= $self->__object_is_self_named() ? |
509
|
|
|
|
|
|
|
[ |
510
|
|
|
|
|
|
|
'Spring equinox', 'Summer solstice', |
511
|
|
|
|
|
|
|
'Fall equinox', 'Winter solstice', |
512
|
|
|
|
|
|
|
] : [ |
513
|
|
|
|
|
|
|
'%s Spring equinox', '%s Summer solstice', |
514
|
|
|
|
|
|
|
'%s Fall equinox', '%s Winter solstice', |
515
|
|
|
|
|
|
|
]; |
516
|
2
|
|
|
|
|
5
|
my $station; |
517
|
|
|
|
|
|
|
$station = $self->get( 'station' ) |
518
|
|
|
|
|
|
|
and ( $station->geodetic() )[0] < 0 |
519
|
2
|
50
|
66
|
|
|
6
|
and $event = ( $event + @{ $tplt } / 2 ) % @{ $tplt }; |
|
0
|
|
|
|
|
0
|
|
|
0
|
|
|
|
|
0
|
|
520
|
2
|
|
|
|
|
10
|
return $self->__event_name( $event, $tplt ); |
521
|
|
|
|
|
|
|
} |
522
|
|
|
|
|
|
|
|
523
|
|
|
|
|
|
|
sub __transit_name_tplt { |
524
|
4
|
|
|
4
|
|
10
|
my ( $self ) = @_; |
525
|
4
|
50
|
|
|
|
13
|
return $self->__object_is_self_named() ? |
526
|
|
|
|
|
|
|
[ 'local midnight', 'local noon' ] : |
527
|
|
|
|
|
|
|
[ '%s transits nadir', '%s transits meridian' ]; |
528
|
|
|
|
|
|
|
} |
529
|
|
|
|
|
|
|
|
530
|
|
|
|
|
|
|
sub __twilight_name { |
531
|
4
|
|
|
4
|
|
12
|
my ( $self, $event, $tplt ) = @_; |
532
|
4
|
50
|
33
|
|
|
27
|
$tplt ||= $self->__object_is_self_named() ? |
533
|
|
|
|
|
|
|
[ 'end twilight', 'begin twilight' ] : |
534
|
|
|
|
|
|
|
[ 'end %s twilight', 'begin %s twilight' ]; |
535
|
4
|
|
|
|
|
16
|
return $self->__event_name( $event, $tplt ); |
536
|
|
|
|
|
|
|
} |
537
|
|
|
|
|
|
|
|
538
|
|
|
|
|
|
|
=begin comment |
539
|
|
|
|
|
|
|
|
540
|
|
|
|
|
|
|
=item $time = $self->season( $year, $season ); |
541
|
|
|
|
|
|
|
|
542
|
|
|
|
|
|
|
This method calculates the time of the given season of the given |
543
|
|
|
|
|
|
|
Gregorian year. The $season is an integer from 0 to 3, with 0 being the |
544
|
|
|
|
|
|
|
astronomical Spring equinox (first point of Aries), and so on. |
545
|
|
|
|
|
|
|
|
546
|
|
|
|
|
|
|
The algorithm comes from Jean Meeus' "Astronomical Algorithms", 2nd |
547
|
|
|
|
|
|
|
Edition, Chapter 27 ("Equinoxes and Solstices), pages 278ff. |
548
|
|
|
|
|
|
|
|
549
|
|
|
|
|
|
|
THIS METHOD IS UNSUPPORTED. I understand the temptation to call it if |
550
|
|
|
|
|
|
|
all you want are the seasons, but if possible I would like to be able to |
551
|
|
|
|
|
|
|
remove it if its use in next_quarter() turns out to e a bad idea. I am |
552
|
|
|
|
|
|
|
not unwilling to support it; if you want me to, please contact me. |
553
|
|
|
|
|
|
|
|
554
|
|
|
|
|
|
|
Because it is unsupported, its name may change without warning if |
555
|
|
|
|
|
|
|
L becomes smart enough to |
556
|
|
|
|
|
|
|
realize that the =begin/end comment markers mean that this method is not |
557
|
|
|
|
|
|
|
documented after all. |
558
|
|
|
|
|
|
|
|
559
|
|
|
|
|
|
|
=end comment |
560
|
|
|
|
|
|
|
|
561
|
|
|
|
|
|
|
=cut |
562
|
|
|
|
|
|
|
|
563
|
|
|
|
|
|
|
sub season { |
564
|
11
|
|
|
11
|
1
|
50
|
my ( $self, $year, $season ) = @_; |
565
|
11
|
50
|
|
|
|
108
|
my ( $Y, $d ) = $year < 1000 ? ( |
566
|
|
|
|
|
|
|
$year / 1000, |
567
|
|
|
|
|
|
|
[ # Meeus table 27 A |
568
|
|
|
|
|
|
|
[ -0.00071, 0.00111, 0.06134, 365242.13740, 1721139.29189 ], |
569
|
|
|
|
|
|
|
[ 0.00025, 0.00907, -0.05323, 365241.72562, 1721233.25401 ], |
570
|
|
|
|
|
|
|
[ 0.00074, -0.00297, -0.11677, 365242.49558, 1721325.70455 ], |
571
|
|
|
|
|
|
|
[ -0.00006, -0.00933, -0.00769, 365242.88257, 1721414.39987 ], |
572
|
|
|
|
|
|
|
]->[ $season ], |
573
|
|
|
|
|
|
|
) : ( |
574
|
|
|
|
|
|
|
( $year - 2000 ) / 1000, |
575
|
|
|
|
|
|
|
[ # Meeus table 27 B |
576
|
|
|
|
|
|
|
[ -0.00057, -0.00411, 0.05169, 365242.37404, 2451623.80984 ], |
577
|
|
|
|
|
|
|
[ -0.00030, 0.00888, 0.00325, 365241.62603, 2451716.56767 ], |
578
|
|
|
|
|
|
|
[ 0.00078, 0.00337, -0.11575, 365242.01767, 2451810.21715 ], |
579
|
|
|
|
|
|
|
[ 0.00032, -0.00823, -0.06223, 365242.74049, 2451900.05952 ], |
580
|
|
|
|
|
|
|
]->[ $season ], |
581
|
|
|
|
|
|
|
); |
582
|
11
|
|
|
|
|
47
|
my $JDE0 = ( ( ( $d->[0] * $Y + $d->[1] ) * $Y + $d->[2] ) * $Y + |
583
|
|
|
|
|
|
|
$d->[3] ) * $Y + $d->[4]; |
584
|
11
|
|
|
|
|
19
|
my $T = ( $JDE0 - 2451545.0 ) / 36525; |
585
|
|
|
|
|
|
|
$self->{debug} |
586
|
11
|
50
|
|
|
|
26
|
and print "Debug - T = $T\n"; |
587
|
11
|
|
|
|
|
33
|
my $W = mod2pi( deg2rad( 35999.373 * $T - 2.47 ) ); |
588
|
11
|
|
|
|
|
34
|
my $delta_lambda = 1 + 0.0334 * cos( $W ) + 0.0007 * cos( 2 * $W ); |
589
|
|
|
|
|
|
|
$self->{debug} |
590
|
11
|
50
|
|
|
|
26
|
and print "Debug - delta lambda = $delta_lambda\n"; |
591
|
11
|
|
|
|
|
16
|
my $S = 0; |
592
|
11
|
|
|
|
|
121
|
foreach my $term ( # Meeus table 27 C |
593
|
|
|
|
|
|
|
[ 485, 324.96, 1934.136 ], |
594
|
|
|
|
|
|
|
[ 203, 337.23, 32964.467 ], |
595
|
|
|
|
|
|
|
[ 199, 342.08, 20.186 ], |
596
|
|
|
|
|
|
|
[ 182, 27.85, 445267.112 ], |
597
|
|
|
|
|
|
|
[ 156, 73.14, 45036.886 ], |
598
|
|
|
|
|
|
|
[ 136, 171.52, 22518.443 ], |
599
|
|
|
|
|
|
|
[ 77, 222.54, 65928.934 ], |
600
|
|
|
|
|
|
|
[ 74, 296.72, 3034.906 ], |
601
|
|
|
|
|
|
|
[ 70, 243.58, 9037.513 ], |
602
|
|
|
|
|
|
|
[ 58, 119.81, 33718.147 ], |
603
|
|
|
|
|
|
|
[ 52, 297.17, 150.678 ], |
604
|
|
|
|
|
|
|
[ 50, 21.02, 2281.226 ], |
605
|
|
|
|
|
|
|
[ 45, 247.54, 29929.562 ], |
606
|
|
|
|
|
|
|
[ 44, 325.15, 31555.956 ], |
607
|
|
|
|
|
|
|
[ 29, 60.93, 4443.417 ], |
608
|
|
|
|
|
|
|
[ 18, 155.12, 67555.328 ], |
609
|
|
|
|
|
|
|
[ 17, 288.79, 4562.452 ], |
610
|
|
|
|
|
|
|
[ 16, 198.04, 62894.029 ], |
611
|
|
|
|
|
|
|
[ 14, 199.76, 31436.921 ], |
612
|
|
|
|
|
|
|
[ 12, 95.39, 14577.848 ], |
613
|
|
|
|
|
|
|
[ 12, 287.11, 31931.756 ], |
614
|
|
|
|
|
|
|
[ 12, 320.81, 34777.259 ], |
615
|
|
|
|
|
|
|
[ 9, 227.73, 1222.114 ], |
616
|
|
|
|
|
|
|
[ 8, 15.45, 16859.074 ], |
617
|
|
|
|
|
|
|
) { |
618
|
264
|
|
|
|
|
556
|
$S += $term->[0] * cos( mod2pi( deg2rad( $term->[1] + $term->[2] |
619
|
|
|
|
|
|
|
* $T ) ) ); |
620
|
|
|
|
|
|
|
} |
621
|
|
|
|
|
|
|
$self->{debug} |
622
|
11
|
50
|
|
|
|
56
|
and print "Debug - S = $S\n"; |
623
|
11
|
|
|
|
|
18
|
my $JDE = 0.00001 * $S / $delta_lambda + $JDE0; |
624
|
|
|
|
|
|
|
$self->{debug} |
625
|
11
|
50
|
|
|
|
23
|
and print "Debug - JDE = $JDE\n"; |
626
|
11
|
|
|
|
|
19
|
my $time = ( $JDE - JD_OF_EPOCH ) * SECSPERDAY; |
627
|
|
|
|
|
|
|
# Note that gmtime() in the following needs the parens because it |
628
|
|
|
|
|
|
|
# might have come from Time::y2038, which appears to take more than |
629
|
|
|
|
|
|
|
# one argument -- even though, as I read it, its prototype is (;$). |
630
|
|
|
|
|
|
|
$self->{debug} |
631
|
11
|
50
|
|
|
|
21
|
and print "Debug - dynamical date is ", scalar gmtime( $time ), "\n"; |
632
|
11
|
|
|
|
|
27
|
return $time - dynamical_delta( $time ); # Not quite right. |
633
|
|
|
|
|
|
|
} |
634
|
|
|
|
|
|
|
|
635
|
|
|
|
|
|
|
=item $sun->set( ... ) |
636
|
|
|
|
|
|
|
|
637
|
|
|
|
|
|
|
This method has been overridden to silently ignore any attempt to set |
638
|
|
|
|
|
|
|
the C<'sun'> attribute. |
639
|
|
|
|
|
|
|
|
640
|
|
|
|
|
|
|
=cut |
641
|
|
|
|
|
|
|
|
642
|
|
|
|
|
|
|
sub set { |
643
|
75
|
|
|
75
|
1
|
215
|
my ( $self, @args ) = @_; |
644
|
75
|
|
|
|
|
200
|
while ( @args ) { |
645
|
188
|
|
|
|
|
472
|
my ( $name, $val ) = splice @args, 0, 2; |
646
|
188
|
50
|
|
|
|
511
|
if ( 'sun' eq $name ) { |
|
|
100
|
|
|
|
|
|
647
|
|
|
|
|
|
|
} elsif ( $attrib{$name} ) { |
648
|
37
|
50
|
|
|
|
104
|
if ( ref $self ) { |
649
|
37
|
|
|
|
|
132
|
$self->{$name} = $val; |
650
|
|
|
|
|
|
|
} else { |
651
|
0
|
|
|
|
|
0
|
$static{$name} = $val; |
652
|
|
|
|
|
|
|
} |
653
|
|
|
|
|
|
|
} else { |
654
|
151
|
|
|
|
|
402
|
$self->SUPER::set( $name, $val ); |
655
|
|
|
|
|
|
|
} |
656
|
|
|
|
|
|
|
} |
657
|
75
|
|
|
|
|
159
|
return $self; |
658
|
|
|
|
|
|
|
} |
659
|
|
|
|
|
|
|
|
660
|
|
|
|
|
|
|
=item $sun->time_set () |
661
|
|
|
|
|
|
|
|
662
|
|
|
|
|
|
|
This method sets coordinates of the object to the coordinates of the |
663
|
|
|
|
|
|
|
Sun at the object's currently-set universal time. The velocity |
664
|
|
|
|
|
|
|
components are arbitrarily set to 0. The 'equinox_dynamical' attribute |
665
|
|
|
|
|
|
|
is set to the object's currently-set dynamical time. |
666
|
|
|
|
|
|
|
|
667
|
|
|
|
|
|
|
Although there's no reason this method can't be called directly, it |
668
|
|
|
|
|
|
|
exists to take advantage of the hook in the B |
669
|
|
|
|
|
|
|
object, to allow the position of the Sun to be computed when the |
670
|
|
|
|
|
|
|
object's time is set. |
671
|
|
|
|
|
|
|
|
672
|
|
|
|
|
|
|
The algorithm comes from Jean Meeus' "Astronomical Algorithms", 2nd |
673
|
|
|
|
|
|
|
Edition, Chapter 25, pages 163ff. |
674
|
|
|
|
|
|
|
|
675
|
|
|
|
|
|
|
=cut |
676
|
|
|
|
|
|
|
|
677
|
|
|
|
|
|
|
# The following constants are used in the time_set() method, |
678
|
|
|
|
|
|
|
# because Meeus' equations are in degrees, I was too lazy |
679
|
|
|
|
|
|
|
# to hand-convert them to radians, but I didn't want to |
680
|
|
|
|
|
|
|
# penalize the user for the conversion every time. |
681
|
|
|
|
|
|
|
|
682
|
16
|
|
|
16
|
|
133
|
use constant SUN_C1_0 => deg2rad (1.914602); |
|
16
|
|
|
|
|
37
|
|
|
16
|
|
|
|
|
76
|
|
683
|
16
|
|
|
16
|
|
106
|
use constant SUN_C1_1 => deg2rad (-0.004817); |
|
16
|
|
|
|
|
36
|
|
|
16
|
|
|
|
|
67
|
|
684
|
16
|
|
|
16
|
|
107
|
use constant SUN_C1_2 => deg2rad (-0.000014); |
|
16
|
|
|
|
|
41
|
|
|
16
|
|
|
|
|
71
|
|
685
|
16
|
|
|
16
|
|
104
|
use constant SUN_C2_0 => deg2rad (0.019993); |
|
16
|
|
|
|
|
43
|
|
|
16
|
|
|
|
|
66
|
|
686
|
16
|
|
|
16
|
|
100
|
use constant SUN_C2_1 => deg2rad (0.000101); |
|
16
|
|
|
|
|
41
|
|
|
16
|
|
|
|
|
62
|
|
687
|
16
|
|
|
16
|
|
98
|
use constant SUN_C3_0 => deg2rad (0.000289); |
|
16
|
|
|
|
|
43
|
|
|
16
|
|
|
|
|
56
|
|
688
|
16
|
|
|
16
|
|
103
|
use constant SUN_LON_2000 => deg2rad (- 0.01397); |
|
16
|
|
|
|
|
30
|
|
|
16
|
|
|
|
|
90
|
|
689
|
|
|
|
|
|
|
|
690
|
|
|
|
|
|
|
sub time_set { |
691
|
11257
|
|
|
11257
|
1
|
16724
|
my $self = shift; |
692
|
11257
|
|
|
|
|
23012
|
my $time = $self->dynamical; |
693
|
|
|
|
|
|
|
|
694
|
|
|
|
|
|
|
# The following algorithm is from Meeus, chapter 25, page, 163 ff. |
695
|
|
|
|
|
|
|
|
696
|
11257
|
|
|
|
|
26697
|
my $T = jcent2000 ($time); # Meeus (25.1) |
697
|
11257
|
|
|
|
|
28045
|
my $L0 = mod2pi(deg2rad((.0003032 * $T + 36000.76983) * $T # Meeus (25.2) |
698
|
|
|
|
|
|
|
+ 280.46646)); |
699
|
11257
|
|
|
|
|
26261
|
my $M = mod2pi(deg2rad(((-.0001537) * $T + 35999.05029) # Meeus (25.3) |
700
|
|
|
|
|
|
|
* $T + 357.52911)); |
701
|
11257
|
|
|
|
|
20817
|
my $e = (-0.0000001267 * $T - 0.000042037) * $T + 0.016708634;# Meeus (25.4) |
702
|
11257
|
|
|
|
|
32573
|
my $C = ((SUN_C1_2 * $T + SUN_C1_1) * $T + SUN_C1_0) * sin ($M) |
703
|
|
|
|
|
|
|
+ (SUN_C2_1 * $T + SUN_C2_0) * sin (2 * $M) |
704
|
|
|
|
|
|
|
+ SUN_C3_0 * sin (3 * $M); |
705
|
11257
|
|
|
|
|
20819
|
my $O = $self->{_sun_geometric_longitude} = $L0 + $C; |
706
|
11257
|
|
|
|
|
24552
|
my $omega = mod2pi (deg2rad (125.04 - 1934.136 * $T)); |
707
|
11257
|
|
|
|
|
26548
|
my $lambda = mod2pi ($O - deg2rad (0.00569 + 0.00478 * sin ($omega))); |
708
|
11257
|
|
|
|
|
19154
|
my $nu = $M + $C; |
709
|
11257
|
|
|
|
|
26109
|
my $R = (1.000_001_018 * (1 - $e * $e)) / (1 + $e * cos ($nu)) |
710
|
|
|
|
|
|
|
* AU; |
711
|
11257
|
50
|
|
|
|
25574
|
$self->{debug} and print <
|
712
|
0
|
|
|
|
|
0
|
Debug sun - @{[strftime '%d-%b-%Y %H:%M:%S', gmtime( $time )]} |
713
|
|
|
|
|
|
|
T = $T |
714
|
0
|
|
|
|
|
0
|
L0 = @{[rad2deg ($L0)]} degrees |
715
|
0
|
|
|
|
|
0
|
M = @{[rad2deg ($M)]} degrees |
716
|
|
|
|
|
|
|
e = $e |
717
|
0
|
|
|
|
|
0
|
C = @{[rad2deg ($C)]} degrees |
718
|
0
|
|
|
|
|
0
|
O = @{[rad2deg ($O)]} degrees |
719
|
0
|
|
|
|
|
0
|
R = @{[$R / AU]} AU |
720
|
0
|
|
|
|
|
0
|
omega = @{[rad2deg ($omega)]} degrees |
721
|
0
|
|
|
|
|
0
|
lambda = @{[rad2deg ($lambda)]} degrees |
722
|
|
|
|
|
|
|
eod |
723
|
|
|
|
|
|
|
|
724
|
11257
|
|
|
|
|
31606
|
$self->ecliptic (0, $lambda, $R); |
725
|
|
|
|
|
|
|
## $self->set (equinox_dynamical => $time); |
726
|
11257
|
|
|
|
|
32926
|
$self->equinox_dynamical ($time); |
727
|
11257
|
|
|
|
|
21470
|
return $self; |
728
|
|
|
|
|
|
|
} |
729
|
|
|
|
|
|
|
|
730
|
|
|
|
|
|
|
# The Sun is normally positioned in inertial coordinates. |
731
|
|
|
|
|
|
|
|
732
|
37
|
|
|
37
|
|
157
|
sub __initial_inertial { return 1 } |
733
|
|
|
|
|
|
|
|
734
|
|
|
|
|
|
|
1; |
735
|
|
|
|
|
|
|
|
736
|
|
|
|
|
|
|
=back |
737
|
|
|
|
|
|
|
|
738
|
|
|
|
|
|
|
=head2 Attributes |
739
|
|
|
|
|
|
|
|
740
|
|
|
|
|
|
|
This class has the following public attributes. The description gives |
741
|
|
|
|
|
|
|
the data type. |
742
|
|
|
|
|
|
|
|
743
|
|
|
|
|
|
|
=over |
744
|
|
|
|
|
|
|
|
745
|
|
|
|
|
|
|
=item iterate_for_quarters (Boolean) |
746
|
|
|
|
|
|
|
|
747
|
|
|
|
|
|
|
If this attribute is true, the C method uses the old |
748
|
|
|
|
|
|
|
(pre-0.088_01) algorithm. |
749
|
|
|
|
|
|
|
|
750
|
|
|
|
|
|
|
If this attribute is false, the new algorithm is used. |
751
|
|
|
|
|
|
|
|
752
|
|
|
|
|
|
|
The default is C, i.e. false, because I believe the new algorithm |
753
|
|
|
|
|
|
|
to be more accurate for reasonably-current times. |
754
|
|
|
|
|
|
|
|
755
|
|
|
|
|
|
|
This attribute is new with version 0.088_01. |
756
|
|
|
|
|
|
|
|
757
|
|
|
|
|
|
|
=back |
758
|
|
|
|
|
|
|
|
759
|
|
|
|
|
|
|
=head2 Historical Calculations |
760
|
|
|
|
|
|
|
|
761
|
|
|
|
|
|
|
This class was written for the purpose of calculating whether the Sun |
762
|
|
|
|
|
|
|
was shining on a given point on the Earth (or in space) at a given time |
763
|
|
|
|
|
|
|
in or reasonably close to the present. I can not say how accurate it is |
764
|
|
|
|
|
|
|
at times far from the present. Those interested in such calculations may |
765
|
|
|
|
|
|
|
want to consider using |
766
|
|
|
|
|
|
|
L |
767
|
|
|
|
|
|
|
instead. |
768
|
|
|
|
|
|
|
|
769
|
|
|
|
|
|
|
Historical calculations will need to use a 64-bit Perl, or at least one |
770
|
|
|
|
|
|
|
with 64-bit integers, to represent times more than about 38 years from |
771
|
|
|
|
|
|
|
the system epoch. |
772
|
|
|
|
|
|
|
|
773
|
|
|
|
|
|
|
In addition, you will need to be careful how you do input and output |
774
|
|
|
|
|
|
|
conversions. |
775
|
|
|
|
|
|
|
|
776
|
|
|
|
|
|
|
L is a core module and an obvious choice, but |
777
|
|
|
|
|
|
|
it only does Gregorian dates. Historical calculations prior to 1587 |
778
|
|
|
|
|
|
|
typically use the Julian calendar. For this you will need to go to |
779
|
|
|
|
|
|
|
something like L, |
780
|
|
|
|
|
|
|
or L which |
781
|
|
|
|
|
|
|
does either Julian or Gregorian as needed. |
782
|
|
|
|
|
|
|
|
783
|
|
|
|
|
|
|
Should you decide to use L, you should be aware |
784
|
|
|
|
|
|
|
that its C and C interpret the year argument |
785
|
|
|
|
|
|
|
strangely: years in the range C<0 - 999> inclusive are not interpreted |
786
|
|
|
|
|
|
|
as Gregorian years, though years outside that range are so interpreted. |
787
|
|
|
|
|
|
|
Beginning with version 1.27 (released July 9 2018), additional |
788
|
|
|
|
|
|
|
subroutines C and C were added. |
789
|
|
|
|
|
|
|
These always interpret the year argument as a Gregorian year. |
790
|
|
|
|
|
|
|
|
791
|
|
|
|
|
|
|
=head1 ACKNOWLEDGMENTS |
792
|
|
|
|
|
|
|
|
793
|
|
|
|
|
|
|
The author wishes to acknowledge Jean Meeus, whose book "Astronomical |
794
|
|
|
|
|
|
|
Algorithms" (second edition) formed the basis for this module. |
795
|
|
|
|
|
|
|
|
796
|
|
|
|
|
|
|
=head1 SEE ALSO |
797
|
|
|
|
|
|
|
|
798
|
|
|
|
|
|
|
The L |
799
|
|
|
|
|
|
|
documentation for a discussion of how the pieces/parts of this |
800
|
|
|
|
|
|
|
distribution go together and how to use them. |
801
|
|
|
|
|
|
|
|
802
|
|
|
|
|
|
|
L by Brett Hamilton, which contains a |
803
|
|
|
|
|
|
|
function-based module to compute the current phase, distance and angular |
804
|
|
|
|
|
|
|
diameter of the Moon, as well as the angular diameter and distance of |
805
|
|
|
|
|
|
|
the Sun. |
806
|
|
|
|
|
|
|
|
807
|
|
|
|
|
|
|
L by Ron Hill and Jean Forget, which |
808
|
|
|
|
|
|
|
contains a function-based module to compute sunrise and sunset for the |
809
|
|
|
|
|
|
|
given day and location. |
810
|
|
|
|
|
|
|
|
811
|
|
|
|
|
|
|
L by Rob Fugina, which provides |
812
|
|
|
|
|
|
|
functionality similar to B. |
813
|
|
|
|
|
|
|
|
814
|
|
|
|
|
|
|
=head1 SUPPORT |
815
|
|
|
|
|
|
|
|
816
|
|
|
|
|
|
|
Support is by the author. Please file bug reports at |
817
|
|
|
|
|
|
|
L, |
818
|
|
|
|
|
|
|
L, or in |
819
|
|
|
|
|
|
|
electronic mail to the author. |
820
|
|
|
|
|
|
|
|
821
|
|
|
|
|
|
|
=head1 AUTHOR |
822
|
|
|
|
|
|
|
|
823
|
|
|
|
|
|
|
Thomas R. Wyant, III (F) |
824
|
|
|
|
|
|
|
|
825
|
|
|
|
|
|
|
=head1 COPYRIGHT AND LICENSE |
826
|
|
|
|
|
|
|
|
827
|
|
|
|
|
|
|
Copyright (C) 2005-2023 by Thomas R. Wyant, III |
828
|
|
|
|
|
|
|
|
829
|
|
|
|
|
|
|
This program is free software; you can redistribute it and/or modify it |
830
|
|
|
|
|
|
|
under the same terms as Perl 5.10.0. For more details, see the full text |
831
|
|
|
|
|
|
|
of the licenses in the directory LICENSES. |
832
|
|
|
|
|
|
|
|
833
|
|
|
|
|
|
|
This program is distributed in the hope that it will be useful, but |
834
|
|
|
|
|
|
|
without any warranty; without even the implied warranty of |
835
|
|
|
|
|
|
|
merchantability or fitness for a particular purpose. |
836
|
|
|
|
|
|
|
|
837
|
|
|
|
|
|
|
=cut |
838
|
|
|
|
|
|
|
|
839
|
|
|
|
|
|
|
# ex: set textwidth=72 : |