% Thermodynamic Engine Modelling - Sigma 1.6L Data Analysis
% Author: Dominic Walton (Original) / Adapted for Sigma 1.6L
% University of Wales Trinity Saint David

clear all; close all; clc;

%% ==================== PART 1: IMPORT SIGMA 1.6L DATA ====================
fprintf('========================================\n');
fprintf('IMPORTING SIGMA 1.6L DATA FROM PAGE 4\n');
fprintf('========================================\n');

excel_file = 'Sigma 1.6L.xlsx';
sheet_number = 4;

% Check if file exists
if ~exist(excel_file, 'file')
    error('File "%s" not found. Please ensure it is in the current directory.', excel_file);
end

% Import data
try
    data_table = readtable(excel_file, 'Sheet', sheet_number);
    fprintf('Successfully imported %d rows and %d columns from Page 4.\n', ...
            height(data_table), width(data_table));
    
    % Display column names
    fprintf('\nAvailable columns on Page 4:\n');
    for col = 1:width(data_table)
        fprintf('  Column %d: %s\n', col, data_table.Properties.VariableNames{col});
    end
    
    % Show first few rows to help identify data
    fprintf('\nFirst 5 rows of data:\n');
    disp(head(data_table, 5));
    
catch ME
    error('Failed to import Excel file: %s', ME.message);
end

% ==================== USER INPUT FOR COLUMN MAPPING ====================
% IMPORTANT: Based on the column display above, you need to identify which
% columns contain your data. Modify these indices accordingly.

fprintf('\n========================================\n');
fprintf('COLUMN MAPPING REQUIRED\n');
fprintf('========================================\n');
fprintf('Please examine the column display above and update the indices below.\n');
fprintf('Typical engine data columns:\n');
fprintf('  - Crank angle (degrees)\n');
fprintf('  - Cylinder pressure (bar or Pa)\n');
fprintf('  - Cylinder volume (cc or m^3)\n');
fprintf('  - Temperature (K)\n');
fprintf('  - Heat release rate (J/deg)\n\n');

% === MODIFY THESE INDICES BASED ON YOUR DATA STRUCTURE ===
col_crank_angle = 1;        % Column containing crank angle [degrees]
col_pressure = 2;            % Column containing cylinder pressure [bar]
col_volume = 3;              % Column containing cylinder volume [cc or m^3]
col_temperature = 4;         % Column containing temperature [K] (if available)
col_heat_release = 5;        % Column containing heat release rate [J/deg] (if available)

% Initialize variables as empty in case columns don't exist
crank_given = [];
pressure_given = [];
volume_given = [];
volume_given_cc = [];
temperature_given = [];
heat_release_given = [];

% Extract data with error checking
try
    if width(data_table) >= col_crank_angle
        crank_given = data_table{:, col_crank_angle};
    end
    if width(data_table) >= col_pressure
        pressure_given = data_table{:, col_pressure};
    end
    if width(data_table) >= col_volume
        volume_given = data_table{:, col_volume};
        
        % Check units and convert if necessary
        if max(volume_given) < 1  % Likely in m^3
            volume_given_cc = volume_given * 1e6;  % Convert to cc for consistency
            fprintf('Volume data appears to be in m^3 - converted to cc\n');
        else
            volume_given_cc = volume_given;  % Already in cc
        end
    end
    
    fprintf('\n=== EXTRACTED DATA SUMMARY ===\n');
    if ~isempty(crank_given)
        fprintf('Crank angle range: %.1f to %.1f degrees\n', min(crank_given), max(crank_given));
    end
    if ~isempty(pressure_given)
        fprintf('Pressure range: %.2f to %.2f bar\n', min(pressure_given), max(pressure_given));
    end
    if ~isempty(volume_given_cc)
        fprintf('Volume range: %.2f to %.2f cc\n', min(volume_given_cc), max(volume_given_cc));
    end
    
    % Extract additional data if available
    if width(data_table) >= col_temperature
        temperature_given = data_table{:, col_temperature};
        fprintf('Temperature range: %.1f to %.1f K\n', min(temperature_given), max(temperature_given));
    end
    
    if width(data_table) >= col_heat_release
        heat_release_given = data_table{:, col_heat_release};
        fprintf('Heat release range: %.2f to %.2f J/deg\n', min(heat_release_given), max(heat_release_given));
    end
    
catch ME
    warning('Could not extract all columns: %s', ME.message);
    warning('Please adjust column indices based on the data structure shown above.');
end

%% ==================== PART 2: ENGINE GEOMETRY PARAMETERS ====================
fprintf('\n========================================\n');
fprintf('ENGINE GEOMETRY CALCULATIONS\n');
fprintf('========================================\n');

% Ford Engine Parameters (from original code)
Bore = 0.079;                    % meters
Stroke = 0.0814;                  % meters
Rod = 0.137;                       % meters
CR = 11.16;                         % Compression ratio

% Swept Volume
Vs = (pi*(Bore^2)/4) * Stroke;
fprintf('Swept Volume: %.6f m^3\n', Vs);

% Clearance Volume
Vc = Vs / (CR - 1);
fprintf('Clearance Volume: %.6f m^3\n', Vc);

% Crank Radius
Rad = Stroke / 2;
fprintf('Crank Radius: %.4f m\n', Rad);

% Crank Angle Range (for simulation)
Crank = -180:1:180;

% Piston Position Calculation
s = Rad * cosd(Crank) + sqrt(Rod^2 - (Rad * sind(Crank)).^2);
fprintf('Piston movement - Min: %.4f m, Max: %.4f m\n', min(s), max(s));

% Cylinder Volume Calculation
CVol = Vc + (pi*(Bore^2)/4) .* (Rad + Rod - s);
CVol_cc = CVol * 1e6;  % Convert to cc for comparison with given data
fprintf('Cylinder Volume - Min: %.2f cc, Max: %.2f cc\n', min(CVol_cc), max(CVol_cc));

%% ==================== PART 3: NON-COMBUSTION PRESSURE ====================
fprintf('\n========================================\n');
fprintf('NON-COMBUSTION PRESSURE CALCULATION\n');
fprintf('========================================\n');

% Parameters for non-combustion pressure
PolytropicCompression = 1.33;
PolytropicExpansion = 1.27;

% Calculate non-combustion pressure
P_nocomb = zeros(size(Crank));
P_nocomb(1) = 1;  % bar (ambient pressure)

for i = 1:length(Crank)-1
    if Crank(i) < 0
        n = PolytropicCompression;
    else
        n = PolytropicExpansion;
    end
    P_nocomb(i+1) = P_nocomb(i) * (CVol_cc(i)/CVol_cc(i+1))^n;
end

fprintf('Max pressure without combustion: %.2f bar\n', max(P_nocomb));

%% ==================== PART 4: IDEAL OTTO CYCLE ====================
fprintf('\n========================================\n');
fprintf('IDEAL OTTO CYCLE CALCULATION\n');
fprintf('========================================\n');

% Cylinder info
Air_Density = 1.2;               % kg/m^3
Gamma = 1.4;                      % Ratio of Specific Heats
Calorific = 4.40E+07;             % J/kg
Combustion_Eff = 0.95;
Specific_Heat_Capacity = 1007;    % J/kg.K
Stoich = 14.7;                     % Stoichiometric Fuel Ratio
AFR = 0.9;                         % Lambda

% Volume array for compression stroke (BDC to TDC)
Cvol_BDC_TDC = CVol_cc(181:-1:1);

% Compression stroke calculation
P1 = zeros(1,181);
T1 = zeros(1,181);
P1(1) = 1;                            % bar
T1(1) = 298;                           % K

for i = 1:180
    P1(i+1) = P1(i) * (Cvol_BDC_TDC(i)/Cvol_BDC_TDC(i+1))^Gamma;
    T1(i+1) = T1(i) * (Cvol_BDC_TDC(i)/Cvol_BDC_TDC(i+1))^(Gamma-1);
end

Temp2 = T1(end);
Pressure2 = P1(end);
fprintf('State 2 (TDC before combustion): P = %.2f bar, T = %.1f K\n', Pressure2, Temp2);

% Air and Fuel Mass
Air_Cvol = Air_Density * (Cvol_BDC_TDC(1) * 1e-6);  % Convert cc to m^3
AFRreal = Stoich * AFR;
FuelMass = Air_Cvol / AFRreal;
fprintf('Air mass: %.6f kg, Fuel mass: %.6f kg\n', Air_Cvol, FuelMass);

% Heat Energy Released
HeatEnergy = FuelMass * Calorific * Combustion_Eff;
fprintf('Heat released: %.2f J\n', HeatEnergy);

% Temperature Rise
CV_air = Specific_Heat_Capacity / Gamma;
DT = HeatEnergy / ((Air_Cvol + FuelMass) * Specific_Heat_Capacity);
fprintf('Temperature rise: %.1f K\n', DT);

% Expansion stroke calculation
Temp3 = Temp2 + DT;
Pressure3 = Pressure2 * (Temp3 / Temp2);
fprintf('State 3 (TDC after combustion): P = %.2f bar, T = %.1f K\n', Pressure3, Temp3);

P3 = zeros(1,181);
T3 = zeros(1,181);
P3(1) = Pressure3;
T3(1) = Temp3;

for i = 1:180
    P3(i+1) = P3(i) * (Cvol_BDC_TDC(i+1)/Cvol_BDC_TDC(i))^Gamma;
    T3(i+1) = T3(i) * (Cvol_BDC_TDC(i+1)/Cvol_BDC_TDC(i))^(Gamma-1);
end

Pressure4 = P3(end);
Temp4 = T3(end);
fprintf('State 4 (BDC before exhaust): P = %.2f bar, T = %.1f K\n', Pressure4, Temp4);

% Full cycle pressure
P_otto_full = [P1, P3(2:end)];
Crank_full = 0:1:360;

%% ==================== PART 5: WIEBE FUNCTION AND MFB ====================
fprintf('\n========================================\n');
fprintf('WIEBE FUNCTION CALCULATION\n');
fprintf('========================================\n');

% Wiebe Function Parameters
a = 5; m = 2;
SoC = -15;                % Start of Combustion
bdur = 40;                 % Burn Duration

% Calculate Wiebe cumulative function
Wiebe = 1 - exp(-a .* ((Crank - SoC)./bdur).^(m+1));

% Calculate Mass Fraction Burned rate
MFB = zeros(1, length(Crank)-1);
for i = 1:length(Crank)-1
    if Crank(i) < SoC
        MFB(i) = 0;
    else
        MFB(i) = (Wiebe(i+1) - Wiebe(i));
    end
end
MFB = [MFB, 0];  % Pad to match array size

fprintf('Peak MFB: %.6f\n', max(MFB));

%% ==================== PART 6: COMBUSTION MODELLING ====================
fprintf('\n========================================\n');
fprintf('COMBUSTION MODELLING\n');
fprintf('========================================\n');

% Parameters for combustion modelling
LHV = 44e6;
a_comb = 5; m_comb = 2;
Ignition = -12;
IgnDel = 5;
gamma_comb = 1.25;
bdur_comb = 30;
R_gas = 287;
rho_air = 1.167;
Combeff = 0.9;
Voleff = 0.9;
AF = 14.7;

% Derived parameters
Vd_comb = (pi/4) * Bore^2 * Stroke;
Air_comb = rho_air * Vd_comb * Voleff;
Fuel_comb = Air_comb / AF;
total_mass = Air_comb + Fuel_comb;
Pcom = 100000;  % Pa

% Non-Combustion Pressure (detailed)
P_detailed = zeros(1, length(Crank));
P_detailed(1) = Pcom;
for i = 1:length(Crank)-1
    if Crank(i) < 0
        n = 1.33;
    else
        n = 1.27;
    end
    P_detailed(i+1) = P_detailed(i) * (CVol(i)/CVol(i+1))^n;
end

% Wiebe Function for combustion
xb = zeros(1, length(Crank));
SoB = Ignition + IgnDel;
EoB = SoB + bdur_comb;
for i = 1:length(Crank)
    if Crank(i) < SoB
        xb(i) = 0;
    elseif Crank(i) <= EoB
        xb(i) = 1 - exp(-a_comb * ((Crank(i)-SoB)/bdur_comb)^(m_comb+1));
    else
        xb(i) = 1;
    end
end

% Heat Released
Qtol = Combeff * Fuel_comb * LHV;
dQr = zeros(1, length(Crank));
dxb = zeros(1, length(Crank));
for i = 1:length(Crank)
    if Crank(i) >= SoB && Crank(i) <= EoB
        theta_term = (Crank(i) - SoB) / bdur_comb;
        dxb(i) = (a_comb*(m_comb+1)/bdur_comb) * (theta_term^m_comb) * exp(-a_comb*(theta_term^(m_comb+1)));
        dQr(i) = Qtol * dxb(i);
    else
        dxb(i) = 0;
        dQr(i) = 0;
    end
end

% Combustion Pressure
Pcomb_sim = zeros(1, length(Crank));
Pcomb_sim(1) = P_detailed(1);
for i = 1:length(Crank)-1
    if Crank(i) < SoB
        Pcomb_sim(i+1) = P_detailed(i+1);
    elseif Crank(i) <= EoB
        dV = CVol(i+1) - CVol(i);
        dQ = dQr(i);
        dP = ((gamma_comb - 1)/CVol(i)) * dQ - (gamma_comb * Pcomb_sim(i)/CVol(i)) * dV;
        Pcomb_sim(i+1) = Pcomb_sim(i) + dP;
    else
        Pcomb_sim(i+1) = Pcomb_sim(i) * (CVol(i)/CVol(i+1))^1.27;
    end
end

% Convert combustion pressure to bar for comparison
Pcomb_sim_bar = Pcomb_sim / 100000;

% Temperature calculations
Temp_noncomb = (P_detailed .* CVol) / (total_mass * R_gas);
Temp_comb = (Pcomb_sim .* CVol) / (total_mass * R_gas);

fprintf('Simulated peak pressure (with combustion): %.2f bar\n', max(Pcomb_sim_bar));
fprintf('Simulated peak temperature: %.2f K\n', max(Temp_comb));

%% ==================== PART 7: DERIVED QUANTITIES ====================
fprintf('\n========================================\n');
fprintf('CALCULATING DERIVED QUANTITIES\n');
fprintf('========================================\n');

% Rate of pressure change
dP_dCrank = gradient(Pcomb_sim_bar, mean(diff(Crank)));

% Polytropic index calculation (from pressure-volume data)
n_polytropic = zeros(1, length(Crank)-1);
for i = 1:length(Crank)-1
    if Crank(i) < 0 && Crank(i) > -170  % Avoid endpoints
        n_polytropic(i) = log(Pcomb_sim_bar(i+1)/Pcomb_sim_bar(i)) / ...
                          log(CVol_cc(i)/CVol_cc(i+1));
    elseif Crank(i) > 0 && Crank(i) < 170
        n_polytropic(i) = log(Pcomb_sim_bar(i+1)/Pcomb_sim_bar(i)) / ...
                          log(CVol_cc(i)/CVol_cc(i+1));
    else
        n_polytropic(i) = NaN;
    end
end
n_polytropic = [n_polytropic, NaN];  % Pad array

% Torque calculation (simplified - assumes pressure acts on piston)
% FIXED: Using element-wise multiplication (.*) instead of matrix multiplication (*)
PistonArea_m2 = pi * (Bore^2) / 4;
rod_angle = asind((Rad/Rod) * sind(Crank));
torque_Nm = Pcomb_sim_bar .* 1e5 .* PistonArea_m2 .* Rad .* sind(Crank + rod_angle);

% For multiple RPMs (simulated - scaling factor)
RPM_range = [1000, 2000, 3000, 4000, 5000, 6000];
torque_rpm = zeros(length(RPM_range), length(Crank));
for rpm_idx = 1:length(RPM_range)
    % Simple scaling - torque reduces at higher RPM (simplified model)
    rpm_factor = 1 - 0.05 * (RPM_range(rpm_idx)/1000 - 1);
    torque_rpm(rpm_idx, :) = torque_Nm .* max(0.7, rpm_factor);
end

fprintf('Peak torque: %.1f Nm\n', max(torque_Nm));

%% ==================== PART 8: CREATE ALL GRAPHS ====================
fprintf('\n========================================\n');
fprintf('GENERATING GRAPHS\n');
fprintf('========================================\n');

% Graph 1: Cylinder Volume vs Crank Angle
figure(1)
clf
plot(Crank, CVol_cc, 'b-', 'LineWidth', 1.5)
hold on
if ~isempty(crank_given) && ~isempty(volume_given_cc)
    plot(crank_given, volume_given_cc, 'r--', 'LineWidth', 1.5)
    legend('Simulated', 'Given Data', 'Location', 'best')
else
    legend('Simulated', 'Location', 'best')
end
xlabel('Crank Angle (degrees)', 'FontSize', 11)
ylabel('Cylinder Volume (cc)', 'FontSize', 11)
title('Cylinder Volume vs Crank Angle', 'FontSize', 12)
grid on
grid minor

% Graph 2: Pressure vs Crank Angle
figure(2)
clf
plot(Crank, Pcomb_sim_bar, 'b-', 'LineWidth', 1.5)
hold on
plot(Crank, P_nocomb, 'g--', 'LineWidth', 1.5)
if ~isempty(crank_given) && ~isempty(pressure_given)
    % Ensure crank_given and pressure_given are same length and compatible
    if length(crank_given) == length(pressure_given)
        plot(crank_given, pressure_given, 'r:', 'LineWidth', 1.5)
        legend('Simulated (with combustion)', 'Without combustion', 'Given Data', 'Location', 'best')
    else
        warning('Given pressure data length mismatch - skipping plot');
        legend('Simulated (with combustion)', 'Without combustion', 'Location', 'best')
    end
else
    legend('Simulated (with combustion)', 'Without combustion', 'Location', 'best')
end
xlabel('Crank Angle (degrees)', 'FontSize', 11)
ylabel('Pressure (bar)', 'FontSize', 11)
title('Cylinder Pressure vs Crank Angle', 'FontSize', 12)
grid on
grid minor

% Graph 3: Pressure vs Volume (PV Diagram)
figure(3)
clf
plot(CVol_cc, Pcomb_sim_bar, 'b-', 'LineWidth', 1.5)
hold on
plot(CVol_cc, P_nocomb, 'g--', 'LineWidth', 1.5)
if ~isempty(volume_given_cc) && ~isempty(pressure_given)
    if length(volume_given_cc) == length(pressure_given)
        plot(volume_given_cc, pressure_given, 'r:', 'LineWidth', 1.5)
        legend('Simulated (with combustion)', 'Without combustion', 'Given Data', 'Location', 'best')
    else
        legend('Simulated (with combustion)', 'Without combustion', 'Location', 'best')
    end
else
    legend('Simulated (with combustion)', 'Without combustion', 'Location', 'best')
end
xlabel('Cylinder Volume (cc)', 'FontSize', 11)
ylabel('Pressure (bar)', 'FontSize', 11)
title('PV Diagram', 'FontSize', 12)
grid on
grid minor

% Graph 4: Log Pressure vs Log Volume
figure(4)
clf
loglog(CVol_cc, Pcomb_sim_bar, 'b-', 'LineWidth', 1.5)
hold on
loglog(CVol_cc, P_nocomb, 'g--', 'LineWidth', 1.5)
if ~isempty(volume_given_cc) && ~isempty(pressure_given)
    if length(volume_given_cc) == length(pressure_given)
        loglog(volume_given_cc, pressure_given, 'r:', 'LineWidth', 1.5)
        legend('Simulated (with combustion)', 'Without combustion', 'Given Data', 'Location', 'best')
    else
        legend('Simulated (with combustion)', 'Without combustion', 'Location', 'best')
    end
else
    legend('Simulated (with combustion)', 'Without combustion', 'Location', 'best')
end
xlabel('Cylinder Volume (cc) - log scale', 'FontSize', 11)
ylabel('Pressure (bar) - log scale', 'FontSize', 11)
title('Log Pressure vs Log Volume', 'FontSize', 12)
grid on

% Graph 5: Polytropic Index vs Crank Angle
figure(5)
clf
plot(Crank, n_polytropic, 'b-', 'LineWidth', 1.5)
xlabel('Crank Angle (degrees)', 'FontSize', 11)
ylabel('Polytropic Index (n)', 'FontSize', 11)
title('Polytropic Index vs Crank Angle', 'FontSize', 12)
ylim([1, 1.5])
grid on
grid minor
% Add reference lines
hold on
plot([-180 180], [1.33 1.33], 'r--', 'LineWidth', 1)
plot([-180 180], [1.27 1.27], 'g--', 'LineWidth', 1)
legend('Calculated n', 'Compression (1.33)', 'Expansion (1.27)', 'Location', 'best')

% Graph 6: Rate of Pressure Change vs Crank Angle
figure(6)
clf
plot(Crank, dP_dCrank, 'b-', 'LineWidth', 1.5)
xlabel('Crank Angle (degrees)', 'FontSize', 11)
ylabel('dP/dθ (bar/deg)', 'FontSize', 11)
title('Rate of Pressure Change vs Crank Angle', 'FontSize', 12)
grid on
grid minor

% Graph 7: Torque vs Crank Angle for different RPMs
figure(7)
clf
colors = lines(length(RPM_range));
for rpm_idx = 1:length(RPM_range)
    plot(Crank, torque_rpm(rpm_idx, :), 'Color', colors(rpm_idx, :), ...
         'LineWidth', 1.5, 'DisplayName', [num2str(RPM_range(rpm_idx)) ' RPM'])
    hold on
end
xlabel('Crank Angle (degrees)', 'FontSize', 11)
ylabel('Torque (Nm)', 'FontSize', 11)
title('Torque vs Crank Angle for Different Engine Speeds', 'FontSize', 12)
legend('Location', 'best')
grid on
grid minor

% Graph 8: MFB vs Crank Angle
figure(8)
clf
plot(Crank, MFB, 'b-', 'LineWidth', 1.5)
xlabel('Crank Angle (degrees)', 'FontSize', 11)
ylabel('Mass Fraction Burned Rate (1/deg)', 'FontSize', 11)
title('Mass Fraction Burned Rate vs Crank Angle', 'FontSize', 12)
grid on
grid minor
xlim([-30, 50])

% Graph 9: MFB vs Wiebe Function
figure(9)
clf
subplot(2,1,1)
plot(Crank, Wiebe, 'b-', 'LineWidth', 1.5)
xlabel('Crank Angle (degrees)', 'FontSize', 11)
ylabel('Wiebe Cumulative', 'FontSize', 11)
title('Wiebe Cumulative Function', 'FontSize', 12)
grid on
xlim([-30, 50])

subplot(2,1,2)
plot(Crank, MFB, 'r-', 'LineWidth', 1.5)
xlabel('Crank Angle (degrees)', 'FontSize', 11)
ylabel('MFB Rate', 'FontSize', 11)
title('Mass Fraction Burned Rate', 'FontSize', 12)
grid on
xlim([-30, 50])

% Graph 10: Burn Rate from Wiebe Function
figure(10)
clf
plot(Crank, dxb, 'b-', 'LineWidth', 1.5)
xlabel('Crank Angle (degrees)', 'FontSize', 11)
ylabel('Burn Rate (dxb/dθ)', 'FontSize', 11)
title('Burn Rate from Wiebe Function', 'FontSize', 12)
grid on
grid minor
xlim([-30, 50])

% Graph 11: Temperature from Combustion vs Crank Angle
figure(11)
clf
plot(Crank, Temp_comb, 'r-', 'LineWidth', 1.5)
hold on
plot(Crank, Temp_noncomb, 'b--', 'LineWidth', 1.5)
if ~isempty(crank_given) && ~isempty(temperature_given)
    if length(crank_given) == length(temperature_given)
        plot(crank_given, temperature_given, 'g:', 'LineWidth', 1.5)
        legend('With Combustion', 'Without Combustion', 'Given Data', 'Location', 'best')
    else
        legend('With Combustion', 'Without Combustion', 'Location', 'best')
    end
else
    legend('With Combustion', 'Without Combustion', 'Location', 'best')
end
xlabel('Crank Angle (degrees)', 'FontSize', 11)
ylabel('Temperature (K)', 'FontSize', 11)
title('Cylinder Temperature vs Crank Angle', 'FontSize', 12)
grid on
grid minor

% Graph 12: Pressure without Combustion
figure(12)
clf
plot(Crank, P_nocomb, 'b-', 'LineWidth', 1.5)
xlabel('Crank Angle (degrees)', 'FontSize', 11)
ylabel('Pressure (bar)', 'FontSize', 11)
title('Pressure without Combustion', 'FontSize', 12)
grid on
grid minor

% Graph 13: Simulated Pressure vs Crank Angle (all versions)
figure(13)
clf
plot(Crank, P_nocomb, 'g--', 'LineWidth', 1.5, 'DisplayName', 'Simple Non-Combustion')
hold on
plot(Crank, Pcomb_sim_bar, 'b-', 'LineWidth', 1.5, 'DisplayName', 'Detailed Combustion')
plot(Crank_full, P_otto_full, 'm:', 'LineWidth', 1.5, 'DisplayName', 'Ideal Otto Cycle')
xlabel('Crank Angle (degrees)', 'FontSize', 11)
ylabel('Pressure (bar)', 'FontSize', 11)
title('Comparison of Simulated Pressures', 'FontSize', 12)
legend('Location', 'best')
grid on
grid minor

% Graph 14: Simulated Pressure with Combustion vs Given Pressure
figure(14)
clf
if ~isempty(crank_given) && ~isempty(pressure_given) && length(crank_given) == length(pressure_given)
    plot(Crank, Pcomb_sim_bar, 'b-', 'LineWidth', 1.5, 'DisplayName', 'Simulated with Combustion')
    hold on
    plot(crank_given, pressure_given, 'r--', 'LineWidth', 1.5, 'DisplayName', 'Given Data')
    
    % Calculate error
    % Interpolate simulated data to match given crank angles
    P_sim_interp = interp1(Crank, Pcomb_sim_bar, crank_given, 'linear', 'extrap');
    error_rms = sqrt(mean((P_sim_interp - pressure_given).^2));
    error_max = max(abs(P_sim_interp - pressure_given));
    
    xlabel('Crank Angle (degrees)', 'FontSize', 11)
    ylabel('Pressure (bar)', 'FontSize', 11)
    title(['Simulated vs Given Pressure (RMS Error: ' num2str(error_rms, '%.2f') ' bar)'], 'FontSize', 12)
    legend('Location', 'best')
    grid on
    grid minor
    
    fprintf('\n=== MODEL VALIDATION ===\n');
    fprintf('RMS Error vs given data: %.2f bar\n', error_rms);
    fprintf('Max Error vs given data: %.2f bar\n', error_max);
else
    plot(0, 0, 'w')
    text(0.5, 0.5, 'Given pressure data not available or length mismatch', ...
         'HorizontalAlignment', 'center', 'FontSize', 14)
    xlim([0 1]); ylim([0 1]);
    axis off
    title('Simulated vs Given Pressure - No Data Available', 'FontSize', 12)
end

% Graph 15: Comparison of all pressure models
figure(15)
clf
subplot(2,2,1)
plot(Crank, P_nocomb, 'b-', 'LineWidth', 1.5)
xlabel('Crank Angle (deg)')
ylabel('Pressure (bar)')
title('Without Combustion')
grid on

subplot(2,2,2)
plot(Crank, Pcomb_sim_bar, 'r-', 'LineWidth', 1.5)
xlabel('Crank Angle (deg)')
ylabel('Pressure (bar)')
title('With Combustion (Simulated)')
grid on

subplot(2,2,3)
if ~isempty(crank_given) && ~isempty(pressure_given) && length(crank_given) == length(pressure_given)
    plot(crank_given, pressure_given, 'g-', 'LineWidth', 1.5)
    xlabel('Crank Angle (deg)')
    ylabel('Pressure (bar)')
    title('Given Data')
    grid on
else
    text(0.5, 0.5, 'No given data or length mismatch', 'HorizontalAlignment', 'center')
    xlim([0 1]); ylim([0 1]);
    axis off
end

subplot(2,2,4)
plot(Crank, P_nocomb, 'b--', 'LineWidth', 1.5)
hold on
plot(Crank, Pcomb_sim_bar, 'r-', 'LineWidth', 1.5)
if ~isempty(crank_given) && ~isempty(pressure_given) && length(crank_given) == length(pressure_given)
    plot(crank_given, pressure_given, 'g:', 'LineWidth', 1.5)
    legend('No Comb', 'Sim Comb', 'Given', 'Location', 'best')
else
    legend('No Comb', 'Sim Comb', 'Location', 'best')
end
xlabel('Crank Angle (deg)')
ylabel('Pressure (bar)')
title('All Comparisons')
grid on

sgtitle('Comprehensive Pressure Model Comparison', 'FontSize', 14, 'FontWeight', 'bold')

%% ==================== PART 9: SET FIGURE NAMES FOR WORD DOCUMENT ====================
% Set figure names for easy identification in the Word document
figure(1); set(gcf, 'Name', 'Cylinder Volume vs Crank Angle');
figure(2); set(gcf, 'Name', 'Pressure vs Crank Angle');
figure(3); set(gcf, 'Name', 'PV Diagram');
figure(4); set(gcf, 'Name', 'Log Pressure vs Log Volume');
figure(5); set(gcf, 'Name', 'Polytropic Index vs Crank Angle');
figure(6); set(gcf, 'Name', 'Rate of Pressure Change');
figure(7); set(gcf, 'Name', 'Torque vs Crank Angle (Multiple RPMs)');
figure(8); set(gcf, 'Name', 'MFB vs Crank Angle');
figure(9); set(gcf, 'Name', 'MFB vs Wiebe Function');
figure(10); set(gcf, 'Name', 'Burn Rate from Wiebe');
figure(11); set(gcf, 'Name', 'Temperature vs Crank Angle');
figure(12); set(gcf, 'Name', 'Pressure without Combustion');
figure(13); set(gcf, 'Name', 'Simulated Pressure Comparison');
figure(14); set(gcf, 'Name', 'Simulated vs Given Pressure');
figure(15); set(gcf, 'Name', 'Comprehensive Pressure Comparison');

%% ==================== CREATE WORD DOCUMENT WITH GRAPHS ONLY ====================
fprintf('\n========================================\n');
fprintf('CREATING WORD DOCUMENT WITH GRAPHS ONLY\n');
fprintf('========================================\n');

try
    % Create Word object
    Word = actxserver('Word.Application');
    Word.Visible = true;  % Set to false if you want it to run in background
    Document = Word.Documents.Add();
    
    % Set page margins
    Selection = Word.Selection;
    Selection.PageSetup.TopMargin = 36;      % 0.5 inch
    Selection.PageSetup.BottomMargin = 36;
    Selection.PageSetup.LeftMargin = 36;
    Selection.PageSetup.RightMargin = 36;
    
    % ==================== TITLE PAGE (SIMPLE) ====================
    Selection.ParagraphFormat.Alignment = 1;  % Center align
    
    % Add some space from top
    for i = 1:5
        Selection.TypeParagraph();
    end
    
    % Main title
    Selection.Font.Size = 24;
    Selection.Font.Bold = 1;
    Selection.Font.Name = 'Arial';
    Selection.TypeText('Sigma 1.6L Engine Analysis');
    Selection.TypeParagraph();
    Selection.TypeParagraph();
    
    % Subtitle
    Selection.Font.Size = 14;
    Selection.Font.Bold = 0;
    Selection.TypeText('Graphical Results');
    Selection.TypeParagraph();
    Selection.TypeParagraph();
    Selection.TypeParagraph();
    
    % Date
    Selection.Font.Size = 11;
    Selection.TypeText(['Generated: ' datestr(now)]);
    Selection.TypeParagraph();
    
    % Reset alignment
    Selection.ParagraphFormat.Alignment = 0;  % Left align
    
    % Page break
    Selection.InsertBreak(7);  % 7 = wdPageBreak
    
    % ==================== GRAPHS ONLY ====================
    graph_titles = {
        'Figure 1: Cylinder Volume vs Crank Angle',
        'Figure 2: Pressure vs Crank Angle',
        'Figure 3: PV Diagram',
        'Figure 4: Log Pressure vs Log Volume',
        'Figure 5: Polytropic Index vs Crank Angle',
        'Figure 6: Rate of Pressure Change',
        'Figure 7: Torque vs Crank Angle (Multiple RPMs)',
        'Figure 8: MFB vs Crank Angle',
        'Figure 9: MFB vs Wiebe Function',
        'Figure 10: Burn Rate from Wiebe',
        'Figure 11: Temperature vs Crank Angle',
        'Figure 12: Pressure without Combustion',
        'Figure 13: Simulated Pressure Comparison',
        'Figure 14: Simulated vs Given Pressure',
        'Figure 15: Comprehensive Pressure Comparison'
    };
    
    for fig_num = 1:15
        % Center the figure heading
        Selection.ParagraphFormat.Alignment = 1;  % Center
        
        % Figure number and title
        Selection.Font.Size = 14;
        Selection.Font.Bold = 1;
        Selection.Font.Name = 'Arial';
        Selection.TypeText(graph_titles{fig_num});
        Selection.Font.Bold = 0;
        Selection.TypeParagraph();
        Selection.TypeParagraph();
        
        % Activate and copy figure
        figure(fig_num)
        set(gcf, 'PaperPositionMode', 'auto')
        set(gcf, 'Position', [100, 100, 800, 600])
        print -clipboard -dmeta
        
        % Paste the figure (centered)
        Selection.Paste();
        Selection.TypeParagraph();
        Selection.TypeParagraph();
        
        % Reset alignment for next heading
        Selection.ParagraphFormat.Alignment = 1;
        
        % Add page break after each graph (except the last one)
        if fig_num < 15
            Selection.InsertBreak(7);  % Page break
        end
    end
    
    % Add a simple footer with page numbers
    for i = 1:Document.Sections.Count
        Footer = Document.Sections.Item(i).Footers.Item(1);
        Footer.Range.Font.Size = 9;
        Footer.Range.Font.Italic = 1;
        Footer.Range.Font.Name = 'Arial';
        Footer.Range.Text = 'Sigma 1.6L Engine Analysis - Graphical Results';
        Footer.PageNumbers.Add();
    end
    
    % Save the document
    filename = fullfile(pwd, 'Sigma_1.6L_Graphs_Only.docx');
    Document.SaveAs2(filename);
    
    fprintf('\n✅ Word document created with 15 graphs only!\n');
    fprintf('📍 File location: %s\n', filename);
    fprintf('📍 Document contains: Title page + 15 full-page graphs\n');
    
catch ME
    fprintf('\n❌ Error creating Word document: %s\n', ME.message);
    fprintf('Make sure Microsoft Word is installed on your computer.\n');
end

%% ==================== CREATE HTML REPORT WITH GRAPHS AND TEXT BOXES ====================
fprintf('\n========================================\n');
fprintf('CREATING HTML REPORT WITH GRAPHS AND TEXT BOXES\n');
fprintf('========================================\n');

try
    % Create a figures directory
    if ~exist('html_figures', 'dir')
        mkdir('html_figures');
    end
    
    % Save all figures as PNG files for the HTML report
    fprintf('Saving figures for HTML report...\n');
    for fig_num = 1:15
        figure(fig_num)
        % Make it look nice
        set(gcf, 'PaperPositionMode', 'auto')
        set(gcf, 'Position', [100, 100, 800, 600])
        % Save as PNG
        filename = sprintf('html_figures\\figure_%02d.png', fig_num);
        saveas(gcf, filename);
        fprintf('  Saved figure %d\n', fig_num);
    end
    
    % Create filename for the HTML report
    html_filename = 'Sigma_1.6L_Engine_Analysis_Report.html';
    
    % Open file for writing
    fileID = fopen(html_filename, 'w');
    
    % Write HTML header with styling
    fprintf(fileID, '<!DOCTYPE html>\n');
    fprintf(fileID, '<html>\n');
    fprintf(fileID, '<head>\n');
    fprintf(fileID, '    <title>Sigma 1.6L Engine Analysis Results</title>\n');
    fprintf(fileID, '    <style>\n');
    fprintf(fileID, '        body { font-family: Arial, sans-serif; margin: 40px; background-color: #f5f5f5; }\n');
    fprintf(fileID, '        .container { max-width: 1200px; margin: 0 auto; background-color: white; padding: 30px; box-shadow: 0 0 20px rgba(0,0,0,0.1); border-radius: 10px; }\n');
    fprintf(fileID, '        h1 { color: #003366; border-bottom: 3px solid #003366; padding-bottom: 10px; }\n');
    fprintf(fileID, '        h2 { color: #006699; margin-top: 40px; background-color: #e6f3ff; padding: 10px; border-radius: 5px; }\n');
    fprintf(fileID, '        .figure { margin: 40px 0; text-align: center; border: 1px solid #ddd; padding: 20px; border-radius: 5px; }\n');
    fprintf(fileID, '        .figure img { max-width: 900px; border: 1px solid #ccc; border-radius: 5px; }\n');
    fprintf(fileID, '        .notes-box { margin-top: 20px; text-align: left; border: 2px dashed #006699; padding: 15px; background-color: #f9f9f9; border-radius: 5px; }\n');
    fprintf(fileID, '        .top-notes { margin: 30px 0; text-align: left; border: 3px solid #003366; padding: 20px; background-color: #ffffff; border-radius: 10px; }\n');
    fprintf(fileID, '        textarea { width: 100%%; padding: 10px; border: 1px solid #ccc; border-radius: 5px; font-family: Arial, sans-serif; font-size: 14px; box-sizing: border-box; }\n');
    fprintf(fileID, '        .top-textarea { width: 100%%; height: 200px; padding: 15px; border: 2px solid #003366; border-radius: 5px; font-family: Arial, sans-serif; font-size: 14px; }\n');
    fprintf(fileID, '        .btn { padding: 15px 30px; margin: 10px; border: none; border-radius: 5px; cursor: pointer; font-size: 16px; font-weight: bold; }\n');
    fprintf(fileID, '        .btn-save { background-color: #28a745; color: white; }\n');
    fprintf(fileID, '        .btn-load { background-color: #17a2b8; color: white; }\n');
    fprintf(fileID, '        .btn-clear { background-color: #dc3545; color: white; }\n');
    fprintf(fileID, '        .btn-print { background-color: #6c757d; color: white; }\n');
    fprintf(fileID, '        .footer { text-align: center; margin-top: 50px; padding-top: 20px; border-top: 1px solid #ddd; color: #666; }\n');
    fprintf(fileID, '    </style>\n');
    
    % Add JavaScript for saving/loading notes
    fprintf(fileID, '    <script>\n');
    fprintf(fileID, '        function saveNotes() {\n');
    fprintf(fileID, '            var notes = {};\n');
    fprintf(fileID, '            var top = document.getElementById("top_notes");\n');
    fprintf(fileID, '            if (top) notes["TopNotes"] = top.value;\n');
    fprintf(fileID, '            for (var i = 1; i <= 15; i++) {\n');
    fprintf(fileID, '                var textarea = document.getElementById("notes_" + i);\n');
    fprintf(fileID, '                if (textarea) notes["Figure " + i] = textarea.value;\n');
    fprintf(fileID, '            }\n');
    fprintf(fileID, '            var dataStr = "data:text/json;charset=utf-8," + encodeURIComponent(JSON.stringify(notes, null, 2));\n');
    fprintf(fileID, '            var a = document.createElement("a"); a.href = dataStr; a.download = "Sigma_1.6L_Notes.json"; a.click();\n');
    fprintf(fileID, '            alert("All notes saved successfully!");\n');
    fprintf(fileID, '        }\n');
    fprintf(fileID, '        \n');
    fprintf(fileID, '        function loadNotes() {\n');
    fprintf(fileID, '            var input = document.createElement("input"); input.type = "file"; input.accept = ".json";\n');
    fprintf(fileID, '            input.onchange = function(e) {\n');
    fprintf(fileID, '                var reader = new FileReader();\n');
    fprintf(fileID, '                reader.onload = function(readerEvent) {\n');
    fprintf(fileID, '                    var notes = JSON.parse(readerEvent.target.result);\n');
    fprintf(fileID, '                    var top = document.getElementById("top_notes");\n');
    fprintf(fileID, '                    if (top && notes["TopNotes"]) top.value = notes["TopNotes"];\n');
    fprintf(fileID, '                    for (var key in notes) {\n');
    fprintf(fileID, '                        if (key.startsWith("Figure ")) {\n');
    fprintf(fileID, '                            var figNum = key.split(" ")[1];\n');
    fprintf(fileID, '                            var textarea = document.getElementById("notes_" + figNum);\n');
    fprintf(fileID, '                            if (textarea) textarea.value = notes[key];\n');
    fprintf(fileID, '                        }\n');
    fprintf(fileID, '                    }\n');
    fprintf(fileID, '                    alert("All notes loaded successfully!");\n');
    fprintf(fileID, '                };\n');
    fprintf(fileID, '                reader.readAsText(e.target.files[0]);\n');
    fprintf(fileID, '            };\n');
    fprintf(fileID, '            input.click();\n');
    fprintf(fileID, '        }\n');
    fprintf(fileID, '        \n');
    fprintf(fileID, '        function clearAllNotes() {\n');
    fprintf(fileID, '            if (confirm("Clear ALL notes?")) {\n');
    fprintf(fileID, '                var top = document.getElementById("top_notes");\n');
    fprintf(fileID, '                if (top) top.value = "";\n');
    fprintf(fileID, '                for (var i = 1; i <= 15; i++) {\n');
    fprintf(fileID, '                    var textarea = document.getElementById("notes_" + i);\n');
    fprintf(fileID, '                    if (textarea) textarea.value = "";\n');
    fprintf(fileID, '                }\n');
    fprintf(fileID, '                alert("All notes cleared!");\n');
    fprintf(fileID, '            }\n');
    fprintf(fileID, '        }\n');
    fprintf(fileID, '        \n');
    fprintf(fileID, '        function printReport() { window.print(); }\n');
    fprintf(fileID, '    </script>\n');
    fprintf(fileID, '</head>\n');
    fprintf(fileID, '<body>\n');
    fprintf(fileID, '    <div class="container">\n');
    
    % Title and control buttons
    fprintf(fileID, '        <h1>Sigma 1.6L Engine Analysis Results</h1>\n');
    fprintf(fileID, '        <p><strong>Generated on:</strong> %s</p>\n', datestr(now));
    fprintf(fileID, '        <div style="text-align: center;">\n');
    fprintf(fileID, '            <button class="btn btn-save" onclick="saveNotes()">💾 Save Notes</button>\n');
    fprintf(fileID, '            <button class="btn btn-load" onclick="loadNotes()">📂 Load Notes</button>\n');
    fprintf(fileID, '            <button class="btn btn-clear" onclick="clearAllNotes()">🗑️ Clear Notes</button>\n');
    fprintf(fileID, '            <button class="btn btn-print" onclick="printReport()">🖨️ Print</button>\n');
    fprintf(fileID, '        </div>\n\n');
    
    % BLANK TOP TEXT BOX
    fprintf(fileID, '        <div class="top-notes">\n');
    fprintf(fileID, '            <textarea id="top_notes" class="top-textarea" placeholder=""></textarea>\n');
    fprintf(fileID, '        </div>\n\n');
    
    % Graphs with text boxes
    graph_titles = {
        'Cylinder Volume vs Crank Angle',
        'Pressure vs Crank Angle',
        'PV Diagram',
        'Log Pressure vs Log Volume',
        'Polytropic Index vs Crank Angle',
        'Rate of Pressure Change',
        'Torque vs Crank Angle (Multiple RPMs)',
        'MFB vs Crank Angle',
        'MFB vs Wiebe Function',
        'Burn Rate from Wiebe',
        'Temperature vs Crank Angle',
        'Pressure without Combustion',
        'Simulated Pressure Comparison',
        'Simulated vs Given Pressure',
        'Comprehensive Pressure Comparison'
    };
    
    for fig_num = 1:15
        fprintf(fileID, '        <div class="figure">\n');
        fprintf(fileID, '            <h3>Figure %d: %s</h3>\n', fig_num, graph_titles{fig_num});
        fprintf(fileID, '            <img src="html_figures/figure_%02d.png" alt="%s">\n', fig_num, graph_titles{fig_num});
        fprintf(fileID, '            <div class="notes-box">\n');
        fprintf(fileID, '                <textarea id="notes_%d" placeholder=""></textarea>\n', fig_num);
        fprintf(fileID, '            </div>\n');
        fprintf(fileID, '        </div>\n');
    end
    
    % Footer
    fprintf(fileID, '        <div class="footer">\n');
    fprintf(fileID, '            <p>Generated by MATLAB on %s</p>\n', datestr(now));
    fprintf(fileID, '        </div>\n');
    fprintf(fileID, '    </div>\n');
    fprintf(fileID, '</body>\n');
    fprintf(fileID, '</html>\n');
    
    fclose(fileID);
    
    fprintf('\n✅ HTML report created: %s\n', html_filename);
    fprintf('✅ 15 graphs saved in "html_figures" folder\n');
    fprintf('✅ Blank text box at top + text boxes under each graph\n');
    
    % Open in browser
    fprintf('\n🚀 Opening HTML report in your browser...\n');
    web(html_filename, '-browser');
    
catch ME
    fprintf('\n❌ Error creating HTML report: %s\n', ME.message);
end