diff --git a/LICENSE/LICENSE b/LICENSE/LICENSE
deleted file mode 100644
index c9663538595f803b0ae7d52c4c959bbb9d683d65..0000000000000000000000000000000000000000
--- a/LICENSE/LICENSE
+++ /dev/null
@@ -1,11 +0,0 @@
-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/print_license.m b/LICENSE/print_license.m
deleted file mode 100644
index d331bb432e6381cdbf759be10f6206459aa89936..0000000000000000000000000000000000000000
--- a/LICENSE/print_license.m
+++ /dev/null
@@ -1,32 +0,0 @@
-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/LICENSE/run.txt b/LICENSE/run.txt
deleted file mode 100644
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0000000000000000000000000000000000000000
diff --git a/Live_DATCOM/datcoming_mac.bat b/Live_DATCOM/datcoming_mac.bat
deleted file mode 100755
index 32778c6c9e5a1f748b770b63ddbdf7f3fe3bbe0e..0000000000000000000000000000000000000000
--- a/Live_DATCOM/datcoming_mac.bat
+++ /dev/null
@@ -1,9 +0,0 @@
-@echo off
-
-pushd \\Mac\Home\Documents\file_skyward\skyward-matlab-rocket-simulator\Live_DATCOM
-
-start datcom.exe
-
-sleep 1
-
-popd
\ No newline at end of file
diff --git a/scratch/CLN.m b/scratch/CLN.m
deleted file mode 100644
index c93a9ca61ec2bdbb03bfbc21cb49c9857c2921a6..0000000000000000000000000000000000000000
--- a/scratch/CLN.m
+++ /dev/null
@@ -1,51 +0,0 @@
-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
deleted file mode 100644
index 0ea93869b61a07100c14365b3503dfea61e8030f..0000000000000000000000000000000000000000
--- a/scratch/THp.m
+++ /dev/null
@@ -1,32 +0,0 @@
-% 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
-
-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;
diff --git a/scratch/ascend_constThrust.m b/scratch/ascend_constThrust.m
deleted file mode 100644
index 306757b24155de038dfff0f3967e59372d5984d3..0000000000000000000000000000000000000000
--- a/scratch/ascend_constThrust.m
+++ /dev/null
@@ -1,454 +0,0 @@
-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
deleted file mode 100644
index 2a390873f7e68849d8c109cf3666efd6096747f3..0000000000000000000000000000000000000000
--- a/scratch/coeffs_data.m
+++ /dev/null
@@ -1,24 +0,0 @@
-% 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
diff --git a/scratch/rotating_axis.m b/scratch/rotating_axis.m
deleted file mode 100644
index 5cae9808d8c02b93bbdcb73ac450cee39c4eb119..0000000000000000000000000000000000000000
--- a/scratch/rotating_axis.m
+++ /dev/null
@@ -1,74 +0,0 @@
-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
deleted file mode 100644
index 942c2b25b240b947e85fb82547862a488def5dd4..0000000000000000000000000000000000000000
--- a/scratch/test_windgen.m
+++ /dev/null
@@ -1,7 +0,0 @@
-% 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
deleted file mode 100644
index 42f8d8a361cc30d9703dbb996b4065b62e4dcc5b..0000000000000000000000000000000000000000
--- a/scratch/vis_coeffs.m
+++ /dev/null
@@ -1,48 +0,0 @@
-% 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')
diff --git a/utils/config_R2A.m b/utils/old_configs/config_R2A.m
similarity index 100%
rename from utils/config_R2A.m
rename to utils/old_configs/config_R2A.m
diff --git a/utils/config_R2A_hermes.m b/utils/old_configs/config_R2A_hermes.m
similarity index 100%
rename from utils/config_R2A_hermes.m
rename to utils/old_configs/config_R2A_hermes.m
diff --git a/utils/config_R2A_hermes_V1.m b/utils/old_configs/config_R2A_hermes_V1.m
similarity index 100%
rename from utils/config_R2A_hermes_V1.m
rename to utils/old_configs/config_R2A_hermes_V1.m