line |
stmt |
bran |
cond |
sub |
pod |
time |
code |
1
|
|
|
|
|
|
|
package DateTime::Astro::Sunrise;
|
2
|
|
|
|
|
|
|
|
3
|
2
|
|
|
2
|
|
1581495
|
use strict;
|
|
2
|
|
|
|
|
7
|
|
|
2
|
|
|
|
|
103
|
|
4
|
|
|
|
|
|
|
require Exporter;
|
5
|
2
|
|
|
2
|
|
11
|
use POSIX;
|
|
2
|
|
|
|
|
4
|
|
|
2
|
|
|
|
|
23
|
|
6
|
2
|
|
|
2
|
|
10051
|
use Math::Trig;
|
|
2
|
|
|
|
|
84906
|
|
|
2
|
|
|
|
|
1397
|
|
7
|
2
|
|
|
2
|
|
31
|
use Carp;
|
|
2
|
|
|
|
|
5
|
|
|
2
|
|
|
|
|
993
|
|
8
|
2
|
|
|
2
|
|
17
|
use DateTime;
|
|
2
|
|
|
|
|
4
|
|
|
2
|
|
|
|
|
90
|
|
9
|
2
|
|
|
2
|
|
12
|
use vars qw( $VERSION $RADEG $DEGRAD @ISA );
|
|
2
|
|
|
|
|
6
|
|
|
2
|
|
|
|
|
7863
|
|
10
|
|
|
|
|
|
|
@ISA = qw( Exporter );
|
11
|
|
|
|
|
|
|
$VERSION = qw($Revision: 0.01_01 $) [1];
|
12
|
|
|
|
|
|
|
$RADEG = ( 180 / pi );
|
13
|
|
|
|
|
|
|
$DEGRAD = ( pi / 180 );
|
14
|
|
|
|
|
|
|
my $INV360 = ( 1.0 / 360.0 );
|
15
|
|
|
|
|
|
|
|
16
|
|
|
|
|
|
|
my $upper_limb = '1';
|
17
|
|
|
|
|
|
|
|
18
|
|
|
|
|
|
|
sub new {
|
19
|
1
|
|
|
1
|
1
|
724
|
my $class = shift;
|
20
|
1
|
|
|
|
|
3
|
my %args;
|
21
|
1
|
|
|
|
|
3
|
$args{LON} = shift;
|
22
|
1
|
|
|
|
|
3
|
$args{LAT} = shift;
|
23
|
1
|
|
|
|
|
2
|
$args{ALT} = shift;
|
24
|
1
|
|
|
|
|
2
|
$args{ITER} = shift;
|
25
|
1
|
50
|
|
|
|
4
|
unless ( $args{LON} ) {
|
26
|
0
|
|
|
|
|
0
|
croak "You need to have a longitude\n";
|
27
|
|
|
|
|
|
|
}
|
28
|
|
|
|
|
|
|
|
29
|
1
|
50
|
|
|
|
4
|
unless ( $args{LAT} ) {
|
30
|
0
|
|
|
|
|
0
|
croak "You need to have a latitude\n";
|
31
|
|
|
|
|
|
|
}
|
32
|
|
|
|
|
|
|
|
33
|
1
|
|
|
|
|
4
|
return bless \%args, $class;
|
34
|
|
|
|
|
|
|
|
35
|
|
|
|
|
|
|
}
|
36
|
|
|
|
|
|
|
|
37
|
|
|
|
|
|
|
sub sunrise {
|
38
|
1
|
|
|
1
|
0
|
6
|
my ( $self, $dt ) = @_;
|
39
|
|
|
|
|
|
|
|
40
|
1
|
|
|
|
|
6
|
my ( $year, $month, $day ) = ( $dt->year, $dt->month, $dt->day );
|
41
|
1
|
|
50
|
|
|
31
|
my $altit = $self->{ALT} || -0.833;
|
42
|
1
|
50
|
|
|
|
5
|
my $iteration = defined( $self->{ITER} ) ? $self->{ITER} : 0;
|
43
|
|
|
|
|
|
|
|
44
|
1
|
50
|
|
|
|
4
|
if ($iteration) {
|
45
|
|
|
|
|
|
|
|
46
|
|
|
|
|
|
|
# This is the initial start
|
47
|
|
|
|
|
|
|
|
48
|
0
|
|
|
|
|
0
|
my $d =
|
49
|
|
|
|
|
|
|
days_since_2000_Jan_0( $year, $month, $day ) + 0.5 - $self->{LON} /
|
50
|
|
|
|
|
|
|
360.0;
|
51
|
0
|
|
|
|
|
0
|
my ( $tmp_rise_1, $tmp_set_1 ) =
|
52
|
|
|
|
|
|
|
sun_rise_set( $d, $self->{LON}, $self->{LAT}, $altit,
|
53
|
|
|
|
|
|
|
15.04107 );
|
54
|
|
|
|
|
|
|
|
55
|
|
|
|
|
|
|
# Now we have the initial rise/set times next recompute d using the exact moment
|
56
|
|
|
|
|
|
|
# recompute sunrise
|
57
|
|
|
|
|
|
|
|
58
|
0
|
|
|
|
|
0
|
my $tmp_rise_2 = 9;
|
59
|
0
|
|
|
|
|
0
|
my $tmp_rise_3 = 0;
|
60
|
0
|
|
|
|
|
0
|
until ( equal( $tmp_rise_2, $tmp_rise_3, 8 ) ) {
|
61
|
|
|
|
|
|
|
|
62
|
0
|
|
|
|
|
0
|
my $d_sunrise_1 = $d + $tmp_rise_1 / 24.0;
|
63
|
0
|
|
|
|
|
0
|
( $tmp_rise_2, undef ) = sun_rise_set(
|
64
|
|
|
|
|
|
|
$d_sunrise_1, $self->{LON}, $self->{LAT}, $altit,
|
65
|
|
|
|
|
|
|
15.04107
|
66
|
|
|
|
|
|
|
);
|
67
|
0
|
|
|
|
|
0
|
$tmp_rise_1 = $tmp_rise_3;
|
68
|
0
|
|
|
|
|
0
|
my $d_sunrise_2 = $d + $tmp_rise_2 / 24.0;
|
69
|
0
|
|
|
|
|
0
|
( $tmp_rise_3, undef ) = sun_rise_set(
|
70
|
|
|
|
|
|
|
$d_sunrise_2, $self->{LON}, $self->{LAT}, $altit,
|
71
|
|
|
|
|
|
|
15.04107
|
72
|
|
|
|
|
|
|
);
|
73
|
|
|
|
|
|
|
|
74
|
|
|
|
|
|
|
#print "tmp_rise2 is: $tmp_rise_2 tmp_rise_3 is:$tmp_rise_3\n";
|
75
|
|
|
|
|
|
|
|
76
|
|
|
|
|
|
|
}
|
77
|
|
|
|
|
|
|
|
78
|
|
|
|
|
|
|
###################################################################################
|
79
|
|
|
|
|
|
|
# end sunrise
|
80
|
|
|
|
|
|
|
###################################################################################
|
81
|
|
|
|
|
|
|
|
82
|
0
|
|
|
|
|
0
|
my $tmp_set_2 = 9;
|
83
|
0
|
|
|
|
|
0
|
my $tmp_set_3 = 0;
|
84
|
|
|
|
|
|
|
|
85
|
0
|
|
|
|
|
0
|
until ( equal( $tmp_set_2, $tmp_set_3, 8 ) ) {
|
86
|
|
|
|
|
|
|
|
87
|
0
|
|
|
|
|
0
|
my $d_sunset_1 = $d + $tmp_set_1 / 24.0;
|
88
|
0
|
|
|
|
|
0
|
( undef, $tmp_set_2 ) = sun_rise_set(
|
89
|
|
|
|
|
|
|
$d_sunset_1, $self->{LON}, $self->{LAT}, $altit,
|
90
|
|
|
|
|
|
|
15.04107
|
91
|
|
|
|
|
|
|
);
|
92
|
0
|
|
|
|
|
0
|
$tmp_set_1 = $tmp_set_3;
|
93
|
0
|
|
|
|
|
0
|
my $d_sunset_2 = $d + $tmp_set_2 / 24.0;
|
94
|
0
|
|
|
|
|
0
|
( undef, $tmp_set_3 ) = sun_rise_set(
|
95
|
|
|
|
|
|
|
$d_sunset_2, $self->{LON}, $self->{LAT}, $altit,
|
96
|
|
|
|
|
|
|
15.04107
|
97
|
|
|
|
|
|
|
);
|
98
|
|
|
|
|
|
|
|
99
|
|
|
|
|
|
|
#print "tmp_set_1 is: $tmp_set_1 tmp_set_3 is:$tmp_set_3\n";
|
100
|
|
|
|
|
|
|
|
101
|
|
|
|
|
|
|
}
|
102
|
|
|
|
|
|
|
|
103
|
0
|
|
|
|
|
0
|
my ( $hour_rise, $min_rise, $hour_set, $min_set ) = convert_hour($tmp_rise_3,$tmp_set_3);
|
104
|
0
|
|
|
|
|
0
|
my $rise_time = DateTime->new(
|
105
|
|
|
|
|
|
|
year => $dt->year,
|
106
|
|
|
|
|
|
|
month => $dt->month,
|
107
|
|
|
|
|
|
|
day => $dt->day,
|
108
|
|
|
|
|
|
|
hour => $hour_rise,
|
109
|
|
|
|
|
|
|
minute => $min_rise,
|
110
|
|
|
|
|
|
|
time_zone => 'UTC'
|
111
|
|
|
|
|
|
|
);
|
112
|
0
|
|
|
|
|
0
|
my $set_time = DateTime->new(
|
113
|
|
|
|
|
|
|
year => $dt->year,
|
114
|
|
|
|
|
|
|
month => $dt->month,
|
115
|
|
|
|
|
|
|
day => $dt->day,
|
116
|
|
|
|
|
|
|
hour => $hour_set,
|
117
|
|
|
|
|
|
|
minute => $min_set,
|
118
|
|
|
|
|
|
|
time_zone => 'UTC'
|
119
|
|
|
|
|
|
|
);
|
120
|
0
|
|
|
|
|
0
|
return ($rise_time, $set_time);
|
121
|
|
|
|
|
|
|
}
|
122
|
|
|
|
|
|
|
else {
|
123
|
1
|
|
|
|
|
5
|
my $d =
|
124
|
|
|
|
|
|
|
days_since_2000_Jan_0( $year, $month, $day ) + 0.5 - $self->{LON} /
|
125
|
|
|
|
|
|
|
360.0;
|
126
|
1
|
|
|
|
|
5
|
my ( $h1, $h2 ) =
|
127
|
|
|
|
|
|
|
sun_rise_set( $d, $self->{LON}, $self->{LAT}, $altit, 15.0 );
|
128
|
1
|
|
|
|
|
5
|
my ( $hour_rise, $min_rise, $hour_set, $min_set ) =
|
129
|
|
|
|
|
|
|
convert_hour( $h1, $h2 );
|
130
|
|
|
|
|
|
|
|
131
|
1
|
|
|
|
|
5
|
my $rise_time = DateTime->new(
|
132
|
|
|
|
|
|
|
year => $dt->year,
|
133
|
|
|
|
|
|
|
month => $dt->month,
|
134
|
|
|
|
|
|
|
day => $dt->day,
|
135
|
|
|
|
|
|
|
hour => $hour_rise,
|
136
|
|
|
|
|
|
|
minute => $min_rise,
|
137
|
|
|
|
|
|
|
time_zone => 'UTC'
|
138
|
|
|
|
|
|
|
);
|
139
|
1
|
|
|
|
|
289
|
my $set_time = DateTime->new(
|
140
|
|
|
|
|
|
|
year => $dt->year,
|
141
|
|
|
|
|
|
|
month => $dt->month,
|
142
|
|
|
|
|
|
|
day => $dt->day,
|
143
|
|
|
|
|
|
|
hour => $hour_set,
|
144
|
|
|
|
|
|
|
minute => $min_set,
|
145
|
|
|
|
|
|
|
time_zone => 'UTC'
|
146
|
|
|
|
|
|
|
);
|
147
|
0
|
|
|
|
|
0
|
return ($rise_time, $set_time);
|
148
|
|
|
|
|
|
|
}
|
149
|
|
|
|
|
|
|
|
150
|
|
|
|
|
|
|
}
|
151
|
|
|
|
|
|
|
|
152
|
|
|
|
|
|
|
sub sun_rise_set {
|
153
|
1
|
|
|
1
|
0
|
3
|
my ( $d, $lon, $lat, $altit, $h ) = @_;
|
154
|
|
|
|
|
|
|
|
155
|
|
|
|
|
|
|
|
156
|
1
|
|
|
|
|
6
|
my $sidtime = revolution( GMST0($d) + 180.0 + $lon );
|
157
|
|
|
|
|
|
|
|
158
|
1
|
|
|
|
|
4
|
my ( $sRA, $sdec ) = sun_RA_dec($d);
|
159
|
1
|
|
|
|
|
4
|
my $tsouth = 12.0 - rev180( $sidtime - $sRA ) /$h ;
|
160
|
1
|
|
|
|
|
2
|
my $sradius = 0.2666 / $sRA;
|
161
|
|
|
|
|
|
|
|
162
|
1
|
50
|
|
|
|
4
|
if ($upper_limb) {
|
163
|
1
|
|
|
|
|
2
|
$altit -= $sradius;
|
164
|
|
|
|
|
|
|
}
|
165
|
|
|
|
|
|
|
|
166
|
|
|
|
|
|
|
# Compute the diurnal arc that the Sun traverses to reach
|
167
|
|
|
|
|
|
|
# the specified altitude altit:
|
168
|
|
|
|
|
|
|
|
169
|
1
|
|
|
|
|
3
|
my $cost =
|
170
|
|
|
|
|
|
|
( sind($altit) - sind($lat) * sind($sdec) ) /
|
171
|
|
|
|
|
|
|
( cosd($lat) * cosd($sdec) );
|
172
|
|
|
|
|
|
|
|
173
|
1
|
|
|
|
|
3
|
my $t;
|
174
|
1
|
50
|
|
|
|
5
|
if ( $cost >= 1.0 ) {
|
|
|
50
|
|
|
|
|
|
175
|
0
|
|
|
|
|
0
|
carp "Sun never rises!!\n";
|
176
|
0
|
|
|
|
|
0
|
$t = 0.0; # Sun always below altit
|
177
|
|
|
|
|
|
|
}
|
178
|
|
|
|
|
|
|
elsif ( $cost <= -1.0 ) {
|
179
|
0
|
|
|
|
|
0
|
carp "Sun never sets!!\n";
|
180
|
0
|
|
|
|
|
0
|
$t = 12.0; # Sun always above altit
|
181
|
|
|
|
|
|
|
}
|
182
|
|
|
|
|
|
|
else {
|
183
|
1
|
|
|
|
|
3
|
$t = acosd($cost) / 15.0; # The diurnal arc, hours
|
184
|
|
|
|
|
|
|
}
|
185
|
|
|
|
|
|
|
|
186
|
|
|
|
|
|
|
# Store rise and set times - in hours UT
|
187
|
|
|
|
|
|
|
|
188
|
1
|
|
|
|
|
11
|
my $hour_rise_ut = $tsouth - $t;
|
189
|
1
|
|
|
|
|
1
|
my $hour_set_ut = $tsouth + $t;
|
190
|
1
|
|
|
|
|
3
|
return ( $hour_rise_ut, $hour_set_ut );
|
191
|
|
|
|
|
|
|
|
192
|
|
|
|
|
|
|
|
193
|
|
|
|
|
|
|
}
|
194
|
|
|
|
|
|
|
|
195
|
|
|
|
|
|
|
#########################################################################################################
|
196
|
|
|
|
|
|
|
sub GMST0 {
|
197
|
|
|
|
|
|
|
|
198
|
|
|
|
|
|
|
#
|
199
|
|
|
|
|
|
|
#
|
200
|
|
|
|
|
|
|
# FUNCTIONAL SEQUENCE for GMST0
|
201
|
|
|
|
|
|
|
#
|
202
|
|
|
|
|
|
|
# _GIVEN
|
203
|
|
|
|
|
|
|
# Day number
|
204
|
|
|
|
|
|
|
#
|
205
|
|
|
|
|
|
|
# _THEN
|
206
|
|
|
|
|
|
|
#
|
207
|
|
|
|
|
|
|
# computes GMST0, the Greenwich Mean Sidereal Time
|
208
|
|
|
|
|
|
|
# at 0h UT (i.e. the sidereal time at the Greenwhich meridian at
|
209
|
|
|
|
|
|
|
# 0h UT). GMST is then the sidereal time at Greenwich at any
|
210
|
|
|
|
|
|
|
# time of the day..
|
211
|
|
|
|
|
|
|
#
|
212
|
|
|
|
|
|
|
#
|
213
|
|
|
|
|
|
|
# _RETURN
|
214
|
|
|
|
|
|
|
#
|
215
|
|
|
|
|
|
|
# Sidtime
|
216
|
|
|
|
|
|
|
#
|
217
|
1
|
|
|
1
|
0
|
2
|
my ($d) = @_;
|
218
|
|
|
|
|
|
|
|
219
|
1
|
|
|
|
|
6
|
my $sidtim0 =
|
220
|
|
|
|
|
|
|
revolution( ( 180.0 + 356.0470 + 282.9404 ) +
|
221
|
|
|
|
|
|
|
( 0.9856002585 + 4.70935E-5 ) * $d );
|
222
|
1
|
|
|
|
|
5
|
return $sidtim0;
|
223
|
|
|
|
|
|
|
|
224
|
|
|
|
|
|
|
}
|
225
|
|
|
|
|
|
|
|
226
|
|
|
|
|
|
|
sub sunpos {
|
227
|
|
|
|
|
|
|
|
228
|
|
|
|
|
|
|
#
|
229
|
|
|
|
|
|
|
#
|
230
|
|
|
|
|
|
|
# FUNCTIONAL SEQUENCE for sunpos
|
231
|
|
|
|
|
|
|
#
|
232
|
|
|
|
|
|
|
# _GIVEN
|
233
|
|
|
|
|
|
|
# day number
|
234
|
|
|
|
|
|
|
#
|
235
|
|
|
|
|
|
|
# _THEN
|
236
|
|
|
|
|
|
|
#
|
237
|
|
|
|
|
|
|
# Computes the Sun's ecliptic longitude and distance */
|
238
|
|
|
|
|
|
|
# at an instant given in d, number of days since */
|
239
|
|
|
|
|
|
|
# 2000 Jan 0.0.
|
240
|
|
|
|
|
|
|
#
|
241
|
|
|
|
|
|
|
#
|
242
|
|
|
|
|
|
|
# _RETURN
|
243
|
|
|
|
|
|
|
#
|
244
|
|
|
|
|
|
|
# ecliptic longitude and distance
|
245
|
|
|
|
|
|
|
# ie. $True_solar_longitude, $Solar_distance
|
246
|
|
|
|
|
|
|
#
|
247
|
1
|
|
|
1
|
0
|
25
|
my ($d) = @_;
|
248
|
|
|
|
|
|
|
|
249
|
|
|
|
|
|
|
# Mean anomaly of the Sun
|
250
|
|
|
|
|
|
|
# Mean longitude of perihelion
|
251
|
|
|
|
|
|
|
# Note: Sun's mean longitude = M + w
|
252
|
|
|
|
|
|
|
# Eccentricity of Earth's orbit
|
253
|
|
|
|
|
|
|
# Eccentric anomaly
|
254
|
|
|
|
|
|
|
# x, y coordinates in orbit
|
255
|
|
|
|
|
|
|
# True anomaly
|
256
|
|
|
|
|
|
|
|
257
|
|
|
|
|
|
|
# Compute mean elements
|
258
|
1
|
|
|
|
|
3
|
my $Mean_anomaly_of_sun = revolution( 356.0470 + 0.9856002585 * $d );
|
259
|
1
|
|
|
|
|
4
|
my $Mean_longitude_of_perihelion = 282.9404 + 4.70935E-5 * $d;
|
260
|
1
|
|
|
|
|
2
|
my $Eccentricity_of_Earth_orbit = 0.016709 - 1.151E-9 * $d;
|
261
|
|
|
|
|
|
|
|
262
|
|
|
|
|
|
|
# Compute true longitude and radius vector
|
263
|
1
|
|
|
|
|
6
|
my $Eccentric_anomaly =
|
264
|
|
|
|
|
|
|
$Mean_anomaly_of_sun + $Eccentricity_of_Earth_orbit * $RADEG *
|
265
|
|
|
|
|
|
|
sind($Mean_anomaly_of_sun) *
|
266
|
|
|
|
|
|
|
( 1.0 + $Eccentricity_of_Earth_orbit * cosd($Mean_anomaly_of_sun) );
|
267
|
|
|
|
|
|
|
|
268
|
1
|
|
|
|
|
3
|
my $x = cosd($Eccentric_anomaly) - $Eccentricity_of_Earth_orbit;
|
269
|
|
|
|
|
|
|
|
270
|
1
|
|
|
|
|
6
|
my $y =
|
271
|
|
|
|
|
|
|
sqrt( 1.0 - $Eccentricity_of_Earth_orbit * $Eccentricity_of_Earth_orbit )
|
272
|
|
|
|
|
|
|
* sind($Eccentric_anomaly);
|
273
|
|
|
|
|
|
|
|
274
|
1
|
|
|
|
|
3
|
my $Solar_distance = sqrt( $x * $x + $y * $y ); # Solar distance
|
275
|
1
|
|
|
|
|
4
|
my $True_anomaly = atan2d( $y, $x ); # True anomaly
|
276
|
|
|
|
|
|
|
|
277
|
1
|
|
|
|
|
3
|
my $True_solar_longitude =
|
278
|
|
|
|
|
|
|
$True_anomaly + $Mean_longitude_of_perihelion; # True solar longitude
|
279
|
|
|
|
|
|
|
|
280
|
1
|
50
|
|
|
|
4
|
if ( $True_solar_longitude >= 360.0 ) {
|
281
|
1
|
|
|
|
|
2
|
$True_solar_longitude -= 360.0; # Make it 0..360 degrees
|
282
|
|
|
|
|
|
|
}
|
283
|
|
|
|
|
|
|
|
284
|
1
|
|
|
|
|
4
|
return ( $Solar_distance, $True_solar_longitude );
|
285
|
|
|
|
|
|
|
}
|
286
|
|
|
|
|
|
|
|
287
|
|
|
|
|
|
|
sub sun_RA_dec {
|
288
|
|
|
|
|
|
|
|
289
|
|
|
|
|
|
|
#
|
290
|
|
|
|
|
|
|
#
|
291
|
|
|
|
|
|
|
# FUNCTIONAL SEQUENCE for sun_RA_dec
|
292
|
|
|
|
|
|
|
#
|
293
|
|
|
|
|
|
|
# _GIVEN
|
294
|
|
|
|
|
|
|
# day number, $r and $lon (from sunpos)
|
295
|
|
|
|
|
|
|
#
|
296
|
|
|
|
|
|
|
# _THEN
|
297
|
|
|
|
|
|
|
#
|
298
|
|
|
|
|
|
|
# compute RA and dec
|
299
|
|
|
|
|
|
|
#
|
300
|
|
|
|
|
|
|
#
|
301
|
|
|
|
|
|
|
# _RETURN
|
302
|
|
|
|
|
|
|
#
|
303
|
|
|
|
|
|
|
# Sun's Right Ascension (RA) and Declination (dec)
|
304
|
|
|
|
|
|
|
#
|
305
|
|
|
|
|
|
|
#
|
306
|
1
|
|
|
1
|
0
|
2
|
my ($d) = @_;
|
307
|
|
|
|
|
|
|
|
308
|
|
|
|
|
|
|
# Compute Sun's ecliptical coordinates
|
309
|
1
|
|
|
|
|
19
|
my ( $r, $lon ) = sunpos($d);
|
310
|
|
|
|
|
|
|
|
311
|
|
|
|
|
|
|
# Compute ecliptic rectangular coordinates (z=0)
|
312
|
1
|
|
|
|
|
2
|
my $x = $r * cosd($lon);
|
313
|
1
|
|
|
|
|
4
|
my $y = $r * sind($lon);
|
314
|
|
|
|
|
|
|
|
315
|
|
|
|
|
|
|
# Compute obliquity of ecliptic (inclination of Earth's axis)
|
316
|
1
|
|
|
|
|
3
|
my $obl_ecl = 23.4393 - 3.563E-7 * $d;
|
317
|
|
|
|
|
|
|
|
318
|
|
|
|
|
|
|
# Convert to equatorial rectangular coordinates - x is unchanged
|
319
|
1
|
|
|
|
|
8
|
my $z = $y * sind($obl_ecl);
|
320
|
1
|
|
|
|
|
3
|
$y = $y * cosd($obl_ecl);
|
321
|
|
|
|
|
|
|
|
322
|
|
|
|
|
|
|
# Convert to spherical coordinates
|
323
|
1
|
|
|
|
|
2
|
my $RA = atan2d( $y, $x );
|
324
|
1
|
|
|
|
|
5
|
my $dec = atan2d( $z, sqrt( $x * $x + $y * $y ) );
|
325
|
|
|
|
|
|
|
|
326
|
1
|
|
|
|
|
3
|
return ( $RA, $dec );
|
327
|
|
|
|
|
|
|
|
328
|
|
|
|
|
|
|
} # sun_RA_dec
|
329
|
|
|
|
|
|
|
|
330
|
|
|
|
|
|
|
sub days_since_2000_Jan_0 {
|
331
|
|
|
|
|
|
|
|
332
|
|
|
|
|
|
|
#
|
333
|
|
|
|
|
|
|
#
|
334
|
|
|
|
|
|
|
# FUNCTIONAL SEQUENCE for days_since_2000_Jan_0
|
335
|
|
|
|
|
|
|
#
|
336
|
|
|
|
|
|
|
# _GIVEN
|
337
|
|
|
|
|
|
|
# year, month, day
|
338
|
|
|
|
|
|
|
#
|
339
|
|
|
|
|
|
|
# _THEN
|
340
|
|
|
|
|
|
|
#
|
341
|
|
|
|
|
|
|
# process the year month and day (counted in days)
|
342
|
|
|
|
|
|
|
# Day 0.0 is at Jan 1 2000 0.0 UT
|
343
|
|
|
|
|
|
|
# Note that ALL divisions here should be INTEGER divisions
|
344
|
|
|
|
|
|
|
#
|
345
|
|
|
|
|
|
|
# _RETURN
|
346
|
|
|
|
|
|
|
#
|
347
|
|
|
|
|
|
|
# day number
|
348
|
|
|
|
|
|
|
#
|
349
|
2
|
|
|
2
|
|
29
|
use integer;
|
|
2
|
|
|
|
|
6
|
|
|
2
|
|
|
|
|
24
|
|
350
|
1
|
|
|
1
|
0
|
2
|
my ( $year, $month, $day ) = @_;
|
351
|
|
|
|
|
|
|
|
352
|
1
|
|
|
|
|
5
|
my $d =
|
353
|
|
|
|
|
|
|
( 367 * ($year) -
|
354
|
|
|
|
|
|
|
int( ( 7 * ( ($year) + ( ( ($month) + 9 ) / 12 ) ) ) / 4 ) +
|
355
|
|
|
|
|
|
|
int( ( 275 * ($month) ) / 9 ) + ($day) - 730530 );
|
356
|
|
|
|
|
|
|
|
357
|
1
|
|
|
|
|
5
|
return $d;
|
358
|
|
|
|
|
|
|
|
359
|
|
|
|
|
|
|
}
|
360
|
|
|
|
|
|
|
|
361
|
|
|
|
|
|
|
sub sind {
|
362
|
7
|
|
|
7
|
0
|
40
|
sin( ( $_[0] ) * $DEGRAD );
|
363
|
|
|
|
|
|
|
}
|
364
|
|
|
|
|
|
|
|
365
|
|
|
|
|
|
|
sub cosd {
|
366
|
6
|
|
|
6
|
0
|
26
|
cos( ( $_[0] ) * $DEGRAD );
|
367
|
|
|
|
|
|
|
}
|
368
|
|
|
|
|
|
|
|
369
|
|
|
|
|
|
|
sub tand {
|
370
|
0
|
|
|
0
|
0
|
0
|
tan( ( $_[0] ) * $DEGRAD );
|
371
|
|
|
|
|
|
|
}
|
372
|
|
|
|
|
|
|
|
373
|
|
|
|
|
|
|
sub atand {
|
374
|
0
|
|
|
0
|
0
|
0
|
( $RADEG * atan( $_[0] ) );
|
375
|
|
|
|
|
|
|
}
|
376
|
|
|
|
|
|
|
|
377
|
|
|
|
|
|
|
sub asind {
|
378
|
0
|
|
|
0
|
0
|
0
|
( $RADEG * asin( $_[0] ) );
|
379
|
|
|
|
|
|
|
}
|
380
|
|
|
|
|
|
|
|
381
|
|
|
|
|
|
|
sub acosd {
|
382
|
1
|
|
|
1
|
0
|
7
|
( $RADEG * acos( $_[0] ) );
|
383
|
|
|
|
|
|
|
}
|
384
|
|
|
|
|
|
|
|
385
|
|
|
|
|
|
|
sub atan2d {
|
386
|
3
|
|
|
3
|
0
|
17
|
( $RADEG * atan2( $_[0], $_[1] ) );
|
387
|
|
|
|
|
|
|
}
|
388
|
|
|
|
|
|
|
|
389
|
|
|
|
|
|
|
sub revolution {
|
390
|
|
|
|
|
|
|
|
391
|
|
|
|
|
|
|
#
|
392
|
|
|
|
|
|
|
#
|
393
|
|
|
|
|
|
|
# FUNCTIONAL SEQUENCE for revolution
|
394
|
|
|
|
|
|
|
#
|
395
|
|
|
|
|
|
|
# _GIVEN
|
396
|
|
|
|
|
|
|
# any angle
|
397
|
|
|
|
|
|
|
#
|
398
|
|
|
|
|
|
|
# _THEN
|
399
|
|
|
|
|
|
|
#
|
400
|
|
|
|
|
|
|
# reduces any angle to within the first revolution
|
401
|
|
|
|
|
|
|
# by subtracting or adding even multiples of 360.0
|
402
|
|
|
|
|
|
|
#
|
403
|
|
|
|
|
|
|
#
|
404
|
|
|
|
|
|
|
# _RETURN
|
405
|
|
|
|
|
|
|
#
|
406
|
|
|
|
|
|
|
# the value of the input is >= 0.0 and < 360.0
|
407
|
|
|
|
|
|
|
#
|
408
|
|
|
|
|
|
|
|
409
|
3
|
|
|
3
|
0
|
5
|
my $x = $_[0];
|
410
|
3
|
|
|
|
|
23
|
return ( $x - 360.0 * floor( $x * $INV360 ) );
|
411
|
|
|
|
|
|
|
}
|
412
|
|
|
|
|
|
|
|
413
|
|
|
|
|
|
|
sub rev180 {
|
414
|
|
|
|
|
|
|
|
415
|
|
|
|
|
|
|
#
|
416
|
|
|
|
|
|
|
#
|
417
|
|
|
|
|
|
|
# FUNCTIONAL SEQUENCE for rev180
|
418
|
|
|
|
|
|
|
#
|
419
|
|
|
|
|
|
|
# _GIVEN
|
420
|
|
|
|
|
|
|
#
|
421
|
|
|
|
|
|
|
# any angle
|
422
|
|
|
|
|
|
|
#
|
423
|
|
|
|
|
|
|
# _THEN
|
424
|
|
|
|
|
|
|
#
|
425
|
|
|
|
|
|
|
# Reduce input to within +180..+180 degrees
|
426
|
|
|
|
|
|
|
#
|
427
|
|
|
|
|
|
|
#
|
428
|
|
|
|
|
|
|
# _RETURN
|
429
|
|
|
|
|
|
|
#
|
430
|
|
|
|
|
|
|
# angle that was reduced
|
431
|
|
|
|
|
|
|
#
|
432
|
1
|
|
|
1
|
0
|
2
|
my ($x) = @_;
|
433
|
|
|
|
|
|
|
|
434
|
1
|
|
|
|
|
6
|
return ( $x - 360.0 * floor( $x * $INV360 + 0.5 ) );
|
435
|
|
|
|
|
|
|
}
|
436
|
|
|
|
|
|
|
|
437
|
|
|
|
|
|
|
sub equal {
|
438
|
0
|
|
|
0
|
0
|
0
|
my ( $A, $B, $dp ) = @_;
|
439
|
|
|
|
|
|
|
|
440
|
0
|
|
|
|
|
0
|
return sprintf( "%.${dp}g", $A ) eq sprintf( "%.${dp}g", $B );
|
441
|
|
|
|
|
|
|
}
|
442
|
|
|
|
|
|
|
|
443
|
|
|
|
|
|
|
sub convert_hour {
|
444
|
|
|
|
|
|
|
|
445
|
|
|
|
|
|
|
#
|
446
|
|
|
|
|
|
|
#
|
447
|
|
|
|
|
|
|
# FUNCTIONAL SEQUENCE for convert_hour
|
448
|
|
|
|
|
|
|
#
|
449
|
|
|
|
|
|
|
# _GIVEN
|
450
|
|
|
|
|
|
|
# Hour_rise, Hour_set
|
451
|
|
|
|
|
|
|
# hours are in UT
|
452
|
|
|
|
|
|
|
#
|
453
|
|
|
|
|
|
|
# _THEN
|
454
|
|
|
|
|
|
|
#
|
455
|
|
|
|
|
|
|
# split out the hours and minites
|
456
|
|
|
|
|
|
|
#
|
457
|
|
|
|
|
|
|
#
|
458
|
|
|
|
|
|
|
# _RETURN
|
459
|
|
|
|
|
|
|
#
|
460
|
|
|
|
|
|
|
# hour:min rise and set
|
461
|
|
|
|
|
|
|
#
|
462
|
|
|
|
|
|
|
|
463
|
1
|
|
|
1
|
0
|
1
|
my ( $hour_rise_ut, $hour_set_ut ) = @_;
|
464
|
|
|
|
|
|
|
|
465
|
1
|
|
|
|
|
4
|
my $min_rise = int( ( $hour_rise_ut - int($hour_rise_ut) ) * 60 );
|
466
|
1
|
|
|
|
|
2
|
my $min_set = int( ( $hour_set_ut - int($hour_set_ut) ) * 60 );
|
467
|
|
|
|
|
|
|
|
468
|
1
|
|
|
|
|
2
|
my $hour_rise = int($hour_rise_ut);
|
469
|
1
|
|
|
|
|
7
|
my $hour_set = int($hour_set_ut);
|
470
|
1
|
50
|
|
|
|
4
|
if ( $min_rise < 10 ) {
|
471
|
0
|
|
|
|
|
0
|
$min_rise = sprintf( "%02d", $min_rise );
|
472
|
|
|
|
|
|
|
}
|
473
|
|
|
|
|
|
|
|
474
|
1
|
50
|
|
|
|
3
|
if ( $min_set < 10 ) {
|
475
|
1
|
|
|
|
|
8
|
$min_set = sprintf( "%02d", $min_set );
|
476
|
|
|
|
|
|
|
}
|
477
|
|
|
|
|
|
|
|
478
|
1
|
|
|
|
|
2
|
return ( $hour_rise, $min_rise, $hour_set, $min_set );
|
479
|
|
|
|
|
|
|
}
|
480
|
|
|
|
|
|
|
=head1 NAME
|
481
|
|
|
|
|
|
|
|
482
|
|
|
|
|
|
|
DateTime::Astro::Sunrise - Perl DateTime extension for computing the sunrise/sunset on a given day
|
483
|
|
|
|
|
|
|
|
484
|
|
|
|
|
|
|
=head1 SYNOPSIS
|
485
|
|
|
|
|
|
|
use DateTime;
|
486
|
|
|
|
|
|
|
use DateTime::Astro::Sunrise;
|
487
|
|
|
|
|
|
|
|
488
|
|
|
|
|
|
|
my $dt = DateTime->new( year => 2000,
|
489
|
|
|
|
|
|
|
month => 6,
|
490
|
|
|
|
|
|
|
day => 20,
|
491
|
|
|
|
|
|
|
);
|
492
|
|
|
|
|
|
|
|
493
|
|
|
|
|
|
|
my $sunrise = DateTime::Astro::Sunrise ->new('-118','33',undef,1);
|
494
|
|
|
|
|
|
|
|
495
|
|
|
|
|
|
|
my ($tmp_rise, $tmp_set) = $sunrise->sunrise($dt);
|
496
|
|
|
|
|
|
|
|
497
|
|
|
|
|
|
|
|
498
|
|
|
|
|
|
|
=head1 DESCRIPTION
|
499
|
|
|
|
|
|
|
|
500
|
|
|
|
|
|
|
This module will return a DateTime Object for sunrise and sunset for a given day.
|
501
|
|
|
|
|
|
|
|
502
|
|
|
|
|
|
|
=head1 USAGE
|
503
|
|
|
|
|
|
|
|
504
|
|
|
|
|
|
|
=over
|
505
|
|
|
|
|
|
|
|
506
|
|
|
|
|
|
|
=item Bnew(longitutide,latatude,ALT,Iteration);>
|
507
|
|
|
|
|
|
|
|
508
|
|
|
|
|
|
|
=over
|
509
|
|
|
|
|
|
|
Eastern longitude is entered as a positive number
|
510
|
|
|
|
|
|
|
Western longitude is entered as a negative number
|
511
|
|
|
|
|
|
|
Northern latitude is entered as a positive number
|
512
|
|
|
|
|
|
|
Southern latitude is entered as a negative number
|
513
|
|
|
|
|
|
|
|
514
|
|
|
|
|
|
|
inter is set to either 0 or 1.
|
515
|
|
|
|
|
|
|
If set to 0 no Iteration will occur.
|
516
|
|
|
|
|
|
|
If set to 1 Iteration will occur.
|
517
|
|
|
|
|
|
|
Default is 0.
|
518
|
|
|
|
|
|
|
|
519
|
|
|
|
|
|
|
There are a number of sun altitides to chose from. The default is
|
520
|
|
|
|
|
|
|
-0.833 because this is what most countries use. Feel free to
|
521
|
|
|
|
|
|
|
specify it if you need to. Here is the list of values to specify
|
522
|
|
|
|
|
|
|
altitude (ALT) with:
|
523
|
|
|
|
|
|
|
|
524
|
|
|
|
|
|
|
=over
|
525
|
|
|
|
|
|
|
|
526
|
|
|
|
|
|
|
=item B<0> degrees
|
527
|
|
|
|
|
|
|
|
528
|
|
|
|
|
|
|
Center of Sun's disk touches a mathematical horizon
|
529
|
|
|
|
|
|
|
|
530
|
|
|
|
|
|
|
=item B<-0.25> degrees
|
531
|
|
|
|
|
|
|
|
532
|
|
|
|
|
|
|
Sun's upper limb touches a mathematical horizon
|
533
|
|
|
|
|
|
|
|
534
|
|
|
|
|
|
|
=item B<-0.583> degrees
|
535
|
|
|
|
|
|
|
|
536
|
|
|
|
|
|
|
Center of Sun's disk touches the horizon; atmospheric refraction accounted for
|
537
|
|
|
|
|
|
|
|
538
|
|
|
|
|
|
|
=item B<-0.833> degrees
|
539
|
|
|
|
|
|
|
|
540
|
|
|
|
|
|
|
Sun's supper limb touches the horizon; atmospheric refraction accounted for
|
541
|
|
|
|
|
|
|
|
542
|
|
|
|
|
|
|
=item B<-6> degrees
|
543
|
|
|
|
|
|
|
|
544
|
|
|
|
|
|
|
Civil twilight (one can no longer read outside without artificial illumination)
|
545
|
|
|
|
|
|
|
|
546
|
|
|
|
|
|
|
=item B<-12> degrees
|
547
|
|
|
|
|
|
|
|
548
|
|
|
|
|
|
|
Nautical twilight (navigation using a sea horizon no longer possible)
|
549
|
|
|
|
|
|
|
|
550
|
|
|
|
|
|
|
=item B<-15> degrees
|
551
|
|
|
|
|
|
|
|
552
|
|
|
|
|
|
|
Amateur astronomical twilight (the sky is dark enough for most astronomical observations)
|
553
|
|
|
|
|
|
|
|
554
|
|
|
|
|
|
|
=item B<-18> degrees
|
555
|
|
|
|
|
|
|
|
556
|
|
|
|
|
|
|
Astronomical twilight (the sky is completely dark)
|
557
|
|
|
|
|
|
|
|
558
|
|
|
|
|
|
|
=item F
|
559
|
|
|
|
|
|
|
|
560
|
|
|
|
|
|
|
=over
|
561
|
|
|
|
|
|
|
|
562
|
|
|
|
|
|
|
The orginal method only gives an approximate value of the Sun's rise/set times.
|
563
|
|
|
|
|
|
|
The error rarely exceeds one or two minutes, but at high latitudes, when the Midnight Sun
|
564
|
|
|
|
|
|
|
soon will start or just has ended, the errors may be much larger. If you want higher accuracy,
|
565
|
|
|
|
|
|
|
you must then use the iteration feature. This feature is new as of version 0.7. Here is
|
566
|
|
|
|
|
|
|
what I have tried to accomplish with this.
|
567
|
|
|
|
|
|
|
|
568
|
|
|
|
|
|
|
a) Compute sunrise or sunset as always, with one exception: to convert LHA from degrees to hours,
|
569
|
|
|
|
|
|
|
divide by 15.04107 instead of 15.0 (this accounts for the difference between the solar day
|
570
|
|
|
|
|
|
|
and the sidereal day.
|
571
|
|
|
|
|
|
|
|
572
|
|
|
|
|
|
|
b) Re-do the computation but compute the Sun's RA and Decl, and also GMST0, for the moment
|
573
|
|
|
|
|
|
|
of sunrise or sunset last computed.
|
574
|
|
|
|
|
|
|
|
575
|
|
|
|
|
|
|
c) Iterate b) until the computed sunrise or sunset no longer changes significantly.
|
576
|
|
|
|
|
|
|
Usually 2 iterations are enough, in rare cases 3 or 4 iterations may be needed.
|
577
|
|
|
|
|
|
|
|
578
|
|
|
|
|
|
|
|
579
|
|
|
|
|
|
|
|
580
|
|
|
|
|
|
|
=back
|
581
|
|
|
|
|
|
|
|
582
|
|
|
|
|
|
|
=back
|
583
|
|
|
|
|
|
|
|
584
|
|
|
|
|
|
|
=head1 ($sunrise, $sunset) = $sunrise->($dt);
|
585
|
|
|
|
|
|
|
|
586
|
|
|
|
|
|
|
|
587
|
|
|
|
|
|
|
Returns two DateTime objects sunrise and sunset.
|
588
|
|
|
|
|
|
|
Please note that the time zone for these objects
|
589
|
|
|
|
|
|
|
is set to UTC. So don't forget to set your timezone!!
|
590
|
|
|
|
|
|
|
|
591
|
|
|
|
|
|
|
|
592
|
|
|
|
|
|
|
|
593
|
|
|
|
|
|
|
=head1 AUTHOR
|
594
|
|
|
|
|
|
|
|
595
|
|
|
|
|
|
|
Ron Hill
|
596
|
|
|
|
|
|
|
rkhill@firstlight.net
|
597
|
|
|
|
|
|
|
|
598
|
|
|
|
|
|
|
=head1 CREDITS
|
599
|
|
|
|
|
|
|
|
600
|
|
|
|
|
|
|
|
601
|
|
|
|
|
|
|
=item Paul Schlyer, Stockholm, Sweden
|
602
|
|
|
|
|
|
|
|
603
|
|
|
|
|
|
|
for his excellent web page on the subject.
|
604
|
|
|
|
|
|
|
|
605
|
|
|
|
|
|
|
=item Rich Bowen (rbowen@rbowen.com)
|
606
|
|
|
|
|
|
|
|
607
|
|
|
|
|
|
|
for suggestions
|
608
|
|
|
|
|
|
|
|
609
|
|
|
|
|
|
|
=head1 COPYRIGHT and LICENSE
|
610
|
|
|
|
|
|
|
|
611
|
|
|
|
|
|
|
Here is the copyright information provided by Paul Schlyer:
|
612
|
|
|
|
|
|
|
|
613
|
|
|
|
|
|
|
Written as DAYLEN.C, 1989-08-16
|
614
|
|
|
|
|
|
|
|
615
|
|
|
|
|
|
|
Modified to SUNRISET.C, 1992-12-01
|
616
|
|
|
|
|
|
|
|
617
|
|
|
|
|
|
|
(c) Paul Schlyter, 1989, 1992
|
618
|
|
|
|
|
|
|
|
619
|
|
|
|
|
|
|
Released to the public domain by Paul Schlyter, December 1992
|
620
|
|
|
|
|
|
|
|
621
|
|
|
|
|
|
|
Permission is hereby granted, free of charge, to any person obtaining a
|
622
|
|
|
|
|
|
|
copy of this software and associated documentation files (the "Software"),
|
623
|
|
|
|
|
|
|
to deal in the Software without restriction, including without limitation
|
624
|
|
|
|
|
|
|
the rights to use, copy, modify, merge, publish, distribute, sublicense,
|
625
|
|
|
|
|
|
|
and/or sell copies of the Software, and to permit persons to whom the
|
626
|
|
|
|
|
|
|
Software is furnished to do so, subject to the following conditions:
|
627
|
|
|
|
|
|
|
|
628
|
|
|
|
|
|
|
The above copyright notice and this permission notice shall be included
|
629
|
|
|
|
|
|
|
in all copies or substantial portions of the Software.
|
630
|
|
|
|
|
|
|
|
631
|
|
|
|
|
|
|
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
|
632
|
|
|
|
|
|
|
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
|
633
|
|
|
|
|
|
|
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
|
634
|
|
|
|
|
|
|
THE AUTHOR BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
|
635
|
|
|
|
|
|
|
WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT
|
636
|
|
|
|
|
|
|
OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
|
637
|
|
|
|
|
|
|
THE SOFTWARE.
|
638
|
|
|
|
|
|
|
|
639
|
|
|
|
|
|
|
=head1 BUGS
|
640
|
|
|
|
|
|
|
|
641
|
|
|
|
|
|
|
=head1 SEE ALSO
|
642
|
|
|
|
|
|
|
|
643
|
|
|
|
|
|
|
perl(1).
|
644
|
|
|
|
|
|
|
|
645
|
|
|
|
|
|
|
=cut
|
646
|
|
|
|
|
|
|
1;
|
647
|
|
|
|
|
|
|
|