lanbote.zip_轨道根数

%大作业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)