poavel.m

这是用于计算天文中的开普勒轨道元素的程序

function PoAVel
%坐标，速度，旋转矩阵的表达式
syms a e I o O E m
polycoord=[a*cos(E)-a*e;sin(E)*a*sqrt(1-e^2);0];
polyvel=(sqrt(m/a)/(1-e*cos(E))^2)*[(e*sin(E))*(cos(o)*cos(E)-e*cos(o)-sin(o)*sin(E)*sqrt(1-e^2))-sqrt(1-e^2)*(sin(o)*cos(E)-e*sin(o)+(1-e^2)*cos(o)*sin(E));(e*sin(E))*(cos(o)*cos(E)-e*cos(o)-sin(o)*sin(E)*sqrt(1-e^2))+sqrt(1-e^2)*(sin(o)*cos(E)-e*sin(o)+(1-e^2)*cos(o)*sin(E));0];
Ro3OMat=[cos(O) sin(O) 0;-sin(O) cos(O) 0;0 0 1];
Ro1iMat=[1 0 0;0 cos(I) sin(I);0 -sin(I) cos(I)];
Ro3oMat=[cos(o) sin(o) 0;-sin(o) cos(o) 0;0 0 1];
RoMat=Ro3oMat*Ro1iMat*Ro3OMat;
%天球赤道坐标系中的坐标，速度表达式
RoMatcoord=RoMat*polycoord;
RoMatvel=RoMat*polyvel;
%读入数据
KepEle=input('输入轨道元素:[a,e,i,Ω,ω,M,μ]');
ValE=EcNeAngle(KepEle(2),degree2rad(KepEle(6)));
%代入坐标，速度，旋转矩阵，及天球坐标系中的坐标，速度的表达式
Valcoord=vpa(subs(polycoord,{a,e,E},{KepEle(1),KepEle(2),ValE}),7);
Valvel=vpa(subs(polyvel,{a,e,E,m},{KepEle(1),KepEle(2),ValE,KepEle(7)}),7);
ValRoMat=vpa(subs(RoMat,{I,o,O},{degree2rad(KepEle(3)),degree2rad(KepEle(5)),degree2rad(KepEle(4))}),7);
%vO=subs(Ro3OMat,O,degree2rad(KepEle(4)))
%vo=subs(Ro3oMat,o,degree2rad(KepEle(5)))
%vi=subs(Ro1iMat,I,degree2rad(KepEle(3)))
coordinates=subs(RoMatcoord,{a,e,I,o,O,E},{KepEle(1),KepEle(2),degree2rad(KepEle(3)),degree2rad(KepEle(5)),degree2rad(KepEle(4)),ValE});
velocity=subs(RoMatvel,{a,e,I,o,O,E,m},{KepEle(1),KepEle(2),degree2rad(KepEle(3)),degree2rad(KepEle(5)),degree2rad(KepEle(4)),ValE,KepEle(7)});
%将天球赤道坐标系中的坐标，速度转化为数值
Coordinates=vpa(coordinates,7)
Velocity=vpa(velocity,7)