【公式推导+matlab代码】太阳位置（太阳方位角和太阳高度角）计算

一、计算公式

二、太阳高度角计算

太阳高度角的求解应先求解太阳赤纬角ds和太阳时角ts。

2.1太阳时角

定义sd为平太阳时；et为真太阳时时差；st为真太阳时；ts为太阳时角。

		  double sd=hour1+(minute1+(120-116.33)*4)/60;
          double et=0.0028-1.9587*sin(b)+9.9059*sin(2*b)-7.0924*cos(b)-0.6882*cos(2*b);
          double st=sd+et/60;
		  double ts=(st-12.0)*PI/12;

2.2太阳赤纬角

定义ds为太阳赤纬角。

	 double d=sumday;
	 double b=2*PI*(d-1)/365.2422;
     double ds=0.006918+0.070257*sin(b)+0.000907*sin(2*b)+0.00148*sin(3*b)-0.39912*cos(b)-0.006758*cos(2*b)-0.002697*cos(3*b);

三、代码示例

clc
clear
%获取当前时间，该时间为世界时，比北京时间早八个小时
currentTime = datetime('now', 'Format', 'yyyy-MM-dd HH:mm:ss', 'TimeZone', 'UTC');
year1 = year(currentTime);
month1 = month(currentTime);
day1 = day(currentTime);
hour1 = hour(currentTime);
minute1 = minute(currentTime);
day_tab = [0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31];
%sumday表示积日，即一年中第几天
sumday = sum(day_tab(1:month1)) + day1;
dn = double(sumday);

%L1为所在地的维度（弧度制）
L1 = 40.60 * pi / 180.0;
%Bourges太阳赤纬角算法
n0=78.801+(0.2422*(year1-1969))-round(0.25*(year1-1969));
%b1为日角
b=2*pi*(dn-1-n0)/365.2422;
%ds为当日的太阳赤纬角（角度制）
% 一年中某天的太阳赤纬角的角度是固定的，范围从（-26.5，26.5）
ds=0.3723+23.2567*sin(b)+0.1149*sin(2.0*b)-0.1712*sin(3.0*b)-0.758*cos(b)...
    +0.3656*cos(2.0*b)+0.02010*cos(3.0*b);
%ds转化为弧度
ds=deg2rad(ds);
%Ts表示真太阳时减去平太阳时的时差
Ts=0.0028-1.9587*sin(b)+9.9059*sin(2*b)-7.0924*cos(b)-0.6882*cos(2*b);
%hour1+8表示北京时间，后面为计算当地时间
Sd = hour1 + 8.0 + (minute1 - (120.0 - 114.683) * 4.0) / 60.0;
%st为真太阳时
st=Sd+Ts/60;
%ts为所在地的太阳时角
ts = (st - 12.0) * pi / 12.0;
%hs为太阳高度角
hs = asin(sin(ds) * sin(L1) + cos(ds) * cos(L1) * cos(ts));
hs1=rad2deg(hs);
%as1为地理正北坐标系下的太阳方位角
as1 = acos((sin(hs) * sin(L1) - sin(ds)) / (cos(hs) * cos(L1)))*180/pi;
%此时方位角是以南为零方位角，转换为以北为零方位角
as1=180-as1;
%进行角度范围转换
if ts > 0
    as1 =360-as1;
end