Sunrise Sunset Algorithm Example.txt

Source: https://edwilliams.org/sunrise_sunset_example.htm

Source:
	Almanac for Computers, 1990
	published by Nautical Almanac Office
	United States Naval Observatory
	Washington, DC 20392

Inputs:
	day, month, year:      date of sunrise/sunset
	latitude, longitude:   location for sunrise/sunset
	zenith:                Sun's zenith for sunrise/sunset
	  offical      = 90 degrees 50'
	  civil        = 96 degrees
	  nautical     = 102 degrees
	  astronomical = 108 degrees
	
	NOTE: longitude is positive for East and negative for West

Worked example (from book):
	June 25, 1990:	25, 6, 1990
	Wayne, NJ:      40.9, -74.3
	Office zenith:  90 50' cos(zenith) = -0.01454


1. first calculate the day of the year

	N1 = floor(275 * month / 9)
	N2 = floor((month + 9) / 12)
	N3 = (1 + floor((year - 4 * floor(year / 4) + 2) / 3))
	N = N1 - (N2 * N3) + day - 30
	
	Example:
	N1 = 183
	N2 = 1
	N3 = 1 + floor((1990 - 4 * 497 + 2) / 3)
	   = 1 + floor((1990 - 1988 + 2) / 3)
	   = 1 + floor((1990 - 1988 + 2) / 3)
	   = 1 + floor(4 / 3)
	   = 2
	N = 183 - 2 + 25 - 30 = 176

2. convert the longitude to hour value and calculate an approximate time

	lngHour = longitude / 15
	
	if rising time is desired:
	  t = N + ((6 - lngHour) / 24)
	if setting time is desired:
	  t = N + ((18 - lngHour) / 24)
	
	Example:
	lngHour = -74.3 / 15 = -4.953
	t = 176 + ((6 - -4.953) / 24)
	  = 176.456
	
3. calculate the Sun's mean anomaly
	
	M = (0.9856 * t) - 3.289
	
	Example:
	M = (0.9856 * 176.456) - 3.289
	  = 170.626

4. calculate the Sun's true longitude
       [Note throughout the arguments of the trig functions
       (sin, tan) are in degrees. It will likely be necessary to
       convert to radians. eg sin(170.626 deg) =sin(170.626*pi/180
       radians)=0.16287]
	
	L = M + (1.916 * sin(M)) + (0.020 * sin(2 * M)) + 282.634
	NOTE: L potentially needs to be adjusted into the range [0,360) by adding/subtracting 360
	
	Example:
	L = 170.626 + (1.916 * sin(170.626)) + (0.020 * sin(2 * 170.626)) + 282.634
	  = 170.626 + (1.916 * 0.16287) + (0.020 * -0.32141) + 282.634
	  = 170.626 + 0.31206 + -0.0064282 + 282.634
	  = 453.566 - 360
	  = 93.566

5a. calculate the Sun's right ascension
	
	RA = atan(0.91764 * tan(L))
	NOTE: RA potentially needs to be adjusted into the range [0,360) by adding/subtracting 360
	
	Example:
	RA = atan(0.91764 * -16.046)
	   = atan(0.91764 * -16.046)
	   = atan(-14.722)
	   = -86.11412

5b. right ascension value needs to be in the same quadrant as L

	Lquadrant  = (floor( L/90)) * 90
	RAquadrant = (floor(RA/90)) * 90
	RA = RA + (Lquadrant - RAquadrant)
	
	Example:
	Lquadrant  = (floor(93.566/90)) * 90
	           = 90
	RAquadrant = (floor(-86.11412/90)) * 90
	           = -90
	RA = -86.11412 + (90 - -90)
	   = -86.11412 + 180
	   = 93.886

5c. right ascension value needs to be converted into hours

	RA = RA / 15
	
	Example:
	RA = 93.886 / 15
	   = 6.259

6. calculate the Sun's declination

	sinDec = 0.39782 * sin(L)
	cosDec = cos(asin(sinDec))
	
	Example:
	sinDec = 0.39782 * sin(93.566)
	       = 0.39782 * 0.99806
	       = 0.39705
	cosDec = cos(asin(0.39705))
	       = cos(asin(0.39705))
	       = cos(23.394)
	       = 0.91780

7a. calculate the Sun's local hour angle
	
	cosH = (cos(zenith) - (sinDec * sin(latitude))) / (cosDec * cos(latitude))
	
	if (cosH >  1) 
	  the sun never rises on this location (on the specified date)
	if (cosH < -1)
	  the sun never sets on this location (on the specified date)
	
	Example:
	cosH = (-0.01454 - (0.39705 * sin(40.9))) / (0.91780 * cos(40.9))
	     = (-0.01454 - (0.39705 * 0.65474)) / (0.91780 * 0.75585)
	     = (-0.01454 - 0.25996) / 0.69372
	     = -0.2745 / 0.69372
	     = -0.39570

7b. finish calculating H and convert into hours
	
	if if rising time is desired:
	  H = 360 - acos(cosH)
	if setting time is desired:
	  H = acos(cosH)
	
	H = H / 15
	
	Example:
	H = 360 - acos(-0.39570)
	  = 360 - 113.310   [ note result of acos converted to degrees]
	  = 246.690
	H = 246.690 / 15
	  = 16.446

8. calculate local mean time of rising/setting
	
	T = H + RA - (0.06571 * t) - 6.622
	
	Example:
	T = 16.446 + 6.259 - (0.06571 * 176.456) - 6.622
	  = 16.446 + 6.259 - 11.595 - 6.622
	  = 4.488

9. adjust back to UTC
	
	UT = T - lngHour
	NOTE: UT potentially needs to be adjusted into the range [0,24) by adding/subtracting 24
	
	Example:
	UT = 4.488 - -4.953
	   = 9.441
	   = 9h 26m

10. convert UT value to local time zone of latitude/longitude
	
	localT = UT + localOffset
	
	Example:
	localT = 9h 26m + -4
	       = 5h 26m
	       = 5:26 am EDT

Sunrise Sunset Algorithm.txt

Source: https://edwilliams.org/sunrise_sunset_algorithm.htm
Archive: https://web.archive.org/web/20161022214335/http://williams.best.vwh.net/sunrise_sunset_algorithm.htm

Sunrise/Sunset Algorithm
Source:
	Almanac for Computers, 1990
	published by Nautical Almanac Office
	United States Naval Observatory
	Washington, DC 20392

Inputs:
	day, month, year:      date of sunrise/sunset
	latitude, longitude:   location for sunrise/sunset
	zenith:                Sun's zenith for sunrise/sunset
	  offical      = 90 degrees 50'
	  civil        = 96 degrees
	  nautical     = 102 degrees
	  astronomical = 108 degrees
	
	NOTE: longitude is positive for East and negative for West
        NOTE: the algorithm assumes the use of a calculator with the
        trig functions in "degree" (rather than "radian") mode. Most
        programming languages assume radian arguments, requiring back
        and forth convertions. The factor is 180/pi. So, for instance,
        the equation RA = atan(0.91764 * tan(L)) would be coded as RA
        = (180/pi)*atan(0.91764 * tan((pi/180)*L)) to give a degree
        answer with a degree input for L.


1. first calculate the day of the year

	N1 = floor(275 * month / 9)
	N2 = floor((month + 9) / 12)
	N3 = (1 + floor((year - 4 * floor(year / 4) + 2) / 3))
	N = N1 - (N2 * N3) + day - 30

2. convert the longitude to hour value and calculate an approximate time

	lngHour = longitude / 15
	
	if rising time is desired:
	  t = N + ((6 - lngHour) / 24)
	if setting time is desired:
	  t = N + ((18 - lngHour) / 24)

3. calculate the Sun's mean anomaly
	
	M = (0.9856 * t) - 3.289

4. calculate the Sun's true longitude
	
	L = M + (1.916 * sin(M)) + (0.020 * sin(2 * M)) + 282.634
	NOTE: L potentially needs to be adjusted into the range [0,360) by adding/subtracting 360

5a. calculate the Sun's right ascension
	
	RA = atan(0.91764 * tan(L))
	NOTE: RA potentially needs to be adjusted into the range [0,360) by adding/subtracting 360

5b. right ascension value needs to be in the same quadrant as L

	Lquadrant  = (floor( L/90)) * 90
	RAquadrant = (floor(RA/90)) * 90
	RA = RA + (Lquadrant - RAquadrant)

5c. right ascension value needs to be converted into hours

	RA = RA / 15

6. calculate the Sun's declination

	sinDec = 0.39782 * sin(L)
	cosDec = cos(asin(sinDec))

7a. calculate the Sun's local hour angle
	
	cosH = (cos(zenith) - (sinDec * sin(latitude))) / (cosDec * cos(latitude))
	
	if (cosH >  1) 
	  the sun never rises on this location (on the specified date)
	if (cosH < -1)
	  the sun never sets on this location (on the specified date)

7b. finish calculating H and convert into hours
	
	if if rising time is desired:
	  H = 360 - acos(cosH)
	if setting time is desired:
	  H = acos(cosH)
	
	H = H / 15

8. calculate local mean time of rising/setting
	
	T = H + RA - (0.06571 * t) - 6.622

9. adjust back to UTC
	
	UT = T - lngHour
	NOTE: UT potentially needs to be adjusted into the range [0,24) by adding/subtracting 24

10. convert UT value to local time zone of latitude/longitude
	
	localT = UT + localOffset