256 SOFTWARE/matlab2tikz/m2tInputParser.m
View file
@@ -0,0 +1,256 @@
function parser = m2tInputParser()
%MATLAB2TIKZINPUTPARSER Input parsing for matlab2tikz..
% This implementation exists because Octave is lacking one.

% Copyright (c) 2008--2014 Nico Schlömer
% All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%
% * Redistributions of source code must retain the above copyright
% notice, this list of conditions and the following disclaimer.
% * Redistributions in binary form must reproduce the above copyright
% notice, this list of conditions and the following disclaimer in
% the documentation and/or other materials provided with the distribution
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.
% =========================================================================
% Initialize the structure.
parser = struct ();
% Public Properties
parser.Results = {};
% Enabel/disable parameters case sensitivity.
parser.CaseSensitive = false;
% Keep parameters not defined by the constructor.
parser.KeepUnmatched = false;
% Enable/disable warning for parameters not defined by the constructor.
parser.WarnUnmatched = true;
% Enable/disable passing arguments in a structure.
parser.StructExpand = true;
% Names of parameters defined in input parser constructor.
parser.Parameters = {};
% Names of parameters not defined in the constructor.
parser.Unmatched = struct ();
% Names of parameters using default values.
parser.UsingDefaults = {};
% Names of deprecated parameters and their alternatives
parser.DeprecatedParameters = struct();

% Handles for functions that act on the object.
parser.addRequired = @addRequired;
parser.addOptional = @addOptional;
parser.addParamValue = @addParamValue;
parser.deprecateParam = @deprecateParam;
parser.parse = @parse;

% Initialize the parser plan
parser.plan = {};
end
% =========================================================================
function p = parser_plan (q, arg_type, name, default, validator)
p = q;
plan = p.plan;
if (isempty (plan))
plan = struct ();
n = 1;
else
n = numel (plan) + 1;
end
plan(n).type = arg_type;
plan(n).name = name;
plan(n).default = default;
plan(n).validator = validator;
p.plan = plan;
end
% =========================================================================
function p = addRequired (p, name, validator)
p = parser_plan (p, 'required', name, [], validator);
end
% =========================================================================
function p = addOptional (p, name, default, validator)
p = parser_plan (p, 'optional', name, default, validator);
end
% =========================================================================
function p = addParamValue (p, name, default, validator)
p = parser_plan (p, 'paramvalue', name, default, validator);
end
% =========================================================================
function p = deprecateParam (p, name, alternatives)
if isempty(alternatives)
alternatives = {};
elseif ischar(alternatives)
alternatives = {alternatives}; % make cellstr
elseif ~iscellstr(alternatives)
error('m2tInputParser:BadAlternatives',...
'Alternatives for a deprecated parameter must be a char or cellstr');
end
p.DeprecatedParameters.(name) = alternatives;
end
% =========================================================================
function p = parse (p, varargin)
plan = p.plan;
results = p.Results;
using_defaults = {};
if (p.CaseSensitive)
name_cmp = @strcmp;
else
name_cmp = @strcmpi;
end
if (p.StructExpand)
k = find (cellfun (@isstruct, varargin));
for m = numel(k):-1:1
n = k(m);
s = varargin{n};
c = [fieldnames(s).'; struct2cell(s).'];
c = c(:).';
if (n > 1 && n < numel (varargin))
varargin = horzcat (varargin(1:n-1), c, varargin(n+1:end));
elseif (numel (varargin) == 1)
varargin = c;
elseif (n == 1);
varargin = horzcat (c, varargin(n+1:end));
else % n == numel (varargin)
varargin = horzcat (varargin(1:n-1), c);
end
end
end
if (isempty (results))
results = struct ();
end
type = {plan.type};
n = find( strcmp( type, 'paramvalue' ) );
m = setdiff (1:numel( plan ), n );
plan = plan ([n,m]);
for n = 1 : numel (plan)
found = false;
results.(plan(n).name) = plan(n).default;
if (~ isempty (varargin))
switch plan(n).type
case 'required'
found = true;
if (strcmpi (varargin{1}, plan(n).name))
varargin(1) = [];
end
value = varargin{1};
varargin(1) = [];
case 'optional'
m = find (cellfun (@ischar, varargin));
k = find (name_cmp (plan(n).name, varargin(m)));
if (isempty (k) && validate_arg (plan(n).validator, varargin{1}))
found = true;
value = varargin{1};
varargin(1) = [];
elseif (~ isempty (k))
m = m(k);
found = true;
value = varargin{max(m)+1};
varargin(union(m,m+1)) = [];
end
case 'paramvalue'
m = find( cellfun (@ischar, varargin) );
k = find (name_cmp (plan(n).name, varargin(m)));
if (~ isempty (k))
found = true;
m = m(k);
value = varargin{max(m)+1};
varargin(union(m,m+1)) = [];
end
otherwise
error( sprintf ('%s:parse', mfilename), ...
'parse (%s): Invalid argument type.', mfilename ...
)
end
end
if (found)
if (validate_arg (plan(n).validator, value))
results.(plan(n).name) = value;
else
error( sprintf ('%s:invalidinput', mfilename), ...
'%s: Input argument ''%s'' has invalid value.\n', mfilename, plan(n).name ...
);
end
p.Parameters = union (p.Parameters, {plan(n).name});
elseif (strcmp (plan(n).type, 'required'))
error( sprintf ('%s:missinginput', mfilename), ...
'%s: input ''%s'' is missing.\n', mfilename, plan(n).name ...
);
else
using_defaults = union (using_defaults, {plan(n).name});
end
end

if ~isempty(varargin)
% Include properties that do not match specified properties
for n = 1:2:numel(varargin)
if ischar(varargin{n})
if p.KeepUnmatched
results.(varargin{n}) = varargin{n+1};
end
if p.WarnUnmatched
warning(sprintf('%s:unmatchedArgument',mfilename), ...
'Ignoring unknown argument "%s"', varargin{n});
end
p.Unmatched.(varargin{n}) = varargin{n+1};
else
error (sprintf ('%s:invalidinput', mfilename), ...
'%s: invalid input', mfilename)
end
end
end

% Store the results of the parsing
p.Results = results;
p.UsingDefaults = using_defaults;

warnForDeprecatedParameters(p);
end
% =========================================================================
function result = validate_arg (validator, arg)
try
result = validator (arg);
catch %#ok
result = false;
end
end
% =========================================================================
function warnForDeprecatedParameters(p)
usedDeprecatedParameters = intersect(p.Parameters, fieldnames(p.DeprecatedParameters));

for iParam = 1:numel(usedDeprecatedParameters)
oldParameter = usedDeprecatedParameters{iParam};
alternatives = p.DeprecatedParameters.(oldParameter);

switch numel(alternatives)
case 0
replacements = '';
case 1
replacements = ['''' alternatives{1} ''''];
otherwise
replacements = deblank(sprintf('''%s'' and ',alternatives{:}));
replacements = regexprep(replacements,' and$','');
end
if ~isempty(replacements)
replacements = sprintf('From now on, please use %s to control the output.\n',replacements);
end

message = ['\n===============================================================================\n', ...
'You are using the deprecated parameter ''%s''.\n', ...
'%s', ...
'==============================================================================='];
warning('matlab2tikz:deprecatedParameter', ...
message, oldParameter, replacements);

end
end
5,555 SOFTWARE/matlab2tikz/matlab2tikz.m
View file

Large diffs are not rendered by default.

142 SOFTWARE/matlab2tikz/private/m2tUpdater.m
View file
@@ -0,0 +1,142 @@
function updater(name, fileExchangeUrl, version, verbose, env)
%UPDATER Auto-update matlab2tikz.
% Only for internal usage.

% Copyright (c) 2012--2014, Nico Schlömer <nico.schloemer@gmail.com>
% All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%
% * Redistributions of source code must retain the above copyright
% notice, this list of conditions and the following disclaimer.
% * Redistributions in binary form must reproduce the above copyright
% notice, this list of conditions and the following disclaimer in
% the documentation and/or other materials provided with the distribution
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.
% =========================================================================
try
html = urlread([fileExchangeUrl, '/all_files']);
catch %#ok
% Couldn't load the URL -- never mind.
html = '';
end
% Search for a string "/version-1.6.3" in the HTML. This assumes
% that the package author has added a file by that name to
% to package. This is a rather dirty hack around FileExchange's
% lack of native versioning information.
mostRecentVersion = regexp(html, '/version-(\d+\.\d+\.\d+)', 'tokens');
if ~isempty(mostRecentVersion)
if isVersionBelow(env, version, mostRecentVersion{1}{1})
userInfo(verbose, '**********************************************\n');
userInfo(verbose, 'New version available! (%s)\n', mostRecentVersion{1}{1});
userInfo(verbose, '**********************************************\n');

reply = input([' *** Would you like ', name, ' to self-upgrade? y/n [n]:'],'s');
if strcmp(reply, 'y')
% Download the files and unzip its contents into the folder
% above the folder that contains the current script.
% This assumes that the file structure is something like
%
% src/matlab2tikz.m
% src/[...]
% AUTHORS
% ChangeLog
% [...]
%
% on the hard drive and the zip file. In particular, this assumes
% that the folder on the hard drive is writable by the user
% and that matlab2tikz.m is not symlinked from some other place.
pathstr = fileparts(mfilename('fullpath'));
targetPath = [pathstr, filesep, '..', filesep];
userInfo(verbose, ['Downloading and unzipping to ', targetPath, '...']);
upgradeSuccess = false;
try
unzippedFiles = unzip([fileExchangeUrl, '?download=true'], targetPath);
upgradeSuccess = true; %~isempty(unzippedFiles);
userInfo(verbose, 'done.');
catch
userInfo(verbose, 'FAILED.');
end
if upgradeSuccess
% TODO explicitly delete all of the old content
% Delete old version number file.
versionFile = [pathstr, filesep, 'version-', version];
if exist(versionFile, 'file') == 2
delete(versionFile);
end
% TODO anything better than error()?
error('Upgrade successful. Please re-execute.');
else
error('Upgrade failed.');
end
end
userInfo(verbose, '');
end
end
end
% =========================================================================
function isBelow = isVersionBelow(env, versionA, versionB)
% Checks if version string or vector versionA is smaller than
% version string or vector versionB.

vA = versionArray(env, versionA);
vB = versionArray(env, versionB);

isBelow = false;
for i = 1:min(length(vA), length(vB))
if vA(i) > vB(i)
isBelow = false;
break;
elseif vA(i) < vB(i)
isBelow = true;
break
end
end

end
% =========================================================================
function arr = versionArray(env, str)
% Converts a version string to an array, e.g.,
% '2.62.8.1' to [2, 62, 8, 1].

if ischar(str)
if strcmpi(env, 'MATLAB')
split = regexp(str, '\.', 'split');
elseif strcmpi(env, 'Octave')
split = strsplit(str, '.');
end
arr = str2num(char(split)); %#ok
else
arr = str;
end

end
% =========================================================================
function userInfo(verbose, message, varargin)
% Display usage information.

if ~verbose
return
end

mess = sprintf(message, varargin{:});

% Replace '\n' by '\n *** ' and print.
mess = strrep( mess, sprintf('\n'), sprintf('\n *** ') );
fprintf( ' *** %s\n', mess );

end
% =========================================================================
BIN +75 KB SOFTWARE/orbital_model/LaTeX/genPDF.pdf
View file
Binary file not shown.
40 SOFTWARE/orbital_model/LaTeX/genPDF.tex
View file
@@ -0,0 +1,40 @@
\documentclass[12pt]{article}

\usepackage{pgfplots}
\usepackage{caption}
\usepackage{amsmath,amssymb,amsthm,multirow,algorithm,algorithmic,amsfonts}



\newlength\figureheight
\newlength\figurewidth

\begin{document}
\thispagestyle{empty}
\begin{figure}
\centering
\setlength\figureheight{5cm}
\setlength\figurewidth{13cm}
\input{n_CL.tikz}
\caption*{The minimal load for entering an orbit around Mars with no active control, a constant $C_D$ ($1.25$) and variable $C_L$ at different initial velocities.}
\end{figure}

\begin{figure}
\centering
\setlength\figureheight{5cm}
\setlength\figurewidth{12.5cm}
\input{n_CD.tikz}
\caption*{The minimal load for entering an orbit around Mars with no active control, a constant $C_L$ ($0.2$) and variable $C_D$ at different initial velocities.}
\end{figure}

\begin{figure}
\centering
\setlength\figureheight{5cm}
\setlength\figurewidth{10cm}
\input{range_plot.tikz}
\caption*{The range of $R_x$ for entering an orbit around Mars with no active control, a constant $C_D$ ($1.25$) and $C_L$ ($-0.4$) at an initial velocity of $7000$ $\left[\frac{m}{s}\right]$.}
\end{figure}



\end{document}
68 SOFTWARE/orbital_model/LaTeX/n_CD.tikz
View file
@@ -0,0 +1,68 @@
% This file was created by matlab2tikz.
% Minimal pgfplots version: 1.3
%
\definecolor{mycolor1}{rgb}{0.20810,0.16630,0.52920}%
\definecolor{mycolor2}{rgb}{0.19855,0.72140,0.63095}%
%
\begin{tikzpicture}

\begin{axis}[%
width=\figurewidth,
height=0.812522\figureheight,
at={(0\figurewidth,0\figureheight)},
scale only axis,
every outer x axis line/.append style={black},
every x tick label/.append style={font=\color{black}},
xmin=0.8,
xmax=1.25,
xlabel={$C_D \left[-\right]$},
xmajorgrids,
every outer y axis line/.append style={black},
every y tick label/.append style={font=\color{black}},
ymin=20,
ymax=90,
ylabel={$a_{min} \left[\frac{m}{s^2}\right]$},
ymajorgrids,
axis x line*=bottom,
axis y line*=left,
legend style={at={(0.5,1.03)},anchor=south,legend cell align=left,align=left,draw=black}
]
\addplot [color=mycolor1,solid,mark=+,mark options={solid}]
table[row sep=crcr]{%
0.8 53.94987918\\
0.85 52.7488605\\
0.9 51.75588249\\
0.95 50.86001367\\
1 50.03539488\\
1.05 49.19858226\\
1.1 48.35557953\\
1.15 47.68183854\\
1.2 47.06358291\\
1.25 46.51698933\\
};
\addlegendentry{6000 m/s};

\addplot [color=mycolor2,solid,mark=o,mark options={solid}]
table[row sep=crcr]{%
0.8 84.16472823\\
0.85 81.43706754\\
0.9 79.40015838\\
0.95 77.66121816\\
1 76.22476929\\
1.05 74.92213863\\
1.1 73.88776242\\
1.15 72.94044015\\
1.2 71.91335277\\
1.25 70.94272194\\
};
\addlegendentry{7000 m/s};

\addplot [color=red,dash pattern=on 1pt off 3pt on 3pt off 3pt]
table[row sep=crcr]{%
0.8 29.43\\
1.25 29.43\\
};
\addlegendentry{$\text{29.43 m/s}^\text{2}\text{ (3g)}$};

\end{axis}
\end{tikzpicture}%
110 SOFTWARE/orbital_model/LaTeX/n_CL.tikz
View file
@@ -0,0 +1,110 @@
% This file was created by matlab2tikz.
% Minimal pgfplots version: 1.3
%
\definecolor{mycolor1}{rgb}{0.20810,0.16630,0.52920}%
\definecolor{mycolor2}{rgb}{0.02650,0.61370,0.81350}%
\definecolor{mycolor3}{rgb}{0.64730,0.74560,0.41880}%
%
\begin{tikzpicture}

\begin{axis}[%
width=0.95092\figurewidth,
height=\figureheight,
at={(0\figurewidth,0\figureheight)},
scale only axis,
every outer x axis line/.append style={black},
every x tick label/.append style={font=\color{black}},
xmin=-0.6,
xmax=0.8,
xlabel={$C_L \left[-\right]$},
xmajorgrids,
every outer y axis line/.append style={black},
every y tick label/.append style={font=\color{black}},
ymin=0,
ymax=58.86,
ylabel={$a_{min} \left[\frac{m}{s^2}\right]$},
ymajorgrids,
axis x line*=bottom,
axis y line*=left,
legend style={at={(0.97,0.03)},anchor=south east,legend cell align=left,align=left,draw=black}
]
\addplot [color=mycolor1,solid,mark=+,mark options={solid}]
table[row sep=crcr]{%
-0.5 14.13288441\\
-0.44 14.91799833\\
-0.38 15.64377156\\
-0.31 16.49330775\\
-0.25 17.43646077\\
-0.19 18.46461744\\
-0.12 19.57588443\\
-0.06 20.92749642\\
0 22.35420396\\
0.06 24.00941583\\
0.12 25.81731054\\
0.19 27.83173518\\
0.25 30.03779817\\
0.31 32.3486712\\
0.38 34.83845901\\
0.44 37.53788652\\
0.5 40.44023388\\
0.56 43.59989754\\
0.63 47.00644947\\
};
\addlegendentry{5000 m/s};

\addplot [color=mycolor2,solid,mark=o,mark options={solid}]
table[row sep=crcr]{%
-0.5 19.63585296\\
-0.44 20.80301733\\
-0.38 22.11463395\\
-0.31 23.54530473\\
-0.25 25.22311884\\
-0.19 27.25214076\\
-0.12 29.47042701\\
-0.06 32.22318168\\
0 35.07137784\\
0.06 38.25335925\\
0.12 41.75607861\\
0.19 45.65975229\\
0.25 49.97838897\\
0.31 54.72550683\\
0.38 59.8857336\\
0.44 65.56112271\\
0.5 71.45858079\\
0.56 77.78018403\\
0.63 84.47813811\\
};
\addlegendentry{6000 m/s};

\addplot [color=mycolor3,solid,mark=asterisk,mark options={solid}]
table[row sep=crcr]{%
-0.44 27.42435531\\
-0.38 29.48981157\\
-0.31 31.69494261\\
-0.25 34.41160629\\
-0.19 37.7688924\\
-0.12 41.526711\\
-0.06 45.8544906\\
0 50.66340165\\
0.06 56.08215135\\
0.12 62.30730267\\
0.19 69.45922431\\
0.25 77.15041146\\
0.31 85.20808059\\
0.38 94.00090131\\
0.44 103.55305527\\
0.5 114.01320321\\
0.56 125.52448284\\
0.63 138.26689785\\
};
\addlegendentry{7000 m/s};

\addplot [color=red,dash pattern=on 1pt off 3pt on 3pt off 3pt]
table[row sep=crcr]{%
-0.6 29.43\\
0.8 29.43\\
};
\addlegendentry{$\text{29.43 m/s}^\text{2}\text{ (3g)}$};

\end{axis}
\end{tikzpicture}%
51 SOFTWARE/orbital_model/LaTeX/range_plot.tikz
View file
@@ -0,0 +1,51 @@
% This file was created by matlab2tikz.
% Minimal pgfplots version: 1.3
%
\definecolor{mycolor1}{rgb}{0.07947,0.51590,0.83283}%
%
\begin{tikzpicture}

\begin{axis}[%
width=0.95092\figurewidth,
height=\figureheight,
at={(0\figurewidth,0\figureheight)},
scale only axis,
every outer x axis line/.append style={black},
every x tick label/.append style={font=\color{black}},
xmin=-4173999,
xmax=-4161999,
xlabel={$R_x \left[-\right]$},
xmajorgrids,
every outer y axis line/.append style={black},
every y tick label/.append style={font=\color{black}},
ymin=0,
ymax=58.86,
ylabel={$a \left[\frac{m}{s^2}\right]$},
ymajorgrids,
axis x line*=bottom,
axis y line*=left
]
\node[right, align=left, inner sep=0mm, text=black]
at (axis cs:-4170999,31.392,0) {$\leftarrow$};
\node[left, align=right, inner sep=0mm, text=black]
at (axis cs:-4164999,31.392,0) {$\rightarrow$};
\node[align=center, inner sep=0mm, text=black]
at (axis cs:-4167999,31.392,0) { $6000\, \left[m\right]$};
\addplot [color=mycolor1,solid,mark=o,mark options={solid},forget plot]
table[row sep=crcr]{%
-4164999 28.90769598\\
-4165499 28.87376319\\
-4165999 28.83987945\\
-4166499 28.80605457\\
-4166999 28.77226893\\
-4167499 28.73854215\\
-4167999 28.70530587\\
-4168499 28.67253066\\
-4168999 28.6398045\\
-4169499 28.60712739\\
-4169999 28.57448952\\
-4170499 28.5419007\\
-4170999 28.50936093\\
};
\end{axis}
\end{tikzpicture}%
2 SOFTWARE/orbital_model/bach_beun.m
View file
@@ -0,0 +1,2 @@
orbitcheck1
orbitcheck2
20 SOFTWARE/orbital_model/compute_orbit.m
View file
@@ -1,12 +1,18 @@
clc
clear all
close all

% load constants
constants

ry = 10* R_m; %[m]
v = 7000; %[m/s]
dt = 1;
CD = 0.8;
rx = -4133000;
CL = 0.2*CD;
tend = 3600 * 24 * 10;
v = 5000; %[m/s]
dt_init = 1;
dt_atmos = 1;
dt_kep_init = 1e-8;
CD = 1.25;
rx = -4682000.0;
CL = 0.56;
tend = 3600 * 3;

[out] = orbit_full(rx,ry,CD,v,dt,CL,tend);
[out] = orbit_full(rx,ry,CD,v,dt_init,dt_atmos,dt_kep_init,CL,tend);
2 SOFTWARE/orbital_model/constants.m
View file
@@ -6,7 +6,7 @@
d = 12; % [m] diameter heatshield
S = 12^2*pi/4; %[m^2]
R_m = 6794000/2; %[m]
h_atm = 104e3; %[m]
h_atm = 400e3; %[m]
g_earth = 9.81; %[m/s^2]
crash_margin = 5000; %[m]
M_mars = 6.419*10^23; %[kg]
56 SOFTWARE/orbital_model/k_orbit_param.m
View file
@@ -0,0 +1,56 @@
function [kepler_param] = k_orbit_param(R1,R2,V,dt,dtheta,G,M)

format long
R1
R2
V
dt
dtheta
G
M

%% Determine kepler orbit parameters
mu = G*M;
% get magintudes
r = norm(R1)
v = norm(V);
r2 = norm(R2)

% calculate the semi major axis with visviva
a = r*mu / (2*mu - r * v^2);
disp(['visviva a: ' num2str(a)]);

% determine e with areal velocity
dA = 1/2 * r * r2 * sin(dtheta);
e = sqrt( 1 - 4 * (dA / dt)^2 / (a * mu) );

% calculate theta
theta = acos( (a*(1-e^2)-r) / (e*r) );

% calculate the orbital period
T = 2*pi * sqrt(a^3 / mu);

% semi-minor axis: b
b = a*sqrt(1-e^2);

% perifocal distance
rp = a*(1-e);

% apofocal distance
ra = 2*a - rp;

%
theta0 = atan2(R1(2),R1(1));
thetap = theta0 - theta;

% output
kepler_param.e = e; %
kepler_param.a = a; %
kepler_param.b = b; %
kepler_param.theta = theta; %
kepler_param.thetap = thetap; %
kepler_param.rp = rp; %
kepler_param.ra = ra; %
kepler_param.T = T; %

end
60 SOFTWARE/orbital_model/make_n_CD_plot.m
View file
@@ -0,0 +1,60 @@
clc
clear all
close all

constants
addpath('..\matlab2tikz')

fid = fopen('orbit_selection_const_CL.txt','r');
C = textscan(fid,'%f %f %f %f %d %d %d %f');
fclose(fid);

V = C{1,4};
CL = C{1,3};
CD = C{1,2};
max_a = C{1,8};
inorbit = C{1,7};

inorbit_index = find( inorbit == 1 );
V_inorbit = V(inorbit_index);
a_max_inorbit = max_a(inorbit_index);
CD_inorbit = CD(inorbit_index);

different_V = unique(V_inorbit);
different_V = sort(different_V);

figure('name','Different orbits for constant CD')
hold on
grid on
xlabel('$C_D \left[-\right]$','interpreter','LaTeX','fontsize',15)
ylabel('$a_{min} \left[\frac{m}{s^2}\right]$','interpreter','LaTeX','fontsize',18)
%ylim([0 6]*g_earth)
cc = parula(length(different_V) + 1);
legend_str = cell(length(different_V)+1,1);
marker = {'-+'; '-o'; '-*'; '-x'; '-s'; '-d'; '-.'; '-^'; '-v'; '->'; '-<'; '-p'; '-h'};
for i = 1:length(different_V)
index_V = find( V_inorbit == different_V(i) );
CD_V = CD_inorbit(index_V);
[CD_V, index_CD_order] = sort(CD_V);
a_max_V = a_max_inorbit(index_V);
a_max_V = a_max_V(index_CD_order);

% remove non min_a terms
different_CD = unique(CD_V);
a_plot = zeros(length(different_CD),1);

for k = 1:length(different_CD)

index_CD = find(CD_V == different_CD(k));
a_max_CD = a_max_V(index_CD);
a_plot(k) = min(a_max_CD)*g_earth;

end

legend_str{i} = [num2str(different_V(i)) ' m/s'];
plot(different_CD,a_plot,marker{i},'color',cc(i,:));
end
plot(xlim,[3,3]*g_earth,'-.','color','r');
legend_str{end} = '29.43 m/s^2 (3g)';
legend(legend_str,'location','northoutside');
matlab2tikz('.\LaTeX\n_CD.tikz','height','\figureheight','width','\figurewidth','showInfo', false,'checkForUpdates',false);
60 SOFTWARE/orbital_model/make_n_CL_plot.m
View file
@@ -0,0 +1,60 @@
clc
clear all
close all

constants
addpath('..\matlab2tikz')

fid = fopen('orbit_selection_const_CD_new.txt','r');
C = textscan(fid,'%f %f %f %f %d %d %d %f');
fclose(fid);

V = C{1,4};
CL = C{1,3};
CD = C{1,2};
max_a = C{1,8};
inorbit = C{1,7};

inorbit_index = find( inorbit == 1 );
V_inorbit = V(inorbit_index);
a_max_inorbit = max_a(inorbit_index);
CL_inorbit = CL(inorbit_index);

different_V = unique(V_inorbit);
different_V = sort(different_V);

figure('name','Different orbits for constant CD')
hold on
grid on
xlabel('$C_L \left[-\right]$','interpreter','LaTeX','fontsize',15)
ylabel('$a_{min} \left[\frac{m}{s^2}\right]$','interpreter','LaTeX','fontsize',18)
ylim([0 6]*g_earth)
cc = parula(length(different_V) + 1);
legend_str = cell(length(different_V)+1,1);
marker = {'-+'; '-o'; '-*'; '-x'; '-s'; '-d'; '-.'; '-^'; '-v'; '->'; '-<'; '-p'; '-h'};
for i = 1:length(different_V)
index_V = find( V_inorbit == different_V(i) );
CL_V = CL_inorbit(index_V);
[CL_V, index_CL_order] = sort(CL_V);
a_max_V = a_max_inorbit(index_V);
a_max_V = a_max_V(index_CL_order);

% remove non min_a terms
different_CL = unique(CL_V);
a_plot = zeros(length(different_CL),1);

for k = 1:length(different_CL)

index_CL = find(CL_V == different_CL(k));
a_max_CL = a_max_V(index_CL);
a_plot(k) = min(a_max_CL)*g_earth;

end

legend_str{i} = [num2str(different_V(i)) ' m/s'];
plot(different_CL,a_plot,marker{i},'color',cc(i,:));
end
plot(xlim,[3,3]*g_earth,'-.','color','r');
legend_str{end} = '29.43 m/s^2 (3g)';
legend(legend_str,'location','southeast');
matlab2tikz('.\LaTeX\n_CL.tikz','height','\figureheight','width','\figurewidth','showInfo', false,'checkForUpdates',false);
44 SOFTWARE/orbital_model/make_range_plot.m
View file
@@ -0,0 +1,44 @@
% read file and plot
clc
clear all
close all

constants
addpath('..\matlab2tikz')

fid = fopen('orbit_selection_range.txt','r');
C = textscan(fid,'%f %f %f %f %d %d %d %f');
fclose(fid);

RX = C{1,1};
V = C{1,4};
CL = C{1,3};
CD = C{1,2};
max_a = C{1,8};
inorbit = C{1,7};

inorbit_index = find( inorbit == 1 );
V_inorbit = V(inorbit_index);
a_max_inorbit = max_a(inorbit_index);
CL_inorbit = CL(inorbit_index);
RX_inorbit = RX(inorbit_index);

distance = max(RX_inorbit) - min(RX_inorbit);


fig = figure('name','results');
hold on
grid on
v_text = 3.2*g_earth;
text(min(RX_inorbit),v_text,'$\leftarrow$','HorizontalAlignment','left','Interpreter','LaTeX','fontsize',15)
text(max(RX_inorbit),v_text,'$\rightarrow$','HorizontalAlignment','right','Interpreter','LaTeX','fontsize',15)
text(min(RX_inorbit) + distance/2,v_text,[' $' num2str(distance) '\, \left[m\right]$'],'HorizontalAlignment','center','Interpreter','LaTeX','fontsize',15)
ylim([0 6]*g_earth)
xlim([min(RX_inorbit)-distance/2 max(RX_inorbit)+distance/2])
cc = parula(5);
plot(RX_inorbit,a_max_inorbit*g_earth,'-o','color',cc(2,:))
xlabel('$R_x \left[-\right]$','interpreter','LaTeX','fontsize',15)
ylabel('$a \left[\frac{m}{s^2}\right]$','interpreter','LaTeX','fontsize',18)
matlab2tikz('.\LaTeX\range_plot.tikz','height','\figureheight','width','\figurewidth','showInfo', false,'checkForUpdates',false);


4 SOFTWARE/orbital_model/orbit.m
View file
@@ -18,8 +18,8 @@
orbit.al = cross(vel_unit,[0,0,1]) * CL * 0.5 * rho * norm(V - Vatm)^2 * S / m;
else
%no lift and drag outside the atmosphere
orbit.ad = 0;
orbit.al = 0;
orbit.ad = [0,0,0];
orbit.al = [0,0,0];
end
%calculate total accaleration
orbit.a = orbit.ag + orbit.ad + orbit.al;
50 SOFTWARE/orbital_model/orbit_full.m
View file
@@ -1,4 +1,4 @@
function [out] = orbit_full(rx,ry,CD,v,dt,CL,tend)
function [out] = orbit_full(rx,ry,CD,v,dt_init,dt_atmos,dt_kep_init,CL,tend);
% load constants
constants

@@ -25,21 +25,53 @@
i = 1;
% waitbar
h = waitbar(0,'Initializing waitbar...');

% kepler orbit
to_kepler = false;

% define time step
dt = dt_init;

while true

% get orbital parameters at next node
orbit_new = orbit(out.R(i,:),out.V(i,:),out.a(i,:),CD,CL,dt,atm,R_m,Omega_m,S,m);
if to_kepler

[orbit_new,t_kep] = orbit_kepler(kepler_param,orbit_new);
to_kepler = false;
out.inorbit = false;

else
% get orbital parameters at next node
orbit_new = orbit(out.R(i,:),out.V(i,:),out.a(i,:),CD,CL,dt,atm,R_m,Omega_m,S,m);
end

out.R(i+1,:) = orbit_new.R;
out.V(i+1,:) = orbit_new.V;
out.a(i+1,:) = orbit_new.a;
out.ad(i+1,:) = orbit_new.ad;
out.al(i+1,:) = orbit_new.al;
out.ag(i+1,:) = orbit_new.ag;

if out.inorbit && (to_kepler == false)

orbit_new = orbit(out.R(i,:),out.V(i,:),out.a(i,:),CD,CL,dt_kep_init,atm,R_m,Omega_m,S,m);

dtheta = atan2(orbit_new.R(2),orbit_new.R(1)) - atan2(out.R(i,2),out.R(i,1));
kepler_param = k_orbit_param(out.R(i,:),orbit_new.R,out.V(i,:),dt_kep_init,dtheta,G,M_mars);
% BEUN stop om te controleren of script klopt
disp(kepler_param);
to_kepler = true;
out.inatmos = false;
% break;
end



% if the s/c has ever been in the atmosphere in_atmus => true
% change time step to time step in atmos!
if (norm(out.R(i+1,:)) < (R_m + h_atm))
out.inatmos = true;
dt = dt_atmos;
end

% if s/c has been in atmos and is leaving atmos again (at next node)
@@ -104,14 +136,24 @@
% circle plot:
theta_plot = 0:0.01:2*pi;
radius_mars = ones(1,length(theta_plot)) * R_m;
radius_mars_atmos = ones(1,length(theta_plot)) * (R_m + 150000);
radius_mars_atmos = ones(1,length(theta_plot)) * (R_m + h_atm);
figure('name','Orbit')
grid on
axis equal
hold on
plot(out.R(:,1),out.R(:,2))
polar(theta_plot,radius_mars,'r');
polar(theta_plot,radius_mars_atmos,'g')
% plot kepler orbit
if exist('kepler_param','var')

rk = kepler_param.a * (1- kepler_param.e^2) ./ (1 + kepler_param.e .* cos(theta_plot));
polar(theta_plot+kepler_param.thetap,rk,'k');
plot(kepler_param.rp*cos(kepler_param.thetap),kepler_param.rp*sin(kepler_param.thetap),'*')
plot(-kepler_param.ra*cos(kepler_param.thetap),-kepler_param.ra*sin(kepler_param.thetap),'d')
plot((R_m + h_atm)*cos(kepler_param.thetap+kepler_param.theta),(R_m + h_atm)*sin(kepler_param.thetap+kepler_param.theta),'p')
plot((R_m + h_atm)*cos(kepler_param.thetap-kepler_param.theta),(R_m + h_atm)*sin(kepler_param.thetap-kepler_param.theta),'p')
end


t = 0:dt:(length(out.R)*dt - dt);
61 SOFTWARE/orbital_model/orbit_kepler.m
View file
@@ -0,0 +1,61 @@
function [orbit,t] = orbit_kepler(kepler_param,orbit_init)

a = kepler_param.a;
e = kepler_param.e;
b = kepler_param.b;
theta = kepler_param.theta;
thetap = kepler_param.thetap;
T = kepler_param.T;

r = norm(orbit_init.R);

% determine time at reentry
% area of ellipse between exit atmos and going back:
c = a*e + r * cos(theta)
Atot = b*a*pi;
A = 1/2*b*a*pi + b*c*sqrt(1-c^2/a^2) - b*a*asin(-c/a) - 1/2*r^2*cos(theta)*sin(theta);

t = A/Atot * T;

% Determine the R, V, and accelerations at the end point:
angle = 2*theta;
Tz = [ cos(angle), sin(angle), 0;...
-sin(angle), cos(angle), 0;...
0, 0, 1];

% rotation and translation matrix for V

% convert to new point..
Ve = Tz * orbit_init.V';
% rotate to another axis system:
angle = thetap;
Tv = [ cos(angle), sin(angle), 0;...
-sin(angle), cos(angle), 0;...
0, 0, 1];
Ve = Tv * Ve;
%reflect around ye axis
Refy = [1, 0, 0;...
0, -1, 0;...
0, 0, 1];
Ve = Refy * Ve;
% rotate back to express in 0 frame
Tv = [ cos(angle), -sin(angle), 0;...
sin(angle), cos(angle), 0;...
0, 0, 1];
Ve = Tv * Ve;
R = (Tz * orbit_init.R')'
V = Ve'
acc = (Tz * orbit_init.a')'
ag = (Tz * orbit_init.ag')'



% output
orbit.R = R;
orbit.V = V;
orbit.a = acc;
orbit.ad = [0,0,0];
orbit.al = [0,0,0];
orbit.ag = ag;

end
9,980 SOFTWARE/orbital_model/orbit_selection_const_CD.txt
View file

Large diffs are not rendered by default.

14,088 SOFTWARE/orbital_model/orbit_selection_const_CD_new.txt
View file

Large diffs are not rendered by default.

2,057 SOFTWARE/orbital_model/orbit_selection_const_CL.txt
View file

Large diffs are not rendered by default.

22 SOFTWARE/orbital_model/orbit_selection_range.txt
View file
@@ -0,0 +1,22 @@
-4161499.0 1.25 -0.40 7000.0 1 1 0 5.515969
-4161999.0 1.25 -0.40 7000.0 1 0 0 2.657648
-4161999.000000 1.250000 -0.400000 7000.000000 1 1 0 5.383808
-4162499.000000 1.250000 -0.400000 7000.000000 1 1 0 5.237833
-4162999.000000 1.250000 -0.400000 7000.000000 1 1 0 5.068766
-4163499.000000 1.250000 -0.400000 7000.000000 1 1 0 4.870082
-4163999.000000 1.250000 -0.400000 7000.000000 1 1 0 4.611254
-4164499.000000 1.250000 -0.400000 7000.000000 1 1 0 4.178437
-4164999.000000 1.250000 -0.400000 7000.000000 1 0 1 2.946758
-4165499.000000 1.250000 -0.400000 7000.000000 1 0 1 2.943299
-4165999.000000 1.250000 -0.400000 7000.000000 1 0 1 2.939845
-4166499.000000 1.250000 -0.400000 7000.000000 1 0 1 2.936397
-4166999.000000 1.250000 -0.400000 7000.000000 1 0 1 2.932953
-4167499.000000 1.250000 -0.400000 7000.000000 1 0 1 2.929515
-4167999.000000 1.250000 -0.400000 7000.000000 1 0 1 2.926127
-4168499.000000 1.250000 -0.400000 7000.000000 1 0 1 2.922786
-4168999.000000 1.250000 -0.400000 7000.000000 1 0 1 2.919450
-4169499.000000 1.250000 -0.400000 7000.000000 1 0 1 2.916119
-4169999.000000 1.250000 -0.400000 7000.000000 1 0 1 2.912792
-4170499.000000 1.250000 -0.400000 7000.000000 1 0 1 2.909470
-4170999.000000 1.250000 -0.400000 7000.000000 1 0 1 2.906153
-4171499.000000 1.250000 -0.400000 7000.000000 1 0 0 2.902840
125 SOFTWARE/orbital_model/orbitcheck.m
View file
@@ -7,22 +7,22 @@

do_plot = false;
write_to_file = true;
file_name = 'orbit_selection_range.txt';


ry = 10* R_m; %[m]
v = 7000; %[m/s]
dt = 1;

CLCD = 0.2; %[-]
refinement_steps = 40;
refinement_steps = 100;


% changing variables
%
rx = -4.05e6:-1e3:-4.2e6;
CD = 0.4:0.1:0.7;
% CD = 1.3;
% rx = -4.162e6;
rx = -4161499:-5e2:-4.8e6;
CD = 1.25;
%CLCD = [-0.25:0.1:0 0:0.05:0.3]; %[-]
CLCD = 0.25;
CL = -0.4;

cc = parula(length(CD)+3);
if do_plot
@@ -31,66 +31,69 @@

rx_first = rx(1);
% file string
filestr = cell(length(CD),1);
parfor k = 1:length(CD)
CL = CLCD * CD(k);
% store first value of crashed..
rx_crashed = rx_first;
for i=1:length(rx)
[out] = orbit_selection(rx(i),ry,CD(k),v,dt,CL,R_m,Omega_m,S,m,G,M_mars,h_atm,crash_margin,g_earth);

if write_to_file
filestr{k} = [filestr{k} sprintf(' %11.1f %4.2f %d %d %d %f \n',rx(i),CD(k),out.inatmos,out.crash,out.inorbit,out.maxaccel)];
end

% store last rx value if crashed
if (out.crash) || (out.inorbit)
rx_crashed = rx(i);
end

% check if not crashed and not in orbit
if (out.crash == false) && (out.inorbit == false)
previous_inorbit = false;
rx_refine = linspace(rx_crashed,rx(i),refinement_steps+2);
for j = 2:(length(rx_refine)-1)
[out] = orbit_selection(rx_refine(j),ry,CD(k),v,dt,CL,R_m,Omega_m,S,m,G,M_mars,h_atm,crash_margin,g_earth);
if write_to_file
filestr{k} = [filestr{k} sprintf(' %11.1f %4.2f %d %d %d %f \n',rx(i),CD(k),out.inatmos,out.crash,out.inorbit,out.maxaccel)];
end

if previous_inorbit && (out.inorbit == 0)
break;
end

previous_inorbit = out.inorbit;

if do_plot && (out.crash == false)
t = 0:dt:(length(out.R)*dt-dt);
subplot(3,1,1)
hold on
grid on
plot(t,sqrt(out.R(:,1).^2 + out.R(:,2).^2 + out.R(:,3).^2),'color',cc(k,:));
subplot(3,1,2)
hold on
grid on
plot(t,sqrt(out.V(:,1).^2 + out.V(:,2).^2 + out.V(:,3).^2),'color',cc(k,:));
subplot(3,1,3)
hold on
grid on
plot(t,sqrt(out.a(:,1).^2 + out.a(:,2).^2 + out.a(:,3).^2),'color',cc(k,:));
filestr = cell(length(v),1);
parfor u = 1:length(v)
for l = 1:length(CL)
for k = 1:length(CD)

% store first value of crashed..
rx_crashed = rx_first;
for i=1:length(rx)
[out] = orbit_selection(rx(i),ry,CD(k),v(u),dt,CL(l),R_m,Omega_m,S,m,G,M_mars,h_atm,crash_margin,g_earth);

if write_to_file
filestr{u} = [filestr{u} sprintf(' %11.1f %4.2f %4.2f %6.1f %d %d %d %f \n',rx(i),CD(k),CL(l),v(u),out.inatmos,out.crash,out.inorbit,out.maxaccel)];
end

% store last rx value if crashed
if (out.crash) || (out.inorbit)
rx_crashed = rx(i);
end

% check if not crashed and not in orbit
if (out.crash == false) && (out.inorbit == false)
previous_inorbit = false;
rx_refine = linspace(rx_crashed,rx(i),refinement_steps+2);
for j = 2:(length(rx_refine)-1)
[out] = orbit_selection(rx_refine(j),ry,CD(k),v(u),dt,CL(l),R_m,Omega_m,S,m,G,M_mars,h_atm,crash_margin,g_earth);
if write_to_file
filestr{u} = [filestr{u} sprintf(' %f %f %f %f %d %d %d %f \n',rx(j),CD(k),CL(l),v(u),out.inatmos,out.crash,out.inorbit,out.maxaccel)];
end

if previous_inorbit && (out.inorbit == 0)
break;
end

previous_inorbit = out.inorbit;

if do_plot && (out.crash == false)
t = 0:dt:(length(out.R)*dt-dt);
subplot(3,1,1)
hold on
grid on
plot(t,sqrt(out.R(:,1).^2 + out.R(:,2).^2 + out.R(:,3).^2),'color',cc(k,:));
subplot(3,1,2)
hold on
grid on
plot(t,sqrt(out.V(:,1).^2 + out.V(:,2).^2 + out.V(:,3).^2),'color',cc(k,:));
subplot(3,1,3)
hold on
grid on
plot(t,sqrt(out.a(:,1).^2 + out.a(:,2).^2 + out.a(:,3).^2),'color',cc(k,:));
end
end
break;
end
break;



end



end
end
end

if write_to_file
fid = fopen('orbit_selection.txt','a');
for k = 1:length(CD)
fid = fopen(file_name,'a');
for k = 1:length(v)
fprintf(fid,'%s',filestr{k});
end
fclose(fid);
11 SOFTWARE/orbital_model/read_outputfile.m
View file
@@ -3,18 +3,19 @@
clear all
close all

fid = fopen('orbit_selection.txt','r');
C = textscan(fid,'%11.1f %4.2f %d %d %d %f');
fid = fopen('orbit_selection_range.txt','r');
C = textscan(fid,'%f %f %f %f %d %d %d %f');
fclose(fid);

%rx(i),CD(k),CL(l),v(u),out.inatmos,out.crash,out.inorbit,out.maxaccel)];

fig = figure('name','results');
hold on
grid on
inorbit = C{5};
inorbit = C{7};
rx = C{1};
CD = C{2};
a = C{6};
a = C{8};

k = 1;
for i=1:length(inorbit)
@@ -29,6 +30,6 @@

end
[Xm,Ym] = meshgrid(rx_orb,CD_orb);
[X, Y,Z] = griddata(rx_orb,CD_orb,a_orb,Xm,Ym);
[X, Y, Z] = griddata(rx_orb,CD_orb,a_orb,Xm,Ym);
plot3(X,Y,Z,'o')
%colormap(fig,jet)
16 SOFTWARE/orbital_model/testbeunacceleratie.m
View file

This file was deleted.