1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
|
/* libhdate - Hebrew calendar library
*
* Copyright (C) 1984-2003 Amos Shapir, 2004-2007 Yaacov Zamir <kzamir@walla.co.il>
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
/*
* Algorithm from http://www.srrb.noaa.gov/highlights/sunrise/calcdetails.html
* The low accuracy solar position equations are used.
* These routines are based on Jean Meeus's book Astronomical Algorithms.
*/
#include <time.h>
#include <stdio.h>
#include <math.h>
#include <errno.h>
#include "hdate.h"
#include "support.h"
/**
@brief days from 1 january
@parm day this day of month
@parm month this month
@parm year this year
@return the days from 1 jan
*/
int
hdate_get_day_of_year (int day, int month, int year)
{
int jd;
/* get todays julian day number */
jd = (1461 * (year + 4800 + (month - 14) / 12)) / 4 +
(367 * (month - 2 - 12 * ((month - 14) / 12))) / 12 -
(3 * ((year + 4900 + (month - 14) / 12) / 100)) / 4 + day;
/* substruct the julian day of 1/1/year and add one */
jd = jd - ((1461 * (year + 4799)) / 4 +
367 * 11 / 12 - (3 * ((year + 4899) / 100)) / 4);
return jd;
}
/**
@brief utc sun times for altitude at a gregorian date
Returns the sunset and sunrise times in minutes from 00:00 (utc time)
if sun altitude in sunrise is deg degries.
This function only works for altitudes sun realy is.
If the sun never get to this altitude, the returned sunset and sunrise values
will be negative. This can happen in low altitude when latitude is
nearing the pols in winter times, the sun never goes very high in
the sky there.
@param day this day of month
@param month this month
@param year this year
@param longitude longitude to use in calculations
@param latitude latitude to use in calculations
@param deg degrees of sun's altitude (0 - Zenith .. 90 - Horizon)
@param sunrise return the utc sunrise in minutes
@param sunset return the utc sunset in minutes
*/
void
hdate_get_utc_sun_time_deg (int day, int month, int year, double latitude, double longitude, double deg, int *sunrise, int *sunset)
{
double gama; /* location of sun in yearly cycle in radians */
double eqtime; /* diffference betwen sun noon and clock noon */
double decl; /* sun declanation */
double ha; /* solar hour engle */
double sunrise_angle = M_PI * deg / 180.0; /* sun angle at sunrise/set */
int day_of_year;
/* get the day of year */
day_of_year = hdate_get_day_of_year (day, month, year);
/* get radians of sun orbit around erth =) */
gama = 2.0 * M_PI * ((double)(day_of_year - 1) / 365.0);
/* get the diff betwen suns clock and wall clock in minutes */
eqtime = 229.18 * (0.000075 + 0.001868 * cos (gama)
- 0.032077 * sin (gama) - 0.014615 * cos (2.0 * gama)
- 0.040849 * sin (2.0 * gama));
/* calculate suns declanation at the equater in radians */
decl = 0.006918 - 0.399912 * cos (gama) + 0.070257 * sin (gama)
- 0.006758 * cos (2.0 * gama) + 0.000907 * sin (2.0 * gama)
- 0.002697 * cos (3.0 * gama) + 0.00148 * sin (3.0 * gama);
/* we use radians, ratio is 2pi/360 */
latitude = M_PI * latitude / 180.0;
/* the sun real time diff from noon at sunset/rise in radians */
errno = 0;
ha = acos (cos (sunrise_angle) / (cos (latitude) * cos (decl)) - tan (latitude) * tan (decl));
/* check for too high altitudes and return negative values */
if (errno == EDOM)
{
*sunrise = -720;
*sunset = -720;
return;
}
/* we use minutes, ratio is 1440min/2pi */
ha = 720.0 * ha / M_PI;
/* get sunset/rise times in utc wall clock in minutes from 00:00 time */
*sunrise = (int)(720.0 - 4.0 * longitude - ha - eqtime);
*sunset = (int)(720.0 - 4.0 * longitude + ha - eqtime);
return;
}
/**
@brief utc sunrise/set time for a gregorian date
@parm day this day of month
@parm month this month
@parm year this year
@parm longitude longitude to use in calculations
@parm latitude latitude to use in calculations
@parm sunrise return the utc sunrise in minutes
@parm sunset return the utc sunset in minutes
*/
void
hdate_get_utc_sun_time (int day, int month, int year, double latitude, double longitude, int *sunrise, int *sunset)
{
hdate_get_utc_sun_time_deg (day, month, year, latitude, longitude, 90.833, sunrise, sunset);
return;
}
/**
@brief utc sunrise/set time for a gregorian date
@parm day this day of month
@parm month this month
@parm year this year
@parm longitude longitude to use in calculations
@parm latitude latitude to use in calculations
@parm sun_hour return the length of shaa zaminit in minutes
@parm first_light return the utc alut ha-shachar in minutes
@parm talit return the utc tphilin and talit in minutes
@parm sunrise return the utc sunrise in minutes
@parm midday return the utc midday in minutes
@parm sunset return the utc sunset in minutes
@parm first_stars return the utc tzeit hacochavim in minutes
@parm three_stars return the utc shlosha cochavim in minutes
*/
void
hdate_get_utc_sun_time_full (int day, int month, int year, double latitude, double longitude,
int *sun_hour, int *first_light, int *talit, int *sunrise,
int *midday, int *sunset, int *first_stars, int *three_stars)
{
int place_holder;
/* sunset and rise time */
hdate_get_utc_sun_time_deg (day, month, year, latitude, longitude, 90.833, sunrise, sunset);
/* shaa zmanit by gara, 1/12 of light time */
*sun_hour = (*sunset - *sunrise) / 12;
*midday = (*sunset + *sunrise) / 2;
/* get times of the different sun angles */
hdate_get_utc_sun_time_deg (day, month, year, latitude, longitude, 106.01, first_light, &place_holder);
hdate_get_utc_sun_time_deg (day, month, year, latitude, longitude, 101.0, talit, &place_holder);
hdate_get_utc_sun_time_deg (day, month, year, latitude, longitude, 96.0, &place_holder, first_stars);
hdate_get_utc_sun_time_deg (day, month, year, latitude, longitude, 98.5, &place_holder, three_stars);
return;
}
|
