vcs.maemo.org Git - livewp/blob - applet/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     fprintf(stderr, "%f, %f, %i.%i.%i %i:%i %i uu=%f", lon, lat, year, month, day, hour, min, zone, uu);
  68     jj = jd(year, month, day, uu);
  69     T = jj / 36525;
  70     k = M_PI / 180;
  71     M = 357.5291 + 35999.0503 * T - 0.0001559 * T * T - 0.00000045 * T * T * T;
  72     M = fmod(M, 360);
  73     Lo = 280.46645 + 36000.76983 * T + 0.0003032 * T * T;
  74     Lo = fmod(Lo, 360);
  75     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);
  76     L=Lo+DL;
  77     eps=23.43999-0.013*T;
  78     delta=(1/k)*asin(sin(L*k)*sin(eps*k));
  79     RA=(1/k)*atan2(cos(eps*k)*sin(L*k),cos(L*k));
  80     RA=fmod((RA+360), 360);
  81     GMST=280.46061837+360.98564736629*jj+0.000387933*T*T-T*T*T/38710000;
  82     GMST=fmod((GMST+360), 360);
  83     LMST=GMST+lon;
  84     H=LMST-RA;
  85     eqt=(Lo-RA)*4;
  86     azm=(1/k)*atan2(-sin(H*k),cos(lat*k)*tan(delta*k)-sin(lat*k)*cos(H*k));
  87     azm=fmod((azm+360), 360);
  88     return azm;
  89 }
  90 double moon_phase(int year, int month, int day)
  91 {
  92     int yy, mm, dd, k1, k2, k3, lg;
  93     double agepart, ageday, ge, moonday, jd, moonfill;
  94     ge = 1721425.5;
  95     moonday = 29.530588853;
  96     if (month <= 2) lg = 0;
  97     else {
  98         if ((fmod(year, 4) == 0) && (!(fmod(year, 100) == 0) && (fmod(year, 400) != 0)))
  99             lg = -1;
 100         else lg = -2;
 101     }
 102     jd = (ge - 1) +
 103         (365 * (year - 1)) +
 104         floor((year - 1) / 4) +
 105         (-floor((year - 1) / 100)) +
 106         floor((year - 1) / 400) +
 107         floor((((367 * month) - 362) / 12) + lg + day);
 108
 109     agepart = (jd - 2451550.1) / moonday;
 110     agepart = agepart - floor(agepart);
 111     if (agepart < 0) agepart + 1;
 112 /*
 113     ageday = agepart * moonday;
 114
 115     moonfill = 1.0 - ((moonday/2 - ageday) / (moonday/2));
 116     fprintf(stderr, "agepart = %f, ageday = %f, moonfill = %f\n", agepart, ageday, moonfill);
 117     */
 118     return agepart;
 119 }
 120
 121 /*
 122 int main(int argc, char * argv)
 123 {
 124     int year, month, day, hour, min, zone;
 125     double alt, azm, lon, lat, moonpart;
 126     lon = 30.133333;
 127     lat = 55.166666;
 128     year = 2010;
 129     month = 4;
 130     day = 13;
 131     hour = 18;
 132     min = 10;
 133     zone = 2;
 134
 135     //moonpart = moon_phase1(year, month, day);
 136     //fprintf(stderr, "moon path = %f\n", moonpart);
 137     alt = altitude(lon, lat, year, month, day, hour, min, zone);
 138     azm = azimuth(lon, lat, year, month, day, hour, min, zone);
 139     fprintf(stderr, "alt=%5.6f azm=%5.6f\n", alt, azm);
 140
 141     return 1;
 142 }
 143 */