%% 1.4 Camshaft and Valve Motion Analysis
% Ford Duratec 2.0 L
% Measured inlet and exhaust cam profiles


%% Import measured data

filename = 'Camshaft - Measurements.xlsx';

data = readtable(filename, ...
    'VariableNamingRule','preserve');

camAngleDeg = data{:,1};
inletRaw = data{:,2};
exhaustRaw = data{:,3};

% Remove missing or invalid rows
validRows = isfinite(camAngleDeg) & ...
            isfinite(inletRaw) & ...
            isfinite(exhaustRaw);

camAngleDeg = camAngleDeg(validRows);
inletRaw = inletRaw(validRows);
exhaustRaw = exhaustRaw(validRows);

%% Correct the angular reference

% Spreadsheet contains readings from 1 to 361 degrees.
% Convert this to a convenient 0 to 360 degree cam-angle range.
camAngleDeg = camAngleDeg - camAngleDeg(1);

camAngleRad = deg2rad(camAngleDeg);

%% Base-circle correction

% Use readings near the beginning and end of the test to estimate
% the zero-lift/base-circle reading.
numberBasePoints = 10;

inletBase = median([inletRaw(1:numberBasePoints); ...
                    inletRaw(end-numberBasePoints+1:end)]);

exhaustBase = median([exhaustRaw(1:numberBasePoints); ...
                      exhaustRaw(end-numberBasePoints+1:end)]);

inletLift = inletRaw - inletBase;
exhaustLift = exhaustRaw - exhaustBase;

% Negative lift values are caused by zero drift and measurement error.
% Physical valve lift cannot be negative.
inletLift(inletLift < 0) = 0;
exhaustLift(exhaustLift < 0) = 0;

%% Smooth measured lift data

% Numerical differentiation greatly amplifies small measurement noise.
% A short Savitzky-Golay smoothing window is therefore applied.
windowLength = 11;

inletLiftSmooth = smoothdata(inletLift, ...
    'sgolay',windowLength);

exhaustLiftSmooth = smoothdata(exhaustLift, ...
    'sgolay',windowLength);

% Prevent smoothing from producing small negative values
inletLiftSmooth(inletLiftSmooth < 0) = 0;
exhaustLiftSmooth(exhaustLiftSmooth < 0) = 0;

%% Engine operating speed

engineSpeedRPM = 6000;

% Camshaft rotates at half crankshaft speed
camSpeedRPM = engineSpeedRPM / 2;
omegaCam = 2*pi*camSpeedRPM/60;       % rad/s

%% Valve velocity

% First derivative of lift with respect to cam angle
dInlet_dTheta = gradient(inletLiftSmooth,camAngleRad);
dExhaust_dTheta = gradient(exhaustLiftSmooth,camAngleRad);

% Convert mm/rad to mm/s
inletVelocity = dInlet_dTheta * omegaCam;
exhaustVelocity = dExhaust_dTheta * omegaCam;

%% Valve acceleration

dInletVelocity_dTheta = gradient(inletVelocity,camAngleRad);
dExhaustVelocity_dTheta = gradient(exhaustVelocity,camAngleRad);

inletAcceleration = dInletVelocity_dTheta * omegaCam;
exhaustAcceleration = dExhaustVelocity_dTheta * omegaCam;

%% Valve jerk

dInletAcceleration_dTheta = ...
    gradient(inletAcceleration,camAngleRad);

dExhaustAcceleration_dTheta = ...
    gradient(exhaustAcceleration,camAngleRad);

inletJerk = dInletAcceleration_dTheta * omegaCam;
exhaustJerk = dExhaustAcceleration_dTheta * omegaCam;

%% Key results

[maxInletLift, inletPeakIndex] = max(inletLiftSmooth);
[maxExhaustLift, exhaustPeakIndex] = max(exhaustLiftSmooth);

inletPeakCamAngle = camAngleDeg(inletPeakIndex);
exhaustPeakCamAngle = camAngleDeg(exhaustPeakIndex);

fprintf('Camshaft profile results\n');
fprintf('-------------------------\n');
fprintf('Maximum inlet lift: %.3f mm\n',maxInletLift);
fprintf('Inlet peak cam angle: %.1f degrees\n', ...
    inletPeakCamAngle);

fprintf('Maximum exhaust lift: %.3f mm\n',maxExhaustLift);
fprintf('Exhaust peak cam angle: %.1f degrees\n', ...
    exhaustPeakCamAngle);

fprintf('\nKinematic calculations performed at %.0f rpm engine speed.\n', ...
    engineSpeedRPM);

fprintf('Peak inlet valve speed: %.1f mm/s\n', ...
    max(abs(inletVelocity)));

fprintf('Peak exhaust valve speed: %.1f mm/s\n', ...
    max(abs(exhaustVelocity)));

fprintf('Peak inlet acceleration: %.0f mm/s^2\n', ...
    max(abs(inletAcceleration)));

fprintf('Peak exhaust acceleration: %.0f mm/s^2\n', ...
    max(abs(exhaustAcceleration)));

%% Figure 7 - Valve lift

figure;
plot(camAngleDeg,inletLiftSmooth, ...
    'LineWidth',2);
hold on;
plot(camAngleDeg,exhaustLiftSmooth, ...
    'LineWidth',2);

grid on;
xlabel('Camshaft Angle (degrees)');
ylabel('Measured Cam Lift (mm)');
title('Inlet and Exhaust Cam Lift Profiles');
legend('Inlet','Exhaust','Location','best');
xlim([0 360]);

%% Figure 8 - Valve velocity

figure;
plot(camAngleDeg,inletVelocity, ...
    'LineWidth',2);
hold on;
plot(camAngleDeg,exhaustVelocity, ...
    'LineWidth',2);
yline(0,'--');

grid on;
xlabel('Camshaft Angle (degrees)');
ylabel('Valve Velocity (mm/s)');
title(sprintf('Valve Velocity at %.0f rpm Engine Speed', ...
    engineSpeedRPM));
legend('Inlet','Exhaust','Location','best');
xlim([0 360]);

%% Figure 9 - Valve acceleration

figure;
plot(camAngleDeg,inletAcceleration, ...
    'LineWidth',2);
hold on;
plot(camAngleDeg,exhaustAcceleration, ...
    'LineWidth',2);
yline(0,'--');

grid on;
xlabel('Camshaft Angle (degrees)');
ylabel('Valve Acceleration (mm/s^2)');
title(sprintf('Valve Acceleration at %.0f rpm Engine Speed', ...
    engineSpeedRPM));
legend('Inlet','Exhaust','Location','best');
xlim([0 360]);

%% Figure 10 - Valve jerk

figure;
plot(camAngleDeg,inletJerk, ...
    'LineWidth',2);
hold on;
plot(camAngleDeg,exhaustJerk, ...
    'LineWidth',2);
yline(0,'--');

grid on;
xlabel('Camshaft Angle (degrees)');
ylabel('Valve Jerk (mm/s^3)');
title(sprintf('Valve Jerk at %.0f rpm Engine Speed', ...
    engineSpeedRPM));
legend('Inlet','Exhaust','Location','best');
xlim([0 360]);