%真近點角的計算Fs：
%Ms：平近點角    Es：偏近點角
DToR=2*pi/360;  %轉換常數從角度到弧度
Ms=0:1:359; %平近點序列（角度）
es=0.01;    %地球偏心率
Ms=DToR*Ms;
%用迭代法求偏近點角
%Es=Ms+es*sin(Es);  開普勒方程，迭代計算公式
Es=zeros(1,360);    %生成偏近點角序列的存儲變量
for i=1:360
    dEs=1;          %賦初值
    TempMs=Ms(i);
    TempEs=TempMs;
    while (abs(dEs)>0.0001)     %計算精度控制循環的循環次數
        TempEsNext=TempMs+es*sin(TempEs);
        dEs=TempEs-TempEsNext;
        TempEs=TempEsNext;
    end
    Es(i)=TempEs;   %存儲計算所得的偏近點角
end
%Es=Es/DTor;
%Ms=Ms/DTor;
Fs=zeros(1,360);    %生成真近點角序列的存儲變量，已知偏近點角求真近點，
v1=sqrt((1+es)/(1-es)); %公式中使用的常數，注意避免將此計算過程放入循環！
for i=1:360
    v2=v1*tan(Es(i)/2);
    Fs(i)=mod((2*atan(v2)+2*pi),2*pi);
end
Es=Es/DToR; %將平近點轉換為角度
Ms=Ms/DToR; %將偏近點轉換為角度
Fs=Fs/DToR; %將真近點轉換為角度