diff --git a/Settings/functions/Display/plot_shock_polar.m b/Settings/functions/Display/plot_shock_polar.m
new file mode 100644
index 00000000..29a80f5b
--- /dev/null
+++ b/Settings/functions/Display/plot_shock_polar.m
@@ -0,0 +1,70 @@
+function [ax1, ax2, ax3] = plot_shock_polar(self, mix1, mix2)
+    % Routine to obtain shock polar plots
+    %   * Plot (pressure, deflection) 
+    %   * Plot (wave angle, deflection)
+    %   * Plot velocity components
+
+    % Abbreviations
+    config = self.Misc.config;
+    % Definitions
+    N = length(mix1);
+    % Plots
+    for i = N:-1:1
+        % Plot (pressure, deflection) - shock polar
+	    ax1 = plot_shock_polar_pressure(mix1{i}, mix2{i}, config);
+        % Plot (wave angle, deflection) - shock polar
+        ax2 = plot_shock_polar_wave(mix1{i}, mix2{i}, config);
+        % Plot velocity components
+        ax3 = plot_shock_polar_velocities(mix1{i}, mix2{i}, config);
+    end
+end
+
+% SUB-PASS FUNCTIONS
+function ax = set_fixed_figure(fixed_number, config)
+    % Generate a figure with a defined object identifier (fixed_number)
+    figure(fixed_number); ax = gca;
+    set(ax, 'LineWidth', config.linewidth, 'FontSize', config.fontsize-2, 'BoxStyle', 'full');
+    grid(ax, 'off'); box(ax, 'off'); hold(ax, 'on'); ax.Layer = 'Top'; axis(ax, 'tight');
+end
+
+function ax = plot_shock_polar_pressure(mix1, mix2, config)
+    % Plot (pressure, deflection) - shock polar
+    ax = set_fixed_figure(1, config);
+    
+    M1 = velocity_relative(mix1) / soundspeed(mix1);
+    plot(ax, [-flip(mix2.polar.theta), mix2.polar.theta], [flip(mix2.polar.p), mix2.polar.p] / pressure(mix1), 'k', 'LineWidth', config.linewidth);
+%     plot(ax, theta_max, mix2.polar.p(mix2.ind_max), 'ko', 'LineWidth', config.linewidth, 'MarkerFaceColor', 'auto');
+%     plot(ax, mix2.polar.theta(mix2.ind_sonic), mix2.polar.p(mix2.ind_sonic), 'ks', 'LineWidth', config.linewidth, 'MarkerFaceColor', 'auto');
+    text(ax, 0, max(mix2.polar.p), strcat('$\mathcal{M}_1 = ', sprintf('%.2g', M1), '$'), 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'center', 'FontSize', config.fontsize-4, 'Interpreter', 'latex');
+%     text(ax, theta_max, str2.polar.p2(str2.ind_max), 'detachment point \rightarrow m', 'HorizontalAlignment', 'right', 'FontSize', config.fontsize-6)
+%     text(ax, theta_sonic, str2.polar.p2(str2.ind_sonic), 'sonic point \rightarrow s', 'HorizontalAlignment', 'right', 'FontSize', config.fontsize-6)
+    xlabel(ax, 'Deflection angle $[^\circ]$', 'FontSize', config.fontsize, 'Interpreter', 'latex');
+    ylabel(ax, '$p_2/p_1$','FontSize', config.fontsize, 'Interpreter', 'latex');
+%     ylim(ax, [1, max(str2.polar.p2)]);
+end
+
+function ax = plot_shock_polar_wave(mix1, mix2, config)
+    % Plot (wave angle, deflection) - shock polar
+    ax = set_fixed_figure(2, config);
+    
+    M1 = velocity_relative(mix1) / soundspeed(mix1);
+    plot(ax, mix2.polar.theta, mix2.polar.beta, 'k', 'LineWidth', config.linewidth);  
+%     plot(ax, mix2.theta_max, mix2.polar.beta(str2.ind_max), 'ko', 'LineWidth', config.linewidth, 'MarkerFaceColor', 'auto');
+%     plot(ax, str2.polar.theta(str2.ind_sonic), mix2.polar.beta(str2.ind_sonic), 'ks', 'LineWidth', config.linewidth, 'MarkerFaceColor', 'auto');
+%     text(ax, mix2.theta_max, mix2.polar.beta(str2.ind_max),
+%     strcat('$\mathcal{M}_1 = ', sprintf('%.2g', M1), '$'), 'VerticalAlignment', 'middle', 'HorizontalAlignment', 'left', 'FontSize', config.fontsize-6, 'Interpreter', 'latex');
+    xlabel(ax, 'Deflection angle $[^\circ]$','FontSize', config.fontsize, 'Interpreter', 'latex');
+    ylabel(ax, 'Wave angle $[^\circ]$','FontSize', config.fontsize, 'Interpreter', 'latex');
+%     xlim(ax, [0, mix2.theta_max]);
+    ylim(ax, [0, 90]);
+end
+
+function ax = plot_shock_polar_velocities(mix1, mix2, config)
+    % Plot velocity components
+    ax = set_fixed_figure(3, config);
+
+    plot(ax, mix2.polar.ux, mix2.polar.uy, 'k', 'LineWidth', config.linewidth);  
+    xlabel(ax, '$u_{2x}$ [m/s]', 'FontSize', config.fontsize, 'Interpreter', 'latex');
+    ylabel(ax, '$u_{2y}$ [m/s]', 'FontSize', config.fontsize, 'Interpreter', 'latex');
+    %     ylim(ax, [0, max(mix2.polar.uy)]);
+end
\ No newline at end of file
diff --git a/Settings/functions/postResults.m b/Settings/functions/postResults.m
index 667a116a..44f27168 100644
--- a/Settings/functions/postResults.m
+++ b/Settings/functions/postResults.m
@@ -68,7 +68,16 @@ function postResults(self)
     plot_figure(self.PD.vP_vR.value, mix2,'v_P/v_R','T',self.Misc.config,self.PD.CompleteOrIncomplete);  
     T = cell2vector(mix2, 'T');
     displaysweepresults(self, mix2, T);
-    
+elseif strcmp(ProblemType,{'SHOCK_POLAR'})
+    % Shock polars
+    plot_shock_polar(self, mix1, mix2)
+%     % Incident velocity [m/s] againts molar fractions [-] 
+%     self.Misc.config.labelx = 'Incident velocity $u_1$ [m/s]';
+%     self.Misc.config.labely = 'Molar fraction $X_i$';
+%     for i = length(mix2):-1:1
+%         beta = cell2vector(mix2{i}.polar, 'beta');
+%         displaysweepresults(self, mix2{i}.polar, beta);
+%     end
 elseif strcmp(ProblemType,{'SHOCK_I'}) && length(phi) > 1
     self.Misc.config.labelx = 'Incident velocity $u_1$ [m/s]';
     self.Misc.config.labely = 'Temperature $T$ [K]';
diff --git a/Solver/Shocks and detonations/shock_polar.m b/Solver/Shocks and detonations/shock_polar.m
index 76abc576..437923ec 100644
--- a/Solver/Shocks and detonations/shock_polar.m	
+++ b/Solver/Shocks and detonations/shock_polar.m	
@@ -1,36 +1,30 @@
-function [str1, str2] = shock_polar(varargin)
+function [mix1, mix2] = shock_polar(varargin)
     % Solve oblique shock
 
     % Unpack input data
-    [self, str1, str2] = unpack(varargin);
-    % Abbreviations
-    config = self.Misc.config;
+    [self, mix1, mix2] = unpack(varargin);
     % Definitions
-    a1 = soundspeed(str1);
-    u1 = str1.u;
-    M1 = u1/a1;
+    a1 = soundspeed(mix1);
+    u1 = mix1.u;
     step = 5; % [m/s]
     u1n = linspace(a1 + 1, u1, (u1 - a1) / step);
     N = length(u1n);
     beta = asin(u1n ./ u1);
     beta_min = asin(a1 / u1);
     un = u1 .* cos(beta);
-
-    % Solve first case for initialization
-    [~, str2] = shock_incident(self, str1, u1n(end), str2);
     % Loop
     for i = N:-1:1
-        [~, str2] = shock_incident(self, str1, u1n(i), str2);
-        a2(i) = soundspeed(str2);
-        u2n(i) = str2.v_shock;
-        p2(i) = pressure(str2);
+        [~, mix2] = shock_incident(self, mix1, u1n(i), mix2);
+        a2(i) = soundspeed(mix2);
+        u2n(i) = mix2.v_shock;
+        p2(i) = pressure(mix2);
         theta(i) = beta(i) - atan(u2n(i) / un(i));
         u2(i) = u2n(i) * csc(beta(i) - theta(i));
+        Xi(:, i) = mix2.Xi;
     end
-
     % Velocity components downstream - laboratory fixed
-    u2x = u2 .* cos(theta);
-    u2y = u2 .* sin(theta);
+    u2x = u2 .* cos(theta); % [m/s]
+    u2y = u2 .* sin(theta); % [m/s]
     % Sonic point
     M2 = u2 ./ a2;
     ind_sonic = find(M2 > 1, 1, 'last');
@@ -38,52 +32,24 @@
     % Max deflection angle
     [theta_max, ind_max] = max(theta * 180/pi); % [deg]
     % Save results
-    str2.beta = beta(end) * 180/pi;             % [deg]
-    str2.theta = theta(end) * 180/pi;           % [deg]
-    str2.theta_max = theta_max;                 % [deg]
-    str2.beta_min = beta_min * 180/pi;          % [deg]
-    str2.beta_max = beta(ind_max) * 180/pi;     % [deg]
-    str2.theta_sonic = theta_sonic;             % [deg]
-    str2.beta_sonic = beta(ind_sonic) * 180/pi; % [deg]
-    
-    str2.theta_range = theta * 180/pi; % [deg]
-    str2.un = u2n;
-    % Plot (pressure, deflection) - shock polar
-	figure(1); ax = gca;
-    set(ax, 'LineWidth', config.linewidth, 'FontSize', config.fontsize-2, 'BoxStyle', 'full')
-    grid(ax, 'off'); box(ax, 'off'); hold(ax, 'on'); ax.Layer = 'Top'; axis(ax, 'tight')
-
-	plot(ax, [-flip(theta), theta] * 180/pi, [flip(p2), p2] / pressure(str1), 'k', 'LineWidth', config.linewidth);
-%     plot(ax, theta_max, p2(ind_max), 'ko', 'LineWidth', config.linewidth, 'MarkerFaceColor', 'auto');
-%     plot(ax, theta(ind_sonic) * 180/pi, p2(ind_sonic), 'ks', 'LineWidth', config.linewidth, 'MarkerFaceColor', 'auto');
-    text(ax, 0, max(p2), strcat('$\mathcal{M}_1 = ', sprintf('%.2g', M1), '$'), 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'center', 'FontSize', config.fontsize-4, 'Interpreter', 'latex')
-%     text(ax, theta_max, p2(ind_max), 'detachment point \rightarrow m', 'HorizontalAlignment', 'right', 'FontSize', config.fontsize-6)
-%     text(ax, theta_sonic * 180/pi, p2(ind_sonic), 'sonic point \rightarrow s', 'HorizontalAlignment', 'right', 'FontSize', config.fontsize-6)
-    xlabel(ax, 'Deflection angle $[^\circ]$', 'FontSize', config.fontsize, 'Interpreter', 'latex');
-    ylabel(ax, '$p_2/p_1$','FontSize', config.fontsize, 'Interpreter', 'latex');
-%     ylim(ax, [1, max(p2)]);
-    % Plot (wave angle, deflection) - shock polar
-    figure(2); ax = gca;
-    set(ax, 'LineWidth', config.linewidth, 'FontSize', config.fontsize-2, 'BoxStyle', 'full')
-    grid(ax, 'off'); box(ax, 'off'); hold(ax, 'on'); ax.Layer = 'Top'; axis(ax, 'tight');
-
-    plot(ax, theta * 180/pi, beta * 180/pi, 'k', 'LineWidth', config.linewidth);  
-%     plot(ax, theta_max, beta(ind_max) * 180/pi, 'ko', 'LineWidth', config.linewidth, 'MarkerFaceColor', 'auto');
-%     plot(ax, theta(ind_sonic) * 180/pi, beta(ind_sonic) * 180/pi, 'ks', 'LineWidth', config.linewidth, 'MarkerFaceColor', 'auto');
-%     text(ax, theta_max, beta(ind_max) * 180/pi, strcat('$\mathcal{M}_1 = ', sprintf('%.2g', M1), '$'), 'VerticalAlignment', 'middle', 'HorizontalAlignment', 'left', 'FontSize', config.fontsize-6, 'Interpreter', 'latex')
-    xlabel(ax, 'Deflection angle $[^\circ]$','FontSize', config.fontsize, 'Interpreter', 'latex');
-    ylabel(ax, 'Wave angle $[^\circ]$','FontSize', config.fontsize, 'Interpreter', 'latex');
-%     xlim(ax, [0, theta_max]);
-    ylim(ax, [0, 90]);
-    % Plot velocity components
-    figure(3); ax = gca;
-    set(ax, 'LineWidth', config.linewidth, 'FontSize', config.fontsize-2, 'BoxStyle', 'full')
-    grid(ax, 'off'); box(ax, 'off'); hold(ax, 'on'); ax.Layer = 'Top'; axis(ax, 'tight');
-
-    plot(ax, u2x, u2y, 'k', 'LineWidth', config.linewidth);  
-    xlabel(ax, '$u_{2x}$ [m/s]', 'FontSize', config.fontsize, 'Interpreter', 'latex');
-    ylabel(ax, '$u_{2y}$ [m/s]', 'FontSize', config.fontsize, 'Interpreter', 'latex');
-%     ylim(ax, [0, max(u2y)]);
+    mix2.beta = beta(end) * 180/pi;             % [deg]
+    mix2.theta = theta(end) * 180/pi;           % [deg]
+    mix2.theta_max = theta_max;                 % [deg]
+    mix2.beta_min = beta_min * 180/pi;          % [deg]
+    mix2.beta_max = beta(ind_max) * 180/pi;     % [deg]
+    mix2.theta_sonic = theta_sonic;             % [deg]
+    mix2.beta_sonic = beta(ind_sonic) * 180/pi; % [deg]
+    mix2.ind_max = ind_max;
+    mix2.ind_sonic = ind_sonic;
+    % Range values for the shock polar
+    mix2.polar.p = p2;    % [bar]
+    mix2.polar.Mach = M2; % [-]
+    mix2.polar.un = u2n;  % [m/s]
+    mix2.polar.ux = u2x;  % [m/s]
+    mix2.polar.uy = u2y;  % [m/s]
+    mix2.polar.beta = beta * 180/pi;   % [deg]
+    mix2.polar.theta = theta * 180/pi; % [deg]
+    mix2.polar.Xi = Xi;   % [-]
 end
 
 % SUB-PASS FUNCTIONS