wmglobe/src/sunpos.c

   1 /*     WMGlobe 1.3  -  All the Earth on a WMaker Icon
   2  *     copyright (C) 1998,99,2000,01 Jerome Dumonteil <jerome.dumonteil@linuxfr.org>
   3  * sunpos.c is taken from Xearth 1.0 (and part of 1.1):
   4  */
   5 /*
   6  * sunpos.c
   7  * kirk johnson
   8  * july 1993
   9  *
  10  * code for calculating the position on the earth's surface for which
  11  * the sun is directly overhead (adapted from _practical astronomy
  12  * with your calculator, third edition_, peter duffett-smith,
  13  * cambridge university press, 1988.)
  14  *
  15  *
  16  * Copyright (C) 1989, 1990, 1993, 1994, 1995 Kirk Lauritz Johnson
  17  *
  18  * Parts of the source code (as marked) are:
  19  *   Copyright (C) 1989, 1990, 1991 by Jim Frost
  20  *   Copyright (C) 1992 by Jamie Zawinski <jwz@lucid.com>
  21  *
  22  * Permission to use, copy, modify and freely distribute xearth for
  23  * non-commercial and not-for-profit purposes is hereby granted
  24  * without fee, provided that both the above copyright notice and this
  25  * permission notice appear in all copies and in supporting
  26  * documentation.
  27  *
  28  * Unisys Corporation holds worldwide patent rights on the Lempel Zev
  29  * Welch (LZW) compression technique employed in the CompuServe GIF
  30  * image file format as well as in other formats. Unisys has made it
  31  * clear, however, that it does not require licensing or fees to be
  32  * paid for freely distributed, non-commercial applications (such as
  33  * xearth) that employ LZW/GIF technology. Those wishing further
  34  * information about licensing the LZW patent should contact Unisys
  35  * directly at (lzw_info@unisys.com) or by writing to
  36  *
  37  *   Unisys Corporation
  38  *   Welch Licensing Department
  39  *   M/S-C1SW19
  40  *   P.O. Box 500
  41  *   Blue Bell, PA 19424
  42  *
  43  * The author makes no representations about the suitability of this
  44  * software for any purpose. It is provided "as is" without express or
  45  * implied warranty.
  46  *
  47  * THE AUTHOR DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
  48  * INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS,
  49  * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY SPECIAL, INDIRECT
  50  * OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM
  51  * LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT,
  52  * NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN
  53  * CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
  54  */
  55
  56 /*************************************************************************/
  57
  58 #include <math.h>
  59 #include <time.h>
  60
  61 #ifndef PI
  62 #define PI 3.141592653
  63 #endif
  64 #define TWOPI (2*PI)
  65 #define DegsToRads(x) ((x)*(TWOPI/360))
  66
  67 /*
  68  * the epoch upon which these astronomical calculations are based is
  69  * 1990 january 0.0, 631065600 seconds since the beginning of the
  70  * "unix epoch" (00:00:00 GMT, Jan. 1, 1970)
  71  *
  72  * given a number of seconds since the start of the unix epoch,
  73  * DaysSinceEpoch() computes the number of days since the start of the
  74  * astronomical epoch (1990 january 0.0)
  75  */
  76
  77 #define EpochStart           (631065600)
  78 #define DaysSinceEpoch(secs) (((secs)-EpochStart)*(1.0/(24*3600)))
  79
  80 /*
  81  * assuming the apparent orbit of the sun about the earth is circular,
  82  * the rate at which the orbit progresses is given by RadsPerDay --
  83  * TWOPI radians per orbit divided by 365.242191 days per year:
  84  */
  85
  86 #define RadsPerDay (TWOPI/365.242191)
  87
  88 /*
  89  * details of sun's apparent orbit at epoch 1990.0 (after
  90  * duffett-smith, table 6, section 46)
  91  *
  92  * Epsilon_g    (ecliptic longitude at epoch 1990.0) 279.403303 degrees
  93  * OmegaBar_g   (ecliptic longitude of perigee)      282.768422 degrees
  94  * Eccentricity (eccentricity of orbit)                0.016713
  95  */
  96
  97 #define Epsilon_g    (DegsToRads(279.403303))
  98 #define OmegaBar_g   (DegsToRads(282.768422))
  99 #define Eccentricity (0.016713)
 100
 101 /*
 102  * MeanObliquity gives the mean obliquity of the earth's axis at epoch
 103  * 1990.0 (computed as 23.440592 degrees according to the method given
 104  * in duffett-smith, section 27)
 105  */
 106 #define MeanObliquity (DegsToRads(23.440592))
 107
 108 /*
 109  * Lunar parameters, epoch January 0, 1990.0 (from Xearth 1.1)
 110  */
 111 #define MoonMeanLongitude        DegsToRads(318.351648)
 112 #define MoonMeanLongitudePerigee DegsToRads( 36.340410)
 113 #define MoonMeanLongitudeNode    DegsToRads(318.510107)
 114 #define MoonInclination          DegsToRads(  5.145396)
 115
 116 #define SideralMonth (27.3217)
 117
 118 /*
 119  * Force an angular value into the range [-PI, +PI]
 120  */
 121 #define Normalize(x)                        \
 122   do {                                      \
 123     if ((x) < -PI)                          \
 124       do (x) += TWOPI; while ((x) < -PI);   \
 125     else if ((x) > PI)                      \
 126       do (x) -= TWOPI; while ((x) > PI);    \
 127   } while (0)
 128
 129 static double solve_keplers_equation(double M);
 130 static double mean_sun(double D);
 131 static double sun_ecliptic_longitude(time_t ssue);
 132 static void ecliptic_to_equatorial(double lambda, double beta,
 133                                    double *alpha, double *delta);
 134 static double julian_date(int y, int m, int d);
 135 static double GST(time_t ssue);
 136
 137 /*
 138  * solve Kepler's equation via Newton's method
 139  * (after duffett-smith, section 47)
 140  */
 141 static double solve_keplers_equation(double M)
 142 {
 143     double E;
 144     double delta;
 145
 146     E = M;
 147     while (1) {
 148         delta = E - Eccentricity * sin(E) - M;
 149         if (fabs(delta) <= 1e-10)
 150             break;
 151         E -= delta / (1 - Eccentricity * cos(E));
 152     }
 153     return E;
 154 }
 155
 156 /*
 157  * Calculate the position of the mean sun: where the sun would
 158  * be if the earth's orbit were circular instead of ellipictal.
 159  */
 160
 161 static double mean_sun(double D)
 162 /*     double D;                   days since ephemeris epoch */
 163 {
 164     double N, M;
 165
 166     N = RadsPerDay * D;
 167     N = fmod(N, TWOPI);
 168     if (N < 0)
 169         N += TWOPI;
 170
 171     M = N + Epsilon_g - OmegaBar_g;
 172     if (M < 0)
 173         M += TWOPI;
 174     return M;
 175 }
 176
 177 /*
 178  * compute ecliptic longitude of sun (in radians)
 179  * (after duffett-smith, section 47)
 180  */
 181 static double sun_ecliptic_longitude(time_t ssue)
 182                                 /* seconds since unix epoch */
 183 {
 184     double D, N;
 185     double M_sun, E;
 186     double v;
 187
 188     D = DaysSinceEpoch(ssue);
 189
 190     N = RadsPerDay * D;
 191     M_sun = mean_sun(D);
 192     E = solve_keplers_equation(M_sun);
 193     v = 2 * atan(sqrt((1 + Eccentricity) / (1 - Eccentricity)) *
 194                  tan(E / 2));
 195
 196     return (v + OmegaBar_g);
 197 }
 198
 199
 200 /*
 201  * convert from ecliptic to equatorial coordinates
 202  * (after duffett-smith, section 27)
 203  */
 204 static void ecliptic_to_equatorial(double lambda, double beta,
 205                                    double *alpha, double *delta)
 206 /*
 207  *    double  lambda;            ecliptic longitude
 208  *    double  beta;               ecliptic latitude
 209  *    double *alpha;              (return) right ascension
 210  *    double *delta;           (return) declination
 211  */
 212 {
 213     double sin_e, cos_e;
 214
 215     sin_e = sin(MeanObliquity);
 216     cos_e = cos(MeanObliquity);
 217
 218     *alpha = atan2(sin(lambda) * cos_e - tan(beta) * sin_e, cos(lambda));
 219     *delta = asin(sin(beta) * cos_e + cos(beta) * sin_e * sin(lambda));
 220 }
 221
 222
 223 /*
 224  * computing julian dates (assuming gregorian calendar, thus this is
 225  * only valid for dates of 1582 oct 15 or later)
 226  * (after duffett-smith, section 4)
 227  */
 228 static double julian_date(int y, int m, int d)
 229 /*
 230  *     int y;                     year (e.g. 19xx)
 231  *     int m;                      month (jan=1, feb=2, ...)
 232  *     int d;                      day of month
 233  */
 234 {
 235     int A, B, C, D;
 236     double JD;
 237
 238     /* lazy test to ensure gregorian calendar */
 239 /*
 240  *   ASSERT(y >= 1583);
 241  */
 242
 243     if ((m == 1) || (m == 2)) {
 244         y -= 1;
 245         m += 12;
 246     }
 247     A = y / 100;
 248     B = 2 - A + (A / 4);
 249     C = (int) (365.25 * y);
 250     D = (int) (30.6001 * (m + 1));
 251
 252     JD = B + C + D + d + 1720994.5;
 253
 254     return JD;
 255 }
 256
 257
 258 /*
 259  * compute greenwich mean sidereal time (GST) corresponding to a given
 260  * number of seconds since the unix epoch
 261  * (after duffett-smith, section 12)
 262  */
 263 static double GST(time_t ssue)
 264    /*time_t ssue;          seconds since unix epoch */
 265 {
 266     double JD;
 267     double T, T0;
 268     double UT;
 269     struct tm *tm;
 270
 271     tm = gmtime(&ssue);
 272
 273     JD = julian_date(tm->tm_year + 1900, tm->tm_mon + 1, tm->tm_mday);
 274     T = (JD - 2451545) / 36525;
 275
 276     T0 = ((T + 2.5862e-5) * T + 2400.051336) * T + 6.697374558;
 277
 278     T0 = fmod(T0, 24.0);
 279     if (T0 < 0)
 280         T0 += 24;
 281
 282     UT = tm->tm_hour + (tm->tm_min + tm->tm_sec / 60.0) / 60.0;
 283
 284     T0 += UT * 1.002737909;
 285     T0 = fmod(T0, 24.0);
 286     if (T0 < 0)
 287         T0 += 24;
 288
 289     return T0;
 290 }
 291
 292
 293 /*
 294  * given a particular time (expressed in seconds since the unix
 295  * epoch), compute position on the earth (lat, lon) such that sun is
 296  * directly overhead.
 297  */
 298 void sun_position(time_t ssue, double *lat, double *lon)
 299 /*     time_t  ssue;               seconds since unix epoch    */
 300 /*     double *lat;                (return) latitude  in rad   */
 301 /*     double *lon;                (return) longitude in rad   */
 302 {
 303     double lambda;
 304     double alpha, delta;
 305     double tmp;
 306
 307     lambda = sun_ecliptic_longitude(ssue);
 308     ecliptic_to_equatorial(lambda, 0.0, &alpha, &delta);
 309
 310     tmp = alpha - (TWOPI / 24) * GST(ssue);
 311     Normalize(tmp);
 312     *lon = tmp;
 313     *lat = delta;
 314 }
 315
 316 /*
 317  * given a particular time (expressed in seconds since the unix
 318  * epoch), compute position on the earth (lat, lon) such that the
 319  * moon is directly overhead.
 320  *
 321  * Based on duffett-smith **2nd ed** section 61; combines some steps
 322  * into single expressions to reduce the number of extra variables.
 323  */
 324 void moon_position(time_t ssue, double *lat, double *lon)
 325 /*    time_t  ssue;              seconds since unix epoch    */
 326 /*   double *lat;                (return) latitude   in ra   */
 327 /*  double *lon;                 (return) longitude   in ra  */
 328 {
 329     double lambda, beta;
 330     double D, L, Ms, Mm, N, Ev, Ae, Ec, alpha, delta;
 331
 332     D = DaysSinceEpoch(ssue);
 333     lambda = sun_ecliptic_longitude(ssue);
 334     Ms = mean_sun(D);
 335
 336     L = fmod(D / SideralMonth, 1.0) * TWOPI + MoonMeanLongitude;
 337     Normalize(L);
 338     Mm = L - DegsToRads(0.1114041 * D) - MoonMeanLongitudePerigee;
 339     Normalize(Mm);
 340     N = MoonMeanLongitudeNode - DegsToRads(0.0529539 * D);
 341     Normalize(N);
 342     Ev = DegsToRads(1.2739) * sin(2.0 * (L - lambda) - Mm);
 343     Ae = DegsToRads(0.1858) * sin(Ms);
 344     Mm += Ev - Ae - DegsToRads(0.37) * sin(Ms);
 345     Ec = DegsToRads(6.2886) * sin(Mm);
 346     L += Ev + Ec - Ae + DegsToRads(0.214) * sin(2.0 * Mm);
 347     L += DegsToRads(0.6583) * sin(2.0 * (L - lambda));
 348     N -= DegsToRads(0.16) * sin(Ms);
 349
 350     L -= N;
 351     lambda = (fabs(cos(L)) < 1e-12) ?
 352         (N + sin(L) * cos(MoonInclination) * PI / 2) :
 353         (N + atan2(sin(L) * cos(MoonInclination), cos(L)));
 354     Normalize(lambda);
 355     beta = asin(sin(L) * sin(MoonInclination));
 356     ecliptic_to_equatorial(lambda, beta, &alpha, &delta);
 357     alpha -= (TWOPI / 24) * GST(ssue);
 358     Normalize(alpha);
 359     *lon = alpha;
 360     *lat = delta;
 361 }