diff --git a/functions/eventFunctions/eventPad.m b/functions/eventFunctions/eventPad.m
index 9048dc893cffc76c77ad2df182c3fc9f6ce3d655..6f9a8800c77ccee67994fa008df7267deedc31c9 100644
--- a/functions/eventFunctions/eventPad.m
+++ b/functions/eventFunctions/eventPad.m
@@ -1,25 +1,25 @@
function [value, isterminal, direction] = eventPad(~, Y, ~, environment, varargin)
% eventPad - Event function to stop simulation at the launchpad exit
-%
+%
% INPUTS:
% - t, double [1, 1], integration time [s];
-% - Y, double [4, 1], integration state, check launchPadFreeDyn for explanation;
-% - settings, struct (motor, CoeffsE, CoeffsF, para, ode, stoch, prob, wind),
+% - Y, double [4, 1], integration state, check launchPadFreeDyn for explanation;
+% - settings, struct (motor, CoeffsE, CoeffsF, para, ode, stoch, prob, wind),
% simulation data;
% - varargin, cell {3, 1}, for additional inputs of the ode function, un-necessary here.
-%
+%
% OUTPUTS:
% - value, double [1, 1], see eventFunction explantion in matlab website
% - isterminal, double [1, 1], see eventFunction explantion in matlab website
% - direction, double [1, 1], see eventFunction explantion in matlab website
-%
+%
% CALLED FUNCTIONS: /
-%
+%
% REVISIONS:
% - 0 21/10/20, release Adriano Filippo Inno
% Copyright © 2021, Skyward Experimental Rocketry, AFD department
% All rights reserved
-%
+%
% SPDX-License-Identifier: GPL-3.0-or-later
altitude = -Y(3);
diff --git a/functions/miscellaneous/delayControl.m b/functions/miscellaneous/delayControl.m
index 96cdfe3132061a186524f9142dec3023bf1a2028..48b52725fb22c63a02b4bc3d332c02264ecf3dbe 100644
--- a/functions/miscellaneous/delayControl.m
+++ b/functions/miscellaneous/delayControl.m
@@ -1,4 +1,7 @@
-function delay = delayControl(pControl, settings)
+function delay = delayControl(rocket)
+arguments
+ rocket Rocket
+end
%{
delayCOntrol - This function calculates the delay to change airbrakes
@@ -26,16 +29,19 @@ A = -9.43386/(3*1000);
B = 19.86779/(3*1000);
alpha = @(S) (-B + sqrt(B^2 + 4*A*S))/(2*A);
-controlSurf = settings.hprot(pControl) * settings.wprot;
+% settings.hprot = rocket.airbrakes.height
+% settings.wprot = rocket.airbrakes.width
+% pControl = rocket.airbrakes.extension
+% settings.vServo = rocket.airbrakes.servoOmega
+controlSurf = rocket.airbrakes.height(rocket.airbrakes.extension) * rocket.airbrakes.width;
servoAngle = alpha(controlSurf);
-Np = length(pControl);
+Np = length(rocket.airbrakes.extension);
delay = zeros(1, Np-1);
for i = 2:Np
- delay(i-1) = abs(servoAngle(i) - servoAngle(i-1))/settings.vServo;
+ delay(i-1) = abs(servoAngle(i) - servoAngle(i-1))/rocket.airbrakes.servoOmega;
end
-
end
\ No newline at end of file
diff --git a/functions/miscellaneous/gaplotshape.m b/functions/miscellaneous/gaplotshape.m
index 4ed8b84ccd6adec1836b3ff94bc215fc6c7295c4..3366665b87edb471fa48e103fb9268bd81cedebf 100644
--- a/functions/miscellaneous/gaplotshape.m
+++ b/functions/miscellaneous/gaplotshape.m
@@ -1,4 +1,11 @@
-function state = gaplotshape(settings, ~, state, flag)
+function state = gaplotshape(rocket, settings, ~, state, flag)
+arguments
+ rocket Rocket
+ settings Settings
+ ~
+ state struct
+ flag char
+end
%{
GAPLOTSHAPE - Genethic algorithm plot function to visualize current best
rocket shape.
@@ -23,15 +30,14 @@ SPDX-License-Identifier: GPL-3.0-or-later
%}
-
switch flag
case 'init'
hold on
%%% NOSECONE
- r = settings.C/2*1000;
+ r = rocket.diameter/2*1000;
Lnos = state.Population(1, 5)*10;
xNos = Lnos - Lnos*cos(linspace(0, pi/2, 50));
- if settings.modHaack
+ if settings.optimization.modHaack
p = state.Population(1, 6);
c = state.Population(1, 7);
xMod = @(x, p) (x/Lnos).^(p)*Lnos;
@@ -39,7 +45,7 @@ switch flag
haackSeriesMod = @(x, p, C) ( r/sqrt(pi) ) * sqrt( thetaMod(x, p) - ( sin(2*thetaMod(x, p))./ 2 ) + C*sin(thetaMod(x, p)).^3 );
yNos = haackSeriesMod(xNos, p, c);
else
- if settings.boatOptimization
+ if settings.optimization.boatOptimization
ogType = state.Population(1, 7);
else
ogType = state.Population(1, 6);
@@ -74,16 +80,16 @@ switch flag
plotNosEnd = plot([Lnos Lnos], [-r r], 'k');
set(plotNosEnd,'Tag','gaplotNosEnd');
%%% CENTERBODY
- if settings.boatOptimization
- if settings.modHaack
+ if settings.optimization.boatOptimization
+ if settings.optimization.modHaack
Laft = state.Population(1, 8);
- Lcent = settings.Lcenter*1000 + settings.boatL*1000 - Laft;
+ Lcent = rocket.lengthCenterNoMot*1000 + rocket.rear.boatLength*1000 - Laft;
else
Laft = state.Population(1, 6);
- Lcent = settings.Lcenter*1000 + settings.boatL*1000 - Laft;
+ Lcent = rocket.lengthCenterNoMot*1000 + rocket.rear.boatLength*1000 - Laft;
end
else
- Lcent = settings.Lcenter*1000;
+ Lcent = rocket.lengthCenterNoMot*1000;
end
plotCent = plot([Lnos Lnos+Lcent], [r r], 'k');
set(plotCent,'Tag','gaplotCent');
@@ -92,18 +98,18 @@ switch flag
plotCentEnd = plot([Lnos+Lcent Lnos+Lcent], [-r r], 'k');
set(plotCentEnd,'Tag','gaplotCentEnd');
%%% BOAT-TAIL
- Daft = settings.optimBoatD*1000;
- if settings.boatOptimization
- if settings.modHaack
+ Daft = settings.optimization.optimBoatD*1000;
+ if settings.optimization.boatOptimization
+ if settings.optimization.modHaack
Laft = state.Population(1, 8);
else
Laft = state.Population(1, 6);
end
else
- Laft = settings.boatL*1000;
+ Laft = rocket.rear.boatLength*1000;
end
- if settings.boatType == 1
+ if strcmp(rocket.rear.boatType, 'OGIVE') % [-] Boat-tail shape. 0: Cone, 1: Tangent Ogive
[xBoat, yBoat] = computeTangentBoatPoints(2*r, Daft, Laft);
else
xBoat = [0 Laft];
@@ -121,7 +127,7 @@ switch flag
C2 = state.Population(1,2)*10;
deltaXLE = state.Population(1,4)*10;
H = state.Population(1,3)*10;
- Xle1 = Lcent + Lnos - settings.d*1000 - C1;
+ Xle1 = Lcent + Lnos - rocket.rear.finsAxialDistance*1000 - C1;
plotFin1 = plot([Xle1 Xle1+deltaXLE], [r r+H],'k');
set(plotFin1,'Tag','gaplotFin1');
plotFin2 = plot([Xle1+C1 Xle1+deltaXLE+C2], [r r+H],'k');
@@ -141,10 +147,10 @@ switch flag
case 'iter'
[~,i] = min(state.Score);
%%% NOSECONE
- r = settings.C/2*1000;
+ r = rocket.diameter/2*1000;
Lnos = state.Population(i, 5)*10;
xNos = Lnos - Lnos*cos(linspace(0, pi/2, 50));
- if settings.modHaack
+ if settings.optimization.modHaack
p = state.Population(i, 6);
c = state.Population(i, 7);
xMod = @(x, p) (x/Lnos).^(p)*Lnos;
@@ -152,7 +158,7 @@ switch flag
haackSeriesMod = @(x, p, C) ( r/sqrt(pi) ) * sqrt( thetaMod(x, p) - ( sin(2*thetaMod(x, p))./ 2 ) + C*sin(thetaMod(x, p)).^3 );
yNos = haackSeriesMod(xNos, p, c);
else
- if settings.boatOptimization
+ if settings.optimization.boatOptimization
ogType = state.Population(1, 7);
else
ogType = state.Population(1, 6);
@@ -189,16 +195,16 @@ switch flag
h = findobj(get(gca,'Children'),'Tag','gaplotNosEnd');
set(h,'Xdata',[Lnos Lnos]);
%%% CENTERBODY
- if settings.boatOptimization
- if settings.modHaack
+ if settings.optimization.boatOptimization
+ if settings.optimization.modHaack
Laft = state.Population(i, 8);
- Lcent = settings.Lcenter*1000 + settings.boatL*1000 - Laft;
+ Lcent = rocket.lengthCenterNoMot*1000 + rocket.rear.boatLength*1000 - Laft;
else
Laft = state.Population(i, 6);
- Lcent = settings.Lcenter*1000 + settings.boatL*1000 - Laft;
+ Lcent = rocket.lengthCenterNoMot*1000 + rocket.rear.boatLength*1000 - Laft;
end
else
- Lcent = settings.Lcenter*1000;
+ Lcent = rocket.lengthCenterNoMot*1000;
end
h = findobj(get(gca,'Children'),'Tag','gaplotCent');
set(h,'Xdata',[Lnos Lnos+Lcent]);
@@ -207,17 +213,17 @@ switch flag
h = findobj(get(gca,'Children'),'Tag','gaplotCentEnd');
set(h,'Xdata',[Lnos+Lcent Lnos+Lcent]);
%%% BOAT-TAIL
- Daft = settings.optimBoatD*1000;
- if settings.boatOptimization
- if settings.modHaack
+ Daft = settings.optimization.optimBoatD*1000;
+ if settings.optimization.boatOptimization
+ if settings.optimization.modHaack
Laft = state.Population(i, 8);
else
Laft = state.Population(i, 6);
end
else
- Laft = settings.boatL*1000;
+ Laft = rocket.rear.boatLength*1000;
end
- if settings.boatType == 1
+ if strcmp(rocket.rear.boatType, 'OGIVE')
[xBoat, yBoat] = computeTangentBoatPoints(2*r, Daft, Laft);
else
xBoat = [0 Laft];
@@ -238,7 +244,7 @@ switch flag
C2 = state.Population(i,2)*10;
deltaXLE = state.Population(i,4)*10;
H = state.Population(i,3)*10;
- Xle1 = Lcent + Lnos - settings.d*1000 - C1;
+ Xle1 = Lcent + Lnos - rocket.rear.finsAxialDistance*1000 - C1;
h = findobj(get(gca,'Children'),'Tag','gaplotFin1');
set(h,'Xdata',[Xle1 Xle1+deltaXLE]);
set(h,'Ydata',[r r+H]);
@@ -258,7 +264,6 @@ switch flag
set(h,'Xdata',[Xle1+deltaXLE Xle1+deltaXLE+C2]);
set(h,'Ydata',[-r-H, -r-H]);
- set(gca, 'xlim', [-100, Lnos+Lcent+Laft+100]);
-
+ set(gca, 'xlim', [-100, Lnos+Lcent+Laft+100]);
end
-
+end
\ No newline at end of file
diff --git a/functions/miscellaneous/singleGArun.m b/functions/miscellaneous/singleGArun.m
index 4d2dc62dfeb264da04e025882689b6654037a11b..7ec2063f2567ac38f4814c59eb3ce0c2e8f921b0 100644
--- a/functions/miscellaneous/singleGArun.m
+++ b/functions/miscellaneous/singleGArun.m
@@ -1,4 +1,15 @@
-function [time, launchpad, apogee, c] = singleGArun(lb, ub, input, settings)
+function [time, launchpad, apogee, c] = singleGArun(lb, ub, input, data_opt, rocket, environment, wind, settings, wrapper)
+arguments
+ lb double
+ ub double
+ input struct
+ data_opt struct
+ rocket Rocket
+ environment Environment
+ wind WindCustom
+ settings Settings
+ wrapper DataWrapper
+end
%{
singelGArun - estimate the time that a single run of GA scripts takes
on your computer. Furthermore, the script is used to
@@ -52,7 +63,7 @@ end
%%% RETRIVING DATCOM VARS FROM THE OPTIMIZATION VARIABLES
% x is looping
-[input] = createDatcomStruct(x, settings, input);
+[input] = createDatcomStruct(x, settings, rocket, input);
%%% AERODYNAMICS COEFFICIENT - LAUNCHPAD
% this section is computed just once, then the resulting axial coefficient
@@ -63,48 +74,61 @@ end
input.fltcon.MACH = 0.1;
input.fltcon.ALPHA = [-1, -0.5, 0, 0.5, 1];
input.fltcon.BETA = 0;
-input.fltcon.ALT = settings.z0;
+input.fltcon.ALT = environment.z0;
input.fltcon.NALPHA = length(input.fltcon.ALPHA); % [-] Number of alphas
input.fltcon.NMACH = length(input.fltcon.MACH);
-xcg = settings.xcg;
+xcg = rocket.xCg;
input.refq.XCG = xcg(1);
[coeff] = dissileMatcom(input);
Coeffs0 = coeff(:,:,:,:,:,1);
CA0 = Coeffs0(1, 3);
-%%% LAUNCHPAD DYNAMICS
+%% LAUNCHPAD DYNAMICS
%%% Initial State
X0pad = [0; 0; 0; 0];
%%% Attitude
-Q0 = angleToQuat(settings.PHI, settings.OMEGA, 180*pi/180)';
+Q0 = angleToQuat(environment.phi, environment.omega, 180*pi/180)';
-[Tpad, Ypad] = ode45(@launchPadFreeDyn, [0, 10], X0pad, settings.ode.optionspad,...
- settings, Q0, CA0);
+[Tpad, Ypad] = ode45(@launchPadFreeDyn, [0, 10], X0pad, settings.ode.optLaunchpad,...
+ rocket, environment, Q0, CA0);
launchpad.Tpad = Tpad;
+%%% overrides wind config
+% Wind directions to be looped symmetrical wrt to the launch
+% note: 180 deg is discarded because 0 deg is ever a worst case scenario)
+azimuthValues = [0 45 60 90 120 135 180]*pi/180;
+altitude = -Ypad(end, 3);
+
+wind.altitudes = [0 altitude];
+wind.magnitudeParameters = ones(2,2)*settings.optimization.maxWindMag;
+wind.magnitudeDistribution = ["u", "u"];
+wind.azimuthParameters = zeros(2,2);
+wind.azimuthDistribution = ["u", "u"];
+
+wind.updateAll();
+
%%% COMPUTING THE LAUNCHPAD STABILITY DATA
% both lateral and longitudinal XCP involved
% the launchpad dynamics is used to write the aerodyanmics states
-T = (288.15 - 0.0065*settings.z0); % temperature
+T = (288.15 - 0.0065*environment.z0); % temperature
a = sqrt(T*1.4*287.0531); % sound speed @launchpad
input.fltcon.MACH = round(Ypad(end, 4)/a, 2); % Mach @launchpad
-input.fltcon.ALT = settings.z0; % Altitude @launchpad
+input.fltcon.ALT = environment.z0; % Altitude @launchpad
launchpad.datcom.Mach = input.fltcon.MACH;
-% Wind directions to be looped symmetrical wrt to the launch
-% note: 180 deg is discarded because 0 deg is ever a worst case scenario)
-Az = [0 45 60 90 120 135 180]*pi/180;
-
-%%% pre-allocation
-nAz = length(Az);
+%% pre-allocation
+nAz = length(azimuthValues);
alphaExit = zeros(1, nAz);
betaExit = zeros(nAz, 1);
for i = 1:nAz
- settings.wind.Az = Az(i); % set the direction
- [uw, vw, ww] = windConstGenerator(settings.wind);
+ % updating parameters by creating new wind type at each iteration
+ wind.azimuthParameters = ones(2,2).*azimuthValues(i);
+ wind.updateAll();
+ [uw, vw, ww] = wind.getVels(altitude);
+
inertialWind = [uw, vw, ww];
bodyWind = quatrotate(Q0', inertialWind);
bodyVelocity = [Ypad(end, 4), 0, 0];
@@ -115,6 +139,18 @@ for i = 1:nAz
end
launchpad.betaExit = unique(betaExit);
+% for i = 1:nAz
+% settings.wind.Az = Az(i); % set the direction
+% [uw, vw, ww] = windConstGenerator(settings.wind);
+% inertialWind = [uw, vw, ww];
+% bodyWind = quatrotate(Q0', inertialWind);
+% bodyVelocity = [Ypad(end, 4), 0, 0];
+% Vr = bodyVelocity - bodyWind;
+% ur = Vr(1); vr = Vr(2); wr = Vr(3);
+% alphaExit(i) = round(atand(wr/ur), 1);
+% betaExit(i) = round(atand(vr/ur), 1);
+% end
+
%%% alpha vector imposed to be symmetric and populated
alphaVectorPositive = [abs(alphaExit) 1];
alphaVectorPositive = sort(alphaVectorPositive);
@@ -134,11 +170,10 @@ for i = 1:nAz
launchpad.indexAlpha(i) = indexAlpha-1;
launchpad.indexBeta(i) = indexBeta;
end
-settings.launchpad = launchpad;
-[c, ~] = XCPcheck(x, input, settings);
+[c, ~] = XCPcheck(x, input, data_opt, launchpad, rocket, environment, settings);
%% OptimizationGA
-[apogee] = OptimizationGA(x, input, settings);
+[apogee] = OptimizationGA(x, input, rocket, environment, wind, settings, wrapper);
time = toc;
end
\ No newline at end of file
diff --git a/functions/odeFunctions/launchPadFreeDyn.m b/functions/odeFunctions/launchPadFreeDyn.m
index 29b992c20c68126ee0f7dba7669ce6571a10a211..fa9d7067955dc93c4c69523266c85ed9dcc0b67b 100644
--- a/functions/odeFunctions/launchPadFreeDyn.m
+++ b/functions/odeFunctions/launchPadFreeDyn.m
@@ -1,11 +1,9 @@
-function dY = launchPadFreeDyn(t, Y, rocket, environment, settings, newSettings, Q0, CA)
+function dY = launchPadFreeDyn(t, Y, rocket, environment, Q0, CA)
arguments
t
Y
rocket Rocket
environment Environment
- settings struct
- newSettings Settings
Q0 double
CA double
end
diff --git a/functions/simulations/OLD/quickApogeeOnly.m b/functions/simulations/OLD/quickApogeeOnly.m
index 29d4d6bd7ae75edaad92e44ffe560943758041e5..8f8dda1dc1d948e784c67475c164bd29f571c948 100644
--- a/functions/simulations/OLD/quickApogeeOnly.m
+++ b/functions/simulations/OLD/quickApogeeOnly.m
@@ -1,5 +1,11 @@
-function [apogee, maxAccel] = quickApogeeOnly(settings)
-
+function [apogee, maxAccel] = quickApogeeOnly(rocket, environment, wind, settings, wrapper)
+arguments
+ rocket Rocket
+ environment Environment
+ wind WindCustom
+ settings Settings
+ wrapper DataWrapper
+end
%{
quickApogeeOnly - This function tests the fins simulating the ascent
@@ -31,64 +37,57 @@ SPDX-License-Identifier: GPL-3.0-or-later
%}
%% ERROR CHECKING
-if not(settings.multipleAB) && length(settings.control) > 1
+% if not(settings.multipleAB) && length(rocket.airbrakes.extension) > 1
+% error('To simulate different airbrakes opening, please set to true settings.multipleAB in config.m');
+% end
+%
+% if settings.multipleAB && length(rocket.airbrakes.extension) < 2
+% error('In airbrakes smooth opening simulations, airbrakes configuration must be always greater than 2 (from launch), check config.m')
+% end
+%
+% if settings.multipleAB && length(rocket.airbrakes.extension) > 1 && length(settings.dtControl) < length(rocket.airbrakes.extension) - 2
+% error('In airbrakes smooth opening simulations, AB configuration usage time vector must be at least of length length(pCOntrol)-2, check config.m')
+% end
+
+if not(rocket.airbrakes.enabled) && length(rocket.airbrakes.extension) > 1
error('To simulate different airbrakes opening, please set to true settings.multipleAB in config.m');
end
-if settings.multipleAB && length(settings.control) < 2
- error('In airbrakes smooth opening simulations, airbrakes configuration must be always greater than 2 (from launch), check config.m')
-end
+% if rocket.airbrakes.enabled && length(rocket.airbrakes.extension) < 2
+% error('In airbrakes smooth opening simulations, airbrakes configuration must be always greater than 2 (from launch), check config.m')
+% end
-if settings.multipleAB && length(settings.control) > 1 && length(settings.dtControl) < length(settings.control) - 2
+if rocket.airbrakes.enabled && length(rocket.airbrakes.extension) > 1 && length(rocket.airbrakes.deltaTime) < length(rocket.airbrakes.extension) - 2
error('In airbrakes smooth opening simulations, AB configuration usage time vector must be at least of length length(pCOntrol)-2, check config.m')
end
-if settings.wind.variable
- error('Wind variable model available just in stoch simulations');
-end
-
-%% STARTING CONDITIONS
+% if settings.wind.variable
+% error('Wind variable model available just in stoch simulations');
+% end
-% Attitude
-Q0 = angleToQuat(settings.PHI, settings.OMEGA, 180*pi/180)';
+%% DATA PREPARATION
+options = odeset(settings.ode.optAscent, 'OutputFcn', @storeData);
+tf = settings.ode.finalTime;
-%% WIND GENERATION
+%% INITIAL CONDITIONS
+%%% Attitude
+Q0 = angleToQuat(environment.phi, environment.omega, 180*pi/180)';
-[uw, vw, ww, ~] = windConstGenerator(settings.wind);
-settings.constWind = [uw, vw, ww];
-tf = settings.ode.finalTime;
+%%% State
+X0 = [0; 0; 0];
+V0 = [0; 0; 0];
+W0 = [0; 0; 0];
-%% ASCENT
-X0 = [0 0 0]';
-V0 = [0 0 0]';
-W0 = [0 0 0]';
+Y0a = [X0; V0; W0; Q0; 0];
-Y0a = [X0; V0; W0; Q0];
+%% WIND GENERATION
+% [uw, vw, ww, ~] = windConstGenerator(settings.wind);
+% settings.constWind = [uw, vw, ww];
-if settings.multipleAB
-
- nAB = length(settings.control);
- t0 = 0;
- Ta = [];
- Ya = [];
-
- settings.delayControl = delayControl(settings.control, settings);
-
- for iAB = 1:nAB
-
- [Tab, Yab] = ode113(@ascentMultipleAB, [t0, tf], Y0a, settings.ode.optionsascMultipleAB, t0, iAB, settings);
-
- % update variables and the state
- t0 = Tab(end);
- Y0a = Yab(end, :);
- Ta = [Ta; Tab];
- Ya = [Ya; Yab];
-
- end
-
-else
- [Ta, Ya] = ode113(@ascent, [0, tf], Y0a, settings.ode.optionsasc1, settings);
-end
+%% ASCENT
+t0 = 0;
+[Ta, Ya] = ode113(@ballistic, [t0, tf], Y0a, options, ...
+ rocket, environment, wind, [], 0, wrapper);
%% CALCULATE OUTPUT QUANTITIES
apogee = -Ya(end, 3);