From 50e556d4887f374406834e7b5a2bb2b3dcf41c1c Mon Sep 17 00:00:00 2001 From: Filippo Inno <filippo.inno@skywarder.eu> Date: Sat, 16 Dec 2017 12:40:26 +0100 Subject: [PATCH] global variables removed; figures and script names changed; some bugs fixed. --- LICENSE/print_license.m | 32 ++ src/ascend.m | 492 +++++++++++++----------------- src/ascend_plot.m | 37 +++ src/ascend_plot.mat | 3 + src/ballistic_descent.m | 322 ++++++------------- src/config.m | 302 +++++++++--------- src/descent_parachute.m | 158 ++++------ src/event_apogee.m | 28 ++ src/event_drg2_opening.m | 14 + src/event_landing.m | 22 ++ src/event_main_opening.m | 16 + src/interp4_easy.m | 18 +- src/parfor_progress.m | 2 +- src/parfor_progress.txt | 50 ++- src/standard_plot.m | 196 ++++++++++++ src/start_simulation.m | 338 +++++--------------- src/std_run.m | 141 +++++---- src/std_run_ballistic.m | 152 ++++----- src/stoch_run.m | 95 +++--- src/stoch_run_bal.m | 69 ++--- src/stoch_run_bal_p.m | 72 ++--- src/stoch_run_p.m | 89 +++--- src/wind_const_generator.m | 35 +++ src/wind_matlab_generator.m | 40 +++ src/wind_vert_generator.m | 20 ++ {src => src_old}/MAIN.m | 0 {src => src_old}/apogee.m | 0 src_old/ascend.m | 456 +++++++++++++++++++++++++++ src_old/ballistic_descent.m | 339 ++++++++++++++++++++ src_old/config.m | 262 ++++++++++++++++ {src => src_old}/crash.m | 0 src_old/descent_parachute.m | 158 ++++++++++ {src => src_old}/drg2_opening.m | 0 src_old/interp4_easy.m | 40 +++ {src => src_old}/launch_pad.m | 0 {src => src_old}/main_opening.m | 0 src_old/parfor_progress.m | 81 +++++ src_old/parfor_progress.txt | 3 + src_old/plot_data.m | 156 ++++++++++ {src => src_old}/print_license.m | 0 src_old/start_simulation.m | 356 +++++++++++++++++++++ src_old/std_run.m | 283 +++++++++++++++++ src_old/std_run_ballistic.m | 229 ++++++++++++++ src_old/stoch_run.m | 165 ++++++++++ src_old/stoch_run_bal.m | 142 +++++++++ src_old/stoch_run_bal_p.m | 139 +++++++++ src_old/stoch_run_p.m | 163 ++++++++++ {src => src_old}/vert_windgen.m | 0 {src => src_old}/wind_generator.m | 0 {src => src_old}/windgen.m | 0 utils/launch_pad.m | 85 ++++++ 51 files changed, 4428 insertions(+), 1372 deletions(-) create mode 100644 LICENSE/print_license.m create mode 100644 src/ascend_plot.m create mode 100644 src/ascend_plot.mat create mode 100644 src/event_apogee.m create mode 100644 src/event_drg2_opening.m create mode 100644 src/event_landing.m create mode 100644 src/event_main_opening.m create mode 100644 src/standard_plot.m create mode 100644 src/wind_const_generator.m create mode 100644 src/wind_matlab_generator.m create mode 100644 src/wind_vert_generator.m rename {src => src_old}/MAIN.m (100%) rename {src => src_old}/apogee.m (100%) create mode 100644 src_old/ascend.m create mode 100644 src_old/ballistic_descent.m create mode 100644 src_old/config.m rename {src => src_old}/crash.m (100%) create mode 100644 src_old/descent_parachute.m rename {src => src_old}/drg2_opening.m (100%) create mode 100644 src_old/interp4_easy.m rename {src => src_old}/launch_pad.m (100%) rename {src => src_old}/main_opening.m (100%) create mode 100644 src_old/parfor_progress.m create mode 100644 src_old/parfor_progress.txt create mode 100644 src_old/plot_data.m rename {src => src_old}/print_license.m (100%) create mode 100644 src_old/start_simulation.m create mode 100644 src_old/std_run.m create mode 100644 src_old/std_run_ballistic.m create mode 100644 src_old/stoch_run.m create mode 100644 src_old/stoch_run_bal.m create mode 100644 src_old/stoch_run_bal_p.m create mode 100644 src_old/stoch_run_p.m rename {src => src_old}/vert_windgen.m (100%) rename {src => src_old}/wind_generator.m (100%) rename {src => src_old}/windgen.m (100%) create mode 100644 utils/launch_pad.m diff --git a/LICENSE/print_license.m b/LICENSE/print_license.m new file mode 100644 index 00000000..d331bb43 --- /dev/null +++ b/LICENSE/print_license.m @@ -0,0 +1,32 @@ +function [ ] = print_license() +% Function that prints at first run the LICENSE saved in the LICENSE folder +% in the same directory of this file + +% Author: Ruben Di Battista +% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu +% email: ruben.dibattista@skywarder.eu +% Website: http://www.skywarder.eu +% April 2014; Last revision: 25.IV.2014 +% License: 2-clause BSD + +if exist('LICENSE','dir') + if not(exist('LICENSE/run.txt','file')) + fh = fopen('LICENSE/LICENSE','r'); + while not(feof(fh)) + l = fgetl(fh); + fprintf(strcat(l,'\n')); + end + fclose(fh); + fh = fopen('LICENSE/run.txt','wt'); + fprintf(fh,strcat('License accepted:',datestr(clock,0))); + fclose(fh); + + agree = input('Write ''AGREE'' to accept the terms\n\n','s'); + while not(strcmpi(agree,'agree')) + agree = input('Write ''AGREE'' to accept the terms\n\n','s'); + end + end +end + +end + diff --git a/src/ascend.m b/src/ascend.m index 662c2b82..87aca60b 100644 --- a/src/ascend.m +++ b/src/ascend.m @@ -1,4 +1,4 @@ -function [ dY ] = ascend( t,Y,settings,uw,vw,ww ) +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 ) % @@ -45,58 +45,28 @@ 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); +Q = [q0 q1 q2 q3]; +Q_conj = [q0 -q1 -q2 -q3]; +normQ = norm(Q); -if abs(normQ-1)>0.1 +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); +%% ADDING WIND (supposed to be added in NED axes); if settings.wind.model -% Wind Model -[uw,vw,ww] = wind_generator(settings,z,t,Q); -wind = [ uw,vw,ww ]; + [uw,vw,ww] = wind_matlab_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 -% -% % 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 = quatrotate(Q, [uw vw ww]); % constant wind + +end % Relative velocities (plus wind); ur = u - wind(1); @@ -105,29 +75,31 @@ 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 +%% CONSTANTS + +S = settings.S; % [m^2] cross surface +C = settings.C; % [m] caliber +CoeffsE = settings.CoeffsE; % Empty Rocket Coefficients +CoeffsF = settings.CoeffsF; % Full Rocket Coefficients +g = 9.80655; % [N/kg] module of gravitational field at zero +tb = settings.tb; % [s] Burning Time +mfr = settings.mfr; % [kg/s] Mass Flow Rate + +% inertias for full configuration (with all the propellant embarqued) obtained with CAD's +Ixxf = settings.Ixxf; % [kg*m^2] Inertia to x-axis +Iyyf = settings.Iyyf; % [kg*m^2] Inertia to y-axis +Izzf = settings.Izzf; % [kg*m^2] Inertia to z-axis -g = 9.81; -%Time-Dependant Properties -m0 = settings.m0; %kg mass -mfr = settings.mfr; -tb = settings.tb; +% inertias for empty configuration (all the propellant consumed) obtained with CAD's +Ixxe = settings.Ixxe; % [kg*m^2] Inertia to x-axis +Iyye = settings.Iyye; % [kg*m^2] Inertia to y-axis +Izze = settings.Izze; % [kg*m^2] Inertia to z-axis -Ixxf= settings.Ixxf; -Ixxe= settings.Ixxe; -Iyyf= settings.Iyyf; -Iyye= settings.Iyye; -Izzf= settings.Izzf; -Izze= settings.Izze; +%% TIME-DEPENDENTS VARIABLES dI = 1/tb*([Ixxf Iyyf Izzf]'-[Ixxe Iyye Izze]'); @@ -136,10 +108,9 @@ if t<tb 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 +else % for t >= tb the fligth condition is the empty one(no interpolation needed) mdot = 0; Ixxdot = 0; Iyydot = 0; @@ -147,294 +118,203 @@ else T=0; end -global T_plot -if bool == 1 - T_plot(contatore) = T; % thrust -end +T_value = T; -%Atmosphere -if -z<0 +%% ATMOSPHERE DATA + +if -z < 0 % z is directed as the gravity vector z = 0; end -[~, a, ~, rho] = atmoscoesa(-z+settings.z0); +[~, a, ~, rho] = atmosisa(-z+settings.z0); M = V_norm/a; +M_value = M; -global M_plot -if bool == 1 - M_plot(contatore) = M; % mach -end -%Aero Angles -if not(u<1e-1 || V_norm<1e-3); + +%% AERODYNAMICS ANGLES + +if not(u < 1e-1 || V_norm < 1e-3) alpha = atan(wr/ur); beta = asin(vr/V_norm); else - alpha =0; + alpha = 0; beta = 0; end -global alpha_plot beta_plot -if bool == 1 - alpha_plot(contatore) = alpha; - beta_plot(contatore) = beta; -end +alpha_value = alpha; +beta_value = beta; -% beta = 0; % prova simulazione con beta imposto uguale a 0 +%% DATCOM COEFFICIENTS -% Coeffs. Interpolation givA = settings.Alphas*pi/180; givB = settings.Betas*pi/180; givH = settings.Altitudes; givM = settings.Machs; -% Interpolation error check +%% INTERPOLATION AT THE BOUNDARIES 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); +if M < givM(1) + 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); + + 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); +elseif alpha < givA(1) + alpha = givA(1); - - end end if beta > givB(end) - beta= givB(end); -else if beta < givB(1) + + beta = givB(end); + +elseif beta < givB(1) + beta = givB(1); - end end -if -z >givH(end) +if -z > givH(end) + z = -givH(end); -else if -z < givH(1) + +elseif -z < givH(1) + z = -givH(1); - end + end -global XCP -% 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); -XCPf =interp4_easy(givA,givM,givB,givH,CoeffsF.X_C_P,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); -XCPe =interp4_easy(givA,givM,givB,givH,CoeffsE.X_C_P,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 + +%% CHOSING THE FULL CONDITION VALUE +% interpolation of the coefficients with the value in the nearest condition of the Coeffs matrix + +CAf = interp4_easy(givA,givM,givB,givH,CoeffsF.CA,alpha,M,beta,-z); +CYBf = interp4_easy(givA,givM,givB,givH,CoeffsF.CYB,alpha,M,beta,-z); +CNAf = interp4_easy(givA,givM,givB,givH,CoeffsF.CNA,alpha,M,beta,-z); +Clf = interp4_easy(givA,givM,givB,givH,CoeffsF.CLL,alpha,M,beta,-z); +Clpf = interp4_easy(givA,givM,givB,givH,CoeffsF.CLLP,alpha,M,beta,-z); +Cmaf = interp4_easy(givA,givM,givB,givH,CoeffsF.CMA,alpha,M,beta,-z); +Cmadf = interp4_easy(givA,givM,givB,givH,CoeffsF.CMAD,alpha,M,beta,-z); +Cmqf = interp4_easy(givA,givM,givB,givH,CoeffsF.CMQ,alpha,M,beta,-z); +Cnbf = interp4_easy(givA,givM,givB,givH,CoeffsF.CLNB,alpha,M,beta,-z); +Cnrf = interp4_easy(givA,givM,givB,givH,CoeffsF.CLNR,alpha,M,beta,-z); +Cnpf = interp4_easy(givA,givM,givB,givH,CoeffsF.CLNP,alpha,M,beta,-z); +CDf = interp4_easy(givA,givM,givB,givH,CoeffsF.CD,alpha,M,beta,-z); +XCPf = interp4_easy(givA,givM,givB,givH,CoeffsF.X_C_P,alpha,M,beta,-z); + +%% CHOSING THE EMPTY CONDITION VALUE +% interpolation of the coefficients with the value in the nearest condition of the Coeffs matrix + +CAe = interp4_easy(givA,givM,givB,givH,CoeffsE.CA,alpha,M,beta,-z); +CYBe = interp4_easy(givA,givM,givB,givH,CoeffsE.CYB,alpha,M,beta,-z); +CNAe = interp4_easy(givA,givM,givB,givH,CoeffsE.CNA,alpha,M,beta,-z); +Cle = interp4_easy(givA,givM,givB,givH,CoeffsE.CLL,alpha,M,beta,-z); +Clpe = interp4_easy(givA,givM,givB,givH,CoeffsE.CLLP,alpha,M,beta,-z); +Cmae = interp4_easy(givA,givM,givB,givH,CoeffsE.CMA,alpha,M,beta,-z); +Cmade = interp4_easy(givA,givM,givB,givH,CoeffsE.CMAD,alpha,M,beta,-z); +Cmqe = interp4_easy(givA,givM,givB,givH,CoeffsE.CMQ,alpha,M,beta,-z); +Cnbe = interp4_easy(givA,givM,givB,givH,CoeffsE.CLNB,alpha,M,beta,-z); +Cnre = interp4_easy(givA,givM,givB,givH,CoeffsE.CLNR,alpha,M,beta,-z); +Cnpe = interp4_easy(givA,givM,givB,givH,CoeffsE.CLNP,alpha,M,beta,-z); +CDe = interp4_easy(givA,givM,givB,givH,CoeffsE.CD,alpha,M,beta,-z); +XCPe = interp4_easy(givA,givM,givB,givH,CoeffsE.X_C_P,alpha,M,beta,-z); + +%% LINEAR INTERPOLATION BETWEEN THE TWO CONDITIONS +% Computing the value of the aerodynamics coefficients at a certain time +% Needed only for t<tb because for t>=tb the condition is the empty one + +if t < tb CA = t/tb*(CAe-CAf)+CAf; - CYB= t/tb*(CYBe-CYBf)+CYBf; - CNA= t/tb*(CNAe-CNAf)+CNAf; + 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; - XCP(contatore)= t/tb*(XCPe-XCPf)+XCPf; + 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; +CD = t/tb*(CDe-CDf)+CDf; +XCP_value = t/tb*(XCPe-XCPf)+XCPf; else CA = CAe; - CYB= CYBe; - CNA= CNAe; + 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; - XCP(contatore) = XCPe; -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; + Clp = Clpe; + Cma = Cmae; + Cmad =Cmade; + Cmq = Cmqe; + Cnb = Cnbe; + Cnr = Cnre; + Cnp = Cnpe; +CD = CDe; +XCP_value = XCPe; 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 +%% FORCES +% first computed in the body-frame reference system +qdyn = 0.5*rho*V_norm^2; %[Pa] dynamics pressure +qdynL_V = 0.5*rho*V_norm*S*C; % -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; +X = qdyn*S*CA; %[N] x-body component of the aerodynamics force +Y = qdyn*S*CYB*beta; %[N] y-body component of the aerodynamics force +Z = qdyn*S*CNA*alpha; %[N] z-body component of the aerodynamics force +Fg = quatrotate(Q,[0 0 m*g])'; %[N] force due to the gravity +F = Fg +[-X+T,+Y,-Z]'; %[N] total forces vector +%% STATE DERIVATIVES +% velocity +du = F(1)/m-q*w+r*v; +dv = F(2)/m-r*u+p*w; +dw = F(3)/m-p*v+q*u; % 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)... +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)... +dr = (Ixx-Iyy)/Izz*p*q + qdynL_V/Izz*(V_norm*Cnb*beta + (Cnr*r+Cnp*p)*C/2)... -Izzdot*r/Izz; +Xb = quatrotate(Q,[x y z]); % Launch Pad-relative coordinates - -% 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 +if Xb(1) < settings.lrampa % No torque on the Launch dp = 0; dq = 0; dr = 0; - - %Ground Reaction. - if T < m*g - du = 0; - dv = 0; - dw = 0; - - end + + if T < m*g % No velocity untill T = Fg + du = 0; + dv = 0; + dw = 0; + + end end +% Quaternion +OM = 1/2* [ 0 -p -q -r ; + p 0 r -q ; + q -r 0 p ; + r q -p 0 ]; +dQQ = OM*Q'; - - - - -OM = 1/2*[ - 0 -p -q -r - p 0 r -q - q -r 0 p - r q -p 0]; - - -dQQ = OM*Q'; +%% FINAL DERIVATIVE STATE ASSEMBLING dY(1:3) = Vels; dY(4) = du; @@ -448,10 +328,56 @@ dY(14) = mdot; dY(15) = Ixxdot; dY(16) = Iyydot; dY(17) = Izzdot; -dY=dY'; +dY = dY'; + + +%% PERSISTENT VARIABLES -if bool == 1 +persistent XCP contatore t_plot ... + T_plot M_plot CA_plot alpha_plot beta_plot Drag_plot Forces_plot + + +%% SAVING THE QUANTITIES FOR THE PLOTS + +if settings.ascend_plot + + if isempty (contatore) + contatore = 1; + t_plot(contatore) = 0; + T_plot(contatore) = 0; + CA_plot(contatore) = 0; + Drag_plot(contatore) = 0; + Forces_plot(contatore) = 0; + M_plot(contatore) = 0; + XCP(contatore) = 0; + alpha_plot(contatore) = 0; + beta_plot(contatore) = 0; + end + t_plot(contatore) = t; + beta_plot(contatore) = beta_value; + alpha_plot(contatore) = alpha_value; + XCP(contatore) = XCP_value; + M_plot(contatore) = M_value; + T_plot(contatore) = T_value; + CA_plot(contatore) = CD; + Drag_plot(contatore) = qdyn*S*CD; + Forces_plot(contatore) = F(1); contatore = contatore + 1; -end + if Vels(3) <= 0 && t > tb + + ascend.XCP = XCP; + ascend.t = t_plot; + ascend.T = T_plot; + ascend.M = M_plot; + ascend.CA = CA_plot; + ascend.alpha = alpha_plot; + ascend.beta = beta_plot; + ascend.Drag = Drag_plot; + ascend.Forces = Forces_plot; + + + + save ('ascend_plot.mat', 'ascend') + end end diff --git a/src/ascend_plot.m b/src/ascend_plot.m new file mode 100644 index 00000000..6d45a554 --- /dev/null +++ b/src/ascend_plot.m @@ -0,0 +1,37 @@ + +load('ascend_plot.mat') + +%% PLOTS + +figure(); +plot(ascend.t, ascend.T, '.'), title('Thrust vs time'), grid on; +xlabel('Time [s]'); ylabel('Thrust [N]'); + +figure(); +plot(ascend.t, ascend.alpha*180/pi), title('alpha vs time'), grid on; +xlabel('Time [s]'); ylabel('alpha [deg]') + +figure(); +plot(ascend.t, ascend.beta*180/pi), title('beta vs time'), grid on; +xlabel('Time [s]'); ylabel('beta [deg]') + +figure(); +plot(ascend.t(1:end-1), ascend.M(1:end-1)), title('mach vs time'), grid on; +xlabel('Time [s]'); ylabel('Mach M [/]') + +figure(); +plot(ascend.t, ascend.CA), title('Aerodyn Coeff vs time'), grid on; +xlabel('Time [s]'); ylabel('Drag Coeff CD [/]') + +figure(); +plot(ascend.t, ascend.Drag), title('Drag vs time'), grid on; +xlabel('Time [s]'); ylabel('Drag D [N]') + +figure(); +plot(ascend.t, ascend.Forces), title('Axial Force vs time'), grid on; +xlabel('Time [s]'); ylabel('Axial force [N]') + +figure(); +plot(ascend.t, ascend.XCP,'.'), title('Stability margin vs time'), grid on; +xlabel('Time [s]'); ylabel('S.M.[/]') + diff --git a/src/ascend_plot.mat b/src/ascend_plot.mat new file mode 100644 index 00000000..29712b8a --- /dev/null +++ b/src/ascend_plot.mat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:f03f4c28b0bcef1c632b6b5055beec0ae50c71a54b995a3ea1c116a81939a09a +size 938421 diff --git a/src/ballistic_descent.m b/src/ballistic_descent.m index a1fbb657..b5840731 100644 --- a/src/ballistic_descent.m +++ b/src/ballistic_descent.m @@ -1,4 +1,4 @@ -function [ dY ] = ballistic_descent( t,Y,settings,uw,vw,ww ) +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 ) % @@ -27,67 +27,44 @@ function [ dY ] = ballistic_descent( t,Y,settings,uw,vw,ww ) % 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); +% 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); +Q_conj = [ q0 -q1 -q2 -q3]; +normQ = norm(Q); -if abs(normQ-1)>0.1 +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); -% 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 ]; +if settings.wind.model + [uw,vw,ww] = wind_matlab_generator(settings,z,t,Q); + wind = [uw,vw,ww]; + +else + +wind = quatrotate(Q, [uw vw ww]); % constant wind + +end % Relative velocities (plus wind); ur = u - wind(1); @@ -96,47 +73,33 @@ wr = w - wind(3); % Body to Inertial velocities Vels = quatrotate(Q_conj,[u v w]); - V_norm = norm([ur vr wr]); +%% CONSTANTS +% Everything related to empty condition (descent-fase) +S = settings.S; % [m^2] cross surface +C = settings.C; % [m] caliber +CoeffsE = settings.CoeffsE; % [/] Empty Rocket Coefficients +g = 9.80655; % [N/kg] module of gravitational field at zero -% 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; +mdot = 0; % Mass flow rate +Ixxdot = 0; % Inertia to x-axis rate +Iyydot = 0; % Inertia to y-axis rate +Izzdot = 0; % Inertia to z-axis rate +T = 0; % Thrust -%Atmosphere -if -z<0 +%% ATMOSPHERE DATA + +if -z<0 % z is directed as the gravity vector z = 0; end [~, a, ~, rho] = atmoscoesa(-z+settings.z0); M = V_norm/a; +%% AERODYNAMICS ANGLES -%Aero Angles -if not(u<1e-1 || V_norm<1e-3); +if not(u<1e-1 || V_norm<1e-3) alpha = atan(wr/ur); beta = asin(vr/V_norm); else @@ -144,183 +107,100 @@ else beta = 0; end -% beta = 0; % prova simulazione con beta imposto uguale a 0 - +%% DATCOM COEFFICIENTS - -% Coeffs. Interpolation givA = settings.Alphas*pi/180; givB = settings.Betas*pi/180; givH = settings.Altitudes; givM = settings.Machs; -% Interpolation error check +%% INTERPOLATION AT THE BOUNDARIES 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); +if M < givM(1) + 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); + + 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); +elseif alpha < givA(1) + alpha = givA(1); - - end end if beta > givB(end) - beta= givB(end); -else if beta < givB(1) + beta = givB(end); +elseif beta < givB(1) beta = givB(1); - end end -if -z >givH(end) +if -z > givH(end) z = -givH(end); -else if -z < givH(1) +elseif -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; +%% CHOSING THE CONDITION VALUE +% interpolation of the coefficients with the value in the nearest condition of the Coeffs matrix + +CA = interp4_easy(givA,givM,givB,givH,CoeffsE.CA,alpha,M,beta,-z); +CYB = interp4_easy(givA,givM,givB,givH,CoeffsE.CYB,alpha,M,beta,-z); +CNA = interp4_easy(givA,givM,givB,givH,CoeffsE.CNA,alpha,M,beta,-z); +Cl = interp4_easy(givA,givM,givB,givH,CoeffsE.CLL,alpha,M,beta,-z); +Clp = interp4_easy(givA,givM,givB,givH,CoeffsE.CLLP,alpha,M,beta,-z); +Cma = interp4_easy(givA,givM,givB,givH,CoeffsE.CMA,alpha,M,beta,-z); +Cmad = interp4_easy(givA,givM,givB,givH,CoeffsE.CMAD,alpha,M,beta,-z); +Cmq = interp4_easy(givA,givM,givB,givH,CoeffsE.CMQ,alpha,M,beta,-z); +Cnb = interp4_easy(givA,givM,givB,givH,CoeffsE.CLNB,alpha,M,beta,-z); +Cnr = interp4_easy(givA,givM,givB,givH,CoeffsE.CLNR,alpha,M,beta,-z); +Cnp = interp4_easy(givA,givM,givB,givH,CoeffsE.CLNP,alpha,M,beta,-z); + +%% FORCES +% first computed in the body-frame reference system +qdyn = 0.5*rho*V_norm^2; % [Pa] dynamics pressure +qdynL_V = 0.5*rho*V_norm*S*C; % +X = qdyn*S*CA; % [N] x-body component of the aerodynamics force +Y = qdyn*S*CYB*beta; % [N] y-body component of the aerodynamics force +Z = qdyn*S*CNA*alpha; % [N] z-body component of the aerodynamics force +Fg = quatrotate(Q,[0 0 m*g])'; % [N] force due to the gravity +F = Fg +[-X+T,+Y,-Z]'; % [N] total forces vector + +%% STATE DERIVATIVES + +% velocity +du = F(1)/m-q*w+r*v; +dv = F(2)/m-r*u+p*w; +dw = F(3)/m-p*v+q*u; % 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)... +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)... +dr = (Ixx-Iyy)/Izz*p*q + qdynL_V/Izz*(V_norm*Cnb*beta + (Cnr*r+Cnp*p)*C/2)... -Izzdot*r/Izz; +% Quaternion +OM = 1/2* [ 0 -p -q -r ; + p 0 r -q ; + q -r 0 p ; + r q -p 0 ]; -% 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; - +dQQ = OM*Q'; -OM = 1/2*[ - 0 -p -q -r - p 0 r -q - q -r 0 p - r q -p 0]; - - -dQQ = OM*Q'; +%% FINAL DERIVATIVE STATE ASSEMBLING dY(1:3) = Vels; dY(4) = du; diff --git a/src/config.m b/src/config.m index 2e00a77e..247d3a49 100644 --- a/src/config.m +++ b/src/config.m @@ -1,8 +1,5 @@ -%CONFIG SCRIPT - This script sets up all the parameters for the simulation (H1 line) +% CONFIG - This script sets up all the parameters for the simulation % 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 @@ -15,29 +12,23 @@ % email: francesco.colombi@skywarder.eu % Release date: 16/04/2016 -% clear all -% close all -% clc +%% LAUNCH SETUP -% ROCKET NAME +% rocket name settings.rocket_name = 'R2A'; -% LAUNCHPAD 6 -settings.z0 = 5; %Launchpad Altitude -settings.lrampa = 6.5; %LaunchPad route (launchpad lenght-distance from ground of the first hook) - +% launchpad 6 +settings.z0 = 5; %[m] Launchpad Altitude +settings.lrampa = 6.5; %[m] LaunchPad route (launchpad length-distance from ground of the first hook) -%STARTING ATTITUDE SETUP % -settings.OMEGA = 80*pi/180; %Elevation Angle -settings.PHI = 90*pi/180; %Azimuth Angle from North Direction -% ENGINE DETAILS % +% starting altitude +settings.OMEGA = 80*pi/180; %[rad] Elevation Angle, user input in degrees (ex. 80) +settings.PHI = 90*pi/180; %[rad] Azimuth Angle from North Direction, user input in degrees (ex. 90) -% % Thrust on time profile (IV Order Interpolation) -% settings.T_coeffs = [8034.5, 0, 0, 0, 0]; % constant thrust -% data from the motor details sheet (experimental plot of the Thrust) +%% ENGINE DETAILS -% SYNTAX: +% sintax: % engine = 1 -> Cesaroni PRO 150 White Thunder % engine = 2 -> Cesaroni PRO 150 SkidMark % engine = 3 -> Cesaroni PRO 150 BlueStreak @@ -46,104 +37,102 @@ switch engine case 1 % Cesaroni PRO 150 White Thunder % Sampling for thrust interpolation + settings.motor.Name = 'Cesaroni PRO 150 White Thunder'; 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 + 4 4.5 4.8 4.9 5 5.05 5.1 5.15 5.2]; % [s] settings.motor.exp_thrust = [8605.1 8900 7900 8400 8400 8250 8200 8300 ... - 8400 8400 8200 7800 7600 7450 7350 7300 4500 500 100 0]; %N - - settings.m0 = 67.761; %kg Overall Mass (Burnout + Propellant) - settings.ms = 43.961; %kg Structural Mass (Burnout - Nosecone) - settings.mp = 18.6; %kg Propellant Mass - settings.tb = 5.12; %sec Burning Time - settings.mfr = settings.mp/settings.tb; %kg/s Mass Flow Rate + 8400 8400 8200 7800 7600 7450 7350 7300 4500 500 100 0]; % [N] + + settings.m0 = 67.761; % [kg] Overall Mass (Burnout + Propellant) + settings.ms = 43.961; % [kg] Structural Mass (Burnout - Nosecone) + settings.mp = 18.6; % [kg] Propellant Mass + settings.tb = 5.12; % [s] 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.Name = 'Cesaroni PRO 150 SkidMark'; 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 + 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 = 64.9; %kg Overall Mass - settings.ms = 46.8; %kg Structural Mass (Burnout) - settings.mp = settings.m0-settings.ms; %kg Propellant Mass - settings.tb = 8.19; %sec Burning Time - settings.mfr = settings.mp/settings.tb; %kg/s Mass Flow Rate + 3800 4000 4081.6 3900 3800 3700 3500 3350 3200 3000 2000 750 0]; % [N] + + settings.m0 = 64.9; % [kg] Overall Mass + settings.ms = 46.8; % [kg] Structural Mass (Burnout) + settings.mp = settings.m0-settings.ms; % [kg] Propellant Mass + settings.tb = 8.19; % [s] 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.Name = 'Cesaroni PRO 150 BlueStreak'; 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 + 4 5 6 6.8 7.05 7.3 7.6 7.8]; % [s] settings.motor.exp_thrust = [0 800 4000 5500 5160 5130 5400 5300 5450 5347 ... - 5160 4950 4700 4400 4400 3800 300 0]; %N - - settings.m0 = 66.2; %kg Overall Mass - settings.ms = 47.3; %kg Structural Mass (Burnout) - settings.mp = settings.m0-settings.ms; %kg Propellant Mass - settings.tb = 7.60; %sec Burning Time - settings.mfr = settings.mp/settings.tb; %kg/s Mass Flow Rate + 5160 4950 4700 4400 4400 3800 300 0]; % [N] + + settings.m0 = 66.2; % [kg] Overall Mass + settings.ms = 47.3; % [kg] Structural Mass (Burnout) + settings.mp = settings.m0-settings.ms; % [kg] Propellant Mass + settings.mnc = 6.13; % [kg] Nosecone Mass + settings.tb = 7.60; % [s] Burning Time + settings.mfr = settings.mp/settings.tb; % [kg/s] Mass Flow Rate otherwise end -% GEOMETRY DETAILS % +%% GEOMETRY DETAILS % This parameters should be the same parameters set up in MISSILE DATCOM % simulation. -settings.C = 0.174; %m Caliber (Fuselage Diameter) -settings.S = 0.02378; %m2 Cross-sectional Surface -L = 4.4; %m Rocket length +settings.C = 0.174; % [m] Caliber (Fuselage Diameter) +settings.S = 0.02378; % [m^2] Cross-sectional Surface +L = 4.4; % [m] Rocket length -% MASS GEOMERTY DETAILS % +%% 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.27; %kg*m2 Inertia to x-axis (Full) -settings.Ixxe = 0.21; %kg*m2 Inertia to x-axis (Empty) -settings.Iyyf = 86.02; %kg*m2 Inertia to y-axis (Full) -settings.Iyye = 66.84; %kg*m2 Inertia to y-axis (Empty) -settings.Izzf = 86.02; %kg*m2 Inertia to z-axis (Full) -settings.Izze = 66.84; %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) -% The files are stored in the ../data folder with a naming convention of +% inertias for full configuration (with all the propellant embarqued) obtained with CAD's +settings.Ixxf = 0.27; % [kg*m^2] Inertia to x-axis +settings.Iyyf = 86.02; % [kg*m^2] Inertia to y-axis +settings.Izzf = 86.02; % [kg*m^2] Inertia to z-axis + +% inertias for empty configuration (all the propellant consumed) obtained with CAD's +settings.Ixxe = 0.21; % [kg*m^2] Inertia to x-axis +settings.Iyye = 66.84; % [kg*m^2] Inertia to y-axis +settings.Izze = 66.84; % [kg*m^2] Inertia to z-axis + +%% AERODYNAMICS DETAILS +% These coefficients are obtained using MISSILE DATCOM +% (after parsing with the tool datcom_parser.py) +% The files are stored in the ../data folder with the following rule: % rocket_name_full.mat | rocket_name_empty.mat -% e.g. R1X_full.mat etc.. +% e.g. R2a_full.mat | R2a_empty.mat +% Relative Path of the data files (default: ../data/). Remember the trailing slash!! -%Relative Path of the data files (default: ../data/). Remember the trailing -% slash!! +% Coeffs is a 4D matrix given by Datcom that contains the aerodynamics +% coefficient computed for the input parameters (AoA,Betas,Altitudes,Machs) +% Note: All the parameters (AoA,Betas,Altitudes,Machs) must be the same for +%empty and full configuration DATA_PATH = '../data/'; filename = strcat(DATA_PATH, settings.rocket_name); -%Coefficients in full configuration (with all the propellant embarqued) +% Coefficients in full configuration filename_full = strcat(filename,'_full.mat'); CoeffsF = load(filename_full,'Coeffs'); settings.CoeffsF = CoeffsF.Coeffs; clear('CoeffsF'); -%Coefficients in empty configuration (all the propellant consumed) +% Coefficients in empty configuration filename_empty = strcat(filename,'_empty.mat'); CoeffsE = load(filename_empty,'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(filename_full,'State'); settings.Alphas = s.State.Alphas'; settings.Betas = s.State.Betas'; @@ -152,109 +141,122 @@ settings.Machs = s.State.Machs'; clear('s'); -%PARACHUTES DETAILS % -%%% DROGUE 1 %%% -settings.para1.S = 1.55; %m2 Surface -settings.para1.mass = 0.577; %kg Parachute Mass -settings.para1.CD = 0.8; %Parachute Drag Coefficient -settings.para1.CL = 0; %Parachute Lift Coefficient +%% PARACHUTES DETAILS -%Altitude of Drogue 2 Opening -settings.zdrg2 = 5000; +% drogue 1 +settings.para1.S = 1.55; % [m^2] Surface +settings.para1.mass = 0.577; % [kg] Parachute Mass +settings.para1.CD = 0.8; % [/] Parachute Drag Coefficient +settings.para1.CL = 0; % [/] Parachute Lift Coefficient -%%% DROGUE 2 %%% -settings.para2.S = 17.5; %m2 Surface -settings.para2.mass = 1.140; %kg Parachute Mass -settings.para2.CD = 0.59; %Parachute Drag Coefficient -settings.para2.CL = 0; %Parachute Lift Coefficient +% drogue 2 +settings.para2.S = 17.5; % [m^2] Surface +settings.para2.mass = 1.140; % [kg] Parachute Mass +settings.para2.CD = 0.59; % [/] Parachute Drag Coefficient +settings.para2.CL = 0; % [/] Parachute Lift Coefficient +settings.zdrg2 = 5000; % [m] Altitude of drogue 2 opening -%Altitude of Main Parachute Opening -settings.zmain = 2000; +% main - rogallo wing +% The drogue parachute effects are neglected +settings.para3.S = 15; % [m^2] Surface +settings.para3.mass = 1.466; % [kg] Parachute Mass +settings.para3.CD = 0.4; % [/] Parachute Drag Coeff +settings.para3.CL = 0.8; % [/] Parachute Lift Coefficient +settings.zmain = 2000; % [m] Altitude of Main Parachute Opening -%%% MAIN - ROGALLO %%% -%The drogue parachute effects are neglected +%% INTEGRATION OPTIONS -settings.para3.S = 15; %m2 Surface -settings.para3.mass = 1.466; %kg Parachute Mass -settings.para3.CD = 0.4; %Parachute Drag Coeff -settings.para3.CL = 0.8; %Parachute Lift Coefficient +settings.ode.timeasc = 0:0.01:2000; % [s] Time span for ascend +settings.ode.timedrg1 = 0:0.01:2000; % [s] Time span for drogue 1 +settings.ode.timedrg2 = 0:0.01:2000; % [s] Time span for drogue 2 +settings.ode.timemain = 0:0.01:2000; % [s] Time span for main (rogallo) +settings.ode.timedesc = 0:0.01:2000; % [s] Time span for ballistic descent -% 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 + + +% create an option structure for the integrations: + +% - AbsTol is the threshold below which the value of the solution becomes unimportant +% - RelTol is the tolerance betweeen two consecutive values +% - Events is the event function that defines when the integration must be +% - stopped (it has to be created) +% - InitialStep is the highest value tried by the solver settings.ode.optionsasc = odeset('AbsTol',1E-3,'RelTol',1E-3,... - 'Events',@apogee,'InitialStep',1); %ODE options for ascend + 'Events',@event_apogee,'InitialStep',1); %ODE options for ascend settings.ode.optionsdrg1 = odeset('AbsTol',1E-3,'RelTol',1E-3,... - 'Events',@drg2_opening); %ODE options for drogue + 'Events',@event_drg2_opening); %ODE options for drogue settings.ode.optionsdrg2 = odeset('AbsTol',1E-3,'RelTol',1E-3,... - 'Events',@main_opening); %ODE options for drogue + 'Events',@event_main_opening); %ODE options for drogue settings.ode.optionsmain = odeset('AbsTol',1E-3,'RelTol',1E-12,... - 'Events',@crash); %ODE options for ascend + 'Events',@event_landing); %ODE options for descent settings.ode.optionsdesc = odeset('AbsTol',1E-3,'RelTol',1E-12,... - 'Events',@crash); %ODE options for ballistic descent + 'Events',@event_landing); %ODE options for ballistic descent -% WIND DETAILS % +%% WIND DETAILS -%Wind calculator model -settings.wind.model = true; -%true for hwsm wind model -%false for constant wind model +settings.wind.model = true; % Wind calculator model +% set to true for hwsm wind model +% set to false for constant wind model -% Wind is randomly generated. Setting the same values for min and max will -% fix the parameters of the wind. - -settings.wind.MagMin = 3; %Minimum Magnitude -settings.wind.MagMax = 3; %Maximum Magnitude -settings.wind.ElMin = 0*pi/180; %Minimum Elevation -settings.wind.ElMax = 0*pi/180; %Maximum Elevation (Max == 90 Deg) -settings.wind.AzMin = (90)*pi/180; %Minimum Azimuth -settings.wind.AzMax = (90)*pi/180; %Maximum Azimuth +% Wind is generated randomly from the minimum to the maximum parameters which defines the wind. +% Setting the same values for min and max will fix the parameters of the wind. +settings.wind.MagMin = 3; % [m/s] Minimum Magnitude +settings.wind.MagMax = 3; % [m/s] Maximum Magnitude +settings.wind.ElMin = 0*pi/180; % [rad] Minimum Elevation, user input in degrees (ex. 0) +settings.wind.ElMax = 0*pi/180; % [rad] Maximum Elevation, user input in degrees (ex. 0) (Max == 90 Deg) +settings.wind.AzMin = (90)*pi/180; % [rad] Minimum Azimuth, user input in degrees (ex. 90) +settings.wind.AzMax = (90)*pi/180; % [rad] Maximum Azimuth, user input in degrees (ex. 90) % NOTE: wind aziumt angle indications (wind directed towards): % 0 deg (use 360 instead of 0) -> North -% 90 deg -> East -% 180 deg -> South -% 270 deg -> West +% 90 deg -> East +% 180 deg -> South +% 270 deg -> West % Settings for the Wind Model -settings.wind.Lat = 39.552709; %Latitude of launching site -settings.wind.Long = 9.652400; %Longitude of launching site -settings.wind.Day = 290; %Day of the launch -settings.wind.Seconds = 13*60*60; %Second of the day (UTM) +settings.wind.Lat = 39.552709; % [deg] Latitude of launching site +settings.wind.Long = 9.652400; % [deg] Longitude of launching site +settings.wind.Day = 290; % [/] Day of the launch +settings.wind.Seconds = 13*60*60; % [/] Second of the day (UTM) +settings.wind.ww = 0; % [m/s] Vertical wind speed -% settings.wind.ww = vert_windgen(settings.wind.MagMin,settings.wind.MagMax); -%Vertical wind speed -settings.wind.ww = 0; % no vertical wind +%% BALLISTIC SIMULATION -% BALLISTIC SIMULATION -settings.ballistic = false; % Set to True to run a standard ballistic simulation +settings.ballistic = true; % Set to True to run a standard ballistic (without drogues) simulation -% LAST DROGUE FAILURE SIMULATION -% simulation in which rogallo wing does not open, therefore landing is +%% LAST DROGUE FAILURE SIMULATION +% simulation in which rogallo wing does not open and thus landing is % achieved thanks to the 2nd parachute + settings.ldf = false; -% APOGEE ONLY -% simulation stopped when reaching the apogee, therefore there is no -% descend phase. Only available for stochastic runs +%% APOGEE ONLY +% simulation stopped when reaching the apogee and thus there is no +% descend phase. Only available for stochastic runs!!! + settings.ao = false; -% STOCHASTIC DETAILS % -%If N>1 the stochastic routine is fired (different standard plots) -settings.stoch.N = 1; % Number of iterations -settings.stoch.parallel = false; % Using parallel or not parallel +%% STOCHASTIC DETAILS +% If N > 1 the stochastic routine is started (different standard plots) + +settings.stoch.N = 50; % Number of cases +settings.stoch.parallel = true; % Using parallel or not parallel tool + +%% PLOT DETAILS + +settings.ascend_plot = true; % Set to True to Plot with 'trend' plots (ascendplot.m) +settings.standard_plot = true; % Set to True to Plot with standard plots (standard_plot.m) +settings.default_plot = true; % Set to True to Plot with default plots (check the chosen simulation function, ex. std_run.m) +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 + +%% CLEARING VARIABLES NOT NEEDED -% 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 +clear('DATA_PATH','engine','filename','filename_empty','filename_full'); diff --git a/src/descent_parachute.m b/src/descent_parachute.m index 62c0c267..34971cd8 100644 --- a/src/descent_parachute.m +++ b/src/descent_parachute.m @@ -1,5 +1,5 @@ -function [ dY ] = descent_parachute( t,Y,settings,uw,vw,ww,para ) -% ODEFUN for Parachute descent +function [dY] = descent_parachute(t,Y,settings,uw,vw,ww,para) +% ODE-Function for Parachute descent % State = ( x y z | u v w ) % Author: Ruben Di Battista @@ -14,139 +14,101 @@ function [ dY ] = descent_parachute( t,Y,settings,uw,vw,ww,para ) % 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); +% 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); -% global bool -% % if bool = 0 the trends during the integration are NOT requested -% % if bool = 1 the trends during the integration are saved and plotted - - -% Adding Wind (supposed to be added in NED axes); if settings.wind.model -% Wind Model -[uw,vw,ww] = wind_generator(settings,z,t); -wind = [ uw,vw,ww ]; + [uw,vw,ww] = wind_matlab_generator(settings,z,t); + wind = [uw,vw,ww]; + else -% constant wind -wind = [uw vw ww]; + +wind = [uw vw ww]; % constant wind + end -% % Wind model -% h = 0; -% if -z > 0 -% h = -z+settings.z0; -% end -% -% % global alt -% % if bool == 1 -% % alt = [alt; h]; -% % end -% -% % wind in NED axis -% [uw,vw] = atmoshwm07(settings.wind.Lat, settings.wind.Long, h, ... -% 'day',settings.wind.Day,'model','quiet'); -% wind = [uw, vw, ww]; -% -% % global WIND -% % if bool == 1 -% % WIND = [WIND; uw, vw, ww]; -% % end - -% Adding wind; +% Relative velocities (plus wind); ur = u - wind(1); vr = v - wind(2); wr = w - wind(3); V_norm = norm([ur vr wr]); -% Constants +%% CONSTANTS + switch para case 1 - S = settings.para1.S; %Parachute Surface - CD = settings.para1.CD; %Parachute CD - CL = settings.para1.CL; %Parachute CL - pmass = 0; %detached parachute mass + S = settings.para1.S; % [m^2] Surface + CD = settings.para1.CD; % [/] Parachute Drag Coefficient + CL = settings.para1.CL; % [/] Parachute Lift Coefficient + pmass = 0; % [kg] detached mass case 2 - S = settings.para2.S; - CD = settings.para2.CD; - CL = settings.para2.CL; - pmass = settings.para1.mass + 6.13; %parachute + nosecone + S = settings.para2.S; % [m^2] Surface + CD = settings.para2.CD; % [/] Parachute Drag Coefficient + CL = settings.para2.CL; % [/] Parachute Lift Coefficient + pmass = settings.para1.mass + settings.mnc; % [kg] detached mass(drogue1 + nosecone) case 3 - S = settings.para3.S; - CD = settings.para3.CD; - CL = settings.para3.CL; - pmass = settings.para1.mass + settings.para2.mass + 6.13; + S = settings.para3.S; % [m^2] Surface + CD = settings.para3.CD; % [/] Parachute Drag Coefficient + CL = settings.para3.CL; % [/] Parachute Lift Coefficient + pmass = settings.para1.mass + settings.para2.mass + settings.mnc; % [kg] detached mass(drogue1/2 + nosecone) otherwise end -rad = sqrt(S/pi); %Equivalent parachute radius; - -%Side Surface -Ss = pi*(rad/2)^2; - -g = 9.81; +g = 9.80655; % [N/kg] magnitude of the gravitational field at zero +m = settings.ms - pmass; % [kg] descend mass -m = settings.ms - pmass; +%% ATMOSPHERE DATA -%Atmosphere -if -z<0 +if -z < 0 z = 0; end [~, ~, ~, rho] = atmoscoesa(-z); -%Center of Mass +%% REFERENCE FRAME +% The parachutes are approximated as rectangular surfaces with the normal +% vector perpendicular to the relative velocity +t_vect = [ur vr wr]; % Tangenzial vector +h_vect = [-vr ur 0]; % horizontal vector -%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); +if (h_vect(1) == 0) && (h_vect(2) == 0) + h_vect = [-v u 0]; end +t_vers = t_vect/norm(t_vect); % Tangenzial versor +h_vers = -h_vect/norm(h_vect); % horizontal versor +n_vect = cross(t_vers, h_vers); % Normal vector +n_vers = n_vect/norm(n_vect); % Normal versor -%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]'; - +if (n_vers(3) > 0) % If the normal vector is downward directed + n_vect = cross(h_vers, t_vers); + n_vers = n_vect/norm(n_vect); +end -du=F(1)/m; -dv=F(2)/m; -dw=F(3)/m; +%% FORCES +D = 0.5*rho*V_norm^2*S*CD*t_vers'; % [N] Drag vector +L = 0.5*rho*V_norm^2*S*CL*n_vers'; % [N] Lift vector +Fg = m*g*[0 0 1]'; % [N] Gravitational Force vector +F = -D+L+Fg; % [N] total forces vector +%% STATE DERIVATIVES +% velocity +du = F(1)/m; +dv = F(2)/m; +dw = F(3)/m; +%% FINAL DERIVATIVE STATE ASSEMBLING dY(1:3) = [u v w]'; dY(4) = du; diff --git a/src/event_apogee.m b/src/event_apogee.m new file mode 100644 index 00000000..edb27fa3 --- /dev/null +++ b/src/event_apogee.m @@ -0,0 +1,28 @@ +function [value,isterminal,direction] = event_apogee(t,Y,settings,varargin) +% Event function to stop simulation at apogee + +% Author: Ruben Di Battista +% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu +% email: ruben.dibattista@skywarder.eu +% Website: http://www.skywarder.eu +% April 2014; Last revision: 25.IV.2014 +% License: 2-clause BSD + +Q = Y(10:13)'; + +%Inertial Frame velocities +vels = quatrotate(quatconj(Q),Y(4:6)'); + +%Stop checking if I'm in Propulsion Phase +if t > settings.tb + value = vels(3); +else + value = 1; +end + +isterminal = 1; +direction = 1; + + +end + diff --git a/src/event_drg2_opening.m b/src/event_drg2_opening.m new file mode 100644 index 00000000..39b3ffe9 --- /dev/null +++ b/src/event_drg2_opening.m @@ -0,0 +1,14 @@ +function [value,isterminal,direction] = event_drg2_opening(t,Y,settings,varargin) +%Event that sets the main parachute opening altitude + +% Author: Francesco Colombi +% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu +% email: francesco.colombi@skywarder.eu +% Release date: 16/04/2016 + +value = Y(3) + settings.zdrg2; +isterminal = 1; +direction = 1; + +end + diff --git a/src/event_landing.m b/src/event_landing.m new file mode 100644 index 00000000..5ef1b381 --- /dev/null +++ b/src/event_landing.m @@ -0,0 +1,22 @@ +function [value,isterminal,direction] = event_landing(t,Y,settings,varargin) +% Event Function for ODE to stop simulation @ landing + +% Author: Ruben Di Battista +% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu +% email: ruben.dibattista@skywarder.eu +% Website: http://www.skywarder.eu +% April 2014; Last revision: 25.IV.2014 +% License: 2-clause BSD + +% Author: Francesco Colombi +% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu +% email: francesco.colombi@skywarder.eu +% Release date: 16/04/2016 + +value = Y(3); +isterminal = 1; +direction = 1; + + +end + diff --git a/src/event_main_opening.m b/src/event_main_opening.m new file mode 100644 index 00000000..aa3078e0 --- /dev/null +++ b/src/event_main_opening.m @@ -0,0 +1,16 @@ +function [value,isterminal,direction] = event_main_opening(t,Y,settings,varargin) +%Event that sets the main parachute opening altitude + + +% Author: Francesco Colombi +% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu +% email: francesco.colombi@skywarder.eu +% Release date: 16/04/2016 + +value = Y(3)+settings.zmain; +isterminal = 1; +direction = 1; + + +end + diff --git a/src/interp4_easy.m b/src/interp4_easy.m index ba3d1845..db1672e1 100644 --- a/src/interp4_easy.m +++ b/src/interp4_easy.m @@ -1,21 +1,13 @@ -function [ V ] = interp4_easy( X1,X2,X3,X4,TABLE,x1,x2,x3,x4 ) -%windgen( AzMin,AzMax,ElMin,ElMax,MagMin,MagMax ) +function [V] = interp4_easy(X1,X2,X3,X4,TABLE,x1,x2,x3,x4) +% interp4_easy: % 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 +% V = 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 @@ -24,7 +16,7 @@ function [ V ] = interp4_easy( X1,X2,X3,X4,TABLE,x1,x2,x3,x4 ) % License: 2-clause BSD M = {X1 X2 X3 X4}; -m = [ x1 x2 x3 x4]; +m = [x1 x2 x3 x4]; index = zeros(4,1); @@ -32,7 +24,7 @@ for i=1:4 [~,index(i)] = min(abs(M{i} - m(i))); end -V =TABLE(index(1),index(2),index(3),index(4)); +V = TABLE(index(1),index(2),index(3),index(4)); diff --git a/src/parfor_progress.m b/src/parfor_progress.m index 71bd7678..60ea01a2 100644 --- a/src/parfor_progress.m +++ b/src/parfor_progress.m @@ -1,5 +1,5 @@ function percent = parfor_progress(N) -%PARFOR_PROGRESS Progress monitor (progress bar) that works with parfor. +% 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 diff --git a/src/parfor_progress.txt b/src/parfor_progress.txt index afe234be..082b2efd 100644 --- a/src/parfor_progress.txt +++ b/src/parfor_progress.txt @@ -1,3 +1,51 @@ -2 +50 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 1 1 diff --git a/src/standard_plot.m b/src/standard_plot.m new file mode 100644 index 00000000..2caa57e9 --- /dev/null +++ b/src/standard_plot.m @@ -0,0 +1,196 @@ +% Author: Francesco Colombi +% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu +% email: francesco.colombi@skywarder.eu +% Release date: 16/04/2016 + + + %% ALTITUDE + + figure(); + plot(T, z), grid on, xlabel('time [s]'), ylabel('altitude [m]'), title('Altitude vs time'); + + %% VELOCITY + + figure(); + plot(T, vz), grid on, hold on; + plot(T, mod_V); + xlabel('time [s]'), ylabel('vz [m/s], |V| [m/s]'), title('Velocities vs time'); + legend('vz', '|V|'); + + %% ACCELERATION + + figure(); + plot(T, mod_A/9.80665), grid on; + xlabel('time [s]'), ylabel('|A| [g]'), title('Acceleration vs time'); + legend('|A|'); + + %% ALTITUDE,VELOCITY,ACCELERATION(subplotted) + + figure() + subplot(3,1,1) + plot(T, z), grid on, xlabel('time [s]'), ylabel('altitude [m]'), title('Altitude vs time'); + + subplot(3,1,2) + plot(T, vz), grid on, hold on; + plot(T, mod_V); + xlabel('time [s]'), ylabel('vz [m/s], |V| [m/s]'); + legend('vz', '|V|'), title('Velocities vs time'); + + subplot(3,1,3) + plot(T, mod_A/9.80665), grid on; + xlabel('time [s]'), ylabel('|A| [g]'); + legend('|A|'); + title('Acceleration vs time') + + %% 3D TRAJECTORY + + 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') + % concentric circles distance + 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]') + + %% TRAECTORY'S TOP VIEW + + figure(); + plot(y, x), axis equal, hold on, grid on; + p1 = plot(0,0,'r.','MarkerSize',21); + p2 = plot(y(end),x(end),'rx','MarkerSize',14); + title('Trajectory Top view') + xlabel('y, East [m]'), ylabel('x, North [m]'); + legend([p1,p2],{'Launch Point','Landing Point'}) + clear p1 p2 + + % 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 + + %% TRAJECTORY PROJECTIONS(subplotted) + + figure() + suptitle('Trajectory projections') + subplot(1,3,1) + plot(y_flight/1000, x_flight/1000), axis equal, hold on, grid on; + p11 = plot(Ya(end,2)/1000, Ya(end,1)/1000, 'r*','markersize',7); + p12 = plot(0, 0, 'r.','markersize',14); + p13 = plot(Y(end,2)/1000, Y(end,1)/1000, 'rx','markersize',7); + 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, 'r*','markersize',7); + plot(Y(end,1)/1000, -Y(end,3)/1000, 'rx','markersize',7); + plot(0, 0, 'r.','markersize',14); + 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, 'r*','markersize',7); + plot(Y(end,2)/1000, -Y(end,3)/1000, 'rx','markersize',7); + plot(0, 0, 'r.','markersize',14); + xlabel('y, East [Km]'), ylabel('z, Altitude [Km]'); + %setting limit if is parallel to north + if sum(y_flight) == 0 + xlim([-1 1]); + end + legend([p11, p12, p13],{'Apogee','Launch Point','Landing Point'}) + + %% ANGULAR RATES (combined) + + 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 [s]'), ylabel('p, q, r [grad/s]') + legend('p: roll rate', 'q: pitch rate', 'r: yaw rate') + title('Angular rates vs time') + + %% ANGULAR RATES (subplotted) + + figure(),suptitle('Angular rates vs time'); + subplot(3,1,1) + plot(Ta, Ya(:, 8)*180/pi), grid on; + xlabel('time [s]'), ylabel('pitch rate q [grad/s]') + + subplot(3,1,2) + plot(Ta, Ya(:, 9)*180/pi), grid on; + xlabel('time [s]'), ylabel('yaw rate r [grad/s]') + + subplot(3,1,3) + plot(Ta, Ya(:, 7)*180/pi), grid on; + xlabel('time [s]'), ylabel('roll rate p [grad/s]') + + + %% EULERIAN ANGLES (3 figures) + + pitch_angle = zeros(length(Ta),1); + yaw_angle = zeros(length(Ta),1); + roll_angle = zeros(length(Ta),1); + + for k = 2:length(Ta) + pitch_angle(k) = pitch_angle(k-1) + (Ya(k, 8) + Ya(k-1, 8))/2*180/pi*(T(k)-T(k-1)); + yaw_angle(k) = yaw_angle(k-1) + (Ya(k, 9) + Ya(k-1, 9))/2*180/pi*(T(k)-T(k-1)); + roll_angle(k) = roll_angle(k-1) + (Ya(k, 7) + Ya(k-1, 7))/2*180/pi*(T(k)-T(k-1)); + end + + figure() + plot(Ta, pitch_angle+settings.OMEGA*180/pi) + title('Pitch angle vs time') + grid on, xlabel('time [s]'), ylabel('pitch angle [deg]'); + + figure() + plot(Ta, yaw_angle) + title('Yaw angle vs time') + grid on, xlabel('time [s]'), ylabel('yaw angle [deg]'); + + figure() + plot(Ta, roll_angle) + title('Roll angle vs time') + grid on, xlabel('time [s]'), ylabel('roll angle [deg]') + + + %% MACH NUMBER + + figure(), plot(T, M), grid on; + xlabel('time [s]'), ylabel('Mach [/]'), title('Mach number vs time'); + + %% STAGNATION TEMPERATURE + + figure(); + plot(T, Tamb-273.15); + hold on, grid on; + plot(T, Ttot-273.15) + title('Temperature profile') + xlabel('time [s]'), ylabel('Temperature [C]'); + legend('Surrounding', 'Total', 'location', 'best') + + clear variable + \ No newline at end of file diff --git a/src/start_simulation.m b/src/start_simulation.m index 215e7a36..b1ceeb47 100644 --- a/src/start_simulation.m +++ b/src/start_simulation.m @@ -4,141 +4,104 @@ % Release date: 16/04/2016 close all -clear, clc - -%% caricamento dati -run('config.m'); - -global bool contatore; -global t_plot T_plot alpha_plot beta_plot M_plot CA_plot Drag_plot Forces_plot; -% global WIND alt; -global XCP; - - -% 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 = []; - XCP = 0; -end +clear all +clc + +%% LOAD DATA + +run('config.m'); + +%% START THE CHOSEN SIMULATION +% T = vector of time used by ODE, [s] also for Tf Ta +% Y = State = ( x y z | u v w | p q r | q0 q1 q2 q3 ) also for Ya,Yf corresponding to T -%% lancio simulatore -[T,Y, Ta,Ya] = MAIN(settings); - -if (bool == 1) - figure(); - plot(t_plot, T_plot, '.'), title('Thrust vs time'), grid on; - ylabel('Thrust [N]') - figure(); - plot(t_plot, alpha_plot*180/pi), title('alpha vs time'), grid on; - ylabel('alpha [deg]') - figure(); - plot(t_plot, beta_plot*180/pi), title('beta vs time'), grid on; - ylabel('beta [deg]') - - figure(); - plot(t_plot, M_plot), title('mach vs time'), grid on; - ylabel('Mach M [-]') - - figure(); - plot(t_plot, CA_plot), title('Aerodyn Coeff vs time'), grid on; - ylabel('Drag Coeff CD [-]') - - figure(); - plot(t_plot, Drag_plot), title('Drag vs time'), grid on; - ylabel('Drag D [N]') - - figure(); - plot(t_plot, Forces_plot), title('Axial Force vs time'), grid on; - ylabel('Axial force [N]') - -% figure() -% title('Wind'), grid on, hold on; -% plot(WIND(:, 1), alt) -% plot(WIND(:, 2), alt) -% plot(WIND(:, 3), alt) -% xlabel('Wind magnitude [m/s]') -% ylabel('Altitude [m]') -% legend('North','East','Down') - - figure(); - plot(t_plot, XCP,'.'), title('Stability margin vs time'), grid on; - ylabel('s.m.') +tic +% Checking if stochastic or standard simulation needed +if settings.ballistic + if settings.stoch.N > 1 + fprintf('Stochastic Ballistic Simulation Started...\n\n'); + if settings.stoch.parallel + [LP,Z] = stoch_run_bal_p(settings); + else + [LP,Z] = stoch_run_bal(settings); + end + else + fprintf('Standard Ballistic Simulation Started...\n\n'); + [T,Y, Ta,Ya] = std_run_ballistic(settings); + end +else + if settings.stoch.N > 1 + fprintf('Stochastic Simulation Started...\n\n'); + if settings.stoch.parallel + [LP,Z] = stoch_run_p(settings); + else + [LP,Z] = stoch_run(settings); + end + else + fprintf('Standard Simulation Started...\n\n'); + [T,Y, Ta,Ya] = std_run(settings); + end end +toc +if settings.stoch.N == 1 -[apogee, i_apogee] = max(-Y(:,3)); +%% POSITIONS -%% posizione x = Y(:,1); y = Y(:,2); z = -Y(:,3); X = [x, y, z]; +[apogee, i_apogee] = max(-Y(:,3)); % position, index of position at apogee + +%% VELOCITIES -%% 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 +%% ACCELERATIONS + +% main derivatives 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 + +% add derivative at the boundaries 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 +clear('ax', 'ay', 'az', 'vx', 'vy') +%% MAXIMUM POSITIONS, VELOCITIES AND ACCELERATIONS + +% pre-allocation mod_X = zeros(length(X), 1); mod_V = mod_X; mod_A = mod_X; +% determine the module of every row element 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); +iexit = find(mod_X <= settings.lrampa); % checking where the missile is undocked from the hook of the launch pad +iexit = iexit(end); + +%% TEMPERATURE AND MACH NUMBER +% pre-allocation M = zeros(length(X), 1); Tamb = M; Ttot = Tamb; @@ -147,19 +110,15 @@ for n = 1:length(M) if (h < 0) h = 0; end - [Tamb(n), a, ~, ~] = atmoscoesa(h); + [Tamb(n), a, ~, ~] = atmosisa(h); M(n) = mod_V(n)/a; Ttot(n) = Tamb(n)*(1+0.2*M(n)^2); end - +% determine the maximum Mach number [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); @@ -184,173 +143,16 @@ 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 +%% PLOT - 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 - p1 = plot(0,0,'r.','MarkerSize',21); - p2 = plot(y(end),x(end),'rx','MarkerSize',14); - title('Trajectory Top view') - xlabel('y, East [m]'), ylabel('x, North [m]'); - legend([p1,p2],{'Launch Point','Landing Point'}) - clear p1 p2 - - % 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 +if settings.standard_plot + run('standard_plot.m') +end - figure(); - subplot(1,3,1); - plot(y_flight/1000, x_flight/1000), axis equal, hold on, grid on; - p11 = plot(Ya(end,2)/1000, Ya(end,1)/1000, 'r*','markersize',7); - p12 = plot(0, 0, 'r.','markersize',14); - p13 = plot(Y(end,2)/1000, Y(end,1)/1000, 'rx','markersize',7); - 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, 'r*','markersize',7); - plot(Y(end,1)/1000, -Y(end,3)/1000, 'rx','markersize',7); - plot(0, 0, 'r.','markersize',14); - 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, 'r*','markersize',7); - plot(Y(end,2)/1000, -Y(end,3)/1000, 'rx','markersize',7); - plot(0, 0, 'r.','markersize',14); - xlabel('y, East [Km]'), ylabel('z, Altitude [Km]'); - %setting limit if is parallel to north - if sum(y_flight) == 0 - xlim([-1 1]); - end - legend([p11, p12, p13],{'Apogee','Launch Point','Landing Point'}) - - % angular rates - figure(); - plot(Ta, Ya(:, 7)*180/pi) - hold on, grid on - plot(Ta, Ya(:, 8)*180/pi) - plot(Ta, Ya(:, 9)*180/pi) - xlabel('time [sec]'), ylabel('p, q, r [grad/sec]') - legend('p: roll rate', 'q: pitch rate', 'r: yaw rate') - title('angular rates vs time') - - figure(); - subplot(3,1,1) - plot(Ta, Ya(:, 8)*180/pi), grid on; - xlabel('time [sec]'), ylabel('pitch rate q [grad/sec]') - subplot(3,1,2) - plot(Ta, Ya(:, 9)*180/pi), grid on; - xlabel('time [sec]'), ylabel('yaw rate r [grad/sec]') - subplot(3,1,3) - plot(Ta, Ya(:, 7)*180/pi), grid on; - xlabel('time [sec]'), ylabel('roll rate p [grad/sec]') - - - % angle - alpha = zeros(length(Ta),1); - beta = zeros(length(Ta),1); - roll = zeros(length(Ta),1); - - for k = 2:length(Ta) - alpha(k) = alpha(k-1) + (Ya(k, 8) + Ya(k-1, 8))/2*180/pi*(T(k)-T(k-1)); - beta(k) = beta(k-1) + (Ya(k, 9) + Ya(k-1, 9))/2*180/pi*(T(k)-T(k-1)); - roll(k) = roll(k-1) + (Ya(k, 7) + Ya(k-1, 7))/2*180/pi*(T(k)-T(k-1)); - end - figure(), plot(Ta, alpha+settings.OMEGA*180/pi), title('Alpha vs time'), grid on; - figure(), plot(Ta, beta), title('Beta vs time'), grid on; - figure(), plot(Ta, roll), title('Roll vs time'), grid on; - - - % Mach - figure(), plot(T, M), grid on; - xlabel('time [sec]'), ylabel('Mach [-]'); - - % Stagnation Temperature - figure(); - plot(T, Tamb-273.15); - hold on, grid on; - plot(T, Ttot-273.15) - title('Temperature profile') - xlabel('time [sec]'), ylabel('Temperature [�C]'); - legend('Surrounding', 'Total', 'location', 'best') +if settings.ascend_plot + run('ascend_plot.m') +end + +clear('ascend_plot.mat') end diff --git a/src/std_run.m b/src/std_run.m index 224b7020..b0efdee0 100644 --- a/src/std_run.m +++ b/src/std_run.m @@ -1,5 +1,5 @@ -function [ Tf,Yf, Ta,Ya ] = std_run( settings ) -%STD RUN - This function runs a standard (non-stochastic) simulation +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 @@ -16,110 +16,109 @@ function [ Tf,Yf, Ta,Ya ] = std_run( settings ) % email: francesco.colombi@skywarder.eu % Release date: 16/04/2016 -%Starting Attitude +%% STARTING CONDITIONS + +% Attitude Q0 = angle2quat(settings.PHI,settings.OMEGA,0*pi/180,'ZYX')'; -%Starting State +% 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 -if settings.wind.model +%% WIND GENERATION + +if settings.wind.model % will be computed inside the integrations uw = 0; vw = 0; ww = 0; else - [uw,vw,ww] = windgen(settings.wind.AzMin,settings.wind.AzMax,... + [uw,vw,ww] = wind_const_generator(settings.wind.AzMin,settings.wind.AzMax,... settings.wind.ElMin,settings.wind.ElMax,settings.wind.MagMin,... settings.wind.MagMax); end -%% ASCEND %% + +%% ASCEND [Ta,Ya] = ode45(@ascend,settings.ode.timeasc,X0a,settings.ode.optionsasc,... settings,uw,vw,ww); -%% DROGUE 1 %% -para = 1; %Flag for Drogue 1 +%% DROGUE 1 +% Initial Condition are the last from ascend (need to rotate because +% velocities are outputted in body axes) -%Initial Condition are the last from ascend (need to rotate because -%velocities are outputted in body axes) +para = 1; % Flag for Drogue 1 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 %% +%% DROGUE 2 +% Initial Condition are the last from drogue 1 descent -para = 2; %Flag for Drogue 2 - -%Initial Condition are the last from drogue 1 descent +para = 2; % Flag for Drogue 2 X0d2 = Yd1(end,:); [Td2,Yd2] = ode45(@descent_parachute,settings.ode.timedrg2,X0d2,... settings.ode.optionsdrg2,settings,uw,vw,ww,para); -%% MAIN %% +%% MAIN - ROGALLO WING +% Initial Condition are the last from drogue 2 descent + if not(settings.ldf) -para = 3; %Flag for Main (Rogall) +para = 3; % Flag for Main (Rogall) end -%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 %% +%% FINAL STATE ASSEMBLING -%Total State +% Total State Yf = [Ya(:,1:3) quatrotate(quatconj(Ya(:,10:13)),Ya(:,4:6));Yd1; Yd2;Yrog]; -%global Par1 Par2 Par3 Asc -%Asc = [Ya(:,1:3) quatrotate(quatconj(Ya(:,10:13)),Ya(:,4:6))]; -%Par1 = Yd1; -%Par2 = Yd2; -%Par3 = Yrog; - -%Total Time +% Total Time Tf = [Ta; Ta(end)+Td1; Ta(end)+Td1(end)+Td2;Ta(end)+Td1(end)+Td2(end)+Trog]; -%Parachute State +% Parachute State Yp = [Yd1;Yd2;Yrog]; Tp = [Ta(end)+Td1;Ta(end)+Td1(end)+Td2;Ta(end)+Td1(end)+Td2(end)+Trog]; -%% PLOTTING THINGS %% +%% DEFAULT PLOTS -if settings.plot == 1 +if settings.default_plot set(0,'DefaultAxesFontSize',settings.DefaultFontSize,... 'DefaultLineLineWidth',settings.DefaultLineWidth); - % ASCENT % + % ascend 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,... + h1 = plot(Ta(end),-Ya(end,3)+settings.z0,'ko','MarkerSize',8,... 'MarkerFaceColor','k'); - h3=plot(Ta(end)+Td1(end)+Td2(end),-Yd2(end,3)+settings.z0,'ks','MarkerSize',8,... + 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 + % Interpolation for less points --> visualization Tinterp = linspace(Tf(1),Tf(end),settings.tSteps); - [Tunique,ia,~]= unique(Tf); + [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 + % Velocities + Va = quatrotate(quatconj(Ya(end,10:13)),Ya(end,4:6)); figure; suptitle('Velocities Profiles on Time') @@ -128,11 +127,11 @@ if settings.plot == 1 hold on xlabel('Time[s]'); ylabel('Velocity-x [m/s]'); - h1=plot(Ta(end),Va(1),'ko','MarkerSize',8,... + h1 = plot(Ta(end),Va(1),'ko','MarkerSize',8,... 'MarkerFaceColor','k'); - h2=plot(Ta(end)+Td1(end),Yd1(end,4),'ks','MarkerSize',8,... + 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,... + 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',... @@ -144,11 +143,11 @@ if settings.plot == 1 hold on xlabel('Time[s]'); ylabel('Velocity-y [m/s]'); - h1=plot(Ta(end),Va(2),'ko','MarkerSize',8,... + h1 = plot(Ta(end),Va(2),'ko','MarkerSize',8,... 'MarkerFaceColor','k'); - h2=plot(Ta(end)+Td1(end),Yd1(end,5),'ks','MarkerSize',8,... + 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,... + 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',... @@ -161,67 +160,66 @@ if settings.plot == 1 hold on xlabel('Time[s]'); ylabel('Velocity-z [m/s]'); - h1=plot(Ta(end),Va(3),'ko','MarkerSize',8,... + h1 = plot(Ta(end),Va(3),'ko','MarkerSize',8,... 'MarkerFaceColor','k'); - h2=plot(Ta(end)+Td1(end),Yd1(end,6),'ks','MarkerSize',8,... + 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,... + 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% + % Trajectory figure; - h0=plot3(Yf(1,2),Yf(1,1),-Yf(1,3)+settings.z0,'k+','MarkerSize',10); + 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,... + 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,... + 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,... + 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 % - + % 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,... + h1 = plot(Ya(end,2),Ya(end,1),'ko','MarkerSize',8,... 'MarkerFaceColor','k'); - h2=plot(Yd1(end,2),Yd1(end,1),'ks','MarkerSize',8,... + h2 = plot(Yd1(end,2),Yd1(end,1),'ks','MarkerSize',8,... 'MarkerFaceColor','k'); - h3=plot(Yd2(end,2),Yd2(end,1),'ks','MarkerSize',8,... + 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 % + % 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,... + 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,... + 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,... + 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',... @@ -230,14 +228,13 @@ if settings.plot == 1 - %Interpolation for less points --> visualization + % 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 % - + % Angular rates figure; suptitle('Angular rates on Time') subplot(3,1,1); @@ -245,7 +242,7 @@ if settings.plot == 1 hold on xlabel('Time[s]'); ylabel('Roll rate p [rad/s]'); - h1=plot(Ta(end),Ya(end,7),'ko','MarkerSize',8,... + h1 = plot(Ta(end),Ya(end,7),'ko','MarkerSize',8,... 'MarkerFaceColor','k'); legend(h1,'Apogee/1st Drogue Deployment','Location','southeast'); grid on @@ -255,7 +252,7 @@ if settings.plot == 1 hold on xlabel('Time[s]'); ylabel('Pitch rate q [rad/s]'); - h1=plot(Ta(end),Ya(end,8),'ko','MarkerSize',8,... + h1 = plot(Ta(end),Ya(end,8),'ko','MarkerSize',8,... 'MarkerFaceColor','k'); legend(h1,'Apogee/1st Drogue Deployment','Location','northeast'); grid on @@ -266,18 +263,18 @@ if settings.plot == 1 hold on xlabel('Time[s]'); ylabel('Yaw rate r [rad/s]'); - h1=plot(Ta(end),Ya(end,9),'ko','MarkerSize',8,... + h1 = plot(Ta(end),Ya(end,9),'ko','MarkerSize',8,... 'MarkerFaceColor','k'); legend(h1,'Apogee/1st Drogue Deployment','Location','southeast'); grid on end -%Resizing +% Resizing h = get(0,'children'); scrsz = get(0,'ScreenSize'); -for i=1:length(h) +for i = 1:length(h) set(h(i),'OuterPosition',[0 0 scrsz(4) scrsz(4)]) - %saveas(h(i), ['figure' num2str(i)], 'fig'); end + end diff --git a/src/std_run_ballistic.m b/src/std_run_ballistic.m index 2bf83adb..7ea8bada 100644 --- a/src/std_run_ballistic.m +++ b/src/std_run_ballistic.m @@ -1,5 +1,5 @@ -function [ Tf,Yf, Ta,Ya ] = std_run_ballistic( settings ) -%STD RUN - This function runs a standard (non-stochastic) simulation +function [Tf,Yf,Ta,Ya] = std_run_ballistic(settings) +% STD RUN BALLISTIC - This function runs a ballistic (non-stochastic) simulation % OUTPUTS % Tf: Time steps % Yf: Final State @@ -16,93 +16,71 @@ function [ Tf,Yf, Ta,Ya ] = std_run_ballistic( settings ) % email: francesco.colombi@skywarder.eu % Release date: 16/04/2016 -%Starting Attitude +%% STARTING CONDITIONS + +% Attitude Q0 = angle2quat(settings.PHI,settings.OMEGA,0*pi/180,'ZYX')'; -%Starting State +% 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,... +%% WIND GENERATION + +if settings.wind.model % will be computed inside the integrations + uw = 0; + vw = 0; + ww = 0; +else + [uw,vw,ww] = wind_const_generator(settings.wind.AzMin,settings.wind.AzMax,... settings.wind.ElMin,settings.wind.ElMax,settings.wind.MagMin,... settings.wind.MagMax); +end -%% ASCEND %% +%% ASCEND [Ta,Ya] = ode45(@ascend,settings.ode.timeasc,X0a,settings.ode.optionsasc,... settings,uw,vw,ww); -%% DESCEND %% +%% DESCEND +% ballistic descend, so no drogue will be used [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 +%% FINAL STATE ASSEMBLING + +% Total State Yf = [Ya;Yd]; -%Total Time +% Total Time Tf = [Ta; Ta(end)+Td]; -%% PLOTTING THINGS %% +%% DEFAULT PLOTS -if settings.plot == 1 +if settings.default_plot == 1 set(0,'DefaultAxesFontSize',settings.DefaultFontSize,... 'DefaultLineLineWidth',settings.DefaultLineWidth); - %Interpolation for less points --> visualization + % 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; - subtitle('Velocities Profiles on Time') + 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'); + h1 = plot(Ta(end),Ya(end,4),'ko','MarkerSize',8,'MarkerFaceColor','k'); + h2 = plot(Tf(end),Yd(end,4),'ks','MarkerSize',8,'MarkerFaceColor','k'); + legend([h1,h2],'Apogee','Landing Point','Location','southeast'); grid on subplot(3,1,2) @@ -110,13 +88,9 @@ if settings.plot == 1 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'); + h1 = plot(Ta(end),Ya(end,5),'ko','MarkerSize',8,'MarkerFaceColor','k'); + h2 = plot(Tf(end),Yd(end,5),'ks','MarkerSize',8,'MarkerFaceColor','k'); + legend([h1,h2],'Apogee','Landing Point','Location','southeast'); grid on @@ -125,62 +99,48 @@ if settings.plot == 1 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'); + h1 = plot(Ta(end),Ya(end,6),'ko','MarkerSize',8,'MarkerFaceColor','k'); + h2 = plot(Tf(end),Yd(end,6),'ks','MarkerSize',8,'MarkerFaceColor','k'); + legend([h1,h2],'Apogee','Landing Point','Location','southeast'); grid on - %COMPLETE TRAJECTORY% + % Trajectory figure; - h0=plot3(Yf(1,2),Yf(1,1),-Yf(1,3)+settings.z0,'k+','MarkerSize',10); + 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,... + 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,... + 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'); + legend([h0,h1,h2],'Launch Pad','Apogee','Landing Point'); grid on - % PLANAR DISPLACEMENT % - + % 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'); + 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'); + legend([h1,h2],'Apogee','Landing Point','Location','southeast'); grid on - - - - %Interpolation for less points --> visualization + % 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 % - + % Angular rates figure; suptitle('Angular rates on Time') subplot(3,1,1); @@ -188,9 +148,8 @@ if settings.plot == 1 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'); + h1 = plot(Ta(end),Ya(end,7),'ko','MarkerSize',8,'MarkerFaceColor','k'); + legend(h1,'Apogee','Location','southeast'); grid on subplot(3,1,2) @@ -198,9 +157,8 @@ if settings.plot == 1 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'); + h1 = plot(Ta(end),Ya(end,8),'ko','MarkerSize',8,'MarkerFaceColor','k'); + legend(h1,'Apogee','Location','northeast'); grid on @@ -209,21 +167,17 @@ if settings.plot == 1 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'); + h1 = plot(Ta(end),Ya(end,9),'ko','MarkerSize',8,'MarkerFaceColor','k'); + legend(h1,'Apogee','Location','southeast'); grid on end -%Resizing +% Resizing h = get(0,'children'); scrsz = get(0,'ScreenSize'); for i=1:length(h) set(h(i),'OuterPosition',[0 0 scrsz(4) scrsz(4)]) - %saveas(h(i), ['figure' num2str(i)], 'fig'); end - - end \ No newline at end of file diff --git a/src/stoch_run.m b/src/stoch_run.m index eb17b5c8..7342b71f 100644 --- a/src/stoch_run.m +++ b/src/stoch_run.m @@ -1,5 +1,5 @@ -function [ LP ,Z ] = stoch_run( settings ) -%STD RUN - This function runs a stochastic simulation +function [LP,Z] = stoch_run(settings) +% STOCH RUN - This function runs a standard stochastic simulation % OUTPUTS % LP: Landing Points % Z: Apogee Altitudes @@ -11,124 +11,125 @@ function [ LP ,Z ] = stoch_run( settings ) % April 2014; Last revision: 29.V.2014 % License: 2-clause BSD -%Starting Attitude +%% STARTING CONDITIONS + +% Attitude Q0 = angle2quat(settings.PHI,settings.OMEGA,0*pi/180,'ZYX')'; -%Starting State +% 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); +%% FOR LOOP -parfor_progress(settings.stoch.N); for i=1:settings.stoch.N - %% Wind Generation - [uw,vw,ww] = windgen(settings.wind.AzMin,settings.wind.AzMax,... + + %% WIND GENERATION + + [uw,vw,ww] = wind_const_generator(settings.wind.AzMin,settings.wind.AzMax,... settings.wind.ElMin,settings.wind.ElMax,settings.wind.MagMin,... settings.wind.MagMax); - %% ASCEND %% + %% ASCEND [Ta,Ya] = ode45(@ascend,settings.ode.timeasc,X0a,settings.ode.optionsasc,... settings,uw,vw,ww); - %% DROGUE 1 %% - para = 1; %Flag for Drogue 1 + %% DROGUE 1 + % Initial Condition are the last from ascend (need to rotate because + % velocities are outputted in body axes) - %Initial Condition are the last from ascend (need to rotate because - %velocities are outputted in body axes) + para = 1; % Flag for Drogue 1 X0d1 = [Ya(end,1:3) quatrotate(quatconj(Ya(end,10:13)),Ya(end,4:6))]; - [Td1,Yd1] = ode45(@descent_parachute,settings.ode.timedrg1,X0d1,... + [~,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 + + para = 2; %Flag for Drogue 2 X0d2 = Yd1(end,:); - [Td2,Yd2] = ode45(@descent_parachute,settings.ode.timedrg2,X0d2,... + [~,Yd2] = ode45(@descent_parachute,settings.ode.timedrg2,X0d2,... settings.ode.optionsdrg2,settings,uw,vw,ww,para); - %% MAIN %% + %% MAIN - ROGALLO WING + % Initial Condition are the last from drogue 2 descent + if not(settings.ldf) - para = 3; %Flag for Main (Rogall) + para = 3; % Flag for Main (Rogall) end - %Initial Condition are the last from drogue 2 descent X0m = Yd2(end,:); - [Tm,Ym] = ode45(@descent_parachute,settings.ode.timemain,X0m,... + [~,Ym] = ode45(@descent_parachute,settings.ode.timemain,X0m,... settings.ode.optionsmain,settings,uw,vw,ww,para); - %% FINAL STATE ASSEMBLING %% + %% FINAL STATE ASSEMBLING - %Total State + % Total State if settings.ao Yf = [Ya(:,1:3) quatrotate(quatconj(Ya(:,10:13)),Ya(:,4:6))]; else Yf = [Ya(:,1:3) quatrotate(quatconj(Ya(:,10:13)),Ya(:,4:6));Yd1; Yd2;Ym]; end - %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,:)) +%% CHECKING BAD SIMULATION + +if numel(LP(LP(:,3) < -10,:)) fprintf(['Some landing points might be incorrect' ... 'please check parameters!\n']); end -% Writing Things +%% PRINTING VALUES -%Mean Landing Point +% Mean Landing Point xm = mean(LP(:,1)); ym = mean(LP(:,2)); -%Mean Apogee Time +% Mean Apogee Time ApoTimem = mean(ApoTime); -%Std. Deviation Apogee Time +% Std. Deviation Apogee Time ApoTimestd = std(ApoTime); -%Mean Altitude +% Mean Altitude Zm = mean(Z); -%Std. Deviation Altitude +% Std. Deviation Altitude Zstd = std(Z); -%Best Fit Ellipse +% 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',... +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 - +%% DEFAULT PLOTS + +if settings.default_plot == 1 plot(xm,ym,'bs','MarkerSize',20,'MarkerFacecolor','b'); hold on - %Point of launch + % Point of launch plot(0,0,'ro','MarkerSize',20,'MarkerFacecolor','r'); - %All the landing points + % All the landing points plot(LP(:,1),LP(:,2),'k+'); title('Landing Points'); @@ -138,7 +139,7 @@ if settings.plot == 1 view(90,270) axis equal - %Histogram + % Histogram [f,x] = hist(Z,10); figure; bar(x,f/settings.stoch.N); @@ -146,17 +147,13 @@ if settings.plot == 1 xlabel('Apogee [m]') ylabel('n_i/n_{tot}') - - - end -%Resizing +% Resizing h = get(0,'children'); scrsz = get(0,'ScreenSize'); -for i=1:length(h) +for i = 1:length(h) set(h(i),'OuterPosition',[0 0 scrsz(4) scrsz(4)]) - %saveas(h(i), ['figure' num2str(i)], 'fig'); end diff --git a/src/stoch_run_bal.m b/src/stoch_run_bal.m index 8292038e..0bc19458 100644 --- a/src/stoch_run_bal.m +++ b/src/stoch_run_bal.m @@ -1,5 +1,5 @@ -function [ LP ,Z ] = stoch_run_bal( settings ) -%STD RUN - This function runs a stochastic simulation +function [LP ,Z] = stoch_run_bal(settings) +% STOCH RUN BAL - This function runs a stochastic ballistic simulation % OUTPUTS % LP: Landing Points % Z: Apogee Altitudes @@ -11,101 +11,101 @@ function [ LP ,Z ] = stoch_run_bal( settings ) % April 2014; Last revision: 29.V.2014 % License: 2-clause BSD -%Starting Attitude +%% STARTING CONDITIONS + +% Attitude Q0 = angle2quat(settings.PHI,settings.OMEGA,0*pi/180,'ZYX')'; -%Starting State +% 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 +% 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 LOOP + for i=1:settings.stoch.N - %% Wind Generation - [uw,vw,ww] = windgen(settings.wind.AzMin,settings.wind.AzMax,... + + %% WIND GENERATION + + [uw,vw,ww] = wind_const_generator(settings.wind.AzMin,settings.wind.AzMax,... settings.wind.ElMin,settings.wind.ElMax,settings.wind.MagMin,... settings.wind.MagMax); - %% ASCEND %% + %% ASCEND [Ta,Ya] = ode45(@ascend,settings.ode.timeasc,X0a,settings.ode.optionsasc,... settings,uw,vw,ww); -%% DESCEND %% + %% DESCEND -[Td,Yd] = ode45(@ballistic_descent,settings.ode.timedesc,Ya(end,:),settings.ode.optionsdesc,... + [~,Yd] = ode45(@ballistic_descent,settings.ode.timedesc,Ya(end,:),settings.ode.optionsdesc,... settings,uw,vw,ww); - %% FINAL STATE ASSEMBLING %% + %% FINAL STATE ASSEMBLING %Total State Yf = [Ya;Yd]; - %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 +%% CHECKING BAD SIMULATION if numel(LP(LP(:,3)<-10,:)) fprintf(['Some landing points might be incorrect' ... 'please check parameters!\n']); end -% Writing Things +%% PRINTING VALUES -%Mean Landing Point +% Mean Landing Point xm = mean(LP(:,1)); ym = mean(LP(:,2)); -%Mean Apogee Time +% Mean Apogee Time ApoTimem = mean(ApoTime); -%Std. Deviation Apogee Time +% Std. Deviation Apogee Time ApoTimestd = std(ApoTime); -%Mean Altitude +% Mean Altitude Zm = mean(Z); -%Std. Deviation Altitude +% Std. Deviation Altitude Zstd = std(Z); -%Best Fit Ellipse +% 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',... +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 - +%% DEFAULT PLOTS + +if settings.default_plot == 1 plot(xm,ym,'bs','MarkerSize',20,'MarkerFacecolor','b'); hold on - %Point of launch + % Point of launch plot(0,0,'ro','MarkerSize',20,'MarkerFacecolor','r'); - %All the landing points + % All the landing points plot(LP(:,1),LP(:,2),'k+'); title('Landing Points'); @@ -115,7 +115,7 @@ if settings.plot == 1 view(90,270) axis equal - %Histogram + % Histogram [f,x] = hist(Z,10); figure; bar(x,f/settings.stoch.N); @@ -124,16 +124,13 @@ if settings.plot == 1 ylabel('n_i/n_{tot}') - - end -%Resizing +% 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 diff --git a/src/stoch_run_bal_p.m b/src/stoch_run_bal_p.m index f3206b74..93d6fef8 100644 --- a/src/stoch_run_bal_p.m +++ b/src/stoch_run_bal_p.m @@ -1,4 +1,4 @@ -function [ LP ,Z ] = stoch_run_bal_p( settings ) +function [LP ,Z] = stoch_run_bal_p(settings) %STD RUN - This function runs a stochastic simulation (parallel) % OUTPUTS % LP: Landing Points @@ -11,47 +11,52 @@ function [ LP ,Z ] = stoch_run_bal_p( settings ) % April 2014; Last revision: 29.V.2014 % License: 2-clause BSD -%Starting Attitude +%% STARTING CONDITIONS + +% Attitude Q0 = angle2quat(settings.PHI,settings.OMEGA,0*pi/180,'ZYX')'; -%Starting State +% 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 +% PreAllocation LP = zeros(settings.stoch.N,3); Z = zeros(settings.stoch.N,1); ApoTime = zeros(settings.stoch.N,1); -parfor_progress(settings.stoch.N); + +%% PARFOR LOOP + +parfor_progress(settings.stoch.N); % initiaize parfor loop parpool; + parfor i=1:settings.stoch.N - %% Wind Generation - [uw,vw,ww] = windgen(settings.wind.AzMin,settings.wind.AzMax,... + + %% WIND GENERATION + + [uw,vw,ww] = wind_const_generator(settings.wind.AzMin,settings.wind.AzMax,... settings.wind.ElMin,settings.wind.ElMax,settings.wind.MagMin,... settings.wind.MagMax); - %% ASCEND %% + %% ASCEND -[Ta,Ya] = ode45(@ascend,settings.ode.timeasc,X0a,settings.ode.optionsasc,... - settings,uw,vw,ww); + [Ta,Ya] = ode45(@ascend,settings.ode.timeasc,X0a,settings.ode.optionsasc,... + settings,uw,vw,ww); -%% DESCEND %% + + %% DESCEND -[Td,Yd] = ode45(@ballistic_descent,settings.ode.timedesc,Ya(end,:),settings.ode.optionsdesc,... + [~,Yd] = ode45(@ballistic_descent,settings.ode.timedesc,Ya(end,:),settings.ode.optionsdesc,... settings,uw,vw,ww); - %% FINAL STATE ASSEMBLING %% + %% FINAL STATE ASSEMBLING %Total State Yf = [Ya;Yd]; - - %Total Time - %Tf = [Ta; Ta(end)+Td;Ta(end)+Td(end)+Tm]; LP(i,:) = Yf(end,1:3); Z(i) = -Ya(end,3); @@ -61,28 +66,28 @@ parfor i=1:settings.stoch.N end -%Checking bad simulations +%% CHECKING BAD SIMULATION if numel(LP(LP(:,3)<-10,:)) fprintf(['Some landing points might be incorrect' ... 'please check parameters!\n']); end -% Writing Things +%% PRINTING VALUES -%Mean Landing Point +% Mean Landing Point xm = mean(LP(:,1)); ym = mean(LP(:,2)); -%Mean Apogee Time +% Mean Apogee Time ApoTimem = mean(ApoTime); -%Std. Deviation Apogee Time +% Std. Deviation Apogee Time ApoTimestd = std(ApoTime); -%Mean Altitude +% Mean Altitude Zm = mean(Z); -%Std. Deviation Altitude +% Std. Deviation Altitude Zstd = std(Z); % Printing to screen @@ -91,17 +96,17 @@ text =['Mean Landing Point:X:%3.3f m, Y:%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 - +%% DEFAULT PLOTS + +if settings.default_plot == 1 plot(xm,ym,'bs','MarkerSize',20,'MarkerFacecolor','b'); hold on - %Point of launch + % Point of launch plot(0,0,'ro','MarkerSize',20,'MarkerFacecolor','r'); - %All the landing points + % All the landing points plot(LP(:,1),LP(:,2),'k+'); title('Landing Points'); @@ -111,7 +116,7 @@ if settings.plot == 1 view(90,270) axis equal - %Histogram + % Histogram [f,x] = hist(Z,10); figure; bar(x,f/settings.stoch.N); @@ -120,20 +125,15 @@ if settings.plot == 1 ylabel('n_i/n_{tot}') - - end delete(gcp); -%Resizing +% Resizing h = get(0,'children'); scrsz = get(0,'ScreenSize'); for i=1:length(h) set(h(i),'OuterPosition',[0 0 scrsz(4) scrsz(4)]) - %saveas(h(i), ['figure' num2str(i)], 'fig'); end - - end diff --git a/src/stoch_run_p.m b/src/stoch_run_p.m index 0e786db0..4de17545 100644 --- a/src/stoch_run_p.m +++ b/src/stoch_run_p.m @@ -1,4 +1,4 @@ -function [ LP ,Z ] = stoch_run_p( settings ) +function [LP,Z] = stoch_run_p(settings) %STD RUN - This function runs a stochastic simulation (parallel) % OUTPUTS % LP: Landing Points @@ -11,71 +11,75 @@ function [ LP ,Z ] = stoch_run_p( settings ) % April 2014; Last revision: 29.V.2014 % License: 2-clause BSD -%Starting Attitude +%% STARTING CONDITIONS + +% Attitude Q0 = angle2quat(settings.PHI,settings.OMEGA,0*pi/180,'ZYX')'; -%Starting State +% 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 LOOP parfor_progress(settings.stoch.N); parpool; + parfor i=1:settings.stoch.N - %% Wind Generation - [uw,vw,ww] = windgen(settings.wind.AzMin,settings.wind.AzMax,... + + %% WIND GENERATION + + [uw,vw,ww] = wind_const_generator(settings.wind.AzMin,settings.wind.AzMax,... settings.wind.ElMin,settings.wind.ElMax,settings.wind.MagMin,... settings.wind.MagMax); - %% ASCEND %% + %% ASCEND [Ta,Ya] = ode45(@ascend,settings.ode.timeasc,X0a,settings.ode.optionsasc,... settings,uw,vw,ww); - %% DROGUE 1 %% - para = 1; %Flag for Drogue 1 + %% DROGUE 1 + % Initial Condition are the last from ascend (need to rotate because + % velocities are outputted in body axes) - %Initial Condition are the last from ascend (need to rotate because - %velocities are outputted in body axes) + para = 1; % Flag for Drogue 1 X0d1 = [Ya(end,1:3) quatrotate(quatconj(Ya(end,10:13)),Ya(end,4:6))]; - [Td1,Yd1] = ode45(@descent_parachute,settings.ode.timedrg1,X0d1,... + [~,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 + + para = 2; %Flag for Drogue 2 X0d2 = Yd1(end,:); - [Td2,Yd2] = ode45(@descent_parachute,settings.ode.timedrg2,X0d2,... + [~,Yd2] = ode45(@descent_parachute,settings.ode.timedrg2,X0d2,... settings.ode.optionsdrg2,settings,uw,vw,ww,para); - %% MAIN %% + %% MAIN - ROGALLO WING + % Initial Condition are the last from drogue 2 descent + if not(settings.ldf) - para = 3; %Flag for Main (Rogall) + para = 3; % Flag for Main (Rogall) end - %Initial Condition are the last from drogue 2 descent X0m = Yd2(end,:); - [Tm,Ym] = ode45(@descent_parachute,settings.ode.timemain,X0m,... + [~,Ym] = ode45(@descent_parachute,settings.ode.timemain,X0m,... settings.ode.optionsmain,settings,uw,vw,ww,para); - %% FINAL STATE ASSEMBLING %% + %% FINAL STATE ASSEMBLING - %Total State + % Total State if settings.ao Yf = [Ya(:,1:3) quatrotate(quatconj(Ya(:,10:13)),Ya(:,4:6))]; else Yf = [Ya(:,1:3) quatrotate(quatconj(Ya(:,10:13)),Ya(:,4:6));Yd1; Yd2;Ym]; end - %Total Time - %Tf = [Ta; Ta(end)+Td;Ta(end)+Td(end)+Tm]; LP(i,:) = Yf(end,1:3); Z(i) = -Ya(end,3); @@ -85,28 +89,29 @@ parfor i=1:settings.stoch.N end -%Checking bad simulations -if numel(LP(LP(:,3)<-10,:)) +%% CHECKING BAD SIMULATION + +if numel(LP(LP(:,3) < -10,:)) fprintf(['Some landing points might be incorrect' ... 'please check parameters!\n']); end -% Writing Things +%% PRINTING VALUES -%Mean Landing Point +% Mean Landing Point xm = mean(LP(:,1)); ym = mean(LP(:,2)); -%Mean Apogee Time +% Mean Apogee Time ApoTimem = mean(ApoTime); -%Std. Deviation Apogee Time +% Std. Deviation Apogee Time ApoTimestd = std(ApoTime); -%Mean Altitude +% Mean Altitude Zm = mean(Z); -%Std. Deviation Altitude +% Std. Deviation Altitude Zstd = std(Z); % Printing to screen @@ -115,17 +120,17 @@ text =['Mean Landing Point:X:%3.3f m, Y:%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 - +%% DEFAULT PLOTS + +if settings.default_plot == 1 plot(xm,ym,'bs','MarkerSize',20,'MarkerFacecolor','b'); hold on - %Point of launch + % Point of launch plot(0,0,'ro','MarkerSize',20,'MarkerFacecolor','r'); - %All the landing points + % All the landing points plot(LP(:,1),LP(:,2),'k+'); title('Landing Points'); @@ -135,7 +140,7 @@ if settings.plot == 1 view(90,270) axis equal - %Histogram + % Histogram [f,x] = hist(Z,10); figure; bar(x,f/settings.stoch.N); @@ -143,21 +148,15 @@ if settings.plot == 1 xlabel('Apogee [m]') ylabel('n_i/n_{tot}') - - - end delete(gcp); -%Resizing +% Resizing h = get(0,'children'); scrsz = get(0,'ScreenSize'); -for i=1:length(h) +for i = 1:length(h) set(h(i),'OuterPosition',[0 0 scrsz(4) scrsz(4)]) - %saveas(h(i), ['figure' num2str(i)], 'fig'); end - - end diff --git a/src/wind_const_generator.m b/src/wind_const_generator.m new file mode 100644 index 00000000..ebe9e7e7 --- /dev/null +++ b/src/wind_const_generator.m @@ -0,0 +1,35 @@ +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 +% +% 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 +R = Mag*angle2dcm(Az,El,0,'ZYX'); +uw = R(1,1); +vw = R(1,2); +ww = R(1,3); + + + + + + +end + diff --git a/src/wind_matlab_generator.m b/src/wind_matlab_generator.m new file mode 100644 index 00000000..f6878f4c --- /dev/null +++ b/src/wind_matlab_generator.m @@ -0,0 +1,40 @@ +function [uw,vw,ww] = wind_generator(settings,z,t,Q) +% wind_generator( settings,z,t,Q ) +% Function that generates wind components in body reference frame +% Based on hwm07 model +% nargin = 4:rotation of wind components in body axis +% nargin = 3: wind components in NED + + +% Author: Gabriele Poiana +% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu +% email: gabriele.poiana@skywarder.eu +% Website: http://www.skywarder.eu +% January 2016; Last revision: 17.I.2016 +% License: 2-clause BSD + + +h = -z+settings.z0; +if h < 0 + h = 0; +end + +%% HORIZONTAL WIND + +[uw,vw] = atmoshwm(settings.wind.Lat,settings.wind.Long,h,'day',settings.wind.Day,... + 'seconds',settings.wind.Seconds+t,'model','quiet','version','14'); %NED reference +ww = settings.wind.ww; + +% Rotation in body reference frame (ascend) +if nargin == 4 + + wind=quatrotate(Q,[uw vw ww]); % wind speed in body reference + + uw = wind(1); + vw = wind(2); + ww = wind(3); +end + +end + + diff --git a/src/wind_vert_generator.m b/src/wind_vert_generator.m new file mode 100644 index 00000000..61444b37 --- /dev/null +++ b/src/wind_vert_generator.m @@ -0,0 +1,20 @@ +function [ww] = vert_windgen(MagMin,MagMax) +% vert_windgen(MagMin,MagMax) +%function that generates vertical wind component + +% Author: Ruben Di Battista & Gabriele Poiana +% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu +% email: ruben.dibattista@skywarder.eu gabriele.poiana@skywarder.eu +% Website: http://www.skywarder.eu +% April 2014; Last revision: 17.I.2016 +% License: 2-clause BSD + +%Generating random value for magnitude + +x = rand; +ww = MagMin+(MagMax-MagMin)*x; +if x < 0.5 + ww = -ww; +end +end + diff --git a/src/MAIN.m b/src_old/MAIN.m similarity index 100% rename from src/MAIN.m rename to src_old/MAIN.m diff --git a/src/apogee.m b/src_old/apogee.m similarity index 100% rename from src/apogee.m rename to src_old/apogee.m diff --git a/src_old/ascend.m b/src_old/ascend.m new file mode 100644 index 00000000..bca5154e --- /dev/null +++ b/src_old/ascend.m @@ -0,0 +1,456 @@ +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.wind.model +% 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 +% +% % 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]); + +% 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 +global XCP +% 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); +XCPf =interp4_easy(givA,givM,givB,givH,CoeffsF.X_C_P,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); +XCPe =interp4_easy(givA,givM,givB,givH,CoeffsE.X_C_P,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; + XCP(contatore)= t/tb*(XCPe-XCPf)+XCPf; +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; + XCP(contatore) = XCPe; +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 diff --git a/src_old/ballistic_descent.m b/src_old/ballistic_descent.m new file mode 100644 index 00000000..a1fbb657 --- /dev/null +++ b/src_old/ballistic_descent.m @@ -0,0 +1,339 @@ +function [ dY ] = ballistic_descent( t,Y,settings,uw,vw,ww ) +% ODE-Function of the 6DOF Rigid Rocket Model +% State = ( x y z | u v w | p q r | q0 q1 q2 q3 ) +% +% (x y z): NED Earth's Surface Centered Frame ("Inertial") coordinates +% (u v w): body frame velocities +% (p q r): body frame angular rates +% (q0 q1 q2 q3): attitude unit quaternion +% +% +% NOTE: To get the NED velocities the body-frame must be multiplied for the +% conjugated of the current attitude quaternion +% E.G. +% +% +% quatrotate(quatconj(Y(:,10:13)),Y(:,4:6)) + +% Author: Ruben Di Battista +% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu +% email: ruben.dibattista@skywarder.eu +% Website: http://www.skywarder.eu +% April 2014; Last revision: 31.XII.2014 +% License: 2-clause BSD + +% Author: Francesco Colombi +% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu +% email: francesco.colombi@skywarder.eu +% Release date: 16/04/2016 + +x = Y(1); +y = Y(2); +z = Y(3); +u = Y(4); +v = Y(5); +w = Y(6); +p = Y(7); +q = Y(8); +r = Y(9); +q0 = Y(10); +q1 = Y(11); +q2 = Y(12); +q3 = Y(13); +m = Y(14); +Ixx = Y(15); +Iyy = Y(16); +Izz = Y(17); + + +Q = [ q0 q1 q2 q3]; +Q_conj= [ q0 -q1 -q2 -q3]; +normQ=norm(Q); + +if abs(normQ-1)>0.1 + Q = Q/normQ; + %frprintf('Norm of the quaternion is (too much!) less than one (%3.2f)\n', ... + %normQ); +end + + +% Adding Wind (supposed to be added in NED axes); + +% constant wind +wind = quatrotate(Q, [uw vw ww]); + +% % Wind model +% h = 0; +% if -z > 0 +% h = -z+settings.z0; +% end +% % global alt +% % if bool == 1 +% % alt = [alt; h]; +% % end +% +% % wind in NED axis +% [uw,vw] = atmoshwm07(settings.wind.Lat, settings.wind.Long, h, ... +% 'day',settings.wind.Day,'model','quiet'); +% +% % global WIND +% % if bool == 1 +% % WIND = [WIND; uw, vw, ww]; +% % end +% +% % wind in body frame +% wind = quatrotate(Q, [uw, vw, ww]); + + +% % Wind Model +% [uw,vw,ww] = wind_generator(settings,z,t,Q); +% wind = [ uw,vw,ww ]; + +% Relative velocities (plus wind); +ur = u - wind(1); +vr = v - wind(2); +wr = w - wind(3); + +% Body to Inertial velocities +Vels = quatrotate(Q_conj,[u v w]); + +V_norm = norm([ur vr wr]); + + + +% Constants +S = settings.S; +C = settings.C; +CoeffsE = settings.CoeffsE; %Empty Rocket Coefficients + +g = 9.81; +%Time-Dependant Properties +m0 = settings.m0; %kg mass +mfr = settings.mfr; +tb = settings.tb; + +Ixxf= settings.Ixxf; +Ixxe= settings.Ixxe; +Iyyf= settings.Iyyf; +Iyye= settings.Iyye; +Izzf= settings.Izzf; +Izze= settings.Izze; + +dI = 1/tb*([Ixxf Iyyf Izzf]'-[Ixxe Iyye Izze]'); + + mdot = 0; + Ixxdot = 0; + Iyydot = 0; + Izzdot = 0; + T=0; + +%Atmosphere +if -z<0 + z = 0; +end +[~, a, ~, rho] = atmoscoesa(-z+settings.z0); +M = V_norm/a; + + +%Aero Angles +if not(u<1e-1 || V_norm<1e-3); + alpha = atan(wr/ur); + beta = asin(vr/V_norm); +else + alpha =0; + beta = 0; +end + +% beta = 0; % prova simulazione con beta imposto uguale a 0 + + + +% Coeffs. Interpolation +givA = settings.Alphas*pi/180; +givB = settings.Betas*pi/180; +givH = settings.Altitudes; +givM = settings.Machs; + +% Interpolation error check + +if M > givM(end) + %frprintf('Mach No. is more than the max provided (M = %3.2f @ t=%2.2f)\n\n', ... + %M,t); + M = givM(end); + +end + +if M<givM(1) + %frprintf('Mach No. is less than the min provided (M = %3.2f @ t=%2.2f)\n\n', ... + % M,t); + M = givM(1); + +end + +if alpha > givA(end) + %frprintf('AoA is more than the max provided (alpha = %3.2f @ t=%2.2f)\n\n', ... + %alpha*180/pi,t) + alpha= givA(end); + +else if alpha < givA(1) + %frprintf('AoA is less than the min provided (alpha = %3.2f @ t=%2.2f)\n\n',... + %alpha*180/pi,t); + alpha = givA(1); + + + end +end + +if beta > givB(end) + beta= givB(end); +else if beta < givB(1) + beta = givB(1); + end +end + +if -z >givH(end) + z = -givH(end); +else if -z < givH(1) + z = -givH(1); + end +end + +% interpolation with the value in the nearest condition: interp_easy +% empty +CAe=interp4_easy(givA,givM,givB,givH,CoeffsE.CA,alpha,M,beta,-z);%,'nearest'); +CYBe=interp4_easy(givA,givM,givB,givH,CoeffsE.CYB,alpha,M,beta,-z);%,'nearest'); +CNAe=interp4_easy(givA,givM,givB,givH,CoeffsE.CNA,alpha,M,beta,-z);%,'nearest'); +Cle=interp4_easy(givA,givM,givB,givH,CoeffsE.CLL,alpha,M,beta,-z);%,'nearest'); +Clpe=interp4_easy(givA,givM,givB,givH,CoeffsE.CLLP,alpha,M,beta,-z);%,'nearest'); +Cmae=interp4_easy(givA,givM,givB,givH,CoeffsE.CMA,alpha,M,beta,-z);%,'nearest'); +Cmade=interp4_easy(givA,givM,givB,givH,CoeffsE.CMAD,alpha,M,beta,-z);%,'nearest'); +Cmqe=interp4_easy(givA,givM,givB,givH,CoeffsE.CMQ,alpha,M,beta,-z);%,'nearest'); +Cnbe=interp4_easy(givA,givM,givB,givH,CoeffsE.CLNB,alpha,M,beta,-z);%,'nearest'); +Cnre=interp4_easy(givA,givM,givB,givH,CoeffsE.CLNR,alpha,M,beta,-z);%,'nearest'); +Cnpe=interp4_easy(givA,givM,givB,givH,CoeffsE.CLNP,alpha,M,beta,-z);%,'nearest'); +% +Clre=interp4_easy(givA,givM,givB,givH,CoeffsE.CLLR,alpha,M,beta,-z);%,'nearest'); +Clbe=interp4_easy(givA,givM,givB,givH,CoeffsE.CLLB,alpha,M,beta,-z);%,'nearest'); +CDe =interp4_easy(givA,givM,givB,givH,CoeffsE.CD,alpha,M,beta,-z); + +% % linear interpolation of the coeff +% %full +% CAf = interpn(givA,givM,givB,givH,CoeffsF.CA,alpha,M,beta,-z); +% CYBf = interpn(givA,givM,givB,givH,CoeffsF.CYB,alpha,M,beta,-z); +% CNAf = interpn(givA,givM,givB,givH,CoeffsF.CNA,alpha,M,beta,-z); +% Clf = interpn(givA,givM,givB,givH,CoeffsF.CLL,alpha,M,beta,-z); +% Clpf = interpn(givA,givM,givB,givH,CoeffsF.CLLP,alpha,M,beta,-z); +% Cmaf = interpn(givA,givM,givB,givH,CoeffsF.CMA,alpha,M,beta,-z); +% Cmadf = interpn(givA,givM,givB,givH,CoeffsF.CMAD,alpha,M,beta,-z); +% Cmqf = interpn(givA,givM,givB,givH,CoeffsF.CMQ,alpha,M,beta,-z); +% Cnbf = interpn(givA,givM,givB,givH,CoeffsF.CLNB,alpha,M,beta,-z); +% Cnrf = interpn(givA,givM,givB,givH,CoeffsF.CLNR,alpha,M,beta,-z); +% Cnpf = interpn(givA,givM,givB,givH,CoeffsF.CLNP,alpha,M,beta,-z); +% % +% Clrf = interpn(givA,givM,givB,givH,CoeffsF.CLLR,alpha,M,beta,-z); +% Clbf = interpn(givA,givM,givB,givH,CoeffsF.CLLB,alpha,M,beta,-z); +% CDf = interpn(givA,givM,givB,givH,CoeffsF.CD,alpha,M,beta,-z); +% +% % empty +% CAe = interpn(givA,givM,givB,givH,CoeffsE.CA,alpha,M,beta,-z); +% CYBe = interpn(givA,givM,givB,givH,CoeffsE.CYB,alpha,M,beta,-z); +% CNAe = interpn(givA,givM,givB,givH,CoeffsE.CNA,alpha,M,beta,-z); +% Cle = interpn(givA,givM,givB,givH,CoeffsE.CLL,alpha,M,beta,-z); +% Clpe = interpn(givA,givM,givB,givH,CoeffsE.CLLP,alpha,M,beta,-z); +% Cmae = interpn(givA,givM,givB,givH,CoeffsE.CMA,alpha,M,beta,-z); +% Cmade = interpn(givA,givM,givB,givH,CoeffsE.CMAD,alpha,M,beta,-z); +% Cmqe = interpn(givA,givM,givB,givH,CoeffsE.CMQ,alpha,M,beta,-z); +% Cnbe = interpn(givA,givM,givB,givH,CoeffsE.CLNB,alpha,M,beta,-z); +% Cnre = interpn(givA,givM,givB,givH,CoeffsE.CLNR,alpha,M,beta,-z); +% Cnpe = interpn(givA,givM,givB,givH,CoeffsE.CLNP,alpha,M,beta,-z); +% % +% Clre = interpn(givA,givM,givB,givH,CoeffsE.CLLR,alpha,M,beta,-z); +% Clbe = interpn(givA,givM,givB,givH,CoeffsE.CLLB,alpha,M,beta,-z); +% CDe =interpn(givA,givM,givB,givH,CoeffsE.CD,alpha,M,beta,-z); + +%Linear interpolation from empty and full configuration coefficients + CA = CAe; + CYB= CYBe; + CNA= CNAe; + Cl = Cle; + Clp= Clpe; + Cma= Cmae; + Cmad=Cmade; + Cmq= Cmqe; + Cnb= Cnbe; + Cnr= Cnre; + Cnp= Cnpe; + % + Clr= Clre; + Clb= Clbe; + CD = CDe; + + +%Center of Mass + + + +%Body Frame + +qdynS=0.5*rho*V_norm^2*S; +qdynL_V = 0.5*rho*V_norm*S*C; + +X=qdynS*CA; +Y=qdynS*CYB*beta; +Z=qdynS*CNA*alpha; + + + +%Forces +F = quatrotate(Q,[0 0 m*g])'; + +F = F +[-X+T,+Y,-Z]'; + + +du=F(1)/m -q*w + r*v;% - mdot/m*u; +dv=F(2)/m -r*u + p*w;% - mdot/m*v; +dw=F(3)/m -p*v + q*u;% - mdot/m*w; + + + + +% Rotation +dp=(Iyy-Izz)/Ixx*q*r + qdynL_V/Ixx*(V_norm*Cl+Clp*p*C/2)-Ixxdot*p/Ixx; +dq=(Izz-Ixx)/Iyy*p*r + qdynL_V/Iyy*(V_norm*Cma*alpha + (Cmad+Cmq)*q*C/2)... + -Iyydot*q/Iyy; +dr=(Ixx-Iyy)/Izz*p*q + qdynL_V/Izz*(V_norm*Cnb*beta + (Cnr*r+Cnp*p)*C/2)... + -Izzdot*r/Izz; + + +% dp=(Iyy-Izz)/Ixx*q*r + qdynL_V/Ixx*(V_norm*(Cl+Clb*beta)+(Clp*p+Clr*r)*C/2)... +% -Ixxdot*p/Ixx; +% dq=(Izz-Ixx)/Iyy*p*r + qdynL_V/Iyy*(V_norm*Cma*alpha + (Cmad+Cmq)*q*C/2)... +% -Iyydot*q/Iyy; +% dr=(Ixx-Iyy)/Izz*p*q + qdynL_V/Izz*(V_norm*Cnb*beta + (Cnr*r+Cnp*p)*C/2)... +% -Izzdot*r/Izz; + + +OM = 1/2*[ + 0 -p -q -r + p 0 r -q + q -r 0 p + r q -p 0]; + + +dQQ = OM*Q'; + +dY(1:3) = Vels; +dY(4) = du; +dY(5) = dv; +dY(6) = dw; +dY(7) = dp; +dY(8) = dq; +dY(9) = dr; +dY(10:13) = dQQ; +dY(14) = mdot; +dY(15) = Ixxdot; +dY(16) = Iyydot; +dY(17) = Izzdot; +dY=dY'; + +end \ No newline at end of file diff --git a/src_old/config.m b/src_old/config.m new file mode 100644 index 00000000..072926f4 --- /dev/null +++ b/src_old/config.m @@ -0,0 +1,262 @@ +%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 + +%% LAUNCH SETUP + +% rocket name +settings.rocket_name = 'R2A'; + +% launchpad 6 +settings.z0 = 5; %[m] Launchpad Altitude +settings.lrampa = 6.5; %[m] LaunchPad route (launchpad length-distance from ground of the first hook) + + +% starting altitude +settings.OMEGA = 80*pi/180; %[rad] Elevation Angle, user input in degrees (ex. 80) +settings.PHI = 90*pi/180; %[rad] Azimuth Angle from North Direction, user input in degrees (ex. 90) + +%% ENGINE DETAILS + +% % Thrust on time profile (IV Order Interpolation) +% settings.T_coeffs = [8034.5, 0, 0, 0, 0]; % constant thrust +% data from the motor details sheet (experimental plot of the Thrust) + +% SYNTAX: +% engine = 1 -> Cesaroni PRO 150 White Thunder +% engine = 2 -> Cesaroni PRO 150 SkidMark +% engine = 3 -> Cesaroni PRO 150 BlueStreak +engine = 3; +switch engine + case 1 + % Cesaroni PRO 150 White Thunder + % Sampling for thrust interpolation + settings.motor.exp_time = [0 0.05 0.15 0.5 0.6 0.74 0.85 1.15 1.7 2.4 3 ... + 4 4.5 4.8 4.9 5 5.05 5.1 5.15 5.2]; %s + settings.motor.exp_thrust = [8605.1 8900 7900 8400 8400 8250 8200 8300 ... + 8400 8400 8200 7800 7600 7450 7350 7300 4500 500 100 0]; %N + + settings.m0 = 67.761; %kg Overall Mass (Burnout + Propellant) + settings.ms = 43.961; %kg Structural Mass (Burnout - Nosecone) + settings.mp = 18.6; %kg Propellant Mass + settings.tb = 5.12; %sec Burning Time + settings.mfr = settings.mp/settings.tb; %kg/s Mass Flow Rate + case 2 + % Cesaroni PRO 150 SkidMark + % Sampling for thrust interpolation + settings.motor.exp_time = [0 0.1 0.2 0.3 0.4 0.5 0.6 0.8 1.2 1.8 3.2 ... + 3.6 4.8 6 7 7.2 7.6 7.8 7.9 8 8.1 8.19]; %s + settings.motor.exp_thrust = [0 3400 3100 3000 3300 3400 3500 3700 3700 ... + 3800 4000 4081.6 3900 3800 3700 3500 3350 3200 3000 2000 750 0]; %N + + settings.m0 = 64.9; %kg Overall Mass + settings.ms = 46.8; %kg Structural Mass (Burnout) + settings.mp = settings.m0-settings.ms; %kg Propellant Mass + settings.tb = 8.19; %sec Burning Time + settings.mfr = settings.mp/settings.tb; %kg/s Mass Flow Rate + case 3 + % Cesaroni PRO 150 BlueStreak + % Sampling for thrust interpolation + settings.motor.exp_time = [0 0.06 0.1 0.15 0.25 0.45 0.8 1 2 3 ... + 4 5 6 6.8 7.05 7.3 7.6 7.8]; %s + settings.motor.exp_thrust = [0 800 4000 5500 5160 5130 5400 5300 5450 5347 ... + 5160 4950 4700 4400 4400 3800 300 0]; %N + + settings.m0 = 66.2; %kg Overall Mass + settings.ms = 47.3; %kg Structural Mass (Burnout) + settings.mp = settings.m0-settings.ms; %kg Propellant Mass + settings.tb = 7.60; %sec Burning Time + settings.mfr = settings.mp/settings.tb; %kg/s Mass Flow Rate + otherwise +end + + +% GEOMETRY DETAILS % +% This parameters should be the same parameters set up in MISSILE DATCOM +% simulation. + +settings.C = 0.174; %m Caliber (Fuselage Diameter) +settings.S = 0.02378; %m2 Cross-sectional Surface +L = 4.4; %m Rocket length + +% MASS GEOMERTY DETAILS % +% x-axis: along the fuselage +% y-axis: right wing +% z-axis: downward + +% % Note: inerzias are used in "apogee_reached.m" +% % HP: rocket inertias = full finite cilinder inertias +% settings.Ixxf=settings.m0*(settings.C/2)^2/2; %Inertia to x-axis (Full) +% settings.Ixxe=settings.ms*(settings.C/2)^2/2; %Inertia to x-axis (Empty) +% settings.Iyyf=settings.m0.*((settings.C/2)^2/4 + L^2/3); %Inertia to y-axis (Full) +% settings.Iyye=settings.ms.*((settings.C/2)^2/4 + L^2/3); %Inertia to y-axis (Empty) +% settings.Izzf=settings.Iyyf; %Inertia to z-axis (Full) +% settings.Izze=settings.Iyye; %Inertia to z-axis (Empty) + + +% inertias first approximation +settings.Ixxf = 0.27; %kg*m2 Inertia to x-axis (Full) +settings.Ixxe = 0.21; %kg*m2 Inertia to x-axis (Empty) +settings.Iyyf = 86.02; %kg*m2 Inertia to y-axis (Full) +settings.Iyye = 66.84; %kg*m2 Inertia to y-axis (Empty) +settings.Izzf = 86.02; %kg*m2 Inertia to z-axis (Full) +settings.Izze = 66.84; %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) +% The files are stored in the ../data folder with a naming convention of +% rocket_name_full.mat | rocket_name_empty.mat +% e.g. R1X_full.mat etc.. + +%Relative Path of the data files (default: ../data/). Remember the trailing +% slash!! + +DATA_PATH = '../data/'; +filename = strcat(DATA_PATH, settings.rocket_name); + +%Coefficients in full configuration (with all the propellant embarqued) +filename_full = strcat(filename,'_full.mat'); +CoeffsF = load(filename_full,'Coeffs'); +settings.CoeffsF = CoeffsF.Coeffs; +clear('CoeffsF'); + +%Coefficients in empty configuration (all the propellant consumed) +filename_empty = strcat(filename,'_empty.mat'); +CoeffsE = load(filename_empty,'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(filename_full,'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.55; %m2 Surface +settings.para1.mass = 0.577; %kg Parachute Mass +settings.para1.CD = 0.8; %Parachute Drag Coefficient +settings.para1.CL = 0; %Parachute Lift Coefficient + +%Altitude of Drogue 2 Opening +settings.zdrg2 = 5000; + +%%% DROGUE 2 %%% +settings.para2.S = 17.5; %m2 Surface +settings.para2.mass = 1.140; %kg Parachute Mass +settings.para2.CD = 0.59; %Parachute Drag Coefficient +settings.para2.CL = 0; %Parachute Lift Coefficient + +%Altitude of Main Parachute Opening +settings.zmain = 2000; + +%%% MAIN - ROGALLO %%% +%The drogue parachute effects are neglected + +settings.para3.S = 15; %m2 Surface +settings.para3.mass = 1.466; %kg Parachute Mass +settings.para3.CD = 0.4; %Parachute Drag Coeff +settings.para3.CL = 0.8; %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 calculator model +settings.wind.model = false; +%true for hwsm wind model +%false for constant wind model + +% Wind is randomly generated. Setting the same values for min and max will +% fix the parameters of the wind. + +settings.wind.MagMin = 3; %Minimum Magnitude +settings.wind.MagMax = 3; %Maximum Magnitude +settings.wind.ElMin = 0*pi/180; %Minimum Elevation +settings.wind.ElMax = 0*pi/180; %Maximum Elevation (Max == 90 Deg) +settings.wind.AzMin = (90)*pi/180; %Minimum Azimuth +settings.wind.AzMax = (90)*pi/180; %Maximum Azimuth + +% NOTE: wind aziumt angle indications (wind directed towards): +% 0 deg (use 360 instead of 0) -> North +% 90 deg -> East +% 180 deg -> South +% 270 deg -> West + +% Settings for the Wind Model +settings.wind.Lat = 39.552709; %Latitude of launching site +settings.wind.Long = 9.652400; %Longitude of launching site +settings.wind.Day = 290; %Day of the launch +settings.wind.Seconds = 13*60*60; %Second of the day (UTM) + +% 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 = true; % Set to True to run a standard ballistic simulation + +% LAST DROGUE FAILURE SIMULATION +% simulation in which rogallo wing does not open, therefore landing is +% achieved thanks to the 2nd parachute +settings.ldf = false; + +% APOGEE ONLY +% simulation stopped when reaching the apogee, therefore there is no +% descend phase. Only available for stochastic runs +settings.ao = false; + +% STOCHASTIC DETAILS % +%If N>1 the stochastic routine is fired (different standard plots) +settings.stoch.N = 1; % Number of iterations +settings.stoch.parallel = false; % Using parallel or not parallel + +% PLOT DETAILS % +settings.plot = true; % 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 diff --git a/src/crash.m b/src_old/crash.m similarity index 100% rename from src/crash.m rename to src_old/crash.m diff --git a/src_old/descent_parachute.m b/src_old/descent_parachute.m new file mode 100644 index 00000000..62c0c267 --- /dev/null +++ b/src_old/descent_parachute.m @@ -0,0 +1,158 @@ +function [ dY ] = descent_parachute( t,Y,settings,uw,vw,ww,para ) +% ODEFUN for Parachute descent +% State = ( x y z | u v w ) + +% Author: Ruben Di Battista +% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu +% email: ruben.dibattista@skywarder.eu +% Website: http://www.skywarder.eu +% April 2014; Last revision: 25.IV.2014 +% License: 2-clause BSD + +% Author: Francesco Colombi +% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu +% email: francesco.colombi@skywarder.eu +% Release date: 16/04/2016 + +x = Y(1); +y = Y(2); +z = Y(3); +u = Y(4); +v = Y(5); +w = Y(6); + + +% global bool +% % if bool = 0 the trends during the integration are NOT requested +% % if bool = 1 the trends during the integration are saved and plotted + + +% Adding Wind (supposed to be added in NED axes); +if settings.wind.model +% 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 +% +% % wind in NED axis +% [uw,vw] = atmoshwm07(settings.wind.Lat, settings.wind.Long, h, ... +% 'day',settings.wind.Day,'model','quiet'); +% wind = [uw, vw, ww]; +% +% % global WIND +% % if bool == 1 +% % WIND = [WIND; uw, vw, ww]; +% % end + +% Adding wind; +ur = u - wind(1); +vr = v - wind(2); +wr = w - wind(3); + +V_norm = norm([ur vr wr]); + +% Constants +switch para + case 1 + S = settings.para1.S; %Parachute Surface + CD = settings.para1.CD; %Parachute CD + CL = settings.para1.CL; %Parachute CL + pmass = 0; %detached parachute mass + case 2 + S = settings.para2.S; + CD = settings.para2.CD; + CL = settings.para2.CL; + pmass = settings.para1.mass + 6.13; %parachute + nosecone + case 3 + S = settings.para3.S; + CD = settings.para3.CD; + CL = settings.para3.CL; + pmass = settings.para1.mass + settings.para2.mass + 6.13; + otherwise +end + +rad = sqrt(S/pi); %Equivalent parachute radius; + +%Side Surface +Ss = pi*(rad/2)^2; + +g = 9.81; + +m = settings.ms - pmass; + +%Atmosphere +if -z<0 + z = 0; +end + +[~, ~, ~, rho] = atmoscoesa(-z); + +%Center of Mass + + + +%Body Frame +D = 0.5*rho*V_norm^2*S*CD; +L = 0.5*rho*V_norm^2*S*CL; + +% HP: il paracadute � approssimato come una superficie rettangolare +% orientata in modo da essere sempre normale alla direzione della velocita +% relativa all'aria e sempre con un asse orizzontale e la forza Lift verso +% l'alto +% in modo da puntare ad atterrare nel punto [x,y] = (0m a nord,500m ad +% ovest) rispetto la postazione di lancio + + +tt = [ur vr wr]; % versore tangenziale +kk = [-vr ur 0]; % versore orizzontale +if (kk(1) == 0) && (kk(2) == 0) + kk = [-v u 0]; +end +tt = tt/norm(tt); +kk = -kk/norm(kk); +nn = cross(tt, kk); nn = nn/norm(nn); % versore normale +if (nn(3) > 0) % nn diretto verso il basso ? + nn = cross(kk, tt); nn = nn/norm(nn); +end + + +%Forces + +% X = D*ur/V_norm; +% Y = D*vr/V_norm; +% Z = D*wr/V_norm -m*g; + +% F = [-X,-Y,-Z]'; +F = -D*tt' + L*nn' + m*g*[0 0 1]'; + + +du=F(1)/m; +dv=F(2)/m; +dw=F(3)/m; + + + + + +dY(1:3) = [u v w]'; +dY(4) = du; +dY(5) = dv; +dY(6) = dw; + +dY = dY'; + +end \ No newline at end of file diff --git a/src/drg2_opening.m b/src_old/drg2_opening.m similarity index 100% rename from src/drg2_opening.m rename to src_old/drg2_opening.m diff --git a/src_old/interp4_easy.m b/src_old/interp4_easy.m new file mode 100644 index 00000000..ba3d1845 --- /dev/null +++ b/src_old/interp4_easy.m @@ -0,0 +1,40 @@ +function [ V ] = interp4_easy( X1,X2,X3,X4,TABLE,x1,x2,x3,x4 ) +%windgen( AzMin,AzMax,ElMin,ElMax,MagMin,MagMax ) +% This function interpolate with nearest-neighbor method a R4->R function +% F(x1...x4), given a [NxMxLxI] Matrix that discretize the function +% +% +% +% (X1...X4): column vectors of the variables the function depends on +% (x1....x4): single values for each variable the interpolation is wanted +% on +% TABLE: the [NxMxLxI] matrix +% +% +% length(X1) = N +% length(X2) = M +% length(X3) = L +% length(X4) = I + +% Author: Ruben Di Battista +% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu +% email: ruben.dibattista@skywarder.eu +% Website: http://www.skywarder.eu +% April 2014; Last revision: 25.IV.2014 +% License: 2-clause BSD + +M = {X1 X2 X3 X4}; +m = [ x1 x2 x3 x4]; + +index = zeros(4,1); + +for i=1:4 + [~,index(i)] = min(abs(M{i} - m(i))); +end + +V =TABLE(index(1),index(2),index(3),index(4)); + + + +end + diff --git a/src/launch_pad.m b/src_old/launch_pad.m similarity index 100% rename from src/launch_pad.m rename to src_old/launch_pad.m diff --git a/src/main_opening.m b/src_old/main_opening.m similarity index 100% rename from src/main_opening.m rename to src_old/main_opening.m diff --git a/src_old/parfor_progress.m b/src_old/parfor_progress.m new file mode 100644 index 00000000..71bd7678 --- /dev/null +++ b/src_old/parfor_progress.m @@ -0,0 +1,81 @@ +function percent = parfor_progress(N) +%PARFOR_PROGRESS Progress monitor (progress bar) that works with parfor. +% PARFOR_PROGRESS works by creating a file called parfor_progress.txt in +% your working directory, and then keeping track of the parfor loop's +% progress within that file. This workaround is necessary because parfor +% workers cannot communicate with one another so there is no simple way +% to know which iterations have finished and which haven't. +% +% PARFOR_PROGRESS(N) initializes the progress monitor for a set of N +% upcoming calculations. +% +% PARFOR_PROGRESS updates the progress inside your parfor loop and +% displays an updated progress bar. +% +% PARFOR_PROGRESS(0) deletes parfor_progress.txt and finalizes progress +% bar. +% +% To suppress output from any of these functions, just ask for a return +% variable from the function calls, like PERCENT = PARFOR_PROGRESS which +% returns the percentage of completion. +% +% Example: +% +% N = 100; +% parfor_progress(N); +% parfor i=1:N +% pause(rand); % Replace with real code +% parfor_progress; +% end +% parfor_progress(0); +% +% See also PARFOR. + +% By Jeremy Scheff - jdscheff@gmail.com - http://www.jeremyscheff.com/ + +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/src_old/parfor_progress.txt b/src_old/parfor_progress.txt new file mode 100644 index 00000000..afe234be --- /dev/null +++ b/src_old/parfor_progress.txt @@ -0,0 +1,3 @@ +2 +1 +1 diff --git a/src_old/plot_data.m b/src_old/plot_data.m new file mode 100644 index 00000000..2ed42bf7 --- /dev/null +++ b/src_old/plot_data.m @@ -0,0 +1,156 @@ +% Author: Francesco Colombi +% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu +% email: francesco.colombi@skywarder.eu +% Release date: 16/04/2016 + +if plots + 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; + p1 = plot(0,0,'r.','MarkerSize',21); + p2 = plot(y(end),x(end),'rx','MarkerSize',14); + title('Trajectory Top view') + xlabel('y, East [m]'), ylabel('x, North [m]'); + legend([p1,p2],{'Launch Point','Landing Point'}) + clear p1 p2 + + % 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; + p11 = plot(Ya(end,2)/1000, Ya(end,1)/1000, 'r*','markersize',7); + p12 = plot(0, 0, 'r.','markersize',14); + p13 = plot(Y(end,2)/1000, Y(end,1)/1000, 'rx','markersize',7); + 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, 'r*','markersize',7); + plot(Y(end,1)/1000, -Y(end,3)/1000, 'rx','markersize',7); + plot(0, 0, 'r.','markersize',14); + 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, 'r*','markersize',7); + plot(Y(end,2)/1000, -Y(end,3)/1000, 'rx','markersize',7); + plot(0, 0, 'r.','markersize',14); + xlabel('y, East [Km]'), ylabel('z, Altitude [Km]'); + %setting limit if is parallel to north + if sum(y_flight) == 0 + xlim([-1 1]); + end + legend([p11, p12, p13],{'Apogee','Launch Point','Landing Point'}) + + % angular rates + figure(); + plot(Ta, Ya(:, 7)*180/pi) + hold on, grid on + plot(Ta, Ya(:, 8)*180/pi) + plot(Ta, Ya(:, 9)*180/pi) + xlabel('time [sec]'), ylabel('p, q, r [grad/sec]') + legend('p: roll rate', 'q: pitch rate', 'r: yaw rate') + title('angular rates vs time') + + figure(); + subplot(3,1,1) + plot(Ta, Ya(:, 8)*180/pi), grid on; + xlabel('time [sec]'), ylabel('pitch rate q [grad/sec]') + subplot(3,1,2) + plot(Ta, Ya(:, 9)*180/pi), grid on; + xlabel('time [sec]'), ylabel('yaw rate r [grad/sec]') + subplot(3,1,3) + plot(Ta, Ya(:, 7)*180/pi), grid on; + xlabel('time [sec]'), ylabel('roll rate p [grad/sec]') + + + % angle + alpha = zeros(length(Ta),1); + beta = zeros(length(Ta),1); + roll = zeros(length(Ta),1); + + for k = 2:length(Ta) + alpha(k) = alpha(k-1) + (Ya(k, 8) + Ya(k-1, 8))/2*180/pi*(T(k)-T(k-1)); + beta(k) = beta(k-1) + (Ya(k, 9) + Ya(k-1, 9))/2*180/pi*(T(k)-T(k-1)); + roll(k) = roll(k-1) + (Ya(k, 7) + Ya(k-1, 7))/2*180/pi*(T(k)-T(k-1)); + end + figure(), plot(Ta, alpha+settings.OMEGA*180/pi), title('Alpha vs time'), grid on; + figure(), plot(Ta, beta), title('Beta vs time'), grid on; + figure(), plot(Ta, roll), title('Roll vs time'), grid on; + + + % Mach + figure(), plot(T, M), grid on; + xlabel('time [sec]'), ylabel('Mach [-]'); + + % Stagnation Temperature + figure(); + plot(T, Tamb-273.15); + hold on, grid on; + plot(T, Ttot-273.15) + title('Temperature profile') + xlabel('time [sec]'), ylabel('Temperature [???C]'); + legend('Surrounding', 'Total', 'location', 'best') + +end \ No newline at end of file diff --git a/src/print_license.m b/src_old/print_license.m similarity index 100% rename from src/print_license.m rename to src_old/print_license.m diff --git a/src_old/start_simulation.m b/src_old/start_simulation.m new file mode 100644 index 00000000..215e7a36 --- /dev/null +++ b/src_old/start_simulation.m @@ -0,0 +1,356 @@ +% Author: Francesco Colombi +% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu +% email: francesco.colombi@skywarder.eu +% Release date: 16/04/2016 + +close all +clear, clc + +%% caricamento dati +run('config.m'); + +global bool contatore; +global t_plot T_plot alpha_plot beta_plot M_plot CA_plot Drag_plot Forces_plot; +% global WIND alt; +global XCP; + + +% 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 = []; + XCP = 0; +end + +%% lancio simulatore +[T,Y, Ta,Ya] = MAIN(settings); + +if (bool == 1) + figure(); + plot(t_plot, T_plot, '.'), title('Thrust vs time'), grid on; + ylabel('Thrust [N]') + figure(); + plot(t_plot, alpha_plot*180/pi), title('alpha vs time'), grid on; + ylabel('alpha [deg]') + figure(); + plot(t_plot, beta_plot*180/pi), title('beta vs time'), grid on; + ylabel('beta [deg]') + + figure(); + plot(t_plot, M_plot), title('mach vs time'), grid on; + ylabel('Mach M [-]') + + figure(); + plot(t_plot, CA_plot), title('Aerodyn Coeff vs time'), grid on; + ylabel('Drag Coeff CD [-]') + + figure(); + plot(t_plot, Drag_plot), title('Drag vs time'), grid on; + ylabel('Drag D [N]') + + figure(); + plot(t_plot, Forces_plot), title('Axial Force vs time'), grid on; + ylabel('Axial force [N]') + +% figure() +% title('Wind'), grid on, hold on; +% plot(WIND(:, 1), alt) +% plot(WIND(:, 2), alt) +% plot(WIND(:, 3), alt) +% xlabel('Wind magnitude [m/s]') +% ylabel('Altitude [m]') +% legend('North','East','Down') + + figure(); + plot(t_plot, XCP,'.'), title('Stability margin vs time'), grid on; + ylabel('s.m.') + +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 + p1 = plot(0,0,'r.','MarkerSize',21); + p2 = plot(y(end),x(end),'rx','MarkerSize',14); + title('Trajectory Top view') + xlabel('y, East [m]'), ylabel('x, North [m]'); + legend([p1,p2],{'Launch Point','Landing Point'}) + clear p1 p2 + + % 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; + p11 = plot(Ya(end,2)/1000, Ya(end,1)/1000, 'r*','markersize',7); + p12 = plot(0, 0, 'r.','markersize',14); + p13 = plot(Y(end,2)/1000, Y(end,1)/1000, 'rx','markersize',7); + 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, 'r*','markersize',7); + plot(Y(end,1)/1000, -Y(end,3)/1000, 'rx','markersize',7); + plot(0, 0, 'r.','markersize',14); + 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, 'r*','markersize',7); + plot(Y(end,2)/1000, -Y(end,3)/1000, 'rx','markersize',7); + plot(0, 0, 'r.','markersize',14); + xlabel('y, East [Km]'), ylabel('z, Altitude [Km]'); + %setting limit if is parallel to north + if sum(y_flight) == 0 + xlim([-1 1]); + end + legend([p11, p12, p13],{'Apogee','Launch Point','Landing Point'}) + + % angular rates + figure(); + plot(Ta, Ya(:, 7)*180/pi) + hold on, grid on + plot(Ta, Ya(:, 8)*180/pi) + plot(Ta, Ya(:, 9)*180/pi) + xlabel('time [sec]'), ylabel('p, q, r [grad/sec]') + legend('p: roll rate', 'q: pitch rate', 'r: yaw rate') + title('angular rates vs time') + + figure(); + subplot(3,1,1) + plot(Ta, Ya(:, 8)*180/pi), grid on; + xlabel('time [sec]'), ylabel('pitch rate q [grad/sec]') + subplot(3,1,2) + plot(Ta, Ya(:, 9)*180/pi), grid on; + xlabel('time [sec]'), ylabel('yaw rate r [grad/sec]') + subplot(3,1,3) + plot(Ta, Ya(:, 7)*180/pi), grid on; + xlabel('time [sec]'), ylabel('roll rate p [grad/sec]') + + + % angle + alpha = zeros(length(Ta),1); + beta = zeros(length(Ta),1); + roll = zeros(length(Ta),1); + + for k = 2:length(Ta) + alpha(k) = alpha(k-1) + (Ya(k, 8) + Ya(k-1, 8))/2*180/pi*(T(k)-T(k-1)); + beta(k) = beta(k-1) + (Ya(k, 9) + Ya(k-1, 9))/2*180/pi*(T(k)-T(k-1)); + roll(k) = roll(k-1) + (Ya(k, 7) + Ya(k-1, 7))/2*180/pi*(T(k)-T(k-1)); + end + figure(), plot(Ta, alpha+settings.OMEGA*180/pi), title('Alpha vs time'), grid on; + figure(), plot(Ta, beta), title('Beta vs time'), grid on; + figure(), plot(Ta, roll), title('Roll vs time'), grid on; + + + % Mach + figure(), plot(T, M), grid on; + xlabel('time [sec]'), ylabel('Mach [-]'); + + % Stagnation Temperature + figure(); + plot(T, Tamb-273.15); + hold on, grid on; + plot(T, Ttot-273.15) + title('Temperature profile') + xlabel('time [sec]'), ylabel('Temperature [�C]'); + legend('Surrounding', 'Total', 'location', 'best') + +end diff --git a/src_old/std_run.m b/src_old/std_run.m new file mode 100644 index 00000000..224b7020 --- /dev/null +++ b/src_old/std_run.m @@ -0,0 +1,283 @@ +function [ Tf,Yf, Ta,Ya ] = std_run( settings ) +%STD RUN - This function runs a standard (non-stochastic) simulation +% OUTPUTS +% Tf: Time steps +% Yf: Final State + +% Author: Ruben Di Battista +% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu +% email: ruben.dibattista@skywarder.eu +% Website: http://www.skywarder.eu +% April 2014; Last revision: 31.XII.2014 +% License: 2-clause BSD + +% Author: Francesco Colombi +% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu +% email: francesco.colombi@skywarder.eu +% Release date: 16/04/2016 + +%Starting Attitude +Q0 = angle2quat(settings.PHI,settings.OMEGA,0*pi/180,'ZYX')'; + +%Starting State +X0 = [0 0 0]'; +V0 = [0 0 0]'; +W0 = [0 0 0]'; +X0a = [X0;V0;W0;Q0;settings.m0;settings.Ixxf;settings.Iyyf;settings.Izzf]; + +%% Wind Generation +if settings.wind.model + 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 %% +if not(settings.ldf) +para = 3; %Flag for Main (Rogall) +end +%Initial Condition are the last from drogue 2 descent +X0m = Yd2(end,:); +[Trog,Yrog] = ode45(@descent_parachute,settings.ode.timemain,X0m,... + settings.ode.optionsmain,settings,uw,vw,ww,para); + + + +%% FINAL STATE ASSEMBLING %% + +%Total State +Yf = [Ya(:,1:3) quatrotate(quatconj(Ya(:,10:13)),Ya(:,4:6));Yd1; Yd2;Yrog]; + +%global Par1 Par2 Par3 Asc +%Asc = [Ya(:,1:3) quatrotate(quatconj(Ya(:,10:13)),Ya(:,4:6))]; +%Par1 = Yd1; +%Par2 = Yd2; +%Par3 = Yrog; + +%Total Time +Tf = [Ta; Ta(end)+Td1; Ta(end)+Td1(end)+Td2;Ta(end)+Td1(end)+Td2(end)+Trog]; + +%Parachute State +Yp = [Yd1;Yd2;Yrog]; +Tp = [Ta(end)+Td1;Ta(end)+Td1(end)+Td2;Ta(end)+Td1(end)+Td2(end)+Trog]; + + +%% PLOTTING THINGS %% + +if settings.plot == 1 + + set(0,'DefaultAxesFontSize',settings.DefaultFontSize,... + 'DefaultLineLineWidth',settings.DefaultLineWidth); + + % ASCENT % + plot(Tf,-Yf(:,3)+settings.z0,'k-','LineWidth',2) + hold on + title('Altitude Profile on Time'); + xlabel('Time [s]') + ylabel('Altitude [m]'); + h1=plot(Ta(end),-Ya(end,3)+settings.z0,'ko','MarkerSize',8,... + 'MarkerFaceColor','k'); + h2=plot(Ta(end)+Td1(end),-Yd1(end,3)+settings.z0,'ks','MarkerSize',8,... + 'MarkerFaceColor','k'); + h3=plot(Ta(end)+Td1(end)+Td2(end),-Yd2(end,3)+settings.z0,'ks','MarkerSize',8,... + 'MarkerFaceColor','k'); + legend([h1 h2 h3],'Apogee/1st Drogue Deployment',... + '2nd Drogue Deployment', 'Main Parachute Deployment') + grid on + + %Interpolation for less points --> visualization + Tinterp = linspace(Tf(1),Tf(end),settings.tSteps); + [Tunique,ia,~]= unique(Tf); + u = interp1(Tunique,Yf(ia,4),Tinterp); + v = interp1(Tunique,Yf(ia,5),Tinterp); + w = interp1(Tunique,Yf(ia,6),Tinterp); + + % VELOCITIES % + Va=quatrotate(quatconj(Ya(end,10:13)),Ya(end,4:6)); %apogee + + figure; + suptitle('Velocities Profiles on Time') + subplot(3,1,1); + plot(Tinterp,u,'k-') + hold on + xlabel('Time[s]'); + ylabel('Velocity-x [m/s]'); + h1=plot(Ta(end),Va(1),'ko','MarkerSize',8,... + 'MarkerFaceColor','k'); + h2=plot(Ta(end)+Td1(end),Yd1(end,4),'ks','MarkerSize',8,... + 'MarkerFaceColor','k'); + h3=plot(Ta(end)+Td1(end)+Td2(end),Yd2(end,4),'ks','MarkerSize',8,... + 'MarkerFaceColor','k'); + legend([h1,h2,h3],'Apogee/1st Drogue Deployment',... + '2nd Drogue Deployment','Main Parachute Deployment',... + 'Location','southeast'); + grid on + + subplot(3,1,2) + plot(Tinterp,v,'k-') + hold on + xlabel('Time[s]'); + ylabel('Velocity-y [m/s]'); + h1=plot(Ta(end),Va(2),'ko','MarkerSize',8,... + 'MarkerFaceColor','k'); + h2=plot(Ta(end)+Td1(end),Yd1(end,5),'ks','MarkerSize',8,... + 'MarkerFaceColor','k'); + h3=plot(Ta(end)+Td1(end)+Td2(end),Yd2(end,5),'ks','MarkerSize',8,... + 'MarkerFaceColor','k'); + legend([h1,h2,h3],'Apogee/1st Drogue Deployment',... + '2nd Drogue Deployment','Main Parachute Deployment',... + 'Location','southeast'); + grid on + + + subplot(3,1,3) + plot(Tinterp,w,'k-') + hold on + xlabel('Time[s]'); + ylabel('Velocity-z [m/s]'); + h1=plot(Ta(end),Va(3),'ko','MarkerSize',8,... + 'MarkerFaceColor','k'); + h2=plot(Ta(end)+Td1(end),Yd1(end,6),'ks','MarkerSize',8,... + 'MarkerFaceColor','k'); + h3=plot(Ta(end)+Td1(end)+Td2(end),Yd2(end,6),'ks','MarkerSize',8,... + 'MarkerFaceColor','k'); + legend([h1,h2,h3],'Apogee/1st Drogue Deployment',... + '2nd Drogue Deployment','Main Parachute Deployment',... + 'Location','southeast'); + grid on + + %COMPLETE TRAJECTORY% + figure; + h0=plot3(Yf(1,2),Yf(1,1),-Yf(1,3)+settings.z0,'k+','MarkerSize',10); + hold on + plot3(Yf(:,2),Yf(:,1),-Yf(:,3)+settings.z0,'k-','LineWidth',2) + title('Complete Trajectory'); + xlabel('East [m]') + ylabel('Noth [m]'); + zlabel('Altitude [m]'); + h1=plot3(Ya(end,2),Ya(end,1),-Ya(end,3)+settings.z0,'ko','MarkerSize',8,... + 'MarkerFaceColor','k'); + h2=plot3(Yd1(end,2),Yd1(end,1),-Yd1(end,3)+settings.z0,'ks','MarkerSize',8,... + 'MarkerFaceColor','k'); + h3=plot3(Yd2(end,2),Yd2(end,1),-Yd2(end,3)+settings.z0,'ks','MarkerSize',8,... + 'MarkerFaceColor','k'); + legend([h0,h1,h2,h3],'Launch Pad','Apogee/1st Drogue Deployment',... + '2nd Drogue Deployment','Main Parachute Deployment'); + grid on + + % PLANAR DISPLACEMENT % + + figure + plot(Yf(:,2),Yf(:,1),'k-','LineWidth',2) + hold on + title('Planar Displacement'); + xlabel('East [m]') + ylabel('North [m]'); + h1=plot(Ya(end,2),Ya(end,1),'ko','MarkerSize',8,... + 'MarkerFaceColor','k'); + h2=plot(Yd1(end,2),Yd1(end,1),'ks','MarkerSize',8,... + 'MarkerFaceColor','k'); + h3=plot(Yd2(end,2),Yd2(end,1),'ks','MarkerSize',8,... + 'MarkerFaceColor','k'); + legend([h1,h2,h3],'Apogee/1st Drogue Deployment',... + '2nd Drogue Deployment','Main Parachute Deployment',... + 'Location','southeast'); + grid on + + % PARACHUTES % + figure + plot(Tp,-Yp(:,3)+settings.z0,'k-','LineWidth',2) + hold on + title('Altitude Profile on Time (parachutes)'); + xlabel('Time [s]') + ylabel('Altitude [m]'); + h1=plot(Ta(end),-Ya(end,3)+settings.z0,'ko','MarkerSize',8,... + 'MarkerFaceColor','k'); + h2=plot(Ta(end)+Td1(end),-Yd1(end,3)+settings.z0,'ks','MarkerSize',8,... + 'MarkerFaceColor','k'); + h3=plot(Ta(end)+Td1(end)+Td2(end),-Yd2(end,3)+settings.z0,'ks','MarkerSize',8,... + 'MarkerFaceColor','k'); + legend([h1,h2,h3],'Apogee/1st Drogue Deployment',... + '2nd Drogue Deployment','Main Parachute Deployment',... + 'Location','southeast'); + grid on + + + + %Interpolation for less points --> visualization + Tinterp = linspace(Ta(1),Ta(end),settings.tSteps); + p = interp1(Ta,Ya(:,7),Tinterp); + q = interp1(Ta,Ya(:,8),Tinterp); + r = interp1(Ta,Ya(:,9),Tinterp); + + % ANGULAR RATES % + + figure; + suptitle('Angular rates on Time') + subplot(3,1,1); + plot(Tinterp,p,'k-') + hold on + xlabel('Time[s]'); + ylabel('Roll rate p [rad/s]'); + h1=plot(Ta(end),Ya(end,7),'ko','MarkerSize',8,... + 'MarkerFaceColor','k'); + legend(h1,'Apogee/1st Drogue Deployment','Location','southeast'); + grid on + + subplot(3,1,2) + plot(Tinterp,q,'k-') + hold on + xlabel('Time[s]'); + ylabel('Pitch rate q [rad/s]'); + h1=plot(Ta(end),Ya(end,8),'ko','MarkerSize',8,... + 'MarkerFaceColor','k'); + legend(h1,'Apogee/1st Drogue Deployment','Location','northeast'); + grid on + + + subplot(3,1,3) + plot(Tinterp,r,'k-') + hold on + xlabel('Time[s]'); + ylabel('Yaw rate r [rad/s]'); + h1=plot(Ta(end),Ya(end,9),'ko','MarkerSize',8,... + 'MarkerFaceColor','k'); + legend(h1,'Apogee/1st Drogue Deployment','Location','southeast'); + grid on + +end + +%Resizing +h = get(0,'children'); +scrsz = get(0,'ScreenSize'); +for i=1:length(h) + set(h(i),'OuterPosition',[0 0 scrsz(4) scrsz(4)]) + %saveas(h(i), ['figure' num2str(i)], 'fig'); +end +end diff --git a/src_old/std_run_ballistic.m b/src_old/std_run_ballistic.m new file mode 100644 index 00000000..2bf83adb --- /dev/null +++ b/src_old/std_run_ballistic.m @@ -0,0 +1,229 @@ +function [ Tf,Yf, Ta,Ya ] = std_run_ballistic( settings ) +%STD RUN - This function runs a standard (non-stochastic) simulation +% OUTPUTS +% Tf: Time steps +% Yf: Final State + +% Author: Ruben Di Battista +% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu +% email: ruben.dibattista@skywarder.eu +% Website: http://www.skywarder.eu +% April 2014; Last revision: 31.XII.2014 +% License: 2-clause BSD + +% Author: Francesco Colombi +% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu +% email: francesco.colombi@skywarder.eu +% Release date: 16/04/2016 + +%Starting Attitude +Q0 = angle2quat(settings.PHI,settings.OMEGA,0*pi/180,'ZYX')'; + +%Starting State +X0 = [0 0 0]'; +V0 = [0 0 0]'; +W0 = [0 0 0]'; +X0a = [X0;V0;W0;Q0;settings.m0;settings.Ixxf;settings.Iyyf;settings.Izzf]; + +% Wind Generation +[uw,vw,ww] = windgen(settings.wind.AzMin,settings.wind.AzMax,... + settings.wind.ElMin,settings.wind.ElMax,settings.wind.MagMin,... + settings.wind.MagMax); + +%% ASCEND %% + +[Ta,Ya] = ode45(@ascend,settings.ode.timeasc,X0a,settings.ode.optionsasc,... + settings,uw,vw,ww); + +%% DESCEND %% + +[Td,Yd] = ode45(@ballistic_descent,settings.ode.timedesc,Ya(end,:),settings.ode.optionsdesc,... + settings,uw,vw,ww); + +% %% DROGUE 1 %% +% para = 1; %Flag for Drogue 1 +% +% %Initial Condition are the last from ascend (need to rotate because +% %velocities are outputted in body axes) +% X0d1 = [Ya(end,1:3) quatrotate(quatconj(Ya(end,10:13)),Ya(end,4:6))]; +% [Td1,Yd1] = ode45(@descent_parachute,settings.ode.timedrg1,X0d1,... +% settings.ode.optionsdrg1,settings,uw,vw,ww,para); +% +% %% DROGUE 2 %% +% para = 2; %Flag for Drogue 2 +% +% %Initial Condition are the last from drogue 1 descent +% X0d2 = Yd1(end,:); +% [Td2,Yd2] = ode45(@descent_parachute,settings.ode.timedrg2,X0d2,... +% settings.ode.optionsdrg2,settings,uw,vw,ww,para); +% +% %% MAIN %% +% para = 3; %Flag for Main (Rogall) +% +% %Initial Condition are the last from drogue 2 descent +% X0m = Yd2(end,:); +% [Trog,Yrog] = ode45(@descent_parachute,settings.ode.timemain,X0m,... +% settings.ode.optionsmain,settings,uw,vw,ww,para); + +%% FINAL STATE ASSEMBLING %% + +%Total State +Yf = [Ya;Yd]; +%Total Time +Tf = [Ta; Ta(end)+Td]; + +%% PLOTTING THINGS %% + +if settings.plot == 1 + + set(0,'DefaultAxesFontSize',settings.DefaultFontSize,... + 'DefaultLineLineWidth',settings.DefaultLineWidth); + + %Interpolation for less points --> visualization + Tinterp = linspace(Tf(1),Tf(end),settings.tSteps); + [Tunique,ia,~]= unique(Tf); + u = interp1(Tunique,Yf(ia,4),Tinterp); + v = interp1(Tunique,Yf(ia,5),Tinterp); + w = interp1(Tunique,Yf(ia,6),Tinterp); + + % VELOCITIES % + Va=quatrotate(quatconj(Ya(end,10:13)),Ya(end,4:6)); %apogee + + figure; + subtitle('Velocities Profiles on Time') + subplot(3,1,1); + plot(Tinterp,u,'k-') + hold on + xlabel('Time[s]'); + ylabel('Velocity-x [m/s]'); + h1=plot(Ta(end),Va(1),'ko','MarkerSize',8,... + 'MarkerFaceColor','k'); + h2=plot(Ta(end)+Td(end),Yd(end,4),'ks','MarkerSize',8,... + 'MarkerFaceColor','k'); + legend([h1,h2],'Apogee/1st Drogue Deployment',... + '2nd Drogue Deployment',... + 'Location','southeast'); + grid on + + subplot(3,1,2) + plot(Tinterp,v,'k-') + hold on + xlabel('Time[s]'); + ylabel('Velocity-y [m/s]'); + h1=plot(Ta(end),Va(2),'ko','MarkerSize',8,... + 'MarkerFaceColor','k'); + h2=plot(Ta(end)+Td(end),Yd(end,5),'ks','MarkerSize',8,... + 'MarkerFaceColor','k'); + legend([h1,h2],'Apogee/1st Drogue Deployment',... + '2nd Drogue Deployment',... + 'Location','southeast'); + grid on + + + subplot(3,1,3) + plot(Tinterp,w,'k-') + hold on + xlabel('Time[s]'); + ylabel('Velocity-z [m/s]'); + h1=plot(Ta(end),Va(3),'ko','MarkerSize',8,... + 'MarkerFaceColor','k'); + h2=plot(Ta(end)+Td(end),Yd(end,6),'ks','MarkerSize',8,... + 'MarkerFaceColor','k'); + legend([h1,h2],'Apogee/1st Drogue Deployment',... + '2nd Drogue Deployment',... + 'Location','southeast'); + grid on + + %COMPLETE TRAJECTORY% + figure; + h0=plot3(Yf(1,2),Yf(1,1),-Yf(1,3)+settings.z0,'k+','MarkerSize',10); + hold on + plot3(Yf(:,2),Yf(:,1),-Yf(:,3)+settings.z0,'k-','LineWidth',2) + title('Complete Trajectory'); + xlabel('East [m]') + ylabel('Noth [m]'); + zlabel('Altitude [m]'); + h1=plot3(Ya(end,2),Ya(end,1),-Ya(end,3)+settings.z0,'ko','MarkerSize',8,... + 'MarkerFaceColor','k'); + h2=plot3(Yd(end,2),Yd(end,1),-Yd(end,3)+settings.z0,'ks','MarkerSize',8,... + 'MarkerFaceColor','k'); + + legend([h0,h1,h2],'Launch Pad','Apogee/1st Drogue Deployment',... + '2nd Drogue Deployment','Main Parachute Deployment'); + grid on + + % PLANAR DISPLACEMENT % + + figure + plot(Yf(:,2),Yf(:,1),'k-','LineWidth',2) + hold on + title('Planar Displacement'); + xlabel('East [m]') + ylabel('North [m]'); + h1=plot(Ya(end,2),Ya(end,1),'ko','MarkerSize',8,... + 'MarkerFaceColor','k'); + h2=plot(Yd(end,2),Yd(end,1),'ks','MarkerSize',8,... + 'MarkerFaceColor','k'); + + legend([h1,h2],'Apogee/1st Drogue Deployment',... + '2nd Drogue Deployment',... + 'Location','southeast'); + grid on + + + + + %Interpolation for less points --> visualization + Tinterp = linspace(Ta(1),Ta(end),settings.tSteps); + p = interp1(Ta,Ya(:,7),Tinterp); + q = interp1(Ta,Ya(:,8),Tinterp); + r = interp1(Ta,Ya(:,9),Tinterp); + + % ANGULAR RATES % + + figure; + suptitle('Angular rates on Time') + subplot(3,1,1); + plot(Tinterp,p,'k-') + hold on + xlabel('Time[s]'); + ylabel('Roll rate p [rad/s]'); + h1=plot(Ta(end),Ya(end,7),'ko','MarkerSize',8,... + 'MarkerFaceColor','k'); + legend(h1,'Apogee/1st Drogue Deployment','Location','southeast'); + grid on + + subplot(3,1,2) + plot(Tinterp,q,'k-') + hold on + xlabel('Time[s]'); + ylabel('Pitch rate q [rad/s]'); + h1=plot(Ta(end),Ya(end,8),'ko','MarkerSize',8,... + 'MarkerFaceColor','k'); + legend(h1,'Apogee/1st Drogue Deployment','Location','northeast'); + grid on + + + subplot(3,1,3) + plot(Tinterp,r,'k-') + hold on + xlabel('Time[s]'); + ylabel('Yaw rate r [rad/s]'); + h1=plot(Ta(end),Ya(end,9),'ko','MarkerSize',8,... + 'MarkerFaceColor','k'); + legend(h1,'Apogee/1st Drogue Deployment','Location','southeast'); + grid on + +end + +%Resizing +h = get(0,'children'); +scrsz = get(0,'ScreenSize'); +for i=1:length(h) + set(h(i),'OuterPosition',[0 0 scrsz(4) scrsz(4)]) + %saveas(h(i), ['figure' num2str(i)], 'fig'); +end + + + +end \ No newline at end of file diff --git a/src_old/stoch_run.m b/src_old/stoch_run.m new file mode 100644 index 00000000..eb17b5c8 --- /dev/null +++ b/src_old/stoch_run.m @@ -0,0 +1,165 @@ +function [ LP ,Z ] = stoch_run( settings ) +%STD RUN - This function runs a stochastic simulation +% OUTPUTS +% LP: Landing Points +% Z: Apogee Altitudes + +% Author: Ruben Di Battista +% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu +% email: ruben.dibattista@skywarder.eu +% Website: http://www.skywarder.eu +% April 2014; Last revision: 29.V.2014 +% License: 2-clause BSD + +%Starting Attitude +Q0 = angle2quat(settings.PHI,settings.OMEGA,0*pi/180,'ZYX')'; + +%Starting State +X0 = [0 0 0]'; +V0 = [0 0 0]'; +W0 = [0 0 0]'; +X0a = [X0;V0;W0;Q0;settings.m0;settings.Ixxf;settings.Iyyf;settings.Izzf]; + + +%PreAllocation +LP = zeros(settings.stoch.N,3); +Z = zeros(settings.stoch.N,1); +ApoTime = zeros(settings.stoch.N,1); + + +parfor_progress(settings.stoch.N); +for i=1:settings.stoch.N + %% Wind Generation + [uw,vw,ww] = windgen(settings.wind.AzMin,settings.wind.AzMax,... + settings.wind.ElMin,settings.wind.ElMax,settings.wind.MagMin,... + settings.wind.MagMax); + + %% ASCEND %% + + [Ta,Ya] = ode45(@ascend,settings.ode.timeasc,X0a,settings.ode.optionsasc,... + settings,uw,vw,ww); + + %% DROGUE 1 %% + para = 1; %Flag for Drogue 1 + + %Initial Condition are the last from ascend (need to rotate because + %velocities are outputted in body axes) + X0d1 = [Ya(end,1:3) quatrotate(quatconj(Ya(end,10:13)),Ya(end,4:6))]; + [Td1,Yd1] = ode45(@descent_parachute,settings.ode.timedrg1,X0d1,... + settings.ode.optionsdrg1,settings,uw,vw,ww,para); + + %% DROGUE 2 %% + para = 2; %Flag for Drogue 2 + + %Initial Condition are the last from drogue 1 descent + X0d2 = Yd1(end,:); + [Td2,Yd2] = ode45(@descent_parachute,settings.ode.timedrg2,X0d2,... + settings.ode.optionsdrg2,settings,uw,vw,ww,para); + + %% MAIN %% + if not(settings.ldf) + para = 3; %Flag for Main (Rogall) + end + + %Initial Condition are the last from drogue 2 descent + X0m = Yd2(end,:); + [Tm,Ym] = ode45(@descent_parachute,settings.ode.timemain,X0m,... + settings.ode.optionsmain,settings,uw,vw,ww,para); + + %% FINAL STATE ASSEMBLING %% + + %Total State + if settings.ao + Yf = [Ya(:,1:3) quatrotate(quatconj(Ya(:,10:13)),Ya(:,4:6))]; + else + Yf = [Ya(:,1:3) quatrotate(quatconj(Ya(:,10:13)),Ya(:,4:6));Yd1; Yd2;Ym]; + end + %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,'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('North [m]'); + ylabel('East [m]'); + legend('Mean Landing Point','Launch Site','Landing Points'); + view(90,270) + axis equal + + %Histogram + [f,x] = hist(Z,10); + figure; + bar(x,f/settings.stoch.N); + title('Apogee Altitudes Distribution') + xlabel('Apogee [m]') + ylabel('n_i/n_{tot}') + + + + +end + +%Resizing +h = get(0,'children'); +scrsz = get(0,'ScreenSize'); +for i=1:length(h) + set(h(i),'OuterPosition',[0 0 scrsz(4) scrsz(4)]) + %saveas(h(i), ['figure' num2str(i)], 'fig'); +end + + + + +end diff --git a/src_old/stoch_run_bal.m b/src_old/stoch_run_bal.m new file mode 100644 index 00000000..8292038e --- /dev/null +++ b/src_old/stoch_run_bal.m @@ -0,0 +1,142 @@ +function [ LP ,Z ] = stoch_run_bal( 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); + + +%% DESCEND %% + +[Td,Yd] = ode45(@ballistic_descent,settings.ode.timedesc,Ya(end,:),settings.ode.optionsdesc,... + settings,uw,vw,ww); + + + %% FINAL STATE ASSEMBLING %% + + %Total State + Yf = [Ya;Yd]; + + %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,'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('North [m]'); + ylabel('East [m]'); + legend('Mean Landing Point','Launch Site','Landing Points'); + view(90,270) + axis equal + + %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 \ No newline at end of file diff --git a/src_old/stoch_run_bal_p.m b/src_old/stoch_run_bal_p.m new file mode 100644 index 00000000..f3206b74 --- /dev/null +++ b/src_old/stoch_run_bal_p.m @@ -0,0 +1,139 @@ +function [ LP ,Z ] = stoch_run_bal_p( settings ) +%STD RUN - This function runs a stochastic simulation (parallel) +% OUTPUTS +% LP: Landing Points +% Z: Apogee Altitudes + +% Author: Ruben Di Battista +% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu +% email: ruben.dibattista@skywarder.eu +% Website: http://www.skywarder.eu +% April 2014; Last revision: 29.V.2014 +% License: 2-clause BSD + +%Starting Attitude +Q0 = angle2quat(settings.PHI,settings.OMEGA,0*pi/180,'ZYX')'; + +%Starting State +X0 = [0 0 0]'; +V0 = [0 0 0]'; +W0 = [0 0 0]'; +X0a = [X0;V0;W0;Q0;settings.m0;settings.Ixxf;settings.Iyyf;settings.Izzf]; + + +%PreAllocation +LP = zeros(settings.stoch.N,3); +Z = zeros(settings.stoch.N,1); +ApoTime = zeros(settings.stoch.N,1); + +parfor_progress(settings.stoch.N); +parpool; +parfor i=1:settings.stoch.N + %% Wind Generation + [uw,vw,ww] = windgen(settings.wind.AzMin,settings.wind.AzMax,... + settings.wind.ElMin,settings.wind.ElMax,settings.wind.MagMin,... + settings.wind.MagMax); + + %% ASCEND %% + +[Ta,Ya] = ode45(@ascend,settings.ode.timeasc,X0a,settings.ode.optionsasc,... + settings,uw,vw,ww); + +%% DESCEND %% + +[Td,Yd] = ode45(@ballistic_descent,settings.ode.timedesc,Ya(end,:),settings.ode.optionsdesc,... + settings,uw,vw,ww); + + + %% FINAL STATE ASSEMBLING %% + + %Total State + Yf = [Ya;Yd]; + + %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('North [m]'); + ylabel('East [m]'); + legend('Mean Landing Point','Launch Site','Landing Points'); + view(90,270) + axis equal + + %Histogram + [f,x] = hist(Z,10); + figure; + bar(x,f/settings.stoch.N); + title('Apogee Altitudes Distribution') + xlabel('Apogee [m]') + ylabel('n_i/n_{tot}') + + + + +end + +delete(gcp); + +%Resizing +h = get(0,'children'); +scrsz = get(0,'ScreenSize'); +for i=1:length(h) + set(h(i),'OuterPosition',[0 0 scrsz(4) scrsz(4)]) + %saveas(h(i), ['figure' num2str(i)], 'fig'); +end + + + +end diff --git a/src_old/stoch_run_p.m b/src_old/stoch_run_p.m new file mode 100644 index 00000000..0e786db0 --- /dev/null +++ b/src_old/stoch_run_p.m @@ -0,0 +1,163 @@ +function [ LP ,Z ] = stoch_run_p( settings ) +%STD RUN - This function runs a stochastic simulation (parallel) +% OUTPUTS +% LP: Landing Points +% Z: Apogee Altitudes + +% Author: Ruben Di Battista +% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu +% email: ruben.dibattista@skywarder.eu +% Website: http://www.skywarder.eu +% April 2014; Last revision: 29.V.2014 +% License: 2-clause BSD + +%Starting Attitude +Q0 = angle2quat(settings.PHI,settings.OMEGA,0*pi/180,'ZYX')'; + +%Starting State +X0 = [0 0 0]'; +V0 = [0 0 0]'; +W0 = [0 0 0]'; +X0a = [X0;V0;W0;Q0;settings.m0;settings.Ixxf;settings.Iyyf;settings.Izzf]; + + +%PreAllocation +LP = zeros(settings.stoch.N,3); +Z = zeros(settings.stoch.N,1); +ApoTime = zeros(settings.stoch.N,1); + +parfor_progress(settings.stoch.N); +parpool; +parfor i=1:settings.stoch.N + %% Wind Generation + [uw,vw,ww] = windgen(settings.wind.AzMin,settings.wind.AzMax,... + settings.wind.ElMin,settings.wind.ElMax,settings.wind.MagMin,... + settings.wind.MagMax); + + %% ASCEND %% + + [Ta,Ya] = ode45(@ascend,settings.ode.timeasc,X0a,settings.ode.optionsasc,... + settings,uw,vw,ww); + + %% DROGUE 1 %% + para = 1; %Flag for Drogue 1 + + %Initial Condition are the last from ascend (need to rotate because + %velocities are outputted in body axes) + X0d1 = [Ya(end,1:3) quatrotate(quatconj(Ya(end,10:13)),Ya(end,4:6))]; + [Td1,Yd1] = ode45(@descent_parachute,settings.ode.timedrg1,X0d1,... + settings.ode.optionsdrg1,settings,uw,vw,ww,para); + + %% DROGUE 2 %% + para = 2; %Flag for Drogue 2 + + %Initial Condition are the last from drogue 1 descent + X0d2 = Yd1(end,:); + [Td2,Yd2] = ode45(@descent_parachute,settings.ode.timedrg2,X0d2,... + settings.ode.optionsdrg2,settings,uw,vw,ww,para); + + %% MAIN %% + if not(settings.ldf) + para = 3; %Flag for Main (Rogall) + end + + %Initial Condition are the last from drogue 2 descent + X0m = Yd2(end,:); + [Tm,Ym] = ode45(@descent_parachute,settings.ode.timemain,X0m,... + settings.ode.optionsmain,settings,uw,vw,ww,para); + + %% FINAL STATE ASSEMBLING %% + + %Total State + if settings.ao + Yf = [Ya(:,1:3) quatrotate(quatconj(Ya(:,10:13)),Ya(:,4:6))]; + else + Yf = [Ya(:,1:3) quatrotate(quatconj(Ya(:,10:13)),Ya(:,4:6));Yd1; Yd2;Ym]; + end + %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('North [m]'); + ylabel('East [m]'); + legend('Mean Landing Point','Launch Site','Landing Points'); + view(90,270) + axis equal + + %Histogram + [f,x] = hist(Z,10); + figure; + bar(x,f/settings.stoch.N); + title('Apogee Altitudes Distribution') + xlabel('Apogee [m]') + ylabel('n_i/n_{tot}') + + + + +end + +delete(gcp); + +%Resizing +h = get(0,'children'); +scrsz = get(0,'ScreenSize'); +for i=1:length(h) + set(h(i),'OuterPosition',[0 0 scrsz(4) scrsz(4)]) + %saveas(h(i), ['figure' num2str(i)], 'fig'); +end + + + +end diff --git a/src/vert_windgen.m b/src_old/vert_windgen.m similarity index 100% rename from src/vert_windgen.m rename to src_old/vert_windgen.m diff --git a/src/wind_generator.m b/src_old/wind_generator.m similarity index 100% rename from src/wind_generator.m rename to src_old/wind_generator.m diff --git a/src/windgen.m b/src_old/windgen.m similarity index 100% rename from src/windgen.m rename to src_old/windgen.m diff --git a/utils/launch_pad.m b/utils/launch_pad.m new file mode 100644 index 00000000..5f71c3f0 --- /dev/null +++ b/utils/launch_pad.m @@ -0,0 +1,85 @@ +% Author: Francesco Colombi +% Skyward Experimental Rocketry | CRD Dept | crd@skywarder.eu +% email: francesco.colombi@skywarder.eu +% Release date: 16/04/2016 + +close all +clear, clc + +% inclinazione rampa +theta = 85*pi/180; + +% [vv]horizon = R'*[vv]corpo +% [vv]corpo = R*[vv]horizon +R = [ cos(theta), 0, sin(theta); ... + 0, 1, 0; ... + -sin(theta), 0, cos(theta)]; + +% velocit� di uscita dalla rampa +Lrampa = 4.5; % m +acc = 5*9.80665; % m/s2 +texit = sqrt(Lrampa/(.5*acc)); % sec +V = acc*texit % m/s +% V = 32;% m/s +vv = V*[1, 0, 0]'; + +% vento +% incidenza vento +dev = [-180:5:180]*pi/180; +W = sqrt(2^2+5^2) * 1.2 % m/s +% W = 10 + +[~, a, ~, rho] = atmoscoesa(0); + +C = 0.174; % m Caliber (Fuselage Diameter) +S = 0.02378; % m2 Cross-sectional Surface +L = 4.32050; % m lunghezza missile + +load for006_full.mat +% Coeffs. Interpolation +givA = State.Alphas; +givB = State.Betas; +givH = State.Altitudes; +givM = State.Machs; + +norm_vr = zeros(length(dev), 1); +M = zeros(length(dev), 1); +alpha = zeros(length(dev), 1); +beta = zeros(length(dev), 1); +XCP = zeros(length(dev), 1); +CM = zeros(length(dev), 1); +Mpitch = zeros(length(dev), 1); + +for n = 1:length(dev) + ww = -W*[cos(dev(n)), sin(dev(n)), 0]'; % assi horizon + ww = R*ww; % assi corpo + + % vel relativa al vento + vvr = vv - ww; % assi corpo + norm_vr(n) = norm(vvr); + M(n) = norm_vr(n)/a; + + ur = vvr(1); + vr = vvr(2); + wr = vvr(3); + + % alpha + alpha(n) = atan(-wr/ur) *180/pi; % deg + beta(n) = asin(vr/norm(vvr)) *180/pi; % deg + + XCP(n) = interpn(givA,givM,givB,givH,Coeffs.X_C_P,alpha(n),M(n),beta(n),0); + CM(n) = interpn(givA,givM,givB,givH,Coeffs.CM,alpha(n),M(n),beta(n),0); + Mpitch(n) = 0.5*rho*norm_vr(n)^2*S*C*CM(n); +end + +figure() +plot(dev*180/pi, XCP), grid on +xlabel('angolo vento [deg]'), ylabel('XCP [ref length]') + +figure() +plot(dev*180/pi, CM), grid on +xlabel('angolo vento [deg]'), ylabel('CM [-]') + +figure() +plot(dev*180/pi, Mpitch), grid on +xlabel('angolo vento [deg]'), ylabel('M pitch [Nm]') -- GitLab