diff --git a/functions/odeFunctions/ballistic.m b/functions/odeFunctions/ballistic.m
index 60f3e92e05341f06f7585af86247e9b90496e940..0194eff2564ba79ba00e30e47181915a47a14891 100644
--- a/functions/odeFunctions/ballistic.m
+++ b/functions/odeFunctions/ballistic.m
@@ -10,7 +10,7 @@ arguments
wrapper = [] %DataWrapper = DataWrapper.empty
end
% ballistic - ode function of the 6DOF Rigid Rocket Model
-%
+%
% INPUTS:
% - t, double [1, 1] integration time [s];
% - Y, double [14, 1] state vector [ x y z | u v w | p q r | q0 q1 q2 q3 | angle]:
@@ -24,16 +24,16 @@ end
% - control, struct (angleRef, timeChangeRef, partialTime) GNC parameters:
% * angleRef, airbrake servo angle reference
% * timeChangeRef, time to change airbrake configuration
-% * partialTime,
+% * partialTime,
% - wrapper, handle to export data
-%
+%
% OUTPUTS:
% - dY, double [14, 1] state derivatives
% - parout, struct, additional flight qunatities
%
% Copyright © 2021, Skyward Experimental Rocketry, AFD, GNC departments
% All rights reserved
-%
+%
% SPDX-License-Identifier: GPL-3.0-or-later
%% RECALLING STATES
@@ -183,14 +183,14 @@ Cn = (Cn0 + Cnb*(beta - beta0));
if abs(alpha) <= pi/180
XCPlon = Cma/CNA;
else
- XCPlon = Cm/CN;
+ XCPlon = Cm/CN;
end
-
+
% XCPlat = Cn/CY;
if abs(beta) <= pi/180
- XCPlat = Cnb/CYB;
+ XCPlat = Cnb/CYB;
else
- XCPlat = Cn/CY;
+ XCPlat = Cn/CY;
end
% if Cn == 0 && CY == 0
@@ -221,10 +221,10 @@ if altitude < environment.effectiveRampAltitude % No
if T < Fg % No velocity untill T = Fg
du = 0;
end
-
- XCPtot = NaN;
+
+ XCPtot = NaN;
else
-%% FORCES
+ %% FORCES
% first computed in the body-frame reference system
qdyn = 0.5*rho*velsNorm^2; % [Pa] dynamics pressure
qdynL_V = 0.5*rho*velsNorm*S*C;
@@ -234,12 +234,12 @@ else
fZ = qdyn*S*CN; % [N] z-body component of the aerodynamics force
Fg = dcm*[0; 0; m*g]; % [N] force due to the gravity in body frame
F = Fg + [-fX+T, fY, -fZ]'; % [N] total forces vector
-
+
%-----------------------------------------------------
%F = Fg + [-X+T*cos(chi), Y+T*sin(chi), -Z]'; % [N] total forces vector
%-----------------------------------------------------
-%% STATE DERIVATIVES
+ %% STATE DERIVATIVES
% velocity
du = F(1)/m - q*w + r*v;
dv = F(2)/m - r*u + p*w;
@@ -252,37 +252,37 @@ else
dr = (Ixx - Iyy)/Izz*p*q + qdynL_V/Izz*(velsNorm*Cn + (Cnr*r+Cnp*p)*C/2)...
- Izzdot*r/Izz;
- % Compute the aerodynamic roll angle
- [~, phi] = getAlphaPhi(alpha, beta);
-
+ % Compute the aerodynamic roll angle
+ [~, phi] = getAlphaPhi(alpha, beta);
+
% Aerodynamic-force coefficient in the alpha-total plane
- CFaTot = sin(phi)*CY + cos(phi)*(-CN);
+ CFaTot = sin(phi)*CY + cos(phi)*(-CN);
% Aerodynanic-moment coefficient in the alpha-total plane
CMaTot = cos(phi)*Cm - sin(phi)*Cn;
if CFaTot ~= 0
- XCPtot = CMaTot/CFaTot;
+ XCPtot = CMaTot/CFaTot;
else
- XCPtot = XCPlon;
+ XCPtot = XCPlon;
end
dAngle = 0;
if isControlActive && ~isAngleSaturated
- angleRef = control.angleRef;
+ angleRefVec = control.angleRef;
if ~rocket.airbrakes.identification
% set velocity of servo (air brakes)
- if length(angleRef) == 2 % for the recallOdeFunction
- if theta < control.timeChangeRef
- angleRef = angleRef(1);
+ if length(angleRefVec) == 2 % for the recallOdeFunction
+ if t < control.timeChangeRef
+ angleRef = angleRefVec(1);
else
- angleRef = angleRef(2);
+ angleRef = angleRefVec(2);
end
else
- [~,idxAirbrake] = min(control.partialTime - theta);
- angleRef = angleRef(idxAirbrake); % don't delete this unless you change how the recall ode works.
+ [~,idxAirbrake] = min(control.partialTime - t);
+ angleRef = angleRefVec(idxAirbrake); % don't delete this unless you change how the recall ode works.
end
else
- [~,idxAirbrake] = min(abs(theta-angleRef(:,1))); % if we are trying to identify we need to have the same input of the flight
- angleRef = angleRef(idxAirbrake,2);
+ [~,idxAirbrake] = min(abs(t-angleRefVec(:,1))); % if we are trying to identify we need to have the same input of the flight
+ angleRef = angleRefVec(idxAirbrake,2);
end
dAngle = (angleRef-angle)/rocket.airbrakes.servoTau;
@@ -293,9 +293,9 @@ else
end
% Quaternions
OM = [ 0 -p -q -r ;
- p 0 r -q ;
- q -r 0 p ;
- r q -p 0 ];
+ p 0 r -q ;
+ q -r 0 p ;
+ r q -p 0 ];
dQQ = 1/2*OM*Q';
@@ -316,26 +316,26 @@ if nargout == 2 || ~isempty(wrapper)
parout.interp.beta = betaOut;
parout.interp.alt = altitude;
parout.interp.mass = m;
- parout.interp.inertias = [Ixx; Iyy; Izz];
+ parout.interp.inertias = [Ixx; Iyy; Izz];
parout.wind.NED = [uw; vw; ww];
parout.wind.body = windVels;
-
+
parout.rotations.dcm = dcm;
-
+
parout.velocities = vels;
-
+
parout.forces.aero = [fX; fY; fZ];
parout.forces.T = T;
-
+
parout.air.rho = rho;
parout.air.P = P;
-
+
parout.accelerations.body = [du; dv; dw];
parout.accelerations.angular = [dp; dq; dr];
parout.accelerometer.body_acc = ([-fX+T, fY, -fZ]')/m;
-
+
parout.coeff.CA = CA;
parout.coeff.CYB = CYB;
parout.coeff.CNA = CNA;
@@ -351,9 +351,9 @@ if nargout == 2 || ~isempty(wrapper)
parout.coeff.XCPlon = XCPlon;
parout.coeff.XCPlat = XCPlat;
parout.coeff.XCPtot = XCPtot;
-
+
parout.uncertanty = [];
if ~isempty(wrapper)
- wrapper.setCache(parout);
+ wrapper.setCache(parout);
end
end
\ No newline at end of file