Why am I getting this error? Unrecognized function or variable 'snopt'. Error in Untitledidk (line 362) [x, f, inform, xmul, fmul] = snopt(Pls can some one show me a solution

12 views (last 30 days)
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

Answers (1)

Lokesh
Lokesh on 21 Mar 2024 at 6:35
Hi Mudrik,
As per my understanding, you are encountering an error related to using 'snopt' in your MATLAB code.
The error message you are encountering, "Unrecognized function or variable 'snopt'," indicates that MATLAB is unable to locate the 'snopt' function. This could be because the function is not in the current working directory or MATLAB's search path, or the 'SNOPT' software package, which includes this function, has not been properly installed or integrated into MATLAB.
Please ensure that the 'SNOPT' software package is correctly installed on your system. Refer to the following links for more information about the MATLAB interface for sparse nonlinear optimization software 'SNOPT'.

Categories

Find more on Environment in Help Center and File Exchange

Products


Release

R2021a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!