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

Projects Programming
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.