How to calculate sunrise and sunset times?
By Michael Henderson •
I need to create a function (in C++) to calculate the sunrise and sunset times, but I am not a mathematician and I cannot find a correct (and easy) way to do that.
I need to get the same results as can be found in:
and
I tried to implement a function based on these articles and but the results are wrong. (maybe I am doing something wrong)
Does anyone know a correct (and easy) formula to calculate it? Maybe the formula used by the websites that I mentioned.
Note: values that I have as input: latitude, longitude, date and UTC offset. (I don't have the altitude)
Thanks
Update:
I developed a new function on Matlab that seems to be more accurate but I still not get the exact sunrise and sunset times:
% Parameters definition
lat = -23.545570; % Latitude
lng = -46.704082; % Longitude
UTCoff = -3; % UTC offset
nDays = daysact('01-jan-2017', '15-mar-2017'); % Number of days since 01/01
% Longitudinal correction
longCorr = 4*(lng - 15*UTCoff);
B = 360*(nDays - 81)/365; % I have no idea
% Equation of Time Correction
EoTCorr = 9.87*sind(2*B) - 7.53*cosd(B) - 1.5*sind(B);
% Solar correction
solarCorr = longCorr - EoTCorr;
% Solar declination
delta = asind(sind(23.45)*sind(360*(nDays - 81)/365));
sunrise = 12 - acosd(-tand(lat)*tand(delta))/15 - solarCorr/60;
sunset = 12 + acosd(-tand(lat)*tand(delta))/15 - solarCorr/60;
sprintf('%2.0f:%2.0f:%2.0f\n', degrees2dms(sunrise))
sprintf('%2.0f:%2.0f:%2.0f\n', degrees2dms(sunset))This function gives me the sunrise at 05:51:25 when it should be 06:09 and the sunset as 18:02:21 when it should be 18:22, according to ESRL (NOAA).
The function was developed based on this:
Any idea what can I do to improve the accuracy?
$\endgroup$ 32 Answers
$\begingroup$Here is a complete routine to that calculates the sunrise (optionally the sunset, see the code) in C++. It only requires latitude and longitude as input, and is accurate to within seconds for the civil sunrise time. This code is running in a UWP C++ app created with Visual Studio 2017. It includes a subroutine to fill the output minutes and seconds with zeros in the case of single digit results.
String^ toplatformstring(bool fill, long ll) { // convert ll to platform string Platform::String^ p_string; std::string doit; if (fill == false) { doit = std::to_string(ll); // convert, don't fill with zeros } else { //convert ll to std doit and fill with zeros std::stringstream ss; ss << std::setw(2) << std::setfill('0') << ll; doit = ss.str(); } //convert doit to platform string char const *pchar = doit.c_str(); std::string s_str = std::string(pchar); std::wstring wid_str = std::wstring(s_str.begin(), s_str.end()); const wchar_t* w_char = wid_str.c_str(); p_string = ref new Platform::String(w_char); return p_string;
}
//double la = 39.299236; // baltimore
//double lo = -76.609383;
//double la = 37.0; // SF California
//double lo = -122.0;
Platform::String^ MainPage::sunrise(double la, double lo) { /*double la = 39.300213; double lo = -76.610516;*/ Platform::String^ sunrisetime; //// get year, month, day integers time_t rawtime; struct tm timeinfo; // get date and time info time(&rawtime); localtime_s(&timeinfo, &rawtime); double xday = timeinfo.tm_mday; double xmonth = timeinfo.tm_mon; xmonth = xmonth + 1; // correct for origin 0 //textblockc->Text = xmonth.ToString(); // for debugging double xyear = timeinfo.tm_year; //double dayofyear = timeinfo.tm_yday; // day of year also // calculate the day of the year // N1 = floor(275 * month / 9) double xxN1 = floor(275 * xmonth / 9); // N2 = floor((month + 9) / 12) double xxN2 = floor((xmonth + 9) / 12); // N3 = (1 + floor((year - 4 * floor(year / 4) + 2) / 3)) double xxN3 = (1 + floor((xyear - 4 * floor(xyear / 4) + 2) / 3)); // N = N1 - (N2 * N3) + day - 30 double day = xxN1 - (xxN2 * xxN3) + xday - 30; double zenith = 90.83333333333333; double D2R = M_PI / 180; double R2D = 180 / M_PI; // convert the longitude to hour value and calculate an approximate time double lnHour = lo / 15; double t; //if (sunrise) { t = day + ((6 - lnHour) / 24); //} else { //t = day + ((18 - lnHour) / 24); //}; //calculate the Sun's mean anomaly double M = (0.9856 * t) - 3.289; //calculate the Sun's true longitude double L = M + (1.916 * sin(M * D2R)) + (0.020 * sin(2 * M * D2R)) + 282.634; if (L > 360) { L = L - 360; } else if (L < 0) { L = L + 360; }; //calculate the Sun's right ascension double RA = R2D * atan(0.91764 * tan(L * D2R)); if (RA > 360) { RA = RA - 360; } else if (RA < 0) { RA = RA + 360; }; //right ascension value needs to be in the same qua double Lquadrant = (floor(L / (90))) * 90; double RAquadrant = (floor(RA / 90)) * 90; RA = RA + (Lquadrant - RAquadrant); //right ascension value needs to be converted into hours RA = RA / 15; //calculate the Sun's declination double sinDec = 0.39782 * sin(L * D2R); double cosDec = cos(asin(sinDec)); //calculate the Sun's local hour angle double cosH = (cos(zenith * D2R) - (sinDec * sin(la * D2R))) / (cosDec * cos(la * D2R)); double H; //if (sunrise) { H = 360 - R2D * acos(cosH); //} else { //H = R2D * Math.acos(cosH) //}; H = H / 15; //calculate local mean time of rising/setting double T = H + RA - (0.06571 * t) - 6.622; //adjust back to UTC double UT = T - lnHour; if (UT > 24) { UT = UT - 24; } else if (UT < 0) { UT = UT + 24; } //convert UT value to local time zone of latitude/longitude int offset = (int)(lo / 15); // estimate utc correction double localT = UT + offset; // -5 for baltimore //convert to seconds int seconds = (int)(localT * 3600); long sec = seconds % 60; long minutes = seconds % 3600 / 60; long hours = seconds % 86400 / 3600; hours = hours % 12; Platform::String^ ssec = toplatformstring(true, sec); Platform::String^ mminutes = toplatformstring(true, minutes); Platform::String^ hhours = toplatformstring(false, hours); sunrisetime = hhours + ":" + mminutes + ":" + ssec; return sunrisetime;
}As @Richard already answered here, I was mixing things.
- My Matlab script is calculating the actual sunrise and sunset (geometrically).
- The NOAA website gives the apparent sunrise and sunset. These values are corrected for atmospheric refraction.
In the glossary to the NOAA website, it is written:
Due to atmospheric refraction, sunrise occurs shortly before the sun crosses above the horizon. Light from the sun is bent, or refracted, as it enters earth's atmosphere. See Apparent Sunrise Figure. This effect causes the apparent sunrise to be earlier than the actual sunrise. Similarly, apparent sunset occurs slightly later than actual sunset.
So this is exactly the effect that is causing the 'calculation error'.
@Richard have reverse engineered the functions from the Excel sheet provided on NOAA's website and created a Matlab function to calculate it:
function [sun_rise_set, varargout] = sunRiseSet( lat, lng, UTCoff, date, PLOT)
%SUNRISESET Compute apparent sunrise and sunset times in seconds.
% sun_rise_set = sunRiseSet( lat, lng, UTCoff, date, PLOT) Computes the *apparent** (refraction
% corrected) sunrise and sunset times in seconds from mignight and returns them as
% sun_rise_set. lat and lng are the latitude (+ to N) and longitude (+ to E), UTCoff is the
% local time offset to UTC in hours and date is the date in format 'dd-mmm-yyyy' ( see below for
% an example). Set PLOT to true to create some plots.
%
% [sun_rise_set, noon] = sunRiseSet( lat, lng, UTCoff, date, PLOT) additionally returns the
% solar noon in seconds from midnight.
%
% [sun_rise_set, noon, opt] = sunRiseSet( lat, lng, UTCoff, date, PLOT) additionally returns the
% information opt, which contains information on every second of the day:
% opt.elev_ang_corr : Apparent (refraction corrected) solar elevation in degrees
% opt.azmt_ang : Solar azimuthal angle (deg cw from N)
% opt.solar_decl : Solar declination in degrees
%
% EXAMPLE:
% lat = -23.545570; % Latitude
% lng = -46.704082; % Longitude
% UTCoff = -3; % UTC offset
% date = '15-mar-2017';
%
% [sun_rise_set, noon, opt] = sunRiseSet( lat, lng, UTCoff, date, 1);
%
%
% Richard Droste
%
% Reverse engineered from the NOAA Excel:
% ()
%
% The formulas are from:
% Meeus, Jean H. Astronomical algorithms. Willmann-Bell, Incorporated, 1991.
% Process input
nDays = daysact('30-dec-1899', date); % Number of days since 01/01
nTimes = 24*3600; % Number of seconds in the day
tArray = linspace(0,1,nTimes);
% Compute
% Letters correspond to colums in the NOAA Excel
E = tArray;
F = nDays+2415018.5+E-UTCoff/24;
G = (F-2451545)/36525;
I = mod(280.46646+G.*(36000.76983+G*0.0003032),360);
J = 357.52911+G.*(35999.05029-0.0001537*G);
K = 0.016708634-G.*(0.000042037+0.0000001267*G);
L = sin(deg2rad(J)).*(1.914602-G.*(0.004817+0.000014*G))+sin(deg2rad(2*J)).* ... (0.019993-0.000101*G)+sin(deg2rad(3*J))*0.000289;
M = I+L;
P = M-0.00569-0.00478*sin(deg2rad(125.04-1934.136*G));
Q = 23+(26+((21.448-G.*(46.815+G.*(0.00059-G*0.001813))))/60)/60;
R = Q+0.00256*cos(deg2rad(125.04-1934.136*G));
T = rad2deg(asin(sin(deg2rad(R)).*sin(deg2rad(P))));
U = tan(deg2rad(R/2)).*tan(deg2rad(R/2));
V = 4*rad2deg(U.*sin(2*deg2rad(I))-2*K.*sin(deg2rad(J))+4*K.*U.*sin(deg2rad(J)).* ... cos(2*deg2rad(I))-0.5.*U.*U.*sin(4*deg2rad(I))-1.25.*K.*K.*sin(2.*deg2rad(J)));
AB = mod(E*1440+V+4*lng-60*UTCoff,1440);
if AB/4 < 0, AC = AB/4+180;else, AC = AB/4-180; end
AD = rad2deg(acos(sin(deg2rad(lat))*sin(deg2rad(T))+cos(deg2rad(lat))*cos(deg2rad(T)).*... cos(deg2rad(AC))));
W = rad2deg(acos(cos(deg2rad(90.833))./(cos(deg2rad(lat))*cos(deg2rad(T))) ... -tan(deg2rad(lat))*tan(deg2rad(T))));
X = (720-4*lng-V+UTCoff*60)*60;
% Results in seconds
[~,noon] = min(abs(X - nTimes*tArray));
[~,sunrise] = min(abs(X-round(W*4*60) - nTimes*tArray));
[~,sunset] = min(abs(X+round(W*4*60) - nTimes*tArray));
% Results in degrees
if nargout > 2 || PLOT solar_decl = T; elev_ang_corr = 90-AD; AC_ind = AC > 0; azmt_ang = mod(rad2deg(acos(((sin(deg2rad(lat))*cos(deg2rad(AD)))-sin(deg2rad(T)))./ ... (cos(deg2rad(lat))*sin(deg2rad(AD)))))+180,360); azmt_ang_2 = mod(540-rad2deg(acos(((sin(deg2rad(lat))*cos(deg2rad(AD)))-sin(deg2rad(T)))./ ... (cos(deg2rad(lat))*sin(deg2rad(AD))))),360); azmt_ang(~AC_ind) = azmt_ang_2(~AC_ind);
end
% Print in hours, minutes and seconds
fprintf('Sunrise: %s \nSunset: %s\n', ... datestr(sunrise/nTimes,'HH:MM:SS'), datestr(sunset/nTimes,'HH:MM:SS'));
sun_rise_set = [sunrise sunset];
if nargout > 1 varargout{1} = noon;
end
if nargout > 2 opt.elev_ang_corr = elev_ang_corr; opt.azmt_ang = azmt_ang; opt.solar_decl = solar_decl; varargout{2} = opt;
end
if PLOT figure; hold on plot(linspace(0,24,nTimes), elev_ang_corr); xlabel('Hour'), ylabel('Angle [Deg]') xlim([0 24]), grid on title('Corrected Elevation Angle') figure; plot(linspace(0,24,nTimes), azmt_ang); xlabel('Hour'), ylabel('Angle [Deg]') xlim([0 24]), grid on title('Azimuthal Angle')
endEdit: Richard's uploaded an extended version on Matlab File Exchange
See the complete discussion here.
Now, I can use this Matlab function to convert it to a C++ function.
$\endgroup$ 3
