%大作业2---兰伯特问题
%问题描述：已知t1,t2时刻的位置r1,r2，求轨道根数
%输入参数：t1,t2,r1,r2,a的精度jd,u
%输出参数：轨道根数 a,e,i,ou,w,f(或E)
%a的第一次估计值a0取为r1,r2的均值
function [a,e,i,ou,w,M0]=lanbote(t1,t2,R1,R2,jd,u)
% r1=sqrt((R1(1))^2+(R1(2))^2+(R1(3))^2);
% r2=sqrt((R2(1))^2+(R2(2))^2+(R2(3))^2);
%第9页PPT第三题---先求轨道根数，再求速度
% clear;clc;
% t1=0;
% t2=1200;
% R1=[5644 -2830 -4170];
% R2=[-2240 7320 -4980];
% jd=1;
% u=398600;
r1=norm(R1);
r2=norm(R2);
%求弦长
R3=R1-R2;
% r3=sqrt((R3(1))^2+(R3(2))^2+(R3(3))^2);
r3=norm(R3);
t3=abs(t2-t1);
%迭代初值a0
a0=(r1+r2)/2;
%迭代求a
aa(1)=a0;
aa(2)=diedai(aa(1),a0,r3,t3,u);
j=1;%注意i已经用过了
while (abs(aa(j+1)-aa(j))>jd && j<5)
j=j+1;
aa(j+1)=diedai(aa(j),a0,r3,t3,u);
end
a=aa(j);
n=sqrt(u/(a^3));
%以上求得a,下面求出ep,dt，然后依次求出轨道根数
ep=2*asind(sqrt((2*a0+r3)/4/a));
dt=2*asind(sqrt((2*a0-r3)/4/a));
%根据（2-178）（2-179）求e 和 E1，然后根据开普勒方程求M
k1=((1-r1/a)*cosd(ep-dt)-(1-r2/a))/(sind(ep-dt));
k2=1-r1/a;
e=sqrt((k1)^2+(k1)^2);
E1=k2/e;
E2=ep-dt+E1;
M0=E1-e*sind(E1)-n*t1;
%下面求i,w,ou
kk1=a*(cosd(E1)-e);
ss1=a*sqrt(1-e*e)*sind(E1);
kk2=a*(cosd(E2)-e);
ss2=a*sqrt(1-e*e)*sind(E2);
% [i,ou]=solve('kk1*sind(ou)*sind(i)+ss1*cosd(ou)*sind(i)-R1(3)','kk2*sind(ou)*sind(i)+ss2*cosd(ou)*sind(i)-R2(3)');
l1=(R1(3)*kk2-kk1*R2(3))/(ss1*kk2-ss2*kk1);
l2=(R1(3)*ss2-R2(3)*ss1)/(kk1*ss2-kk2*ss1);
i=asind((l1)^2+(l2)^2);
w=atand(l2/l1);
%求ou
ll1=kk1*cosd(w)-ss1*sind(w);
ll2=kk1*sind(w)*cosd(i)+ss1*cosd(w)*cosd(i);
ll3=kk2*cosd(w)-ss2*sind(w);
ll4=kk2*sind(w)*cosd(i)+ss2*cosd(w)*cosd(i);
ou=acosd((R1(1)*ll4-R2(1)*ll2)/(ll1*ll4-ll2*ll3));
jieguo=[a,e,i,ou,w,M0];
function y=diedai(x,a0,r3,t3,u)
ep=2*asind(sqrt((2*a0+r3)/4/x));
dt=2*asind(sqrt((2*a0-r3)/4/x));
n=(ep-dt-sind(ep)+sind(dt))/t3;
y=power(u/n/n,1/3);
%nthroot(x,3)或power(x,n)
评论0