yevent是空的——它没有元素。调试程序找出原因。
我不知道为什么我得到错误毕竟努力工作
1视图(30天)
显示旧的评论
% patched-conic重力帮助星际
%轨迹设计和优化
%喷气推进实验室星历表和64位SNOPT算法
%与MATLAB轨道力学
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
清晰;
全球鸸鹋smu xmu aunit ip1 ip2 ip3 jdtdb0 jdtdbi1
全球iephem公里ephname eq2000 pmu代表rsoi fbalt
全球dv_launch dvm_launch dv_arrival dvm_arrival vinfm_in vinfm_out
全球fbrp c3launch jdtdb1 jdtdb2
全球otype rito1 vito1 rito2 vito2逆向
全球vinf_in vinf_out jdtdb_soi rp2sc vp2sc
% J2000 ecliptic-to-equatorial变换矩阵
eq2000 = [(1.000000000000000 0 0);…
[0 0.917482062069182 - -0.397777155931914);…
0.397777155931914 - 0.917482062069182 [0]];
%定义天文单位(公里)
aunit = 149597870.691;
%弧度度转换因子
rtd = 180.0 /π;
%定义“参考”tdb朱利安日期
jdtdb0 =朱利安(1,1,2000);
%初始化喷气推进实验室星历表
ephname =“de421.bin”;
iephem = 1;
公里= 1;
%定义行星的名字向量
pdata = [“水星”;“金星”;“地球”;“火星”;…
“木星”;“土星”;“天王星”;“海王星”;“冥王星”];
pname = cellstr (pdata);
%的地球引力常量(公里/秒^ ^ 3 2)
pmu (1) = 22032.08;
pmu (2) = 324858.599;
pmu (3) = 398600.436;
pmu (4) = 42828.314;
pmu (5) = 126712767.863;
pmu (6) = 37940626.063;
pmu (7) = 5794549.007;
pmu (8) = 6836534.064;
pmu (9) = 981.601;
xmu = 0.295912208285591149 e 03;
smu = xmu * aunit ^ 3/86400.0 ^ 2;
xmu = 0.899701152970881167 e-09;
鸸鹋= xmu * aunit ^ 3/86400.0 ^ 2;
%的地球赤道半径(公里)
代表(1)= 2439.7;
代表(2)= 6051.9;
代表(3)= 6378.14;
代表(4)= 3397.0;
代表(5)= 71492.0;
代表(6)= 60268.0;
代表(7)= 25559.0;
代表(8)= 24764.0;
代表(9)= 1151.0;
%的势力范围为每个行星(公里)
rsoi (1) = 112409.906759936;
rsoi (2) = 616277.129297850;
rsoi (3) = 924647.107795632;
rsoi (4) = 577231.618115568;
rsoi (5) = 48204698.6886979;
rsoi (6) = 54553723.6086354;
rsoi (7) = 51771106.3741412;
rsoi (8) = 86696170.2674129;
rsoi (9) = 15110628.1503479;
%开始模拟
clc;回家;
流(“\ nprogram flyby_matlab \ n”);
流(“\ npatched-conic重力辅助轨迹分析\ n \ n”);
%请求离开tdb日历日期
流(“\ ndeparture日历猜\ n”);
(' 1 < < = 12日1 < = =月天< = 31年=所有数字!”);
3;5;2030;
% (\ n \ 03、05年2030年)=获取当前日期;
%计算离职tdb朱利安日期
jdtdbi1 =朱利安(5 5 2030);
而(1)
流(‘\ nplease输入离职日期搜索边界天\ n”);
ddays1 =输入(“30”);
如果(ddays1 > = 0.0)
打破;
结束
结束
%请求飞越tdb日历日期
流(' \ n \ nflyby——日历猜\ n”);
(' 1 < < = 12日1 < = =月天< = 31年=所有数字!”);
3;5;2037;
%(月、日、年)=获取当前日期;
%计算飞越tdb朱利安日期
jdtdbi2 =朱利安(2037年03、05);
% jdtdbi2 =朱利安(月、日、年);
而(1)
流(输入“\ nplease飞越日期搜索边界天\ n”);
ddays2 =输入(“?”);
如果(ddays2 > = 0.0)
打破;
结束
结束
%请求到来tdb日历日期
流(' \ n \ narrival——日历猜\ n”);
(' 1 < < = 12日1 < = =月天< = 31年=所有数字!”);
8;12;2050;
%(月;一天;年)=获取当前日期;
%计算到达tdb朱利安日期
jdtdbi3 =朱利安(8,2050);
% jdtdbi3 =朱利安(月、日、年);
而(1)
流(‘\ nplease输入到达日期搜索边界天\ n”);
ddays3 =输入(“?”);
如果(ddays3 > = 0.0)
打破;
结束
结束
%请求离开,飞越和到达行星
为i = 1:1:3
流(“\ n \ n星球菜单”);
流(“\ n < 1 >汞”);
流(“\ n < 2 >金星”);
流(“\ n < 3 >地球”);
流(' \ n < 4 >火星的);
流(“\ n < 5 >木星”);
流(“土星\ n < 6 >”);
流(“\ n < 7 >天王星”);
流(“\ n < 8 >海王星”);
流(“\ n < 9 >冥王星”);
如果(i = = 1)
而(1)
流(' \ n \ nplease选择离开地球\ n”);
ip1 =输入(“?”);
如果(ip1 > = 1 & & ip1 < = 9)
打破;
结束
结束
结束
如果(我= = 2)
而(1)
流(' \ n \ nplease选择飞越地球\ n”);
ip2 =输入(“?”);
如果(ip2 > = 1 & & ip2 < = 9)
打破;
结束
结束
结束
如果(我= = 3)
而(1)
流(' \ n \ nplease选择到达行星\ n”);
ip3 =输入(“?”);
如果(ip3 > = 1 & & ip3 < = 9)
打破;
结束
结束
结束
结束
%的要求低,上界的飞越高度(公里)
而(1)
流(' \ n \ nplease输入的下界飞越高度(公里)\ n ');
fbalt_lwr =输入(“?”);
如果(fbalt_lwr > = 0.0)
打破;
结束
结束
而(1)
流(' \ n \ nplease输入飞越高度的上限(公里)\ n ');
fbalt_upr =输入(“?”);
如果(fbalt_upr > = fbalt_lwr)
打破;
结束
结束
%的请求类型的优化
而(1)
流(“\ n \ n \ n优化菜单”);
流(' \ n < 1 >最小化出发δv \ n ');
流(' \ n < 2 >最小化到来δv \ n”);
流(“总δv \ n \ n < 3 >最小化”);
流(' \ n选择(1、2或3)\ n ');
otype =输入(“?”);
如果(otype = = 1 | | otype = = 2 | | otype = = 3)
打破;
结束
结束
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% %解决patched-conic重力辅助问题
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
clc;回家;
xg = 0 (3,1);
%控制变量初始猜测
%(离职,飞越,到达tdb朱利安日期)
xg (1) = jdtdbi1 - jdtdb0;
xg (2) = jdtdbi2 - jdtdb0;
xg (3) = jdtdbi3 - jdtdb0;
%控制变量界限
xlwr (1) = xg (1) - ddays1;
xupr (1) = xg (1) + ddays1;
xlwr (2) = xg (2) - ddays2;
xupr (2) = xg (2) + ddays2;
xlwr (3) = xg (3) - ddays3;
xupr (3) = xg (3) + ddays3;
xlwr = xlwr ';
xupr = xupr ';
在目标函数%范围
流(1)=无穷;
fupp(1) = +正;
%飞越界限高度
流(2)= fbalt_lwr;
fupp (2) = fbalt_upr;
%界限v-infinity匹配(平等)约束
流(3)= 0.0;
fupp (3) = 0.0;
流=流”;
fupp = fupp ';
xmul = 0 (3,1);
xstate = 0 (3,1);
fmul = 0 (3,1);
fstate = 0 (3,1);
%使用snopt找到最优解
f (x),通知、xmul fmul] = snopt (xg、xlwr xupr, xmul, xstate,…
流、fupp fmul fstate,“fbyfunc”);
%计算最优解
(~ g) = fbyfunc (x);
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
%计算出发的双曲线定位
% (J2000地球赤道和equinox)
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
dveq1 = eq2000 * dv_launch;
dveqm1 =规范(dveq1);
decl1 = 90.0 - rtd *这些“可信赖医疗组织”(dveq1 (3) / dveqm1);
rasc1 = rtd * atan3 (dveq1 (2), dveq1 (1));
%计算tdb朱利安日期最优转移
jdtdb1 = jdtdb0 + x (1);
jdtdb2 = jdtdb0 + x (2);
jdtdb3 = jdtdb0 + x (3);
%转换解决方案朱利安日期日历日期和tdb倍
[cdstr1, tdbstr1] = jd2str (jdtdb1);
[cdstr2, tdbstr2] = jd2str (jdtdb2);
[cdstr3, tdbstr3] = jd2str (jdtdb3);
% % % % % % % % % % % % % % % % %
% %打印结果
% % % % % % % % % % % % % % % % %
流(“\ nprogram flyby_matlab”);
流(' \ n - - - - - - - - - - - - - - - - - - - - - - \ n”);
流(' \ n \ npatched-conic重力辅助轨迹分析);
如果(otype = = 1)
流(“\ nminimize离开δv \ n”);
结束
如果(otype = = 2)
流(“δv \ n \ nminimize到来”);
结束
如果(otype = = 3)
流(“总δv \ n \ nminimize”);
结束
流(“\ ndeparture星球”);
disp (pname (ip1));
流(“飞越地球”);
disp (pname (ip2));
流(到达地球的);
disp (pname (ip3));
%离职条件
流(“\ nLAUNCH条件”);
流(' \ n = = = = = = = = = = = = = = = = = \ n ');
流(“\ ndeparture日程表日期”);
disp (cdstr1);
流(“\ ndeparture TDB时间”);
disp (tdbstr1);
流(“\ % 12.6 f ndeparture TDB朱利安日期”,jdtdb1);
流(' \ n \ nheliocentric发射δv向量和大小的);
流(' \ n(地球黄道平和J2000 equinox)”);
流(' \ n - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - \ n”);
流(“\ nlaunch delta-vx f % 12.6米/秒的,1000.0 * dv_launch (1));
流(“\ nlaunch delta-vy f % 12.6米/秒的,1000.0 * dv_launch (2));
流(“\ nlaunch delta-vz f % 12.6米/秒的,1000.0 * dv_launch (3));
流(' \ n \ nlaunchδv % 12.6 f m / s \ n '1000.0 * dvm_launch);
流(“\ nlaunch双曲线特征”);
流(' \ n(地球赤道的意思和J2000 equinox)”);
流(' \ n - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - \ n”);
流(“\ nasymptote对提升% 12.6 f度\ n ',rasc1);
流(“\ nasymptote赤纬% 12.6 f度\ n ',decl1);
流(' \ nlaunch能源% 12.6 f公里^ 2 /秒^ 2 \ n”,c3launch);
流(“\ npost-impulse日心轨道航天器的元素和状态向量);
流(' \ n(地球黄道平和J2000 equinox)”);
流(' \ n - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - \ n”);
oev1 = eci2orb1 (smu rito1 vito1);
oeprint1 (smu oev1 3);
svprint (rito1 vito1);
流(“日心轨道参数和状态向量离开地球的);
流(' \ n(地球黄道平和J2000 equinox)”);
流(' \ n - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - \ n”);
[一国,vp1] = p2000 (ip1 jdtdb1);
oev1 = eci2orb1 (vp1一国,smu);
oeprint1 (smu oev1 3);
svprint(一国,vp1);
%飞越条件
流(“\ nPATCHED-CONIC飞越条件”);
流(' \ n = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = \ n ');
流(“\ nflyby日程表日期”);
disp (cdstr2);
流(“\ nflyby TDB时间”);
disp (tdbstr2);
流(“\ nflyby TDB朱利安日期% 12.6 f \ n ',jdtdb2);
流(' \ n \ f nlaunch-to-flyby时间% 12.6天的,jdtdb2 jdtdb1);
流(“\ nv-infinity f % 12.6米/秒的1000.0 * vinfm_in);
流(“\ nv-infinity % 12.6 f m / s \ n '1000.0 * vinfm_out);
流(“\ nflyby高度% 12.6 f公里\ n ',fbalt);
%将角度
tmp =代表(ip2) * vinfm_in * vinfm_in / pmu (ip2);
tangle1 = 2.0 *印度历的7月(1.0 / (1.0 + tmp));
流(“\ nmaximum把% 12.6 f度”tangle1 * rtd);
tmp = (fbalt +代表(ip2)) * vinfm_in * vinfm_in / pmu (ip2);
tangle2 = 2.0 *印度历的7月(1.0 / (1.0 + tmp));
流(“\ nactual转角% 12.6 f度\ n 'tangle2 * rtd);
流(“\ nheliocentricδv % 12.6 f m / s \ n '1000.0 *逆向);
%日心均
dvh_max =√pmu (ip2) /代表(ip2));
流(“\ nmax日心δv % 12.6 f m / s \ n '1000.0 * dvh_max);
在近拱点% planet-centered轨道要素
oev_periapsis = fbhyper (pmu (ip2)、vinf_in vinf_out, fbrp);
[rp2sc1, vp2sc1] = orb2eci (pmu (ip2) oev_periapsis);
流(“\ nplanet-centered轨道航天器在近拱点的元素和状态向量);
流(' \ n(地球黄道平和J2000 equinox)”);
流(' \ n - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - \ n”);
oev_periapsis oeprint1 (pmu (ip2), 1);
svprint (rp2sc1 vp2sc1);
[bplane, ~, ~, ~] = rv2bplane (pmu (ip2)、rp2sc1 vp2sc1);
流(“b-plane航天器在近拱点的坐标);
流(' \ n(地球黄道平和J2000 equinox)”);
流(' \ n - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - \ n”);
bpprint (rp2sc1 vp2sc1 bplane);
cdecl_asy = cos (bplane (2));
crasc_asy = cos (bplane (3));
srasc_asy =罪(bplane (3));
sdecl_asy =罪(bplane (2));
流(“\ nheliocentric航天器的轨道参数和状态向量);
流(' \ n(地球黄道平和J2000 equinox)”);
流(' \ n - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - \ n”);
τ= 86400.0 * (jdtdb2 - jdtdb1);
副总裁(rp) = twobody2 (smu,τ,rito1 vito1);
项目= eci2orb1(副总裁smu, rp);
oeprint1 (smu项目3);
副总裁svprint (rp);
流(“\ nheliocentric飞越行星的轨道参数和状态向量);
流(' \ n(地球黄道平和J2000 equinox)”);
流(' \ n - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - \ n”);
副总裁(rp) = p2000 (ip2 jdtdb2);
项目= eci2orb1(副总裁smu, rp);
oeprint1 (smu项目3);
副总裁svprint (rp);
%到达条件
流(“\ nARRIVAL条件”);
流(' \ n = = = = = = = = = = = = = = = = = = \ n ');
流(“\ narrival日程表日期”);
disp (cdstr3);
流(“\ narrival TDB时间”);
disp (tdbstr3);
流(“\ narrival TDB朱利安日期% 12.6 f \ n ',jdtdb3);
流(' \ n \ f nflyby-to-arrival时间% 12.6天的,jdtdb3 jdtdb2);
流(“\ nheliocentric到来δv向量和级”);
流(' \ n(地球黄道平和J2000 equinox)”);
流(' \ n - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - \ n”);
流(“\ narrival delta-vx f % 12.6米/秒的,1000.0 * dv_arrival (1));
流(“\ narrival delta-vy f % 12.6米/秒的,1000.0 * dv_arrival (2));
流(“\ narrival delta-vz f % 12.6米/秒的,1000.0 * dv_arrival (3));
流(' \ n \ narrivalδv % 12.6 f m / s \ n '1000.0 * dvm_arrival);
%双体传播第二回合的转移轨道
[rf3, vf3] = twobody2 (smu 86400 * (jdtdb3 - jdtdb2), rito2, vito2);
流(“\ npre-impulse日心轨道航天器的元素和状态向量);
流(' \ n(地球黄道平和J2000 equinox)”);
流(' \ n - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - \ n”);
oev3 = eci2orb1 (smu rf3 vf3);
oeprint1 (smu oev3 3);
svprint (rf3 vf3);
流(“\ npost-impulse日心轨道航天器的元素和状态向量);
流(' \ n(地球黄道平和J2000 equinox)”);
流(' \ n - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - \ n”);
rf4 = rf3;
vf4 (1) = vf3 (1) + dv_arrival (1);
vf4 (2) = vf3 (2) + dv_arrival (2);
vf4 (3) = vf3 (3) + dv_arrival (3);
oev4 = eci2orb1 (smu rf4 vf4);
oeprint1 (smu oev4 3);
svprint (rf4 vf4);
流(“日心轨道参数和状态向量到达地球的);
流(' \ n(地球黄道平和J2000 equinox)”);
流(' \ n - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - \ n”);
副总裁(rp) = p2000 (ip3 jdtdb3);
项目= eci2orb1(副总裁smu, rp);
oeprint1 (smu项目3);
副总裁svprint (rp);
流(“\ nMISSION总结”);
流(' \ n = = = = = = = = = = = = = = = \ n ');
流(“\ ntotalδv % 12.6 f m / s \ n ',…
1000.0 * (dvm_launch + dvm_arrival));
流(' \ ntotal能源% 12.6 f公里^ 2 /秒^ 2 \ n”,dvm_launch ^ 2 + dvm_arrival ^ 2);
流(‘\ ntotal任务持续时间% 12.6 f天\ n \ n ',jdtdb3 jdtdb1);
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
%双体集成的轨迹
%的第一个冲动在飞越行星SOI入口
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
%设置数值
选择= odeset (“RelTol”1.0 e-6“AbsTol”1.0 e-8“事件”,@soi_event);
%定义最大搜索时间(秒)
tof = 86400.0 * (jdtdb2 - jdtdb1);
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
%求飞越地球SOI入口条件
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
[tevent ~, ~, ~, ~] =数值(@twobody_heqm [0 tof], [rito1 vito1],选项);
jdtdb_soi = jdtdb1 + tevent / 86400.0;
[cdstr_ca, utstr_ca] = jd2str (jdtdb_soi);
流(“\ nPATCHED-CONIC SOI入口条件”);
流(' \ n = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = \ n ');
流(“\ ncalendar日期”);
disp (cdstr_ca);
流(“\ nTDB时间”);
disp (utstr_ca);
流(' % \ nTDB朱利安日期12.6 f \ n”,jdtdb_soi);
流(“\ nplanet-centered航天器的轨道参数和状态向量);
流(' \ n(地球黄道平和J2000 equinox)”);
流(' \ n - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - \ n”);
oev_soi = eci2orb1 (pmu (ip2)、rp2sc vp2sc);
oev_soi oeprint1 (pmu (ip2), 1);
svprint (rp2sc vp2sc);
[bplane, ~, ~, ~] = rv2bplane (pmu (ip2)、rp2sc vp2sc);
流(“b-plane坐标在势力范围”);
流(' \ n(地球黄道平和J2000 equinox)”);
流(' \ n - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - \ n”);
bpprint (rp2sc vp2sc bplane);
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
%集成轨迹从SOI入口
%向飞越最接近地球
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
%设置数值
选择= odeset (“RelTol”1.0平台以及“AbsTol”1.0平台以及“事件”,@fpa_event);
%定义最大搜索时间(秒)
tof = 10.0 * 86400.0;
[~,~,tevent yevent() =数值(@peqm [0 tof], [rp2sc vp2sc],选项);
在最接近%提取状态向量
rsc_ca = yevent (1:3);
vsc_ca = yevent (6);
在最接近% B-plane飞越双曲线的坐标
[bplane、房车、电视,ibperr] = rv2bplane (pmu (ip2)、rsc_ca vsc_ca);
在最接近% tdb朱利安日期
jdtdb_ca = jdtdb_soi + tevent / 86400.0;
[cdstr_ca, utstr_ca] = jd2str (jdtdb_ca);
流(“\ n \ nNUMERICALLY集成最接近条件”);
流(' \ n = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = \ n ');
流(“\ ncalendar日期”);
disp (cdstr_ca);
流(“\ nTDB时间”);
disp (utstr_ca);
流(' % \ nTDB朱利安日期12.6 f \ n”,jdtdb_ca);
流(“\ nplanet-centered航天器的轨道参数和状态向量);
流(' \ n(地球黄道平和J2000 equinox)”);
流(' \ n - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - \ n”);
oev_ca = eci2orb1 (pmu (ip2)、rsc_ca vsc_ca);
oev_ca oeprint1 (pmu (ip2), 1);
svprint (rsc_ca vsc_ca);
流(b-plane坐标在最接近的);
流(' \ n(地球黄道平和J2000 equinox)”);
流(' \ n - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - \ n”);
bpprint (rsc_ca vsc_ca bplane);
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
%日心和飞越图形%
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
而(1)
流(' \ n \ nwould你想为这个任务创建图形(y =是的,n = no) \ n”);
slct =输入(“?”,“年代”);
如果(slct = =“y”| | slct = =“n”)
打破;
结束
结束
如果(slct = =“y”)
而(1)
流(' \ n \ nplease输入情节步长(天)\ n ');
deltat =输入(“?”);
如果(deltat > 0)
打破;
结束
结束
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
%情节日心行星轨道和轨道转移
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
图(1);
%离开行星半轴长和周期
(r1, v1) = p2000 (ip1 jdtdb1);
oev1 = eci2orb1 (smu, r1, v1);
period1 = 2.0 *π* oev1 (1) * sqrt (oev1 (1) / smu) / 86400.0;
%飞越行星半轴长和周期
(r2, v2) = p2000 (ip2 jdtdb1);
oev2 = eci2orb1 (smu, r2, v2);
period2 = 2.0 *π* oev2 (1) * sqrt (oev2 (1) / smu) / 86400.0;
%到达行星半轴长和周期
(r3, v3) = p2000 (ip3 jdtdb1);
oev3 = eci2orb1 (smu, r3, v3);
period3 = 2.0 *π* oev3 (1) * sqrt (oev3 (1) / smu) / 86400.0;
%的“行星”情节
npts1 =修复(period1 / deltat);
npts2 =修复(period2 / deltat);
npts3 =修复(period3 / deltat);
%离开地球轨道
为我= 0:1:npts1
jdtdb = jdtdb1 + i * deltat;
(r1 ~) = p2000 (ip1 jdtdb);
x1 (i + 1) = r1 (1) / aunit;
日元(i + 1) = r1 (2) / aunit;
结束
%计算数据点
(r1, v1) = p2000 (ip1 jdtdb1 + period1);
x1 (npts1 + 1) = r1 (1) / aunit;
日元(npts1 + 1) = r1 (2) / aunit;
%飞越地球日心轨道
为我= 0:1:npts2
jdtdb = jdtdb1 + i * deltat;
(r2, ~) = p2000 (ip2 jdtdb);
r2 x2 (i + 1) = (1) / aunit;
y2 r2 (i + 1) = (2) / aunit;
结束
%计算数据点
(r2, v2) = p2000 (ip2 jdtdb1 + period2);
r2 x2 (npts2 + 1) = (1) / aunit;
y2 r2 (npts2 + 1) = (2) / aunit;
%到达地球日心轨道
为我= 0:1:npts3
jdtdb = jdtdb1 + i * deltat;
(r3, ~) = p2000 (ip3 jdtdb);
r3 x3 (i + 1) = (1) / aunit;
r3 y3 (i + 1) = (2) / aunit;
结束
%计算数据点
(r3, v3) = p2000 (ip3 jdtdb1 + period3);
r3 x3 (npts3 + 1) = (1) / aunit;
r3 y3 (npts3 + 1) = (2) / aunit;
%转移轨道的第一站
npts4 =修复((jdtdb2 - jdtdb1) / deltat);
为我= 0:1:npts4
τ= 86400.0 *我* deltat;
[rft1 ~] = twobody2 (smu,τ,rito1 vito1);
x4 (i + 1) = rft1 (1) / aunit;
y4 (i + 1) = rft1 (2) / aunit;
结束
%计算首回合的最后数据点
τ= 86400.0 * (jdtdb2 - jdtdb1);
[rft1, vft1] = twobody2 (smu,τ,rito1 vito1);
x4 (npts4 + 1) = rft1 (1) / aunit;
y4 (npts4 + 1) = rft1 (2) / aunit;
%日心转移轨道的第二站
npts5 =修复((jdtdb3 - jdtdb2) / deltat);
为我= 0:1:npts5
τ= 86400.0 *我* deltat;
[rft2, vft2] = twobody2 (smu,τ,rito2 vito2);
x5 (i + 1) = rft2 (1) / aunit;
日元(i + 1) = rft2 (2) / aunit;
结束
%计算第二回合的最后数据点
τ= 86400.0 * (jdtdb3 - jdtdb2);
[rfto2, vfto2] = twobody2 (smu,τ,rito2 vito2);
x5 (npts5 + 1) = rfto2 (1) / aunit;
日元(npts5 + 1) = rfto2 (2) / aunit;
%确定春分“大小”
xve = oev1 (1) / aunit;
如果(oev2 (1) > oev1 (1))
xve = oev2 (1) / aunit;
结束
如果(oev3 (1) > oev2 (1))
xve = oev3 (1) / aunit;
结束
%情节日心行星轨道和patched-conic转移轨迹
持有在;
情节(x1, y1,“。b”);
情节(x1, y1,“- b”);
情节(x2, y2,“.g”);
情节(x2, y2,“g”);
情节(x3, y3,“r”);
情节(x3, y3,“- r”);
情节(x4、y4“同意”);
情节(x4、y4“- k”);
情节(x5,日元,“同意”);
情节(x5,日元,“- k”);
情节(0,0,“衔接”,“MarkerSize”10);
%的阴谋和标签春分方向
线([0.05,xve], (0, 0),“颜色”,“黑”);
文本(1.05 * xve 0“\ Upsilon”);
%标签发射,飞越和到达地点
情节(x4 (1)、y4 (1),‘* k”);
文本(x4 (1) + 0.06, y4 (1) - 0.01,' D ');
情节(x5(1),日元(1),‘* k”);
文本(x5(1) + 0.06,日元(1)+ 0.01,“F”);
情节(x5 (npts5 + 1),日元(npts5 + 1),‘* k”);
文本(x5 (npts5 + 1) + 0.06,日元(npts5 + 1) - 0.01,“一个”);
轴平等的;
%显示发射,飞越和到达日期
文本(0.85 * xve -xve + 0.8,“离开”,“字形大小”8);
文本(0.875 * xve -xve + 0.7, pname (ip1),“字形大小”8);
文本(0.875 * xve -xve + 0.6, cdstr1,“字形大小”8);
文本(0.85 * xve -xve + 0.4,“到达”,“字形大小”8);
文本(0.875 * xve -xve + 0.3, pname (ip3),“字形大小”8);
文本(0.875 * xve -xve + 0.2, cdstr3,“字形大小”8);
文本(0.85 * xve xve - 0.2,“飞越”,“字形大小”8);
文本(0.875 * xve xve - 0.3, pname (ip2),“字形大小”8);
文本(0.875 * xve xve - 0.4, cdstr2,“字形大小”8);
文本(0.875 * xve xve - 0.5, tdbstr2,“字形大小”8);
%标签图和轴在天文单位(AU)
包含(“X坐标(AU)”,“字形大小”12);
ylabel (“Y坐标(AU)”,“字形大小”12);
标题(“Patched-conic重力辅助轨迹”,“字形大小”16);
变焦在;
%下一行创建一个颜色用tiff预览eps图形文件
打印-depsc tiff r300 flyby_matlab1.eps
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
%创建三维图形的飞越
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
%单位指向向量从飞越地球到太阳
副总裁(rp) = p2000 (ip2 jdtdb2);
呼吁(1)= rp(1) /规范(rp);
呼吁(2)= rp(2) /规范(rp);
呼吁(3)= rp(3) /规范(rp);
sun_axisx =[0 5 *呼吁(1)];
sun_axisy =[0 5 *呼吁(2)];
sun_axisz =[0 5 *呼吁(3)];
%单位飞越地球的速度矢量
紫外线(1)=(1)副总裁/规范(vp);
紫外线(2)=(2)副总裁/规范(vp);
紫外线(3)=(3)副总裁/规范(vp);
uv_axisx =[0 5 *紫外线(1)];
uv_axisy =[0 5 *紫外线(2)];
uv_axisz =[0 5 *紫外线(3)];
dtof = 5000.0;
[ri_ho, vi_ho] = twobody2 (pmu (ip2)、-dtof rp2sc1, vp2sc1);
deltat1 = 2.0 * dtof / 300;
simtime1 = -deltat1;
为i = 1:1:301
simtime1 = simtime1 + deltat1;
%双曲线轨道位置向量“规范化”
(射频,vf) = twobody2 (pmu (ip2)、simtime1 ri_ho, vi_ho);
rp1_x (i) =射频(1)/代表(ip2);
rp1_y (i) =射频(2)/代表(ip2);
rp1_z (i) =射频(3)/代表(ip2);
结束
%建立轴向量
xaxisx = 1.5 [1];
xaxisy = [0 0];
xaxisz = [0 0];
yaxisx = [0 0];
yaxisy = 1.5 [1];
yaxisz = [0 0];
zaxisx = [0 0];
zaxisy = [0 0];
zaxisz = 1.5 [1];
图(2);
持有在;
网格在;
%的阴谋飞越地球
[x, y, z] =球(24);
h =冲浪(x, y, z);
colormap ((127/255 1 222/255));
集(h,“edgecolor”,(1 1 1));
%情节坐标系坐标轴
plot3 (xaxisx xaxisy xaxisz,“- r”,“线宽”1);
plot3 (yaxisx yaxisy yaxisz,“g”,“线宽”1);
plot3 (zaxisx zaxisy zaxisz,“- b”,“线宽”1);
%的情节单元向量指向太阳
plot3 (sun_axisx sun_axisy sun_axisz,“- y”,“线宽”1);
plot3 (uv_axisx uv_axisy uv_axisz,“- k”,“线宽”1);
%绘制planet-centered飞越双曲线轨道
plot3 (rp1_x rp1_y rp1_z,“- m”,“线宽”1);
%飞越双曲线的标签传入站
plot3 (rp1_x (1) rp1_y (1) rp1_z (1),“* m”);
%的标签近拱点飞越双曲线
plot3 (rp2sc1(1) /代表(ip2) rp2sc1(2) /代表(ip2) rp2sc1(3) /代表(ip2),“om”);
%标签情节在飞越地球半径(PR)
包含(“X坐标(PR)”,“字形大小”12);
ylabel (“Y坐标(PR)”,“字形大小”12);
zlabel (“Z坐标(PR)”,“字形大小”12);
标题(“Patched-conic重力辅助轨迹”,“字形大小”16);
轴平等的;
38岁的视图(30);
rotate3d在;
%下一行创建一个颜色用tiff预览eps图形文件
打印-depsc tiff r300 flyby_matlab2.eps
结束
错误:
索引超出了数组元素的数量(0)。
错误Untitled2(第742行)
rsc_ca = yevent (1:3);
