出处:张闯, 吕东辉, 顼超静. 太阳实时位置计算及在图像光照方向中的应用[J]. 电子测量技术, 2010, 33(11):87-89.
算法描述如上,不过有两处错误:
L = (D+E/60)/15; %经度订正
omega = (TT-12)*15; %太阳时角
代码如下:
NF = 2010; %年份
Y = 3; %月
%R = 25; %日
D = 121; %观测处经度的度值
E = 0.4*60; %观测处经度的分度值
phi = 31;%观测地的纬度
%S=16; %观测时的时值
%F=32; %观测时的分值
R = [17 17 19 19 25 25 25];
S = [14 14 9 10 15 15 16];
F = [20 30 56 11 41 47 32];
%% 计算t
N0 = 79.6764+0.2422*(NF-1985)-floor((NF-1985)/4);
%积日N
A = NF/4;
B = A-floor(A);
C = 32.8;
if Y <= 2
C = 30.6;
end
if (B == 0) && (Y > 2)
C = 31.8;
end
N = floor(30.6*Y-C+0.5)+R;
%积日修正值dN
L = (D+E/60)/15; %经度订正
W = (S+F)/60; %时刻订正
dN = (W-L)/24;
t = N-N0+dN
%t = (floor(30.6*Y-C+0.5)+R)+((S+F/60)-(D+E/60)*15)/24-(79.6764+0.2422*(NF-1985)-floor((NF-1985)/4));
%% 日角
theta = 2*pi*t/365.2422;
theta = theta*(pi/180);
%% 太阳赤纬
delta = 0.3723+23.2567*sin(theta)+0.1149*sin(2*theta)-0.1712*sin(3*theta)-0.758*cos(theta)+0.3656*cos(2*theta)+0.0201*cos(3*theta)
%% 时差
Eq = 0.0028-1.9857*sin(theta)+9.9059*sin(2*theta)-7.0924*cos(theta)-0.6882*cos(2*theta);
%% 太阳时角
LC = 4*(D-120);
TT = (S+F/60+LC/60+Eq/60);
omega = (TT-12)*15;
%% 转换为弧度制
delta = delta*(pi/180);
phi = phi*(pi/180);
omega = omega*(pi/180);
%% 太阳高度角
hs = asind(sin(phi)*sin(delta)+cos(phi)*cos(delta).*cos(omega))
验证数据: