如何计算日出和日落时间（matlab）？

Question

我需要在 Matlab 中计算日出和日落时间，但我找不到正确（且简单）的方法。

我需要获得与以下内容相同的结果：

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

我已经尝试根据https://en.wikipedia.org/wiki/Sunrise_equation和http://www.wikihow.com/Estimate-the-Time-of-Sunrise-or-Sunset这几篇文章实现一个功能，但是结果是错误的。 （也许我做错了什么）

我还在 Matlab 中开发了一个似乎更准确的脚本，但我仍然没有得到确切的日出和日落时间：

% 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))

根据ESRL (NOAA)，这个函数给了我 05:51:25 的日出，它应该是 06:09，日落是 18:02:21，它应该是 18:22。

函数就是基于这个开发的：https://www.mathworks.com/matlabcentral/fileexchange/55509-sunrise-sunset/content/SunriseSunset.mlx

我可以做些什么来提高准确性并从 ESRL (NOAA) 获得相同的值？

Answer 1

你在这里把苹果和橙子混在一起！

• 您使用的公式是计算实际日出和日落（几何）。

• NOAA website 提供明显日出和日落。这些值已针对大气折射进行了校正！

the glossary to the NOAA website中写着：

由于大气折射，日出会在太阳穿过地平线之前不久出现。来自太阳的光在进入地球大气层时会发生弯曲或折射。参见明显的日出图。这种效应导致明显的日出早于实际日出。同样，明显的日落比实际的日落稍晚。

所以这正是您在计算错误时观察到的效果。

如果您真的想计算明显的日出和日落，请参考 NOAA 本身的Solar Calculation Details 或this SO answer。但请注意：“……这很复杂！”

编辑：请参阅我的其他答案，了解在 MatLab 中计算明显日出和日落的精确函数

Answer 2

所以，我对NOAA's website 上提供的 Excel 表中的函数进行了逆向工程。

给你。它计算出表观（经折射校正的）日出和日落，精确得像瑞士手表：

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

编辑：我在Matlab File Exchange上传了一个扩展版本