git.maemo.org Git - navit-package/blob - navit/sunriset.c

   1 /*
   2
   3 SUNRISET.C - computes Sun rise/set times, start/end of twilight, and
   4              the length of the day at any date and latitude
   5
   6 Written as DAYLEN.C, 1989-08-16
   7
   8 Modified to SUNRISET.C, 1992-12-01
   9
  10 (c) Paul Schlyter, 1989, 1992
  11
  12 Released to the public domain by Paul Schlyter, December 1992
  13
  14 */
  15
  16
  17 #include <stdio.h>
  18 #include <math.h>
  19
  20 #include "sunriset.h"
  21
  22 /* The "workhorse" function for sun rise/set times */
  23
  24 int __sunriset__( int year, int month, int day, double lon, double lat,
  25                   double altit, int upper_limb, double *trise, double *tset )
  26 /***************************************************************************/
  27 /* Note: year,month,date = calendar date, 1801-2099 only.             */
  28 /*       Eastern longitude positive, Western longitude negative       */
  29 /*       Northern latitude positive, Southern latitude negative       */
  30 /*       The longitude value IS critical in this function!            */
  31 /*       altit = the altitude which the Sun should cross              */
  32 /*               Set to -35/60 degrees for rise/set, -6 degrees       */
  33 /*               for civil, -12 degrees for nautical and -18          */
  34 /*               degrees for astronomical twilight.                   */
  35 /*         upper_limb: non-zero -> upper limb, zero -> center         */
  36 /*               Set to non-zero (e.g. 1) when computing rise/set     */
  37 /*               times, and to zero when computing start/end of       */
  38 /*               twilight.                                            */
  39 /*        *rise = where to store the rise time                        */
  40 /*        *set  = where to store the set  time                        */
  41 /*                Both times are relative to the specified altitude,  */
  42 /*                and thus this function can be used to comupte       */
  43 /*                various twilight times, as well as rise/set times   */
  44 /* Return value:  0 = sun rises/sets this day, times stored at        */
  45 /*                    *trise and *tset.                               */
  46 /*               +1 = sun above the specified "horizon" 24 hours.     */
  47 /*                    *trise set to time when the sun is at south,    */
  48 /*                    minus 12 hours while *tset is set to the south  */
  49 /*                    time plus 12 hours. "Day" length = 24 hours     */
  50 /*               -1 = sun is below the specified "horizon" 24 hours   */
  51 /*                    "Day" length = 0 hours, *trise and *tset are    */
  52 /*                    both set to the time when the sun is at south.  */
  53 /*                                                                    */
  54 /**********************************************************************/
  55 {
  56       double  d,  /* Days since 2000 Jan 0.0 (negative before) */
  57       sr,         /* Solar distance, astronomical units */
  58       sRA,        /* Sun's Right Ascension */
  59       sdec,       /* Sun's declination */
  60       sradius,    /* Sun's apparent radius */
  61       t,          /* Diurnal arc */
  62       tsouth,     /* Time when Sun is at south */
  63       sidtime;    /* Local sidereal time */
  64
  65       int rc = 0; /* Return cde from function - usually 0 */
  66
  67       /* Compute d of 12h local mean solar time */
  68       d = days_since_2000_Jan_0(year,month,day) + 0.5 - lon/360.0;
  69
  70       /* Compute local sideral time of this moment */
  71       sidtime = revolution( GMST0(d) + 180.0 + lon );
  72
  73       /* Compute Sun's RA + Decl at this moment */
  74       sun_RA_dec( d, &sRA, &sdec, &sr );
  75
  76       /* Compute time when Sun is at south - in hours UT */
  77       tsouth = 12.0 - rev180(sidtime - sRA)/15.0;
  78
  79       /* Compute the Sun's apparent radius, degrees */
  80       sradius = 0.2666 / sr;
  81
  82       /* Do correction to upper limb, if necessary */
  83       if ( upper_limb )
  84             altit -= sradius;
  85
  86       /* Compute the diurnal arc that the Sun traverses to reach */
  87       /* the specified altitide altit: */
  88       {
  89             double cost;
  90             cost = ( sind(altit) - sind(lat) * sind(sdec) ) /
  91                   ( cosd(lat) * cosd(sdec) );
  92             if ( cost >= 1.0 )
  93                   rc = -1, t = 0.0;       /* Sun always below altit */
  94             else if ( cost <= -1.0 )
  95                   rc = +1, t = 12.0;      /* Sun always above altit */
  96             else
  97                   t = acosd(cost)/15.0;   /* The diurnal arc, hours */
  98       }
  99
 100       /* Store rise and set times - in hours UT */
 101       *trise = tsouth - t;
 102       *tset  = tsouth + t;
 103
 104       return rc;
 105 }  /* __sunriset__ */
 106
 107
 108
 109 /* The "workhorse" function */
 110
 111
 112 double __daylen__( int year, int month, int day, double lon, double lat,
 113                    double altit, int upper_limb )
 114 /**********************************************************************/
 115 /* Note: year,month,date = calendar date, 1801-2099 only.             */
 116 /*       Eastern longitude positive, Western longitude negative       */
 117 /*       Northern latitude positive, Southern latitude negative       */
 118 /*       The longitude value is not critical. Set it to the correct   */
 119 /*       longitude if you're picky, otherwise set to to, say, 0.0     */
 120 /*       The latitude however IS critical - be sure to get it correct */
 121 /*       altit = the altitude which the Sun should cross              */
 122 /*               Set to -35/60 degrees for rise/set, -6 degrees       */
 123 /*               for civil, -12 degrees for nautical and -18          */
 124 /*               degrees for astronomical twilight.                   */
 125 /*         upper_limb: non-zero -> upper limb, zero -> center         */
 126 /*               Set to non-zero (e.g. 1) when computing day length   */
 127 /*               and to zero when computing day+twilight length.      */
 128 /**********************************************************************/
 129 {
 130       double  d,  /* Days since 2000 Jan 0.0 (negative before) */
 131       obl_ecl,    /* Obliquity (inclination) of Earth's axis */
 132       sr,         /* Solar distance, astronomical units */
 133       slon,       /* True solar longitude */
 134       sin_sdecl,  /* Sine of Sun's declination */
 135       cos_sdecl,  /* Cosine of Sun's declination */
 136       sradius,    /* Sun's apparent radius */
 137       t;          /* Diurnal arc */
 138
 139       /* Compute d of 12h local mean solar time */
 140       d = days_since_2000_Jan_0(year,month,day) + 0.5 - lon/360.0;
 141
 142       /* Compute obliquity of ecliptic (inclination of Earth's axis) */
 143       obl_ecl = 23.4393 - 3.563E-7 * d;
 144
 145       /* Compute Sun's position */
 146       sunpos( d, &slon, &sr );
 147
 148       /* Compute sine and cosine of Sun's declination */
 149       sin_sdecl = sind(obl_ecl) * sind(slon);
 150       cos_sdecl = sqrt( 1.0 - sin_sdecl * sin_sdecl );
 151
 152       /* Compute the Sun's apparent radius, degrees */
 153       sradius = 0.2666 / sr;
 154
 155       /* Do correction to upper limb, if necessary */
 156       if ( upper_limb )
 157             altit -= sradius;
 158
 159       /* Compute the diurnal arc that the Sun traverses to reach */
 160       /* the specified altitide altit: */
 161       {
 162             double cost;
 163             cost = ( sind(altit) - sind(lat) * sin_sdecl ) /
 164                   ( cosd(lat) * cos_sdecl );
 165             if ( cost >= 1.0 )
 166                   t = 0.0;                      /* Sun always below altit */
 167             else if ( cost <= -1.0 )
 168                   t = 24.0;                     /* Sun always above altit */
 169             else  t = (2.0/15.0) * acosd(cost); /* The diurnal arc, hours */
 170       }
 171       return t;
 172 }  /* __daylen__ */
 173
 174
 175 /* This function computes the Sun's position at any instant */
 176
 177 void sunpos( double d, double *lon, double *r )
 178 /******************************************************/
 179 /* Computes the Sun's ecliptic longitude and distance */
 180 /* at an instant given in d, number of days since     */
 181 /* 2000 Jan 0.0.  The Sun's ecliptic latitude is not  */
 182 /* computed, since it's always very near 0.           */
 183 /******************************************************/
 184 {
 185       double M,         /* Mean anomaly of the Sun */
 186              w,         /* Mean longitude of perihelion */
 187                         /* Note: Sun's mean longitude = M + w */
 188              e,         /* Eccentricity of Earth's orbit */
 189              E,         /* Eccentric anomaly */
 190              x, y,      /* x, y coordinates in orbit */
 191              v;         /* True anomaly */
 192
 193       /* Compute mean elements */
 194       M = revolution( 356.0470 + 0.9856002585 * d );
 195       w = 282.9404 + 4.70935E-5 * d;
 196       e = 0.016709 - 1.151E-9 * d;
 197
 198       /* Compute true longitude and radius vector */
 199       E = M + e * RADEG * sind(M) * ( 1.0 + e * cosd(M) );
 200       x = cosd(E) - e;
 201       y = sqrt( 1.0 - e*e ) * sind(E);
 202       *r = sqrt( x*x + y*y );              /* Solar distance */
 203       v = atan2d( y, x );                  /* True anomaly */
 204       *lon = v + w;                        /* True solar longitude */
 205       if ( *lon >= 360.0 )
 206             *lon -= 360.0;                   /* Make it 0..360 degrees */
 207 }
 208
 209 void sun_RA_dec( double d, double *RA, double *dec, double *r )
 210 {
 211   double lon, obl_ecl;
 212   double xs, ys, zs;
 213   double xe, ye, ze;
 214
 215   /* Compute Sun's ecliptical coordinates */
 216   sunpos( d, &lon, r );
 217
 218   /* Compute ecliptic rectangular coordinates */
 219   xs = *r * cosd(lon);
 220   ys = *r * sind(lon);
 221   zs = 0; /* because the Sun is always in the ecliptic plane! */
 222
 223   /* Compute obliquity of ecliptic (inclination of Earth's axis) */
 224   obl_ecl = 23.4393 - 3.563E-7 * d;
 225
 226   /* Convert to equatorial rectangular coordinates - x is unchanged */
 227   xe = xs;
 228   ye = ys * cosd(obl_ecl);
 229   ze = ys * sind(obl_ecl);
 230
 231   /* Convert to spherical coordinates */
 232   *RA = atan2d( ye, xe );
 233   *dec = atan2d( ze, sqrt(xe*xe + ye*ye) );
 234
 235 }  /* sun_RA_dec */
 236
 237
 238 /******************************************************************/
 239 /* This function reduces any angle to within the first revolution */
 240 /* by subtracting or adding even multiples of 360.0 until the     */
 241 /* result is >= 0.0 and < 360.0                                   */
 242 /******************************************************************/
 243
 244 #define INV360    ( 1.0 / 360.0 )
 245
 246 double revolution( double x )
 247 /*****************************************/
 248 /* Reduce angle to within 0..360 degrees */
 249 /*****************************************/
 250 {
 251       return( x - 360.0 * floor( x * INV360 ) );
 252 }  /* revolution */
 253
 254 double rev180( double x )
 255 /*********************************************/
 256 /* Reduce angle to within -180..+180 degrees */
 257 /*********************************************/
 258 {
 259       return( x - 360.0 * floor( x * INV360 + 0.5 ) );
 260 }  /* revolution */
 261
 262
 263 /*******************************************************************/
 264 /* This function computes GMST0, the Greenwhich Mean Sidereal Time */
 265 /* at 0h UT (i.e. the sidereal time at the Greenwhich meridian at  */
 266 /* 0h UT).  GMST is then the sidereal time at Greenwich at any     */
 267 /* time of the day.  I've generelized GMST0 as well, and define it */
 268 /* as:  GMST0 = GMST - UT  --  this allows GMST0 to be computed at */
 269 /* other times than 0h UT as well.  While this sounds somewhat     */
 270 /* contradictory, it is very practical:  instead of computing      */
 271 /* GMST like:                                                      */
 272 /*                                                                 */
 273 /*  GMST = (GMST0) + UT * (366.2422/365.2422)                      */
 274 /*                                                                 */
 275 /* where (GMST0) is the GMST last time UT was 0 hours, one simply  */
 276 /* computes:                                                       */
 277 /*                                                                 */
 278 /*  GMST = GMST0 + UT                                              */
 279 /*                                                                 */
 280 /* where GMST0 is the GMST "at 0h UT" but at the current moment!   */
 281 /* Defined in this way, GMST0 will increase with about 4 min a     */
 282 /* day.  It also happens that GMST0 (in degrees, 1 hr = 15 degr)   */
 283 /* is equal to the Sun's mean longitude plus/minus 180 degrees!    */
 284 /* (if we neglect aberration, which amounts to 20 seconds of arc   */
 285 /* or 1.33 seconds of time)                                        */
 286 /*                                                                 */
 287 /*******************************************************************/
 288
 289 double GMST0( double d )
 290 {
 291       double sidtim0;
 292       /* Sidtime at 0h UT = L (Sun's mean longitude) + 180.0 degr  */
 293       /* L = M + w, as defined in sunpos().  Since I'm too lazy to */
 294       /* add these numbers, I'll let the C compiler do it for me.  */
 295       /* Any decent C compiler will add the constants at compile   */
 296       /* time, imposing no runtime or code overhead.               */
 297       sidtim0 = revolution( ( 180.0 + 356.0470 + 282.9404 ) +
 298                           ( 0.9856002585 + 4.70935E-5 ) * d );
 299       return sidtim0;
 300 }  /* GMST0 */