已知轨道根数：半长轴、偏心率、倾角、近地点幅角、近地点时刻、生交点赤经，求解目标（卫星）的位置矢量
clc clear all close all miu = 3.9860047e14; %地球引力常数(m^3/s^2) M = @(t,tao,a)(t-tao)*(sqrt(miu/a^3)); v = @(M,a,e)M+e*(2-e^2/4+5*e^4/96)*sin(M)+e^2*(5/4-11*e^2/24)*sin(2*M) ... +e^3*(13/12-43*e^2/64)*sin(3*M)+103*e^4*sin(4*M)/96+1097*e^5*sin(5*M)/960; r = @(v,a,e)a*(1-e.^2)./(1+e.*cos(v)); %三轴转换矩阵 Rx = @(i)[1 0 0 0 cos(i) sin(i) 0 -sin(i) cos(i)]; %轨道倾角 Ry = @(theta)[cos(theta) 0 sin(theta) 0 1 0 -sin(theta) 0 cos(theta)]; Rz = @(w)[cos(w) sin(w) 0 -sin(w) cos(w) -2 0 0 1]; %升交点赤经或者近地点幅角 % p0 = [r*cos(v),r*sin(v),0]; % p = Rz(-Omiga)*Rx(-i)*Rz(-w)*p0; %各量在J2000惯性坐标系下的矢量坐标 %% 目标卫星轨道 a_o = 6862.8 ; %轨道半长轴 (km) e_o = 0.001884; %偏心率 i_o = 98.79*pi/180 ; %轨道倾角 (度) w_o = 60*pi/180 ; %近地点幅角（度） Omiga_o = 30*pi/180; %升交点赤经（度） tao_o = 2.452685925e6 ; %过近地点时刻：2003-2-5-10.12.00 %% 距离解算 Position = []; Velocity = []; count = 0; for t = 2456961:0.01:2456961.5 %探测时刻从20014年10月30日12：00至24:00，为儒略日计时 count = count+1; M_o = M(t,tao_o,a_o); v_o = v(M_o,a_o,e_o); r_o = r(v_o,a_o,e_o); p0_o = [r_o*cos(v_o);r_o*sin(v_o);0]; %位置矢量 p_o = Rz(-Omiga_o)*Rx(-i_o)*Rz(-w_o)*p0_o; %目标卫星在J2000惯性坐标系下的矢量坐标 Position(count) = p_o; %写入数据存入向量 c = sqrt(miu*a_o/(1-e^2)); V0 = [-c*sin(v_o);c*(e.o + cos(v_o)); 0]; %速度 V_o = Rz(-Omiga_o)*Rx(-i_o)*Rz(-w_o)*V0; %目标卫星在J2000惯性坐标系下的速度矢量 Velocity (count) = V_o; end有了位置矢量可以对轨道进行轨迹描绘
plot3(Position(1,:),Position(2,:),Position(3,:))
儒略日的求解
Year=2014; Month=10; Day=30; Hour=12; Min=0; Sec=0; if Month == 1 || Month ==2 f=Year-1; g=Month+12; end if Month>=3 f=Year; g=Month; end mid1=floor(365.25*f); mid2=floor(30.6001*(g+1)); A=2-floor(f/100)+floor(f/400); J=mid1+mid2+Day+A+1720994.5; disp('儒略日：') JDE = double(J+Hour/24+Min/1440+Sec/86400)2000年1月1日12：00的儒略日为2451545
参考文献：汪洪源, 陈赟. 天基空间目标红外动态辐射特性建模与仿真[J]. 红外与激光工程, 2016(5):11-17.
