diff --git a/classes/@Rocket/Rocket.m b/classes/@Rocket/Rocket.m
index a7796d40e0e63d2cb72130d62aa8ca29810c5b12..f050235147b644d3a07242764a10445361213e54 100644
--- a/classes/@Rocket/Rocket.m
+++ b/classes/@Rocket/Rocket.m
@@ -52,7 +52,8 @@ classdef Rocket < Config
         parachutes                  Para            % [-]       (nParachutes, nStages) Parachutes onboard
 
         dynamicDerivatives  (1, 1)  logical         % [-]       True if dynamic derivatives will be loaded
-        
+        highAOACoefficients (1, 1)  logical         % [-]       True if highAOACoefficients are loaded
+
         coefficients                Coefficient     % [-]       Aerodynamic coefficients
         coefficientsHighAOA         Coefficient     % [-]       Aerodynamic coefficients at high angle of attack
     end
@@ -149,13 +150,16 @@ classdef Rocket < Config
                     fullfile(mission.dataPath, 'aeroCoefficients.mat'), ...
                     coeffName);
 
-                obj.coefficientsHighAOA = Coefficient( ...
-                    fullfile(mission.dataPath, 'aeroCoefficientsHighAOA.mat'), ...
-                    coeffName);
+                if obj.highAOACoefficients
+                    obj.coefficientsHighAOA = Coefficient( ...
+                        fullfile(mission.dataPath, 'aeroCoefficientsHighAOA.mat'), ...
+                        coeffName);
+                end
 
                 answer = '';
                 
-                if isempty(obj.coefficients.static) || isempty(obj.coefficientsHighAOA.static)
+                if isempty(obj.coefficients.static) || ...
+                        (obj.highAOACoefficients && isempty(obj.coefficientsHighAOA.static))
                     answer = questdlg(['Coefficient matrices not found. ' ...
                         'Do you want to create new matrices?']);
                 elseif options.checkGeometry
@@ -194,7 +198,7 @@ classdef Rocket < Config
                         parserPath = fullfile(mission.msaPath, 'autoMatricesProtub');
                         addpath(genpath(parserPath));
                         [obj.coefficients, obj.coefficientsHighAOA] = ...
-                            mainAutoMatProtub(obj);
+                            mainAutoMatProtub(obj, 'computeHighAOA', obj.highAOACoefficients);
                     case 'Cancel'
                         error('Rocket creation aborted')
                     otherwise
diff --git a/classes/enums/CoeffType.m b/classes/enums/CoeffType.m
new file mode 100644
index 0000000000000000000000000000000000000000..85e051cba449a53e81cbc4661ab2025987e83284
--- /dev/null
+++ b/classes/enums/CoeffType.m
@@ -0,0 +1,7 @@
+classdef CoeffType
+    enumeration
+        AlphaTotPhi
+        AlphaBeta
+    end
+end
+
diff --git a/classes/bays/MotorType.m b/classes/enums/MotorType.m
similarity index 100%
rename from classes/bays/MotorType.m
rename to classes/enums/MotorType.m
diff --git a/classes/misc/Coefficient.asv b/classes/misc/Coefficient.asv
new file mode 100644
index 0000000000000000000000000000000000000000..924dc074dcde46149da25f7dbe5052a1e5eecfc7
--- /dev/null
+++ b/classes/misc/Coefficient.asv
@@ -0,0 +1,322 @@
+classdef Coefficient
+    % Coefficient: Manages the coefficients associated with a rocket
+    %
+    %   Constructor:
+    %       - Coefficient: Creates an instance of the Coefficient class.
+    %           Loaded config: -
+    %           Loaded data: aeroCoefficients.mat
+    %           Arguments:
+    %               - filePath: char, path to coefficient file to load
+    %               - name: (optional) coefficient name. Used when dynamic
+    %               derivatives are also needed
+    %
+    %   Coefficients are retrieved from DATCOM in the following format:
+    %       Coefficient matrix: double (15, nAlpha, nMach, nBeta, nAlt, nAbk, nXcg) 
+    %       Where the 15 coefficients are 
+    %           - CA        Axial force 
+    %           - CYB       Side force derivative wrt beta
+    %           - CY0       Side force
+    %           - CNA       Normal force derivative wrt alpha
+    %           - CN0       Normal force
+    %           - Cl        Rolling moment
+    %           - Clp       Rolling moment derivative wrt roll rate
+    %           - Cma       Pitching moment derivative wrt alpha
+    %           - Cm0       Pitching moment
+    %           - Cmad      Pitching moment derivative wrt alpha dot
+    %           - Cmq       Pitching moment derivative wrt pitch rate
+    %           - Cnb       Yawing moment derivative wrt beta
+    %           - Cn0       Yawing moment
+    %           - Cnr       Yawing moment derivative wrt yaw rate
+    %           - Cnp       Yawing moment derivative wrt roll rate
+    %
+    %   Coefficients are then stored in the following format:
+    %       static:  double (6, nAlpha, nMach, nBeta, nAlt, nAbk)
+    %       Where the 6 coefficients are: 
+    %           - [CA, CY, CN, Cl, Cm, Cn]
+    %       dynamic: double (5, nAlpha, nMach, nBeta, nAlt, nAbk, nXcg)
+    %           - [Clp, Cmad, Cmq, Cnr, Cnp]
+    %
+    %   NOTE: When dynamic derivatives are not loaded, the matrix will be
+    %   6-D, as there will be no need to interpolate wrt to the xcg: moment
+    %   transport formulas will be used instead
+    %
+    %   NOTE: Coefficients in an interpolants are stored on the last
+    %   dimension, instead of the first, to improve performance
+
+    properties(Dependent)
+        total                       double  % Coefficients stored in DATCOM format
+        geometry                    struct  % Reference geometry
+        state                       struct  % Flight envelope
+    end
+
+    properties
+        finsCN                      double  % Fins-only CN
+    end
+
+    properties(Dependent, SetAccess = private)
+        static                      double  % Static coefficients: [CA, CY, CN, Cl, Cm, Cn]
+        dynamic                     double  % Dynamic derivatives: [Clp, Cmad, Cmq, Cnr, Cnp]
+    end
+
+    properties(SetAccess = private)
+        type                        CoeffType   % Wheter based on alphaTot-phi or alpha-beta
+        symmetry            (1, 1)  logical     % Apply axial simmetry: valid for alphaTot-phi
+
+        isReady             (1, 1)  logical     % Whether all coefficients are loaded and of the right size
+        isDynamic           (1, 1)  logical     % Whether to load dynamic derivatives. Adds XCG dependece
+    end
+
+    properties(Access = private)
+        TOTAL                               % Cached variable for total
+        GEOMETRY                            % Cached variable for geometry
+        STATE                               % Cached variable for state
+    end
+
+    properties(Access = private)
+        staticInterpolant   (1, 1)  griddedInterpolant
+        dynamicInterpolant  (1, 1)  griddedInterpolant
+    end
+
+    methods
+        function obj = Coefficient(filePath, name)
+            arguments
+                filePath    (1, :)  char = ''
+                name        (1, :)  char = 'generic'
+            end
+
+            if isempty(filePath), return; end
+            
+            obj = obj.loadData(filePath, name);
+            obj.isReady = obj.checkProperties();
+            obj         = obj.updateInterpolants();
+        end
+
+        function [coefficients] = get(obj, alpha, mach, beta, altitude, airbakes, xcg)
+            % GET: Retrieve aerodynamic coefficients for specified flight conditions
+            %
+            %   This method interpolates the aerodynamic coefficients based on
+            %   the provided flight conditions and adjusts them for the given
+            %   center of gravity (xcg).
+            %
+            %   Arguments:
+            %       - alpha     double, angle of attack             [rad]
+            %       - mach      double, Mach number                 [-]
+            %       - beta      double, sideslip angle              [rad]
+            %       - altitude  double, altitude                    [m]
+            %       - airbakes  double, airbrake deployment         [-]
+            %       - xcg       double, center of gravity position  [m]
+            %
+            %   Returns:
+            %       - coefficients: double (11, 1), aerodynamic coefficients
+            %           [CA, CY, CN, Cl, Cm, Cn, Clp, Cmad, Cmq, Cnr, Cnp]
+            %
+            %   NOTE: Static coefficients are interpolated and adjusted for
+            %   the specified xcg. Dynamic derivatives are included only if
+            %   loadDynamic is true.
+            
+            if ~obj.isReady, error('Cannot interpolate coefficients: check that all dimensions match'); end
+
+            % Interpolating coefficients
+            coefficients = zeros(11, 1);
+            if obj.type == CoeffType.AlphaBeta
+                coefficients(1:6) = obj.staticInterpolant(alpha, mach, beta, altitude, airbakes);
+            else
+                [alphaTot, phi] = getAlphaPhi(alpha, beta);
+                phi = wrapTo2Pi(phi);
+
+                finAngle = 2*pi / obj.GEOMETRY.nPanel;                           % Angle between two fin panels
+                n = floor(phi / finAngle);
+                deltaPhi = n*finAngle;  
+                psi = phi - deltaPhi;                                            % Angle wrapped to finAngle
+                
+                coefficients(1:6) = obj.staticInterpolant(alphaTot, mach, psi, altitude, airbakes); % Get coeffs in limited range
+
+                if n > 0 % If necessary, perform rotation on different frame
+                    R = [cos(deltaPhi), -sin(deltaPhi);
+                        sin(deltaPhi), cos(deltaPhi)];
+    
+                    coefficients([2, 3]) = R  * coefficients([2, 3]);
+                    coefficients([5, 6]) = R' * coefficients([5, 6]);
+                end
+            end
+
+            % Transporting static coefficients to new xcg
+            % C_B = C_A + d / c * [0; -CN; CY]_A 
+            d = xcg - obj.GEOMETRY.xcg(1);
+            l = obj.GEOMETRY.diameter;
+
+            CY = coefficients(2); CN = coefficients(3);
+            forceCoeffs = [0; CN; CY];
+            
+            coefficients(4:6) = coefficients(4:6) + d/l * forceCoeffs;
+
+            if ~obj.isDynamic, return; end
+            if obj.type == CoeffType.AlphaTotPhi
+                error('Alpha-Phi coefficients do not support dynamic derivatives yet')
+            end
+            coefficients(7:11) = obj.dynamicInterpolant(alpha, mach, beta, altitude, airbakes, xcg);
+        end
+    end
+
+    methods % Getters
+        function value = get.total(obj)
+            value = obj.TOTAL;
+        end
+
+        function value = get.static(obj)
+            if isempty(obj.TOTAL), value = []; return; end
+            value = obj.TOTAL([1, 3, 5, 6, 9, 13], :, :, :, :, :, 1);
+        end
+
+        function value = get.dynamic(obj)
+            if ~obj.isDynamic || isempty(obj.TOTAL), value = []; return; end
+            value = obj.TOTAL([7, 10, 11, 14, 15], :, :, :, :, :, :);
+        end
+
+        function value = get.geometry(obj)
+            value = obj.GEOMETRY;
+        end
+
+        function value = get.state(obj)
+            value = obj.STATE;
+        end
+    end
+
+    methods % Setters
+        function obj = set.total(obj, value)
+            obj.TOTAL = value;
+            obj.isDynamic   = size(value, 7) ~= 1;
+            obj.isReady     = obj.checkProperties();
+            obj             = obj.updateInterpolants();
+        end
+
+        function obj = set.geometry(obj, value)
+            if ~isempty(obj.STATE) && isfield(obj.STATE, 'phis') && ~isempty(obj.STATE.phis)
+                obj.symmetry = (360/obj.STATE.phis(end) == value.nPanel);
+            end
+
+            obj.GEOMETRY = value;
+            obj.isReady  = obj.checkProperties();
+            obj          = obj.updateInterpolants();
+        end
+
+        function obj = set.state(obj, value)
+            if isfield(value, 'phis') && ~isempty(value.phis)
+                if ~isempty(value.betas)
+                    error('Cannot set both alpha, beta and phi'); 
+                end
+                if ~isempty(obj.GEOMETRY)
+                    obj.symmetry = (360/value.phis(end) == obj.GEOMETRY.nPanel);
+                end
+                obj.type = CoeffType.AlphaTotPhi;
+            else
+                obj.type = CoeffType.AlphaBeta;
+            end
+
+            obj.STATE   = value;
+            obj.isReady = obj.checkProperties();
+            obj         = obj.updateInterpolants();
+        end
+    end
+
+    methods(Access = private)
+        function obj = loadData(obj, filePath, name)
+            % Load the coefficients from the specified file
+            % filePath: char, path to coefficient file to load
+            % name: (optional) coefficient name. Used when dynamic
+            % derivatives are also needed
+            %
+            % NOTE: Coefficient alpha and beta derivatives are not loaded
+            %       as coefficients are already interpolated linearly
+            
+
+            if ~exist(filePath, 'file'), return; end
+            warning('off'); % Disable warning in case coefficients are not present
+            dataCoeffs = load(filePath, name);
+            warning('on');
+            
+            if ~isfield(dataCoeffs, name), return; end
+            dataCoeffs = dataCoeffs.(name);
+
+            obj.finsCN = dataCoeffs.finsCN;
+            obj.GEOMETRY = dataCoeffs.geometry;
+            obj.state = dataCoeffs.state;
+            
+            % Load coefficients
+            obj.total = dataCoeffs.total;
+        end
+
+        function ready = checkProperties(obj)
+            % Check if STATIC, DYNAMIC, GEOMETRY, and STATE are non-empty
+            ready = ~isempty(obj.static) && ...
+                    ~isempty(obj.GEOMETRY) && ~isempty(obj.STATE) && ...
+                    ~(isempty(obj.dynamic) && obj.isDynamic);
+
+            if ~ready, return; end
+
+            alpha = obj.STATE.alphas;
+            mach = obj.STATE.machs;
+            altitude = obj.STATE.altitudes;
+            airbakes = obj.STATE.hprot;
+            xcg = obj.GEOMETRY.xcg;
+ 
+            if obj.type == CoeffType.AlphaBeta
+                theta = obj.STATE.betas;                             % Converting to rad
+            else
+                theta = obj.STATE.phis;                              % Converting to rad
+            end
+
+            gridVecs = {alpha, mach, theta, altitude, airbakes, xcg};
+            dims = cellfun(@(x) length(x), gridVecs);
+            dims(dims == 0) = 1;                                            % Empty case is the same as scalar case
+
+            ready = all(size(obj.total, 2:7) == dims);
+        end
+
+        function obj = updateInterpolants(obj)
+            % Update static and dynamic interpolants
+            if ~obj.isReady, return; end
+        
+            %% Retrieve flight conditions
+            alpha = obj.STATE.alphas*pi/180;                                % Converting to rad
+            mach = obj.STATE.machs;
+            altitude = obj.STATE.altitudes;
+            xcg = obj.GEOMETRY.xcg;
+
+            if isempty(obj.STATE.hprot), airbakes = 0;
+                else, airbakes = obj.STATE.hprot/obj.STATE.hprot(end);      % Normalizing on height
+            end
+            
+            if obj.type == CoeffType.AlphaBeta
+                theta = obj.STATE.betas*pi/180;                             % Converting to rad
+            else
+                theta = obj.STATE.phis*pi/180;                              % Converting to rad
+            end
+
+            gridVecs = {alpha, mach, theta, altitude, airbakes, xcg};
+
+            % Find singleton dims (last dimension is coefficients and will not be of length 1)
+            singletonDims = cellfun(@(x) isscalar(x) || isempty(x), gridVecs);
+            gridVecs(singletonDims) = deal({[0, 1]});
+
+            %% Create interpolants
+            staticCoeffs = permute(obj.static, [2, 3, 4, 5, 6, 1]);
+            staticCoeffs = repmat(staticCoeffs, singletonDims(1:end-1) + 1); % Exclude xcg for static coefficients
+
+            obj.staticInterpolant = griddedInterpolant(gridVecs(1:end-1), staticCoeffs, 'linear', 'nearest');
+
+            if ~obj.isDynamic, return; end
+
+            %% Handle dynamic derivatives
+            dynamicCoeffs = permute(obj.dynamic, [2, 3, 4, 5, 6, 7, 1]);
+
+            % Flip xcg in case of descending order
+            if ~isscalar(xcg) && xcg(2) < xcg(1)
+                gridVecs{end} = flip(xcg);
+                dynamicCoeffs = flip(dynamicCoeffs, 7);
+            end
+
+            obj.dynamicInterpolant = griddedInterpolant(gridVecs, dynamicCoeffs, 'linear', 'nearest');
+        end
+    end
+end
\ No newline at end of file
diff --git a/classes/misc/Coefficient.m b/classes/misc/Coefficient.m
index 02f721db9297a34c3ce20681ab86d28b1ca4ad4c..2573c6d0cd59a63132d03d87f50dbc751c32d30f 100644
--- a/classes/misc/Coefficient.m
+++ b/classes/misc/Coefficient.m
@@ -59,10 +59,12 @@ classdef Coefficient < Structured
         dynamic                     double  % Dynamic derivatives: [Clp, Cmad, Cmq, Cnr, Cnp]
     end
 
-
     properties(SetAccess = private)
-        isReady             (1, 1)  logical % Whether all coefficients are loaded and of the right size
-        isDynamic           (1, 1)  logical % Whether to load dynamic derivatives. Adds XCG dependece
+        type                        CoeffType   % Wheter based on alphaTot-phi or alpha-beta
+        symmetry            (1, 1)  logical     % Apply axial simmetry: valid for alphaTot-phi
+
+        isReady             (1, 1)  logical     % Whether all coefficients are loaded and of the right size
+        isDynamic           (1, 1)  logical     % Whether to load dynamic derivatives. Adds XCG dependece
     end
 
     properties(Access = private)
@@ -91,7 +93,7 @@ classdef Coefficient < Structured
             obj         = obj.updateInterpolants();
         end
 
-        function coefficients = get(obj, alpha, mach, beta, altitude, airbakes, xcg)
+        function [coefficients] = get(obj, alpha, mach, beta, altitude, airbakes, xcg)
             % GET: Retrieve aerodynamic coefficients for specified flight conditions
             %
             %   This method interpolates the aerodynamic coefficients based on
@@ -116,11 +118,32 @@ classdef Coefficient < Structured
             
             if ~obj.isReady, error('Cannot interpolate coefficients: check that all dimensions match'); end
 
+            % Interpolating coefficients
             coefficients = zeros(11, 1);
-            coefficients(1:6) = obj.staticInterpolant(alpha, mach, beta, altitude, airbakes);
-            % Transporting static coefficients to new xcg
+            if obj.type == CoeffType.AlphaBeta
+                coefficients(1:6) = obj.staticInterpolant(alpha, mach, beta, altitude, airbakes);
+            else
+                [alphaTot, phi] = getAlphaPhi(alpha, beta);
+                phi = wrapTo2Pi(phi);
+
+                finAngle = 2*pi / obj.GEOMETRY.nPanel;                           % Angle between two fin panels
+                n = floor(phi / finAngle);
+                deltaPhi = n*finAngle;  
+                psi = phi - deltaPhi;                                            % Angle wrapped to finAngle
+                
+                coefficients(1:6) = obj.staticInterpolant(alphaTot, mach, psi, altitude, airbakes); % Get coeffs in limited range
+
+                if n > 0 % If necessary, perform rotation on different frame
+                    R = [cos(deltaPhi), -sin(deltaPhi);
+                        sin(deltaPhi), cos(deltaPhi)];
+    
+                    coefficients([2, 3]) = R  * coefficients([2, 3]);
+                    coefficients([5, 6]) = R' * coefficients([5, 6]);
+                end
+            end
 
-            % C_B = C_A + d / c * [0; -CN; CY]_A <- NOTE: Non torna il meno su CN
+            % Transporting static coefficients to new xcg
+            % C_B = C_A + d / c * [0; -CN; CY]_A 
             d = xcg - obj.GEOMETRY.xcg(1);
             l = obj.GEOMETRY.diameter;
 
@@ -130,6 +153,9 @@ classdef Coefficient < Structured
             coefficients(4:6) = coefficients(4:6) + d/l * forceCoeffs;
 
             if ~obj.isDynamic, return; end
+            if obj.type == CoeffType.AlphaTotPhi
+                error('Alpha-Phi coefficients do not support dynamic derivatives yet')
+            end
             coefficients(7:11) = obj.dynamicInterpolant(alpha, mach, beta, altitude, airbakes, xcg);
         end
 
@@ -190,6 +216,10 @@ classdef Coefficient < Structured
         end
 
         function obj = set.geometry(obj, value)
+            if ~isempty(obj.STATE) && isfield(obj.STATE, 'phis') && ~isempty(obj.STATE.phis)
+                obj.symmetry = (360/obj.STATE.phis(end) == value.nPanel);
+            end
+
             obj.GEOMETRY = value;
             obj.isReady  = obj.checkProperties();
             obj          = obj.updateChad();
@@ -197,6 +227,18 @@ classdef Coefficient < Structured
         end
 
         function obj = set.state(obj, value)
+            if isfield(value, 'phis') && ~isempty(value.phis)
+                if ~isempty(value.betas)
+                    error('Cannot set both alpha, beta and phi'); 
+                end
+                if ~isempty(obj.GEOMETRY)
+                    obj.symmetry = (360/value.phis(end) == obj.GEOMETRY.nPanel);
+                end
+                obj.type = CoeffType.AlphaTotPhi;
+            else
+                obj.type = CoeffType.AlphaBeta;
+            end
+
             obj.STATE   = value;
             obj.isReady = obj.checkProperties();
             obj         = obj.updateChad();
@@ -241,12 +283,17 @@ classdef Coefficient < Structured
 
             alpha = obj.STATE.alphas;
             mach = obj.STATE.machs;
-            beta = obj.STATE.betas;
             altitude = obj.STATE.altitudes;
             airbakes = obj.STATE.hprot;
             xcg = obj.GEOMETRY.xcg;
  
-            gridVecs = {alpha, mach, beta, altitude, airbakes, xcg};
+            if obj.type == CoeffType.AlphaBeta
+                theta = obj.STATE.betas;                             % Converting to rad
+            else
+                theta = obj.STATE.phis;                              % Converting to rad
+            end
+
+            gridVecs = {alpha, mach, theta, altitude, airbakes, xcg};
             dims = cellfun(@(x) length(x), gridVecs);
             dims(dims == 0) = 1;                                            % Empty case is the same as scalar case
 
@@ -260,15 +307,20 @@ classdef Coefficient < Structured
             %% Retrieve flight conditions
             alpha = obj.STATE.alphas*pi/180;                                % Converting to rad
             mach = obj.STATE.machs;
-            beta = obj.STATE.betas*pi/180;                                  % Converting to rad
             altitude = obj.STATE.altitudes;
             xcg = obj.GEOMETRY.xcg;
 
             if isempty(obj.STATE.hprot), airbakes = 0;
                 else, airbakes = obj.STATE.hprot/obj.STATE.hprot(end);      % Normalizing on height
             end
+            
+            if obj.type == CoeffType.AlphaBeta
+                theta = obj.STATE.betas*pi/180;                             % Converting to rad
+            else
+                theta = obj.STATE.phis*pi/180;                              % Converting to rad
+            end
 
-            gridVecs = {alpha, mach, beta, altitude, airbakes, xcg};
+            gridVecs = {alpha, mach, theta, altitude, airbakes, xcg};
 
             % Find singleton dims (last dimension is coefficients and will not be of length 1)
             singletonDims = cellfun(@(x) isscalar(x) || isempty(x), gridVecs);
diff --git a/functions/miscellaneous/createDissileInput.m b/functions/miscellaneous/createDissileInput.m
index 39bc2b0b88f53030d739d657bdb7668c3fae7f6f..1c3f5eee6da918539abc97779c31585e5bc65d75 100644
--- a/functions/miscellaneous/createDissileInput.m
+++ b/functions/miscellaneous/createDissileInput.m
@@ -49,12 +49,12 @@ if not(any(vars.alpha == 1))
     error('vars.alpha does not contains 1');
 end
 
-if not(any(vars.alpha == -1))
-    error('vars.alpha does not contains -1');
+if isempty(vars.phi) && not(any(vars.alpha == -1))
+    error('vars.alpha does not contain -1 in alpha-beta configuration');
 end
 
-if not(isequal(vars.alpha, -fliplr(vars.alpha)))
-    error('vars.alpha is not symmetric');
+if isempty(vars.phi) && not(isequal(vars.alpha, -fliplr(vars.alpha)))
+    error('vars.alpha is not symmetric in alpha-beta configuration');
 end
 
 % if any(vars.abk > 1)
@@ -69,6 +69,7 @@ input.fltcon.about = 'Flight conditions quantities';
 input.fltcon.MACH = vars.mach;
 input.fltcon.ALPHA = vars.alpha;
 input.fltcon.BETA = vars.beta;
+input.fltcon.PHI = vars.phi;
 input.fltcon.ALT = vars.alt;
 
 %% REFQ
diff --git a/functions/odeFunctions/ballistic.m b/functions/odeFunctions/ballistic.m
index d08b3069d7274836e0f686c19963919facaec233..fed5600a547105015d3835c998120ef8accf7f44 100644
--- a/functions/odeFunctions/ballistic.m
+++ b/functions/odeFunctions/ballistic.m
@@ -134,7 +134,7 @@ end
 if not(abs(ur) < 1e-9 || velsNorm < 1e-9)
     alpha = atan(wr/ur);
     beta = atan(vr/ur);                                                     % beta = asin(vr/V_norm) is the classical notation, Datcom uses this one though.
-    % alpha_tot = atan(sqrt(wr^2 + vr^2)/ur);   % datcom 97' definition
+    alphaTot = atan(sqrt(wr^2 + vr^2)/ur);   % datcom 97' definition
 else
     alpha = 0;
     beta = 0;
@@ -144,14 +144,17 @@ alphaOut = alpha;
 betaOut = beta;
 
 %% INTERPOLATE AERODYNAMIC COEFFICIENTS:
-
-if abs(alpha) > rocket.coefficients.state.alphas(end)*pi/180 || ...
-        abs(beta) > rocket.coefficients.state.betas(end)*pi/180
-    coeffsValues = rocket.coefficientsHighAOA.get(alpha, mach, beta, absoluteAltitude, extension, xcg);
+if rocket.coefficients.type == CoeffType.AlphaBeta
+    if abs(alpha) > rocket.coefficients.state.alphas(end)*pi/180 || abs(beta) > rocket.coefficients.state.betas(end)*pi/180
+        coeffsValues = rocket.coefficientsHighAOA.get(alpha, mach, beta, absoluteAltitude, extension, xcg);
+    end
+elseif abs(alphaTot) > rocket.coefficients.state.alphas(end)*pi/180
+        coeffsValues = rocket.coefficientsHighAOA.get(alpha, mach, beta, absoluteAltitude, extension, xcg);
 else
     coeffsValues = rocket.coefficients.get(alpha, mach, beta, absoluteAltitude, extension, xcg);
 end
 
+
 % Retrieve Coefficients
 CA = coeffsValues(1); CY = coeffsValues(2); CN = coeffsValues(3);
 Cl = coeffsValues(4); Cm = coeffsValues(5); Cn = coeffsValues(6);
diff --git a/missions/2021_Lynx_Portugal_October/config/rocketConfig.m b/missions/2021_Lynx_Portugal_October/config/rocketConfig.m
index 160c5b09b850ef81406e13cf89dfcda9005c2e9c..cb264015ae392fd9ed08b661398f6b5efb8f1d2e 100644
--- a/missions/2021_Lynx_Portugal_October/config/rocketConfig.m
+++ b/missions/2021_Lynx_Portugal_October/config/rocketConfig.m
@@ -20,6 +20,7 @@ rocket.lengthCenterNoMot = 1.7640;                                     % [m]
 %   When false, coefficients are saved with current motor's name
 %   When true,  coefficients are saved as 'generic'
 rocket.dynamicDerivatives = false;                                    % [-]      True if dynamic derivatives will be loaded
+rocket.highAOACoefficients = false;                                     % [-]      True if separate matrix for high AOA will be loaded
 
 %% PLD - Includes Parafoil + Nose
 parafoil = Parafoil();
diff --git a/missions/2021_Lynx_Roccaraso_September/config/rocketConfig.m b/missions/2021_Lynx_Roccaraso_September/config/rocketConfig.m
index d0dc3df55356e3aebc59d51214a30794f8668dc1..0851f94e15aa603fff564f4dd2d8ed15c0e56bf8 100644
--- a/missions/2021_Lynx_Roccaraso_September/config/rocketConfig.m
+++ b/missions/2021_Lynx_Roccaraso_September/config/rocketConfig.m
@@ -20,6 +20,7 @@ rocket.lengthCenterNoMot = 1.7840;                                     % [m]
 %   When false, coefficients are saved with current motor's name
 %   When true,  coefficients are saved as 'generic'
 rocket.dynamicDerivatives = false;                                    % [-]      True if dynamic derivatives will be loaded
+rocket.highAOACoefficients = false;                                     % [-]      True if separate matrix for high AOA will be loaded
 
 %% PLD - Includes Parafoil + Nose
 parafoil = Parafoil();
diff --git a/missions/2022_Pyxis_Portugal_October/config/rocketConfig.m b/missions/2022_Pyxis_Portugal_October/config/rocketConfig.m
index 2f949bfb46c9b09dd02fea6a1c87253e9c6c5b20..17d63ca520136746a7f874cb8e2d9c34c726670b 100644
--- a/missions/2022_Pyxis_Portugal_October/config/rocketConfig.m
+++ b/missions/2022_Pyxis_Portugal_October/config/rocketConfig.m
@@ -20,6 +20,7 @@ rocket.lengthCenterNoMot = 1.4470;                                     % [m]
 %   When false, coefficients are saved with current motor's name
 %   When true,  coefficients are saved as 'generic'
 rocket.dynamicDerivatives = false;                                    % [-]      True if dynamic derivatives will be loaded
+rocket.highAOACoefficients = false;                                     % [-]      True if separate matrix for high AOA will be loaded
 
 %% PLD - Includes Parafoil + Nose
 parafoil = Parafoil();
diff --git a/missions/2022_Pyxis_Roccaraso_September/config/rocketConfig.m b/missions/2022_Pyxis_Roccaraso_September/config/rocketConfig.m
index ee3a3c71bd73ab263db92f3e929db488a6a5735c..5b7406a3e5a8362c54c1caddde3a35c9fa1798cc 100644
--- a/missions/2022_Pyxis_Roccaraso_September/config/rocketConfig.m
+++ b/missions/2022_Pyxis_Roccaraso_September/config/rocketConfig.m
@@ -22,6 +22,7 @@ rocket.lengthCenterNoMot = 1.61;                                     % [m]
 %   When false, coefficients are saved with current motor's name
 %   When true,  coefficients are saved as 'generic'
 rocket.dynamicDerivatives = false;                                    % [-]      True if dynamic derivatives will be loaded
+rocket.highAOACoefficients = false;                                     % [-]      True if separate matrix for high AOA will be loaded
 
 %% PLD - Includes Parafoil + Nose
 parafoil = Parafoil();
diff --git a/missions/2023_Gemini_Portugal_October/config/rocketConfig.m b/missions/2023_Gemini_Portugal_October/config/rocketConfig.m
index 401c415be5e782af00b4318e62db4d13ea6ef5ca..f9d6d4c6f95e6d93445e6289d5115b9bba35837a 100644
--- a/missions/2023_Gemini_Portugal_October/config/rocketConfig.m
+++ b/missions/2023_Gemini_Portugal_October/config/rocketConfig.m
@@ -19,7 +19,7 @@ rocket.lengthCenterNoMot = 1.517;                                     % [m]
 %   When false, coefficients are saved with current motor's name
 %   When true,  coefficients are saved as 'generic'
 rocket.dynamicDerivatives = false;                                    % [-]      True if dynamic derivatives will be loaded
-
+rocket.highAOACoefficients = false;                                     % [-]      True if separate matrix for high AOA will be loaded
 
 %% PLD - Includes Parafoil + Nose
 parafoil = Parafoil();
diff --git a/missions/2023_Gemini_Roccaraso_September/config/rocketConfig.m b/missions/2023_Gemini_Roccaraso_September/config/rocketConfig.m
index fc2715bc0dd7d61dd740e20b3138c8552440dcf1..cadbc7ef2b84f7f504eae85f73766ece8a83ec9c 100644
--- a/missions/2023_Gemini_Roccaraso_September/config/rocketConfig.m
+++ b/missions/2023_Gemini_Roccaraso_September/config/rocketConfig.m
@@ -19,6 +19,7 @@ rocket.lengthCenterNoMot = 1.517;                                     % [m]
 %   When false, coefficients are saved with current motor's name
 %   When true,  coefficients are saved as 'generic'
 rocket.dynamicDerivatives = false;                                    % [-]      True if dynamic derivatives will be loaded
+rocket.highAOACoefficients = false;                                     % [-]      True if separate matrix for high AOA will be loaded
 
 %% PLD - Includes Parafoil + Nose
 parafoil = Parafoil();
diff --git a/missions/2024_Lyra_Portugal_October/config/rocketConfig.m b/missions/2024_Lyra_Portugal_October/config/rocketConfig.m
index f87640c176ca52f0ab3274f1b32602462a1489ff..1d90d15565d46515b3ca9ef5bcc7677b369fd810 100644
--- a/missions/2024_Lyra_Portugal_October/config/rocketConfig.m
+++ b/missions/2024_Lyra_Portugal_October/config/rocketConfig.m
@@ -13,6 +13,7 @@ rocket.lengthCenterNoMot = [];                                          % [m]
 %   When false, coefficients are saved with current motor's name
 %   When true,  coefficients are saved as 'generic'
 rocket.dynamicDerivatives = false;                                    % [-]      True if dynamic derivatives will be loaded
+rocket.highAOACoefficients = false;                                     % [-]      True if separate matrix for high AOA will be loaded
 
 %% PLD - Includes Parafoil + Nose
 parafoil = Parafoil();
diff --git a/missions/2024_Lyra_Roccaraso_September/config/rocketConfig.m b/missions/2024_Lyra_Roccaraso_September/config/rocketConfig.m
index 2939db0ce6550575817cdd355d8b2f6ca2e68943..133be74e22c422360b80e924c6fb291a147c61db 100644
--- a/missions/2024_Lyra_Roccaraso_September/config/rocketConfig.m
+++ b/missions/2024_Lyra_Roccaraso_September/config/rocketConfig.m
@@ -19,6 +19,7 @@ rocket.lengthCenterNoMot = [];                                          % [m]
 %   When false, coefficients are saved with current motor's name
 %   When true,  coefficients are saved as 'generic'
 rocket.dynamicDerivatives = false;                                    % [-]      True if dynamic derivatives will be loaded
+rocket.highAOACoefficients = false;                                     % [-]      True if separate matrix for high AOA will be loaded
 
 %% PLD - Includes Parafoil + Nose
 parafoil = Parafoil();
diff --git a/missions/2025_Orion_Portugal_October/config/rocketConfig.m b/missions/2025_Orion_Portugal_October/config/rocketConfig.m
index e73d2a22bc555fc4cefd08600e01497798ccb620..6e2c6c63d6d7c2f6d22e1d2a05f6e40766512558 100644
--- a/missions/2025_Orion_Portugal_October/config/rocketConfig.m
+++ b/missions/2025_Orion_Portugal_October/config/rocketConfig.m
@@ -12,7 +12,8 @@ rocket.diameter = 0.15;                                                   % [m]
 % If dynamic derivatives are loaded, coefficient will depend on rocket xcg
 %   When false, coefficients are saved with current motor's name
 %   When true,  coefficients are saved as 'generic'
-rocket.dynamicDerivatives = false;                                    % [-]      True if dynamic derivatives will be loaded
+rocket.dynamicDerivatives = false;                                      % [-]      True if dynamic derivatives will be loaded
+rocket.highAOACoefficients = false;                                     % [-]      True if separate matrix for high AOA will be loaded
 
 %% PRF - Includes Parafoil + Nose
 parafoil = Parafoil();
diff --git a/missions/2025_Orion_Portugal_October/data/aeroCoefficients.mat b/missions/2025_Orion_Portugal_October/data/aeroCoefficients.mat
index 1df11798e89a63e54468544fd12b9fb0f90002e5..e7019e85c88b0ac771bc85569864afe0be9298dd 100644
--- a/missions/2025_Orion_Portugal_October/data/aeroCoefficients.mat
+++ b/missions/2025_Orion_Portugal_October/data/aeroCoefficients.mat
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:fb0ccbfea51866683e3cdaf1b7aa21c42ca3d387575aaa40df3694557c009fd4
-size 14789031
+oid sha256:f544fdb75b3dcc7180915ef4b0263969a17444b7609cb7e8947fc482d5cffa14
+size 300751768
diff --git a/missions/2025_Orion_Portugal_October/data/aeroCoefficientsHighAOA.mat b/missions/2025_Orion_Portugal_October/data/aeroCoefficientsHighAOA.mat
index 965f46e73771a4b87becdd6b43fd7e2c83959090..00d9ce5686896e88b5ebfa31fd23bb911eaa6da2 100644
--- a/missions/2025_Orion_Portugal_October/data/aeroCoefficientsHighAOA.mat
+++ b/missions/2025_Orion_Portugal_October/data/aeroCoefficientsHighAOA.mat
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:f6f664625606e3fdfcbf0de0d8a45b8f49a43be69d8a879104f1315d4d9f5e53
-size 5776386
+oid sha256:1659d1dc6b19ea13c3cdd75b8f9627d42123c4cb6dff535d93618247a17ab9d1
+size 46224320
diff --git a/missions/2025_Orion_Roccaraso_September/config/rocketConfig.m b/missions/2025_Orion_Roccaraso_September/config/rocketConfig.m
index c4e1ffcee38a736e09fe06543010f7699fe3a553..e320c1583c941f7fe5b25b5c54d0a4ee114c5e15 100644
--- a/missions/2025_Orion_Roccaraso_September/config/rocketConfig.m
+++ b/missions/2025_Orion_Roccaraso_September/config/rocketConfig.m
@@ -19,6 +19,7 @@ rocket.lengthCenterNoMot = [];                                          % [m]
 %   When false, coefficients are saved with current motor's name
 %   When true,  coefficients are saved as 'generic'
 rocket.dynamicDerivatives = false;                                    % [-]      True if dynamic derivatives will be loaded
+rocket.highAOACoefficients = false;                                     % [-]      True if separate matrix for high AOA will be loaded
 
 %% PLD - Includes Parafoil + Nose
 parafoil = Parafoil();