How to calculate sunrise and sunset times?

Question

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:

https://www.esrl.noaa.gov/gmd/grad/solcalc/ and http://sunrise-sunset.org/api

I tried to implement a function based on these articles https://en.wikipedia.org/wiki/Sunrise_equation and http://www.wikihow.com/Estimate-the-Time-of-Sunrise-or-Sunset 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: https://www.mathworks.com/matlabcentral/fileexchange/55509-sunrise-sunset/content/SunriseSunset.mlx

Any idea what can I do to improve the accuracy?

Answer 1

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:
% (https://www.esrl.noaa.gov/gmd/grad/solcalc/calcdetails.html)
% 
% 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')
end

Edit: 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.

Answer 2

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;
}