Date formula hurts my head, have you seen it before and how does it work?

Projects Programming
1

Here is the formula;

d = (int)(30.6*((month + 9) % 12) + 58.5 + date) % 365;

All the variables are int.
Month is the month of the year 1...12
Date is the day of the month 1...31

The formula sets 'd' to be the day of the year 0...364
where 0 is the 1/Jan and 364 is the 31/Dec.
The formula does not include the year and so it cannot compensate for leap years.
For 29/Feb the formula returns 59 which is the same as 1/Mar so in leap years day 59 occurs twice.

I was unaware of this formula until recently.
My questions are;

  • Is it well known
  • Can anybody explain simply how it works

The formula does work. Below is some output which I hope is self explanatory.
I printed the returned values as if each month had 31 days.

      31   28   31   30   31   30   31   31   30   31   30   31
     Jan  Feb  Mar  Apr  May  Jun  Jul  Aug  Sep  Oct  Nov  Dec
01     0   31   59   90  120  151  181  212  243  273  304  334
02     1   32   60   91  121  152  182  213  244  274  305  335
03     2   33   61   92  122  153  183  214  245  275  306  336
04     3   34   62   93  123  154  184  215  246  276  307  337
05     4   35   63   94  124  155  185  216  247  277  308  338
06     5   36   64   95  125  156  186  217  248  278  309  339
07     6   37   65   96  126  157  187  218  249  279  310  340
08     7   38   66   97  127  158  188  219  250  280  311  341
09     8   39   67   98  128  159  189  220  251  281  312  342
10     9   40   68   99  129  160  190  221  252  282  313  343
11    10   41   69  100  130  161  191  222  253  283  314  344
12    11   42   70  101  131  162  192  223  254  284  315  345
13    12   43   71  102  132  163  193  224  255  285  316  346
14    13   44   72  103  133  164  194  225  256  286  317  347
15    14   45   73  104  134  165  195  226  257  287  318  348
16    15   46   74  105  135  166  196  227  258  288  319  349
17    16   47   75  106  136  167  197  228  259  289  320  350
18    17   48   76  107  137  168  198  229  260  290  321  351
19    18   49   77  108  138  169  199  230  261  291  322  352
20    19   50   78  109  139  170  200  231  262  292  323  353
21    20   51   79  110  140  171  201  232  263  293  324  354
22    21   52   80  111  141  172  202  233  264  294  325  355
23    22   53   81  112  142  173  203  234  265  295  326  356
24    23   54   82  113  143  174  204  235  266  296  327  357
25    24   55   83  114  144  175  205  236  267  297  328  358
26    25   56   84  115  145  176  206  237  268  298  329  359
27    26   57   85  116  146  177  207  238  269  299  330  360
28    27   58   86  117  147  178  208  239  270  300  331  361
29    28   59   87  118  148  179  209  240  271  301  332  362
30    29   60   88  119  149  180  210  241  272  302  333  363
31    30   61   89  120  150  181  211  242  273  303  334  364
2

Hey guys I can see people reading the post, but nobody is responding.
If you have or have not seen this formula before then post even if you don't know how it works.
I am really surprised that I had not seen it before.

3

Not sure I've seen that exact formula. I use the function below, which does allow for leap years. Looks like it has at least some similarities.

//Return ordinal day of year for the given time_t
int ordinalDate(time_t t)
{
    int m = month(t);
    int d = day(t);
    
    if (m == 1)
        return d;
    else if (m == 2)
        return d + 31;
    else {
        int n = floor(30.6 * (m + 1)) + d - 122;
        return n + (isLeap(t) ? 60 : 59);
    }
}

//Leap years are those divisible by 4, but not those divisible by 100,
//except that those divisble by 400 *are* leap years.
//See Kernighan & Ritchie, 2nd edition, section 2.5.
boolean isLeap(time_t t)
{
    int y = year(t);
    return (y % 4 == 0 && y % 100 != 0) || y % 400 == 0;
}
4

Hello Jack,

I can see the similarities, but you have several lines of code and a lot of magic numbers.
The original formula has the magic numbers but it is just one line.

Where does the 30.6 come from?

5

Dunno. Not my formula. 365/12 is 30.4167 and 366/12 is 30.5.

6

nobody is responding.

Probably because we don't understand it either :slight_smile:

I suspect the 30.6 is the average days per month bumped by .1 to force a rounding upwards or some such??

Rob

7

I can't quite get my head around it but I partly follow it.

It's a way of fudging the number of days in a month.

30.6 is the average number of days in a month excluding Jan and Feb.
((month + 9) % 12) is doing two things:

It converts the month number from 1-12 to 0-11.
It adds ten to the month. Multiply that by the average-days-per-month for Mar-Dec, add in 58.5 representing the number of days in Jan+Feb, and you have added a year's worth of days to the date. Since the final answer is mod 365, that doesn't affect the answer.

I can kinda vaguely see how that might be producing the right answer, but it still feels like a most extraordinary cheat to me.

8

This formula seems to be related to or is the same as Zeller's Congruence:

Maybe if you remove the "mod 7" modulo operation in Zeller's formula you obtain a similar result as yours.

Zeller's Congruence does take into account leap years.

:wink:

9

OK I am slightly less brain-dead than when I replied before, here is an explanation of the code I posted above:

10

I can kinda vaguely see how that might be producing the right answer, but it still feels like a most extraordinary cheat to me.

Exactly! When I first saw it I thought that cannot possibly work, but checking it it does. My next thought was "Why have I never seen this before, is it just me, has everybody else been using this equation for years?".

Zeller's Congruence and the Ordinal Date code are making my head hurt even more, but I'll try and work through them -thanks.

I would still be interested if people have or have not seen the formula before. Is it a little gem that has been lying lost in cyberspace or is it really quite well know?

11

I wouldn't call this an explanation, but it's something like one.

Starting with March, the number of days in each month follows this pattern: 31, 30, 31, 30, 31 / 31, 30, 31, 30, 31 / 31, (28 or 29). You could look at it as three sets of (31, 30, 31, 30, 31), with the last set truncated at the end of February. Advantages of looking at the calendar this way, in terms of calculating which day of the year corresponds to a particular month and day, are:- The number of days in each month follows a simple pattern, and- The oddity, February, is at the end of the list, where it requires no special treatment beyond "% 365."

This code(month + 9) % 12yields a number between 0 and 11, inclusive, with zero for March, 9 for December, 10 for January, and 11 for February, so it's function is to realign the months starting with March, to take advantage of 11 months of regularity. 30.6 is the average number of days in each group of five months, so multiplying by that number yields an estimate of the "ordinal date" of the first day of each month. That figure isn't precise, though, and it needs some adjustment. Looking at adding 58.5 as two separate operations - adding 0.5 and adding 58 - then adding 0.5 adjusts the number so that taking the int() function yields the proper ordinal date for the first day of the month. Adding 58 - the ordinal date of February 28 - adjusts for the fact that we temporarily treated March 1 as the first day of the year. Taking modulo 365 accounts for the fact that the calculation yields ordinal dates for January and February that are bigger than 365.

Voila!

You'll note that the calculation doesn't do anything for leap year - the year doesn't figure in to the formula, and it gives the same result for Feb 29 as for March 1. To be general, it needs to know the year and how to tell if it's a leap year, and make an adjustment for ordinal dates bigger than 58.

There's some hand-waving in that description above - particularly with regard to adjusting an estimate of the ordinal date by adding 0.5. I don't know why that works, but I suspect that it's largely because the size of the pattern - five elements - is small, and because the elements only differ by 1. But I'm not sure, and I don't want to do the work to figure it out. Calendar calculations are can be ferociously non-intuitive - otherwise, the Gregorian calendar would have been developed in ancient Mesopotamia, and we'd call it the Babylonian calendar. That, and they've usually been worked out long ago by someone who cared a lot more about them than I do. It's worth noting that the formula works with constants in a small range around 30.6, and around 58.5.

If you're hoping that you'll be able to see how this formula was derived, you may want to give that up. My guess is that this formula was developed at least in part by laborious trial-and-error. Unless you're eager to learn how to manipulate the calendar for your own edification, you'll be well-advised to be satisfied with a general understanding of why this calculation works, and move on to something you really like.

12

move on to something you really like.

:slight_smile:

13

Where did you find this? It looks like something I wrote a long time ago for a VB program at an old employer (except I don't think I was calculating a julian date - it's been a long time). I just remember thinking at the time that it was pretty neat that I was able to cram into a short line of code something that would have been spread out more; one reason I did it that way was because the code ran faster (I think it was part of some reporting on a bunch of rows from a DB or something). I am pretty sure I commented the code, though, so that others would understand what was happening (even so, I remember thinking that someone in the future would either think I was mad, insane, a genious, or some combination of all three - and that they'd be cursing my name). I still get a smile from it, even now (evil little part of me, I guess). :smiley:

14

Where did you find this?

A person called David Williams wrote a program to control a heliostat in 1987.
He produced a QBasic program in 2007 based on his original program.
Somebody called Marius Gundersen converted this to C in 2009.
The formula was uncommented and I am tending to tmd3's thought that it was developed by trial and error.

Gundersen's code gave a clean compile immediately. However I think I may not have the the final working version as I have spotted a couple problem areas.

The QBasic formula was;

D = INT(30.6 * ((Mth% + 9) MOD 12) + 58.5 + Day%) MOD 365 ' day of year (D = 0 on Jan 1)

Gundersen converted this to;

uint8_t d = 30.6*((month + 9) % 12) + 58.5 + date % 365;

I think the conversion is a bit suspect because the month and date are uint8_t and he has got rid of the explicit casting to an int and setting the precedence of the %365. So it might work on some compilers but not all.
That is why I put the cast back in and used int throughout;

d = (int)(30.6*((month + 9) % 12) + 58.5 + date) % 365;

What makes me really suspect the code is this function where *az is set but az is tested;

void c2p(double x, double y, double z, double *az, double *el)
{
	//Cartesian to Polar. Convert from X,Y,Z to AZ,EL 

	/*
		Ang is used twice, first to calculate the angle of elevation, and then the azimuth.
		Ang can produce a negative output, down to -pi/2, and this is acceptable as an angle of elevation. 
		The direction in which the mirror must be aimed, for example, can be below the horizontal. 
		But negative azimuths are not used. If Ang produces a negative value for the azimuth,
		the constant PY2 (2pi) is added to produce a positive equivalent. 
		AZ is therefore in the range 0 <= AZ < 2pi.
	*/
 
	*el = ang(sqrt(x*x + y*y), z);
	*az = ang(y, x);
	if(az<0) // !!!!! *az
	{
		*az+=pi2;
	}
	return;	  
}

The comments were not in the code I bolted them in from a description of the program given by David Williams. Maybe -Ve azimuths don't occur very often so the problem was not spotted?

So far that is all I can see that seems wrong. I would quite like to get the code working and have ordered a book containg tables of time and solar declination to try to check the output.