CombustionToolbox · AlbertoCuadra · Oct 19, 2025 · Oct 19, 2025
diff --git a/+combustiontoolbox/+core/@CaloricGasModel/CaloricGasModel.m b/+combustiontoolbox/+core/@CaloricGasModel/CaloricGasModel.m
@@ -0,0 +1,126 @@
+classdef CaloricGasModel < uint8
+    % The :mat:class:`CaloricGasModel` class is an enumeration that defines
+    % the available caloric gas models for thermodynamic property calculations.
+    %
+    % The available caloric gas models are:
+    %     * `perfect`: Calorically perfect gas model
+    %     * `thermallyPerfect`: Thermally perfect gas model
+    %     * `imperfect`: Calorically imperfect gas model
+    %
+    % Example:
+    %     model = CaloricGasModel.perfect;
+
+    enumeration
+        perfect          (0)
+        thermallyPerfect (1)
+        imperfect        (2)
+    end
+
+    methods
+
+        function obj = setPerfect(obj)
+            % Set the caloric gas model to perfect.
+            %
+            % Example:
+            %     model = model.setPerfect()
+
+            obj = combustiontoolbox.core.CaloricGasModel.perfect;
+        end
+
+        function obj = setThermallyPerfect(obj)
+            % Set the caloric gas model to thermally perfect.
+            %
+            % Example:
+            %     model = model.setThermallyPerfect()
+
+            obj = combustiontoolbox.core.CaloricGasModel.thermallyPerfect;
+        end
+
+        function obj = setImperfect(obj)
+            % Set the caloric gas model to imperfect.
+            %
+            % Example:
+            %     model = model.setImperfect()
+
+            obj = combustiontoolbox.core.CaloricGasModel.imperfect;
+        end
+
+        function value = isPerfect(obj)
+            % Check if the caloric gas model is perfect.
+            %
+            % Returns:
+            %     value (logical): True if the model is perfect, false otherwise.
+            %
+            % Example:
+            %     isPerfect = model.isPerfect()
+
+            value = (obj == combustiontoolbox.core.CaloricGasModel.perfect);
+        end
+
+        function value = isThermallyPerfect(obj)
+            % Check if the caloric gas model is thermally perfect.
+            %
+            % Returns:
+            %     value (logical): True if the model is thermally perfect, false otherwise.
+            %
+            % Example:
+            %     isThermallyPerfect = model.isThermallyPerfect()
+
+            value = (obj == combustiontoolbox.core.CaloricGasModel.thermallyPerfect);
+        end
+
+        function value = isImperfect(obj)
+            % Check if the caloric gas model is imperfect.
+            %
+            % Returns:
+            %     value (logical): True if the model is imperfect, false otherwise.
+            %
+            % Example:
+            %     isImperfect = model.isImperfect()
+
+            value = (obj == combustiontoolbox.core.CaloricGasModel.imperfect);
+        end
+
+    end
+
+    methods (Static)
+
+        function model = fromFlag(FLAG_TCHEM_FROZEN, FLAG_FROZEN)
+            % Convert legacy solver flags to CaloricGasModel enumeration.
+            %
+            % Args:
+            %     FLAG_TCHEM_FROZEN (logical): True if calorically perfect gas
+            %     FLAG_FROZEN (logical): True if thermally perfect gas
+            %
+            % Returns:
+            %     model (CaloricGasModel): Corresponding caloric gas model
+            %
+            % Legacy compatibility:
+            %     * FLAG_TCHEM_FROZEN = true,  FLAG_FROZEN = false  --- perfect
+            %     * FLAG_TCHEM_FROZEN = false, FLAG_FROZEN = true   --- thermallyPerfect
+            %     * FLAG_TCHEM_FROZEN = false, FLAG_FROZEN = false  --- imperfect
+            %
+            % Example:
+            %     model = CaloricGasModel.fromFlag(true, false)
+
+            % Check inputs
+            if ~islogical(FLAG_TCHEM_FROZEN) || ~islogical(FLAG_FROZEN)
+                error('Input flags must be logical values.');
+            end
+
+            % Map legacy flags to caloric model
+            if FLAG_TCHEM_FROZEN && FLAG_FROZEN
+                error('Incompatible flags: Both FLAG_TCHEM_FROZEN and FLAG_FROZEN cannot be true simultaneously.');
+            elseif FLAG_TCHEM_FROZEN
+                model = combustiontoolbox.core.CaloricGasModel.perfect;
+            elseif FLAG_FROZEN
+                model = combustiontoolbox.core.CaloricGasModel.thermallyPerfect;
+            else
+                model = combustiontoolbox.core.CaloricGasModel.imperfect;
+            end
+
+        end
+
+    end
+
+end
diff --git a/+combustiontoolbox/+equilibrium/@EquilibriumSolver/EquilibriumSolver.m b/+combustiontoolbox/+equilibrium/@EquilibriumSolver/EquilibriumSolver.m
@@ -24,6 +24,7 @@
 
     properties
         problemType                % Problem type [TP, TV, HP, EV, SP, SV]
+        caloricGasModel            % Caloric gas model (perfect, thermallyPerfect, imperfect)
         tolGibbs = 1e-6            % Tolerance of the Gibbs/Helmholtz minimization method
         tolE = 1e-6                % Tolerance of the mass balance
         tolMoles = 1e-14           % Tolerance of the composition of the mixture                       
@@ -63,10 +64,12 @@
             % Constructor
             defaultProblemType = 'TP';
             defaultPlotConfig = combustiontoolbox.utils.display.PlotConfig();
+            defaultCaloricGasModel = combustiontoolbox.core.CaloricGasModel.imperfect;
 
             % Parse input arguments
             p = inputParser;
             addOptional(p, 'problemType', defaultProblemType, @(x) ischar(x) && any(strcmpi(x, {'TP', 'TV', 'HP', 'EV', 'SP', 'SV'})));
+            addParameter(p, 'caloricGasModel', defaultCaloricGasModel, @(x) isa(x, 'combustiontoolbox.core.CaloricGasModel'));
             addParameter(p, 'tolGibbs', obj.tolGibbs, @(x) isnumeric(x) && x > 0);
             addParameter(p, 'tolE', obj.tolE, @(x) isnumeric(x) && x > 0);
             addParameter(p, 'tolMoles', obj.tolMoles, @(x) isnumeric(x) && x > 0);
@@ -97,6 +100,7 @@
 
             % Set properties
             obj.problemType = upper(p.Results.problemType);
+            obj.caloricGasModel = p.Results.caloricGasModel;
             obj.tolGibbs = p.Results.tolGibbs;
             obj.tolE = p.Results.tolE;
             obj.tolMoles = p.Results.tolMoles;
@@ -123,6 +127,15 @@
             obj.FLAG_REPORT = p.Results.FLAG_REPORT;
             obj.FLAG_CACHE = p.Results.FLAG_CACHE;
             obj.plotConfig = p.Results.plotConfig;
+
+            % Display warning if deprecated flags are used
+            if ~ismember('FLAG_TCHEM_FROZEN', p.UsingDefaults) || ~ismember('FLAG_FROZEN', p.UsingDefaults)
+                warning(['The flags ''FLAG_TCHEM_FROZEN'' and ''FLAG_FROZEN'' are deprecated. ', ...
+                         'Please use the ''caloricGasModel'' parameter with values from the CaloricGasModel enumeration instead.']);
+
+                obj.caloricGasModel = obj.caloricGasModel.fromFlag(obj.FLAG_TCHEM_FROZEN, obj.FLAG_FROZEN);
+            end
+
         end
 
         function obj = set(obj, property, value, varargin)

diff --git a/+combustiontoolbox/+equilibrium/@EquilibriumSolver/equilibrateT.m b/+combustiontoolbox/+equilibrium/@EquilibriumSolver/equilibrateT.m
@@ -20,20 +20,21 @@
 
     % Import packages
     import combustiontoolbox.utils.findIndex
-    
-    % Check if calculations are for a thermochemical frozen gas (calorically perfect gas)
-    if obj.FLAG_TCHEM_FROZEN
+
+    % Check if calculations are for a calorically perfect gas
+    if obj.caloricGasModel.isPerfect()
         mix2 = equilibrateTPerfect(mix1, mix2, T);
         return
     end
 
-    % Check if calculations are for a calorically imperfect gas with frozen chemistry (thermally perfect gas)
-    if obj.FLAG_FROZEN
+    % Check if calculations are for a thermally perfect gas
+    if obj.caloricGasModel.isThermallyPerfect()
         mix2 = equilibrateTFrozen(mix2, T);
         return
     end
 
     % Definitions
+    molesGuess = [];
     N_mix0 = moles(mix1); % Get moles of inert species
     system = mix2.chemicalSystem;
     systemProducts = mix2.chemicalSystemProducts;

diff --git a/+combustiontoolbox/+rocket/@RocketSolver/RocketSolver.m b/+combustiontoolbox/+rocket/@RocketSolver/RocketSolver.m
@@ -35,6 +35,7 @@
         function obj = RocketSolver(varargin)
             % Constructor
             defaultProblemType = 'ROCKET_IAC';
+            defaultCaloricGasModel = combustiontoolbox.core.CaloricGasModel.imperfect;
             defaultEquilibriumSolver = combustiontoolbox.equilibrium.EquilibriumSolver();
             defaultPlotConfig = combustiontoolbox.utils.display.PlotConfig();
             defaultFLAG_TCHEM_FROZEN = false;
@@ -43,6 +44,7 @@
             % Parse input arguments
             p = inputParser;
             addOptional(p, 'problemType', defaultProblemType, @(x) ischar(x) && any(strcmpi(x, {'ROCKET_IAC', 'ROCKET_FAC'})));
+            addParameter(p, 'caloricGasModel', defaultCaloricGasModel, @(x) isa(x, 'combustiontoolbox.core.CaloricGasModel'));
             addParameter(p, 'equilibriumSolver', defaultEquilibriumSolver);
             addParameter(p, 'FLAG_TCHEM_FROZEN', defaultFLAG_TCHEM_FROZEN, @(x) islogical(x))
             addParameter(p, 'FLAG_FROZEN', defaultFLAG_FROZEN, @(x) islogical(x));
@@ -64,6 +66,7 @@
             obj.FLAG_CACHE = p.Results.FLAG_CACHE;
 
             if sum(contains(p.UsingDefaults, 'equilibriumSolver'))
+                obj.equilibriumSolver.caloricGasModel = p.Results.caloricGasModel;
                 obj.equilibriumSolver.FLAG_TCHEM_FROZEN = p.Results.FLAG_TCHEM_FROZEN;
                 obj.equilibriumSolver.FLAG_FROZEN = p.Results.FLAG_FROZEN;
                 obj.equilibriumSolver.tolMoles = p.Results.tolMoles;
@@ -75,6 +78,15 @@
             obj.equilibriumSolver.FLAG_CACHE = false;
             obj.plotConfig.plotProperties = [obj.plotConfig.plotProperties, {'u', 'I_sp', 'I_vac'}];
             obj.plotConfig.plotPropertiesBasis = [obj.plotConfig.plotPropertiesBasis, {[], [], []}];
+
+            % Display warning if deprecated flags are used
+            if ~ismember('FLAG_TCHEM_FROZEN', p.UsingDefaults) || ~ismember('FLAG_FROZEN', p.UsingDefaults)
+                warning(['The flags ''FLAG_TCHEM_FROZEN'' and ''FLAG_FROZEN'' are deprecated. ', ...
+                         'Please use the ''caloricGasModel'' parameter with values from the CaloricGasModel enumeration instead.']);
+
+                obj.equilibriumSolver.caloricGasModel = obj.equilibriumSolver.caloricGasModel.fromFlag(obj.equilibriumSolver.FLAG_TCHEM_FROZEN, obj.equilibriumSolver.FLAG_FROZEN);
+            end
+
         end
 
         function obj = set(obj, property, value, varargin)

diff --git a/+combustiontoolbox/+rocket/@RocketSolver/rocketChamberIAC.m b/+combustiontoolbox/+rocket/@RocketSolver/rocketChamberIAC.m
@@ -17,11 +17,11 @@
     %     mix = rocketChamberIAC(obj, mix, mix_guess)
 
     % Temporal value
-    TEMP_FLAG_FROZEN = obj.equilibriumSolver.FLAG_FROZEN;
+    tempCaloricGasModel = obj.equilibriumSolver.caloricGasModel;
 
     % Definitions
     obj.equilibriumSolver.problemType = 'HP';
-    obj.equilibriumSolver.FLAG_FROZEN = false;
+    obj.equilibriumSolver.caloricGasModel = obj.equilibriumSolver.caloricGasModel.setImperfect;
 
     % Compute chemical equilibria at the exit of the chamber (HP)
     if isempty(mix_guess)
@@ -35,5 +35,5 @@
 
     % Restore problemType and FLAG_FROZEN values;
     mix.problemType = 'ROCKET_IAC';
-    obj.equilibriumSolver.FLAG_FROZEN = TEMP_FLAG_FROZEN;
+    obj.equilibriumSolver.caloricGasModel = tempCaloricGasModel;
 end
diff --git a/+combustiontoolbox/+shockdetonation/@DetonationSolver/DetonationSolver.m b/+combustiontoolbox/+shockdetonation/@DetonationSolver/DetonationSolver.m
@@ -61,6 +61,7 @@
         function obj = DetonationSolver(varargin)
             % Constructor
             defaultProblemType = 'DET';
+            defaultCaloricGasModel = combustiontoolbox.core.CaloricGasModel.imperfect;
             defaultPlotConfig = combustiontoolbox.utils.display.PlotConfig();
             defaultEquilibriumSolver = combustiontoolbox.equilibrium.EquilibriumSolver('plotConfig', defaultPlotConfig);
             defaultShockSolver = combustiontoolbox.shockdetonation.ShockSolver('plotConfig', defaultPlotConfig);
@@ -70,6 +71,7 @@
             % Parse input arguments
             p = inputParser;
             addOptional(p, 'problemType', defaultProblemType, @(x) ischar(x) && any(strcmpi(x, {'DET', 'DET_R', 'DET_OVERDRIVEN', 'DET_OVERDRIVEN_R', 'DET_UNDERDRIVEN', 'DET_UNDERDRIVEN_R', 'DET_OBLIQUE', 'DET_OBLIQUE_R', 'DET_POLAR', 'DET_POLAR_R'})));
+            addParameter(p, 'caloricGasModel', defaultCaloricGasModel, @(x) isa(x, 'combustiontoolbox.core.CaloricGasModel'));
             addParameter(p, 'equilibriumSolver', defaultEquilibriumSolver);
             addParameter(p, 'shockSolver', defaultShockSolver);
             addParameter(p, 'tol0', obj.tol0, @(x) isnumeric(x) && x > 0);
@@ -109,6 +111,7 @@
             obj.plotConfig = p.Results.plotConfig;
 
             if sum(contains(p.UsingDefaults, 'equilibriumSolver'))
+                obj.equilibriumSolver.caloricGasModel = p.Results.caloricGasModel;
                 obj.equilibriumSolver.FLAG_TCHEM_FROZEN = p.Results.FLAG_TCHEM_FROZEN;
                 obj.equilibriumSolver.FLAG_FROZEN = p.Results.FLAG_FROZEN;
                 obj.equilibriumSolver.tolMoles = p.Results.tolMoles;
@@ -123,6 +126,15 @@
             obj.shockSolver.FLAG_CACHE = false;
             obj.plotConfig.plotProperties{end + 1} = 'uShock';
             obj.plotConfig.plotPropertiesBasis{end + 1} = [];
+
+            % Display warning if deprecated flags are used
+            if ~ismember('FLAG_TCHEM_FROZEN', p.UsingDefaults) || ~ismember('FLAG_FROZEN', p.UsingDefaults)
+                warning(['The flags ''FLAG_TCHEM_FROZEN'' and ''FLAG_FROZEN'' are deprecated. ', ...
+                         'Please use the ''caloricGasModel'' parameter with values from the CaloricGasModel enumeration instead.']);
+
+                obj.equilibriumSolver.caloricGasModel = obj.equilibriumSolver.caloricGasModel.fromFlag(obj.equilibriumSolver.FLAG_TCHEM_FROZEN, obj.equilibriumSolver.FLAG_FROZEN);
+            end
+
         end
 
         function obj = set(obj, property, value, varargin)

diff --git a/+combustiontoolbox/+shockdetonation/@JumpConditionsSolver/JumpConditionsSolver.m b/+combustiontoolbox/+shockdetonation/@JumpConditionsSolver/JumpConditionsSolver.m
@@ -95,6 +95,7 @@
             %     obj = JumpConditionsSolver('param1', value1, 'param2', value2, ...);
 
             % Default values
+            defaultCaloricGasModel = combustiontoolbox.core.CaloricGasModel.imperfect;
             defaultEquilibriumSolver = combustiontoolbox.equilibrium.EquilibriumSolver();
             defaultShockSolver = combustiontoolbox.shockdetonation.ShockSolver('equilibriumSolver', defaultEquilibriumSolver);
             defaultFLAG_TCHEM_FROZEN = false;
@@ -106,6 +107,7 @@
 
             % Parse input arguments
             p = inputParser;
+            addParameter(p, 'caloricGasModel', defaultCaloricGasModel, @(x) isa(x, 'combustiontoolbox.core.CaloricGasModel'));
             addParameter(p, 'equilibriumSolver', defaultEquilibriumSolver, @(x) isa(x, 'combustiontoolbox.equilibrium.EquilibriumSolver'));
             addParameter(p, 'shockSolver', defaultShockSolver, @(x) isa(x, 'combustiontoolbox.shockdetonation.ShockSolver'));
             addParameter(p, 'tolGammas1', obj.tolGammas1, @(x) isnumeric(x) && isscalar(x) && x > 0);
@@ -127,6 +129,7 @@
             parse(p, varargin{:});
 
             % Set parameters
+            obj.caloricGasModel = p.Results.caloricGasModel;
             obj.equilibriumSolver = p.Results.equilibriumSolver;
             obj.shockSolver = p.Results.shockSolver;
             obj.tolGammas1 = p.Results.tolGammas1;
@@ -158,6 +161,15 @@
             obj.equilibriumSolver.FLAG_TIME = false;
             obj.equilibriumSolver.FLAG_CACHE = false;
             obj.shockSolver.FLAG_TIME = false;
+
+            % Display warning if deprecated flags are used
+            if ~ismember('FLAG_TCHEM_FROZEN', p.UsingDefaults) || ~ismember('FLAG_FROZEN', p.UsingDefaults)
+                warning(['The flags ''FLAG_TCHEM_FROZEN'' and ''FLAG_FROZEN'' are deprecated. ', ...
+                         'Please use the ''caloricGasModel'' parameter with values from the CaloricGasModel enumeration instead.']);
+
+                obj.equilibriumSolver.caloricGasModel = obj.equilibriumSolver.caloricGasModel.fromFlag(obj.equilibriumSolver.FLAG_TCHEM_FROZEN, obj.equilibriumSolver.FLAG_FROZEN);
+            end
+
         end
 
         function value = get.MIN_MACH(obj)
@@ -266,7 +278,7 @@
             jumpConditions = getJumpData(obj);
 
             % Solve Gammas (calorically perfect gas)
-            if obj.equilibriumSolver.FLAG_TCHEM_FROZEN
+            if obj.caloricGasModel.isPerfect()
                 [jumpConditions.Gammas1, jumpConditions.Gammas2, jumpConditions.Gammas3] = obj.getGammasPerfect(jumpConditions.gamma1, jumpConditions.M1);
 
                 % Interpolate results into a smaller grid

diff --git a/+combustiontoolbox/+shockdetonation/@ShockSolver/ShockSolver.m b/+combustiontoolbox/+shockdetonation/@ShockSolver/ShockSolver.m
@@ -57,6 +57,7 @@
         function obj = ShockSolver(varargin)
             % Constructor
             defaultProblemType = 'SHOCK_I';
+            defaultCaloricGasModel = combustiontoolbox.core.CaloricGasModel.imperfect;
             defaultEquilibriumSolver = combustiontoolbox.equilibrium.EquilibriumSolver();
             defaultPlotConfig = combustiontoolbox.utils.display.PlotConfig();
             defaultFLAG_TCHEM_FROZEN = false;
@@ -65,6 +66,7 @@
             % Parse input arguments
             p = inputParser;
             addOptional(p, 'problemType', defaultProblemType, @(x) ischar(x) && any(strcmpi(x, {'SHOCK_I', 'SHOCK_R', 'SHOCK_OBLIQUE', 'SHOCK_OBLIQUE_R', 'SHOCK_POLAR', 'SHOCK_POLAR_R', 'SHOCK_POLAR_LIMITRR', 'SHOCK_PRANDTL_MEYER'})));
+            addParameter(p, 'caloricGasModel', defaultCaloricGasModel, @(x) isa(x, 'combustiontoolbox.core.CaloricGasModel'));
             addParameter(p, 'equilibriumSolver', defaultEquilibriumSolver);
             addParameter(p, 'tol0', obj.tol0, @(x) isnumeric(x) && x > 0);
             addParameter(p, 'itMax', obj.itMax, @(x) isnumeric(x) && x > 0);
@@ -103,6 +105,7 @@
             obj.plotConfig = p.Results.plotConfig;
 
             if sum(contains(p.UsingDefaults, 'equilibriumSolver'))
+                obj.equilibriumSolver.caloricGasModel = p.Results.caloricGasModel;
                 obj.equilibriumSolver.FLAG_TCHEM_FROZEN = p.Results.FLAG_TCHEM_FROZEN;
                 obj.equilibriumSolver.FLAG_FROZEN = p.Results.FLAG_FROZEN;
             end
@@ -111,6 +114,15 @@
             obj.equilibriumSolver.FLAG_RESULTS = false;
             obj.equilibriumSolver.FLAG_TIME = false;
             obj.equilibriumSolver.FLAG_CACHE = false;
+
+            % Display warning if deprecated flags are used
+            if ~ismember('FLAG_TCHEM_FROZEN', p.UsingDefaults) || ~ismember('FLAG_FROZEN', p.UsingDefaults)
+                warning(['The flags ''FLAG_TCHEM_FROZEN'' and ''FLAG_FROZEN'' are deprecated. ', ...
+                         'Please use the ''caloricGasModel'' parameter with values from the CaloricGasModel enumeration instead.']);
+
+                obj.equilibriumSolver.caloricGasModel = obj.equilibriumSolver.caloricGasModel.fromFlag(obj.equilibriumSolver.FLAG_TCHEM_FROZEN, obj.equilibriumSolver.FLAG_FROZEN);
+            end
+
         end
 
         function obj = set(obj, property, value, varargin)