根据观测时间，经纬度，求太阳高度角

出处：张闯, 吕东辉, 顼超静. 太阳实时位置计算及在图像光照方向中的应用[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))

验证数据：