%% TASK 3 - Model Refinement and Validation
% Refined combustion model calibrated against measured pressure data

clear
clc
close all

%% --------------------------------------------------------
% IMPORT MEASURED DATA
%% --------------------------------------------------------

filename = 'Sigma 1.6L.xls';

opts = detectImportOptions(filename,'Sheet','6500');
opts.DataRange = 'A12';

T = readtable(filename,opts);

theta = T{:,1};
V_meas = T{:,2}*1e-6;
P_meas = T{:,3};

valid = isfinite(theta) & isfinite(V_meas) & isfinite(P_meas);

theta = theta(valid);
V_meas = V_meas(valid);
P_meas = P_meas(valid);

%% --------------------------------------------------------
% USE ONE ENGINE CYCLE
%% --------------------------------------------------------

cycle = theta >= -180 & theta <= 180;

Crank = theta(cycle);
Vol = V_meas(cycle);
P_exp = P_meas(cycle);

N = length(Crank);

%% --------------------------------------------------------
% ENGINE PARAMETERS
%% --------------------------------------------------------

gamma = 1.30;

nComp = 1.33;
nExp = 1.27;

LHV = 44e6;
rho = 1.18;
VE = 0.90;
AF = 14.7;
CombEff = 0.90;

R = 287;

%% --------------------------------------------------------
% CYLINDER MASS
%% --------------------------------------------------------

Vd = max(Vol)-min(Vol);

AirMass = rho*Vd*VE;
FuelMass = AirMass/AF;
TotalMass = AirMass + FuelMass;

%% --------------------------------------------------------
% INITIAL WIEBE PARAMETERS
%% --------------------------------------------------------

Spark = -12;
IgnDelay = 5;

SOC0 = Spark + IgnDelay;
Burn0 = 30;

a = 5;
m = 2;

%% --------------------------------------------------------
% CALIBRATE SOC AND BURN DURATION
%% --------------------------------------------------------

x0 = [SOC0 Burn0];

objective = @(x) pressureError(...
    x,Crank,Vol,P_exp,...
    FuelMass,LHV,CombEff,...
    gamma,nComp,nExp,a,m);

options = optimset('Display','off');

x = fminsearch(objective,x0,options);

SOC = x(1);
BurnDuration = x(2);

%% --------------------------------------------------------
% FINAL WIEBE FUNCTION
%% --------------------------------------------------------

Wiebe = zeros(N,1);

for i = 1:N

    if Crank(i) < SOC

        Wiebe(i)=0;

    elseif Crank(i) > SOC+BurnDuration

        Wiebe(i)=1;

    else

        y=(Crank(i)-SOC)/BurnDuration;

        Wiebe(i)=1-exp(-a*y^(m+1));

    end

end

MFB=[0;diff(Wiebe)];

%% --------------------------------------------------------
% CYLINDER PRESSURE
%% --------------------------------------------------------

P_model=zeros(N,1);

% Use measured initial pressure
P_model(1)=P_exp(1);

for i=1:N-1

    if Crank(i)<=0
        n=nComp;
    else
        n=nExp;
    end

    dQ=FuelMass*LHV*CombEff*MFB(i);

    Ppoly=P_model(i)*(Vol(i)/Vol(i+1))^n;

    Pheat=(gamma-1)*dQ/(Vol(i+1)*1e5);

    P_model(i+1)=Ppoly+Pheat;

end

%% --------------------------------------------------------
% CYLINDER TEMPERATURE
%% --------------------------------------------------------

Temp=(P_model*1e5.*Vol)/(TotalMass*R);

%% --------------------------------------------------------
% VALIDATION
%% --------------------------------------------------------

RMSE=sqrt(mean((P_model-P_exp).^2));

PeakMeasured=max(P_exp);
PeakModel=max(P_model);

PeakError=abs(PeakMeasured-PeakModel);
PeakErrorPercent=PeakError/PeakMeasured*100;

[~,i1]=max(P_exp);
[~,i2]=max(P_model);

PeakAngleMeasured=Crank(i1);
PeakAngleModel=Crank(i2);

TimingError=PeakAngleModel-PeakAngleMeasured;

Rvalue=corrcoef(P_model,P_exp);
Correlation=Rvalue(1,2);

%% --------------------------------------------------------
% RESULTS
%% --------------------------------------------------------

fprintf('\n');
fprintf('------------ TASK 3 RESULTS ------------\n');

fprintf('Optimised SOC            %.2f deg\n',SOC);
fprintf('Optimised Burn Duration  %.2f deg\n',BurnDuration);

fprintf('\n');

fprintf('Measured Peak Pressure   %.2f bar\n',PeakMeasured);
fprintf('Model Peak Pressure      %.2f bar\n',PeakModel);

fprintf('Peak Error               %.2f %%\n',PeakErrorPercent);

fprintf('Measured Peak Angle      %.1f deg\n',PeakAngleMeasured);
fprintf('Model Peak Angle         %.1f deg\n',PeakAngleModel);

fprintf('Timing Error             %.1f deg\n',TimingError);

fprintf('RMSE                     %.3f bar\n',RMSE);

fprintf('Correlation              %.4f\n',Correlation);

%% --------------------------------------------------------
% PLOTS
%% --------------------------------------------------------

figure
plot(Crank,P_exp,'k','LineWidth',2)
hold on
plot(Crank,P_model,'r','LineWidth',2)
grid on
xlabel('Crank Angle (deg)')
ylabel('Pressure (bar)')
title('Measured vs Refined Model')
legend('Measured','Refined Model','Location','Best')

figure
plot(Crank,P_model-P_exp,'LineWidth',2)
grid on
xlabel('Crank Angle (deg)')
ylabel('Pressure Error (bar)')
title('Pressure Difference')

figure
plot(Crank,Wiebe,'LineWidth',2)
grid on
xlabel('Crank Angle (deg)')
ylabel('Mass Fraction Burned')
title('Optimised Wiebe Function')

figure
plot(Crank,Temp,'LineWidth',2)
grid on
xlabel('Crank Angle (deg)')
ylabel('Temperature (K)')
title('Cylinder Temperature')

%% --------------------------------------------------------
% LOCAL FUNCTION
%% --------------------------------------------------------

function err = pressureError(x,Crank,Vol,Pexp,...
    FuelMass,LHV,CombEff,...
    gamma,nComp,nExp,a,m)

SOC=x(1);
Burn=x(2);

N=length(Crank);

Wiebe=zeros(N,1);

for i=1:N

    if Crank(i)<SOC

        Wiebe(i)=0;

    elseif Crank(i)>SOC+Burn

        Wiebe(i)=1;

    else

        y=(Crank(i)-SOC)/Burn;

        Wiebe(i)=1-exp(-a*y^(m+1));

    end

end

MFB=[0;diff(Wiebe)];

P=zeros(N,1);

P(1)=Pexp(1);

for i=1:N-1

    if Crank(i)<=0
        n=nComp;
    else
        n=nExp;
    end

    dQ=FuelMass*LHV*CombEff*MFB(i);

    P(i+1)=P(i)*(Vol(i)/Vol(i+1))^n + ...
        (gamma-1)*dQ/(Vol(i+1)*1e5);

end

err=sqrt(mean((P-Pexp).^2));

end