我需要在 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) 获得相同的值？
最佳答案
因此，我对 NOAA's website 上提供的 Excel 工作表中的函数进行了逆向工程。
干得好。它计算 表观 (折射校正)日出和日落，像瑞士 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];
编辑: 我已经在 Matlab File Exchange 上上传了一个扩展版本