计算观测点对太阳指向的方位角度和俯仰角度

在计算观测点对太阳的指向时，网上相关的公式五花八门，在百度的过程中，我发现很多都是错误的，经过多番的验证，最终参考了两位大佬的博客，也对程序做了一下改动。 

【公式推导+matlab代码】太阳位置（太阳方位角和太阳高度角）计算_太阳方位角计算公式详解-优快云博客【C语言】计算实时太阳角度（高度角、方位角），以及使用stm32单片机实时获取时间戳_输入当地经纬度信息 计算得出当地太阳高度角和方位角 c语言代码实现-优快云博客

首先是观测点对于太阳指向方位和俯仰的计算公式：

%%计算公式
%%太阳方位角相对于正北方向的夹角为：
%%%As=arccos((sinhssinL-sindeta)/coshscosL)  Ts<0
%%%As=360-arccos((sinhssinL-sindeta)/coshscosL)  Ts>0
%%%太阳高度角为
%%%hs=arcsin(sindetashinL+cosdetacosLcosTs)
%%其中L为观测点的地理纬度，deta和Ts分别为太阳的赤纬角和太阳时角

然后是时间的获取，这里面有两种方式，一个是通过matlab自动获取目前的系统时间，一个是如果有外接的时统端机，可以通过时统端机不断地获取最新时间，时间点是从2000/1/1 0:00:00.000开始的UTC时间秒。

首先是通过matlab获取时间

%获取当前时间，该时间为世界时，比北京时间早八个小时
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, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31];
%sumday表示积日，即一年中第几天
sumday = sum(day_tab(1:month1)) + day1;
dn = double(sumday); 

然后是通过时统获取UTC时间秒

%%
%%求解年份
Timestamp=766162398+12*3600;   %%这里以北京时间2024/4/12 10:53:18.000为例子
Year_data = floor(Timestamp/31556736+2000);					%%用时间戳获取年份 
DayData =floor(mod(Timestamp,31556736)/86400+1);		%%获取日期 
Hour_data = floor(mod((Timestamp/3600),24));				   %获取小时数 
if (Hour_data==24)
    Hour_data = 0;										%%当时间为24点时归零 
end
minute=mod((Timestamp/60),60);
time_hour = Hour_data+mod((Timestamp/60),60)/60.0;					%%获取分钟信息并且转为小数形式加入小时数中

输入当地的纬度

%L1为所在地的纬度（弧度制）
L1 = 43.976850 * pi / 180.0;  %%此处为长春的纬度

计算当日的赤纬角，这里面参考了三个公式，测试发现均可用

% %%
% %%%ds为当日的太阳赤纬角（角度制）-----采用第三个公式
% ds=23.45*sin((sumday-80)/370.0*360*pi/180);

%ds为当日的太阳赤纬角（角度制）-----采用第一个公式
%Bourges太阳赤纬角算法
%n0=78.801+(0.2422*(year1-1969))-round(0.25*(year1-1969));
 n0=79.676+(0.2422*(year1-1985))-round(0.25*(year1-1985));
%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为当日的太阳赤纬角（角度制）-----采用第二个公式
% % 一年中某天的太阳赤纬角的角度是固定的，范围从（-26.5，26.5）
% d=sumday;  %%一年中天数的累加和
% b=2*pi*(d-1)/365.2422;
% 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)+0.00148*sin(3*b))*180/pi;

赤纬角转换为弧度

%%
%ds转化为弧度
ds=deg2rad(ds);
%ds转换为度
ds_deg=rad2deg(ds);

考虑时间差和不考虑时间差的计算，这里面要注意的是，经度需要修改为当地的经度，否则会引起计算偏差，作者的位置是长春，经纬度为（125.394173,43.986850）

%%
%考虑时间差的计算
%%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 - 125.394173) * 4.0) / 60.0;
%st为真太阳时
st=Sd+Ts/60;

% % %%
% %%不考虑时间差的计算，目前和软件比较接近的结果
% st=hour1+minute1/60+8;   %%太阳时间

观测点对太阳指向的方位和俯仰

%%
%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

验证阶段，作者使用的是orbitron软件，平时使用这个软件计算对卫星指向，特别实用，强烈推荐！！！

验证之后的误差在0.5度以内。