diff --git a/AUTHORS.md b/AUTHORS.md
new file mode 100644
index 0000000000000000000000000000000000000000..505a24659e59ab83bc0a6fdf5bd227b724dc73e1
--- /dev/null
+++ b/AUTHORS.md
@@ -0,0 +1,14 @@
+# Skyward Experimental Rocketry
+website: http://www.skywarder.eu
+
+## Ruben Di Battista
+email: ruben.dibattista@skywarder.eu
+
+## Francesco Colombi
+email: francesco.colombi@skywarder.eu
+
+## Gabriele Poiana
+email: gabriele.poiana@skywarder.eu
+
+## Salvatore Andrea Bella
+email: salvatore.andrea.bella@skywarder.eu
\ No newline at end of file
diff --git a/LICENSE/LICENSE b/LICENSE/LICENSE
new file mode 100644
index 0000000000000000000000000000000000000000..c9663538595f803b0ae7d52c4c959bbb9d683d65
--- /dev/null
+++ b/LICENSE/LICENSE
@@ -0,0 +1,11 @@
+Copyright (c) 2014, Skyward Experimental Rocketry <www.skywarder.eu>
+Author: Ruben Di Battista <ruben.dibattista@skywarder.eu> | CRD Department <crd@skywarder.eu>
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
+
+1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
+
+2. 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 HOLDER 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.
diff --git a/LICENSE/run.txt b/LICENSE/run.txt
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/data/R1X_empty.mat b/data/R1X_empty.mat
new file mode 100644
index 0000000000000000000000000000000000000000..e7cdab6b5f1a1d2b2dac0238b4e06ab5b2ba85aa
--- /dev/null
+++ b/data/R1X_empty.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:04d546630ff2edcd4e5a80135510f684369b2f5afc1161a922c2ff643a8d57db
+size 180528
diff --git a/data/R1X_full.mat b/data/R1X_full.mat
new file mode 100644
index 0000000000000000000000000000000000000000..6f97ea490a7aa94de1f458380561052abdaec5aa
--- /dev/null
+++ b/data/R1X_full.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:245f2d7e21bfcf139b47b40a4d700131fe600c861031af8c6448a1d7d1935955
+size 180528
diff --git a/data/R2A_empty.mat b/data/R2A_empty.mat
new file mode 100644
index 0000000000000000000000000000000000000000..45806efcab6491bf095005e29a19d293a4b78a90
--- /dev/null
+++ b/data/R2A_empty.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:aca8fe82cda534ed2e378eb1a026a983391218882ace088874653fbc9b27e3d0
+size 1583632
diff --git a/data/R2A_full.mat b/data/R2A_full.mat
new file mode 100644
index 0000000000000000000000000000000000000000..cee99450e79c57988cf610e6901a454135d070a6
--- /dev/null
+++ b/data/R2A_full.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:582428d14fc90a340b2d9e361b4225b4f25b9b9743c16328426ee4bd40dbd23d
+size 1583632
diff --git a/scratch/CLN.m b/scratch/CLN.m
new file mode 100644
index 0000000000000000000000000000000000000000..c93a9ca61ec2bdbb03bfbc21cb49c9857c2921a6
--- /dev/null
+++ b/scratch/CLN.m
@@ -0,0 +1,51 @@
+close all
+clear, clc
+
+% AERODYNAMICS DETAILS %
+% This coefficients are intended to be obtained through MISSILE DATCOM
+% (than parsed with the tool datcom_parser.py)
+CoeffsF = load('for006_full.mat','Coeffs');
+settings.CoeffsF = CoeffsF.Coeffs;
+clear('CoeffsF');
+
+%Note: All the parameters (AoA,Betas,Altitudes,Machs) must be the same for
+%empty and full configuration
+s = load('for006_full.mat','State');
+settings.Alphas = s.State.Alphas';
+settings.Betas = s.State.Betas';
+settings.Altitudes = s.State.Altitudes';
+settings.Machs = s.State.Machs';
+clear('s');
+
+% Coeffs. Interpolation
+givA = settings.Alphas*pi/180;
+givB = settings.Betas*pi/180;
+givH = settings.Altitudes;
+givM = settings.Machs;
+
+alpha = 0;
+M = 1.7;
+beta = [-5:0.1:5]*pi/180;
+z = 0;
+
+Cln = zeros(length(beta),2);
+Cnb = zeros(length(beta),2);
+Cnr = zeros(length(beta),2);
+Cnp = zeros(length(beta),2);
+
+for k = 1:length(beta)
+ Cln(k,1)=interp4_easy(givA,givM,givB,givH,settings.CoeffsF.CLN,alpha,M,beta(k),-z);%,'nearest');
+ Cnb(k,1)=interp4_easy(givA,givM,givB,givH,settings.CoeffsF.CLNB,alpha,M,beta(k),-z);%,'nearest');
+ Cnr(k,1)=interp4_easy(givA,givM,givB,givH,settings.CoeffsF.CLNR,alpha,M,beta(k),-z);%,'nearest');
+ Cnp(k,1)=interp4_easy(givA,givM,givB,givH,settings.CoeffsF.CLNP,alpha,M,beta(k),-z);%,'nearest');
+
+ Cln(k,2)=interpn(givA,givM,givB,givH,settings.CoeffsF.CLN,alpha,M,beta(k),-z);%,'nearest');
+ Cnb(k,2)=interpn(givA,givM,givB,givH,settings.CoeffsF.CLNB,alpha,M,beta(k),-z);%,'nearest');
+ Cnr(k,2)=interpn(givA,givM,givB,givH,settings.CoeffsF.CLNR,alpha,M,beta(k),-z);%,'nearest');
+ Cnp(k,2)=interpn(givA,givM,givB,givH,settings.CoeffsF.CLNP,alpha,M,beta(k),-z);%,'nearest');
+end
+
+figure(), plot(beta*180/pi, Cln), title('CLN'), grid on;
+figure(), plot(beta*180/pi, Cnb), title('CLN/B'), grid on;
+figure(), plot(beta*180/pi, Cnr), title('CLN/r'), grid on;
+figure(), plot(beta*180/pi, Cnp), title('CLN/p'), grid on;
\ No newline at end of file
diff --git a/scratch/THp.m b/scratch/THp.m
new file mode 100644
index 0000000000000000000000000000000000000000..616769783dd26bc2970e450ad138180c57190261
--- /dev/null
+++ b/scratch/THp.m
@@ -0,0 +1,38 @@
+% Author: Francesco Colombi
+% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu
+% email: francesco.colombi@skywarder.eu
+% Release date: 16/04/2016
+
+clear, clc
+
+% %Cesaroni PRO 150 White Thunder
+% sTH = [0 0.05 0.15 0.5 0.6 0.74 0.85 1.15 1.35 2.45 ...
+% 3 3.7 4 4.5 4.8 4.9 5 5.05 5.1 5.15 5.2]; %s
+% THs = [0 9200 7900 8500 8500 8350 8300 8400 8500 8500 ...
+% 8300 8000 7800 7600 7450 7300 7300 4500 500 100 0]; %N
+
+%Cesaroni PRO 150 SkidMark
+sTH = [0 0.1 0.2 0.3 0.4 0.5 0.6 0.8 1.2 1.8 3.2 3.6 ...
+ 4.8 6 7 7.2 7.6 7.8 7.9 8 8.1 8.19]; %s
+THs = [0 3400 3100 3000 3300 3400 3500 3700 3700 3800 ...
+ 4000 4081.6 3900 3800 3700 3500 3350 3200 3000 2000 750 0]; %N
+
+
+%sTH = [0 0.74 2.56 4.37 5.2];
+%THs = [0 8500 8500 7600 0];
+%sTH = [0 0.15 0.59 1.28 2.115 3 3.84 4.52 4.97 5.12];
+%THs = [0 5500 6100 6100 6250 6050 5800 5200 5050 0];
+
+
+s=linspace(0,sTH(end));
+dt = s(2);
+TH = interp1(sTH,THs,s);
+
+sum = 0;
+for j = 1:length(s)
+ sum = sum + TH(j)*dt;
+end
+media = sum/5.12 %N
+
+plot(s,TH);
+grid on;
\ No newline at end of file
diff --git a/scratch/ascend_constThrust.m b/scratch/ascend_constThrust.m
new file mode 100644
index 0000000000000000000000000000000000000000..306757b24155de038dfff0f3967e59372d5984d3
--- /dev/null
+++ b/scratch/ascend_constThrust.m
@@ -0,0 +1,454 @@
+function [ dY ] = ascend( t,Y,settings,uw,vw,ww )
+% ODE-Function of the 6DOF Rigid Rocket Model
+% State = ( x y z | u v w | p q r | q0 q1 q2 q3 )
+%
+% (x y z): NED Earth's Surface Centered Frame ("Inertial") coordinates
+% (u v w): body frame velocities
+% (p q r): body frame angular rates
+% (q0 q1 q2 q3): attitude unit quaternion
+%
+%
+% NOTE: To get the NED velocities the body-frame must be multiplied for the
+% conjugated of the current attitude quaternion
+% E.G.
+%
+%
+% quatrotate(quatconj(Y(:,10:13)),Y(:,4:6))
+
+% Author: Ruben Di Battista
+% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu
+% email: ruben.dibattista@skywarder.eu
+% Website: http://www.skywarder.eu
+% April 2014; Last revision: 31.XII.2014
+% License: 2-clause BSD
+
+% Author: Francesco Colombi
+% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu
+% email: francesco.colombi@skywarder.eu
+% Release date: 16/04/2016
+
+x = Y(1);
+y = Y(2);
+z = Y(3);
+u = Y(4);
+v = Y(5);
+w = Y(6);
+p = Y(7);
+q = Y(8);
+r = Y(9);
+q0 = Y(10);
+q1 = Y(11);
+q2 = Y(12);
+q3 = Y(13);
+m = Y(14);
+Ixx = Y(15);
+Iyy = Y(16);
+Izz = Y(17);
+
+global bool
+% if bool = 0 the trends during the integration are NOT requested
+% if bool = 1 the trends during the integration are saved and plotted
+global t_plot contatore
+if bool == 1
+ t_plot(contatore) = t;
+end
+
+
+Q = [ q0 q1 q2 q3];
+Q_conj= [ q0 -q1 -q2 -q3];
+normQ=norm(Q);
+
+if abs(normQ-1)>0.1
+ Q = Q/normQ;
+ %frprintf('Norm of the quaternion is (too much!) less than one (%3.2f)\n', ...
+ %normQ);
+end
+
+
+% Adding Wind (supposed to be added in NED axes);
+
+% constant wind
+wind = quatrotate(Q, [uw vw ww]);
+
+% % Wind model
+% h = 0;
+% if -z > 0
+% h = -z+settings.z0;
+% end
+% % global alt
+% % if bool == 1
+% % alt = [alt; h];
+% % end
+%
+% % wind in NED axis
+% [uw,vw] = atmoshwm07(settings.wind.Lat, settings.wind.Long, h, ...
+% 'day',settings.wind.Day,'model','quiet');
+%
+% % global WIND
+% % if bool == 1
+% % WIND = [WIND; uw, vw, ww];
+% % end
+%
+% % wind in body frame
+% wind = quatrotate(Q, [uw, vw, ww]);
+
+
+% % Wind Model
+% [uw,vw,ww] = wind_generator(settings,z,t,Q);
+% wind = [ uw,vw,ww ];
+
+% Relative velocities (plus wind);
+ur = u - wind(1);
+vr = v - wind(2);
+wr = w - wind(3);
+
+% Body to Inertial velocities
+Vels = quatrotate(Q_conj,[u v w]);
+
+V_norm = norm([ur vr wr]);
+
+
+
+% Constants
+S = settings.S;
+C = settings.C;
+CoeffsE = settings.CoeffsE; %Empty Rocket Coefficients
+CoeffsF = settings.CoeffsF; % Full Rocket Coefficients
+
+g = 9.81;
+%Time-Dependant Properties
+m0 = settings.m0; %kg mass
+mfr = settings.mfr;
+tb = settings.tb;
+
+Ixxf= settings.Ixxf;
+Ixxe= settings.Ixxe;
+Iyyf= settings.Iyyf;
+Iyye= settings.Iyye;
+Izzf= settings.Izzf;
+Izze= settings.Izze;
+
+dI = 1/tb*([Ixxf Iyyf Izzf]'-[Ixxe Iyye Izze]');
+
+if t<tb
+ mdot = -mfr;
+ Ixxdot = -dI(1);
+ Iyydot = -dI(2);
+ Izzdot = -dI(3);
+ T = settings.T_coeffs*[t^0 t^1 t^2 t^3 t^4]';
+
+
+else
+ mdot = 0;
+ Ixxdot = 0;
+ Iyydot = 0;
+ Izzdot = 0;
+ T=0;
+end
+
+global T_plot
+if bool == 1
+ T_plot(contatore) = T; % thrust
+end
+
+%Atmosphere
+if -z<0
+ z = 0;
+end
+[~, a, ~, rho] = atmoscoesa(-z+settings.z0);
+M = V_norm/a;
+
+global M_plot
+if bool == 1
+ M_plot(contatore) = M; % mach
+end
+
+%Aero Angles
+if not(u<1e-1 || V_norm<1e-3);
+ alpha = atan(wr/ur);
+ beta = asin(vr/V_norm);
+else
+ alpha =0;
+ beta = 0;
+end
+
+global alpha_plot beta_plot
+if bool == 1
+ alpha_plot(contatore) = alpha;
+ beta_plot(contatore) = beta;
+end
+
+% beta = 0; % prova simulazione con beta imposto uguale a 0
+
+
+
+% Coeffs. Interpolation
+givA = settings.Alphas*pi/180;
+givB = settings.Betas*pi/180;
+givH = settings.Altitudes;
+givM = settings.Machs;
+
+% Interpolation error check
+
+if M > givM(end)
+ %frprintf('Mach No. is more than the max provided (M = %3.2f @ t=%2.2f)\n\n', ...
+ %M,t);
+ M = givM(end);
+
+end
+
+if M<givM(1)
+ %frprintf('Mach No. is less than the min provided (M = %3.2f @ t=%2.2f)\n\n', ...
+ % M,t);
+ M = givM(1);
+
+end
+
+if alpha > givA(end)
+ %frprintf('AoA is more than the max provided (alpha = %3.2f @ t=%2.2f)\n\n', ...
+ %alpha*180/pi,t)
+ alpha= givA(end);
+
+else if alpha < givA(1)
+ %frprintf('AoA is less than the min provided (alpha = %3.2f @ t=%2.2f)\n\n',...
+ %alpha*180/pi,t);
+ alpha = givA(1);
+
+
+ end
+end
+
+if beta > givB(end)
+ beta= givB(end);
+else if beta < givB(1)
+ beta = givB(1);
+ end
+end
+
+if -z >givH(end)
+ z = -givH(end);
+else if -z < givH(1)
+ z = -givH(1);
+ end
+end
+
+% interpolation with the value in the nearest condition: interp_easy
+% full
+CAf=interp4_easy(givA,givM,givB,givH,CoeffsF.CA,alpha,M,beta,-z);%,'nearest');
+CYBf=interp4_easy(givA,givM,givB,givH,CoeffsF.CYB,alpha,M,beta,-z);%,'nearest');
+CNAf=interp4_easy(givA,givM,givB,givH,CoeffsF.CNA,alpha,M,beta,-z);%,'nearest');
+Clf=interp4_easy(givA,givM,givB,givH,CoeffsF.CLL,alpha,M,beta,-z);%,'nearest');
+Clpf=interp4_easy(givA,givM,givB,givH,CoeffsF.CLLP,alpha,M,beta,-z);%,'nearest');
+Cmaf=interp4_easy(givA,givM,givB,givH,CoeffsF.CMA,alpha,M,beta,-z);%,'nearest');
+Cmadf=interp4_easy(givA,givM,givB,givH,CoeffsF.CMAD,alpha,M,beta,-z);%,'nearest');
+Cmqf=interp4_easy(givA,givM,givB,givH,CoeffsF.CMQ,alpha,M,beta,-z);%,'nearest');
+Cnbf=interp4_easy(givA,givM,givB,givH,CoeffsF.CLNB,alpha,M,beta,-z);%,'nearest');
+Cnrf=interp4_easy(givA,givM,givB,givH,CoeffsF.CLNR,alpha,M,beta,-z);%,'nearest');
+Cnpf=interp4_easy(givA,givM,givB,givH,CoeffsF.CLNP,alpha,M,beta,-z);%,'nearest');
+%
+Clrf=interp4_easy(givA,givM,givB,givH,CoeffsF.CLLR,alpha,M,beta,-z);%,'nearest');
+Clbf=interp4_easy(givA,givM,givB,givH,CoeffsF.CLLB,alpha,M,beta,-z);%,'nearest');
+CDf =interp4_easy(givA,givM,givB,givH,CoeffsF.CD,alpha,M,beta,-z);
+
+% empty
+CAe=interp4_easy(givA,givM,givB,givH,CoeffsE.CA,alpha,M,beta,-z);%,'nearest');
+CYBe=interp4_easy(givA,givM,givB,givH,CoeffsE.CYB,alpha,M,beta,-z);%,'nearest');
+CNAe=interp4_easy(givA,givM,givB,givH,CoeffsE.CNA,alpha,M,beta,-z);%,'nearest');
+Cle=interp4_easy(givA,givM,givB,givH,CoeffsE.CLL,alpha,M,beta,-z);%,'nearest');
+Clpe=interp4_easy(givA,givM,givB,givH,CoeffsE.CLLP,alpha,M,beta,-z);%,'nearest');
+Cmae=interp4_easy(givA,givM,givB,givH,CoeffsE.CMA,alpha,M,beta,-z);%,'nearest');
+Cmade=interp4_easy(givA,givM,givB,givH,CoeffsE.CMAD,alpha,M,beta,-z);%,'nearest');
+Cmqe=interp4_easy(givA,givM,givB,givH,CoeffsE.CMQ,alpha,M,beta,-z);%,'nearest');
+Cnbe=interp4_easy(givA,givM,givB,givH,CoeffsE.CLNB,alpha,M,beta,-z);%,'nearest');
+Cnre=interp4_easy(givA,givM,givB,givH,CoeffsE.CLNR,alpha,M,beta,-z);%,'nearest');
+Cnpe=interp4_easy(givA,givM,givB,givH,CoeffsE.CLNP,alpha,M,beta,-z);%,'nearest');
+%
+Clre=interp4_easy(givA,givM,givB,givH,CoeffsE.CLLR,alpha,M,beta,-z);%,'nearest');
+Clbe=interp4_easy(givA,givM,givB,givH,CoeffsE.CLLB,alpha,M,beta,-z);%,'nearest');
+CDe =interp4_easy(givA,givM,givB,givH,CoeffsE.CD,alpha,M,beta,-z);
+
+% % linear interpolation of the coeff
+% %full
+% CAf = interpn(givA,givM,givB,givH,CoeffsF.CA,alpha,M,beta,-z);
+% CYBf = interpn(givA,givM,givB,givH,CoeffsF.CYB,alpha,M,beta,-z);
+% CNAf = interpn(givA,givM,givB,givH,CoeffsF.CNA,alpha,M,beta,-z);
+% Clf = interpn(givA,givM,givB,givH,CoeffsF.CLL,alpha,M,beta,-z);
+% Clpf = interpn(givA,givM,givB,givH,CoeffsF.CLLP,alpha,M,beta,-z);
+% Cmaf = interpn(givA,givM,givB,givH,CoeffsF.CMA,alpha,M,beta,-z);
+% Cmadf = interpn(givA,givM,givB,givH,CoeffsF.CMAD,alpha,M,beta,-z);
+% Cmqf = interpn(givA,givM,givB,givH,CoeffsF.CMQ,alpha,M,beta,-z);
+% Cnbf = interpn(givA,givM,givB,givH,CoeffsF.CLNB,alpha,M,beta,-z);
+% Cnrf = interpn(givA,givM,givB,givH,CoeffsF.CLNR,alpha,M,beta,-z);
+% Cnpf = interpn(givA,givM,givB,givH,CoeffsF.CLNP,alpha,M,beta,-z);
+% %
+% Clrf = interpn(givA,givM,givB,givH,CoeffsF.CLLR,alpha,M,beta,-z);
+% Clbf = interpn(givA,givM,givB,givH,CoeffsF.CLLB,alpha,M,beta,-z);
+% CDf = interpn(givA,givM,givB,givH,CoeffsF.CD,alpha,M,beta,-z);
+%
+% % empty
+% CAe = interpn(givA,givM,givB,givH,CoeffsE.CA,alpha,M,beta,-z);
+% CYBe = interpn(givA,givM,givB,givH,CoeffsE.CYB,alpha,M,beta,-z);
+% CNAe = interpn(givA,givM,givB,givH,CoeffsE.CNA,alpha,M,beta,-z);
+% Cle = interpn(givA,givM,givB,givH,CoeffsE.CLL,alpha,M,beta,-z);
+% Clpe = interpn(givA,givM,givB,givH,CoeffsE.CLLP,alpha,M,beta,-z);
+% Cmae = interpn(givA,givM,givB,givH,CoeffsE.CMA,alpha,M,beta,-z);
+% Cmade = interpn(givA,givM,givB,givH,CoeffsE.CMAD,alpha,M,beta,-z);
+% Cmqe = interpn(givA,givM,givB,givH,CoeffsE.CMQ,alpha,M,beta,-z);
+% Cnbe = interpn(givA,givM,givB,givH,CoeffsE.CLNB,alpha,M,beta,-z);
+% Cnre = interpn(givA,givM,givB,givH,CoeffsE.CLNR,alpha,M,beta,-z);
+% Cnpe = interpn(givA,givM,givB,givH,CoeffsE.CLNP,alpha,M,beta,-z);
+% %
+% Clre = interpn(givA,givM,givB,givH,CoeffsE.CLLR,alpha,M,beta,-z);
+% Clbe = interpn(givA,givM,givB,givH,CoeffsE.CLLB,alpha,M,beta,-z);
+% CDe =interpn(givA,givM,givB,givH,CoeffsE.CD,alpha,M,beta,-z);
+
+%Linear interpolation from empty and full configuration coefficients
+if t<tb
+ CA = t/tb*(CAe-CAf)+CAf;
+ CYB= t/tb*(CYBe-CYBf)+CYBf;
+ CNA= t/tb*(CNAe-CNAf)+CNAf;
+ Cl = t/tb*(Cle-Clf)+Clf;
+ Clp= t/tb*(Clpe-Clpf)+Clpf;
+ Cma= t/tb*(Cmae-Cmaf)+Cmaf;
+ Cmad=t/tb*(Cmade-Cmadf)+Cmadf;
+ Cmq= t/tb*(Cmqe-Cmqf)+Cmqf;
+ Cnb= t/tb*(Cnbe-Cnbf)+Cnbf;
+ Cnr= t/tb*(Cnre-Cnrf)+Cnrf;
+ Cnp= t/tb*(Cnpe-Cnpf)+Cnpf;
+ %
+ Clr= t/tb*(Clre-Clrf)+Clrf;
+ Clb= t/tb*(Clbe-Clbf)+Clbf;
+ CD= t/tb*(CDe-CDf)+CDf;
+else
+ CA = CAe;
+ CYB= CYBe;
+ CNA= CNAe;
+ Cl = Cle;
+ Clp= Clpe;
+ Cma= Cmae;
+ Cmad=Cmade;
+ Cmq= Cmqe;
+ Cnb= Cnbe;
+ Cnr= Cnre;
+ Cnp= Cnpe;
+ %
+ Clr= Clre;
+ Clb= Clbe;
+ CD = CDe;
+end
+
+
+global CA_plot % coefficiente aerodinamico
+if bool == 1
+ CA_plot(contatore) = CD; % Axial force coeff
+end
+
+
+%Center of Mass
+
+
+
+%Body Frame
+
+qdynS=0.5*rho*V_norm^2*S;
+qdynL_V = 0.5*rho*V_norm*S*C;
+
+X=qdynS*CA;
+Y=qdynS*CYB*beta;
+Z=qdynS*CNA*alpha;
+
+global Drag_plot
+if bool == 1
+ Drag_plot(contatore) = qdynS*CD;
+end
+
+
+%Forces
+F = quatrotate(Q,[0 0 m*g])';
+
+F = F +[-X+T,+Y,-Z]';
+
+% global Forces_plot
+% if bool == 1
+% Forces_plot(contatore) = F(1);
+% end
+
+
+du=F(1)/m -q*w + r*v;% - mdot/m*u;
+dv=F(2)/m -r*u + p*w;% - mdot/m*v;
+dw=F(3)/m -p*v + q*u;% - mdot/m*w;
+
+
+
+
+% Rotation
+dp=(Iyy-Izz)/Ixx*q*r + qdynL_V/Ixx*(V_norm*Cl+Clp*p*C/2)-Ixxdot*p/Ixx;
+dq=(Izz-Ixx)/Iyy*p*r + qdynL_V/Iyy*(V_norm*Cma*alpha + (Cmad+Cmq)*q*C/2)...
+ -Iyydot*q/Iyy;
+dr=(Ixx-Iyy)/Izz*p*q + qdynL_V/Izz*(V_norm*Cnb*beta + (Cnr*r+Cnp*p)*C/2)...
+ -Izzdot*r/Izz;
+
+
+% dp=(Iyy-Izz)/Ixx*q*r + qdynL_V/Ixx*(V_norm*(Cl+Clb*beta)+(Clp*p+Clr*r)*C/2)...
+% -Ixxdot*p/Ixx;
+% dq=(Izz-Ixx)/Iyy*p*r + qdynL_V/Iyy*(V_norm*Cma*alpha + (Cmad+Cmq)*q*C/2)...
+% -Iyydot*q/Iyy;
+% dr=(Ixx-Iyy)/Izz*p*q + qdynL_V/Izz*(V_norm*Cnb*beta + (Cnr*r+Cnp*p)*C/2)...
+% -Izzdot*r/Izz;
+
+
+
+% Launch Pad-relative coordinates
+Xb = quatrotate(Q,[x y z]);
+
+%
+
+%On Launch Pad No Torque
+ if Xb(1) < settings.lrampa
+ dp = 0;
+ dq = 0;
+ dr = 0;
+
+ %Ground Reaction.
+ if T < m*g
+ du = 0;
+ dv = 0;
+ dw = 0;
+
+ end
+end
+
+
+
+
+
+
+
+OM = 1/2*[
+ 0 -p -q -r
+ p 0 r -q
+ q -r 0 p
+ r q -p 0];
+
+
+dQQ = OM*Q';
+
+dY(1:3) = Vels;
+dY(4) = du;
+dY(5) = dv;
+dY(6) = dw;
+dY(7) = dp;
+dY(8) = dq;
+dY(9) = dr;
+dY(10:13) = dQQ;
+dY(14) = mdot;
+dY(15) = Ixxdot;
+dY(16) = Iyydot;
+dY(17) = Izzdot;
+dY=dY';
+
+if bool == 1
+ contatore = contatore + 1;
+end
+
+end
\ No newline at end of file
diff --git a/scratch/coeffs_data.m b/scratch/coeffs_data.m
new file mode 100644
index 0000000000000000000000000000000000000000..889999251b2ffefa20ea7ee497e40a6557e8309e
--- /dev/null
+++ b/scratch/coeffs_data.m
@@ -0,0 +1,23 @@
+% Author: Francesco Colombi
+% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu
+% email: francesco.colombi@skywarder.eu
+% Release date: 16/04/2016
+
+
+load for006_empty.mat
+
+iAlpha0 = find(State.Alphas == 0);
+iBeta0 = find(State.Betas == 0);
+
+M = State.Machs';
+altitudes = State.Altitudes';
+
+CD = [];
+for i = 1:length(altitudes)
+ CD(:, i) = Coeffs.CD(iAlpha0, :, iBeta0, i)';
+end
+
+CA = [];
+for i = 1:length(altitudes)
+ CA(:, i) = Coeffs.CA(iAlpha0, :, iBeta0, i)';
+end
\ No newline at end of file
diff --git a/scratch/rotating_axis.m b/scratch/rotating_axis.m
new file mode 100644
index 0000000000000000000000000000000000000000..5cae9808d8c02b93bbdcb73ac450cee39c4eb119
--- /dev/null
+++ b/scratch/rotating_axis.m
@@ -0,0 +1,74 @@
+function rotating_axis(R,limits)
+
+
+
+% limits: vector containing 6 boundaries for axis and as last term the time
+% for representation
+
+if nargin==1
+ x1= -1;
+ x2= 1;
+ y1= -1;
+ y2= 1;
+ z1= -1;
+ z2= 1;
+ time=0.01;
+else
+ x1=limits(1);
+ x2=limits(2);
+ y1=limits(3);
+ y2=limits(4);
+ z1=limits(5);
+ z2=limits(6);
+ time=limits(7);
+end
+
+if size(R,3)>1
+ [R]=square_to_line(R);
+end
+
+u1=R(:,1);
+u2=R(:,2);
+u3=R(:,3);
+v1=R(:,4);
+v2=R(:,5);
+v3=R(:,6);
+w1=R(:,7);
+w2=R(:,8);
+w3=R(:,9);
+
+figure
+
+for i=1:length(u1)
+
+ clf('reset')
+ grid on
+ axis on
+ axis([x1 x2 y1 y2 z1 z2])
+ hold on
+ quiver3(0,0,0,u1(i),v1(i),w1(i))
+ quiver3(0,0,0,u2(i),v2(i),w2(i))
+ quiver3(0,0,0,u3(i),v3(i),w3(i))
+ hold off
+
+ legend('x','y','z')
+
+
+ pause(time)
+
+end
+
+end
+
+function[R_l]=square_to_line(R_s)
+
+n=size(R_s,3);
+R_l=zeros(n,9);
+for i=1:n
+
+ R_l(i,:)=[R_s(1,1,i),R_s(2,1,i),R_s(3,1,i),R_s(1,2,i),R_s(2,2,i),R_s(3,2,i),R_s(1,3,i),R_s(2,3,i),R_s(3,3,i)];
+
+end
+
+
+end
\ No newline at end of file
diff --git a/scratch/test_windgen.m b/scratch/test_windgen.m
new file mode 100644
index 0000000000000000000000000000000000000000..942c2b25b240b947e85fb82547862a488def5dd4
--- /dev/null
+++ b/scratch/test_windgen.m
@@ -0,0 +1,7 @@
+% Test windgen
+
+for i=1:50
+ [u,v,w]=windgen(pi/2,pi/2,0,0,0,20);
+ quiver(u,v);
+ hold on
+end
diff --git a/scratch/vis_coeffs.m b/scratch/vis_coeffs.m
new file mode 100644
index 0000000000000000000000000000000000000000..a0961ac89563d7a7220ab6e4b3865103b839b633
--- /dev/null
+++ b/scratch/vis_coeffs.m
@@ -0,0 +1,48 @@
+% Author: Francesco Colombi
+% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu
+% email: francesco.colombi@skywarder.eu
+% Release date: 16/04/2016
+
+close all
+clear, clc
+
+load for006_empty.mat
+
+iAlpha0 = find(State.Alphas == 0);
+iBeta0 = find(State.Betas == 0);
+
+% for iMach = 1:length(State.Machs)
+% figure();
+% for iAlt = 1:length(State.Altitudes)
+% subplot(3, 1, 1), hold on, grid on;
+% title(['CL vs Alpha for Beta = 0, Mach = ', num2str(State.Machs(iMach))]);
+% plot(State.Alphas, Coeffs.CL(:, iMach, iBeta0, iAlt));
+% end
+%
+% for iAlt = 1:length(State.Altitudes)
+% subplot(3, 1, 2), hold on, grid on;
+% title(['CD vs Alpha for Beta = 0, Mach = ', num2str(State.Machs(iMach))]);
+% plot(State.Alphas, Coeffs.CD(:, iMach, iBeta0, iAlt));
+% end
+%
+% for iAlt = 1:length(State.Altitudes)
+% subplot(3, 1, 3), hold on, grid on;
+% title(['Cm vs Alpha for Beta = 0, Mach = ', num2str(State.Machs(iMach))]);
+% plot(State.Alphas, Coeffs.CM(:, iMach, iBeta0, iAlt));
+% end
+% end
+
+figure(), grid on, hold on;
+for iMach = 1:length(State.Machs)
+ for iAlt = 1
+ title('X_{CP} vs Alpha for Beta = 0, different Mach');
+ plot(State.Alphas, Coeffs.X_C_P(:, iMach, iBeta0, iAlt));
+ end
+end
+legend(num2str(State.Machs(1:end)'))
+
+figure();
+plot(State.Machs, Coeffs.X_C_P(iAlpha0, :, iBeta0, 1)), grid on;
+title('X_{CP} vs Mach, with alpha and beta equal to 0');
+xlabel('Mach [-]')
+ylabel('X_{CP} from moment center in ref. lenght, neg. AFT of moment center')
\ No newline at end of file
diff --git a/src/MAIN.m b/src/MAIN.m
new file mode 100644
index 0000000000000000000000000000000000000000..55799930ed70faf3327271c3f3c4e075f0b80841
--- /dev/null
+++ b/src/MAIN.m
@@ -0,0 +1,46 @@
+function [T,Y, Ta,Ya] = MAIN(settings)
+ %MAIN SCRIPT - This function retrieves all the parameters and run the
+ %simulation
+ %
+ % All the parameters are configured in config.m
+
+ % Author: Ruben Di Battista
+ % Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu
+ % email: ruben.dibattista@skywarder.eu
+ % Website: http://www.skywarder.eu
+ % April 2014; Last revision: 31.XII.2014
+ % License: 2-clause BSD
+
+ % Author: Francesco Colombi
+ % Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu
+ % email: francesco.colombi@skywarder.eu
+ % Release date: 16/04/2016
+
+ %Printing License
+ print_license();
+ tic
+
+ %Retrieving Parameters
+% run('config.m');
+
+
+ %Checking if stochastic or standard simulation needed
+ if settings.ballistic
+ fprintf('Standard Ballistic Simulation Fired...\n\n');
+ [T,Y, Ta,Ya]=std_run_ballistic(settings);
+ else
+ if settings.stoch.N > 1
+ fprintf('Stochastic Simulation Fired...\n\n');
+ if settings.stoch.parallel
+ [LP,Z]=stoch_run_p(settings);
+ else
+ [LP,Z]=stoch_run(settings);
+ end
+ else
+ fprintf('Standard Simulation Fired...\n\n');
+ [T,Y, Ta,Ya]=std_run(settings);
+ end
+ end
+
+ toc
+end
diff --git a/src/apogee.m b/src/apogee.m
new file mode 100644
index 0000000000000000000000000000000000000000..9b7830561cc9861617d6344c11b1afcf5d027669
--- /dev/null
+++ b/src/apogee.m
@@ -0,0 +1,28 @@
+function [ value,isterminal,direction ] = apogee( t,Y,settings,varargin )
+% Event function to stop simulation at apogee
+
+% Author: Ruben Di Battista
+% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu
+% email: ruben.dibattista@skywarder.eu
+% Website: http://www.skywarder.eu
+% April 2014; Last revision: 25.IV.2014
+% License: 2-clause BSD
+
+Q = Y(10:13)';
+
+%Inertial Frame velocities
+vels = quatrotate(quatconj(Q),Y(4:6)');
+
+%Stop checking if I'm in Propulsion Phase
+if t>settings.tb
+ value = vels(3);
+else
+ value = 1;
+end
+
+isterminal = 1;
+direction = 1;
+
+
+end
+
diff --git a/src/ascend.m b/src/ascend.m
new file mode 100644
index 0000000000000000000000000000000000000000..7af3dfc507d99bc2417a5e20991e38837f3ac353
--- /dev/null
+++ b/src/ascend.m
@@ -0,0 +1,454 @@
+function [ dY ] = ascend( t,Y,settings,uw,vw,ww )
+% ODE-Function of the 6DOF Rigid Rocket Model
+% State = ( x y z | u v w | p q r | q0 q1 q2 q3 )
+%
+% (x y z): NED Earth's Surface Centered Frame ("Inertial") coordinates
+% (u v w): body frame velocities
+% (p q r): body frame angular rates
+% (q0 q1 q2 q3): attitude unit quaternion
+%
+%
+% NOTE: To get the NED velocities the body-frame must be multiplied for the
+% conjugated of the current attitude quaternion
+% E.G.
+%
+%
+% quatrotate(quatconj(Y(:,10:13)),Y(:,4:6))
+
+% Author: Ruben Di Battista
+% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu
+% email: ruben.dibattista@skywarder.eu
+% Website: http://www.skywarder.eu
+% April 2014; Last revision: 31.XII.2014
+% License: 2-clause BSD
+
+% Author: Francesco Colombi
+% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu
+% email: francesco.colombi@skywarder.eu
+% Release date: 16/04/2016
+
+x = Y(1);
+y = Y(2);
+z = Y(3);
+u = Y(4);
+v = Y(5);
+w = Y(6);
+p = Y(7);
+q = Y(8);
+r = Y(9);
+q0 = Y(10);
+q1 = Y(11);
+q2 = Y(12);
+q3 = Y(13);
+m = Y(14);
+Ixx = Y(15);
+Iyy = Y(16);
+Izz = Y(17);
+
+global bool
+% if bool = 0 the trends during the integration are NOT requested
+% if bool = 1 the trends during the integration are saved and plotted
+global t_plot contatore
+if bool == 1
+ t_plot(contatore) = t;
+end
+
+
+Q = [ q0 q1 q2 q3];
+Q_conj= [ q0 -q1 -q2 -q3];
+normQ=norm(Q);
+
+if abs(normQ-1)>0.1
+ Q = Q/normQ;
+ %frprintf('Norm of the quaternion is (too much!) less than one (%3.2f)\n', ...
+ %normQ);
+end
+
+
+% Adding Wind (supposed to be added in NED axes);
+
+% constant wind
+wind = quatrotate(Q, [uw vw ww]);
+
+% % Wind model
+% h = 0;
+% if -z > 0
+% h = -z+settings.z0;
+% end
+% % global alt
+% % if bool == 1
+% % alt = [alt; h];
+% % end
+%
+% % wind in NED axis
+% [uw,vw] = atmoshwm07(settings.wind.Lat, settings.wind.Long, h, ...
+% 'day',settings.wind.Day,'model','quiet');
+%
+% % global WIND
+% % if bool == 1
+% % WIND = [WIND; uw, vw, ww];
+% % end
+%
+% % wind in body frame
+% wind = quatrotate(Q, [uw, vw, ww]);
+
+
+% % Wind Model
+% [uw,vw,ww] = wind_generator(settings,z,t,Q);
+% wind = [ uw,vw,ww ];
+
+% Relative velocities (plus wind);
+ur = u - wind(1);
+vr = v - wind(2);
+wr = w - wind(3);
+
+% Body to Inertial velocities
+Vels = quatrotate(Q_conj,[u v w]);
+
+V_norm = norm([ur vr wr]);
+
+
+
+% Constants
+S = settings.S;
+C = settings.C;
+CoeffsE = settings.CoeffsE; %Empty Rocket Coefficients
+CoeffsF = settings.CoeffsF; % Full Rocket Coefficients
+
+g = 9.81;
+%Time-Dependant Properties
+m0 = settings.m0; %kg mass
+mfr = settings.mfr;
+tb = settings.tb;
+
+Ixxf= settings.Ixxf;
+Ixxe= settings.Ixxe;
+Iyyf= settings.Iyyf;
+Iyye= settings.Iyye;
+Izzf= settings.Izzf;
+Izze= settings.Izze;
+
+dI = 1/tb*([Ixxf Iyyf Izzf]'-[Ixxe Iyye Izze]');
+
+if t<tb
+ mdot = -mfr;
+ Ixxdot = -dI(1);
+ Iyydot = -dI(2);
+ Izzdot = -dI(3);
+% T = settings.T_coeffs*[t^0 t^1 t^2 t^3 t^4]';
+ T = interp1(settings.motor.exp_time, settings.motor.exp_thrust, t);
+
+else
+ mdot = 0;
+ Ixxdot = 0;
+ Iyydot = 0;
+ Izzdot = 0;
+ T=0;
+end
+
+global T_plot
+if bool == 1
+ T_plot(contatore) = T; % thrust
+end
+
+%Atmosphere
+if -z<0
+ z = 0;
+end
+[~, a, ~, rho] = atmoscoesa(-z+settings.z0);
+M = V_norm/a;
+
+global M_plot
+if bool == 1
+ M_plot(contatore) = M; % mach
+end
+
+%Aero Angles
+if not(u<1e-1 || V_norm<1e-3);
+ alpha = atan(wr/ur);
+ beta = asin(vr/V_norm);
+else
+ alpha =0;
+ beta = 0;
+end
+
+global alpha_plot beta_plot
+if bool == 1
+ alpha_plot(contatore) = alpha;
+ beta_plot(contatore) = beta;
+end
+
+% beta = 0; % prova simulazione con beta imposto uguale a 0
+
+
+
+% Coeffs. Interpolation
+givA = settings.Alphas*pi/180;
+givB = settings.Betas*pi/180;
+givH = settings.Altitudes;
+givM = settings.Machs;
+
+% Interpolation error check
+
+if M > givM(end)
+ %frprintf('Mach No. is more than the max provided (M = %3.2f @ t=%2.2f)\n\n', ...
+ %M,t);
+ M = givM(end);
+
+end
+
+if M<givM(1)
+ %frprintf('Mach No. is less than the min provided (M = %3.2f @ t=%2.2f)\n\n', ...
+ % M,t);
+ M = givM(1);
+
+end
+
+if alpha > givA(end)
+ %frprintf('AoA is more than the max provided (alpha = %3.2f @ t=%2.2f)\n\n', ...
+ %alpha*180/pi,t)
+ alpha= givA(end);
+
+else if alpha < givA(1)
+ %frprintf('AoA is less than the min provided (alpha = %3.2f @ t=%2.2f)\n\n',...
+ %alpha*180/pi,t);
+ alpha = givA(1);
+
+
+ end
+end
+
+if beta > givB(end)
+ beta= givB(end);
+else if beta < givB(1)
+ beta = givB(1);
+ end
+end
+
+if -z >givH(end)
+ z = -givH(end);
+else if -z < givH(1)
+ z = -givH(1);
+ end
+end
+
+% interpolation with the value in the nearest condition: interp_easy
+% full
+CAf=interp4_easy(givA,givM,givB,givH,CoeffsF.CA,alpha,M,beta,-z);%,'nearest');
+CYBf=interp4_easy(givA,givM,givB,givH,CoeffsF.CYB,alpha,M,beta,-z);%,'nearest');
+CNAf=interp4_easy(givA,givM,givB,givH,CoeffsF.CNA,alpha,M,beta,-z);%,'nearest');
+Clf=interp4_easy(givA,givM,givB,givH,CoeffsF.CLL,alpha,M,beta,-z);%,'nearest');
+Clpf=interp4_easy(givA,givM,givB,givH,CoeffsF.CLLP,alpha,M,beta,-z);%,'nearest');
+Cmaf=interp4_easy(givA,givM,givB,givH,CoeffsF.CMA,alpha,M,beta,-z);%,'nearest');
+Cmadf=interp4_easy(givA,givM,givB,givH,CoeffsF.CMAD,alpha,M,beta,-z);%,'nearest');
+Cmqf=interp4_easy(givA,givM,givB,givH,CoeffsF.CMQ,alpha,M,beta,-z);%,'nearest');
+Cnbf=interp4_easy(givA,givM,givB,givH,CoeffsF.CLNB,alpha,M,beta,-z);%,'nearest');
+Cnrf=interp4_easy(givA,givM,givB,givH,CoeffsF.CLNR,alpha,M,beta,-z);%,'nearest');
+Cnpf=interp4_easy(givA,givM,givB,givH,CoeffsF.CLNP,alpha,M,beta,-z);%,'nearest');
+%
+Clrf=interp4_easy(givA,givM,givB,givH,CoeffsF.CLLR,alpha,M,beta,-z);%,'nearest');
+Clbf=interp4_easy(givA,givM,givB,givH,CoeffsF.CLLB,alpha,M,beta,-z);%,'nearest');
+CDf =interp4_easy(givA,givM,givB,givH,CoeffsF.CD,alpha,M,beta,-z);
+
+% empty
+CAe=interp4_easy(givA,givM,givB,givH,CoeffsE.CA,alpha,M,beta,-z);%,'nearest');
+CYBe=interp4_easy(givA,givM,givB,givH,CoeffsE.CYB,alpha,M,beta,-z);%,'nearest');
+CNAe=interp4_easy(givA,givM,givB,givH,CoeffsE.CNA,alpha,M,beta,-z);%,'nearest');
+Cle=interp4_easy(givA,givM,givB,givH,CoeffsE.CLL,alpha,M,beta,-z);%,'nearest');
+Clpe=interp4_easy(givA,givM,givB,givH,CoeffsE.CLLP,alpha,M,beta,-z);%,'nearest');
+Cmae=interp4_easy(givA,givM,givB,givH,CoeffsE.CMA,alpha,M,beta,-z);%,'nearest');
+Cmade=interp4_easy(givA,givM,givB,givH,CoeffsE.CMAD,alpha,M,beta,-z);%,'nearest');
+Cmqe=interp4_easy(givA,givM,givB,givH,CoeffsE.CMQ,alpha,M,beta,-z);%,'nearest');
+Cnbe=interp4_easy(givA,givM,givB,givH,CoeffsE.CLNB,alpha,M,beta,-z);%,'nearest');
+Cnre=interp4_easy(givA,givM,givB,givH,CoeffsE.CLNR,alpha,M,beta,-z);%,'nearest');
+Cnpe=interp4_easy(givA,givM,givB,givH,CoeffsE.CLNP,alpha,M,beta,-z);%,'nearest');
+%
+Clre=interp4_easy(givA,givM,givB,givH,CoeffsE.CLLR,alpha,M,beta,-z);%,'nearest');
+Clbe=interp4_easy(givA,givM,givB,givH,CoeffsE.CLLB,alpha,M,beta,-z);%,'nearest');
+CDe =interp4_easy(givA,givM,givB,givH,CoeffsE.CD,alpha,M,beta,-z);
+
+% % linear interpolation of the coeff
+% %full
+% CAf = interpn(givA,givM,givB,givH,CoeffsF.CA,alpha,M,beta,-z);
+% CYBf = interpn(givA,givM,givB,givH,CoeffsF.CYB,alpha,M,beta,-z);
+% CNAf = interpn(givA,givM,givB,givH,CoeffsF.CNA,alpha,M,beta,-z);
+% Clf = interpn(givA,givM,givB,givH,CoeffsF.CLL,alpha,M,beta,-z);
+% Clpf = interpn(givA,givM,givB,givH,CoeffsF.CLLP,alpha,M,beta,-z);
+% Cmaf = interpn(givA,givM,givB,givH,CoeffsF.CMA,alpha,M,beta,-z);
+% Cmadf = interpn(givA,givM,givB,givH,CoeffsF.CMAD,alpha,M,beta,-z);
+% Cmqf = interpn(givA,givM,givB,givH,CoeffsF.CMQ,alpha,M,beta,-z);
+% Cnbf = interpn(givA,givM,givB,givH,CoeffsF.CLNB,alpha,M,beta,-z);
+% Cnrf = interpn(givA,givM,givB,givH,CoeffsF.CLNR,alpha,M,beta,-z);
+% Cnpf = interpn(givA,givM,givB,givH,CoeffsF.CLNP,alpha,M,beta,-z);
+% %
+% Clrf = interpn(givA,givM,givB,givH,CoeffsF.CLLR,alpha,M,beta,-z);
+% Clbf = interpn(givA,givM,givB,givH,CoeffsF.CLLB,alpha,M,beta,-z);
+% CDf = interpn(givA,givM,givB,givH,CoeffsF.CD,alpha,M,beta,-z);
+%
+% % empty
+% CAe = interpn(givA,givM,givB,givH,CoeffsE.CA,alpha,M,beta,-z);
+% CYBe = interpn(givA,givM,givB,givH,CoeffsE.CYB,alpha,M,beta,-z);
+% CNAe = interpn(givA,givM,givB,givH,CoeffsE.CNA,alpha,M,beta,-z);
+% Cle = interpn(givA,givM,givB,givH,CoeffsE.CLL,alpha,M,beta,-z);
+% Clpe = interpn(givA,givM,givB,givH,CoeffsE.CLLP,alpha,M,beta,-z);
+% Cmae = interpn(givA,givM,givB,givH,CoeffsE.CMA,alpha,M,beta,-z);
+% Cmade = interpn(givA,givM,givB,givH,CoeffsE.CMAD,alpha,M,beta,-z);
+% Cmqe = interpn(givA,givM,givB,givH,CoeffsE.CMQ,alpha,M,beta,-z);
+% Cnbe = interpn(givA,givM,givB,givH,CoeffsE.CLNB,alpha,M,beta,-z);
+% Cnre = interpn(givA,givM,givB,givH,CoeffsE.CLNR,alpha,M,beta,-z);
+% Cnpe = interpn(givA,givM,givB,givH,CoeffsE.CLNP,alpha,M,beta,-z);
+% %
+% Clre = interpn(givA,givM,givB,givH,CoeffsE.CLLR,alpha,M,beta,-z);
+% Clbe = interpn(givA,givM,givB,givH,CoeffsE.CLLB,alpha,M,beta,-z);
+% CDe =interpn(givA,givM,givB,givH,CoeffsE.CD,alpha,M,beta,-z);
+
+%Linear interpolation from empty and full configuration coefficients
+if t<tb
+ CA = t/tb*(CAe-CAf)+CAf;
+ CYB= t/tb*(CYBe-CYBf)+CYBf;
+ CNA= t/tb*(CNAe-CNAf)+CNAf;
+ Cl = t/tb*(Cle-Clf)+Clf;
+ Clp= t/tb*(Clpe-Clpf)+Clpf;
+ Cma= t/tb*(Cmae-Cmaf)+Cmaf;
+ Cmad=t/tb*(Cmade-Cmadf)+Cmadf;
+ Cmq= t/tb*(Cmqe-Cmqf)+Cmqf;
+ Cnb= t/tb*(Cnbe-Cnbf)+Cnbf;
+ Cnr= t/tb*(Cnre-Cnrf)+Cnrf;
+ Cnp= t/tb*(Cnpe-Cnpf)+Cnpf;
+ %
+ Clr= t/tb*(Clre-Clrf)+Clrf;
+ Clb= t/tb*(Clbe-Clbf)+Clbf;
+ CD= t/tb*(CDe-CDf)+CDf;
+else
+ CA = CAe;
+ CYB= CYBe;
+ CNA= CNAe;
+ Cl = Cle;
+ Clp= Clpe;
+ Cma= Cmae;
+ Cmad=Cmade;
+ Cmq= Cmqe;
+ Cnb= Cnbe;
+ Cnr= Cnre;
+ Cnp= Cnpe;
+ %
+ Clr= Clre;
+ Clb= Clbe;
+ CD = CDe;
+end
+
+
+global CA_plot % coefficiente aerodinamico
+if bool == 1
+ CA_plot(contatore) = CD; % Axial force coeff
+end
+
+
+%Center of Mass
+
+
+
+%Body Frame
+
+qdynS=0.5*rho*V_norm^2*S;
+qdynL_V = 0.5*rho*V_norm*S*C;
+
+X=qdynS*CA;
+Y=qdynS*CYB*beta;
+Z=qdynS*CNA*alpha;
+
+global Drag_plot
+if bool == 1
+ Drag_plot(contatore) = qdynS*CD;
+end
+
+
+%Forces
+F = quatrotate(Q,[0 0 m*g])';
+
+F = F +[-X+T,+Y,-Z]';
+
+global Forces_plot
+if bool == 1
+ Forces_plot(contatore) = F(1);
+end
+
+
+du=F(1)/m -q*w + r*v;% - mdot/m*u;
+dv=F(2)/m -r*u + p*w;% - mdot/m*v;
+dw=F(3)/m -p*v + q*u;% - mdot/m*w;
+
+
+
+
+% Rotation
+dp=(Iyy-Izz)/Ixx*q*r + qdynL_V/Ixx*(V_norm*Cl+Clp*p*C/2)-Ixxdot*p/Ixx;
+dq=(Izz-Ixx)/Iyy*p*r + qdynL_V/Iyy*(V_norm*Cma*alpha + (Cmad+Cmq)*q*C/2)...
+ -Iyydot*q/Iyy;
+dr=(Ixx-Iyy)/Izz*p*q + qdynL_V/Izz*(V_norm*Cnb*beta + (Cnr*r+Cnp*p)*C/2)...
+ -Izzdot*r/Izz;
+
+
+% dp=(Iyy-Izz)/Ixx*q*r + qdynL_V/Ixx*(V_norm*(Cl+Clb*beta)+(Clp*p+Clr*r)*C/2)...
+% -Ixxdot*p/Ixx;
+% dq=(Izz-Ixx)/Iyy*p*r + qdynL_V/Iyy*(V_norm*Cma*alpha + (Cmad+Cmq)*q*C/2)...
+% -Iyydot*q/Iyy;
+% dr=(Ixx-Iyy)/Izz*p*q + qdynL_V/Izz*(V_norm*Cnb*beta + (Cnr*r+Cnp*p)*C/2)...
+% -Izzdot*r/Izz;
+
+
+
+% Launch Pad-relative coordinates
+Xb = quatrotate(Q,[x y z]);
+
+%
+
+%On Launch Pad No Torque
+ if Xb(1) < settings.lrampa
+ dp = 0;
+ dq = 0;
+ dr = 0;
+
+ %Ground Reaction.
+ if T < m*g
+ du = 0;
+ dv = 0;
+ dw = 0;
+
+ end
+end
+
+
+
+
+
+
+
+OM = 1/2*[
+ 0 -p -q -r
+ p 0 r -q
+ q -r 0 p
+ r q -p 0];
+
+
+dQQ = OM*Q';
+
+dY(1:3) = Vels;
+dY(4) = du;
+dY(5) = dv;
+dY(6) = dw;
+dY(7) = dp;
+dY(8) = dq;
+dY(9) = dr;
+dY(10:13) = dQQ;
+dY(14) = mdot;
+dY(15) = Ixxdot;
+dY(16) = Iyydot;
+dY(17) = Izzdot;
+dY=dY';
+
+if bool == 1
+ contatore = contatore + 1;
+end
+
+end
\ No newline at end of file
diff --git a/src/ballistic_descent.m b/src/ballistic_descent.m
new file mode 100644
index 0000000000000000000000000000000000000000..a1fbb65748ebb1e97a59020bf48addf3ed33e649
--- /dev/null
+++ b/src/ballistic_descent.m
@@ -0,0 +1,339 @@
+function [ dY ] = ballistic_descent( t,Y,settings,uw,vw,ww )
+% ODE-Function of the 6DOF Rigid Rocket Model
+% State = ( x y z | u v w | p q r | q0 q1 q2 q3 )
+%
+% (x y z): NED Earth's Surface Centered Frame ("Inertial") coordinates
+% (u v w): body frame velocities
+% (p q r): body frame angular rates
+% (q0 q1 q2 q3): attitude unit quaternion
+%
+%
+% NOTE: To get the NED velocities the body-frame must be multiplied for the
+% conjugated of the current attitude quaternion
+% E.G.
+%
+%
+% quatrotate(quatconj(Y(:,10:13)),Y(:,4:6))
+
+% Author: Ruben Di Battista
+% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu
+% email: ruben.dibattista@skywarder.eu
+% Website: http://www.skywarder.eu
+% April 2014; Last revision: 31.XII.2014
+% License: 2-clause BSD
+
+% Author: Francesco Colombi
+% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu
+% email: francesco.colombi@skywarder.eu
+% Release date: 16/04/2016
+
+x = Y(1);
+y = Y(2);
+z = Y(3);
+u = Y(4);
+v = Y(5);
+w = Y(6);
+p = Y(7);
+q = Y(8);
+r = Y(9);
+q0 = Y(10);
+q1 = Y(11);
+q2 = Y(12);
+q3 = Y(13);
+m = Y(14);
+Ixx = Y(15);
+Iyy = Y(16);
+Izz = Y(17);
+
+
+Q = [ q0 q1 q2 q3];
+Q_conj= [ q0 -q1 -q2 -q3];
+normQ=norm(Q);
+
+if abs(normQ-1)>0.1
+ Q = Q/normQ;
+ %frprintf('Norm of the quaternion is (too much!) less than one (%3.2f)\n', ...
+ %normQ);
+end
+
+
+% Adding Wind (supposed to be added in NED axes);
+
+% constant wind
+wind = quatrotate(Q, [uw vw ww]);
+
+% % Wind model
+% h = 0;
+% if -z > 0
+% h = -z+settings.z0;
+% end
+% % global alt
+% % if bool == 1
+% % alt = [alt; h];
+% % end
+%
+% % wind in NED axis
+% [uw,vw] = atmoshwm07(settings.wind.Lat, settings.wind.Long, h, ...
+% 'day',settings.wind.Day,'model','quiet');
+%
+% % global WIND
+% % if bool == 1
+% % WIND = [WIND; uw, vw, ww];
+% % end
+%
+% % wind in body frame
+% wind = quatrotate(Q, [uw, vw, ww]);
+
+
+% % Wind Model
+% [uw,vw,ww] = wind_generator(settings,z,t,Q);
+% wind = [ uw,vw,ww ];
+
+% Relative velocities (plus wind);
+ur = u - wind(1);
+vr = v - wind(2);
+wr = w - wind(3);
+
+% Body to Inertial velocities
+Vels = quatrotate(Q_conj,[u v w]);
+
+V_norm = norm([ur vr wr]);
+
+
+
+% Constants
+S = settings.S;
+C = settings.C;
+CoeffsE = settings.CoeffsE; %Empty Rocket Coefficients
+
+g = 9.81;
+%Time-Dependant Properties
+m0 = settings.m0; %kg mass
+mfr = settings.mfr;
+tb = settings.tb;
+
+Ixxf= settings.Ixxf;
+Ixxe= settings.Ixxe;
+Iyyf= settings.Iyyf;
+Iyye= settings.Iyye;
+Izzf= settings.Izzf;
+Izze= settings.Izze;
+
+dI = 1/tb*([Ixxf Iyyf Izzf]'-[Ixxe Iyye Izze]');
+
+ mdot = 0;
+ Ixxdot = 0;
+ Iyydot = 0;
+ Izzdot = 0;
+ T=0;
+
+%Atmosphere
+if -z<0
+ z = 0;
+end
+[~, a, ~, rho] = atmoscoesa(-z+settings.z0);
+M = V_norm/a;
+
+
+%Aero Angles
+if not(u<1e-1 || V_norm<1e-3);
+ alpha = atan(wr/ur);
+ beta = asin(vr/V_norm);
+else
+ alpha =0;
+ beta = 0;
+end
+
+% beta = 0; % prova simulazione con beta imposto uguale a 0
+
+
+
+% Coeffs. Interpolation
+givA = settings.Alphas*pi/180;
+givB = settings.Betas*pi/180;
+givH = settings.Altitudes;
+givM = settings.Machs;
+
+% Interpolation error check
+
+if M > givM(end)
+ %frprintf('Mach No. is more than the max provided (M = %3.2f @ t=%2.2f)\n\n', ...
+ %M,t);
+ M = givM(end);
+
+end
+
+if M<givM(1)
+ %frprintf('Mach No. is less than the min provided (M = %3.2f @ t=%2.2f)\n\n', ...
+ % M,t);
+ M = givM(1);
+
+end
+
+if alpha > givA(end)
+ %frprintf('AoA is more than the max provided (alpha = %3.2f @ t=%2.2f)\n\n', ...
+ %alpha*180/pi,t)
+ alpha= givA(end);
+
+else if alpha < givA(1)
+ %frprintf('AoA is less than the min provided (alpha = %3.2f @ t=%2.2f)\n\n',...
+ %alpha*180/pi,t);
+ alpha = givA(1);
+
+
+ end
+end
+
+if beta > givB(end)
+ beta= givB(end);
+else if beta < givB(1)
+ beta = givB(1);
+ end
+end
+
+if -z >givH(end)
+ z = -givH(end);
+else if -z < givH(1)
+ z = -givH(1);
+ end
+end
+
+% interpolation with the value in the nearest condition: interp_easy
+% empty
+CAe=interp4_easy(givA,givM,givB,givH,CoeffsE.CA,alpha,M,beta,-z);%,'nearest');
+CYBe=interp4_easy(givA,givM,givB,givH,CoeffsE.CYB,alpha,M,beta,-z);%,'nearest');
+CNAe=interp4_easy(givA,givM,givB,givH,CoeffsE.CNA,alpha,M,beta,-z);%,'nearest');
+Cle=interp4_easy(givA,givM,givB,givH,CoeffsE.CLL,alpha,M,beta,-z);%,'nearest');
+Clpe=interp4_easy(givA,givM,givB,givH,CoeffsE.CLLP,alpha,M,beta,-z);%,'nearest');
+Cmae=interp4_easy(givA,givM,givB,givH,CoeffsE.CMA,alpha,M,beta,-z);%,'nearest');
+Cmade=interp4_easy(givA,givM,givB,givH,CoeffsE.CMAD,alpha,M,beta,-z);%,'nearest');
+Cmqe=interp4_easy(givA,givM,givB,givH,CoeffsE.CMQ,alpha,M,beta,-z);%,'nearest');
+Cnbe=interp4_easy(givA,givM,givB,givH,CoeffsE.CLNB,alpha,M,beta,-z);%,'nearest');
+Cnre=interp4_easy(givA,givM,givB,givH,CoeffsE.CLNR,alpha,M,beta,-z);%,'nearest');
+Cnpe=interp4_easy(givA,givM,givB,givH,CoeffsE.CLNP,alpha,M,beta,-z);%,'nearest');
+%
+Clre=interp4_easy(givA,givM,givB,givH,CoeffsE.CLLR,alpha,M,beta,-z);%,'nearest');
+Clbe=interp4_easy(givA,givM,givB,givH,CoeffsE.CLLB,alpha,M,beta,-z);%,'nearest');
+CDe =interp4_easy(givA,givM,givB,givH,CoeffsE.CD,alpha,M,beta,-z);
+
+% % linear interpolation of the coeff
+% %full
+% CAf = interpn(givA,givM,givB,givH,CoeffsF.CA,alpha,M,beta,-z);
+% CYBf = interpn(givA,givM,givB,givH,CoeffsF.CYB,alpha,M,beta,-z);
+% CNAf = interpn(givA,givM,givB,givH,CoeffsF.CNA,alpha,M,beta,-z);
+% Clf = interpn(givA,givM,givB,givH,CoeffsF.CLL,alpha,M,beta,-z);
+% Clpf = interpn(givA,givM,givB,givH,CoeffsF.CLLP,alpha,M,beta,-z);
+% Cmaf = interpn(givA,givM,givB,givH,CoeffsF.CMA,alpha,M,beta,-z);
+% Cmadf = interpn(givA,givM,givB,givH,CoeffsF.CMAD,alpha,M,beta,-z);
+% Cmqf = interpn(givA,givM,givB,givH,CoeffsF.CMQ,alpha,M,beta,-z);
+% Cnbf = interpn(givA,givM,givB,givH,CoeffsF.CLNB,alpha,M,beta,-z);
+% Cnrf = interpn(givA,givM,givB,givH,CoeffsF.CLNR,alpha,M,beta,-z);
+% Cnpf = interpn(givA,givM,givB,givH,CoeffsF.CLNP,alpha,M,beta,-z);
+% %
+% Clrf = interpn(givA,givM,givB,givH,CoeffsF.CLLR,alpha,M,beta,-z);
+% Clbf = interpn(givA,givM,givB,givH,CoeffsF.CLLB,alpha,M,beta,-z);
+% CDf = interpn(givA,givM,givB,givH,CoeffsF.CD,alpha,M,beta,-z);
+%
+% % empty
+% CAe = interpn(givA,givM,givB,givH,CoeffsE.CA,alpha,M,beta,-z);
+% CYBe = interpn(givA,givM,givB,givH,CoeffsE.CYB,alpha,M,beta,-z);
+% CNAe = interpn(givA,givM,givB,givH,CoeffsE.CNA,alpha,M,beta,-z);
+% Cle = interpn(givA,givM,givB,givH,CoeffsE.CLL,alpha,M,beta,-z);
+% Clpe = interpn(givA,givM,givB,givH,CoeffsE.CLLP,alpha,M,beta,-z);
+% Cmae = interpn(givA,givM,givB,givH,CoeffsE.CMA,alpha,M,beta,-z);
+% Cmade = interpn(givA,givM,givB,givH,CoeffsE.CMAD,alpha,M,beta,-z);
+% Cmqe = interpn(givA,givM,givB,givH,CoeffsE.CMQ,alpha,M,beta,-z);
+% Cnbe = interpn(givA,givM,givB,givH,CoeffsE.CLNB,alpha,M,beta,-z);
+% Cnre = interpn(givA,givM,givB,givH,CoeffsE.CLNR,alpha,M,beta,-z);
+% Cnpe = interpn(givA,givM,givB,givH,CoeffsE.CLNP,alpha,M,beta,-z);
+% %
+% Clre = interpn(givA,givM,givB,givH,CoeffsE.CLLR,alpha,M,beta,-z);
+% Clbe = interpn(givA,givM,givB,givH,CoeffsE.CLLB,alpha,M,beta,-z);
+% CDe =interpn(givA,givM,givB,givH,CoeffsE.CD,alpha,M,beta,-z);
+
+%Linear interpolation from empty and full configuration coefficients
+ CA = CAe;
+ CYB= CYBe;
+ CNA= CNAe;
+ Cl = Cle;
+ Clp= Clpe;
+ Cma= Cmae;
+ Cmad=Cmade;
+ Cmq= Cmqe;
+ Cnb= Cnbe;
+ Cnr= Cnre;
+ Cnp= Cnpe;
+ %
+ Clr= Clre;
+ Clb= Clbe;
+ CD = CDe;
+
+
+%Center of Mass
+
+
+
+%Body Frame
+
+qdynS=0.5*rho*V_norm^2*S;
+qdynL_V = 0.5*rho*V_norm*S*C;
+
+X=qdynS*CA;
+Y=qdynS*CYB*beta;
+Z=qdynS*CNA*alpha;
+
+
+
+%Forces
+F = quatrotate(Q,[0 0 m*g])';
+
+F = F +[-X+T,+Y,-Z]';
+
+
+du=F(1)/m -q*w + r*v;% - mdot/m*u;
+dv=F(2)/m -r*u + p*w;% - mdot/m*v;
+dw=F(3)/m -p*v + q*u;% - mdot/m*w;
+
+
+
+
+% Rotation
+dp=(Iyy-Izz)/Ixx*q*r + qdynL_V/Ixx*(V_norm*Cl+Clp*p*C/2)-Ixxdot*p/Ixx;
+dq=(Izz-Ixx)/Iyy*p*r + qdynL_V/Iyy*(V_norm*Cma*alpha + (Cmad+Cmq)*q*C/2)...
+ -Iyydot*q/Iyy;
+dr=(Ixx-Iyy)/Izz*p*q + qdynL_V/Izz*(V_norm*Cnb*beta + (Cnr*r+Cnp*p)*C/2)...
+ -Izzdot*r/Izz;
+
+
+% dp=(Iyy-Izz)/Ixx*q*r + qdynL_V/Ixx*(V_norm*(Cl+Clb*beta)+(Clp*p+Clr*r)*C/2)...
+% -Ixxdot*p/Ixx;
+% dq=(Izz-Ixx)/Iyy*p*r + qdynL_V/Iyy*(V_norm*Cma*alpha + (Cmad+Cmq)*q*C/2)...
+% -Iyydot*q/Iyy;
+% dr=(Ixx-Iyy)/Izz*p*q + qdynL_V/Izz*(V_norm*Cnb*beta + (Cnr*r+Cnp*p)*C/2)...
+% -Izzdot*r/Izz;
+
+
+OM = 1/2*[
+ 0 -p -q -r
+ p 0 r -q
+ q -r 0 p
+ r q -p 0];
+
+
+dQQ = OM*Q';
+
+dY(1:3) = Vels;
+dY(4) = du;
+dY(5) = dv;
+dY(6) = dw;
+dY(7) = dp;
+dY(8) = dq;
+dY(9) = dr;
+dY(10:13) = dQQ;
+dY(14) = mdot;
+dY(15) = Ixxdot;
+dY(16) = Iyydot;
+dY(17) = Izzdot;
+dY=dY';
+
+end
\ No newline at end of file
diff --git a/src/config.m b/src/config.m
new file mode 100644
index 0000000000000000000000000000000000000000..6d0898ec4ed9c56ba5592c66e52318fd20e63bdd
--- /dev/null
+++ b/src/config.m
@@ -0,0 +1,225 @@
+%CONFIG SCRIPT - This script sets up all the parameters for the simulation (H1 line)
+% All the parameters are stored in the "settings" structure.
+%
+
+% TODO: GUI to fill automatically this configuration script
+
+% Author: Ruben Di Battista
+% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu
+% email: ruben.dibattista@skywarder.eu
+% Website: http://www.skywarder.eu
+% License: 2-clause BSD
+
+% Author: Francesco Colombi
+% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu
+% email: francesco.colombi@skywarder.eu
+% Release date: 16/04/2016
+
+% clear all
+% close all
+% clc
+
+% LAUNCHPAD %
+settings.z0 = 100; %Launchpad Altitude
+settings.lrampa = 4.443; %LaunchPad route (launchpad lenght-distance from ground of the first hook)
+
+
+%STARTING ATTITUDE SETUP %
+settings.OMEGA = 85*pi/180; %Elevation Angle
+settings.PHI = 270*pi/180; %Azimuth Angle from North Direction
+
+% ENGINE DETAILS %
+
+% % Thrust on time profile (IV Order Interpolation)
+% settings.T_coeffs = [8034.5, 0, 0, 0, 0]; % constant thrust
+% data from the motor details sheet (experimental plot of the Thrust)
+
+% SYNTAX:
+% engine = 1 -> Cesaroni PRO 150 White Thunder
+% engine = 2 -> Cesaroni PRO 150 SkidMark
+% engine = 3 -> Cesaroni PRO 150 BlueStreak
+engine = 3;
+switch engine
+ case 1
+ % Cesaroni PRO 150 White Thunder
+ % Sampling for thrust interpolation
+ settings.motor.exp_time = [0 0.05 0.15 0.5 0.6 0.74 0.85 1.15 1.7 2.4 3 ...
+ 4 4.5 4.8 4.9 5 5.05 5.1 5.15 5.2]; %s
+ settings.motor.exp_thrust = [8605.1 8900 7900 8400 8400 8250 8200 8300 ...
+ 8400 8400 8200 7800 7600 7450 7350 7300 4500 500 100 0]; %N
+
+ settings.m0 = 67.761; %kg Overall Mass (Burnout + Propellant)
+ settings.ms = 43.961; %kg Structural Mass (Burnout - Nosecone)
+ settings.mp = 18.6; %kg Propellant Mass
+ settings.tb = 5.12; %sec Burning Time
+ settings.mfr = settings.mp/settings.tb; %kg/s Mass Flow Rate
+ case 2
+ % Cesaroni PRO 150 SkidMark
+ % Sampling for thrust interpolation
+ settings.motor.exp_time = [0 0.1 0.2 0.3 0.4 0.5 0.6 0.8 1.2 1.8 3.2 ...
+ 3.6 4.8 6 7 7.2 7.6 7.8 7.9 8 8.1 8.19]; %s
+ settings.motor.exp_thrust = [0 3400 3100 3000 3300 3400 3500 3700 3700 ...
+ 3800 4000 4081.6 3900 3800 3700 3500 3350 3200 3000 2000 750 0]; %N
+
+ settings.m0 = 66.052; %kg Overall Mass (Burnout + Propellant)
+ settings.ms = 43.695; %kg Structural Mass (Burnout - Nosecone)
+ settings.mp = 17.157; %kg Propellant Mass
+ settings.tb = 8.19; %sec Burning Time
+ settings.mfr = settings.mp/settings.tb; %kg/s Mass Flow Rate
+ case 3
+ % Cesaroni PRO 150 BlueStreak
+ % Sampling for thrust interpolation
+ settings.motor.exp_time = [0 0.06 0.1 0.15 0.25 0.45 0.8 1 2 3 ...
+ 4 5 6 6.8 7.05 7.3 7.6 7.8]; %s
+ settings.motor.exp_thrust = [0 800 4000 5500 5160 5130 5400 5300 5450 5347 ...
+ 5160 4950 4700 4400 4400 3800 300 0]; %N
+
+ settings.m0 = 67.111; %kg Overall Mass (Burnout + Propellant)
+ settings.ms = 44.211; %kg Structural Mass (Burnout - Nosecone)
+ settings.mp = 17.7; %kg Propellant Mass
+ settings.tb = 7.60; %sec Burning Time
+ settings.mfr = settings.mp/settings.tb; %kg/s Mass Flow Rate
+ otherwise
+end
+
+
+% GEOMETRY DETAILS %
+% This parameters should be the same parameters set up in MISSILE DATCOM
+% simulation.
+
+settings.C = 0.174; %m Caliber (Fuselage Diameter)
+settings.S = 0.02378; %m2 Cross-sectional Surface
+L = 4.4; %m Rocket length
+
+% MASS GEOMERTY DETAILS %
+% x-axis: along the fuselage
+% y-axis: right wing
+% z-axis: downward
+
+% % Note: inerzias are used in "apogee_reached.m"
+% % HP: rocket inertias = full finite cilinder inertias
+% settings.Ixxf=settings.m0*(settings.C/2)^2/2; %Inertia to x-axis (Full)
+% settings.Ixxe=settings.ms*(settings.C/2)^2/2; %Inertia to x-axis (Empty)
+% settings.Iyyf=settings.m0.*((settings.C/2)^2/4 + L^2/3); %Inertia to y-axis (Full)
+% settings.Iyye=settings.ms.*((settings.C/2)^2/4 + L^2/3); %Inertia to y-axis (Empty)
+% settings.Izzf=settings.Iyyf; %Inertia to z-axis (Full)
+% settings.Izze=settings.Iyye; %Inertia to z-axis (Empty)
+
+
+% inertias first approximation
+settings.Ixxf = 0.14412762; %kg*m2 Inertia to x-axis (Full)
+settings.Ixxe = 0.14412762; %kg*m2 Inertia to x-axis (Empty)
+settings.Iyyf = 30.16891336; %kg*m2 Inertia to y-axis (Full)
+settings.Iyye = 30.16891335; %kg*m2 Inertia to y-axis (Empty)
+settings.Izzf = 30.16930792; %kg*m2 Inertia to z-axis (Full)
+settings.Izze = 30.16930792; %kg*m2 Inertia to z-axis (Empty)
+
+% AERODYNAMICS DETAILS %
+% This coefficients are intended to be obtained through MISSILE DATCOM
+% (than parsed with the tool datcom_parser.py)
+CoeffsF = load('for006_full.mat','Coeffs');
+settings.CoeffsF = CoeffsF.Coeffs;
+clear('CoeffsF');
+CoeffsE = load('for006_empty.mat','Coeffs');
+settings.CoeffsE = CoeffsE.Coeffs;
+clear('CoeffsE');
+
+% Note: All the parameters (AoA,Betas,Altitudes,Machs) must be the same for
+% empty and full configuration
+s = load('for006_full.mat','State');
+settings.Alphas = s.State.Alphas';
+settings.Betas = s.State.Betas';
+settings.Altitudes = s.State.Altitudes';
+settings.Machs = s.State.Machs';
+clear('s');
+
+%PARACHUTES DETAILS %
+%%% DROGUE 1 %%%
+settings.para1.S = 5; %m2 Surface
+settings.para1.mass = 0.577; %kg Parachute Mass
+settings.para1.CD = 0.8; %Parachute Drag Coefficient
+settings.para1.CL = 0; %Parachute Lift Coefficient
+
+%Altitude of Drogue 2 Opening
+settings.zdrg2 = 2000;
+
+%%% DROGUE 2 %%%
+settings.para2.S = 17.5; %m2 Surface
+settings.para2.mass = 1.140; %kg Parachute Mass
+settings.para2.CD = 0.59; %Parachute Drag Coefficient
+settings.para2.CL = 0; %Parachute Lift Coefficient
+
+%Altitude of Main Parachute Opening
+settings.zmain = 1000;
+
+%%% MAIN - ROGALLO %%%
+%The drogue parachute effects are neglected
+
+settings.para3.S = 15; %m2 Surface
+settings.para3.mass = 1.466; %kg Parachute Mass
+settings.para3.CD = 0.4; %Parachute Drag Coeff
+settings.para3.CL = 0.9; %Parachute Lift Coefficient
+
+% INTEGRATION OPTIONS %
+settings.ode.timeasc = 0:0.01:2000; %sec %Time span for ascend
+settings.ode.timedrg1 = 0:0.01:2000; %sec %Time span for drogue 1
+settings.ode.timedrg2 = 0:0.01:2000; %sec %Time span for drogue 2
+settings.ode.timemain = 0:0.01:2000; %sec %Time span for main (rogallo)
+settings.ode.timedesc = 0:0.01:2000; %sec %Time span for ballistic descent
+
+settings.ode.optionsasc = odeset('AbsTol',1E-3,'RelTol',1E-3,...
+ 'Events',@apogee,'InitialStep',1); %ODE options for ascend
+
+settings.ode.optionsdrg1 = odeset('AbsTol',1E-3,'RelTol',1E-3,...
+ 'Events',@drg2_opening); %ODE options for drogue
+
+settings.ode.optionsdrg2 = odeset('AbsTol',1E-3,'RelTol',1E-3,...
+ 'Events',@main_opening); %ODE options for drogue
+
+settings.ode.optionsmain = odeset('AbsTol',1E-3,'RelTol',1E-12,...
+ 'Events',@crash); %ODE options for ascend
+
+settings.ode.optionsdesc = odeset('AbsTol',1E-3,'RelTol',1E-12,...
+ 'Events',@crash); %ODE options for ballistic descent
+
+
+% WIND DETAILS %
+
+% Wind is randomly generated. Setting the same values for min and max will
+% fix the parameters of the wind.
+
+settings.wind.MagMin = 6; %Minimum Magnitude
+settings.wind.MagMax = 6; %Maximum Magnitude
+settings.wind.ElMin = 0*pi/180; %Minimum Elevation
+settings.wind.ElMax = 0*pi/180; %Maximum Elevation (Max == 90 Deg)
+settings.wind.AzMin = (180 + (135))*pi/180; %Minimum Azimuth
+settings.wind.AzMax = (180 + (135))*pi/180; %Maximum Azimuth
+
+% NOTE: wind aziumt angle = 180 + ...
+% means that I'm upwind with an angle of ...
+% 180+(..) heading sud-ovest
+% 180-(..) heading sud-est
+
+% Settings for the Wind Model
+settings.wind.Lat = 52.85; %Latitude of launching site
+settings.wind.Long = 16.033333; %Longitude of launching site
+settings.wind.Day = 180; %Day of the launch
+settings.wind.Seconds = 36000; %Second of the day
+
+% settings.wind.ww = vert_windgen(settings.wind.MagMin,settings.wind.MagMax);
+%Vertical wind speed
+settings.wind.ww = 0; % no vertical wind
+
+% BALLISTIC SIMULATION
+settings.ballistic = false; % Set to True to run a standard ballistic simulation
+
+% STOCHASTIC DETAILS %
+%If N>1 the stochastic routine is fired (different standard plots)
+settings.stoch.N = 1; % Number of iterations
+settings.stoch.parallel = false; % Using parallel or not parallel
+
+% PLOT DETAILS %
+settings.plot = false; % Set to True to Plot with default plots
+settings.tSteps = 250; % Set the number of time steps to visualize
+settings.DefaultFontSize = 10; % Default font size for plot
+settings.DefaultLineWidth = 1; % Default Line Width for plot
\ No newline at end of file
diff --git a/src/config_old.m b/src/config_old.m
new file mode 100644
index 0000000000000000000000000000000000000000..de14d2d4393e3ad2b57e2f9887664d904c9844a8
--- /dev/null
+++ b/src/config_old.m
@@ -0,0 +1,227 @@
+%CONFIG SCRIPT - This script sets up all the parameters for the simulation (H1 line)
+% All the parameters are stored in the "settings" structure.
+%
+
+% TODO: GUI to fill automatically this configuration script
+
+% Author: Ruben Di Battista
+% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu
+% email: ruben.dibattista@skywarder.eu
+% Website: http://www.skywarder.eu
+% License: 2-clause BSD
+
+% Author: Francesco Colombi
+% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu
+% email: francesco.colombi@skywarder.eu
+% Release date: 16/04/2016
+
+% clear all
+% close all
+% clc
+
+% LAUNCHPAD %
+settings.z0 = 100; %Launchpad Altitude
+settings.lrampa = 4.443; %LaunchPad route (launchpad lenght-distance from ground of the first hook)
+
+
+%STARTING ATTITUDE SETUP %
+settings.OMEGA = 85*pi/180; %Elevation Angle
+settings.PHI = 270*pi/180; %Azimuth Angle from North Direction
+
+% ENGINE DETAILS %
+
+% % Thrust on time profile (IV Order Interpolation)
+% settings.T_coeffs = [8034.5, 0, 0, 0, 0]; % constant thrust
+% data from the motor details sheet (experimental plot of the Thrust)
+
+% SYNTAX:
+% engine = 1 -> Cesaroni PRO 150 White Thunder
+% engine = 2 -> Cesaroni PRO 150 SkidMark
+% engine = 3 -> Cesaroni PRO 150 BlueStreak
+engine = 3;
+switch engine
+ case 1
+ % Cesaroni PRO 150 White Thunder
+ % Sampling for thrust interpolation
+ settings.motor.exp_time = [0 0.05 0.15 0.5 0.6 0.74 0.85 1.15 1.7 2.4 3 ...
+ 4 4.5 4.8 4.9 5 5.05 5.1 5.15 5.2]; %s
+ settings.motor.exp_thrust = [8605.1 8900 7900 8400 8400 8250 8200 8300 ...
+ 8400 8400 8200 7800 7600 7450 7350 7300 4500 500 100 0]; %N
+
+ settings.m0 = 61.8; %kg Overall Mass
+ settings.mp = 18.6; %kg Propellant Mass
+ settings.ms = settings.m0 - settings.mp; %kg Structural Mass
+ settings.tb = 5.12; %sec Burning Time
+ settings.mfr = settings.mp/settings.tb; %kg/s Mass Flow Rate
+ case 2
+ % Cesaroni PRO 150 SkidMark
+ % Sampling for thrust interpolation
+ settings.motor.exp_time = [0 0.1 0.2 0.3 0.4 0.5 0.6 0.8 1.2 1.8 3.2 ...
+ 3.6 4.8 6 7 7.2 7.6 7.8 7.9 8 8.1 8.19]; %s
+ settings.motor.exp_thrust = [0 3400 3100 3000 3300 3400 3500 3700 3700 ...
+ 3800 4000 4081.6 3900 3800 3700 3500 3350 3200 3000 2000 750 0]; %N
+
+ settings.m0 = 61.8; %kg Overall Mass
+ settings.mp = 17.157; %kg Propellant Mass
+ settings.ms = settings.m0 - settings.mp; %kg Structural Mass
+ settings.tb = 8.19; %sec Burning Time
+ settings.mfr = settings.mp/settings.tb; %kg/s Mass Flow Rate
+ case 3
+ % Cesaroni PRO 150 BlueStreak
+ % Sampling for thrust interpolation
+ settings.motor.exp_time = [0 0.06 0.1 0.15 0.25 0.45 0.8 1 2 3 ...
+ 4 5 6 6.8 7.05 7.3 7.6 7.8]; %s
+ settings.motor.exp_thrust = [0 800 4000 5500 5160 5130 5400 5300 5450 5347 ...
+ 5160 4950 4700 4400 4400 3800 300 0]; %N
+
+
+ settings.mp = 17.7; %kg Propellant Mass
+ settings.ms = 39 + 3.9598; %kg Structural Mass
+ settings.m0 = settings.mp + settings.ms; %kg Overall Mass without parachutes
+ %considering 3~5 kg DPL compartment
+ settings.tb = 7.60; %sec Burning Time
+ settings.mfr = settings.mp/settings.tb; %kg/s Mass Flow Rate
+ otherwise
+end
+
+
+% GEOMETRY DETAILS %
+% This parameters should be the same parameters set up in MISSILE DATCOM
+% simulation.
+
+settings.C = 0.174; %m Caliber (Fuselage Diameter)
+settings.S = 0.02378; %m2 Cross-sectional Surface
+L = 4.396; %m Rocket length
+
+% MASS GEOMERTY DETAILS %
+% x-axis: along the fuselage
+% y-axis: right wing
+% z-axis: downward
+
+% % Note: inerzias are used in "apogee_reached.m"
+% % HP: rocket inertias = full finite cilinder inertias
+% settings.Ixxf=settings.m0*(settings.C/2)^2/2; %Inertia to x-axis (Full)
+% settings.Ixxe=settings.ms*(settings.C/2)^2/2; %Inertia to x-axis (Empty)
+% settings.Iyyf=settings.m0.*((settings.C/2)^2/4 + L^2/3); %Inertia to y-axis (Full)
+% settings.Iyye=settings.ms.*((settings.C/2)^2/4 + L^2/3); %Inertia to y-axis (Empty)
+% settings.Izzf=settings.Iyyf; %Inertia to z-axis (Full)
+% settings.Izze=settings.Iyye; %Inertia to z-axis (Empty)
+
+
+% inertias first approximation
+settings.Ixxf = 0.14412762; %kg*m2 Inertia to x-axis (Full)
+settings.Ixxe = 0.14412762; %kg*m2 Inertia to x-axis (Empty)
+settings.Iyyf = 30.16891336; %kg*m2 Inertia to y-axis (Full)
+settings.Iyye = 30.16891335; %kg*m2 Inertia to y-axis (Empty)
+settings.Izzf = 30.16930792; %kg*m2 Inertia to z-axis (Full)
+settings.Izze = 30.16930792; %kg*m2 Inertia to z-axis (Empty)
+
+% AERODYNAMICS DETAILS %
+% This coefficients are intended to be obtained through MISSILE DATCOM
+% (than parsed with the tool datcom_parser.py)
+CoeffsF = load('for006_full.mat','Coeffs');
+settings.CoeffsF = CoeffsF.Coeffs;
+clear('CoeffsF');
+CoeffsE = load('for006_empty.mat','Coeffs');
+settings.CoeffsE = CoeffsE.Coeffs;
+clear('CoeffsE');
+
+% Note: All the parameters (AoA,Betas,Altitudes,Machs) must be the same for
+% empty and full configuration
+s = load('for006_full.mat','State');
+settings.Alphas = s.State.Alphas';
+settings.Betas = s.State.Betas';
+settings.Altitudes = s.State.Altitudes';
+settings.Machs = s.State.Machs';
+clear('s');
+
+%PARACHUTES DETAILS %
+%%% DROGUE 1 %%%
+settings.para1.S = 1.153; %m2 Surface
+settings.para1.mass = 0.0692; %kg Parachute Mass
+settings.para1.CD = 0.9; %Parachute Drag Coefficient
+settings.para1.CL = 0; %Parachute Lift Coefficient
+
+%Altitude of Drogue 2 Opening
+settings.zdrg2 = 2000;
+
+%%% DROGUE 2 %%%
+settings.para2.S = 11.520; %m2 Surface
+settings.para2.mass = 0.691; %kg Parachute Mass
+settings.para2.CD = 0.59; %Parachute Drag Coefficient
+settings.para2.CL = 0; %Parachute Lift Coefficient
+
+%Altitude of Main Parachute Opening
+settings.zmain = 1000;
+
+%%% MAIN - ROGALLO %%%
+%The drogue parachute effects are neglected
+
+settings.para3.S = 6.327; %m2 Surface
+settings.para3.mass = 0.380; %kg Parachute Mass
+settings.para3.CD = 0.4; %Parachute Drag Coeff
+settings.para3.CL = 0.9; %Parachute Lift Coefficient
+
+% INTEGRATION OPTIONS %
+settings.ode.timeasc = 0:0.01:2000; %sec %Time span for ascend
+settings.ode.timedrg1 = 0:0.01:2000; %sec %Time span for drogue 1
+settings.ode.timedrg2 = 0:0.01:2000; %sec %Time span for drogue 2
+settings.ode.timemain = 0:0.01:2000; %sec %Time span for main (rogallo)
+settings.ode.timedesc = 0:0.01:2000; %sec %Time span for ballistic descent
+
+settings.ode.optionsasc = odeset('AbsTol',1E-3,'RelTol',1E-3,...
+ 'Events',@apogee,'InitialStep',1); %ODE options for ascend
+
+settings.ode.optionsdrg1 = odeset('AbsTol',1E-3,'RelTol',1E-3,...
+ 'Events',@drg2_opening); %ODE options for drogue
+
+settings.ode.optionsdrg2 = odeset('AbsTol',1E-3,'RelTol',1E-3,...
+ 'Events',@main_opening); %ODE options for drogue
+
+settings.ode.optionsmain = odeset('AbsTol',1E-3,'RelTol',1E-12,...
+ 'Events',@crash); %ODE options for ascend
+
+settings.ode.optionsdesc = odeset('AbsTol',1E-3,'RelTol',1E-12,...
+ 'Events',@crash); %ODE options for ballistic descent
+
+
+% WIND DETAILS %
+
+% Wind is randomly generated. Setting the same values for min and max will
+% fix the parameters of the wind.
+
+settings.wind.MagMin = 0; %Minimum Magnitude
+settings.wind.MagMax = 0; %Maximum Magnitude
+settings.wind.ElMin = 0*pi/180; %Minimum Elevation
+settings.wind.ElMax = 0*pi/180; %Maximum Elevation (Max == 90 Deg)
+settings.wind.AzMin = (180 + (135))*pi/180; %Minimum Azimuth
+settings.wind.AzMax = (180 + (135))*pi/180; %Maximum Azimuth
+
+% NOTE: wind aziumt angle = 180 + ...
+% means that I'm upwind with an angle of ...
+% 180+(..) heading sud-ovest
+% 180-(..) heading sud-est
+
+% Settings for the Wind Model
+settings.wind.Lat = 52.85; %Latitude of launching site
+settings.wind.Long = 16.033333; %Longitude of launching site
+settings.wind.Day = 180; %Day of the launch
+settings.wind.Seconds = 36000; %Second of the day
+
+% settings.wind.ww = vert_windgen(settings.wind.MagMin,settings.wind.MagMax);
+%Vertical wind speed
+settings.wind.ww = 0; % no vertical wind
+
+% BALLISTIC SIMULATION
+settings.ballistic = false; % Set to True to run a ballistic simulation
+
+% STOCHASTIC DETAILS %
+%If N>1 the stochastic routine is fired (different standard plots)
+settings.stoch.N = 1; % Number of iterations
+settings.stoch.parallel = false; % Using parallel or not parallel
+
+% PLOT DETAILS %
+settings.plot = false; % Set to True to Plot with default plots
+settings.tSteps = 250; % Set the number of time steps to visualize
+settings.DefaultFontSize = 10; % Default font size for plot
+settings.DefaultLineWidth = 1; % Default Line Width for plot
\ No newline at end of file
diff --git a/src/crash.m b/src/crash.m
new file mode 100644
index 0000000000000000000000000000000000000000..4a0c7584c9ec21d036ac898314637162ee806e72
--- /dev/null
+++ b/src/crash.m
@@ -0,0 +1,22 @@
+function [ value,isterminal,direction ] = crash( t,Y,settings,varargin )
+%Event Function for ODE to stop simulation @ landing
+
+% Author: Ruben Di Battista
+% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu
+% email: ruben.dibattista@skywarder.eu
+% Website: http://www.skywarder.eu
+% April 2014; Last revision: 25.IV.2014
+% License: 2-clause BSD
+
+% Author: Francesco Colombi
+% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu
+% email: francesco.colombi@skywarder.eu
+% Release date: 16/04/2016
+
+value = Y(3);
+isterminal = 1;
+direction = 1;
+
+
+end
+
diff --git a/src/descent_parachute.m b/src/descent_parachute.m
new file mode 100644
index 0000000000000000000000000000000000000000..f8c1d44d7454c4c9d97ce01082701d27d0a6129a
--- /dev/null
+++ b/src/descent_parachute.m
@@ -0,0 +1,157 @@
+function [ dY ] = descent_parachute( t,Y,settings,uw,vw,ww,para )
+% ODEFUN for Parachute descent
+% State = ( x y z | u v w )
+
+% Author: Ruben Di Battista
+% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu
+% email: ruben.dibattista@skywarder.eu
+% Website: http://www.skywarder.eu
+% April 2014; Last revision: 25.IV.2014
+% License: 2-clause BSD
+
+% Author: Francesco Colombi
+% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu
+% email: francesco.colombi@skywarder.eu
+% Release date: 16/04/2016
+
+x = Y(1);
+y = Y(2);
+z = Y(3);
+u = Y(4);
+v = Y(5);
+w = Y(6);
+
+
+% global bool
+% % if bool = 0 the trends during the integration are NOT requested
+% % if bool = 1 the trends during the integration are saved and plotted
+
+
+% Adding Wind (supposed to be added in NED axes);
+
+% constant wind
+wind = [uw vw ww];
+
+% % Wind model
+% h = 0;
+% if -z > 0
+% h = -z+settings.z0;
+% end
+%
+% % global alt
+% % if bool == 1
+% % alt = [alt; h];
+% % end
+%
+% % wind in NED axis
+% [uw,vw] = atmoshwm07(settings.wind.Lat, settings.wind.Long, h, ...
+% 'day',settings.wind.Day,'model','quiet');
+% wind = [uw, vw, ww];
+%
+% % global WIND
+% % if bool == 1
+% % WIND = [WIND; uw, vw, ww];
+% % end
+
+% % Wind Model
+% [uw,vw,ww] = wind_generator(settings,z,t,Q);
+% wind = [ uw,vw,ww ];
+
+% Adding wind;
+ur = u - wind(1);
+vr = v - wind(2);
+wr = w - wind(3);
+
+V_norm = norm([ur vr wr]);
+
+% Constants
+switch para
+ case 1
+ S = settings.para1.S; %Parachute Surface
+ CD = settings.para1.CD; %Parachute CD
+ CL = settings.para1.CL; %Parachute CL
+ pmass = 0; %detached parachute mass
+ case 2
+ S = settings.para2.S;
+ CD = settings.para2.CD;
+ CL = settings.para2.CL;
+ pmass = settings.para1.mass;
+ case 3
+ S = settings.para3.S;
+ CD = settings.para3.CD;
+ CL = settings.para3.CL;
+ pmass = settings.para1.mass + settings.para2.mass;
+ otherwise
+end
+
+rad = sqrt(S/pi); %Equivalent parachute radius;
+
+%Side Surface
+Ss = pi*(rad/2)^2;
+
+g = 9.81;
+
+m = settings.ms - pmass;
+
+%Atmosphere
+if -z<0
+ z = 0;
+end
+
+[~, ~, ~, rho] = atmoscoesa(-z);
+
+%Center of Mass
+
+
+
+%Body Frame
+D = 0.5*rho*V_norm^2*S*CD;
+L = 0.5*rho*V_norm^2*S*CL;
+
+% HP: il paracadute � approssimato come una superficie rettangolare
+% orientata in modo da essere sempre normale alla direzione della velocita
+% relativa all'aria e sempre con un asse orizzontale e la forza Lift verso
+% l'alto
+% in modo da puntare ad atterrare nel punto [x,y] = (0m a nord,500m ad
+% ovest) rispetto la postazione di lancio
+
+
+tt = [ur vr wr]; % versore tangenziale
+kk = [-vr ur 0]; % versore orizzontale
+if (kk(1) == 0) && (kk(2) == 0)
+ kk = [-v u 0];
+end
+tt = tt/norm(tt);
+kk = -kk/norm(kk);
+nn = cross(tt, kk); nn = nn/norm(nn); % versore normale
+if (nn(3) > 0) % nn diretto verso il basso ?
+ nn = cross(kk, tt); nn = nn/norm(nn);
+end
+
+
+%Forces
+
+% X = D*ur/V_norm;
+% Y = D*vr/V_norm;
+% Z = D*wr/V_norm -m*g;
+
+% F = [-X,-Y,-Z]';
+F = -D*tt' + L*nn' + m*g*[0 0 1]';
+
+
+du=F(1)/m;
+dv=F(2)/m;
+dw=F(3)/m;
+
+
+
+
+
+dY(1:3) = [u v w]';
+dY(4) = du;
+dY(5) = dv;
+dY(6) = dw;
+
+dY = dY';
+
+end
\ No newline at end of file
diff --git a/src/drg2_opening.m b/src/drg2_opening.m
new file mode 100644
index 0000000000000000000000000000000000000000..9e7b31af1ec3afd0371a6bd3b5b0642bc2161ace
--- /dev/null
+++ b/src/drg2_opening.m
@@ -0,0 +1,15 @@
+function [ value,isterminal,direction ] = drg2_opening( t,Y,settings,varargin)
+%Event that sets the main parachute opening altitude
+
+% Author: Francesco Colombi
+% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu
+% email: francesco.colombi@skywarder.eu
+% Release date: 16/04/2016
+
+value = Y(3)+settings.zdrg2;
+isterminal = 1;
+direction = 1;
+
+
+end
+
diff --git a/src/interp4_easy.m b/src/interp4_easy.m
new file mode 100644
index 0000000000000000000000000000000000000000..ba3d18455f0159cfa8194f43904b7c6ac57ec782
--- /dev/null
+++ b/src/interp4_easy.m
@@ -0,0 +1,40 @@
+function [ V ] = interp4_easy( X1,X2,X3,X4,TABLE,x1,x2,x3,x4 )
+%windgen( AzMin,AzMax,ElMin,ElMax,MagMin,MagMax )
+% This function interpolate with nearest-neighbor method a R4->R function
+% F(x1...x4), given a [NxMxLxI] Matrix that discretize the function
+%
+%
+%
+% (X1...X4): column vectors of the variables the function depends on
+% (x1....x4): single values for each variable the interpolation is wanted
+% on
+% TABLE: the [NxMxLxI] matrix
+%
+%
+% length(X1) = N
+% length(X2) = M
+% length(X3) = L
+% length(X4) = I
+
+% Author: Ruben Di Battista
+% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu
+% email: ruben.dibattista@skywarder.eu
+% Website: http://www.skywarder.eu
+% April 2014; Last revision: 25.IV.2014
+% License: 2-clause BSD
+
+M = {X1 X2 X3 X4};
+m = [ x1 x2 x3 x4];
+
+index = zeros(4,1);
+
+for i=1:4
+ [~,index(i)] = min(abs(M{i} - m(i)));
+end
+
+V =TABLE(index(1),index(2),index(3),index(4));
+
+
+
+end
+
diff --git a/src/launch_pad.m b/src/launch_pad.m
new file mode 100644
index 0000000000000000000000000000000000000000..8928c14a0e9dd70e50295752649964f515651394
--- /dev/null
+++ b/src/launch_pad.m
@@ -0,0 +1,85 @@
+% Author: Francesco Colombi
+% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu
+% email: francesco.colombi@skywarder.eu
+% Release date: 16/04/2016
+
+close all
+clear, clc
+
+% inclinazione rampa
+theta = 85*pi/180;
+
+% [vv]horizon = R'*[vv]corpo
+% [vv]corpo = R*[vv]horizon
+R = [ cos(theta), 0, sin(theta); ...
+ 0, 1, 0; ...
+ -sin(theta), 0, cos(theta)];
+
+% velocit� di uscita dalla rampa
+Lrampa = 4.5; % m
+acc = 5*9.80665; % m/s2
+texit = sqrt(Lrampa/(.5*acc)); % sec
+V = acc*texit % m/s
+% V = 32;% m/s
+vv = V*[1, 0, 0]';
+
+% vento
+% incidenza vento
+dev = [-180:5:180]*pi/180;
+W = sqrt(2^2+5^2) * 1.2 % m/s
+% W = 10
+
+[~, a, ~, rho] = atmoscoesa(0);
+
+C = 0.174; % m Caliber (Fuselage Diameter)
+S = 0.02378; % m2 Cross-sectional Surface
+L = 4.32050; % m lunghezza missile
+
+load for006_full.mat
+% Coeffs. Interpolation
+givA = State.Alphas;
+givB = State.Betas;
+givH = State.Altitudes;
+givM = State.Machs;
+
+norm_vr = zeros(length(dev), 1);
+M = zeros(length(dev), 1);
+alpha = zeros(length(dev), 1);
+beta = zeros(length(dev), 1);
+XCP = zeros(length(dev), 1);
+CM = zeros(length(dev), 1);
+Mpitch = zeros(length(dev), 1);
+
+for n = 1:length(dev)
+ ww = -W*[cos(dev(n)), sin(dev(n)), 0]'; % assi horizon
+ ww = R*ww; % assi corpo
+
+ % vel relativa al vento
+ vvr = vv - ww; % assi corpo
+ norm_vr(n) = norm(vvr);
+ M(n) = norm_vr(n)/a;
+
+ ur = vvr(1);
+ vr = vvr(2);
+ wr = vvr(3);
+
+ % alpha
+ alpha(n) = atan(-wr/ur) *180/pi; % deg
+ beta(n) = asin(vr/norm(vvr)) *180/pi; % deg
+
+ XCP(n) = interpn(givA,givM,givB,givH,Coeffs.X_C_P,alpha(n),M(n),beta(n),0);
+ CM(n) = interpn(givA,givM,givB,givH,Coeffs.CM,alpha(n),M(n),beta(n),0);
+ Mpitch(n) = 0.5*rho*norm_vr(n)^2*S*C*CM(n);
+end
+
+figure()
+plot(dev*180/pi, XCP), grid on
+xlabel('angolo vento [deg]'), ylabel('XCP [ref length]')
+
+figure()
+plot(dev*180/pi, CM), grid on
+xlabel('angolo vento [deg]'), ylabel('CM [-]')
+
+figure()
+plot(dev*180/pi, Mpitch), grid on
+xlabel('angolo vento [deg]'), ylabel('M pitch [Nm]')
\ No newline at end of file
diff --git a/src/main_opening.m b/src/main_opening.m
new file mode 100644
index 0000000000000000000000000000000000000000..18fd536e62201512bbb5f3de4295140fa2ff733a
--- /dev/null
+++ b/src/main_opening.m
@@ -0,0 +1,16 @@
+function [ value,isterminal,direction ] = main_opening( t,Y,settings,varargin)
+%Event that sets the main parachute opening altitude
+
+
+% Author: Francesco Colombi
+% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu
+% email: francesco.colombi@skywarder.eu
+% Release date: 16/04/2016
+
+value = Y(3)+settings.zmain;
+isterminal = 1;
+direction = 1;
+
+
+end
+
diff --git a/src/parfor_progress.m b/src/parfor_progress.m
new file mode 100644
index 0000000000000000000000000000000000000000..d22b9f74211d725d6d3c8e99dad6f0d4715edf0a
--- /dev/null
+++ b/src/parfor_progress.m
@@ -0,0 +1,81 @@
+function percent = parfor_progress(N)
+%PARFOR_PROGRESS Progress monitor (progress bar) that works with parfor.
+% PARFOR_PROGRESS works by creating a file called parfor_progress.txt in
+% your working directory, and then keeping track of the parfor loop's
+% progress within that file. This workaround is necessary because parfor
+% workers cannot communicate with one another so there is no simple way
+% to know which iterations have finished and which haven't.
+%
+% PARFOR_PROGRESS(N) initializes the progress monitor for a set of N
+% upcoming calculations.
+%
+% PARFOR_PROGRESS updates the progress inside your parfor loop and
+% displays an updated progress bar.
+%
+% PARFOR_PROGRESS(0) deletes parfor_progress.txt and finalizes progress
+% bar.
+%
+% To suppress output from any of these functions, just ask for a return
+% variable from the function calls, like PERCENT = PARFOR_PROGRESS which
+% returns the percentage of completion.
+%
+% Example:
+%
+% N = 100;
+% parfor_progress(N);
+% parfor i=1:N
+% pause(rand); % Replace with real code
+% parfor_progress;
+% end
+% parfor_progress(0);
+%
+% See also PARFOR.
+
+% By Jeremy Scheff - jdscheff@gmail.com - http://www.jeremyscheff.com/
+
+error(nargchk(0, 1, nargin, 'struct'));
+
+if nargin < 1
+ N = -1;
+end
+
+percent = 0;
+w = 50; % Width of progress bar
+
+if N > 0
+ f = fopen('parfor_progress.txt', 'w');
+ if f<0
+ error('Do you have write permissions for %s?', pwd);
+ end
+ fprintf(f, '%d\n', N); % Save N at the top of progress.txt
+ fclose(f);
+
+ if nargout == 0
+ disp([' 0%[>', repmat(' ', 1, w), ']']);
+ end
+elseif N == 0
+ delete('parfor_progress.txt');
+ percent = 100;
+
+ if nargout == 0
+ disp([repmat(char(8), 1, (w+9)), char(10), '100%[', repmat('=', 1, w+1), ']']);
+ end
+else
+ if ~exist('parfor_progress.txt', 'file')
+ error('parfor_progress.txt not found. Run PARFOR_PROGRESS(N) before PARFOR_PROGRESS to initialize parfor_progress.txt.');
+ end
+
+ f = fopen('parfor_progress.txt', 'a');
+ fprintf(f, '1\n');
+ fclose(f);
+
+ f = fopen('parfor_progress.txt', 'r');
+ progress = fscanf(f, '%d');
+ fclose(f);
+ percent = (length(progress)-1)/progress(1)*100;
+
+ if nargout == 0
+ perc = sprintf('%3.0f%%', percent); % 4 characters wide, percentage
+ disp([repmat(char(8), 1, (w+9)), char(10), perc, '[', repmat('=', 1, round(percent*w/100)), '>', repmat(' ', 1, w - round(percent*w/100)), ']']);
+ end
+end
diff --git a/src/parfor_progress.txt b/src/parfor_progress.txt
new file mode 100644
index 0000000000000000000000000000000000000000..7d2156005f7e255337fab9500fb5cbc881bf6930
--- /dev/null
+++ b/src/parfor_progress.txt
@@ -0,0 +1,11 @@
+10
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
diff --git a/src/print_license.m b/src/print_license.m
new file mode 100644
index 0000000000000000000000000000000000000000..51c086d073d45824625d0c2c8edf9123c37d72ad
--- /dev/null
+++ b/src/print_license.m
@@ -0,0 +1,33 @@
+function [ ] = print_license()
+% Function that prints at first run the LICENSE saved in the LICENSE folder
+% in the same directory of this file
+
+% Author: Ruben Di Battista
+% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu
+% email: ruben.dibattista@skywarder.eu
+% Website: http://www.skywarder.eu
+% April 2014; Last revision: 25.IV.2014
+% License: 2-clause BSD
+
+if exist('LICENSE','dir')
+ if not(exist('LICENSE/run.txt','file'))
+ fh = fopen('LICENSE/LICENSE','r');
+ while not(feof(fh))
+ l = fgetl(fh);
+ fprintf(strcat(l,'\n'));
+ end
+ fclose(fh);
+ fh = fopen('LICENSE/run.txt','wt');
+ fprintf(fh,strcat('License accepted:',datestr(clock,0)));
+ fclose(fh);
+
+ agree = input('Write ''AGREE'' to accept the terms\n\n','s');
+ while not(strcmpi(agree,'agree'))
+ agree = input('Write ''AGREE'' to accept the terms\n\n','s');
+ end
+ end
+end
+
+
+end
+
diff --git a/src/start_simulation.m b/src/start_simulation.m
new file mode 100644
index 0000000000000000000000000000000000000000..e0bdddfa1c5ac6a7739f2c70957c3d8c9f6e7380
--- /dev/null
+++ b/src/start_simulation.m
@@ -0,0 +1,339 @@
+% Author: Francesco Colombi
+% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu
+% email: francesco.colombi@skywarder.eu
+% Release date: 16/04/2016
+
+close all
+clear, clc
+
+%% caricamento dati
+run('config.m');
+
+global bool contatore;
+global t_plot T_plot alpha_plot beta_plot M_plot CA_plot Drag_plot Forces_plot;
+% global WIND alt;
+
+
+% if bool = 0 the trends during the integration are NOT requested
+% if bool = 1 the trends during the integration are saved and plotted
+bool = 0;
+if bool == 1
+ contatore = 1;
+ t_plot = 0;
+ T_plot = 0;
+ alpha_plot = 0;
+ beta_plot = 0;
+ M_plot = 0;
+ CA_plot = 0;
+ Drag_plot = 0;
+ Forces_plot = 0;
+% WIND = [];
+% alt = [];
+end
+
+%% lancio simulatore
+[T,Y, Ta,Ya] = MAIN(settings);
+
+if (bool == 1)
+ figure();
+ plot(t_plot, T_plot, '.'), title('Thrust vs time'), grid on;
+ ylabel('Thrust [N]')
+ figure();
+ plot(t_plot, alpha_plot*180/pi), title('alpha vs time'), grid on;
+ ylabel('alpha [deg]')
+ figure();
+ plot(t_plot, beta_plot*180/pi), title('beta vs time'), grid on;
+ ylabel('beta [deg]')
+
+ figure();
+ plot(t_plot, M_plot), title('mach vs time'), grid on;
+ ylabel('Mach M [-]')
+
+ figure();
+ plot(t_plot, CA_plot), title('Aerodyn Coeff vs time'), grid on;
+ ylabel('Drag Coeff CD [-]')
+
+ figure();
+ plot(t_plot, Drag_plot), title('Drag vs time'), grid on;
+ ylabel('Drag D [N]')
+
+ figure();
+ plot(t_plot, Forces_plot), title('Axial Force vs time'), grid on;
+ ylabel('Axial force [N]')
+
+% figure()
+% title('Wind'), grid on, hold on;
+% plot(WIND(:, 1), alt)
+% plot(WIND(:, 2), alt)
+% plot(WIND(:, 3), alt)
+% xlabel('Wind magnitude [m/s]')
+% ylabel('Altitude [m]')
+% legend('North','Eas','Down')
+end
+
+
+[apogee, i_apogee] = max(-Y(:,3));
+
+%% posizione
+x = Y(:,1);
+y = Y(:,2);
+z = -Y(:,3);
+X = [x, y, z];
+
+%% velocit�
+N = length(x);
+% % derivata centrale
+% vx = (x(3:N)-x(1:N-2))./(T(3:N)-T(1:N-2));
+% vy = (y(3:N)-y(1:N-2))./(T(3:N)-T(1:N-2));
+% vz = (z(3:N)-z(1:N-2))./(T(3:N)-T(1:N-2));
+% % aggiungo derivata estremi
+% vx = [0; vx; (x(end)-x(end-1))/(T(end)-T(end-1))];
+% vy = [0; vy; (y(end)-y(end-1))/(T(end)-T(end-1))];
+% vz = [0; vz; (z(end)-z(end-1))/(T(end)-T(end-1))];
+vx = Y(:,4);
+vy = Y(:,5);
+vz = -Y(:,6);
+V = [vx, vy, vz];
+
+%% accelerazione
+% derivata centrale
+ax = (vx(3:N)-vx(1:N-2))./(T(3:N)-T(1:N-2));
+ay = (vy(3:N)-vy(1:N-2))./(T(3:N)-T(1:N-2));
+az = (vz(3:N)-vz(1:N-2))./(T(3:N)-T(1:N-2));
+% aggiungo derivata estremi
+ax = [vx(2)/T(2); ax; (vx(end)-vx(end-1))/(T(end)-T(end-1))];
+ay = [vy(2)/T(2); ay; (vy(end)-vy(end-1))/(T(end)-T(end-1))];
+az = [vz(2)/T(2); az; (vz(end)-vz(end-1))/(T(end)-T(end-1))];
+
+% % derivata in avanti
+% ax = (vx(2:N)-vx(1:N-1))./(T(2:N)-T(1:N-1));
+% ay = (vy(2:N)-vy(1:N-1))./(T(2:N)-T(1:N-1));
+% az = (vz(2:N)-vz(1:N-1))./(T(2:N)-T(1:N-1));
+% % aggiungo derivata estremi
+% ax = [0; ax];
+% ay = [0; ay];
+% az = [0; az];
+
+A = [ax, ay, az];
+
+clear ax ay az
+clear vx vy
+
+mod_X = zeros(length(X), 1);
+mod_V = mod_X;
+mod_A = mod_X;
+
+for k = 1:length(X)
+ mod_X(k) = norm(X(k,:));
+ mod_V(k) = norm(V(k,:));
+ mod_A(k) = norm(A(k,:));
+end
+[max_dist, imax_dist] = max(mod_X);
+[max_v, imax_v] = max(mod_V);
+[max_a, imax_a] = max(mod_A);
+
+M = zeros(length(X), 1);
+Tamb = M;
+Ttot = Tamb;
+for n = 1:length(M)
+ h = z(n);
+ if (h < 0)
+ h = 0;
+ end
+ [Tamb(n), a, ~, ~] = atmoscoesa(h);
+ M(n) = mod_V(n)/a;
+ Ttot(n) = Tamb(n)*(1+0.2*M(n)^2);
+end
+
+[max_M, imax_M] = max(M);
+
+%% launch pad exit index
+iexit = find(mod_X <= settings.lrampa);
+iexit = iexit(end);
+
+
+%% DATA RECORD (display)
+disp(' ')
+disp('DATA RECORD:')
+fprintf('apogee reached: %g [m] \n', apogee);
+fprintf('apogee velocity: %g [m/s] \n', mod_V(i_apogee));
+fprintf('time: %g [sec] \n\n', Ta(end))
+
+fprintf('max speed reached: %g [m/s] \n', max_v)
+fprintf('altitude: %g [m] \n', z(imax_v))
+fprintf('Mach: %g [-] \n', M(imax_v))
+fprintf('time: %g [sec] \n\n', T(imax_v))
+
+fprintf('max Mach reached: %g [-] \n', max_M)
+fprintf('altitude: %g [m] \n', z(imax_M))
+fprintf('veocity: %g [m/s] \n', mod_V(imax_M))
+fprintf('time: %g [sec] \n\n', T(imax_M))
+
+fprintf('max acceleration reached: %g [m/s2] = %g [g] \n', max_a, max_a/9.80665)
+fprintf('veocity: %g [m/s] \n', mod_V(imax_a))
+fprintf('time: %g [sec] \n\n', T(imax_a))
+
+fprintf('run on launch pad: %g [m] \n', mod_X(iexit))
+fprintf('speed at launch pad exit: %g [m/s] \n', mod_V(iexit))
+fprintf('time: %g [sec] \n\n', T(iexit))
+
+%% plot some interesting stuff
+plots = true;
+if plots
+ % figure(), hold on, grid on;
+ % plot(T, z);
+ % plot(T, vz);
+ % plot(T, mod_V);
+ % plot(T, mod_A/9.80665)
+ % legend('z', 'vz', 'mod V', 'mod A');
+ % xlabel('Time [sec]');
+ % ylabel('z [m], vz [m/s], mod V [m/s], mod A [g]')
+
+ figure();
+ subplot(3,1,1)
+ plot(T, z), grid on, xlabel('time [sec]'), ylabel('altitude [m]');
+ subplot(3,1,2)
+ plot(T, vz), grid on, hold on;
+ plot(T, mod_V);
+ xlabel('time [sec]'), ylabel('vz [m/s], |V| [m/s]');
+ legend('vz', '|V|');
+ subplot(3,1,3)
+ plot(T, mod_A/9.80665), grid on;
+ xlabel('time [sec]'), ylabel('|A| [g]');
+ legend('|A|');
+
+ figure();
+ plot(T, z), grid on, xlabel('time [sec]'), ylabel('altitude [m]');
+
+ figure();
+ plot(T, vz), grid on, hold on;
+ plot(T, mod_V);
+ xlabel('time [sec]'), ylabel('vz [m/s], |V| [m/s]');
+ legend('vz', '|V|');
+
+ figure();
+ plot(T, mod_A/9.80665), grid on;
+ xlabel('time [sec]'), ylabel('|A| [g]');
+ legend('|A|');
+
+ figure();
+ plot3(y, x, z), axis equal, hold on, grid on;
+ plot3(Ya(end,2), Ya(end,1), -Ya(end,3), '*')
+ plot3(0, 0, 0, 'or')
+ % visualizzazione circonferenze distanze
+ theta_plot = linspace(0,2*pi);
+ R_plot = [1, 2, 3, 4, 5]*1000;
+ for j = 1:length(R_plot)
+ x_plot = R_plot(j)*cos(theta_plot');
+ y_plot = R_plot(j)*sin(theta_plot');
+ z_plot = zeros(length(theta_plot), 1);
+ plot3(y_plot, x_plot, z_plot, '--r')
+ end
+
+ title('Trajectory')
+ xlabel('y, East [m]'), ylabel('x, North [m]'), zlabel('Altitude [m]')
+
+ figure();
+ plot(y, x), axis equal, hold on, grid on;
+ % visualizzazione circonferenze distanze
+ theta_plot = linspace(0,2*pi);
+ R_plot = [1, 2, 3, 4, 5]*1000;
+ for j = 1:length(R_plot)
+ x_plot = R_plot(j)*cos(theta_plot');
+ y_plot = R_plot(j)*sin(theta_plot');
+ plot(y_plot, x_plot, '--r')
+ end
+ xlabel('y, East [m]'), ylabel('x, North [m]');
+
+ % lengths to 0 if they are less than 1 meter
+ x_flight=x;
+ y_flight=y;
+ z_flight=z;
+ for ii = 1 : length(x)
+ if norm(x_flight(ii))<1
+ x_flight(ii) = 0;
+ end
+ if norm(y_flight(ii))<1
+ y_flight(ii) = 0;
+ end
+ if norm(z_flight(ii))<1
+ z_flight(ii) = 0;
+ end
+ end
+
+ figure();
+ subplot(1,3,1);
+ plot(y_flight/1000, x_flight/1000), axis equal, hold on, grid on;
+ plot(Ya(end,2)/1000, Ya(end,1)/1000, '*');
+ plot(0, 0, 'or')
+ xlabel('y, East [Km]'), ylabel('x, North [Km]');
+ subplot(1,3,2);
+ plot(x_flight/1000, z_flight/1000), hold on, grid on;
+ plot(Ya(end,1)/1000, -Ya(end,3)/1000, '*');
+ plot(0, 0, 'or')
+ xlabel('x, North [Km]'), ylabel('z, Altitude [Km]');
+ %setting limit if is parallel to east
+ if sum(x_flight) == 0
+ xlim([-1 1]);
+ end
+ subplot(1,3,3);
+ plot(y_flight/1000, z_flight/1000), hold on, grid on;
+ plot(Ya(end,2)/1000, -Ya(end,3)/1000, '*');
+ plot(0, 0, 'or')
+ xlabel('y, East [Km]'), ylabel('z, Altitude [Km]');
+ %setting limit if is parallel to north
+ if sum(y_flight) == 0
+ xlim([-1 1]);
+ end
+
+ % angular rates
+ figure();
+ plot(Ta, Ya(:, 7)*180/pi)
+ hold on, grid on
+ plot(Ta, Ya(:, 8)*180/pi)
+ plot(Ta, Ya(:, 9)*180/pi)
+ xlabel('time [sec]'), ylabel('p, q, r [grad/sec]')
+ legend('p: roll rate', 'q: pitch rate', 'r: yaw rate')
+ title('angular rates vs time')
+
+ figure();
+ subplot(3,1,1)
+ plot(Ta, Ya(:, 8)*180/pi), grid on;
+ xlabel('time [sec]'), ylabel('pitch rate q [grad/sec]')
+ subplot(3,1,2)
+ plot(Ta, Ya(:, 9)*180/pi), grid on;
+ xlabel('time [sec]'), ylabel('yaw rate r [grad/sec]')
+ subplot(3,1,3)
+ plot(Ta, Ya(:, 7)*180/pi), grid on;
+ xlabel('time [sec]'), ylabel('roll rate p [grad/sec]')
+
+
+ % angle
+ alpha = zeros(length(Ta),1);
+ beta = zeros(length(Ta),1);
+ roll = zeros(length(Ta),1);
+
+ for k = 2:length(Ta)
+ alpha(k) = alpha(k-1) + (Ya(k, 8) + Ya(k-1, 8))/2*180/pi*(T(k)-T(k-1));
+ beta(k) = beta(k-1) + (Ya(k, 9) + Ya(k-1, 9))/2*180/pi*(T(k)-T(k-1));
+ roll(k) = roll(k-1) + (Ya(k, 7) + Ya(k-1, 7))/2*180/pi*(T(k)-T(k-1));
+ end
+ figure(), plot(Ta, alpha+settings.OMEGA*180/pi), title('Alpha vs time'), grid on;
+ figure(), plot(Ta, beta), title('Beta vs time'), grid on;
+ figure(), plot(Ta, roll), title('Roll vs time'), grid on;
+
+
+ % Mach
+ figure(), plot(T, M), grid on;
+ xlabel('time [sec]'), ylabel('Mach [-]');
+
+ % Stagnation Temperature
+ figure();
+ plot(T, Tamb-273.15);
+ hold on, grid on;
+ plot(T, Ttot-273.15)
+ title('Temperature profile')
+ xlabel('time [sec]'), ylabel('Temperature [�C]');
+ legend('Surrounding', 'Total', 'location', 'best')
+
+end
diff --git a/src/std_run.m b/src/std_run.m
new file mode 100644
index 0000000000000000000000000000000000000000..5cbe2d875fa5ff60147fd95ab32e7c9acbff3afa
--- /dev/null
+++ b/src/std_run.m
@@ -0,0 +1,278 @@
+function [ Tf,Yf, Ta,Ya ] = std_run( settings )
+%STD RUN - This function runs a standard (non-stochastic) simulation
+% OUTPUTS
+% Tf: Time steps
+% Yf: Final State
+
+% Author: Ruben Di Battista
+% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu
+% email: ruben.dibattista@skywarder.eu
+% Website: http://www.skywarder.eu
+% April 2014; Last revision: 31.XII.2014
+% License: 2-clause BSD
+
+% Author: Francesco Colombi
+% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu
+% email: francesco.colombi@skywarder.eu
+% Release date: 16/04/2016
+
+%Starting Attitude
+Q0 = angle2quat(settings.PHI,settings.OMEGA,0*pi/180,'ZYX')';
+
+%Starting State
+X0 = [0 0 0]';
+V0 = [0 0 0]';
+W0 = [0 0 0]';
+X0a = [X0;V0;W0;Q0;settings.m0;settings.Ixxf;settings.Iyyf;settings.Izzf];
+
+%% Wind Generation
+[uw,vw,ww] = windgen(settings.wind.AzMin,settings.wind.AzMax,...
+ settings.wind.ElMin,settings.wind.ElMax,settings.wind.MagMin,...
+ settings.wind.MagMax);
+
+%% ASCEND %%
+
+[Ta,Ya] = ode45(@ascend,settings.ode.timeasc,X0a,settings.ode.optionsasc,...
+ settings,uw,vw,ww);
+
+%% DROGUE 1 %%
+para = 1; %Flag for Drogue 1
+
+%Initial Condition are the last from ascend (need to rotate because
+%velocities are outputted in body axes)
+X0d1 = [Ya(end,1:3) quatrotate(quatconj(Ya(end,10:13)),Ya(end,4:6))];
+[Td1,Yd1] = ode45(@descent_parachute,settings.ode.timedrg1,X0d1,...
+ settings.ode.optionsdrg1,settings,uw,vw,ww,para);
+
+%% DROGUE 2 %%
+para = 2; %Flag for Drogue 2
+
+%Initial Condition are the last from drogue 1 descent
+X0d2 = Yd1(end,:);
+[Td2,Yd2] = ode45(@descent_parachute,settings.ode.timedrg2,X0d2,...
+ settings.ode.optionsdrg2,settings,uw,vw,ww,para);
+
+%% MAIN %%
+para = 3; %Flag for Main (Rogall)
+
+%Initial Condition are the last from drogue 2 descent
+X0m = Yd2(end,:);
+[Trog,Yrog] = ode45(@descent_parachute,settings.ode.timemain,X0m,...
+ settings.ode.optionsmain,settings,uw,vw,ww,para);
+
+
+
+%% FINAL STATE ASSEMBLING %%
+
+%Total State
+Yf = [Ya(:,1:3) quatrotate(quatconj(Ya(:,10:13)),Ya(:,4:6));Yd1; Yd2;Yrog];
+
+global Par1 Par2 Par3 Asc
+Asc = [Ya(:,1:3) quatrotate(quatconj(Ya(:,10:13)),Ya(:,4:6))];
+Par1 = Yd1;
+Par2 = Yd2;
+Par3 = Yrog;
+%Total Time
+Tf = [Ta; Ta(end)+Td1; Ta(end)+Td1(end)+Td2;Ta(end)+Td1(end)+Td2(end)+Trog];
+
+%Parachute State
+Yp = [Yd1;Yd2;Yrog];
+Tp = [Ta(end)+Td1;Ta(end)+Td1(end)+Td2;Ta(end)+Td1(end)+Td2(end)+Trog];
+
+
+%% PLOTTING THINGS %%
+
+if settings.plot == 1
+
+ set(0,'DefaultAxesFontSize',settings.DefaultFontSize,...
+ 'DefaultLineLineWidth',settings.DefaultLineWidth);
+
+ % ASCENT %
+ plot(Tf,-Yf(:,3)+settings.z0,'k-','LineWidth',2)
+ hold on
+ title('Altitude Profile on Time');
+ xlabel('Time [s]')
+ ylabel('Altitude [m]');
+ h1=plot(Ta(end),-Ya(end,3)+settings.z0,'ko','MarkerSize',8,...
+ 'MarkerFaceColor','k');
+ h2=plot(Ta(end)+Td1(end),-Yd1(end,3)+settings.z0,'ks','MarkerSize',8,...
+ 'MarkerFaceColor','k');
+ h3=plot(Ta(end)+Td1(end)+Td2(end),-Yd2(end,3)+settings.z0,'ks','MarkerSize',8,...
+ 'MarkerFaceColor','k');
+ legend([h1 h2 h3],'Apogee/1st Drogue Deployment',...
+ '2nd Drogue Deployment', 'Main Parachute Deployment')
+ grid on
+
+ %Interpolation for less points --> visualization
+ Tinterp = linspace(Tf(1),Tf(end),settings.tSteps);
+ [Tunique,ia,~]= unique(Tf);
+ u = interp1(Tunique,Yf(ia,4),Tinterp);
+ v = interp1(Tunique,Yf(ia,5),Tinterp);
+ w = interp1(Tunique,Yf(ia,6),Tinterp);
+
+ % VELOCITIES %
+ Va=quatrotate(quatconj(Ya(end,10:13)),Ya(end,4:6)); %apogee
+
+ figure;
+ suptitle('Velocities Profiles on Time')
+ subplot(3,1,1);
+ plot(Tinterp,u,'k-')
+ hold on
+ xlabel('Time[s]');
+ ylabel('Velocity-x [m/s]');
+ h1=plot(Ta(end),Va(1),'ko','MarkerSize',8,...
+ 'MarkerFaceColor','k');
+ h2=plot(Ta(end)+Td1(end),Yd1(end,4),'ks','MarkerSize',8,...
+ 'MarkerFaceColor','k');
+ h3=plot(Ta(end)+Td1(end)+Td2(end),Yd2(end,4),'ks','MarkerSize',8,...
+ 'MarkerFaceColor','k');
+ legend([h1,h2,h3],'Apogee/1st Drogue Deployment',...
+ '2nd Drogue Deployment','Main Parachute Deployment',...
+ 'Location','southeast');
+ grid on
+
+ subplot(3,1,2)
+ plot(Tinterp,v,'k-')
+ hold on
+ xlabel('Time[s]');
+ ylabel('Velocity-y [m/s]');
+ h1=plot(Ta(end),Va(2),'ko','MarkerSize',8,...
+ 'MarkerFaceColor','k');
+ h2=plot(Ta(end)+Td1(end),Yd1(end,5),'ks','MarkerSize',8,...
+ 'MarkerFaceColor','k');
+ h3=plot(Ta(end)+Td1(end)+Td2(end),Yd2(end,5),'ks','MarkerSize',8,...
+ 'MarkerFaceColor','k');
+ legend([h1,h2,h3],'Apogee/1st Drogue Deployment',...
+ '2nd Drogue Deployment','Main Parachute Deployment',...
+ 'Location','southeast');
+ grid on
+
+
+ subplot(3,1,3)
+ plot(Tinterp,w,'k-')
+ hold on
+ xlabel('Time[s]');
+ ylabel('Velocity-z [m/s]');
+ h1=plot(Ta(end),Va(3),'ko','MarkerSize',8,...
+ 'MarkerFaceColor','k');
+ h2=plot(Ta(end)+Td1(end),Yd1(end,6),'ks','MarkerSize',8,...
+ 'MarkerFaceColor','k');
+ h3=plot(Ta(end)+Td1(end)+Td2(end),Yd2(end,6),'ks','MarkerSize',8,...
+ 'MarkerFaceColor','k');
+ legend([h1,h2,h3],'Apogee/1st Drogue Deployment',...
+ '2nd Drogue Deployment','Main Parachute Deployment',...
+ 'Location','southeast');
+ grid on
+
+ %COMPLETE TRAJECTORY%
+ figure;
+ h0=plot3(Yf(1,2),Yf(1,1),-Yf(1,3)+settings.z0,'k+','MarkerSize',10);
+ hold on
+ plot3(Yf(:,2),Yf(:,1),-Yf(:,3)+settings.z0,'k-','LineWidth',2)
+ title('Complete Trajectory');
+ xlabel('East [m]')
+ ylabel('Noth [m]');
+ zlabel('Altitude [m]');
+ h1=plot3(Ya(end,2),Ya(end,1),-Ya(end,3)+settings.z0,'ko','MarkerSize',8,...
+ 'MarkerFaceColor','k');
+ h2=plot3(Yd1(end,2),Yd1(end,1),-Yd1(end,3)+settings.z0,'ks','MarkerSize',8,...
+ 'MarkerFaceColor','k');
+ h3=plot3(Yd2(end,2),Yd2(end,1),-Yd2(end,3)+settings.z0,'ks','MarkerSize',8,...
+ 'MarkerFaceColor','k');
+ legend([h0,h1,h2,h3],'Launch Pad','Apogee/1st Drogue Deployment',...
+ '2nd Drogue Deployment','Main Parachute Deployment');
+ grid on
+
+ % PLANAR DISPLACEMENT %
+
+ figure
+ plot(Yf(:,2),Yf(:,1),'k-','LineWidth',2)
+ hold on
+ title('Planar Displacement');
+ xlabel('East [m]')
+ ylabel('North [m]');
+ h1=plot(Ya(end,2),Ya(end,1),'ko','MarkerSize',8,...
+ 'MarkerFaceColor','k');
+ h2=plot(Yd1(end,2),Yd1(end,1),'ks','MarkerSize',8,...
+ 'MarkerFaceColor','k');
+ h3=plot(Yd2(end,2),Yd2(end,1),'ks','MarkerSize',8,...
+ 'MarkerFaceColor','k');
+ legend([h1,h2,h3],'Apogee/1st Drogue Deployment',...
+ '2nd Drogue Deployment','Main Parachute Deployment',...
+ 'Location','southeast');
+ grid on
+
+ % PARACHUTES %
+ figure
+ plot(Tp,-Yp(:,3)+settings.z0,'k-','LineWidth',2)
+ hold on
+ title('Altitude Profile on Time (parachutes)');
+ xlabel('Time [s]')
+ ylabel('Altitude [m]');
+ h1=plot(Ta(end),-Ya(end,3)+settings.z0,'ko','MarkerSize',8,...
+ 'MarkerFaceColor','k');
+ h2=plot(Ta(end)+Td1(end),-Yd1(end,3)+settings.z0,'ks','MarkerSize',8,...
+ 'MarkerFaceColor','k');
+ h3=plot(Ta(end)+Td1(end)+Td2(end),-Yd2(end,3)+settings.z0,'ks','MarkerSize',8,...
+ 'MarkerFaceColor','k');
+ legend([h1,h2,h3],'Apogee/1st Drogue Deployment',...
+ '2nd Drogue Deployment','Main Parachute Deployment',...
+ 'Location','southeast');
+ grid on
+
+
+
+ %Interpolation for less points --> visualization
+ Tinterp = linspace(Ta(1),Ta(end),settings.tSteps);
+ p = interp1(Ta,Ya(:,7),Tinterp);
+ q = interp1(Ta,Ya(:,8),Tinterp);
+ r = interp1(Ta,Ya(:,9),Tinterp);
+
+ % ANGULAR RATES %
+
+ figure;
+ suptitle('Angular rates on Time')
+ subplot(3,1,1);
+ plot(Tinterp,p,'k-')
+ hold on
+ xlabel('Time[s]');
+ ylabel('Roll rate p [rad/s]');
+ h1=plot(Ta(end),Ya(end,7),'ko','MarkerSize',8,...
+ 'MarkerFaceColor','k');
+ legend(h1,'Apogee/1st Drogue Deployment','Location','southeast');
+ grid on
+
+ subplot(3,1,2)
+ plot(Tinterp,q,'k-')
+ hold on
+ xlabel('Time[s]');
+ ylabel('Pitch rate q [rad/s]');
+ h1=plot(Ta(end),Ya(end,8),'ko','MarkerSize',8,...
+ 'MarkerFaceColor','k');
+ legend(h1,'Apogee/1st Drogue Deployment','Location','northeast');
+ grid on
+
+
+ subplot(3,1,3)
+ plot(Tinterp,r,'k-')
+ hold on
+ xlabel('Time[s]');
+ ylabel('Yaw rate r [rad/s]');
+ h1=plot(Ta(end),Ya(end,9),'ko','MarkerSize',8,...
+ 'MarkerFaceColor','k');
+ legend(h1,'Apogee/1st Drogue Deployment','Location','southeast');
+ grid on
+
+end
+
+%Resizing
+h = get(0,'children');
+scrsz = get(0,'ScreenSize');
+for i=1:length(h)
+ set(h(i),'OuterPosition',[0 0 scrsz(4) scrsz(4)])
+ %saveas(h(i), ['figure' num2str(i)], 'fig');
+end
+
+
+
+end
\ No newline at end of file
diff --git a/src/std_run_ballistic.m b/src/std_run_ballistic.m
new file mode 100644
index 0000000000000000000000000000000000000000..7d1cb4e3d254a217a85d36900e397b8ae82c424c
--- /dev/null
+++ b/src/std_run_ballistic.m
@@ -0,0 +1,229 @@
+function [ Tf,Yf, Ta,Ya ] = std_run_test( settings )
+%STD RUN - This function runs a standard (non-stochastic) simulation
+% OUTPUTS
+% Tf: Time steps
+% Yf: Final State
+
+% Author: Ruben Di Battista
+% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu
+% email: ruben.dibattista@skywarder.eu
+% Website: http://www.skywarder.eu
+% April 2014; Last revision: 31.XII.2014
+% License: 2-clause BSD
+
+% Author: Francesco Colombi
+% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu
+% email: francesco.colombi@skywarder.eu
+% Release date: 16/04/2016
+
+%Starting Attitude
+Q0 = angle2quat(settings.PHI,settings.OMEGA,0*pi/180,'ZYX')';
+
+%Starting State
+X0 = [0 0 0]';
+V0 = [0 0 0]';
+W0 = [0 0 0]';
+X0a = [X0;V0;W0;Q0;settings.m0;settings.Ixxf;settings.Iyyf;settings.Izzf];
+
+% Wind Generation
+[uw,vw,ww] = windgen(settings.wind.AzMin,settings.wind.AzMax,...
+ settings.wind.ElMin,settings.wind.ElMax,settings.wind.MagMin,...
+ settings.wind.MagMax);
+
+%% ASCEND %%
+
+[Ta,Ya] = ode45(@ascend,settings.ode.timeasc,X0a,settings.ode.optionsasc,...
+ settings,uw,vw,ww);
+
+%% DESCEND %%
+
+[Td,Yd] = ode45(@ballistic_descent,settings.ode.timedesc,Ya(end,:),settings.ode.optionsdesc,...
+ settings,uw,vw,ww);
+
+% %% DROGUE 1 %%
+% para = 1; %Flag for Drogue 1
+%
+% %Initial Condition are the last from ascend (need to rotate because
+% %velocities are outputted in body axes)
+% X0d1 = [Ya(end,1:3) quatrotate(quatconj(Ya(end,10:13)),Ya(end,4:6))];
+% [Td1,Yd1] = ode45(@descent_parachute,settings.ode.timedrg1,X0d1,...
+% settings.ode.optionsdrg1,settings,uw,vw,ww,para);
+%
+% %% DROGUE 2 %%
+% para = 2; %Flag for Drogue 2
+%
+% %Initial Condition are the last from drogue 1 descent
+% X0d2 = Yd1(end,:);
+% [Td2,Yd2] = ode45(@descent_parachute,settings.ode.timedrg2,X0d2,...
+% settings.ode.optionsdrg2,settings,uw,vw,ww,para);
+%
+% %% MAIN %%
+% para = 3; %Flag for Main (Rogall)
+%
+% %Initial Condition are the last from drogue 2 descent
+% X0m = Yd2(end,:);
+% [Trog,Yrog] = ode45(@descent_parachute,settings.ode.timemain,X0m,...
+% settings.ode.optionsmain,settings,uw,vw,ww,para);
+
+%% FINAL STATE ASSEMBLING %%
+
+%Total State
+Yf = [Ya;Yd];
+%Total Time
+Tf = [Ta; Ta(end)+Td];
+
+%% PLOTTING THINGS %%
+
+if settings.plot == 1
+
+ set(0,'DefaultAxesFontSize',settings.DefaultFontSize,...
+ 'DefaultLineLineWidth',settings.DefaultLineWidth);
+
+ %Interpolation for less points --> visualization
+ Tinterp = linspace(Tf(1),Tf(end),settings.tSteps);
+ [Tunique,ia,~]= unique(Tf);
+ u = interp1(Tunique,Yf(ia,4),Tinterp);
+ v = interp1(Tunique,Yf(ia,5),Tinterp);
+ w = interp1(Tunique,Yf(ia,6),Tinterp);
+
+ % VELOCITIES %
+ Va=quatrotate(quatconj(Ya(end,10:13)),Ya(end,4:6)); %apogee
+
+ figure;
+ suptitle('Velocities Profiles on Time')
+ subplot(3,1,1);
+ plot(Tinterp,u,'k-')
+ hold on
+ xlabel('Time[s]');
+ ylabel('Velocity-x [m/s]');
+ h1=plot(Ta(end),Va(1),'ko','MarkerSize',8,...
+ 'MarkerFaceColor','k');
+ h2=plot(Ta(end)+Td(end),Yd(end,4),'ks','MarkerSize',8,...
+ 'MarkerFaceColor','k');
+ legend([h1,h2],'Apogee/1st Drogue Deployment',...
+ '2nd Drogue Deployment',...
+ 'Location','southeast');
+ grid on
+
+ subplot(3,1,2)
+ plot(Tinterp,v,'k-')
+ hold on
+ xlabel('Time[s]');
+ ylabel('Velocity-y [m/s]');
+ h1=plot(Ta(end),Va(2),'ko','MarkerSize',8,...
+ 'MarkerFaceColor','k');
+ h2=plot(Ta(end)+Td(end),Yd(end,5),'ks','MarkerSize',8,...
+ 'MarkerFaceColor','k');
+ legend([h1,h2],'Apogee/1st Drogue Deployment',...
+ '2nd Drogue Deployment',...
+ 'Location','southeast');
+ grid on
+
+
+ subplot(3,1,3)
+ plot(Tinterp,w,'k-')
+ hold on
+ xlabel('Time[s]');
+ ylabel('Velocity-z [m/s]');
+ h1=plot(Ta(end),Va(3),'ko','MarkerSize',8,...
+ 'MarkerFaceColor','k');
+ h2=plot(Ta(end)+Td(end),Yd(end,6),'ks','MarkerSize',8,...
+ 'MarkerFaceColor','k');
+ legend([h1,h2],'Apogee/1st Drogue Deployment',...
+ '2nd Drogue Deployment',...
+ 'Location','southeast');
+ grid on
+
+ %COMPLETE TRAJECTORY%
+ figure;
+ h0=plot3(Yf(1,2),Yf(1,1),-Yf(1,3)+settings.z0,'k+','MarkerSize',10);
+ hold on
+ plot3(Yf(:,2),Yf(:,1),-Yf(:,3)+settings.z0,'k-','LineWidth',2)
+ title('Complete Trajectory');
+ xlabel('East [m]')
+ ylabel('Noth [m]');
+ zlabel('Altitude [m]');
+ h1=plot3(Ya(end,2),Ya(end,1),-Ya(end,3)+settings.z0,'ko','MarkerSize',8,...
+ 'MarkerFaceColor','k');
+ h2=plot3(Yd(end,2),Yd(end,1),-Yd(end,3)+settings.z0,'ks','MarkerSize',8,...
+ 'MarkerFaceColor','k');
+
+ legend([h0,h1,h2],'Launch Pad','Apogee/1st Drogue Deployment',...
+ '2nd Drogue Deployment','Main Parachute Deployment');
+ grid on
+
+ % PLANAR DISPLACEMENT %
+
+ figure
+ plot(Yf(:,2),Yf(:,1),'k-','LineWidth',2)
+ hold on
+ title('Planar Displacement');
+ xlabel('East [m]')
+ ylabel('North [m]');
+ h1=plot(Ya(end,2),Ya(end,1),'ko','MarkerSize',8,...
+ 'MarkerFaceColor','k');
+ h2=plot(Yd(end,2),Yd(end,1),'ks','MarkerSize',8,...
+ 'MarkerFaceColor','k');
+
+ legend([h1,h2],'Apogee/1st Drogue Deployment',...
+ '2nd Drogue Deployment',...
+ 'Location','southeast');
+ grid on
+
+
+
+
+ %Interpolation for less points --> visualization
+ Tinterp = linspace(Ta(1),Ta(end),settings.tSteps);
+ p = interp1(Ta,Ya(:,7),Tinterp);
+ q = interp1(Ta,Ya(:,8),Tinterp);
+ r = interp1(Ta,Ya(:,9),Tinterp);
+
+ % ANGULAR RATES %
+
+ figure;
+ suptitle('Angular rates on Time')
+ subplot(3,1,1);
+ plot(Tinterp,p,'k-')
+ hold on
+ xlabel('Time[s]');
+ ylabel('Roll rate p [rad/s]');
+ h1=plot(Ta(end),Ya(end,7),'ko','MarkerSize',8,...
+ 'MarkerFaceColor','k');
+ legend(h1,'Apogee/1st Drogue Deployment','Location','southeast');
+ grid on
+
+ subplot(3,1,2)
+ plot(Tinterp,q,'k-')
+ hold on
+ xlabel('Time[s]');
+ ylabel('Pitch rate q [rad/s]');
+ h1=plot(Ta(end),Ya(end,8),'ko','MarkerSize',8,...
+ 'MarkerFaceColor','k');
+ legend(h1,'Apogee/1st Drogue Deployment','Location','northeast');
+ grid on
+
+
+ subplot(3,1,3)
+ plot(Tinterp,r,'k-')
+ hold on
+ xlabel('Time[s]');
+ ylabel('Yaw rate r [rad/s]');
+ h1=plot(Ta(end),Ya(end,9),'ko','MarkerSize',8,...
+ 'MarkerFaceColor','k');
+ legend(h1,'Apogee/1st Drogue Deployment','Location','southeast');
+ grid on
+
+end
+
+%Resizing
+h = get(0,'children');
+scrsz = get(0,'ScreenSize');
+for i=1:length(h)
+ set(h(i),'OuterPosition',[0 0 scrsz(4) scrsz(4)])
+ %saveas(h(i), ['figure' num2str(i)], 'fig');
+end
+
+
+
+end
\ No newline at end of file
diff --git a/src/stoch_run.m b/src/stoch_run.m
new file mode 100644
index 0000000000000000000000000000000000000000..f4dad8e2bbb6dd790a0a25cfa993fb22fd14641a
--- /dev/null
+++ b/src/stoch_run.m
@@ -0,0 +1,158 @@
+function [ LP ,Z ] = stoch_run( settings )
+%STD RUN - This function runs a stochastic simulation
+% OUTPUTS
+% LP: Landing Points
+% Z: Apogee Altitudes
+
+% Author: Ruben Di Battista
+% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu
+% email: ruben.dibattista@skywarder.eu
+% Website: http://www.skywarder.eu
+% April 2014; Last revision: 29.V.2014
+% License: 2-clause BSD
+
+%Starting Attitude
+Q0 = angle2quat(settings.PHI,settings.OMEGA,0*pi/180,'ZYX')';
+
+%Starting State
+X0 = [0 0 0]';
+V0 = [0 0 0]';
+W0 = [0 0 0]';
+X0a = [X0;V0;W0;Q0;settings.m0;settings.Ixxf;settings.Iyyf;settings.Izzf];
+
+
+%PreAllocation
+LP = zeros(settings.stoch.N,3);
+Z = zeros(settings.stoch.N,1);
+ApoTime = zeros(settings.stoch.N,1);
+
+
+parfor_progress(settings.stoch.N);
+for i=1:settings.stoch.N
+ %% Wind Generation
+ [uw,vw,ww] = windgen(settings.wind.AzMin,settings.wind.AzMax,...
+ settings.wind.ElMin,settings.wind.ElMax,settings.wind.MagMin,...
+ settings.wind.MagMax);
+
+ %% ASCEND %%
+
+ [Ta,Ya] = ode45(@ascend,settings.ode.timeasc,X0a,settings.ode.optionsasc,...
+ settings,uw,vw,ww);
+
+ %% DROGUE 1 %%
+ para = 1; %Flag for Drogue 1
+
+ %Initial Condition are the last from ascend (need to rotate because
+ %velocities are outputted in body axes)
+ X0d1 = [Ya(end,1:3) quatrotate(quatconj(Ya(end,10:13)),Ya(end,4:6))];
+ [Td1,Yd1] = ode45(@descent_parachute,settings.ode.timedrg1,X0d1,...
+ settings.ode.optionsdrg1,settings,uw,vw,ww,para);
+
+ %% DROGUE 2 %%
+ para = 2; %Flag for Drogue 2
+
+ %Initial Condition are the last from drogue 1 descent
+ X0d2 = Yd1(end,:);
+ [Td2,Yd2] = ode45(@descent_parachute,settings.ode.timedrg2,X0d2,...
+ settings.ode.optionsdrg2,settings,uw,vw,ww,para);
+
+ %% MAIN %%
+ para = 3; %Flag for Main (Rogall)
+
+ %Initial Condition are the last from drogue 2 descent
+ X0m = Yd2(end,:);
+ [Tm,Ym] = ode45(@descent_parachute,settings.ode.timemain,X0m,...
+ settings.ode.optionsmain,settings,uw,vw,ww,para);
+
+ %% FINAL STATE ASSEMBLING %%
+
+ %Total State
+ Yf = [Ya(:,1:3) quatrotate(quatconj(Ya(:,10:13)),Ya(:,4:6));Yd1; Yd2;Ym];
+ %Total Time
+ %Tf = [Ta; Ta(end)+Td;Ta(end)+Td(end)+Tm];
+
+ LP(i,:) = Yf(end,1:3);
+ Z(i) = -Ya(end,3);
+ ApoTime(i) = Ta(end);
+
+ parfor_progress;
+
+end
+
+%Checking bad simulations
+if numel(LP(LP(:,3)<-10,:))
+ fprintf(['Some landing points might be incorrect' ...
+ 'please check parameters!\n']);
+end
+
+% Writing Things
+
+%Mean Landing Point
+xm = mean(LP(:,1));
+ym = mean(LP(:,2));
+
+%Mean Apogee Time
+ApoTimem = mean(ApoTime);
+
+%Std. Deviation Apogee Time
+ApoTimestd = std(ApoTime);
+
+%Mean Altitude
+Zm = mean(Z);
+
+%Std. Deviation Altitude
+Zstd = std(Z);
+
+%Best Fit Ellipse
+plot(LP(:,1),LP(:,2),'k*');
+
+% Printing to screen
+text =['Mean Landing Point:X:%3.3f m, Y:%3.3f m\n',...
+ 'Mean Altitude: %3.3f m || STD: %3.3f m\n',...
+ 'Mean Apogee Time: %3.3f s || STD: %3.3f s\n'];
+fprintf(text,xm,ym,Zm,Zstd,ApoTimem,ApoTimestd);
+
+if settings.plot == 1
+ %% PLOTTING THINGS
+
+
+ plot(xm,ym,'ks','MarkerSize',20);
+ hold on
+
+ %Point of launch
+ plot(0,0,'ro','MarkerSize',20);
+
+ %All the landing points
+ plot(LP(:,1),LP(:,2),'k+');
+
+ title('Landing Points');
+ xlabel('X [m]');
+ ylabel('Y [m]');
+ legend('Mean Landing Point','Launch Site','Landing Points');
+
+
+ %Histogram
+ [f,x] = hist(Z,10);
+ figure;
+ bar(x,f/settings.stoch.N);
+ title('Apogee Altitudes Distribution')
+ xlabel('Apogee [m]')
+ ylabel('n_i/n_{tot}')
+
+
+
+
+end
+
+%Resizing
+h = get(0,'children');
+scrsz = get(0,'ScreenSize');
+for i=1:length(h)
+ set(h(i),'OuterPosition',[0 0 scrsz(4) scrsz(4)])
+ %saveas(h(i), ['figure' num2str(i)], 'fig');
+end
+
+
+
+
+end
diff --git a/src/stoch_run_p.m b/src/stoch_run_p.m
new file mode 100644
index 0000000000000000000000000000000000000000..e32d2c7d141580108a93488a0caa910d644ccba1
--- /dev/null
+++ b/src/stoch_run_p.m
@@ -0,0 +1,156 @@
+function [ LP ,Z ] = stoch_run_p( settings )
+%STD RUN - This function runs a stochastic simulation (parallel)
+% OUTPUTS
+% LP: Landing Points
+% Z: Apogee Altitudes
+
+% Author: Ruben Di Battista
+% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu
+% email: ruben.dibattista@skywarder.eu
+% Website: http://www.skywarder.eu
+% April 2014; Last revision: 29.V.2014
+% License: 2-clause BSD
+
+%Starting Attitude
+Q0 = angle2quat(settings.PHI,settings.OMEGA,0*pi/180,'ZYX')';
+
+%Starting State
+X0 = [0 0 0]';
+V0 = [0 0 0]';
+W0 = [0 0 0]';
+X0a = [X0;V0;W0;Q0;settings.m0;settings.Ixxf;settings.Iyyf;settings.Izzf];
+
+
+%PreAllocation
+LP = zeros(settings.stoch.N,3);
+Z = zeros(settings.stoch.N,1);
+ApoTime = zeros(settings.stoch.N,1);
+
+parfor_progress(settings.stoch.N);
+parpool;
+parfor i=1:settings.stoch.N
+ %% Wind Generation
+ [uw,vw,ww] = windgen(settings.wind.AzMin,settings.wind.AzMax,...
+ settings.wind.ElMin,settings.wind.ElMax,settings.wind.MagMin,...
+ settings.wind.MagMax);
+
+ %% ASCEND %%
+
+ [Ta,Ya] = ode45(@ascend,settings.ode.timeasc,X0a,settings.ode.optionsasc,...
+ settings,uw,vw,ww);
+
+ %% DROGUE 1 %%
+ para = 1; %Flag for Drogue 1
+
+ %Initial Condition are the last from ascend (need to rotate because
+ %velocities are outputted in body axes)
+ X0d1 = [Ya(end,1:3) quatrotate(quatconj(Ya(end,10:13)),Ya(end,4:6))];
+ [Td1,Yd1] = ode45(@descent_parachute,settings.ode.timedrg1,X0d1,...
+ settings.ode.optionsdrg1,settings,uw,vw,ww,para);
+
+ %% DROGUE 2 %%
+ para = 2; %Flag for Drogue 2
+
+ %Initial Condition are the last from drogue 1 descent
+ X0d2 = Yd1(end,:);
+ [Td2,Yd2] = ode45(@descent_parachute,settings.ode.timedrg2,X0d2,...
+ settings.ode.optionsdrg2,settings,uw,vw,ww,para);
+
+ %% MAIN %%
+ para = 3; %Flag for Main (Rogall)
+
+ %Initial Condition are the last from drogue 2 descent
+ X0m = Yd2(end,:);
+ [Tm,Ym] = ode45(@descent_parachute,settings.ode.timemain,X0m,...
+ settings.ode.optionsmain,settings,uw,vw,ww,para);
+
+ %% FINAL STATE ASSEMBLING %%
+
+ %Total State
+ Yf = [Ya(:,1:3) quatrotate(quatconj(Ya(:,10:13)),Ya(:,4:6));Yd1; Yd2;Ym];
+ %Total Time
+ %Tf = [Ta; Ta(end)+Td;Ta(end)+Td(end)+Tm];
+
+ LP(i,:) = Yf(end,1:3);
+ Z(i) = -Ya(end,3);
+ ApoTime(i) = Ta(end);
+
+ parfor_progress;
+
+end
+
+%Checking bad simulations
+if numel(LP(LP(:,3)<-10,:))
+ fprintf(['Some landing points might be incorrect' ...
+ 'please check parameters!\n']);
+end
+
+% Writing Things
+
+%Mean Landing Point
+xm = mean(LP(:,1));
+ym = mean(LP(:,2));
+
+%Mean Apogee Time
+ApoTimem = mean(ApoTime);
+
+%Std. Deviation Apogee Time
+ApoTimestd = std(ApoTime);
+
+%Mean Altitude
+Zm = mean(Z);
+
+%Std. Deviation Altitude
+Zstd = std(Z);
+
+% Printing to screen
+text =['Mean Landing Point:X:%3.3f m, Y:%3.3f m\n',...
+ 'Mean Altitude: %3.3f m || STD: %3.3f m\n',...
+ 'Mean Apogee Time: %3.3f s || STD: %3.3f s\n'];
+fprintf(text,xm,ym,Zm,Zstd,ApoTimem,ApoTimestd);
+
+if settings.plot == 1
+ %% PLOTTING THINGS
+
+
+ plot(xm,ym,'bs','MarkerSize',20,'MarkerFacecolor','b');
+ hold on
+
+ %Point of launch
+ plot(0,0,'ro','MarkerSize',20,'MarkerFacecolor','r');
+
+ %All the landing points
+ plot(LP(:,1),LP(:,2),'k+');
+
+ title('Landing Points');
+ xlabel('X [m]');
+ ylabel('Y [m]');
+ legend('Mean Landing Point','Launch Site','Landing Points');
+
+
+ %Histogram
+ [f,x] = hist(Z,10);
+ figure;
+ bar(x,f/settings.stoch.N);
+ title('Apogee Altitudes Distribution')
+ xlabel('Apogee [m]')
+ ylabel('n_i/n_{tot}')
+
+
+
+
+end
+
+delete(gcp);
+
+%Resizing
+h = get(0,'children');
+scrsz = get(0,'ScreenSize');
+for i=1:length(h)
+ set(h(i),'OuterPosition',[0 0 scrsz(4) scrsz(4)])
+ %saveas(h(i), ['figure' num2str(i)], 'fig');
+end
+
+
+
+end
diff --git a/src/vert_windgen.m b/src/vert_windgen.m
new file mode 100644
index 0000000000000000000000000000000000000000..5765f4cf3d49007d4e2ec467c9700246f1b464c9
--- /dev/null
+++ b/src/vert_windgen.m
@@ -0,0 +1,23 @@
+function [ ww ] = vert_windgen( MagMin,MagMax )
+%windgen( MagMin,MagMax )
+%function that generates vertical wind component
+%
+
+
+
+% Author: Ruben Di Battista & Gabriele Poiana
+% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu
+% email: ruben.dibattista@skywarder.eu gabriele.poiana@skywarder.eu
+% Website: http://www.skywarder.eu
+% April 2014; Last revision: 17.I.2016
+% License: 2-clause BSD
+
+%Generating random value for magnitude
+
+x=rand;
+ww = MagMin + (MagMax-MagMin)*x;
+if x<0.5
+ ww=-ww;
+end
+end
+
diff --git a/src/wind_generator.m b/src/wind_generator.m
new file mode 100644
index 0000000000000000000000000000000000000000..bc8da33db1889c8abe83beaea9214213010bc69e
--- /dev/null
+++ b/src/wind_generator.m
@@ -0,0 +1,39 @@
+function [ uw,vw,ww ] = wind_generator( settings,z,t,Q )
+%wind_generator( settings,z,t,Q )
+%Function that generates wind components in body reference frame
+%Based on hwm07 model
+%nargin=4:rotation of wind components in body axis
+%nargin=3: wind components in NED
+%
+
+% Author: Gabriele Poiana
+% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu
+% email: gabriele.poiana@skywarder.eu
+% Website: http://www.skywarder.eu
+% January 2016; Last revision: 17.I.2016
+% License: 2-clause BSD
+
+
+h=-z+settings.z0;
+%global vento
+
+%% Horizontal Wind Components
+[uw,vw]=atmoshwm07(settings.wind.Lat,settings.wind.Long,h,'day',settings.wind.Day,...
+ 'seconds',settings.wind.Seconds+t,'model','quiet'); %NED reference
+
+%vento=[vento;uw,vw,ww,h];
+
+%% Wind Components
+if nargin==4
+ % Rotation in body reference frame (ascend)
+ wind=quatrotate(Q,[uw vw ww]); % wind speed in body reference
+
+ uw=wind(1);
+ vw=wind(2);
+ ww=wind(3);
+
+end
+
+end
+
+
diff --git a/src/windgen.m b/src/windgen.m
new file mode 100644
index 0000000000000000000000000000000000000000..ea003d70f05a659965cc6fbdff815cfdac1f25da
--- /dev/null
+++ b/src/windgen.m
@@ -0,0 +1,44 @@
+function [ uw,vw,ww ] = windgen( AzMin,AzMax,ElMin,ElMax,MagMin,MagMax )
+%windgen( AzMin,AzMax,ElMin,ElMax,MagMin,MagMax )
+%function that generates wind components in NED axes based on altitude
+%
+%Very basic implementation: it has to be improved (Power Spectral Density
+%etc..)
+%
+%
+%Vector Orientation
+%AzMin = 0; Minimum angle of Azimuth from North
+%AzMax = 2*pi; Maximum angle of Azimuth from North
+%ElMin = 0; Minimum angle of Elevation
+%ElMax = pi/2; Maximum angle of Elevatiom
+
+% Author: Ruben Di Battista
+% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu
+% email: ruben.dibattista@skywarder.eu
+% Website: http://www.skywarder.eu
+% April 2014; Last revision: 25.IV.2014
+% License: 2-clause BSD
+
+%Generating random values for orientation and magnitude
+Az = AzMin + (AzMax-AzMin)*rand;
+El = ElMin + (ElMax-ElMin)*rand;
+Mag = MagMin + (MagMax-MagMin)*rand;
+
+%Random Wind Vector
+
+%uw = Mag*cos(Az)*cos(El);
+%vw = Mag*cos(El)*sin(Az);
+%ww = -Mag*sin(El);
+
+R = Mag*angle2dcm(Az,El,0,'ZYX');
+uw = R(1,1);
+vw = R(1,2);
+ww = R(1,3);
+
+
+
+
+
+
+end
+
diff --git a/utils/Drag_curve.m b/utils/Drag_curve.m
new file mode 100644
index 0000000000000000000000000000000000000000..3895cd78f520687d28afe3d969274969f3ffe07c
--- /dev/null
+++ b/utils/Drag_curve.m
@@ -0,0 +1,48 @@
+% Author: Francesco Colombi
+% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu
+% email: francesco.colombi@skywarder.eu
+% Release date: 16/04/2016
+
+close all
+clear, clc
+
+[~, a, ~, rho] = atmoscoesa(0);
+
+% DATCOM DATA
+load for006_full.mat;
+CoeffsF = Coeffs;
+load for006_empty.mat;
+CoeffsE = Coeffs;
+
+givA = State.Alphas;
+givB = State.Betas;
+givH = State.Altitudes;
+givM = State.Machs;
+
+% GEOMETRIC DATA
+C = 0.174; %m Reference Length (Fuselage Diameter)
+S = 0.02378; %m2 Cross-sectional Surface
+L = 4.30; %m lunghezza missile
+
+% DRAG CURVE
+M = linspace(0, 3);
+V = a*M;
+CD = zeros(1, length(M));
+D = zeros(1, length(M));
+for n = 1:length(M)
+ CD(n) = interpn(givA,givM,givB,givH,CoeffsF.CD,0,M(n),0,0);
+ D(n) = 1/2 * rho * V(n)^2 * S * CD(n);
+end
+
+figure();
+subplot(2,1,1)
+plot(V, CD), grid on;
+ylabel('CD [-]')
+xlabel('V [m/s]')
+ylim([0, 0.7])
+
+% figure();
+subplot(2,1,2)
+plot(V, D), grid on;
+ylabel('D [N]')
+xlabel('V [m/s]')
\ No newline at end of file
diff --git a/utils/apogee_reached.m b/utils/apogee_reached.m
new file mode 100644
index 0000000000000000000000000000000000000000..dff278557ae8a3e640b6ca600dc42385ea8a00a4
--- /dev/null
+++ b/utils/apogee_reached.m
@@ -0,0 +1,175 @@
+% TARGET ??
+% usare il simulatore a 6dof per simulare diverse configuarzioni del
+% missile e vedere la quota massima raggiunta per poter fare un confronto e
+% decidere la soluzione pi� perormante
+
+%% COSA MI SERVE ??
+% 1) i vari file di DATCOM per le diverse configurazioni (geometrie e CG)
+% 2) le inerzie del missile per le diverse configurazioni
+% 3) il profilo della spinta (rocket engine data)
+
+%% COME OTTENERE LA QUOTA DELL'APOGEO DAL SIMULATORE ??
+
+% ODE-Function of the 6DOF Rigid Rocket Model
+% State = ( x y z | u v w | p q r | q0 q1 q2 q3 )
+%
+% (x y z): NED Earth's Surface Centered Frame ("Inertial") coordinates
+% (u v w): body frame velocities
+% (p q r): body frame angular rates
+% (q0 q1 q2 q3): attitude unit quaternion
+
+% % [T,Y]=MAIN();
+% %
+% % apogeo = max(-Y(:,3));
+% % fprintf('Apogee: %g [m] \n', apogeo);
+
+%% CODICE
+% close all
+clear, clc
+
+global bool
+bool = 0;
+% if bool = 0 the trends during the integration are NOT requested
+% if bool = 1 the trends during the integration are saved and plotted
+
+
+% se la geometria non cambia faccio caricare i dati con i coeff
+% aerodinamici da "config.m" solo una volta e poi cambio le inerzie dentro
+% al ciclo "for"
+run('config.m');
+
+% masse
+m_fuel = 18.6; %kg
+% m_empty = 40.4;
+% m_empty = [44.6:]; %kg
+% m_full = m_empty + m_fuel;
+m_full = linspace(44.6, 84.6, 30)';
+m_empty = m_full - m_fuel;
+
+
+% dimensioni missile
+D = 0.174; %m
+% D = 0.214; %m
+R = D/2; %m
+S = pi*D^2/4; %m2
+L = 4.30; %m
+
+% inerzie
+Ixx_full = m_full.*R^2/2;
+Ixx_empty = m_empty.*R^2/2;
+Iyy_full = m_full.*(R^2/4 + L^2/3);
+Iyy_empty = m_empty.*(R^2/4 + L^2/3);
+Izz_full = Iyy_full;
+Izz_empty = Iyy_empty;
+
+n_casi = length(Ixx_full);
+
+apogeo = zeros(n_casi, 1);
+max_dist = zeros(n_casi, 1);
+max_v = zeros(n_casi, 1);
+max_a = zeros(n_casi, 1);
+
+for n = 1:n_casi
+ % carica i dati dai file di DATCOM per ogni caso e salvali nel file
+ % .mat per il run del simulatore
+
+% % se invece per ogni caso cambia anche la geometria
+% file_name = ['for006_empty(',num2str(n),').mat'];
+% load(file_name);
+% save('for006_empty.mat', 'Coeffs', 'State');
+%
+% file_name = ['for006_full(',num2str(n),').mat'];
+% save('for006_full.mat', 'Coeffs', 'State');
+%
+% run('config.m');
+
+ % setto la massa e le inerzie del caso n-th
+ settings.ms = m_empty(n);
+ settings.m0 = m_full(n);
+
+ settings.Ixxf = Ixx_full(n); %Inertia to x-axis (Full)
+ settings.Ixxe = Ixx_empty(n); %Inertia to x-axis (Empty)
+ settings.Iyyf = Iyy_full(n); %Inertia to y-axis (Full)
+ settings.Iyye = Iyy_empty(n); %Inertia to y-axis (Empty)
+ settings.Izzf = Izz_full(n); %Inertia to z-axis (Full)
+ settings.Izze = Izz_empty(n); %Inertia to z-axis (Empty)
+
+% % prova senza inerzie
+% settings.Ixxf=2.1; %Inertia to x-axis (Full)
+% settings.Ixxe=1.4; %Inertia to x-axis (Empty)
+% settings.Iyyf=20; %Inertia to y-axis (Full)
+% settings.Iyye=15; %Inertia to y-axis (Empty)
+% settings.Izzf=20; %Inertia to z-axis (Full)
+% settings.Izze=15; %Inertia to z-axis (Empty)
+
+% save('settings.mat', 'settings')
+
+ % faccio girare il simulatore
+ [T,Y, Ta,Ya] = MAIN(settings);
+
+ [apogeo(n), iapogeo] = max(-Y(:,3));
+ fprintf('apogee is equal to %g m \n', apogeo(n));
+ fprintf('apogeo raggiunto in %g sec \n\n', Ta(end))
+
+% % plot vari
+% % posizione
+% x = Y(:,1);
+% y = Y(:,2);
+% z = -Y(:,3);
+% X = [x, y, z];
+%
+% % velocit�
+% N = length(x);
+% % % derivata centrale
+% % vx = (x(3:N)-x(1:N-2))./(T(3:N)-T(1:N-2));
+% % vy = (y(3:N)-y(1:N-2))./(T(3:N)-T(1:N-2));
+% % vz = (z(3:N)-z(1:N-2))./(T(3:N)-T(1:N-2));
+% % % aggiungo derivata estremi
+% % vx = [0; vx; (x(end)-x(end-1))/(T(end)-T(end-1))];
+% % vy = [0; vy; (y(end)-y(end-1))/(T(end)-T(end-1))];
+% % vz = [0; vz; (z(end)-z(end-1))/(T(end)-T(end-1))];
+% vx = Y(:,4);
+% vy = Y(:,5);
+% vz = -Y(:,6);
+% V = [vx, vy, vz];
+%
+% % accelerazione
+% % derivata centrale
+% ax = (vx(3:N)-vx(1:N-2))./(T(3:N)-T(1:N-2));
+% ay = (vy(3:N)-vy(1:N-2))./(T(3:N)-T(1:N-2));
+% az = (vz(3:N)-vz(1:N-2))./(T(3:N)-T(1:N-2));
+% % aggiungo derivata estremi
+% ax = [vx(2)/T(2); ax; (vx(end)-vx(end-1))/(T(end)-T(end-1))];
+% ay = [vy(2)/T(2); ay; (vy(end)-vy(end-1))/(T(end)-T(end-1))];
+% az = [vz(2)/T(2); az; (vz(end)-vz(end-1))/(T(end)-T(end-1))];
+% A = [ax, ay, az];
+%
+% mod_X = zeros(length(X), 1);
+% mod_V = mod_X;
+% mod_A = mod_X;
+%
+% for k = 1:length(X)
+% mod_X(k) = norm(X(k,:));
+% mod_V(k) = norm(V(k,:));
+% mod_A(k) = norm(A(k,:));
+% end
+% [max_dist, imax_dist] = max(mod_X);
+% [max_v, imax_v] = max(mod_V);
+%
+% % plot
+% figure(), hold on, grid on;
+% plot(T, z);
+% plot(T, vz);
+% plot(T, mod_V);
+% legend('z', 'vz', 'mod V');
+end
+
+% istogramma apogeo vs massa struttura
+figure(), grid on;
+stem(m_full, apogeo);
+xlabel('overall mass = structural mass + fuel(=18.6kg) [kg]');
+ylabel('apogee [m]')
+
+[maxApogeo, imaxApogeo] = max(apogeo);
+fprintf('massimo apogeo raggiunto: %g m \n', maxApogeo);
+fprintf('ottenuto con massa strutturale di %g kg \n\n', m_empty(imaxApogeo));
\ No newline at end of file