It's been a couple months since I heard anything, so I wanted to check the
status of my request and update you on my efforts.  Things got pretty busy
so I didn't have much time to work on this in August, but the last couple
weeks I've been able to really dive into it.

The latest patch is attached (based off 5.3.3 r303757 now).  As you can see,
it's no longer just a hacky proof-of-concept like the previous patch sent
out.  This actually makes full use of the formulas and has been thoroughly
tested.

The only issue I'm running into is precision.  You'll notice I've got a ton
of type casting all over the place (including a lot of unnecessary long
doubles lol) which I'll have cleaned-up in the next patch.  These formulas
are accurate to within a few seconds (over the next several thousand years
at least), but the level of precision required is a bit much for standard
floating points to handle.  As a result, the equinox/solstice start dates
are off by about 6-14 minutes when compiled on stock Ubuntu 10.04 and 0-7
minutes when compiled on Windows (VC9).  Is the PHP source currently making
use of any arbitrary precision libraries (like MPFR, etc) anywhere?  That
would easily solve this rounding issue and enable me to put together some
PHPT tests with predictable results.


Also, I would like to ping for an update on my SVN account request.  If
there is anything further you need me to do for you in order to complete
this process, I will be more than happy to do it.


Thanks!

--Kris


On Tue, Aug 10, 2010 at 1:36 PM, Kris Craig <kris.cr...@gmail.com> wrote:

> Woops, sorry.  Here's the file renamed to .txt.  Thanks for the tip!
>
>
> --Kris
>
>
>
> On Tue, Aug 10, 2010 at 12:50 PM, Michael Maclean <m...@php.net> wrote:
>
>> On 10/08/10 20:28, Kris Craig wrote:
>>
>>> Sorry, I guess it would help if I actually attached the patch.....  Here
>>> it is.
>>>
>>
>> The list strips attachments with filenames ending in something other than
>> .txt - resend or perhaps put it online somewhere?
>>
>> --
>> Cheers,
>> Michael
>>
>
>
--- php_date.c_diffbase.c       Fri Sep 24 18:50:09 2010
+++ php_date.c  Fri Sep 24 18:17:06 2010
@@ -967,6 +967,10 @@
        "Sun", "Mon", "Tue", "Wed", "Thu", "Fri", "Sat"
 };
 
+static char *season_names[] = {
+       "Winter", "Spring", "Summer", "Fall"
+};
+
 static char *english_suffix(timelib_sll number)
 {
        if (number >= 10 && number <= 19) {
@@ -982,6 +986,162 @@
 }
 /* }}} */
 
+/* {{{ Astronomical season equinox/solstice calculations.  --Kris */
+long double php_date_get_equinox_jde( int y, int season )
+{
+       long double jde, y2;
+       
+       //Equations derived from "Astronomical Algorithms" by Jean Meeus.  
--Kris
+       
+       y2 = (long double) ((int) y - 2000) / 1000;
+       
+       switch ( (int) season )
+       {
+               default:
+               case 0:
+                       //December solstice, start of winter.  --Kris
+                       jde = 2451900.05952 + (365242.74049 * y2) - (0.06223 * 
pow( y2, 2 )) - (0.00823 * pow( y2, 3 )) + (0.00032 * pow( y2, 4 ));
+                       break;
+               case 1:
+                       //March equinox, start of spring.  --Kris
+                       jde = 2451623.80984 + (365242.37404 * y2) + (0.05169 * 
pow( y2, 2 )) - (0.00411 * pow( y2, 3 )) - (0.00057 * pow( y2, 4 ));
+                       break;
+               case 2:
+                       //June solstice, start of summer.  --Kris
+                       jde = 2451716.56767 + (365241.62603 * y2) + (0.00325 * 
pow( y2, 2 )) + (0.00888 * pow( y2, 3 )) - (0.00030 * pow( y2, 4 ));
+                       break;
+               case 3:
+                       //September equinox, start of fall.  --Kris
+                       jde = 2451810.21715 + (365242.01767 * y2) - (0.11575 * 
pow( y2, 2 )) + (0.00337 * pow( y2, 3 )) + (0.00078 * pow( y2, 4 ));
+                       break;
+       }
+       
+       return jde;
+}
+
+int gregorian_base_from_jd( double jd )
+{
+       double Z, A, a;
+       
+       jd += 0.5;
+       
+       Z = (int) jd;
+       
+       if ( Z < 2299161 )
+       {
+               A = Z;
+       }
+       else
+       {
+               a = (int) ((Z - 1867216.25) / 36524.25);
+               A = Z + 1 + a - (int) (a / 4);
+       }
+       
+       return (A + 1524);
+}
+
+double dayofmonth_from_jd( double jd )
+{
+       double B, C, D, E;  //Ambiguity aside, I'm just keeping the variable 
names consistent with the original formulas.  --Kris
+       double F;
+       
+       B = (int) gregorian_base_from_jd( jd );
+       C = (int) ((B - 122.1) / 365.25);
+       D = (int) (365.25 * C);
+       E = (int) ((B - D) / 30.6001);
+       F = (double) (jd + 0.05) - (int) (jd + 0.05);
+       
+       return (B - D - (int) (30.6001 * E) + F);
+}
+
+int month_from_jd( double jd )
+{
+       int B, C, D, E;  //Ambiguity aside, I'm just keeping the variable names 
consistent with the original formulas.  --Kris
+       double F;
+       
+       B = (int) gregorian_base_from_jd( jd );
+       C = (int) ((B - 122.1) / 365.25);
+       D = (int) (365.25 * C);
+       E = (int) ((B - D) / 30.6001);
+       F = (double) (jd + 0.05) - (int) (jd + 0.05);
+       
+       if ( E < 14 )
+       {
+               return (E - 1);
+       }
+       else
+       {
+               return (E - 13);
+       }
+}
+
+int year_from_jd( double jd )
+{
+       double B, C;
+       
+       B = (int) gregorian_base_from_jd( jd );
+       C = (int) ((B - 122.1) / 365.25);
+       
+       if ( month_from_jd( jd ) > 2 )
+       {
+               return (C - 4716);
+       }
+       else
+       {
+               return (C - 4715);
+       }
+}
+
+double day_exact( double d, double h, double i, double s )
+{
+       return (double) (d + (h / 24) + (i / (24 * 60)) + (s / (24 * pow( 60, 2 
))));
+}
+
+double jd_from_date( double y, double m, double d )
+{
+       /* Our working formula variables.  --Kris */
+       double A, B;
+       
+       /* If Jan or Feb, it's considered to be the 13th or 14th month of the 
preceding year.  --Kris */
+       if ( m <= 2 )
+       {
+               m += 12;
+               y--;
+       }
+       
+       A = (int) (y / 100);
+       B = 2 - A + (int) (A / 4);
+       
+       return (double) ((int) (365.25 * (y + 4716)) + (int) (30.6001 * (m + 
1)) + d + B - 1524.5);
+}
+
+/* }}} */
+
+/* {{{ Get the current season.  If hemisphere cannot be determined, assume 
north.  --Kris */
+int php_date_get_season( double y, double m, double d, double h, double i, 
double s, char *timezone )
+{
+       double hemisphere, season, sloop;
+       double jd;
+       
+       //TODO - return (Season# + Hemisphere#) % 4;  //Where Hemisphere# is 0 
for North, 2 for South.  --Kris
+       hemisphere = 0;
+       
+       jd = jd_from_date( y, m, day_exact( (double) d, (double) h, (double) i, 
(double) s ) );
+       
+       season = 0;
+       for ( sloop = 4; sloop > 0; sloop-- )
+       {
+               if ( jd >= php_date_get_equinox_jde( y, ((int) sloop % 4) ) )
+               {
+                       season = ((int) sloop % 4);
+                       break;
+               }
+       }
+       
+       return (int) (season + hemisphere) % 4;
+}
+/* }}} */
+
 /* {{{ day of week helpers */
 char *php_date_full_day_name(timelib_sll y, timelib_sll m, timelib_sll d)
 {
@@ -1067,6 +1227,12 @@
                        case 'L': length = slprintf(buffer, 32, "%d", 
timelib_is_leap((int) t->y)); break;
                        case 'y': length = slprintf(buffer, 32, "%02d", (int) 
t->y % 100); break;
                        case 'Y': length = slprintf(buffer, 32, "%s%04lld", 
t->y < 0 ? "-" : "", php_date_llabs((timelib_sll) t->y)); break;
+                       
+                       /* season */
+                       case 'v': length = slprintf(buffer, 32, "%s", 
season_names[(int) php_date_get_season( (double) t->y, (double) t->m, (double) 
t->d, 
+                               (double) t->h, (double) t->i, (double) t->s, 
+                               localtime && t->zone_type == 
TIMELIB_ZONETYPE_ID ? t->tz_info->name : "" )]); break;
 
                        /* time */
                        case 'a': length = slprintf(buffer, 32, "%s", t->h >= 
12 ? "pm" : "am"); break;
-- 
PHP Internals - PHP Runtime Development Mailing List
To unsubscribe, visit: http://www.php.net/unsub.php

Reply via email to