clear
clc
%JD=1.8849556;%108經度?
JD=108.56
FA=0.593412;%34緯度

?i=1;
pi?=?3.14159265;
for?j?=?1:1:365
???N=1:365;%積日
%??N=111;%積日
??for?i=1:47;
???????S?=?1?:0.5:24;
??????%N0?
????%??i=25*j+k;
????????N0=79.6764+0.2422*（2016-1985）-7;?
????????%日角·弧度
????????seita=（2*pi*（N（1j）-N0））/365.2422;%日角·弧度
????????%赤緯角??弧度
????????ED=pi*（0.3723+23.2567*sin（seita）+0.1149*sin（2*seita）-0.1712*sin（3*seita）-0.758*cos（seita）+0.3656*cos（2*seita）+0.0201*cos（3*seita））/180;
????????%地方時計算//小時
????????SD（ji）=S（1i）-（120-JD）/15;
????????%時差??小時
????????ET=（0.0028-1.9857*sin（seita）+9.9057*sin（2*seita）-7.0924*cos（seita）-0.6882*cos（2*seita））/60;
????????%真太陽時/小時
????????S0（ji）=SD（ji）+ET;
????????%太陽時角/弧度
????????TAO（ji）=pi*（S0（ji）-12）*15/180;
????????%高度角/弧度
????????H0（ji）=asin（sin（ED）*sin（FA）+cos（ED）*cos（FA）*cos（TAO（ji）））;
????????%方位角??弧度
????????A（ji）=acos（（sin（H0（ji））*sin（FA）-sin（ED））/（cos（H0（ji））*cos（FA）））;

????????%高度角弧度轉為角度
????????H1（ji）=（180*H0（ji））/pi;
????????%方位角弧度轉為角度
????