1
2# Copyright (c) 1999-2017 Rob Fugina <robf@fugina.com>
3# Distributed under the terms of the GNU Public License, Version 3.0
4
5package Astro::SunTime;
6use vars qw(@ISA @EXPORT $VERSION);
7$VERSION = 0.06;
8@ISA = qw(Exporter);
9@EXPORT = qw(sun_time);
10
11#  Results can be checked with: http://aa.usno.navy.mil/data/docs/RS_OneYear.php
12
13# 09/03/00 :: winter Make ParseDate optional.  It is overkill and I could not get it to
14#                    compile in perl2exe.  It gave runaway comment errors :(
15# 10/12/00 :: winter Change time_zone check to defined, to allow for time_zone 0
16
17use POSIX;
18
19use strict;
20
21# sun_time takes:
22# type => 'rise' | 'set'
23# latitude
24# longitude
25# time_zone => hours from GMT
26# date => date parsable by Time::ParseDate::parsedate()
27# time => to feed to localtime
28
29sub sun_time
30{
31   my %params = @_;
32
33   my $type = $params{type} || 'rise';
34   my $latitude  = (defined $params{latitude}) ? $params{latitude} : 38.74274;
35   my $longitude = (defined $params{longitude}) ? $params{longitude} : -90.560143;
36   my $time_zone = (defined $params{time_zone}) ? $params{time_zone} : -6;
37
38   my $time;
39   if ($params{date}) {
40       require Time::ParseDate;
41       $time = Time::ParseDate::parsedate($params{date});
42   }
43   elsif ($params{time}) {
44       $time = $params{time};
45   }
46   else {
47       $time = time;
48   }
49   my @suntime = localtime($time);
50
51   my $yday = $suntime[7] + 1;
52
53   my $A = 1.5708;
54   my $B = 3.14159;
55   my $C = 4.71239;
56   my $D = 6.28319;
57   my $E = 0.0174533 * $latitude;
58   my $F = 0.0174533 * $longitude;
59   my $G = 0.261799  * $time_zone;
60
61   # For astronomical twilight, use R = -.309017
62   # For     nautical twilight, use R = -.207912
63   # For        civil twilight, use R = -.104528
64   # For     sunrise or sunset, use R = -.0145439
65
66   my $R = -.0145439;
67   if ($params{twilight}) {
68       if($params{twilight} eq 'astronomical') {
69           $R = -.309017;
70       }
71       elsif($params{twilight} eq 'nautical') {
72           $R = -.207912;
73       }
74       elsif($params{twilight} eq 'civil') {
75           $R = -.104528;
76       }
77   }
78
79
80   my $J = ($type eq 'rise') ? $A : $C;
81   my $K = $yday + (($J - $F) / $D);
82   my $L = ($K * .017202) - .0574039;   # Solar Mean Anomoly
83   my $M = $L + .0334405 * sin($L);     # Solar True Longitude
84   $M += 4.93289 + (3.49066E-04) * sin(2 * $L);
85   $M = &normalize($M, $D);              # Quadrant Determination
86   $M += 4.84814E-06 if ($M / $A) - int($M / $A) == 0;
87   my $P = sin($M) / cos($M);           # Solar Right Ascension
88   $P = atan2(.91746 * $P, 1);
89
90   # Quadrant Adjustment
91   if ($M > $C)
92   {
93      $P += $D;
94   }
95   elsif ($M > $A)
96   {
97      $P += $B;
98   }
99
100   my $Q = .39782 * sin($M);            # Solar Declination
101   $Q = $Q / sqrt(-$Q * $Q + 1);     # This is how the original author wrote it!
102   $Q = atan2($Q, 1);
103
104   my $S = $R - (sin($Q) * sin($E));
105   $S = $S / (cos($Q) * cos($E));
106
107   return 'none' if abs($S) > 1;      # Null phenomenon
108
109   $S = $S / sqrt(-$S * $S + 1);
110   $S = $A - atan2($S, 1);
111   $S = $D - $S if $type eq 'rise';
112
113   my $T = $S + $P - 0.0172028 * $K - 1.73364; # Local apparent time
114   my $U = $T - $F;                     # Universal timer
115   my $V = $U + $G;                     # Wall clock time
116   $V = &normalize($V, $D);
117   $V = $V * 3.81972;
118
119   my $hour = int($V);
120   my $min  = int(($V - $hour) * 60 + 0.5);
121
122   @suntime[2,1,0] = ($hour, $min, 0);
123
124   @suntime = localtime(mktime(@suntime)); # normalize date structure
125
126   return sprintf("%d:%02d", @suntime[2,1]);
127}
128
129sub normalize
130{
131   my $Z = shift;
132   my $D = shift;
133
134   die "Trying to normalize with zero offset..." if ($D == 0);
135
136   while ($Z < 0)   {$Z = $Z + $D}
137   while ($Z >= $D) {$Z = $Z - $D}
138
139   return $Z;
140}
141
142
143
1441;
145
146__END__
147
148
149