2 * This file is part of Live Wallpaper (livewp)
4 * Copyright (C) 2010 Vlad Vasiliev
5 * Copyright (C) 2010 Tanya Makova
8 * This software is free software; you can redistribute it and/or
9 * modify it under the terms of the GNU Lesser General Public License
10 * as published by the Free Software Foundation; either version 2.1 of
11 * the License, or (at your option) any later version.
13 * This software is distributed in the hope that it will be useful, but
14 * WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 * Lesser General Public License for more details.
18 * You should have received a copy of the GNU Lesser General Public
19 * License along with this software; if not, write to the Free Software
20 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA
23 /*******************************************************************************/
25 #include "livewp-astro.h"
27 double ut(int hour, int min, int zone)
29 return (hour - zone + (double)min/60);
32 double jd(int year, int month, int day, double uu)
34 return (367*year)-floor((floor((month+9)/12)+year)*7/4)+floor(275*month/9)+day-730531.5+(uu/24);
37 double altitude(double lon, double lat, int year, int month, int day, int hour, int min, int zone)
39 double uu, jj, T, k, M, Lo, DL, L, eps, delta, RA, GMST, LMST, H, eqt, alt;
40 uu = ut(hour, min, zone);
41 jj = jd(year, month, day, uu);
44 M = 357.5291 + 35999.0503 * T - 0.0001559 * T * T - 0.00000045 * T * T * T;
46 Lo = 280.46645 + 36000.76983 * T + 0.0003032 * T * T;
48 DL = (1.9146-0.004817*T-0.000014*T*T)*sin(k*M)+(0.019993-0.000101*T)*sin(k*2*M)+0.00029*sin(k*3*M);
51 delta=(1/k)*asin(sin(L*k)*sin(eps*k));
52 RA=(1/k)*atan2(cos(eps*k)*sin(L*k),cos(L*k));
53 RA=fmod((RA+360), 360);
54 GMST=280.46061837+360.98564736629*jj+0.000387933*T*T-T*T*T/38710000;
55 GMST=fmod((GMST+360), 360);
59 alt=(1/k)*asin(sin(lat*k)*sin(delta*k)+cos(lat*k)*cos(delta*k)*cos(H*k));
63 double azimuth(double lon, double lat, int year, int month, int day, int hour, int min, int zone)
65 double uu, jj, T, k, M, Lo, DL, L, eps, delta, RA, GMST, LMST, H, eqt, azm;
66 uu = ut(hour, min, zone);
67 jj = jd(year, month, day, uu);
70 M = 357.5291 + 35999.0503 * T - 0.0001559 * T * T - 0.00000045 * T * T * T;
72 Lo = 280.46645 + 36000.76983 * T + 0.0003032 * T * T;
74 DL = (1.9146-0.004817*T-0.000014*T*T)*sin(k*M)+(0.019993-0.000101*T)*sin(k*2*M)+0.00029*sin(k*3*M);
77 delta=(1/k)*asin(sin(L*k)*sin(eps*k));
78 RA=(1/k)*atan2(cos(eps*k)*sin(L*k),cos(L*k));
79 RA=fmod((RA+360), 360);
80 GMST=280.46061837+360.98564736629*jj+0.000387933*T*T-T*T*T/38710000;
81 GMST=fmod((GMST+360), 360);
85 azm=(1/k)*atan2(-sin(H*k),cos(lat*k)*tan(delta*k)-sin(lat*k)*cos(H*k));
86 azm=fmod((azm+360), 360);
89 int day_of_year(int year, int month, int day)
92 n1 = floor(275 * month / 9);
93 n2 = floor((month + 9) / 12);
94 n3 = (1 + floor((year - 4 * floor(year / 4) + 2) / 3));
95 n = n1 - (n2 * n3) + day - 30;
98 double deg2rad(double deg)
100 return M_PI * deg / 180;
102 double truel(double m)
105 l = m + (1.916 * sin(deg2rad(m))) + (0.020 * sin(deg2rad(2*m))) + 282.634;
112 double rightasc(double l)
115 ra = 180 * atan(0.91764 * tan(l * M_PI/180)) / M_PI;
116 lq = 90 * floor(l/90);
117 rq = 90 * floor(ra / 90);
121 double locangle(double lat, double sindec, double cosdec)
123 double cosh = (-0.01454 - sindec * sin(deg2rad(lat))) / (cosdec * cos(deg2rad(lat)));
126 double tout(double t, double lng_hour)
128 double ut = t - lng_hour;
136 sun_rise_set(double lon, double lat,
137 int year, int month, int day,
139 int * hourrise, int *minrise, int *hourset, int *minset)
142 double lng_hour, trise, tset, mrise, mset, lrise, lset,
143 rarise, raset, sindecrise, cosdecrise, sindecset, cosdecset,
145 hrise, hset, ttrise, ttset, utrise, utset;
146 int utrisehour, utrisemin, utsethour, utsetmin;
147 n = day_of_year(year, month, day);
149 trise = (double)n + (double)((6 - lng_hour) / 24);
150 tset = (double)n + (double)((18 - lng_hour) / 24);
151 mrise = (0.9856 * trise) - 3.289;
152 mset = (0.9856 * tset) - 3.289;
153 lrise = truel(mrise);
155 rarise = rightasc(lrise);
156 raset = rightasc(lset);
157 sindecrise = 0.39782 * sin(deg2rad(lrise));
158 cosdecrise = cos(asin(sindecrise));
159 sindecset = 0.39782 * sin(deg2rad(lset));
160 cosdecset = cos(asin(sindecset));
161 coshrise = locangle(lat, sindecrise, cosdecrise);
162 /* if (coshrise > 1) sun never rises
163 * if (coshrise < -1) sun never sets
165 coshset = locangle(lat, sindecset, cosdecset);
166 hrise = (360 - 180 * acos(coshrise) / M_PI) / 15;
167 hset = (180 * acos(coshrise) / M_PI) / 15;
168 ttrise = hrise + rarise - (0.06571 * trise) - 6.622;
169 ttset = hset + raset - (0.06571 * tset) - 6.622;
170 utrise = tout(ttrise, lng_hour);
171 utset = tout(ttset, lng_hour);
173 utrisehour = floor(utrise);
174 utrisemin = floor((utrise - utrisehour) * 60);
175 utsethour = floor(utset);
176 utsetmin = floor((utset - utsethour) * 60);
178 *hourrise = utrisehour + zone;
181 else if (*hourrise < 0)
183 *minrise = utrisemin;
184 *hourset = utsethour + zone;
187 else if (*hourset < 0)
194 double moon_phase(int year, int month, int day)
197 double agepart, ge, jd, moonday;
199 moonday = 29.530588853;
200 if (month <= 2) lg = 0;
202 if ((fmod(year, 4) == 0) && (!(fmod(year, 100) == 0) && (fmod(year, 400) != 0)))
208 floor((year - 1) / 4) +
209 (-floor((year - 1) / 100)) +
210 floor((year - 1) / 400) +
211 floor((((367 * month) - 362) / 12) + lg + day);
213 agepart = (jd - 2451550.1) / moonday;
214 agepart = agepart - floor(agepart);
215 if (agepart < 0) agepart ++;
217 ageday = agepart * moonday;
219 moonfill = 1.0 - ((moonday/2 - ageday) / (moonday/2));
220 fprintf(stderr, "agepart = %f, ageday = %f, moonfill = %f\n", agepart, ageday, moonfill);
226 int main(int argc, char * argv)
228 int year, month, day, hour, min, zone;
229 double alt, azm, lon, lat, moonpart;
239 //moonpart = moon_phase1(year, month, day);
240 //fprintf(stderr, "moon path = %f\n", moonpart);
241 alt = altitude(lon, lat, year, month, day, hour, min, zone);
242 azm = azimuth(lon, lat, year, month, day, hour, min, zone);
243 fprintf(stderr, "alt=%5.6f azm=%5.6f\n", alt, azm);