Sunset/Sunrise Calculator

/*****************************************************
*This is an implementation of the Sunrise/Sunset Algorithm
*found at 
* http://williams.best.vwh.net/sunrise_sunset_algorithm.htm
*from the Almanac for Computers, 1990
*Published by Nautical Almanac Office
*Washington, DC 20392
*Implemented by Chris Snyder
*****************************************************/
/*****************************************************
*Function name:calcSunset
*Params: $date(The day to calculate sunset or sunrise for)
*	$lat(The latitude for the location to calculate sunrise or sunset for)
*	$lon(The longitude for the location to calculate sunrise or sunset for west is negative)
*	$sunset(if set to 1 calculates sunset if set to 0 sunrise)
*	$GMToffset(difference in hours from GMT)
*****************************************************/
function calcSunset($date, $lat=41.57, $lon=-86.2, $sunset=1, $GMToffset=-4)
{
//The sun's Zenith
$zenith = 90.83;

//The sun's Zenith to calculate civil twilight
//$zenith=96;

//need the day of year
$day = getdate($date);
$N = $day['yday']+1;
        
//Convert the longitude to hour value and calculate and
//Approximate time.
	$lonHour = ($lon/15);
	if($sunset)
		$t=$N+((18-$lonHour)/24);
	else
		$t=$N+((6-$lonHour)/24);

//Calculate the Sun's mean anomaly
	$M = (0.9856*$t)-3.289;

//Calculate the Sun's true longitude
	$sinM = sin(deg2rad($M));
	$sin2M = sin(2*deg2rad($M));
	
	$L=$M +(1.916*$sinM) + (0.02*$sin2M) + 282.634;

	if($L>360)
		$L-=360;
	if($L<0)
		$L+=360;

//Calculate the Sun's right ascension(RA)
	$tanL = 0.91764 * tan(deg2rad($L));
	$RA = rad2deg(atan($tanL));
	if($RA>360)
		$RA-=360;
	if($RA<0)
		$RA+=360;

//Putting the RA value into the same quadrant as L
	$LQ = (floor($L/90))*90;
	$RAQ = (floor($RA/90))*90;
	$RA = $RA + ($LQ - $RAQ);

//Convert RA values to hours
	$RA /= 15;

//calculate the Sun's declination
	$sinDec = 0.39782 * sin(deg2rad($L));
	$cosDec = cos(asin($sinDec));

//calculate the Sun's local hour angle
//if cosH > 1 the sun never rises on this date at this location
//if cosH < -1 the sun never sets on this date at this location
	$cosH = (cos(deg2rad($zenith)) - ($sinDec * sin(deg2rad($lat)))) / ($cosDec * cos(deg2rad($lat)));

//finish calculating H and convert into hours
	if($sunset)
		$H = rad2deg(acos($cosH));
	else
		$H = 360 - rad2deg(acos($cosH));
	$H /= 15;

//Calculate local mean time of rising/setting
	$T = $H + $RA - (0.06571 * $t) - 6.622;

//Adjust back to UTC
	$UT = $T - $lonHour;
	if($UT>24)
		$UT -= 24;
	if($UT<0)
		$UT += 24;
//Adjust for current time zone
	$UT=$UT+$GMToffset;
	if($UT<0)
		$UT+=24;

//format to readable form
	$hour = floor($UT);
	$minutes = round(60*($UT - $hour));
	$minutes = str_pad($minutes, 3, ":0", STR_PAD_LEFT);

	$UT = $hour.$minutes;
	return $UT;
}