%% 1.3 Piston Motion Analysis
% Ford Duratec 2.0 L
% Measured geometry from workshop data


%% Engine geometry

bore = 87.5e-3;          % Bore [m]
stroke = 83.4e-3;        % Stroke [m]
rodLength = 146.25e-3;   % Connecting rod length [m]
compressionRatio = 11.16;
numberOfCylinders = 4;

crankRadius = stroke / 2;

%% Crank angle range

thetaDeg = 0:1:360;
thetaRad = deg2rad(thetaDeg);

%% Piston position

% Piston displacement measured downwards from TDC
pistonPosition = crankRadius .* (1 - cos(thetaRad)) + ...
    rodLength - sqrt(rodLength.^2 - ...
    (crankRadius .* sin(thetaRad)).^2);

% Convert to mm
pistonPosition_mm = pistonPosition * 1000;

%% Cylinder volume

sweptVolume = (pi / 4) * bore^2 * stroke;
clearanceVolume = sweptVolume / (compressionRatio - 1);

cylinderVolume = clearanceVolume + ...
    (pi / 4) * bore^2 .* pistonPosition;

cylinderVolume_cc = cylinderVolume * 1e6;

%% Engine speed for velocity and acceleration calculations

engineSpeedRPM = 6000;
omega = 2 * pi * engineSpeedRPM / 60;

%% Piston velocity

% Differentiate position with respect to crank angle
dPosition_dTheta = gradient(pistonPosition, thetaRad);

% Convert angular derivative to time derivative
pistonVelocity = dPosition_dTheta * omega;

%% Piston acceleration

dVelocity_dTheta = gradient(pistonVelocity, thetaRad);
pistonAcceleration = dVelocity_dTheta * omega;

%% Mean and peak piston speed

meanPistonSpeed = 2 * stroke * engineSpeedRPM / 60;
peakPistonSpeed = max(abs(pistonVelocity));

fprintf('Engine speed: %.0f rpm\n', engineSpeedRPM);
fprintf('Mean piston speed: %.2f m/s\n', meanPistonSpeed);
fprintf('Peak piston speed: %.2f m/s\n', peakPistonSpeed);

%% Figure 1 - Piston position

figure;
plot(thetaDeg, pistonPosition_mm, 'LineWidth', 1.5);
grid on;
xlabel('Crank Angle (degrees)');
ylabel('Piston Position from TDC (mm)');
title('Piston Position vs Crank Angle');
xlim([0 360]);

%% Figure 2 - Piston velocity

figure;
plot(thetaDeg, pistonVelocity, 'LineWidth', 1.5);
grid on;
xlabel('Crank Angle (degrees)');
ylabel('Piston Velocity (m/s)');
title(sprintf('Piston Velocity vs Crank Angle at %.0f rpm', ...
    engineSpeedRPM));
xlim([0 360]);
yline(0, '--');

%% Figure 3 - Piston acceleration

figure;
plot(thetaDeg, pistonAcceleration, 'LineWidth', 1.5);
grid on;
xlabel('Crank Angle (degrees)');
ylabel('Piston Acceleration (m/s^2)');
title(sprintf('Piston Acceleration vs Crank Angle at %.0f rpm', ...
    engineSpeedRPM));
xlim([0 360]);
yline(0, '--');

%% Figure 4 - Cylinder volume

figure;
plot(thetaDeg, cylinderVolume_cc, 'LineWidth', 1.5);
grid on;
xlabel('Crank Angle (degrees)');
ylabel('Cylinder Volume (cm^3)');
title('Cylinder Volume vs Crank Angle');
xlim([0 360]);

%% Key calculated values

rodRatio = rodLength / crankRadius;
displacementTotal = sweptVolume * numberOfCylinders;

fprintf('\nEngine geometry results:\n');
fprintf('Crank radius: %.2f mm\n', crankRadius * 1000);
fprintf('Rod-to-crank ratio: %.3f\n', rodRatio);
fprintf('Swept volume per cylinder: %.2f cm^3\n', ...
    sweptVolume * 1e6);
fprintf('Total calculated displacement: %.1f cm^3\n', ...
    displacementTotal * 1e6);
fprintf('Clearance volume per cylinder: %.2f cm^3\n', ...
    clearanceVolume * 1e6);