PHP  
 PHP: Test and Code Coverage Analysis
downloads | QA | documentation | faq | getting help | mailing lists | reporting bugs | php.net sites | links | my php.net 
 

LCOV - code coverage report
Current view: top level - ext/date/lib - astro.c (source / functions) Hit Total Coverage
Test: PHP Code Coverage Lines: 75 75 100.0 %
Date: 2014-12-15 Functions: 7 7 100.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /*
       2             :    +----------------------------------------------------------------------+
       3             :    | PHP Version 5                                                        |
       4             :    +----------------------------------------------------------------------+
       5             :    | Copyright (c) 1997-2014 The PHP Group                                |
       6             :    +----------------------------------------------------------------------+
       7             :    | This source file is subject to version 3.01 of the PHP license,      |
       8             :    | that is bundled with this package in the file LICENSE, and is        |
       9             :    | available through the world-wide-web at the following url:           |
      10             :    | http://www.php.net/license/3_01.txt                                  |
      11             :    | If you did not receive a copy of the PHP license and are unable to   |
      12             :    | obtain it through the world-wide-web, please send a note to          |
      13             :    | license@php.net so we can mail you a copy immediately.               |
      14             :    +----------------------------------------------------------------------+
      15             :    | Algorithms are taken from a public domain source by Paul             |
      16             :    | Schlyter, who wrote this in December 1992                            |
      17             :    +----------------------------------------------------------------------+
      18             :    | Authors: Derick Rethans <derick@derickrethans.nl>                    |
      19             :    +----------------------------------------------------------------------+
      20             :  */
      21             : 
      22             : /* $Id$ */
      23             : 
      24             : #include <stdio.h>
      25             : #include <math.h>
      26             : #include "timelib.h"
      27             : 
      28             : #define days_since_2000_Jan_0(y,m,d) \
      29             :         (367L*(y)-((7*((y)+(((m)+9)/12)))/4)+((275*(m))/9)+(d)-730530L)
      30             : 
      31             : #ifndef PI
      32             :  #define PI        3.1415926535897932384
      33             : #endif
      34             : 
      35             : #define RADEG     ( 180.0 / PI )
      36             : #define DEGRAD    ( PI / 180.0 )
      37             : 
      38             : /* The trigonometric functions in degrees */
      39             : 
      40             : #define sind(x)  sin((x)*DEGRAD)
      41             : #define cosd(x)  cos((x)*DEGRAD)
      42             : #define tand(x)  tan((x)*DEGRAD)
      43             : 
      44             : #define atand(x)    (RADEG*atan(x))
      45             : #define asind(x)    (RADEG*asin(x))
      46             : #define acosd(x)    (RADEG*acos(x))
      47             : #define atan2d(y,x) (RADEG*atan2(y,x))
      48             : 
      49             : 
      50             : /* Following are some macros around the "workhorse" function __daylen__ */
      51             : /* They mainly fill in the desired values for the reference altitude    */
      52             : /* below the horizon, and also selects whether this altitude should     */
      53             : /* refer to the Sun's center or its upper limb.                         */
      54             : 
      55             : 
      56             : #include "astro.h"
      57             : 
      58             : /******************************************************************/
      59             : /* This function reduces any angle to within the first revolution */
      60             : /* by subtracting or adding even multiples of 360.0 until the     */
      61             : /* result is >= 0.0 and < 360.0                                   */
      62             : /******************************************************************/
      63             : 
      64             : #define INV360    (1.0 / 360.0)
      65             : 
      66             : /*****************************************/
      67             : /* Reduce angle to within 0..360 degrees */
      68             : /*****************************************/
      69        2334 : static double astro_revolution(double x)
      70             : {
      71        2334 :         return (x - 360.0 * floor(x * INV360));
      72             : }
      73             : 
      74             : /*********************************************/
      75             : /* Reduce angle to within +180..+180 degrees */
      76             : /*********************************************/
      77         778 : static double astro_rev180( double x )
      78             : {
      79         778 :         return (x - 360.0 * floor(x * INV360 + 0.5));
      80             : }
      81             : 
      82             : /*******************************************************************/
      83             : /* This function computes GMST0, the Greenwich Mean Sidereal Time  */
      84             : /* at 0h UT (i.e. the sidereal time at the Greenwhich meridian at  */
      85             : /* 0h UT).  GMST is then the sidereal time at Greenwich at any     */
      86             : /* time of the day.  I've generalized GMST0 as well, and define it */
      87             : /* as:  GMST0 = GMST - UT  --  this allows GMST0 to be computed at */
      88             : /* other times than 0h UT as well.  While this sounds somewhat     */
      89             : /* contradictory, it is very practical:  instead of computing      */
      90             : /* GMST like:                                                      */
      91             : /*                                                                 */
      92             : /*  GMST = (GMST0) + UT * (366.2422/365.2422)                      */
      93             : /*                                                                 */
      94             : /* where (GMST0) is the GMST last time UT was 0 hours, one simply  */
      95             : /* computes:                                                       */
      96             : /*                                                                 */
      97             : /*  GMST = GMST0 + UT                                              */
      98             : /*                                                                 */
      99             : /* where GMST0 is the GMST "at 0h UT" but at the current moment!   */
     100             : /* Defined in this way, GMST0 will increase with about 4 min a     */
     101             : /* day.  It also happens that GMST0 (in degrees, 1 hr = 15 degr)   */
     102             : /* is equal to the Sun's mean longitude plus/minus 180 degrees!    */
     103             : /* (if we neglect aberration, which amounts to 20 seconds of arc   */
     104             : /* or 1.33 seconds of time)                                        */
     105             : /*                                                                 */
     106             : /*******************************************************************/
     107             : 
     108         778 : static double astro_GMST0(double d)
     109             : {
     110             :         double sidtim0;
     111             :         /* Sidtime at 0h UT = L (Sun's mean longitude) + 180.0 degr  */
     112             :         /* L = M + w, as defined in sunpos().  Since I'm too lazy to */
     113             :         /* add these numbers, I'll let the C compiler do it for me.  */
     114             :         /* Any decent C compiler will add the constants at compile   */
     115             :         /* time, imposing no runtime or code overhead.               */
     116         778 :         sidtim0 = astro_revolution((180.0 + 356.0470 + 282.9404) + (0.9856002585 + 4.70935E-5) * d);
     117         778 :         return sidtim0;
     118             : } 
     119             : 
     120             : /* This function computes the Sun's position at any instant */
     121             : 
     122             : /******************************************************/
     123             : /* Computes the Sun's ecliptic longitude and distance */
     124             : /* at an instant given in d, number of days since     */
     125             : /* 2000 Jan 0.0.  The Sun's ecliptic latitude is not  */
     126             : /* computed, since it's always very near 0.           */
     127             : /******************************************************/
     128         778 : static void astro_sunpos(double d, double *lon, double *r)
     129             : {
     130             :         double M,         /* Mean anomaly of the Sun */
     131             :                w,         /* Mean longitude of perihelion */
     132             :                           /* Note: Sun's mean longitude = M + w */
     133             :                e,         /* Eccentricity of Earth's orbit */
     134             :                E,         /* Eccentric anomaly */
     135             :                x, y,      /* x, y coordinates in orbit */
     136             :                v;         /* True anomaly */
     137             : 
     138             :         /* Compute mean elements */
     139         778 :         M = astro_revolution(356.0470 + 0.9856002585 * d);
     140         778 :         w = 282.9404 + 4.70935E-5 * d;
     141         778 :         e = 0.016709 - 1.151E-9 * d;
     142             : 
     143             :         /* Compute true longitude and radius vector */
     144         778 :         E = M + e * RADEG * sind(M) * (1.0 + e * cosd(M));
     145         778 :         x = cosd(E) - e;
     146         778 :         y = sqrt(1.0 - e*e) * sind(E);
     147         778 :         *r = sqrt(x*x + y*y);              /* Solar distance */
     148         778 :         v = atan2d(y, x);                  /* True anomaly */
     149         778 :         *lon = v + w;                        /* True solar longitude */
     150         778 :         if (*lon >= 360.0) {
     151          28 :                 *lon -= 360.0;                   /* Make it 0..360 degrees */
     152             :         }
     153         778 : }
     154             : 
     155         778 : static void astro_sun_RA_dec(double d, double *RA, double *dec, double *r)
     156             : {
     157             :         double lon, obl_ecl, x, y, z;
     158             : 
     159             :         /* Compute Sun's ecliptical coordinates */
     160         778 :         astro_sunpos(d, &lon, r);
     161             : 
     162             :         /* Compute ecliptic rectangular coordinates (z=0) */
     163         778 :         x = *r * cosd(lon);
     164         778 :         y = *r * sind(lon);
     165             : 
     166             :         /* Compute obliquity of ecliptic (inclination of Earth's axis) */
     167         778 :         obl_ecl = 23.4393 - 3.563E-7 * d;
     168             : 
     169             :         /* Convert to equatorial rectangular coordinates - x is unchanged */
     170         778 :         z = y * sind(obl_ecl);
     171         778 :         y = y * cosd(obl_ecl);
     172             : 
     173             :         /* Convert to spherical coordinates */
     174         778 :         *RA = atan2d(y, x);
     175         778 :         *dec = atan2d(z, sqrt(x*x + y*y));
     176         778 : }
     177             : 
     178             : /**
     179             :  * Note: timestamp = unixtimestamp (NEEDS to be 00:00:00 UT)
     180             :  *       Eastern longitude positive, Western longitude negative       
     181             :  *       Northern latitude positive, Southern latitude negative       
     182             :  *       The longitude value IS critical in this function!            
     183             :  *       altit = the altitude which the Sun should cross              
     184             :  *               Set to -35/60 degrees for rise/set, -6 degrees       
     185             :  *               for civil, -12 degrees for nautical and -18          
     186             :  *               degrees for astronomical twilight.                   
     187             :  *         upper_limb: non-zero -> upper limb, zero -> center         
     188             :  *               Set to non-zero (e.g. 1) when computing rise/set     
     189             :  *               times, and to zero when computing start/end of       
     190             :  *               twilight.                                            
     191             :  *        *rise = where to store the rise time                        
     192             :  *        *set  = where to store the set  time                        
     193             :  *                Both times are relative to the specified altitude,  
     194             :  *                and thus this function can be used to compute       
     195             :  *                various twilight times, as well as rise/set times   
     196             :  * Return value:  0 = sun rises/sets this day, times stored at        
     197             :  *                    *trise and *tset.                               
     198             :  *               +1 = sun above the specified "horizon" 24 hours.     
     199             :  *                    *trise set to time when the sun is at south,    
     200             :  *                    minus 12 hours while *tset is set to the south  
     201             :  *                    time plus 12 hours. "Day" length = 24 hours     
     202             :  *               -1 = sun is below the specified "horizon" 24 hours   
     203             :  *                    "Day" length = 0 hours, *trise and *tset are    
     204             :  *                    both set to the time when the sun is at south.  
     205             :  *                                                                    
     206             :  */
     207         778 : 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)
     208             : {
     209             :         double  d,  /* Days since 2000 Jan 0.0 (negative before) */
     210             :         sr,         /* Solar distance, astronomical units */
     211             :         sRA,        /* Sun's Right Ascension */
     212             :         sdec,       /* Sun's declination */
     213             :         sradius,    /* Sun's apparent radius */
     214             :         t,          /* Diurnal arc */
     215             :         tsouth,     /* Time when Sun is at south */
     216             :         sidtime;    /* Local sidereal time */
     217             :         timelib_time *t_utc;
     218             :         timelib_sll   timestamp, old_sse;
     219             : 
     220         778 :         int rc = 0; /* Return cde from function - usually 0 */
     221             : 
     222             :         /* Normalize time */
     223         778 :         old_sse = t_loc->sse;
     224         778 :         t_loc->h = 12;
     225         778 :         t_loc->i = t_loc->s = 0;
     226         778 :         timelib_update_ts(t_loc, NULL);
     227             : 
     228             :         /* Calculate TS belonging to UTC 00:00 of the current day */
     229         778 :         t_utc = timelib_time_ctor();
     230         778 :         t_utc->y = t_loc->y;
     231         778 :         t_utc->m = t_loc->m;
     232         778 :         t_utc->d = t_loc->d;
     233         778 :         t_utc->h = t_utc->i = t_utc->s = 0;
     234         778 :         timelib_update_ts(t_utc, NULL);
     235             : 
     236             :         /* Compute d of 12h local mean solar time */
     237         778 :         timestamp = t_loc->sse;
     238         778 :         d = timelib_ts_to_juliandate(timestamp) - lon/360.0;
     239             : 
     240             :         /* Compute local sidereal time of this moment */
     241         778 :         sidtime = astro_revolution(astro_GMST0(d) + 180.0 + lon);
     242             : 
     243             :         /* Compute Sun's RA + Decl at this moment */
     244         778 :         astro_sun_RA_dec( d, &sRA, &sdec, &sr );
     245             : 
     246             :         /* Compute time when Sun is at south - in hours UT */
     247         778 :         tsouth = 12.0 - astro_rev180(sidtime - sRA) / 15.0;
     248             : 
     249             :         /* Compute the Sun's apparent radius, degrees */
     250         778 :         sradius = 0.2666 / sr;
     251             : 
     252             :         /* Do correction to upper limb, if necessary */
     253         778 :         if (upper_limb) {
     254         619 :                 altit -= sradius;
     255             :         }
     256             : 
     257             :         /* Compute the diurnal arc that the Sun traverses to reach */
     258             :         /* the specified altitude altit: */
     259             :         {
     260             :                 double cost;
     261         778 :                 cost = (sind(altit) - sind(lat) * sind(sdec)) / (cosd(lat) * cosd(sdec));
     262         778 :                 *ts_transit = t_utc->sse + (tsouth * 3600);
     263         778 :                 if (cost >= 1.0) {
     264          90 :                         rc = -1;
     265          90 :                         t = 0.0;       /* Sun always below altit */
     266             : 
     267          90 :                         *ts_rise = *ts_set = t_utc->sse + (tsouth * 3600);
     268         688 :                 } else if (cost <= -1.0) {
     269          32 :                         rc = +1;
     270          32 :                         t = 12.0;      /* Sun always above altit */
     271             : 
     272          32 :                         *ts_rise = t_loc->sse - (12 * 3600);
     273          32 :                         *ts_set  = t_loc->sse + (12 * 3600);
     274             :                 } else {
     275         656 :                         t = acosd(cost) / 15.0;   /* The diurnal arc, hours */
     276             : 
     277             :                         /* Store rise and set times - as Unix Timestamp */
     278         656 :                         *ts_rise = ((tsouth - t) * 3600) + t_utc->sse;
     279         656 :                         *ts_set  = ((tsouth + t) * 3600) + t_utc->sse;
     280             : 
     281         656 :                         *h_rise = (tsouth - t);
     282         656 :                         *h_set  = (tsouth + t);
     283             :                 }
     284             :         }
     285             : 
     286             :         /* Kill temporary time and restore original sse */
     287         778 :         timelib_time_dtor(t_utc);
     288         778 :         t_loc->sse = old_sse;
     289             : 
     290         778 :         return rc;
     291             : }
     292             : 
     293         778 : double timelib_ts_to_juliandate(timelib_sll ts)
     294             : {
     295             :         double tmp;
     296             : 
     297         778 :         tmp = ts;
     298         778 :         tmp /= 86400;
     299         778 :         tmp += 2440587.5;
     300         778 :         tmp -= 2451543;
     301             : 
     302         778 :         return tmp;
     303             : }

Generated by: LCOV version 1.10

Generated at Mon, 15 Dec 2014 17:02:40 +0000 (2 days ago)

Copyright © 2005-2014 The PHP Group
All rights reserved.