clear;
global emu smu xmu aunit ip1 ip2 ip3 jdtdb0 jdtdbi1
global iephem km ephname eq2000 pmu rep rsoi fbalt
global dv_launch dvm_launch dv_arrival dvm_arrival vinfm_in vinfm_out
global fbrp c3launch jdtdb1 jdtdb2
global otype rito1 vito1 rito2 vito2 dvh
global vinf_in vinf_out jdtdb_soi rp2sc vp2sc
% J2000 ecliptic-to-equatorial transformation matrix
eq2000 = [[1.000000000000000 0 0]; ...
[0 0.917482062069182 -0.397777155931914]; ...
[0 0.397777155931914 0.917482062069182]];
% define Astronomical Unit (kilometers)
aunit = 149597870.691;
% radians to degrees conversion factor
rtd = 180.0 / pi;
% define "reference" tdb julian date
jdtdb0 = julian(1, 1, 2000);
% initialize jpl ephemeris
ephname = 'de421.bin';
iephem = 1;
km = 1;
% define planet name vector
pdata = ['Mercury'; 'Venus '; 'Earth '; 'Mars '; ...
'Jupiter'; 'Saturn '; 'Uranus '; 'Neptune'; 'Pluto '];
pname = cellstr(pdata);
% planet gravitational constants (kilometer^3/second^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.295912208285591149e-03;
smu = xmu * aunit^3 / 86400.0^2;
xmu = 0.899701152970881167e-09;
emu = xmu * aunit^3 / 86400.0^2;
% planet equatorial radius (kilometers)
rep(1) = 2439.7;
rep(2) = 6051.9;
rep(3) = 6378.14;
rep(4) = 3397.0;
rep(5) = 71492.0;
rep(6) = 60268.0;
rep(7) = 25559.0;
rep(8) = 24764.0;
rep(9) = 1151.0;
% sphere-of-influence for each planet (kilometers)
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;
% begin simulation
clc; home;
fprintf('\nprogram flyby_matlab\n');
fprintf('\npatched-conic gravity-assist trajectory analysis\n\n');
% request departure tdb calendar date
fprintf('\ndeparture - calendar date guess\n');
('1 <= month <= 12, 1 <= day <= 31, year = all digits!');
3; 5; 2030;
%['\03, 05, 2030\n'] = getdate;
% compute departure tdb julian date
jdtdbi1 = julian(5, 5, 2030);
while(1)
fprintf('\nplease input the departure date search boundary in days\n');
ddays1 = input('30');
if (ddays1 >= 0.0)
break;
end
end
% request flyby tdb calendar date
fprintf('\n\nflyby - calendar date guess\n');
('1 <= month <= 12, 1 <= day <= 31, year = all digits!');
3; 5; 2037;
%[month, day, year] = getdate;
% compute flyby tdb julian date
jdtdbi2 = julian(03,05,2037);
%jdtdbi2 = julian(month, day, year);
while(1)
fprintf('\nplease input the flyby date search boundary in days\n');
ddays2 = input('?');
if (ddays2 >= 0.0)
break;
end
end
% request arrival tdb calendar date
fprintf('\n\narrival - calendar date guess\n');
('1 <= month <= 12, 1 <= day <= 31, year = all digits!');
8; 12; 2050;
%(month; day; year) = getdate;
% compute arrival tdb julian date
jdtdbi3 = julian(8, 12, 2050);
%jdtdbi3 = julian(month, day, year);
while(1)
fprintf('\nplease input the arrival date search boundary in days\n');
ddays3 = input('?');
if (ddays3 >= 0.0)
break;
end
end
% request departure, flyby and arrival planets
for i = 1:1:3
fprintf('\n planet menu\n');
fprintf('\n <1> Mercury');
fprintf('\n <2> Venus');
fprintf('\n <3> Earth');
fprintf('\n <4> Mars');
fprintf('\n <5> Jupiter');
fprintf('\n <6> Saturn');
fprintf('\n <7> Uranus');
fprintf('\n <8> Neptune');
fprintf('\n <9> Pluto');
if (i == 1)
while(1)
fprintf('\n\nplease select the departure planet\n');
ip1 = input('?');
if (ip1 >= 1 && ip1 <= 9)
break;
end
end
end
if (i == 2)
while(1)
fprintf('\n\nplease select the flyby planet\n');
ip2 = input('?');
if (ip2 >= 1 && ip2 <= 9)
break;
end
end
end
if (i == 3)
while(1)
fprintf('\n\nplease select the arrival planet\n');
ip3 = input('?');
if (ip3 >= 1 && ip3 <= 9)
break;
end
end
end
end
% request lower and upper bounds for the flyby altitude (kilometers)
while(1)
fprintf('\n\nplease input the lower bound for the flyby altitude (kilometers)\n');
fbalt_lwr = input('?');
if (fbalt_lwr >= 0.0)
break;
end
end
while(1)
fprintf('\n\nplease input the upper bound for the flyby altitude (kilometers)\n');
fbalt_upr = input('?');
if (fbalt_upr >= fbalt_lwr)
break;
end
end
% request type of optimization
while(1)
fprintf('\n\n optimization menu\n');
fprintf('\n <1> minimize departure delta-v\n');
fprintf('\n <2> minimize arrival delta-v\n');
fprintf('\n <3> minimize total delta-v\n');
fprintf('\n selection (1, 2 or 3)\n');
otype = input('?');
if (otype == 1 || otype == 2 || otype == 3)
break;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solve the patched-conic gravity assist problem %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc; home;
xg = zeros(3, 1);
% control variables initial guesses
% (departure, flyby, and arrival tdb julian dates)
xg(1) = jdtdbi1 - jdtdb0;
xg(2) = jdtdbi2 - jdtdb0;
xg(3) = jdtdbi3 - jdtdb0;
% bounds on control variables
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';
% bounds on objective function
flow(1) = -inf;
fupp(1) = +inf;
% bounds on flyby altitude
flow(2) = fbalt_lwr;
fupp(2) = fbalt_upr;
% bounds on v-infinity matching (equality) constraint
flow(3) = 0.0;
fupp(3) = 0.0;
flow = flow';
fupp = fupp';
xmul = zeros(3, 1);
xstate = zeros(3, 1);
fmul = zeros(3, 1);
fstate = zeros(3, 1);
% use snopt to find optimal solution
[x, f, inform, xmul, fmul] = snopt(xg, xlwr, xupr, xmul, xstate, ...
flow, fupp, fmul, fstate, 'fbyfunc');
% evaluate optimal solution
[~, g] = fbyfunc(x);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute orientation of the departure hyperbola
% (Earth mean equator and equinox of J2000)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dveq1 = eq2000 * dv_launch;
dveqm1 = norm(dveq1);
decl1 = 90.0 - rtd * acos(dveq1(3) / dveqm1);
rasc1 = rtd * atan3(dveq1(2), dveq1(1));
% compute tdb julian dates of optimal transfer
jdtdb1 = jdtdb0 + x(1);
jdtdb2 = jdtdb0 + x(2);
jdtdb3 = jdtdb0 + x(3);
% convert solution julian dates to calendar dates and tdb times
[cdstr1, tdbstr1] = jd2str(jdtdb1);
[cdstr2, tdbstr2] = jd2str(jdtdb2);
[cdstr3, tdbstr3] = jd2str(jdtdb3);
%%%%%%%%%%%%%%%%%
% print results %
%%%%%%%%%%%%%%%%%
fprintf('\nprogram flyby_matlab');
fprintf('\n--------------------\n');
fprintf('\npatched-conic gravity assist trajectory analysis\n');
if (otype == 1)
fprintf('\nminimize departure delta-v\n');
end
if (otype == 2)
fprintf('\nminimize arrival delta-v\n');
end
if (otype == 3)
fprintf('\nminimize total delta-v\n');
end
fprintf('\ndeparture planet ');
disp(pname(ip1));
fprintf('flyby planet ');
disp(pname(ip2));
fprintf('arrival planet ');
disp(pname(ip3));
% departure conditions
fprintf('\nLAUNCH CONDITIONS');
fprintf('\n=================\n');
fprintf('\ndeparture calendar date ');
disp(cdstr1);
fprintf('\ndeparture TDB time ');
disp(tdbstr1);
fprintf('\ndeparture TDB julian date %12.6f', jdtdb1);
fprintf('\n\nheliocentric launch delta-v vector and magnitude');
fprintf('\n(Earth mean ecliptic and equinox of J2000)');
fprintf('\n------------------------------------------\n');
fprintf('\nlaunch delta-vx %12.6f meters/second', 1000.0 * dv_launch(1));
fprintf('\nlaunch delta-vy %12.6f meters/second', 1000.0 * dv_launch(2));
fprintf('\nlaunch delta-vz %12.6f meters/second', 1000.0 * dv_launch(3));
fprintf('\n\nlaunch delta-v %12.6f meters/second\n', 1000.0 * dvm_launch);
fprintf('\nlaunch hyperbola characteristics');
fprintf('\n(Earth mean equator and equinox of J2000)');
fprintf('\n-----------------------------------------\n');
fprintf('\nasymptote right ascension %12.6f degrees\n', rasc1);
fprintf('\nasymptote declination %12.6f degrees\n', decl1);
fprintf('\nlaunch energy %12.6f kilometer^2/second^2\n', c3launch);
fprintf('\npost-impulse heliocentric orbital elements and state vector of the spacecraft');
fprintf('\n(Earth mean ecliptic and equinox of J2000)');
fprintf('\n------------------------------------------\n');
oev1 = eci2orb1(smu, rito1, vito1);
oeprint1(smu, oev1, 3);
svprint(rito1, vito1);
fprintf('heliocentric orbital elements and state vector of the departure planet');
fprintf('\n(Earth mean ecliptic and equinox of J2000)');
fprintf('\n------------------------------------------\n');
[rp1, vp1] = p2000(ip1, jdtdb1);
oev1 = eci2orb1(smu, rp1, vp1);
oeprint1(smu, oev1, 3);
svprint(rp1, vp1);
% flyby conditions
fprintf('\nPATCHED-CONIC FLYBY CONDITIONS');
fprintf('\n==============================\n');
fprintf('\nflyby calendar date ');
disp(cdstr2);
fprintf('\nflyby TDB time ');
disp(tdbstr2);
fprintf('\nflyby TDB julian date %12.6f\n', jdtdb2);
fprintf('\nlaunch-to-flyby time %12.6f days\n', jdtdb2 - jdtdb1);
fprintf('\nv-infinity in %12.6f meters/second', 1000.0 * vinfm_in);
fprintf('\nv-infinity out %12.6f meters/second\n', 1000.0 * vinfm_out);
fprintf('\nflyby altitude %12.6f kilometers\n', fbalt);
% turn angles
tmp = rep(ip2) * vinfm_in * vinfm_in / pmu(ip2);
tangle1 = 2.0 * asin(1.0 / (1.0 + tmp));
fprintf('\nmaximum turn angle %12.6f degrees', tangle1 * rtd);
tmp = (fbalt + rep(ip2)) * vinfm_in * vinfm_in / pmu(ip2);
tangle2 = 2.0 * asin(1.0 / (1.0 + tmp));
fprintf('\nactual turn angle %12.6f degrees\n', tangle2 * rtd);
fprintf('\nheliocentric delta-v %12.6f meters/second\n', 1000.0 * dvh);
% heliocentric deltavs
dvh_max = sqrt(pmu(ip2) / rep(ip2));
fprintf('\nmax heliocentric delta-v %12.6f meters/second\n', 1000.0 * dvh_max);
% planet-centered orbital elements at periapsis
oev_periapsis = fbhyper(pmu(ip2), vinf_in, vinf_out, fbrp);
[rp2sc1, vp2sc1] = orb2eci(pmu(ip2), oev_periapsis);
fprintf('\nplanet-centered orbital elements and state vector of the spacecraft at periapsis');
fprintf('\n(Earth mean ecliptic and equinox of J2000)');
fprintf('\n------------------------------------------\n');
oeprint1(pmu(ip2), oev_periapsis, 1);
svprint(rp2sc1, vp2sc1);
[bplane, ~, ~, ~] = rv2bplane(pmu(ip2), rp2sc1, vp2sc1);
fprintf('b-plane coordinates of the spacecraft at periapsis');
fprintf('\n(Earth mean ecliptic and equinox of J2000)');
fprintf('\n------------------------------------------\n');
bpprint(rp2sc1, vp2sc1, bplane);
cdecl_asy = cos(bplane(2));
crasc_asy = cos(bplane(3));
srasc_asy = sin(bplane(3));
sdecl_asy = sin(bplane(2));
fprintf('\nheliocentric orbital elements and state vector of the spacecraft');
fprintf('\n(Earth mean ecliptic and equinox of J2000)');
fprintf('\n------------------------------------------\n');
tau = 86400.0 * (jdtdb2 - jdtdb1);
[rp, vp] = twobody2(smu, tau, rito1, vito1);
oev = eci2orb1(smu, rp, vp);
oeprint1(smu, oev, 3);
svprint(rp, vp);
fprintf('\nheliocentric orbital elements and state vector of the flyby planet');
fprintf('\n(Earth mean ecliptic and equinox of J2000)');
fprintf('\n------------------------------------------\n');
[rp, vp] = p2000(ip2, jdtdb2);
oev = eci2orb1(smu, rp, vp);
oeprint1(smu, oev, 3);
svprint(rp, vp);
% arrival conditions
fprintf('\nARRIVAL CONDITIONS');
fprintf('\n==================\n');
fprintf('\narrival calendar date ');
disp(cdstr3);
fprintf('\narrival TDB time ');
disp(tdbstr3);
fprintf('\narrival TDB julian date %12.6f\n', jdtdb3);
fprintf('\nflyby-to-arrival time %12.6f days\n', jdtdb3 - jdtdb2);
fprintf('\nheliocentric arrival delta-v vector and magnitude');
fprintf('\n(Earth mean ecliptic and equinox of J2000)');
fprintf('\n------------------------------------------\n');
fprintf('\narrival delta-vx %12.6f meters/second', 1000.0 * dv_arrival(1));
fprintf('\narrival delta-vy %12.6f meters/second', 1000.0 * dv_arrival(2));
fprintf('\narrival delta-vz %12.6f meters/second', 1000.0 * dv_arrival(3));
fprintf('\n\narrival delta-v %12.6f meters/second\n', 1000.0 * dvm_arrival);
% two-body propagation of the second leg of the transfer orbit
[rf3, vf3] = twobody2 (smu, 86400 * (jdtdb3 - jdtdb2), rito2, vito2);
fprintf('\npre-impulse heliocentric orbital elements and state vector of the spacecraft');
fprintf('\n(Earth mean ecliptic and equinox of J2000)');
fprintf('\n------------------------------------------\n');
oev3 = eci2orb1(smu, rf3, vf3);
oeprint1(smu, oev3, 3);
svprint(rf3, vf3);
fprintf('\npost-impulse heliocentric orbital elements and state vector of the spacecraft');
fprintf('\n(Earth mean ecliptic and equinox of J2000)');
fprintf('\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);
fprintf('heliocentric orbital elements and state vector of the arrival planet');
fprintf('\n(Earth mean ecliptic and equinox of J2000)');
fprintf('\n------------------------------------------\n');
[rp, vp] = p2000(ip3, jdtdb3);
oev = eci2orb1(smu, rp, vp);
oeprint1(smu, oev, 3);
svprint(rp, vp);
fprintf('\nMISSION SUMMARY');
fprintf('\n===============\n');
fprintf('\ntotal delta-v %12.6f meters/second\n', ...
1000.0 * (dvm_launch + dvm_arrival));
fprintf('\ntotal energy %12.6f kilometer^2/second^2\n', dvm_launch^2 + dvm_arrival^2);
fprintf('\ntotal mission duration %12.6f days\n\n', jdtdb3 - jdtdb1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% two-body integration of the trajectory from
% first impulse to SOI entrance at flyby planet
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% set up for ode45
options = odeset('RelTol', 1.0e-6, 'AbsTol', 1.0e-8, 'Events', @soi_event);
% define maximum search time (seconds)
tof = 86400.0 * (jdtdb2 - jdtdb1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solve for flyby planet SOI entrance conditions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[~, ~, tevent, ~, ~] = ode45(@twobody_heqm, [0 tof], [rito1 vito1], options);
jdtdb_soi = jdtdb1 + tevent / 86400.0;
[cdstr_ca, utstr_ca] = jd2str(jdtdb_soi);
fprintf('\nPATCHED-CONIC SOI ENTRANCE CONDITIONS');
fprintf('\n=====================================\n');
fprintf('\ncalendar date ');
disp(cdstr_ca);
fprintf('\nTDB time ');
disp(utstr_ca);
fprintf('\nTDB julian date %12.6f\n', jdtdb_soi);
fprintf('\nplanet-centered orbital elements and state vector of the spacecraft');
fprintf('\n(Earth mean ecliptic and equinox of J2000)');
fprintf('\n------------------------------------------\n');
oev_soi = eci2orb1(pmu(ip2), rp2sc, vp2sc);
oeprint1(pmu(ip2), oev_soi, 1);
svprint(rp2sc, vp2sc);
[bplane, ~, ~, ~] = rv2bplane(pmu(ip2), rp2sc, vp2sc);
fprintf('b-plane coordinates at sphere-of-influence');
fprintf('\n(Earth mean ecliptic and equinox of J2000)');
fprintf('\n------------------------------------------\n');
bpprint(rp2sc, vp2sc, bplane);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% integrate trajectory from SOI entrance
% to closest approach to the flyby planet
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% set up for ode45
options = odeset('RelTol', 1.0e-10, 'AbsTol', 1.0e-10, 'Events', @fpa_event);
% define maximum search time (seconds)
tof = 10.0 * 86400.0;
[~, ~, tevent, yevent, ie] = ode45(@peqm, [0 tof], [rp2sc vp2sc], options);
% extract state vector at closest approach
rsc_ca = yevent(1:3);
vsc_ca = yevent(4:6);
% B-plane coordinates of flyby hyperbola at closest approach
[bplane, rv, tv, ibperr] = rv2bplane(pmu(ip2), rsc_ca, vsc_ca);
% tdb julian date at closest approach
jdtdb_ca = jdtdb_soi + tevent / 86400.0;
[cdstr_ca, utstr_ca] = jd2str(jdtdb_ca);
fprintf('\n\nNUMERICALLY INTEGRATED CLOSEST APPROACH CONDITIONS');
fprintf('\n==================================================\n');
fprintf('\ncalendar date ');
disp(cdstr_ca);
fprintf('\nTDB time ');
disp(utstr_ca);
fprintf('\nTDB julian date %12.6f\n', jdtdb_ca);
fprintf('\nplanet-centered orbital elements and state vector of the spacecraft');
fprintf('\n(Earth mean ecliptic and equinox of J2000)');
fprintf('\n------------------------------------------\n');
oev_ca = eci2orb1(pmu(ip2), rsc_ca, vsc_ca);
oeprint1(pmu(ip2), oev_ca, 1);
svprint(rsc_ca, vsc_ca);
fprintf('b-plane coordinates at closest approach');
fprintf('\n(Earth mean ecliptic and equinox of J2000)');
fprintf('\n------------------------------------------\n');
bpprint(rsc_ca, vsc_ca, bplane);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% heliocentric and flyby graphics %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
while(1)
fprintf('\n\nwould you like to create graphics for this mission (y = yes, n = no)\n');
slct = input('?', 's');
if (slct == 'y' || slct == 'n')
break;
end
end
if (slct == 'y')
while(1)
fprintf('\n\nplease input the plot step size (days)\n');
deltat = input('?');
if (deltat > 0)
break;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot heliocentric planetary orbits and transfer trajectory
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1);
% departure planet semimajor axis and period
[r1, v1] = p2000(ip1, jdtdb1);
oev1 = eci2orb1(smu, r1, v1);
period1 = 2.0 * pi * oev1(1) * sqrt(oev1(1) / smu) / 86400.0;
% flyby planet semimajor axis and period
[r2, v2] = p2000(ip2, jdtdb1);
oev2 = eci2orb1(smu, r2, v2);
period2 = 2.0 * pi * oev2(1) * sqrt(oev2(1) / smu) / 86400.0;
% arrival planet semimajor axis and period
[r3, v3] = p2000(ip3, jdtdb1);
oev3 = eci2orb1(smu, r3, v3);
period3 = 2.0 * pi * oev3(1) * sqrt(oev3(1) / smu) / 86400.0;
% number of "planet" points to plot
npts1 = fix(period1 / deltat);
npts2 = fix(period2 / deltat);
npts3 = fix(period3 / deltat);
% departure planet orbit
for i = 0:1:npts1
jdtdb = jdtdb1 + i * deltat;
[r1, ~] = p2000(ip1, jdtdb);
x1(i + 1) = r1(1) / aunit;
y1(i + 1) = r1(2) / aunit;
end
% compute last data point
[r1, v1] = p2000(ip1, jdtdb1 + period1);
x1(npts1 + 1) = r1(1) / aunit;
y1(npts1 + 1) = r1(2) / aunit;
% flyby planet heliocentric orbit
for i = 0:1:npts2
jdtdb = jdtdb1 + i * deltat;
[r2, ~] = p2000(ip2, jdtdb);
x2(i + 1) = r2(1) / aunit;
y2(i + 1) = r2(2) / aunit;
end
% compute last data point
[r2, v2] = p2000(ip2, jdtdb1 + period2);
x2(npts2 + 1) = r2(1) / aunit;
y2(npts2 + 1) = r2(2) / aunit;
% arrival planet heliocentric orbit
for i = 0:1:npts3
jdtdb = jdtdb1 + i * deltat;
[r3, ~] = p2000(ip3, jdtdb);
x3(i + 1) = r3(1) / aunit;
y3(i + 1) = r3(2) / aunit;
end
% compute last data point
[r3, v3] = p2000(ip3, jdtdb1 + period3);
x3(npts3 + 1) = r3(1) / aunit;
y3(npts3 + 1) = r3(2) / aunit;
% first leg of transfer orbit
npts4 = fix((jdtdb2 - jdtdb1) / deltat);
for i = 0:1:npts4
tau = 86400.0 * i * deltat;
[rft1, ~] = twobody2(smu, tau, rito1, vito1);
x4(i + 1) = rft1(1) / aunit;
y4(i + 1) = rft1(2) / aunit;
end
% compute last data point of first leg
tau = 86400.0 * (jdtdb2 - jdtdb1);
[rft1, vft1] = twobody2(smu, tau, rito1, vito1);
x4(npts4 + 1) = rft1(1) / aunit;
y4(npts4 + 1) = rft1(2) / aunit;
% second leg of heliocentric transfer orbit
npts5 = fix((jdtdb3 - jdtdb2) / deltat);
for i = 0:1:npts5
tau = 86400.0 * i * deltat;
[rft2, vft2] = twobody2(smu, tau, rito2, vito2);
x5(i + 1) = rft2(1) / aunit;
y5(i + 1) = rft2(2) / aunit;
end
% compute last data point of second leg
tau = 86400.0 * (jdtdb3 - jdtdb2);
[rfto2, vfto2] = twobody2(smu, tau, rito2, vito2);
x5(npts5 + 1) = rfto2(1) / aunit;
y5(npts5 + 1) = rfto2(2) / aunit;
% determine vernal equinox "size"
xve = oev1(1) / aunit;
if (oev2(1) > oev1(1))
xve = oev2(1) / aunit;
end
if (oev3(1) > oev2(1))
xve = oev3(1) / aunit;
end
% plot heliocentric planet orbits and patched-conic transfer trajectory
hold on;
plot(x1, y1, '.b');
plot(x1, y1, '-b');
plot(x2, y2, '.g');
plot(x2, y2, '-g');
plot(x3, y3, '.r');
plot(x3, y3, '-r');
plot(x4, y4, '.k');
plot(x4, y4, '-k');
plot(x5, y5, '.k');
plot(x5, y5, '-k');
plot(0, 0, 'hy', 'MarkerSize', 10);
% plot and label vernal equinox direction
line ([0.05, xve], [0, 0], 'Color', 'black');
text(1.05 * xve, 0, '\Upsilon');
% label launch, flyby and arrival locations
plot(x4(1), y4(1), '*k');
text(x4(1) + 0.06, y4(1) - 0.01, 'D');
plot(x5(1), y5(1), '*k');
text(x5(1) + 0.06, y5(1) + 0.01, 'F');
plot(x5(npts5 + 1), y5(npts5 + 1), '*k');
text(x5(npts5 + 1) + 0.06, y5(npts5 + 1) - 0.01, 'A');
axis equal;
% display launch, flyby and arrival dates
text(0.85 * xve, -xve + 0.8, 'Departure', 'FontSize', 8);
text(0.875 * xve, -xve + 0.7, pname(ip1), 'FontSize', 8);
text(0.875 * xve, -xve + 0.6, cdstr1, 'FontSize', 8);
text(0.85 * xve, -xve + 0.4, 'Arrival', 'FontSize', 8);
text(0.875 * xve, -xve + 0.3, pname(ip3), 'FontSize', 8);
text(0.875 * xve, -xve + 0.2, cdstr3, 'FontSize', 8);
text(0.85 * xve, xve - 0.2, 'Flyby', 'FontSize', 8);
text(0.875 * xve, xve - 0.3, pname(ip2), 'FontSize', 8);
text(0.875 * xve, xve - 0.4, cdstr2, 'FontSize', 8);
text(0.875 * xve, xve - 0.5, tdbstr2, 'FontSize', 8);
% label plot and axes in astronomical units (AU)
xlabel('X coordinate (AU)', 'FontSize', 12);
ylabel('Y coordinate (AU)', 'FontSize', 12);
title('Patched-conic Gravity Assist Trajectory', 'FontSize', 16);
zoom on;
% the next line creates a color eps graphics file with tiff preview
print -depsc -tiff -r300 flyby_matlab1.eps
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% create three-dimensional graphics of the flyby
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% unit pointing vector from the flyby planet to the sun
[rp, vp] = p2000(ip2, jdtdb2);
up2s(1) = -rp(1) / norm(rp);
up2s(2) = -rp(2) / norm(rp);
up2s(3) = -rp(3) / norm(rp);
sun_axisx = [0 5 * up2s(1)];
sun_axisy = [0 5 * up2s(2)];
sun_axisz = [0 5 * up2s(3)];
% unit velocity vector of the flyby planet
uv(1) = vp(1) / norm(vp);
uv(2) = vp(2) / norm(vp);
uv(3) = vp(3) / norm(vp);
uv_axisx = [0 5 * uv(1)];
uv_axisy = [0 5 * uv(2)];
uv_axisz = [0 5 * uv(3)];
dtof = 5000.0;
[ri_ho, vi_ho] = twobody2 (pmu(ip2), -dtof, rp2sc1, vp2sc1);
deltat1 = 2.0 * dtof / 300;
simtime1 = -deltat1;
for i = 1:1:301
simtime1 = simtime1 + deltat1;
% hyperbolic orbit "normalized" position vector
[rf, vf] = twobody2 (pmu(ip2), simtime1, ri_ho, vi_ho);
rp1_x(i) = rf(1) / rep(ip2);
rp1_y(i) = rf(2) / rep(ip2);
rp1_z(i) = rf(3) / rep(ip2);
end
% create axes vectors
xaxisx = [1 1.5];
xaxisy = [0 0];
xaxisz = [0 0];
yaxisx = [0 0];
yaxisy = [1 1.5];
yaxisz = [0 0];
zaxisx = [0 0];
zaxisy = [0 0];
zaxisz = [1 1.5];
figure(2);
hold on;
grid on;
% plot flyby planet
[x, y, z] = sphere(24);
h = surf(x, y, z);
colormap([127/255 1 222/255]);
set (h, 'edgecolor', [1 1 1]);
% plot coordinate system axes
plot3(xaxisx, xaxisy, xaxisz, '-r', 'LineWidth', 1);
plot3(yaxisx, yaxisy, yaxisz, '-g', 'LineWidth', 1);
plot3(zaxisx, zaxisy, zaxisz, '-b', 'LineWidth', 1);
% plot unit pointing vector to the sun
plot3(sun_axisx, sun_axisy, sun_axisz, '-y', 'LineWidth', 1);
plot3(uv_axisx, uv_axisy, uv_axisz, '-k', 'LineWidth', 1);
% plot planet-centered flyby hyperbolic orbit
plot3(rp1_x, rp1_y, rp1_z, '-m', 'LineWidth', 1);
% label incoming leg of flyby hyperbola
plot3(rp1_x(1), rp1_y(1), rp1_z(1), '*m');
% label periapsis of flyby hyperbola
plot3(rp2sc1(1) / rep(ip2), rp2sc1(2) / rep(ip2), rp2sc1(3) / rep(ip2), 'om');
% label plot in flyby planet radii (PR)
xlabel('X coordinate (PR)', 'FontSize', 12);
ylabel('Y coordinate (PR)', 'FontSize', 12);
zlabel('Z coordinate (PR)', 'FontSize', 12);
title('Patched-conic Gravity Assist Trajectory', 'FontSize', 16);
axis equal;
view(38, 30);
rotate3d on;
% the next line creates a color eps graphics file with tiff preview
print -depsc -tiff -r300 flyby_matlab2.eps
end
