<?php
//=============================================================================
// Copyright (C) 2005 Marc A. Murison
//
// This program is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public License
// (see http://www.gnu.org/copyleft/gpl.html) as published by the
// Free Software Foundation; either version 2 of the License, or (at
// your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License (http://www.gnu.org/copyleft/gpl.html)
// for more details.
//=============================================================================
//=============================================================================
// Calculate sunrise, sunset, and twilights for the current day
// for a given latitude & longitude.
//
// Not accurate for high latitudes (> 60 deg) for certain times of year.
// Also, you can calculate rise/set/twilight times and print them accurate
// to less than a minute if you wish (default accuracies for internal
// calculations here are 1 second), but that would be pretty silly.
//
// Author: Marc A. Murison
// U.S. Naval Observatory
// Washington, DC 20392
// murisonATusno.navy.mil
// http://www.alpheratz.net/murison/
//
// You may use this code in your own code as long as you (or your
// company) agree that you (1) will not profit monetarily from it in
// any way, (2) will retain the Author attribution shown above, and
// (3) abide by the terms of the GNU General Public License.
//=============================================================================
require_once("math.php");
//-----------------------------------------------------------------------------
// main function to calculate twilight and rise/set times for the Sun
//
// lat and lon must be decimal degrees
//
// optional arguments (include either all three or none):
// integer day of month, month number (January=1), 4-digit year
//
// A typical call might look like this (defaults to today's date):
// echo printRiseSets( SunRiseSet( $lat_DC, $lon_DC ) );
// or this (for June 3, 2005):
// echo printRiseSets( SunRiseSet( $lat_DC, $lon_DC, 3, 6, 2005 ) );
//-----------------------------------------------------------------------------
function SunRiseSet( $lat /*decimal degrees North*/,
$lon /*decimal degrees West longitude*/ )
{
global $torad, $tzone, $ratio, $terror;
//constants
$ratio = 1.00273790935; //sidereal days per solar day
$torad = M_PI/180.0; //convert degrees to radians
$zset = 90.83327*$torad; //rise/set zenith angle
$zcivil = 96*$torad; //civil twilight zenith angle
$zastro = 108*$torad; //astronomical twilight zenith angle
$terror = 0.1/60.0; //tolerance in decimal hours
//if no date provided, assume today
$nargs = func_num_args();
$args = func_get_args();
if( $nargs > 2 ) {
if( $nargs != 5 ) {
print("<p>SunRiseSet error: wrong number of arguments!</p>");
return;
}
$day = $args[2];
$month = $args[3];
$year = $args[4];
} else {
$today = getdate(time()); //use PHP function to get current local date
$day = $today['mday']; //day of month
$month = $today['mon']; //month number
$year = $today['year'];
}
//day numbers
$jd0 = JD0hUT($lon,$tzone,$day,$month,$year); //JD at 0h UT
$ndays = $jd0 - 2451545.0; //number of days from J2000.0 = Jan. 1.5 UT 2000
//calculate Greenwich mean sideral time at 0h UT in decimal hours of time
$GMST0 = GMST0hUT($jd0);
//local parameters
// $refraction = 34.0/60.0; //average refraction at horizon is 34 arcmin
// $zset = (90 + $sd + $refraction)*$torad; //sunrise/set zenith angle
//Loop to calculate the rise/set times. For latitudes far from the Arctic
//Circle, just one iteration is usually sufficient to provide accuracies
//under a minute. The following loop always does at least two iterations.
$t_rise = 6; //initial guess for rise times
$t_set = 18; //initial guess for set times
$times[2] = $t_rise; //sunrise
$times[3] = $t_set; //sunset
$times[1] = $t_rise; //civil twilight -- rise
$times[4] = $t_set; //civil twilight -- set
$times[0] = $t_rise; //astronomical twilight -- rise
$times[5] = $t_set; //astronomical twilight -- set
$zvals = array($zastro,$zcivil,$zset,$zset,$zcivil,$zastro);
$signh = array(-1,-1,-1,1,1,1);
$countmax = 15; //max number of iterations
for( $i=0; $i < 6; ++$i ) {
$tnew = $times[$i]; $told = 0; $count = 0;
while( abs($tnew - $told) > $terror && $count < $countmax ) {
$told = $tnew;
$tnew = improvelocaltime( $told, $zvals[$i], $GMST0, $ndays,
$lat, $lon, $signh[$i], $azimuths[$i] );
++$count;
}
$times[$i] = $tnew;
$azimuths[$i] /= $torad;
}
//calculate the local meridian transit time of the Sun
$ttrans = Sun_transit_time( $lon, $ndays );
//calculate the altitude of the Sun at transit
SunEqCoords( $ttrans+$tzone, $ndays, $alpha /*decimal hours*/, $dec /*radians*/ );
$alt = 90 - $lat + $dec/$torad;
//output an array of values that can be used $vals[$i]["t"] for times
//and $vals[$i]["a"] for corresponding azimuths (in the case of rise
//and set times) or altitude (in the case of meridian transit), with
//the order being [astro twilight rise, civil twilight rise, sunrise,
//meridian transit, set, civil twilight set, astro twilight set] as
//$i runs from 0-6
for( $i=0; $i < 7; ++$i ) {
if( $i < 3 ) {
$vals[$i]["t"] = $times[$i];
$vals[$i]["a"] = $azimuths[$i];
} else if( $i == 3 ) {
$vals[$i]["t"] = $ttrans;
$vals[$i]["a"] = $alt;
} else {
$vals[$i]["t"] = $times[$i-1];
$vals[$i]["a"] = $azimuths[$i-1];
}
}
return $vals;
}
//-----------------------------------------------------------------------------
// given an initial guess for the local time t (decimal hours) corresponding
// to a sunrise or sunset, calculate the rise or set time (decimal hours)
// based on that initial guess
//
// the sun's azimuth is put into $az, in radians
//-----------------------------------------------------------------------------
function improvelocaltime( $t /*decimal hours*/, $z /*radians*/,
$GMST0 /*decimal hours*/, $ndays,
$lat /*decimal degrees*/, $lon /*decimal degrees*/,
$signh /*+1 for set, -1 for rise*/,
&$az /* radians */ ) {
global $torad, $tzone;
//calculate equatorial coords for Sun
//RA and DEC are returned with decimal hours and radians
SunEqCoords( $t+$tzone, $ndays, $alpha /*decimal hours*/, $dec /*radians*/ );
//calculate rise/set hour angle (decimal hours) and azimuth (radians)
hour_az( $z, $dec, $lat*$torad, $signh, $h, $az );
//get local time in decimal hours from hour angle
$tlocal = getlocaltime( $h, $alpha, $GMST0, $lon/15.0 );
return $tlocal;
}
//-----------------------------------------------------------------------------
// calculate local transit time in decimal hours of the Sun on
// date ndays after J2000.0
//-----------------------------------------------------------------------------
function Sun_transit_time( $lon /*decimal degrees*/, $ndays ) {
global $tzone, $ratio, $terror;
$countmax = 15; //max number of iterations
$tnew = 12; //initial guess (local noon)
$told = 0;
$count = 0;
while( abs($tnew - $told) > $terror && $count < $countmax ) {
$told = $tnew;
//get the ecliptic coords of the Sun
SunEqCoords( $tnew+$tzone, $ndays, $alpha, $dec );
//calculate the local equation of time (hour angle of true Sun at local noon)
//L = solar mean longitude in deg
$L = 280.460 + 0.9856474*($ndays+$tnew/24.0);
$L = fmod($L,360.0);
$EOT = $L/15.0-$alpha;
if( $EOT > 1 ) $EOT -= 24; //spike occurs near RA = L = 0 and/or 24
//calculate the local transit time
$tnew = 12.0 + ($lon/15.0-$tzone) - $EOT;
++$count;
}
if( $tnew < 0 ) $tnew += 24;
return $tnew;
}
//-----------------------------------------------------------------------------
// return hour angle (decimal hours) and azimuth (radians)
// input args z, dec, and lat must all be in radians
// signh = {-1 for rise, +1 for set}
//-----------------------------------------------------------------------------
function hour_az( $z, $dec, $lat, $signh, &$h, &$az ) {
$torad = M_PI/180.0;
$sz = sin($z);
$cz = cos($z);
$sd = sin($dec);
$cd = cos($dec);
$sphi = sin($lat);
$cphi = cos($lat);
$ch = ($cz-$sd*$sphi)/($cd*$cphi); //cos(h*)
$h = acos($ch); //h* in radians
$saz = -$cd*sin($h)/$sz*$signh;
$caz = ($sd-$sphi*$cz)/($cphi*$sz);
$az = angle($caz,$saz);
$h = $signh*$h/$torad/15; //h = sign(h)*(h*) in decimal hours
}
//-----------------------------------------------------------------------------
// calculate local time in decimal hours from hour angle
// (all args in decimal hours)
//-----------------------------------------------------------------------------
function getlocaltime( $h, $alpha, $GMST0, $lon ) {
global $ratio, $tzone;
$LST = $h + $alpha;
if( $LST < 0 ) $LST += 24;
$dt = ($LST + $lon - $GMST0)/$ratio; //time elapsed since UT midnight
//if the time zone correction to get local time sends the latter
//negative, then the elapsed time since UT midnight is actually
//dt + 24 *sidereal* hours
$localt = $dt - $tzone;
if( $localt < 0 ) $localt += 24/$ratio;
return $localt;
}
//-----------------------------------------------------------------------------
// Print a string containing SunRiseSet() output that looks like this:
// astro az civil az rise az trans alt set az civil az astro az
// 03:59 43 05:19 58 05:50 63 13:06 72 20:20 298 20:51 303 22:12 318
//
// Optional arguments: 'heading_only' and 'skipheading'
//
// A typical call might look like this:
// echo printRiseSets( SunRiseSet( $lat_DC, $lon_DC ) );
//-----------------------------------------------------------------------------
function printRiseSets( $A ) {
$nargs = func_num_args();
$args = func_get_args();
$h_only = false;
$skip_h = false;
if( $nargs > 1 ) {
for( $i=1; $i < $nargs; ++$i ) {
if( $args[$i] == 'heading_only' ) $h_only = true;
if( $args[$i] == 'skipheading' ) $skip_h = true;
}
}
$tstyle = "<span style=\"font-family: monospace; font-size: 11px; line-height: 1.2;\">";
$hstyle = "<span style=\"font-family: monospace; font-size: 11px; line-height: 1.7;\">";
$gstyle = "<span style=\"color:#999999\">";
$bgcolor = "<span style=\"background-color: #f0f0f0;\">";
$endspan = "</span>";
$enddiv = "</div>";
$sp = " ";
$heading = "astro".$sp.$sp.$gstyle."az".$endspan.$sp.$sp.
"civil".$sp.$sp.$gstyle."az".$endspan.$sp.$sp.
"rise".$sp.$sp.$sp.$gstyle."az".$endspan.$sp.$sp.
"trans".$sp.$gstyle."alt".$endspan.$sp.
$sp."set".$sp.$sp.$sp.$gstyle."az".$endspan.$sp.$sp.
"civil".$sp.$sp.$gstyle."az".$endspan.$sp.$sp.
"astro".$sp.$sp.$gstyle."az".$endspan."<br>";
if( $h_only ) return $hstyle.$heading.$endspan;
$times = "";
for( $i=0; $i < 7; ++$i ) {
$tstr = printtime($A[$i]["t"],false);
$az = (int)($A[$i]["a"]+0.5);
if( $az == 0 ) $az = "---";
if( $i != 3 ) {
if( strlen($az) == 2 ) $az = $sp.$az;
if( strlen($az) == 1 ) $az = $sp.$sp.$az;
}
if( $i == 2 || $i == 4 ) {
$times = $times.$bgcolor.$tstr.$endspan.$sp.$gstyle.$az.$endspan;
} else {
$times = $times.$tstr.$sp.$gstyle.$az.$endspan;
}
if( $i == 3 ) {
if( strlen($az) == 1 ) $times = $times.$sp;
if( strlen($az) == 3 ) $times = $times.$sp;
}
if( $i != 3 || strlen($az) != 3 ) {
$times = $times.$sp.$sp;
}
}
$nargs = func_num_args();
$args = func_get_args();
if( $skip_h ) {
$msg = $tstyle.$times.$endspan;
} else {
$msg = $hstyle.$heading.$endspan.$tstyle.$times.$endspan;
}
return $msg;
}
//=============================================================================
// other (independent) utility functions
//=============================================================================
//-----------------------------------------------------------------------------
// calculate equatorial coordinates of the Sun for UT time t on
// date ndays after J2000.0 (See 2006 Astronomical Almanac, p. C24)
//-----------------------------------------------------------------------------
function SunEqCoords( $UT /*decimal hours*/, $ndays /*whole days from JD2451545.0*/,
&$alpha /*decimal hours*/, &$dec /*radians*/ )
{
$torad = M_PI/180.0;
$dt = $ndays + $UT/24.0;
$L = 280.460 + 0.9856474*$dt; //solar mean longitude in deg
$M = 357.528 + 0.9856003*$dt; //solar mean anomaly in deg
$L = fmod($L,360.0);
$M = fmod($M,360.0);
$eL = $L + 1.915*sin($M*$torad)
+ 0.020*sin(2*$M*$torad); //solar ecliptic longitude in deg
if( $eL < 0 ) $eL += 360;
$T = $dt/36525.0;
$eps = 23.43929111 - //obliquity of the ecliptic (deg)
(0.01301118 + (0.00000016393 - 0.0000005036*$T)*$T)*$T;
// $R = 1.00014 - 0.01671*cos($M*$torad) - 0.00014*cos(2*$M*$torad); //distance of Sun in AU
// $sd = 0.2666/$R; //apparent solar radius in degrees
//RA and DEC of Sun at time t
$dec = asin( sin($eL*$torad)*sin($eps*$torad) ); //radians
$ca = cos($eL*$torad)/cos($dec);
$sa = sin($eL*$torad)*cos($eps*$torad)/cos($dec);
$alpha = angle($ca,$sa)/$torad/15; //decimal hours
}
//-----------------------------------------------------------------------------
// calculate zenith angle and azimuth
//-----------------------------------------------------------------------------
function zaz( $h /*decimal hours*/, $dec /*radians*/, $lat /*decimal degrees*/,
&$z /*radians*/, &$az /*radians*/ )
{
$torad = M_PI/180.0;
$slat = sin($lat*$torad);
$clat = cos($lat*$torad);
$sdec = sin($dec);
$cdec = cos($dec);
$hrad = $h*15.0*$torad;
$sh = sin($hrad);
$ch = cos($hrad);
$cz = $slat*$sdec + $clat*$cdec*$ch;
$z = acos($cz);
$sz = sin($z);
// $sz = sqrt(sqr($cdec*$sh)+sqr($sdec*$clat-$slat*$cdec*$ch));
// $z = asin($sz);
$saz = -$cdec*$sh/$sz;
$caz = ($sdec-$slat*$cz)/($clat*$sz);
$az = angle($caz,$saz);
}
//-----------------------------------------------------------------------------
// return the hour angle in decimal hours
//
// lon must be decimal degrees, UT and RA in decimal hourse
//
// optional arguments (include either all three or none):
// integer day of month, month number (January=1), 4-digit year
//-----------------------------------------------------------------------------
function hourangle( $UT /*decimal hours*/,
$RA /*decimal hours*/,
$lon /*decimal degrees West longitude*/ )
{
//if no date provided, assume today
$nargs = func_num_args();
$args = func_get_args();
if( $nargs > 3 ) {
if( $nargs != 6 ) {
print("<p>hourangle error: wrong number of arguments!</p>");
return;
}
$day = $args[3];
$month = $args[4];
$year = $args[5];
} else {
$today = getdate(time()); //use PHP function to get current local date
$day = $today['mday']; //day of month
$month = $today['mon']; //month number
$year = $today['year'];
}
//JD at 0h UT
$jd0 = JD0hUT($lon,$tzone,$day,$month,$year);
//Greenwich mean sideral time at 0h UT
$GMST0 = GMST0hUT($jd0);
//hour angle
$ratio = 1.00273790935; //sidereal days per solar day
$h = $GMST0 + $ratio*$UT - $lon/15.0 - $RA;
while( $h < 0 ) $h += 24.0;
while( $h > 0 ) $h -= 24.0;
if( $h < -12 ) $h += 24;
if( $h > 12 ) $h -= 24;
return $h;
}
//-----------------------------------------------------------------------------
// return the Greenwich mean sidereal time (decimal hours) at 0h UT on
// a given Julian Date
//-----------------------------------------------------------------------------
function GMST0hUT( $jd0h )
{
$ndays = $jd0h - 2451545.0; //number of days from J2000.0 = Jan. 1.5 UT 2000
$T = $ndays/36525.0; //Julian centuries from J2000.0
$GMST0 = 24110.54841 + (8640184.812866 + (0.093104 - 0.0000062*$T)*$T)*$T;
$GMST0 /= 3600.0;
while( $GMST0 < 0 ) $GMST0 += 24.0;
while( $GMST0 > 24 ) $GMST0 -= 24.0;
return $GMST0;
}
//-----------------------------------------------------------------------------
// return the Julian Date at 0h UT
// also calculates the time zone (an integer, positive for West longitude)
//
// lon must be decimal degrees
//
// optional arguments (include either all three or none):
// integer day of month, month number (January=1), 4-digit year
//-----------------------------------------------------------------------------
function JD0hUT( $lon /*decimal degrees West longitude*/, &$tzone )
{
//if no date provided, assume today
$nargs = func_num_args();
$args = func_get_args();
if( $nargs > 2 ) {
if( $nargs != 5 ) {
print("<p>JD0hUT error: wrong number of arguments!</p>");
return;
}
$day = $args[2];
$month = $args[3];
$year = $args[4];
} else {
$today = getdate(time()); //use PHP function to get current local date
$day = $today['mday']; //day of month
$month = $today['mon']; //month number
$year = $today['year'];
}
//determine the local time zone
$tzone = get_time_zone( $lon, $day, $month, $year );
//now that time zone is determined, we can get the UT date
$utime = mktime(1,1,1,$month,$day,$year);
$UT = getdate($utime+$tzone*3600); //use PHP function to get current UT date
$day = $UT['mday']; //day of month
$month = $UT['mon']; //month number
$year = $UT['year'];
//finally, calculate the JD at 0h UT
$jd0 = gregoriantojd($month,$day,$year)-0.5;
return $jd0;
}
//-----------------------------------------------------------------------------
// determine the time zone
//-----------------------------------------------------------------------------
function get_time_zone( $lon /*decimal degrees west longitude*/,
$day /*integer day of month*/,
$month /*month number (January=1)*/,
$year /*4-digit year*/ )
{
$tzone = (int)($lon/15.0+0.5);
if( is_dst( $lon, $day, $month, $year ) ) $tzone -= 1;
return $tzone;
}
//-----------------------------------------------------------------------------
// determine if daylight savings time or not
//-----------------------------------------------------------------------------
function is_dst( $lon /*decimal degrees west longitude*/,
$day /*integer day of month*/,
$month /*month number (January=1)*/,
$year /*4-digit year*/ )
{
if( $month > 3 && $month < 11 ) { //adjust tzone for daylight savings
if( $month == 4 and $day < 7 ) { //check for first Sunday in April
$utime = mktime (1,1,1,$month,1,$year);
$day1 = date('w',$utime); //day of week of first day of April
$sunday = 1 + (7-$day1); //day of month of 1st Sunday
if( $day >= $sunday ) {
return true;
}
} else if( $month == 10 and $day > 24 ) { //check for last Sunday in October
$utime = mktime (1,1,1,$month,31,$year);
$day31 = date('w',$utime); //day of week of last day of October
$sunday = 31 - $day31; //day of month of last Sunday
if( $day < $sunday ) {
return true;
}
} else {
return true;
}
}
return false;
}
?>
