vcs.maemo.org Git - livewp/blob - applet/astrolib.c

   1 #include <stdio.h>
   2 #include <math.h>
   3
   4 double ut(int hour, int min, int zone)
   5 {
   6     return (hour - zone + (double)min/60);
   7 }
   8
   9 double jd(int year, int month, int day, double uu)
  10 {
  11     return (367*year)-floor((floor((month+9)/12)+year)*7/4)+floor(275*month/9)+day-730531.5+(uu/24);
  12 }
  13
  14 double altitude(double lon, double lat, int year, int month, int day, int hour, int min, int zone)
  15 {
  16     double uu, jj, T, k, M, Lo, DL, L, eps, delta, RA, GMST, LMST, H, eqt, alt;
  17     uu = ut(hour, min, zone);
  18     jj = jd(year, month, day, uu);
  19     T = jj / 36525;
  20     k = M_PI / 180;
  21     M = 357.5291 + 35999.0503 * T - 0.0001559 * T * T - 0.00000045 * T * T * T;
  22     M = fmod(M, 360);
  23     Lo = 280.46645 + 36000.76983 * T + 0.0003032 * T * T;
  24     Lo = fmod(Lo, 360);
  25     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);
  26     L=Lo+DL;
  27     eps=23.43999-0.013*T;
  28     delta=(1/k)*asin(sin(L*k)*sin(eps*k));
  29     RA=(1/k)*atan2(cos(eps*k)*sin(L*k),cos(L*k));
  30     RA=fmod((RA+360), 360);
  31     GMST=280.46061837+360.98564736629*jj+0.000387933*T*T-T*T*T/38710000;
  32     GMST=fmod((GMST+360), 360);
  33     LMST=GMST+lon;
  34     H=LMST-RA;
  35     eqt=(Lo-RA)*4;
  36     alt=(1/k)*asin(sin(lat*k)*sin(delta*k)+cos(lat*k)*cos(delta*k)*cos(H*k));
  37     return alt;
  38 }
  39
  40 double azimuth(double lon, double lat, int year, int month, int day, int hour, int min, int zone)
  41 {
  42     double uu, jj, T, k, M, Lo, DL, L, eps, delta, RA, GMST, LMST, H, eqt, azm;
  43     uu = ut(hour, min, zone);
  44     fprintf(stderr, "%f, %f, %i.%i.%i %i:%i %i uu=%f", lon, lat, year, month, day, hour, min, zone, uu);
  45     jj = jd(year, month, day, uu);
  46     T = jj / 36525;
  47     k = M_PI / 180;
  48     M = 357.5291 + 35999.0503 * T - 0.0001559 * T * T - 0.00000045 * T * T * T;
  49     M = fmod(M, 360);
  50     Lo = 280.46645 + 36000.76983 * T + 0.0003032 * T * T;
  51     Lo = fmod(Lo, 360);
  52     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);
  53     L=Lo+DL;
  54     eps=23.43999-0.013*T;
  55     delta=(1/k)*asin(sin(L*k)*sin(eps*k));
  56     RA=(1/k)*atan2(cos(eps*k)*sin(L*k),cos(L*k));
  57     RA=fmod((RA+360), 360);
  58     GMST=280.46061837+360.98564736629*jj+0.000387933*T*T-T*T*T/38710000;
  59     GMST=fmod((GMST+360), 360);
  60     LMST=GMST+lon;
  61     H=LMST-RA;
  62     eqt=(Lo-RA)*4;
  63     azm=(1/k)*atan2(-sin(H*k),cos(lat*k)*tan(delta*k)-sin(lat*k)*cos(H*k));
  64     azm=fmod((azm+360), 360);
  65     return azm;
  66 }
  67 double moon_phase(int year, int month, int day)
  68 {
  69     int yy, mm, dd, k1, k2, k3, jd;
  70     double ip, v, ag;
  71     yy = year - floor((12 - month) / 10);
  72     mm = month + 9;
  73     if (mm >= 12) mm = mm - 12;
  74
  75     k1 = floor(365.35 * (yy + 4712));
  76     k2 = floor(30.6 * mm + 0.05);
  77     k3 = floor(floor(yy/100+49) * 0.75) - 38;
  78
  79     jd = k1 + k2 + day + 59;
  80     if (jd > 2299160) jd = jd - k3;
  81
  82     v = (jd - 2451550.1) / 29.530588853;
  83     ip = v - floor(v);
  84     if (ip < 0) ip = ip + 1;
  85     ag = ip * 29.53;
  86     fprintf(stderr, "ip = %f, ag = %f\n", ip, ag);
  87 }
  88
  89 double moon_phase1(int year, int month, int day)
  90 {
  91     int yy, mm, dd, k1, k2, k3, lg;
  92     double agepart, ageday, ge, moonday, jd, moonfill;
  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 }*/