我需要在 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.0fn', degrees2dms(sunrise))
sprintf('%2.0f:%2.0f:%2.0fn', degrees2dms(sunset))
根据 ESRL (NOAA)，这个函数给我的日出时间是 05：51：25，应该是 06：09，日落应该是 18：22 时是 18：02：21。
该功能是基于这个开发的：https://www.mathworks.com/matlabcentral/fileexchange/55509-sunrise-sunset/content/SunriseSunset.mlx
我可以做些什么来提高精度并从 ESRL (NOAA) 获得相同的值？
你在这里混合苹果和橙子！
-
您使用的公式是计算实际的日出和日落(几何)。
-
NOAA网站给出了明显的日出和日落。这些值针对大气折射进行了校正！
在NOAA网站的词汇表中，它写道：
由于大气折射，日出发生在太阳穿过地平线之前不久。来自太阳的光在进入地球大气层时会弯曲或折射。见表观日出图。这种效应导致表观日出早于实际日出。同样，明显的日落比实际日落稍晚。
所以这正是你观察到的计算误差的效果。
如果您真的想计算日出和日落，请参阅 NOAA 本身的太阳计算细节或这个 SO 答案。但请注意："...太复杂了！
编辑：请参阅我的另一个答案，以获取在MatLab中计算表观日出和日落的精确函数
因此，我对NOAA网站上提供的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: %sn', ...
datestr(sunrise/nTimes,'HH:MM:SS'), datestr(sunset/nTimes,'HH:MM:SS'));
sun_rise_set = [sunrise sunset];
编辑：我已经在Matlab文件交换上上传了一个扩展版本