matlab视日轨迹跟踪算法仿真程序

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;
        %方位角弧度转为角度