vcs.maemo.org Git - livewp/blob - applet/src/livewp-astro.c

   1 /*
   2  * This file is part of Live Wallpaper (livewp)
   3  *
   4  * Copyright (C) 2010 Vlad Vasiliev
   5  * Copyright (C) 2010 Tanya Makova
   6  *       for the code
   7  *
   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.
  12  *
  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.
  17  *
  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
  21  * 02110-1301 USA
  22 */
  23 /*******************************************************************************/
  24 #include <stdio.h>
  25 #include "livewp-astro.h"
  26
  27 double ut(int hour, int min, int zone)
  28 {
  29     return (hour - zone + (double)min/60);
  30 }
  31
  32 double jd(int year, int month, int day, double uu)
  33 {
  34     return (367*year)-floor((floor((month+9)/12)+year)*7/4)+floor(275*month/9)+day-730531.5+(uu/24);
  35 }
  36
  37 double altitude(double lon, double lat, int year, int month, int day, int hour, int min, int zone)
  38 {
  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);
  42     T = jj / 36525;
  43     k = M_PI / 180;
  44     M = 357.5291 + 35999.0503 * T - 0.0001559 * T * T - 0.00000045 * T * T * T;
  45     M = fmod(M, 360);
  46     Lo = 280.46645 + 36000.76983 * T + 0.0003032 * T * T;
  47     Lo = fmod(Lo, 360);
  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);
  49     L=Lo+DL;
  50     eps=23.43999-0.013*T;
  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);
  56     LMST=GMST+lon;
  57     H=LMST-RA;
  58     eqt=(Lo-RA)*4;
  59     alt=(1/k)*asin(sin(lat*k)*sin(delta*k)+cos(lat*k)*cos(delta*k)*cos(H*k));
  60     return alt;
  61 }
  62
  63 double azimuth(double lon, double lat, int year, int month, int day, int hour, int min, int zone)
  64 {
  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);
  68     T = jj / 36525;
  69     k = M_PI / 180;
  70     M = 357.5291 + 35999.0503 * T - 0.0001559 * T * T - 0.00000045 * T * T * T;
  71     M = fmod(M, 360);
  72     Lo = 280.46645 + 36000.76983 * T + 0.0003032 * T * T;
  73     Lo = fmod(Lo, 360);
  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);
  75     L=Lo+DL;
  76     eps=23.43999-0.013*T;
  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);
  82     LMST=GMST+lon;
  83     H=LMST-RA;
  84     eqt=(Lo-RA)*4;
  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);
  87     return azm;
  88 }
  89 double moon_phase(int year, int month, int day)
  90 {
  91     int lg;
  92     double agepart, ge, jd, moonday;
  93     ge = 1721425.5;
  94     moonday = 29.530588853;
  95     if (month <= 2) lg = 0;
  96     else {
  97         if ((fmod(year, 4) == 0) && (!(fmod(year, 100) == 0) && (fmod(year, 400) != 0)))
  98             lg = -1;
  99         else lg = -2;
 100     }
 101     jd = (ge - 1) +
 102         (365 * (year - 1)) +
 103         floor((year - 1) / 4) +
 104         (-floor((year - 1) / 100)) +
 105         floor((year - 1) / 400) +
 106         floor((((367 * month) - 362) / 12) + lg + day);
 107
 108     agepart = (jd - 2451550.1) / moonday;
 109     agepart = agepart - floor(agepart);
 110     if (agepart < 0) agepart + 1;
 111 /*
 112     ageday = agepart * moonday;
 113
 114     moonfill = 1.0 - ((moonday/2 - ageday) / (moonday/2));
 115     fprintf(stderr, "agepart = %f, ageday = %f, moonfill = %f\n", agepart, ageday, moonfill);
 116     */
 117     return agepart;
 118 }
 119
 120 /*
 121 int main(int argc, char * argv)
 122 {
 123     int year, month, day, hour, min, zone;
 124     double alt, azm, lon, lat, moonpart;
 125     lon = 30.133333;
 126     lat = 55.166666;
 127     year = 2010;
 128     month = 4;
 129     day = 13;
 130     hour = 18;
 131     min = 10;
 132     zone = 2;
 133
 134     //moonpart = moon_phase1(year, month, day);
 135     //fprintf(stderr, "moon path = %f\n", moonpart);
 136     alt = altitude(lon, lat, year, month, day, hour, min, zone);
 137     azm = azimuth(lon, lat, year, month, day, hour, min, zone);
 138     fprintf(stderr, "alt=%5.6f azm=%5.6f\n", alt, azm);
 139
 140     return 1;
 141 }
 142 */