How to calculate sunrise and sunset times (matlab)?

Question

I need to calculate the sunrise and sunset times in Matlab, but I a 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 already 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)

I also developed a script in 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

What can I do to improve the accuracy and get the same values from the ESRL (NOAA)?

Answer 1

So, I have reverse engineered the functions from the Excel sheet provided on NOAA's website.

Here you go. It computes the apparent (refraction-corrected) sunrise and sunset, accurate like a Swiss watch:

function sun_rise_set = sunRiseSet( lat, lng, UTCoff, date)
%SUNRISESET Compute apparent sunrise and sunset times in seconds.
%     sun_rise_set = sunRiseSet( lat, lng, UTCoff, date) 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).
% 
% EXAMPLE:
%     lat = -23.545570;     % Latitude
%     lng = -46.704082;     % Longitude
%     UTCoff = -3;          % UTC offset
%     date = '15-mar-2017';
% 
%     sun_rise_set = sunRiseSet( lat, lng, UTCoff, date);
% 
%     [sr_h, sr_m, sr_s] = hms(sun_rise_set(1));
%     [ss_h, ss_m, ss_s] = hms(sun_rise_set(2));
%
% 
% 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)));
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
[~,sunrise] = min(abs(X-round(W*4*60) - nTimes*tArray));
[~,sunset] = min(abs(X+round(W*4*60) - nTimes*tArray));

% 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];

Edit: I have uploaded an extended version on Matlab File Exchange

Answer 2

You are mixing apples with oranges here!

• The formula you are using 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 you are observing with your calculation error.

If you really want to calculate the apparent sunrise and sunset, refer to the Solar Calculation Details from NOAA itself or this SO answer. But be aware: "... it's complicated!"

EDIT: See my other answer for a precise function to compute the apparent sunrise and sunset in MatLab