+/*
+ * This file is part of Live Wallpaper (livewp)
+ *
+ * Copyright (C) 2010 Vlad Vasiliev
+ * Copyright (C) 2010 Tanya Makova
+ * for the code
+ *
+ * This software is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1 of
+ * the License, or (at your option) any later version.
+ *
+ * This software 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
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this software; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA
+ * 02110-1301 USA
+*/
+/*******************************************************************************/
+#include <stdio.h>
+#include "livewp-astro.h"
+
+double ut(int hour, int min, int zone)
+{
+ return (hour - zone + (double)min/60);
+}
+
+double jd(int year, int month, int day, double uu)
+{
+ return (367*year)-floor((floor((month+9)/12)+year)*7/4)+floor(275*month/9)+day-730531.5+(uu/24);
+}
+
+double altitude(double lon, double lat, int year, int month, int day, int hour, int min, int zone)
+{
+ double uu, jj, T, k, M, Lo, DL, L, eps, delta, RA, GMST, LMST, H, eqt, alt;
+ uu = ut(hour, min, zone);
+ jj = jd(year, month, day, uu);
+ T = jj / 36525;
+ k = M_PI / 180;
+ M = 357.5291 + 35999.0503 * T - 0.0001559 * T * T - 0.00000045 * T * T * T;
+ M = fmod(M, 360);
+ Lo = 280.46645 + 36000.76983 * T + 0.0003032 * T * T;
+ Lo = fmod(Lo, 360);
+ 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);
+ L=Lo+DL;
+ eps=23.43999-0.013*T;
+ delta=(1/k)*asin(sin(L*k)*sin(eps*k));
+ RA=(1/k)*atan2(cos(eps*k)*sin(L*k),cos(L*k));
+ RA=fmod((RA+360), 360);
+ GMST=280.46061837+360.98564736629*jj+0.000387933*T*T-T*T*T/38710000;
+ GMST=fmod((GMST+360), 360);
+ LMST=GMST+lon;
+ H=LMST-RA;
+ eqt=(Lo-RA)*4;
+ alt=(1/k)*asin(sin(lat*k)*sin(delta*k)+cos(lat*k)*cos(delta*k)*cos(H*k));
+ return alt;
+}
+
+double azimuth(double lon, double lat, int year, int month, int day, int hour, int min, int zone)
+{
+ double uu, jj, T, k, M, Lo, DL, L, eps, delta, RA, GMST, LMST, H, eqt, azm;
+ uu = ut(hour, min, zone);
+ fprintf(stderr, "%f, %f, %i.%i.%i %i:%i %i uu=%f", lon, lat, year, month, day, hour, min, zone, uu);
+ jj = jd(year, month, day, uu);
+ T = jj / 36525;
+ k = M_PI / 180;
+ M = 357.5291 + 35999.0503 * T - 0.0001559 * T * T - 0.00000045 * T * T * T;
+ M = fmod(M, 360);
+ Lo = 280.46645 + 36000.76983 * T + 0.0003032 * T * T;
+ Lo = fmod(Lo, 360);
+ 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);
+ L=Lo+DL;
+ eps=23.43999-0.013*T;
+ delta=(1/k)*asin(sin(L*k)*sin(eps*k));
+ RA=(1/k)*atan2(cos(eps*k)*sin(L*k),cos(L*k));
+ RA=fmod((RA+360), 360);
+ GMST=280.46061837+360.98564736629*jj+0.000387933*T*T-T*T*T/38710000;
+ GMST=fmod((GMST+360), 360);
+ LMST=GMST+lon;
+ H=LMST-RA;
+ eqt=(Lo-RA)*4;
+ azm=(1/k)*atan2(-sin(H*k),cos(lat*k)*tan(delta*k)-sin(lat*k)*cos(H*k));
+ azm=fmod((azm+360), 360);
+ return azm;
+}
+double moon_phase(int year, int month, int day)
+{
+ int yy, mm, dd, k1, k2, k3, lg;
+ double agepart, ageday, ge, moonday, jd, moonfill;
+ ge = 1721425.5;
+ moonday = 29.530588853;
+ if (month <= 2) lg = 0;
+ else {
+ if ((fmod(year, 4) == 0) && (!(fmod(year, 100) == 0) && (fmod(year, 400) != 0)))
+ lg = -1;
+ else lg = -2;
+ }
+ jd = (ge - 1) +
+ (365 * (year - 1)) +
+ floor((year - 1) / 4) +
+ (-floor((year - 1) / 100)) +
+ floor((year - 1) / 400) +
+ floor((((367 * month) - 362) / 12) + lg + day);
+
+ agepart = (jd - 2451550.1) / moonday;
+ agepart = agepart - floor(agepart);
+ if (agepart < 0) agepart + 1;
+/*
+ ageday = agepart * moonday;
+
+ moonfill = 1.0 - ((moonday/2 - ageday) / (moonday/2));
+ fprintf(stderr, "agepart = %f, ageday = %f, moonfill = %f\n", agepart, ageday, moonfill);
+ */
+ return agepart;
+}
+
+/*
+int main(int argc, char * argv)
+{
+ int year, month, day, hour, min, zone;
+ double alt, azm, lon, lat, moonpart;
+ lon = 30.133333;
+ lat = 55.166666;
+ year = 2010;
+ month = 4;
+ day = 13;
+ hour = 18;
+ min = 10;
+ zone = 2;
+
+ //moonpart = moon_phase1(year, month, day);
+ //fprintf(stderr, "moon path = %f\n", moonpart);
+ alt = altitude(lon, lat, year, month, day, hour, min, zone);
+ azm = azimuth(lon, lat, year, month, day, hour, min, zone);
+ fprintf(stderr, "alt=%5.6f azm=%5.6f\n", alt, azm);
+
+ return 1;
+}
+*/