diff --git a/AUTHORS.md b/AUTHORS.md deleted file mode 100644 index 505a24659e59ab83bc0a6fdf5bd227b724dc73e1..0000000000000000000000000000000000000000 --- a/AUTHORS.md +++ /dev/null @@ -1,14 +0,0 @@ -# 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 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/run.txt b/LICENSE/run.txt deleted file mode 100644 index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0000000000000000000000000000000000000000 diff --git a/README.md b/README.md deleted file mode 100644 index c6b798cea4a8893c7cc7ceb89fbe46a48f9381ae..0000000000000000000000000000000000000000 --- a/README.md +++ /dev/null @@ -1,12 +0,0 @@ -# skyward-matlab-rocket-simulator -This is a script developed in MATLAB for the simulation of rocket dynamics. - -## How to use it -Clone the repository. Once in MATLAB, right click on the folder containig all the repository. -Click on "Add to Path". Click on "Selected Folder and Subfolders". - -## Info -For more info, read the README file contained in each folder. - -## License -See "LICENSE/LICENSE" file. diff --git a/Resources/README.md b/Resources/README.md deleted file mode 100644 index 18aa7da3033df2f8243aad3ab74cddae86ab2d6c..0000000000000000000000000000000000000000 --- a/Resources/README.md +++ /dev/null @@ -1,16 +0,0 @@ -# Resources -This folder contains all the input data needed to run the simulator - -## missileDATCOM data -This folder must cointain missileDATCOM data for empty and full states. -Saved as "for006_empty.mat" and "for006_full.mat". -In order to obtain those data from missileDATCOM output file you need to use the parser. - -## "config.m" -The config file is where you set up all the relevant paramenters about the simulation, such as geometry, engine and wind options. - -**STOCHASTIC SIMULATION** -The stochastic simulation only returns a plot of landing points and a histogram of the apogee altitudes distribution. - -## "start_simulation.m" -Run this file to start the simulation. \ No newline at end of file diff --git a/Resources/config.m b/Resources/config.m deleted file mode 100644 index 5a86a86eca5a0e7a96aad7d22c4562a276360172..0000000000000000000000000000000000000000 --- a/Resources/config.m +++ /dev/null @@ -1,227 +0,0 @@ -%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 % - -% 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 Type -settings.hwm = false; %Set to true to use HWM wind model -%otherwise the wind model will be the one described in "windgen.m" - -% Wind is randomly generated but constant during the simulation. -% 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; % Set to true to use parallel simulation - -% 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/Resources/for006_empty.mat b/Resources/for006_empty.mat deleted file mode 100644 index e69a45c67daa289ef496f68d9ea573bc31902058..0000000000000000000000000000000000000000 --- a/Resources/for006_empty.mat +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:68b67e58eedc8b0cb219a180227e7aace61760de1eb11340c5627cfea43a765c -size 5930920 diff --git a/Resources/for006_full.mat b/Resources/for006_full.mat deleted file mode 100644 index bd35de97093b0d8a30cc49804063e6014d317843..0000000000000000000000000000000000000000 --- a/Resources/for006_full.mat +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:65fc1a2656a94a7791ae10774e46666103c082949a846184559a74b44f75d83b -size 5930920 diff --git a/Resources/parfor_progress.txt b/Resources/parfor_progress.txt deleted file mode 100644 index afe234be093e7ad8424af41b8a2162a42d19b3e9..0000000000000000000000000000000000000000 --- a/Resources/parfor_progress.txt +++ /dev/null @@ -1,3 +0,0 @@ -2 -1 -1 diff --git a/Resources/start_simulation.m b/Resources/start_simulation.m deleted file mode 100644 index 0476e213318ebd2c0b6979ed29101c65428a429b..0000000000000000000000000000000000000000 --- a/Resources/start_simulation.m +++ /dev/null @@ -1,350 +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 - -%% caricamento dati -run('config.m'); - -global bool contatore; -global t_plot T_plot alpha_plot beta_plot M_plot CA_plot Drag_plot; -global 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 = 1; -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 - if settings.stoch.N > 1 - fprintf('Stochastic Simulation Fired...\n\n'); - if settings.stoch.parallel == 1 - [LP,Z]=stoch_run_p(settings); - else - [LP,Z]=stoch_run(settings); - end - else -[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('Tamb', 'Ttot', 'location', 'best') - -end - end \ No newline at end of file diff --git a/Sources/Ascent/apogee.m b/Sources/Ascent/apogee.m deleted file mode 100644 index 9b7830561cc9861617d6344c11b1afcf5d027669..0000000000000000000000000000000000000000 --- a/Sources/Ascent/apogee.m +++ /dev/null @@ -1,28 +0,0 @@ -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/Sources/Ascent/apogee_reached.m b/Sources/Ascent/apogee_reached.m deleted file mode 100644 index dff278557ae8a3e640b6ca600dc42385ea8a00a4..0000000000000000000000000000000000000000 --- a/Sources/Ascent/apogee_reached.m +++ /dev/null @@ -1,175 +0,0 @@ -% 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 diff --git a/Sources/Ascent/ascend.m b/Sources/Ascent/ascend.m deleted file mode 100644 index 9d2ebc8444396288626b0d9e3489830f8894f71e..0000000000000000000000000000000000000000 --- a/Sources/Ascent/ascend.m +++ /dev/null @@ -1,456 +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); - -if settings.hwm -% Wind Model -[uw vw ww] = wind_generator(settings,z,t,Q); -wind = [uw vw ww]; -else -% constant wind -wind = quatrotate(Q, [uw vw ww]); -end - -% % Wind model -% h = 0; -% if -z > 0 -% h = -z+settings.z0; -% end -% % -% % global alt -% % if bool == 1 -% % alt = [alt; h]; -% % end -% % -% % -% % global WIND -% % if bool == 1 -% % WIND = [WIND; uw, vw, ww]; -% % end - - - - -% 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/Sources/Ascent/ascend_constThrust.m b/Sources/Ascent/ascend_constThrust.m deleted file mode 100644 index 306757b24155de038dfff0f3967e59372d5984d3..0000000000000000000000000000000000000000 --- a/Sources/Ascent/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/Sources/Common/MAIN.m b/Sources/Common/MAIN.m deleted file mode 100644 index a6d000eb61f49c9477c4d91a705ad060bc4215cf..0000000000000000000000000000000000000000 --- a/Sources/Common/MAIN.m +++ /dev/null @@ -1,38 +0,0 @@ -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 - fprintf('Standard Simulation Fired...\n\n'); - [T,Y, Ta,Ya]=std_run(settings); - end - - - toc -end diff --git a/Sources/Common/interp4_easy.m b/Sources/Common/interp4_easy.m deleted file mode 100644 index ba3d18455f0159cfa8194f43904b7c6ac57ec782..0000000000000000000000000000000000000000 --- a/Sources/Common/interp4_easy.m +++ /dev/null @@ -1,40 +0,0 @@ -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/Sources/Common/parfor_progress.m b/Sources/Common/parfor_progress.m deleted file mode 100644 index 2f053a3723b5ecf5db33cd9006040ad94c838da0..0000000000000000000000000000000000000000 --- a/Sources/Common/parfor_progress.m +++ /dev/null @@ -1,81 +0,0 @@ -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/ - -narginchk(0, 1) - -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/Sources/Common/parfor_progress.txt b/Sources/Common/parfor_progress.txt deleted file mode 100644 index afe234be093e7ad8424af41b8a2162a42d19b3e9..0000000000000000000000000000000000000000 --- a/Sources/Common/parfor_progress.txt +++ /dev/null @@ -1,3 +0,0 @@ -2 -1 -1 diff --git a/Sources/Common/print_license.m b/Sources/Common/print_license.m deleted file mode 100644 index 51c086d073d45824625d0c2c8edf9123c37d72ad..0000000000000000000000000000000000000000 --- a/Sources/Common/print_license.m +++ /dev/null @@ -1,33 +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/Sources/Common/std_run.m b/Sources/Common/std_run.m deleted file mode 100644 index d8ed1c56c8f33fe8e96398a838c77f13e4eb2606..0000000000000000000000000000000000000000 --- a/Sources/Common/std_run.m +++ /dev/null @@ -1,279 +0,0 @@ -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]'; -m_tot = settings.m0 + settings.para1.mass + settings.para2.mass + settings.para3.mass; %Overall Mass -X0a = [X0;V0;W0;Q0;m_tot;settings.Ixxf;settings.Iyyf;settings.Izzf]; - -%% Wind Generation -if settings.hwm - uw = 0; - vw = 0; - ww = 0; -else - [uw,vw,ww] = windgen(settings.wind.AzMin,settings.wind.AzMax,... - settings.wind.ElMin,settings.wind.ElMax,settings.wind.MagMin,... - settings.wind.MagMax); -end -%% 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]; - -%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/Sources/Common/std_run_ballistic.m b/Sources/Common/std_run_ballistic.m deleted file mode 100644 index a5ef5d938b20624855ed2343632e0273ac2f0bcd..0000000000000000000000000000000000000000 --- a/Sources/Common/std_run_ballistic.m +++ /dev/null @@ -1,230 +0,0 @@ -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]'; -m_tot = settings.m0 + settings.para1.mass + settings.para2.mass + settings.para3.mass; %Overall Mass -X0a = [X0;V0;W0;Q0;m_tot;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/Sources/Common/stoch_run.m b/Sources/Common/stoch_run.m deleted file mode 100644 index 6c28af5a68e84451e7ac7a75fdbd9b03e368d7f7..0000000000000000000000000000000000000000 --- a/Sources/Common/stoch_run.m +++ /dev/null @@ -1,161 +0,0 @@ -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]'; -m_tot = settings.m0 + settings.para1.mass + settings.para2.mass + settings.para3.mass; %Overall Mass -X0a = [X0;V0;W0;Q0;m_tot;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 -if settings.hwm - uw = 0; - vw = 0; - ww = 0; -else - [uw,vw,ww] = windgen(settings.wind.AzMin,settings.wind.AzMax,... - settings.wind.ElMin,settings.wind.ElMax,settings.wind.MagMin,... - settings.wind.MagMax); -end -%% 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]; - %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,'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/Sources/Common/stoch_run_p.m b/Sources/Common/stoch_run_p.m deleted file mode 100644 index ba2bbb09ed1cada1afb5c220f0fd2ccbc8794848..0000000000000000000000000000000000000000 --- a/Sources/Common/stoch_run_p.m +++ /dev/null @@ -1,162 +0,0 @@ -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]'; -m_tot = settings.m0 + settings.para1.mass + settings.para2.mass + settings.para3.mass; %Overall Mass -X0a = [X0;V0;W0;Q0;m_tot;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 -if settings.hwm - uw = 0; - vw = 0; - ww = 0; -else - [uw,vw,ww] = windgen(settings.wind.AzMin,settings.wind.AzMax,... - settings.wind.ElMin,settings.wind.ElMax,settings.wind.MagMin,... - settings.wind.MagMax); -end -%% 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]; - %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/Sources/Descent/ballistic_descent.m b/Sources/Descent/ballistic_descent.m deleted file mode 100644 index 7547c0a61de80de059b0aa7701df070c69d3216a..0000000000000000000000000000000000000000 --- a/Sources/Descent/ballistic_descent.m +++ /dev/null @@ -1,332 +0,0 @@ -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); - -if settings.hwm -% Wind Model -[uw vw ww] = wind_generator(settings,z,t,Q); -wind = [uw vw ww]; -else -% constant wind -wind = quatrotate(Q, [uw vw ww]); -end - -% % Wind model -% h = 0; -% if -z > 0 -% h = -z+settings.z0; -% end -% % -% % if trends.bool == 1 -% % trends.alt = [trends.alt; h]; -% % -% % -% % if trends.bool == 1 -% % trends.WIND = [trends.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/Sources/Descent/crash.m b/Sources/Descent/crash.m deleted file mode 100644 index 4a0c7584c9ec21d036ac898314637162ee806e72..0000000000000000000000000000000000000000 --- a/Sources/Descent/crash.m +++ /dev/null @@ -1,22 +0,0 @@ -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/Sources/Descent/descent_parachute.m b/Sources/Descent/descent_parachute.m deleted file mode 100644 index 1cfc78b5ea960f28d1cd71cd19970cf09ee49b15..0000000000000000000000000000000000000000 --- a/Sources/Descent/descent_parachute.m +++ /dev/null @@ -1,150 +0,0 @@ -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); - - -% Adding Wind (supposed to be added in NED axes); - -if settings.hwm -% Wind Model -[uw, vw, ww] = wind_generator(settings,z,t); -wind = [uw vw ww]; -else -% constant wind -wind = [uw vw ww]; -end - -% % Wind model -% h = 0; -% if -z > 0 -% h = -z+settings.z0; -% end -% % -% % global alt -% % if bool == 1 -% % alt = [alt; h]; -% % end -% % -% % -% % global WIND -% % if bool == 1 -% % WIND = [WIND; uw, vw, ww]; -% % end - -% 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 = settings.para1.mass; %parachute mass - case 2 - S = settings.para2.S; - CD = settings.para2.CD; - CL = settings.para2.CL; - pmass = settings.para2.mass; - case 3 - S = settings.para3.S; - CD = settings.para3.CD; - CL = settings.para3.CL; - pmass = settings.para3.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/Sources/Descent/drg2_opening.m b/Sources/Descent/drg2_opening.m deleted file mode 100644 index 9e7b31af1ec3afd0371a6bd3b5b0642bc2161ace..0000000000000000000000000000000000000000 --- a/Sources/Descent/drg2_opening.m +++ /dev/null @@ -1,15 +0,0 @@ -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/Sources/Descent/main_opening.m b/Sources/Descent/main_opening.m deleted file mode 100644 index 18fd536e62201512bbb5f3de4295140fa2ff733a..0000000000000000000000000000000000000000 --- a/Sources/Descent/main_opening.m +++ /dev/null @@ -1,16 +0,0 @@ -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/Sources/Old/MAIN_old.m b/Sources/Old/MAIN_old.m deleted file mode 100644 index 2d55ef68aa5ee33f18c613e9ddce51d5e6627f3b..0000000000000000000000000000000000000000 --- a/Sources/Old/MAIN_old.m +++ /dev/null @@ -1,36 +0,0 @@ -function [T,Y] = MAIN() - %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 - - %Printing License - print_license(); - tic - - %Retrieving Parameters - run('config.m'); - - - %Checking if stochastic or standard simulation needed - if settings.stoch.N > 1 - fprintf('Stochastic Simulation Fired...\n\n'); - if settings.stoch.parallel == 1 - [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 - - toc -end diff --git a/Sources/Old/config_old.m b/Sources/Old/config_old.m deleted file mode 100644 index e5c91ced70d408a51662adbb59eeae11c5f7a814..0000000000000000000000000000000000000000 --- a/Sources/Old/config_old.m +++ /dev/null @@ -1,137 +0,0 @@ -%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 - -% clear all -% close all -% clc - -% LAUNCHPAD % -settings.z0 = 1200; %Launchpad Altitude -settings.lrampa = 3.5; %LaunchPad Length - - -%STARTING ATTITUDE SETUP % -settings.OMEGA = 89*pi/180; %Elevation Angle -settings.PHI = 0*pi/180; %Azimuth Angle from North Direction - -% ENGINE DETAILS % -% settings.T_coeffs = fliplr([-332.799894219468, 2209.93106349784, -5265.48638128018, 4608.82884815720, -133.133036494531]); - %[0 438.918868 -185.329879 33.642643 -2.517265]; -%Thrust on time profile (IV Order Interpolation) -settings.T_coeffs = [8034.5, 0, 0, 0, 0]; - -settings.ms = 40.4; %kg Dry Mass (No Propellant) -settings.tb = 5.12; %sec Burning Time -settings.mp = 18.6; %kg Mass of Propellant -settings.mfr = settings.mp/settings.tb; %kg/sec Mass Flow Rate -settings.m0 = settings.ms+settings.mp; %kg Total Starting Mass - -% 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.250; %m lunghezza missile - -% MASS GEOMERTY DETAILS % -% x-axis: along the fuselage -% y-axis: right wing -% z-axis: downward - -% NB: le inerzie vengono inserite nel codice di "apogee_reached.m" -% HP: inerzie razzo = inerzie cilindro pieno circoscritto -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) - -% % prova senza inerzie -% settings.Ixxf=0.21; %Inertia to x-axis (Full) -% settings.Ixxe=0.14; %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) - -% 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 %%% -settings.para1.S = 0.52; %Surface -settings.para1.mass = 0.02; %Parachute Mass -settings.para1.CD = 0.9; %Parachute Drag Coefficient - - -%Altitude of Main Parachute Opening -settings.zmain = 250; - -%%% MAIN %%% -%The drogue parachute effects are neglected - -settings.para2.S = 11.65; %Surface -settings.para2.mass = 0.43; %Parachute Mass -settings.para2.CD = 0.59; %Parachute Mass - - -% INTEGRATION OPTIONS % -settings.ode.timeasc = 0:0.01:2000; %Time span for ascend -settings.ode.timedrg = 0:0.01:2000; %Time span for drogue -settings.ode.timemain = 0:0.01:2000; %Time span for main - -settings.ode.optionsasc = odeset('AbsTol',1E-3,'RelTol',1E-3,... - 'Events',@apogee,'InitialStep',1); %ODE options for ascend -settings.ode.optionsdrg = 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 - -% WIND DETAILS % -% Wind is randomly generated. Setting the same values for min and max will -% fix the parameters of the wind. - -settings.wind.MagMin = 10; %Minimum Magnitude -settings.wind.MagMax = 10; %Maximum Magnitude -settings.wind.ElMin = 0*pi/180; %Minimum Elevation -settings.wind.ElMax = 0*pi/180; %Maximum Elevation (Max == 90 Deg) -settings.wind.AzMin = 45*pi/180; %Minimum Azimuth -settings.wind.AzMax = 45*pi/180; %Maximum Azimuth - -% STOCHASTIC DETAILS % -%If N>1 the stochastic routine is fired (different standard plots) -settings.stoch.N = 1; %Number of iterations -settings.stoch.parallel = true; %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/Sources/Old/descent_parachute_old.m b/Sources/Old/descent_parachute_old.m deleted file mode 100644 index 9fdd27bd292e7ab02d2b7d9df2864880c00e9b59..0000000000000000000000000000000000000000 --- a/Sources/Old/descent_parachute_old.m +++ /dev/null @@ -1,86 +0,0 @@ -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 - -x = Y(1); -y = Y(2); -z = Y(3); -u = Y(4); -v = Y(5); -w = Y(6); - -%Adding wind - -ur = u - uw; -vr = v - vw; -wr = w - ww; - -V_norm = norm([ur vr wr]); - -% Constants -if para == 1 - S = settings.para1.S; %Parachute Surface - CD = settings.para1.CD; %Parachute CD - pmass = settings.para1.mass; -else - S = settings.para2.S; - CD = settings.para2.CD; - pmass = settings.para2.mass; -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; - -X=D*ur/V_norm; -Y=D*vr/V_norm; -Z=D*wr/V_norm -m*g; - -%Forces -F = [-X,-Y,-Z]'; - - - -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 diff --git a/Sources/Old/std_run_old.m b/Sources/Old/std_run_old.m deleted file mode 100644 index 6d7386533e3b7848db6dacfcfe801ef6bc488a67..0000000000000000000000000000000000000000 --- a/Sources/Old/std_run_old.m +++ /dev/null @@ -1,237 +0,0 @@ -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 - -%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 %% -para = 1; %Flag for Drogue - -%Initial Condition are the last from ascend (need to rotate because -%velocities are outputted in body axes) -X0d = [Ya(end,1:3) quatrotate(quatconj(Ya(end,10:13)),Ya(end,4:6))]; -[Td,Yd] = ode45(@descent_parachute,settings.ode.timedrg,X0d,... - settings.ode.optionsdrg,settings,uw,vw,ww,para); - -%% MAIN %% -para = 2; %Flag for Main - - -%Initial Condition are the last from drogue descent -X0m = Yd(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));Yd;Ym]; -%Total Time -Tf = [Ta; Ta(end)+Td;Ta(end)+Td(end)+Tm]; - -%Parachute State -Yp = [Yd;Ym]; -Tp = [Ta(end)+Td;Ta(end)+Td(end)+Tm]; - -if settings.plot == 1 - %% PLOTTING THINGS - 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)+Td(end),-Yd(end,3)+settings.z0,'ks','MarkerSize',8,... - 'MarkerFaceColor','k'); - legend([h1 h2],'Apogee/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)+Td(end),Yd(end,4),'ks','MarkerSize',8,... - 'MarkerFaceColor','k'); - legend([h1,h2],'Apogee/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)+Td(end),Yd(end,5),'ks','MarkerSize',8,... - 'MarkerFaceColor','k'); - legend([h1,h2],'Apogee/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)+Td(end),Yd(end,6),'ks','MarkerSize',8,... - 'MarkerFaceColor','k'); - legend([h1,h2],'Apogee/Drogue Deployment','Main Parachute Deployment',... - 'Location','southeast'); - grid on - - %COMPLETE TRAJECTORY% - figure; - h0=plot3(Yf(1,1),Yf(1,2),-Yf(1,3)+settings.z0,'k+','MarkerSize',10); - hold on - plot3(Yf(:,1),Yf(:,2),-Yf(:,3)+settings.z0,'k-','LineWidth',2) - title('Complete Trajectory'); - xlabel('X [m]') - ylabel('Y [m]'); - zlabel('Z [m]'); - h1=plot3(Ya(end,1),Ya(end,2),-Ya(end,3)+settings.z0,'ko','MarkerSize',8,... - 'MarkerFaceColor','k'); - h2=plot3(Yd(end,1),Yd(end,2),-Yd(end,3)+settings.z0,'ks','MarkerSize',8,... - 'MarkerFaceColor','k'); - legend([h0,h1,h2],'Launch Pad','Apogee/Drogue Deployment',... - 'Main Parachute Deployment'); - grid on - - % PLANAR DISPLACEMENT % - - figure - plot(Yf(:,1),Yf(:,2),'k-','LineWidth',2) - hold on - title('Planar Displacement'); - xlabel('X [m]') - ylabel('Y [m]'); - h1=plot(Ya(end,1),Ya(end,2),'ko','MarkerSize',8,... - 'MarkerFaceColor','k'); - h2=plot(Yd(end,1),Yd(end,2),'ks','MarkerSize',8,... - 'MarkerFaceColor','k'); - legend([h1,h2],'Apogee/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)+Td(end),-Yd(end,3)+settings.z0,'ks','MarkerSize',8,... - 'MarkerFaceColor','k'); - legend([h1 h2],'Apogee/Drogue Deployment','Main Parachute Deployment') - 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('p [rad/s]'); - h1=plot(Ta(end),Ya(end,7),'ko','MarkerSize',8,... - 'MarkerFaceColor','k'); - legend(h1,'Apogee/Drogue Deployment','Location','southeast'); - grid on - - subplot(3,1,2) - plot(Tinterp,q,'k-') - hold on - xlabel('Time[s]'); - ylabel('q [rad/s]'); - h1=plot(Ta(end),Ya(end,8),'ko','MarkerSize',8,... - 'MarkerFaceColor','k'); - legend(h1,'Apogee/Drogue Deployment','Location','northeast'); - grid on - - - subplot(3,1,3) - plot(Tinterp,r,'k-') - hold on - xlabel('Time[s]'); - ylabel('r [rad/s]'); - h1=plot(Ta(end),Ya(end,9),'ko','MarkerSize',8,... - 'MarkerFaceColor','k'); - legend(h1,'Apogee/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 - diff --git a/Sources/Old/stoch_run_old.m b/Sources/Old/stoch_run_old.m deleted file mode 100644 index 3b3be34a976544476a8256eb787d4b3a12e97a91..0000000000000000000000000000000000000000 --- a/Sources/Old/stoch_run_old.m +++ /dev/null @@ -1,151 +0,0 @@ -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 %% - para = 1; %Flag for Drogue - - %Initial Condition are the last from ascend (need to rotate because - %velocities are outputted in body axes) - X0d = [Ya(end,1:3) quatrotate(quatconj(Ya(end,10:13)),Ya(end,4:6))]; - [Td,Yd] = ode45(@descent_parachute,settings.ode.timedrg,X0d,... - settings.ode.optionsdrg,settings,uw,vw,ww,para); - - %% MAIN %% - para = 2; %Flag for Main - - - %Initial Condition are the last from drogue descent - X0m = Yd(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));Yd;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/Sources/README.md b/Sources/README.md deleted file mode 100644 index 59890d7006c645baf7b3cfd688764145892c8077..0000000000000000000000000000000000000000 --- a/Sources/README.md +++ /dev/null @@ -1,22 +0,0 @@ -# Sources -This folder contains the source code. - -## Ascent -Contains the scripts used for the ascent phase in 6 DoF. - -## Common -Contains the scripts used during any time of the simulation. - -## Descent -Contains the scripts used for the descet phase in 3 DoF. - -Beacause of the 3 DoF, the pendulum effect after the apogee is not considered. - -## Old -Contains old scripts that have been modified or replaced. - -## Utilities -Contains scripts useful to test and analyze particular situations (Devs). - -## Wind -Contains the scripts used for implementing the wind model. diff --git a/Sources/Utilities/CLN.m b/Sources/Utilities/CLN.m deleted file mode 100644 index c93a9ca61ec2bdbb03bfbc21cb49c9857c2921a6..0000000000000000000000000000000000000000 --- a/Sources/Utilities/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/Sources/Utilities/Drag_curve.m b/Sources/Utilities/Drag_curve.m deleted file mode 100644 index 3895cd78f520687d28afe3d969274969f3ffe07c..0000000000000000000000000000000000000000 --- a/Sources/Utilities/Drag_curve.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 - -[~, 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/Sources/Utilities/THp.m b/Sources/Utilities/THp.m deleted file mode 100644 index 616769783dd26bc2970e450ad138180c57190261..0000000000000000000000000000000000000000 --- a/Sources/Utilities/THp.m +++ /dev/null @@ -1,38 +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 - - -%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/Sources/Utilities/coeffs_data.m b/Sources/Utilities/coeffs_data.m deleted file mode 100644 index 889999251b2ffefa20ea7ee497e40a6557e8309e..0000000000000000000000000000000000000000 --- a/Sources/Utilities/coeffs_data.m +++ /dev/null @@ -1,23 +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 \ No newline at end of file diff --git a/Sources/Utilities/launch_pad.m b/Sources/Utilities/launch_pad.m deleted file mode 100644 index 8928c14a0e9dd70e50295752649964f515651394..0000000000000000000000000000000000000000 --- a/Sources/Utilities/launch_pad.m +++ /dev/null @@ -1,85 +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 - -% 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/Sources/Utilities/rotating_axis.m b/Sources/Utilities/rotating_axis.m deleted file mode 100644 index 5cae9808d8c02b93bbdcb73ac450cee39c4eb119..0000000000000000000000000000000000000000 --- a/Sources/Utilities/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/Sources/Utilities/test_windgen.m b/Sources/Utilities/test_windgen.m deleted file mode 100644 index 942c2b25b240b947e85fb82547862a488def5dd4..0000000000000000000000000000000000000000 --- a/Sources/Utilities/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/Sources/Utilities/vis_coeffs.m b/Sources/Utilities/vis_coeffs.m deleted file mode 100644 index a0961ac89563d7a7220ab6e4b3865103b839b633..0000000000000000000000000000000000000000 --- a/Sources/Utilities/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') \ No newline at end of file diff --git a/Sources/Wind/vert_windgen.m b/Sources/Wind/vert_windgen.m deleted file mode 100644 index 5765f4cf3d49007d4e2ec467c9700246f1b464c9..0000000000000000000000000000000000000000 --- a/Sources/Wind/vert_windgen.m +++ /dev/null @@ -1,23 +0,0 @@ -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/Sources/Wind/wind_generator.m b/Sources/Wind/wind_generator.m deleted file mode 100644 index 55d1bba7b264a7bde2dcb70f57936fe007b751fa..0000000000000000000000000000000000000000 --- a/Sources/Wind/wind_generator.m +++ /dev/null @@ -1,39 +0,0 @@ -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 hwm 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]=atmoshwm(settings.wind.Lat,settings.wind.Long,h,'day',settings.wind.Day,... - 'seconds',settings.wind.Seconds+t,'model','quiet'); %NED reference -ww=settings.wind.ww; -%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/Sources/Wind/windgen.m b/Sources/Wind/windgen.m deleted file mode 100644 index ea003d70f05a659965cc6fbdff815cfdac1f25da..0000000000000000000000000000000000000000 --- a/Sources/Wind/windgen.m +++ /dev/null @@ -1,44 +0,0 @@ -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 -