Back to index

php5  5.3.10
astro.c
Go to the documentation of this file.
00001 /*
00002    +----------------------------------------------------------------------+
00003    | PHP Version 5                                                        |
00004    +----------------------------------------------------------------------+
00005    | Copyright (c) 1997-2010 The PHP Group                                |
00006    +----------------------------------------------------------------------+
00007    | This source file is subject to version 3.01 of the PHP license,      |
00008    | that is bundled with this package in the file LICENSE, and is        |
00009    | available through the world-wide-web at the following url:           |
00010    | http://www.php.net/license/3_01.txt                                  |
00011    | If you did not receive a copy of the PHP license and are unable to   |
00012    | obtain it through the world-wide-web, please send a note to          |
00013    | license@php.net so we can mail you a copy immediately.               |
00014    +----------------------------------------------------------------------+
00015    | Algorithms are taken from a public domain source by Paul             |
00016    | Schlyter, who wrote this in December 1992                            |
00017    +----------------------------------------------------------------------+
00018    | Authors: Derick Rethans <derick@derickrethans.nl>                    |
00019    +----------------------------------------------------------------------+
00020  */
00021 
00022 /* $Id: astro.c 293036 2010-01-03 09:23:27Z sebastian $ */
00023 
00024 #include <stdio.h>
00025 #include <math.h>
00026 #include "timelib.h"
00027 
00028 #define days_since_2000_Jan_0(y,m,d) \
00029        (367L*(y)-((7*((y)+(((m)+9)/12)))/4)+((275*(m))/9)+(d)-730530L)
00030 
00031 #ifndef PI
00032  #define PI        3.1415926535897932384
00033 #endif
00034 
00035 #define RADEG     ( 180.0 / PI )
00036 #define DEGRAD    ( PI / 180.0 )
00037 
00038 /* The trigonometric functions in degrees */
00039 
00040 #define sind(x)  sin((x)*DEGRAD)
00041 #define cosd(x)  cos((x)*DEGRAD)
00042 #define tand(x)  tan((x)*DEGRAD)
00043 
00044 #define atand(x)    (RADEG*atan(x))
00045 #define asind(x)    (RADEG*asin(x))
00046 #define acosd(x)    (RADEG*acos(x))
00047 #define atan2d(y,x) (RADEG*atan2(y,x))
00048 
00049 
00050 /* Following are some macros around the "workhorse" function __daylen__ */
00051 /* They mainly fill in the desired values for the reference altitude    */
00052 /* below the horizon, and also selects whether this altitude should     */
00053 /* refer to the Sun's center or its upper limb.                         */
00054 
00055 
00056 #include "astro.h"
00057 
00058 /******************************************************************/
00059 /* This function reduces any angle to within the first revolution */
00060 /* by subtracting or adding even multiples of 360.0 until the     */
00061 /* result is >= 0.0 and < 360.0                                   */
00062 /******************************************************************/
00063 
00064 #define INV360    (1.0 / 360.0)
00065 
00066 /*****************************************/
00067 /* Reduce angle to within 0..360 degrees */
00068 /*****************************************/
00069 static double astro_revolution(double x)
00070 {
00071        return (x - 360.0 * floor(x * INV360));
00072 }
00073 
00074 /*********************************************/
00075 /* Reduce angle to within +180..+180 degrees */
00076 /*********************************************/
00077 static double astro_rev180( double x )
00078 {
00079        return (x - 360.0 * floor(x * INV360 + 0.5));
00080 }
00081 
00082 /*******************************************************************/
00083 /* This function computes GMST0, the Greenwich Mean Sidereal Time  */
00084 /* at 0h UT (i.e. the sidereal time at the Greenwhich meridian at  */
00085 /* 0h UT).  GMST is then the sidereal time at Greenwich at any     */
00086 /* time of the day.  I've generalized GMST0 as well, and define it */
00087 /* as:  GMST0 = GMST - UT  --  this allows GMST0 to be computed at */
00088 /* other times than 0h UT as well.  While this sounds somewhat     */
00089 /* contradictory, it is very practical:  instead of computing      */
00090 /* GMST like:                                                      */
00091 /*                                                                 */
00092 /*  GMST = (GMST0) + UT * (366.2422/365.2422)                      */
00093 /*                                                                 */
00094 /* where (GMST0) is the GMST last time UT was 0 hours, one simply  */
00095 /* computes:                                                       */
00096 /*                                                                 */
00097 /*  GMST = GMST0 + UT                                              */
00098 /*                                                                 */
00099 /* where GMST0 is the GMST "at 0h UT" but at the current moment!   */
00100 /* Defined in this way, GMST0 will increase with about 4 min a     */
00101 /* day.  It also happens that GMST0 (in degrees, 1 hr = 15 degr)   */
00102 /* is equal to the Sun's mean longitude plus/minus 180 degrees!    */
00103 /* (if we neglect aberration, which amounts to 20 seconds of arc   */
00104 /* or 1.33 seconds of time)                                        */
00105 /*                                                                 */
00106 /*******************************************************************/
00107 
00108 static double astro_GMST0(double d)
00109 {
00110        double sidtim0;
00111        /* Sidtime at 0h UT = L (Sun's mean longitude) + 180.0 degr  */
00112        /* L = M + w, as defined in sunpos().  Since I'm too lazy to */
00113        /* add these numbers, I'll let the C compiler do it for me.  */
00114        /* Any decent C compiler will add the constants at compile   */
00115        /* time, imposing no runtime or code overhead.               */
00116        sidtim0 = astro_revolution((180.0 + 356.0470 + 282.9404) + (0.9856002585 + 4.70935E-5) * d);
00117        return sidtim0;
00118 } 
00119 
00120 /* This function computes the Sun's position at any instant */
00121 
00122 /******************************************************/
00123 /* Computes the Sun's ecliptic longitude and distance */
00124 /* at an instant given in d, number of days since     */
00125 /* 2000 Jan 0.0.  The Sun's ecliptic latitude is not  */
00126 /* computed, since it's always very near 0.           */
00127 /******************************************************/
00128 static void astro_sunpos(double d, double *lon, double *r)
00129 {
00130        double M,         /* Mean anomaly of the Sun */
00131               w,         /* Mean longitude of perihelion */
00132                          /* Note: Sun's mean longitude = M + w */
00133               e,         /* Eccentricity of Earth's orbit */
00134               E,         /* Eccentric anomaly */
00135               x, y,      /* x, y coordinates in orbit */
00136               v;         /* True anomaly */
00137 
00138        /* Compute mean elements */
00139        M = astro_revolution(356.0470 + 0.9856002585 * d);
00140        w = 282.9404 + 4.70935E-5 * d;
00141        e = 0.016709 - 1.151E-9 * d;
00142 
00143        /* Compute true longitude and radius vector */
00144        E = M + e * RADEG * sind(M) * (1.0 + e * cosd(M));
00145        x = cosd(E) - e;
00146        y = sqrt(1.0 - e*e) * sind(E);
00147        *r = sqrt(x*x + y*y);              /* Solar distance */
00148        v = atan2d(y, x);                  /* True anomaly */
00149        *lon = v + w;                        /* True solar longitude */
00150        if (*lon >= 360.0) {
00151               *lon -= 360.0;                   /* Make it 0..360 degrees */
00152        }
00153 }
00154 
00155 static void astro_sun_RA_dec(double d, double *RA, double *dec, double *r)
00156 {
00157        double lon, obl_ecl, x, y, z;
00158 
00159        /* Compute Sun's ecliptical coordinates */
00160        astro_sunpos(d, &lon, r);
00161 
00162        /* Compute ecliptic rectangular coordinates (z=0) */
00163        x = *r * cosd(lon);
00164        y = *r * sind(lon);
00165 
00166        /* Compute obliquity of ecliptic (inclination of Earth's axis) */
00167        obl_ecl = 23.4393 - 3.563E-7 * d;
00168 
00169        /* Convert to equatorial rectangular coordinates - x is unchanged */
00170        z = y * sind(obl_ecl);
00171        y = y * cosd(obl_ecl);
00172 
00173        /* Convert to spherical coordinates */
00174        *RA = atan2d(y, x);
00175        *dec = atan2d(z, sqrt(x*x + y*y));
00176 }
00177 
00207 int timelib_astro_rise_set_altitude(timelib_time *t_loc, double lon, double lat, double altit, int upper_limb, double *h_rise, double *h_set, timelib_sll *ts_rise, timelib_sll *ts_set, timelib_sll *ts_transit)
00208 {
00209        double  d,  /* Days since 2000 Jan 0.0 (negative before) */
00210        sr,         /* Solar distance, astronomical units */
00211        sRA,        /* Sun's Right Ascension */
00212        sdec,       /* Sun's declination */
00213        sradius,    /* Sun's apparent radius */
00214        t,          /* Diurnal arc */
00215        tsouth,     /* Time when Sun is at south */
00216        sidtime;    /* Local sidereal time */
00217        timelib_time *t_utc;
00218        timelib_sll   timestamp, old_sse;
00219 
00220        int rc = 0; /* Return cde from function - usually 0 */
00221 
00222        /* Normalize time */
00223        old_sse = t_loc->sse;
00224        t_loc->h = 12;
00225        t_loc->i = t_loc->s = 0;
00226        timelib_update_ts(t_loc, NULL);
00227 
00228        /* Calculate TS belonging to UTC 00:00 of the current day */
00229        t_utc = timelib_time_ctor();
00230        t_utc->y = t_loc->y;
00231        t_utc->m = t_loc->m;
00232        t_utc->d = t_loc->d;
00233        t_utc->h = t_utc->i = t_utc->s = 0;
00234        timelib_update_ts(t_utc, NULL);
00235 
00236        /* Compute d of 12h local mean solar time */
00237        timestamp = t_loc->sse;
00238        d = timelib_ts_to_juliandate(timestamp) - lon/360.0;
00239 
00240        /* Compute local sidereal time of this moment */
00241        sidtime = astro_revolution(astro_GMST0(d) + 180.0 + lon);
00242 
00243        /* Compute Sun's RA + Decl at this moment */
00244        astro_sun_RA_dec( d, &sRA, &sdec, &sr );
00245 
00246        /* Compute time when Sun is at south - in hours UT */
00247        tsouth = 12.0 - astro_rev180(sidtime - sRA) / 15.0;
00248 
00249        /* Compute the Sun's apparent radius, degrees */
00250        sradius = 0.2666 / sr;
00251 
00252        /* Do correction to upper limb, if necessary */
00253        if (upper_limb) {
00254               altit -= sradius;
00255        }
00256 
00257        /* Compute the diurnal arc that the Sun traverses to reach */
00258        /* the specified altitude altit: */
00259        {
00260               double cost;
00261               cost = (sind(altit) - sind(lat) * sind(sdec)) / (cosd(lat) * cosd(sdec));
00262               *ts_transit = t_utc->sse + (tsouth * 3600);
00263               if (cost >= 1.0) {
00264                      rc = -1;
00265                      t = 0.0;       /* Sun always below altit */
00266 
00267                      *ts_rise = *ts_set = t_utc->sse + (tsouth * 3600);
00268               } else if (cost <= -1.0) {
00269                      rc = +1;
00270                      t = 12.0;      /* Sun always above altit */
00271 
00272                      *ts_rise = t_loc->sse - (12 * 3600);
00273                      *ts_set  = t_loc->sse + (12 * 3600);
00274               } else {
00275                      t = acosd(cost) / 15.0;   /* The diurnal arc, hours */
00276 
00277                      /* Store rise and set times - as Unix Timestamp */
00278                      *ts_rise = ((tsouth - t) * 3600) + t_utc->sse;
00279                      *ts_set  = ((tsouth + t) * 3600) + t_utc->sse;
00280 
00281                      *h_rise = (tsouth - t);
00282                      *h_set  = (tsouth + t);
00283               }
00284        }
00285 
00286        /* Kill temporary time and restore original sse */
00287        timelib_time_dtor(t_utc);
00288        t_loc->sse = old_sse;
00289 
00290        return rc;
00291 }
00292 
00293 double timelib_ts_to_juliandate(timelib_sll ts)
00294 {
00295        double tmp;
00296 
00297        tmp = ts;
00298        tmp /= 86400;
00299        tmp += 2440587.5;
00300        tmp -= 2451543;
00301 
00302        return tmp;
00303 }