From 58afe444cc83619e9b9c3f257a60b6025ba22a31 Mon Sep 17 00:00:00 2001
From: Mauco03 <marco.gaibotti@skywarder.eu>
Date: Mon, 22 Apr 2024 09:29:07 +0200
Subject: [PATCH] [refactoring-ode][classes, missions] Bugfix

- Rocket will not load coefficients if matrices do not exist
- Fixed a bug in Wind where it would not extrapolate altitudes
- Added deg2Rad conversion in Environment
- Updated configs/overrides
---
 classes/Rocket.m                              | 14 ++++++++-----
 classes/components/Environment.m              |  4 ++++
 classes/components/WindCustom.m               |  4 ++--
 functions/miscellaneous/createDissileInput.m  |  7 +++++++
 functions/odeFunctions/ballistic.m            | 20 +++++++++----------
 .../config/rocketConfig.m                     |  9 ++++++++-
 .../data/CAinterpCoeffs.mat                   |  2 +-
 .../data/rocketStress.mat                     |  4 ++--
 .../data/rocketStressMain.mat                 |  4 ++--
 .../config/environmentConfig.m                |  6 +++---
 .../config/rocketConfig.m                     | 11 ++++++----
 .../config/windConfig.m                       |  2 +-
 12 files changed, 56 insertions(+), 31 deletions(-)

diff --git a/classes/Rocket.m b/classes/Rocket.m
index 2ad97d1..8968662 100644
--- a/classes/Rocket.m
+++ b/classes/Rocket.m
@@ -125,7 +125,6 @@ classdef Rocket < Component
             if ~isempty(obj.xCg)
                 return;
             end
-            % Motor mass is 5th element of "bays"
             obj.xCg = (obj.xCgNoMotor*obj.massNoMotor + ...
                 obj.motor.mass.*(obj.motor.xCg+obj.lengthCenterNoMot))./ ...
                 obj.mass;
@@ -188,10 +187,15 @@ classdef Rocket < Component
                 return;
             end
             varNames = {'total', 'geometry', 'state', 'finsCN'};
-            dataCoeffs = load(fullfile(obj.mission.dataPath, 'aeroCoefficients.mat'), ...
-                varNames{:});
-            dataCoeffsHighAOA = load(fullfile(obj.mission.dataPath, 'aeroCoefficientsHighAOA.mat'), ...
-                varNames{:});
+            aeroPath = fullfile(obj.mission.dataPath, 'aeroCoefficients.mat');
+            aeroHighAOAPath = fullfile(obj.mission.dataPath, 'aeroCoefficientsHighAOA.mat');
+
+            if ~(exist(aeroPath, 'file') && exist(aeroHighAOAPath,'file'))
+                return;
+            end
+
+            dataCoeffs = load(aeroPath, varNames{:});
+            dataCoeffsHighAOA = load(aeroHighAOAPath, varNames{:});
             
             if isfield(dataCoeffs, "finsCN") && isfield(dataCoeffsHighAOA, "finsCN")
                 obj.coefficients.finsCN = dataCoeffs.finsCN;
diff --git a/classes/components/Environment.m b/classes/components/Environment.m
index 25e3d77..2f76a64 100644
--- a/classes/components/Environment.m
+++ b/classes/components/Environment.m
@@ -47,6 +47,10 @@ classdef Environment < Component
             end
             obj@Component(mission, varIn);
             obj.motor = motor;
+
+            %% Update angles
+            obj.omega = deg2rad(obj.omega);
+            obj.phi = deg2rad(obj.phi);
         end
     end
 
diff --git a/classes/components/WindCustom.m b/classes/components/WindCustom.m
index 69aa45f..af8f0f8 100644
--- a/classes/components/WindCustom.m
+++ b/classes/components/WindCustom.m
@@ -48,8 +48,8 @@ classdef WindCustom < Component
                 vw = round( - obj.magnitude * sin(obj.azimuth), 6);
                 ww = 0;
             else
-                mag = interpLinear(obj.altitudes, obj.magnitude, z);
-                az = interpLinear(obj.altitudes, obj.azimuth, z);
+                mag = interpLinear(obj.altitudes, obj.magnitude, z, true);
+                az = interpLinear(obj.altitudes, obj.azimuth, z, true);
                 
                 uw = round( - mag * cos(az), 6);
                 vw = round( - mag * sin(az), 6);
diff --git a/functions/miscellaneous/createDissileInput.m b/functions/miscellaneous/createDissileInput.m
index a5216fe..ebc2818 100644
--- a/functions/miscellaneous/createDissileInput.m
+++ b/functions/miscellaneous/createDissileInput.m
@@ -98,6 +98,13 @@ if ~isempty(rocket.pitot.initialConeLength)
     input.axibod.PITOTCDB = rocket.pitot.finalConeDiameter;                                      % [m] Pitot final conic section diameter
     input.axibod.PITOTL = rocket.pitot.length;                                                   % [m] Pitot tube length
     input.axibod.PITOTD = rocket.PITOTD;                                                         % [m] Pitot tube diameter
+else
+    input.axibod.PITOTCLI = 0;                                                                  % [m] Pitot initial conic section length
+    input.axibod.PITOTCLB = 0;                                                                  % [m] Pitot final conic section length
+    input.axibod.PITOTCDI = 0;                                                                  % [m] Pitot initial conic section diameter
+    input.axibod.PITOTCDB = 0;                                                                  % [m] Pitot final conic section diameter
+    input.axibod.PITOTL = 0;                                                                    % [m] Pitot tube length
+    input.axibod.PITOTD = 0;                                                                    % [m] Pitot tube diameter
 end
 
 %%% centerbody
diff --git a/functions/odeFunctions/ballistic.m b/functions/odeFunctions/ballistic.m
index b435a4a..3baa2f2 100644
--- a/functions/odeFunctions/ballistic.m
+++ b/functions/odeFunctions/ballistic.m
@@ -78,7 +78,7 @@ g = environment.g0/(1 + (-z*1e-3/6371))^2;              % [N/kg]  module of grav
 tb = rocket.motor.cutoffTime;                           % [s]     Burning Time
 local = [environment.z0, environment.temperature, ...   % vector containing inputs for atmosphereData
     environment.pressure, environment.rho];
-omega = environment.omega*pi/180;
+omega = environment.omega;
 
 %% QUATERION ATTITUDE
 Q = [q0 q1 q2 q3];
@@ -213,8 +213,8 @@ end
 %% RAMP / FREE FLIGHT
 if -z < environment.effectiveRampLength*sin(omega)      % No torque on the launchpad
     Fg = m*g*sin(omega);                % [N] force due to the gravity
-    X = 0.5*rho*V_norm^2*S*CA;
-    F = -Fg +T -X;
+    fX = 0.5*rho*V_norm^2*S*CA;
+    F = -Fg +T -fX;
     du = F/m;
 
     dv = 0;
@@ -225,8 +225,8 @@ if -z < environment.effectiveRampLength*sin(omega)      % No torque on the launc
 
     alpha_value = NaN;
     beta_value = NaN;
-    Y = 0;
-    Z = 0;
+    fY = 0;
+    fZ = 0;
     XCPlon = NaN;
     XCPlat = NaN;
 
@@ -241,11 +241,11 @@ else
     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*CY;                      % [N] y-body component of the aerodynamics force
-    Z = qdyn*S*CN;                      % [N] z-body component of the aerodynamics force
+    fX = qdyn*S*CA;                      % [N] x-body component of the aerodynamics force
+    fY = qdyn*S*CY;                      % [N] y-body component of the aerodynamics force
+    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 + [-X+T, Y, -Z]';             % [N] total forces vector
+    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
@@ -312,7 +312,7 @@ if nargout == 2
     
     parout.velocities = Vels;
     
-    parout.forces.AeroDyn_Forces = [X, Y, Z];
+    parout.forces.AeroDyn_Forces = [fX, fY, fZ];
     parout.forces.T = T;
     
     parout.air.rho = rho;
diff --git a/missions/2024_Lyra_Portugal_October/config/rocketConfig.m b/missions/2024_Lyra_Portugal_October/config/rocketConfig.m
index c039bc4..f2274b5 100644
--- a/missions/2024_Lyra_Portugal_October/config/rocketConfig.m
+++ b/missions/2024_Lyra_Portugal_October/config/rocketConfig.m
@@ -14,7 +14,7 @@ rocket.lengthCenterNoMot = 1.778;                          % [m]      OVERRIDE C
 %% PLD - Includes Payload + Nose
 payload = Payload();
 
-payload.length = 0.32;                                     % [m] Total bay length
+payload.length = 0.6;                                      % [m] Total bay length
 payload.mass = 3.38512;                                    % [kg] Total bay mass
 payload.inertia = [];                                      % [kg*m^2] Total bay inertia (Body reference)
 payload.xCg = [];                                          % [m] Cg relative to bay upper side
@@ -49,12 +49,16 @@ airbrakes.mass = [];                                  % [kg] Total bay mass
 airbrakes.inertia = [];                               % [kg*m^2] Total bay inertia (Body reference)
 airbrakes.xCg = [];                                   % [m] Cg relative to bay upper side
 
+airbrakes.multipleAB = false;                              % If true, multiple and smooth airbrakes opening will be simulated
+airbrakes.opening = [0 1 0.5];                                 % aerobrakes, 1-2-3 for 0%, 50% or 100% opened
+airbrakes.deltaTime = [5 5 5];                                % aerobrakes, configurations usage time
 airbrakes.multipleAB = false;                            % If true, multiple and smooth airbrakes opening will be simulated
 airbrakes.opening = [1];                                 % aerobrakes, 1-2-3 for 0%, 50% or 100% opened
 airbrakes.deltaTime = [];                                % aerobrakes, configurations usage time
 
 airbrakes.n = 3;                                           % [-] number of brakes
 airbrakes.height = linspace(0, 0.0363, 3);                 % [m] Block airbrakes opening coordinate ( First entry must be 0! )
+airbrakes.angles = [0, 45, 68];
 airbrakes.width = 0.1002754821;                            % [m] brakes width (normal)
 airbrakes.thickness = 0.008;                               % [m] brakes thickness
 airbrakes.xDistance = 1.517;                               % [m] axial position from nosecone base
@@ -66,6 +70,9 @@ airbrakes.servoOmega = 150*pi/180;                         % [rad/s] Servo-motor
 %% MOTOR
 motor = Motor();
 
+motor.name = 'HRE_ARM_30_inj-T03T03';                       % [-] Motor name
+motor.cutoffTime = [];                                     % [s] OVERRIDE Cutoff time
+motor.ignitionTransient = 0.4;                             % [s] Ignition transient
 motor.name = 'HRE_ARM_30_inj-T03T03';                      % [-] Motor name
 motor.cutoffTime = inf;                                    % [s] Cutoff time
 motor.ignitionTransient = 0.3;                             % [s] Ignition transient
diff --git a/missions/2024_Lyra_Portugal_October/data/CAinterpCoeffs.mat b/missions/2024_Lyra_Portugal_October/data/CAinterpCoeffs.mat
index 5302536..97df552 100644
--- a/missions/2024_Lyra_Portugal_October/data/CAinterpCoeffs.mat
+++ b/missions/2024_Lyra_Portugal_October/data/CAinterpCoeffs.mat
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:d5d212e8d6501ea7d1e4572c84de2a59887f7a87e965dcbcf189fd37142fffdd
+oid sha256:33c07484db34ceb100aaeb9ea817d40471a25657e4b63395ae3d40244f44691f
 size 457
diff --git a/missions/2024_Lyra_Portugal_October/data/rocketStress.mat b/missions/2024_Lyra_Portugal_October/data/rocketStress.mat
index 7f20d2d..2c82625 100644
--- a/missions/2024_Lyra_Portugal_October/data/rocketStress.mat
+++ b/missions/2024_Lyra_Portugal_October/data/rocketStress.mat
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:a2af604b52fb12df3ffc96825ccb04e4ec3fa9bce8a0f70c3aa2cb2ae4100177
-size 1138617
+oid sha256:284b485d8f2c3fc1ac51dabdff1966fddf8961bf6cfcbc6fac6e0d0cc57da511
+size 636771
diff --git a/missions/2024_Lyra_Portugal_October/data/rocketStressMain.mat b/missions/2024_Lyra_Portugal_October/data/rocketStressMain.mat
index 9e62ba7..b33c119 100644
--- a/missions/2024_Lyra_Portugal_October/data/rocketStressMain.mat
+++ b/missions/2024_Lyra_Portugal_October/data/rocketStressMain.mat
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:50626a92cd9ddc4603b035c1b472966ef0608558fcf1e9073555f9cc2a19b0c8
-size 420614
+oid sha256:1ba65a59b508f8e02c202b1145f3af61e6455b64a80de60039625bc024c9e359
+size 214027
diff --git a/missions/2024_Lyra_Roccaraso_September/config/environmentConfig.m b/missions/2024_Lyra_Roccaraso_September/config/environmentConfig.m
index e90dd15..8f05f5a 100644
--- a/missions/2024_Lyra_Roccaraso_September/config/environmentConfig.m
+++ b/missions/2024_Lyra_Roccaraso_September/config/environmentConfig.m
@@ -9,9 +9,9 @@ environment.lon0 = 14.0546408;          % [deg] Launchpad longitude
 environment.z0 = 1414;                  % [m] Launchpad Altitude
 environment.omega = 85;                 % [deg] Launchpad Elevation
 environment.phi = 0;                    % [deg] Launchpad Azimuth
-environment.pin1Length = 000;           % [m] Distance from the upper pin to the upper tank cap
-environment.pin2Length = 000;           % [m] Distance from the lower pin to the lower tank cap  
-environment.rampLength = 000;           % [m] Total launchpad length 
+environment.pin1Length = 0.5603;        % [m] Distance from the upper pin to the upper tank cap
+environment.pin2Length = 0.2055;        % [m] Distance from the lower pin to the lower tank cap  
+environment.rampLength = 10;            % [m] Total launchpad length 
 
 environment.temperature = [];           % [deg] Ground temperature correction
 environment.pressure = [];              % [Pa] Ground pressure correction
diff --git a/missions/2024_Lyra_Roccaraso_September/config/rocketConfig.m b/missions/2024_Lyra_Roccaraso_September/config/rocketConfig.m
index c619623..c0e80a3 100644
--- a/missions/2024_Lyra_Roccaraso_September/config/rocketConfig.m
+++ b/missions/2024_Lyra_Roccaraso_September/config/rocketConfig.m
@@ -6,11 +6,14 @@ rocket = Rocket();
 rocket.diameter = 0.15;                                    % [m]        Rocket diameter
 rocket.length = [];                                        % [m]        OVERRIDE total length
 rocket.mass = [];                                          % [kg]       OVERRIDE total mass
-rocket.massNoMotor = [];                                   % [kg]       OVERRIDE mass without motor
+rocket.massNoMotor = [15.6348];                            % [kg]       OVERRIDE mass without motor
 rocket.inertia = [];                                       % [kg*m^2]   OVERRIDE total inertia - Axibody reference
-rocket.inertiaNoMotor = [];                                % [kg*m^2]   OVERRIDE inertia without motor
+rocket.inertiaNoMotor = [
+    0.06535397;
+    12.07664659;
+    12.07701314];                                          % [kg*m^2]   OVERRIDE inertia without motor
 rocket.xCg = [];                                           % [m]        OVERRIDE total xCg
-rocket.xCgNoMotor = [];                                    % [m]        OVERRIDE xCg without motor
+rocket.xCgNoMotor = [1.28];                                    % [m]        OVERRIDE xCg without motor
 
 %% PLD - Includes Payload + Nose
 payload = Payload();
@@ -102,7 +105,7 @@ airbrakes.servoOmega = 150*pi/180;                         % [rad/s] Servo-motor
 %% MOTOR
 motor = Motor();
 
-motor.name = 'HRE_FURIA-Rv2-T04T03';                       % [-] Motor name
+motor.name = 'HRE_ARM_Rocc_U06-T03T03';                       % [-] Motor name
 motor.cutoffTime = [];                                     % [s] OVERRIDE Cutoff time
 motor.ignitionTransient = 0.4;                             % [s] Ignition transient
 motor.cutoffTransient = 0.3;                               % [s] Cut-off transient
diff --git a/missions/2024_Lyra_Roccaraso_September/config/windConfig.m b/missions/2024_Lyra_Roccaraso_September/config/windConfig.m
index 332c74e..5392f93 100644
--- a/missions/2024_Lyra_Roccaraso_September/config/windConfig.m
+++ b/missions/2024_Lyra_Roccaraso_September/config/windConfig.m
@@ -16,7 +16,7 @@ windCustom.magnitudeDistribution = ["g", "u", "u"];     % [-] Distribution type:
 windCustom.magnitudeParameters = [7 2 10;               % [m/s] Distribution parameters: "u" - [min; max], "g" - [mu; sigma]
                             0.5 9 20];
 windCustom.azimuthDistribution = ["u", "u", "u"];       % [-] Distribution type: "u" - uniform, "g" - gaussian
-windCustom.azimuthParameters = 0*pi/180 * ones(2,3);    % [deg] Distribution parameters: "u" - [min; max], "g" - [mu; sigma]
+windCustom.azimuthParameters = 90*pi/180 * ones(2,3);    % [deg] Distribution parameters: "u" - [min; max], "g" - [mu; sigma]
 
 %% MATLAB WIND MODEL
 
-- 
GitLab