%% 1.5 Dynamic System Assessment

%% Effective reciprocating masses

mValveInlet = 48e-3;          % kg
mValveExhaust = 41e-3;        % kg

mFollowerInlet = 32.90e-3;    % kg
mFollowerExhaust = 32.41e-3;  % kg

mColletsInlet = 0.92e-3;      % kg
mColletsExhaust = 0.93e-3;    % kg

mSpring = 35e-3;              % kg

% One-third of the spring mass is assumed dynamically effective
mEffInlet = mValveInlet + mFollowerInlet + ...
    mColletsInlet + mSpring/3;

mEffExhaust = mValveExhaust + mFollowerExhaust + ...
    mColletsExhaust + mSpring/3;

%% Representative valve acceleration at 6000 rpm

% Use the 99th percentile to reduce the influence of isolated
% noise spikes created by numerical differentiation
aInletReference = prctile(abs(inletAcceleration),99) / 1000;
aExhaustReference = prctile(abs(exhaustAcceleration),99) / 1000;

referenceRPM = 6000;

%% Assumed available spring forces

springForce = [250 350 450];  % N

floatRPMInlet = referenceRPM .* sqrt( ...
    springForce ./ (mEffInlet .* aInletReference));

floatRPMExhaust = referenceRPM .* sqrt( ...
    springForce ./ (mEffExhaust .* aExhaustReference));

floatTable = table( ...
    springForce', ...
    floatRPMInlet', ...
    floatRPMExhaust', ...
    'VariableNames', ...
    {'AssumedSpringForce_N', ...
    'InletFloatRPM', ...
    'ExhaustFloatRPM'});

disp(floatTable);

fprintf('Effective inlet mass: %.2f g\n',mEffInlet*1000);
fprintf('Effective exhaust mass: %.2f g\n',mEffExhaust*1000);

%% Piston-to-valve timing assessment

crankAngle720 = 0:1:720;

% Repeat piston motion over the full four-stroke cycle
thetaPistonRad = deg2rad(mod(crankAngle720,360));

pistonFromTDC = crankRadius .* ...
    (1-cos(thetaPistonRad)) + ...
    rodLength - sqrt(rodLength.^2 - ...
    (crankRadius.*sin(thetaPistonRad)).^2);

pistonFromTDC_mm = pistonFromTDC * 1000;

%% Convert cam angle to crank angle

measuredCrankAngle = camAngleDeg * 2;

% These offsets must be changed once true valve timing is known
inletPhaseOffset = 0;     % degrees crank angle
exhaustPhaseOffset = 0;   % degrees crank angle

inletCrankAngle = mod(measuredCrankAngle + inletPhaseOffset,720);
exhaustCrankAngle = mod(measuredCrankAngle + exhaustPhaseOffset,720);

% Sort profiles after phase shifting
[inletCrankAngle, idxIn] = sort(inletCrankAngle);
inletLiftPhased = inletLiftSmooth(idxIn);

[exhaustCrankAngle, idxEx] = sort(exhaustCrankAngle);
exhaustLiftPhased = exhaustLiftSmooth(idxEx);

%% Interpolate valve lift onto the 0-720 degree crank grid

inletLift720 = interp1( ...
    inletCrankAngle,inletLiftPhased, ...
    crankAngle720,'linear',0);

exhaustLift720 = interp1( ...
    exhaustCrankAngle,exhaustLiftPhased, ...
    crankAngle720,'linear',0);

%% Plot piston position and valve lift

figure;
yyaxis left
plot(crankAngle720,pistonFromTDC_mm,'LineWidth',2);
ylabel('Piston Position from TDC (mm)');

yyaxis right
plot(crankAngle720,inletLift720,'LineWidth',2);
hold on;
plot(crankAngle720,exhaustLift720,'LineWidth',2);
ylabel('Valve Lift (mm)');

grid on;
xlabel('Crank Angle (degrees)');
title('Piston Position and Valve Lift over a Four-Stroke Cycle');
legend('Piston','Inlet Valve','Exhaust Valve','Location','best');
xlim([0 720]);

xline(0,'--','TDC');
xline(180,'--','BDC');
xline(360,'--','TDC');
xline(540,'--','BDC');
xline(720,'--','TDC');


%% Mass reduction required for 25 percent higher float RPM

targetRPMIncrease = 1.25;

massRatioRequired = 1 / targetRPMIncrease^2;
massReductionFraction = 1 - massRatioRequired;

newMassInlet = mEffInlet * massRatioRequired;
newMassExhaust = mEffExhaust * massRatioRequired;

massReductionInlet = mEffInlet - newMassInlet;
massReductionExhaust = mEffExhaust - newMassExhaust;

fprintf('\nMass reduction analysis\n');
fprintf('Required percentage reduction: %.1f %%\n', ...
    massReductionFraction*100);

fprintf('Original inlet effective mass: %.2f g\n', ...
    mEffInlet*1000);
fprintf('Target inlet effective mass: %.2f g\n', ...
    newMassInlet*1000);
fprintf('Required inlet reduction: %.2f g\n', ...
    massReductionInlet*1000);

fprintf('Original exhaust effective mass: %.2f g\n', ...
    mEffExhaust*1000);
fprintf('Target exhaust effective mass: %.2f g\n', ...
    newMassExhaust*1000);
fprintf('Required exhaust reduction: %.2f g\n', ...
    massReductionExhaust*1000);