%% Passive quarter car assignment - SRN 23103907
% AI-assisted code. Review and run in MATLAB before submission.
% State order: [xs xu vs vu]. SI units throughout.
clear; close all; clc;
p = [303.3 43.7 17628.6 178932.3 1311.7 5.2];
t = (0:0.001:8)';
[M,K,C,A] = matrices(p);
[V,D] = eig(K,M); [wn,j] = sort(sqrt(diag(D))); V = V(:,j);
V = V ./ V(1,:);
fprintf('Natural frequencies (Hz): %.6f %.6f\n',wn/(2*pi));
disp('Mode shapes normalised to body displacement = 1'); disp(V);
disp('Damped poles'); disp(eig(A));
[zf,af,lf,df,mf] = simulate(p,t,[0.029;0.022;0;0],0.5,true);
[zc,ac,lc,dc,mc] = simulate(p,t,[0.029;0.022;0;0],0.5,false);
[zb,ab,lb,db,mb] = simulate(p,t,zeros(4,1),0.5,false);
[r,vr] = road(t,0.5);
figure; tiledlayout(3,1);
nexttile; plot(t,1000*r); ylabel('Road (mm)'); grid on;
nexttile; plot(t,1000*zc(:,1:2)); ylabel('Displacement (mm)'); legend('Body','Wheel'); grid on;
nexttile; plot(t,1000*[dc zc(:,2)-r]); ylabel('Deflection (mm)'); xlabel('Time (s)'); legend('Suspension','Tyre'); grid on;
figure; tiledlayout(2,1);
nexttile; plot(t,1000*zf(:,1:2)); ylabel('Displacement (mm)'); legend('Body','Wheel'); grid on;
E = zeros(size(t));
for i=1:numel(t)
 E(i)=0.5*zf(i,3:4)*M*zf(i,3:4)'+0.5*zf(i,1:2)*K*zf(i,1:2)';
end
nexttile; plot(t,E); ylabel('Energy (J)'); xlabel('Time (s)'); grid on;
% One factor at a time sensitivity: all six parameters +/-20%.
sensitivity = zeros(6,1); low=zeros(6,7); high=low;
for j=1:6
 pl=p; ph=p; pl(j)=0.8*p(j); ph(j)=1.2*p(j);
 [~,~,~,~,low(j,:)]=simulate(pl,t,zeros(4,1),0.5,false);
 [~,~,~,~,high(j,:)]=simulate(ph,t,zeros(4,1),0.5,false);
 sensitivity(j)=(high(j,2)-low(j,2))/(0.4*mb(2));
end
names={'Ms';'Mu';'ks';'kt';'cs';'ct'};
disp(table(names,sensitivity));
figure; bar(sensitivity); xticks(1:6); xticklabels(names); ylabel('Normalised RMS sensitivity'); grid on;
% 144 point design grid. Minimise body RMS displacement with constraints.
best=inf; pm=[];
for k=linspace(0.6,1,9)
 for c=linspace(0.5,2,16)
  pp=p; pp(3)=k*p(3); pp(5)=c*p(5);
  [~,~,~,~,m]=simulate(pp,t,zeros(4,1),0.5,false);
  feasible=m(3)<=mb(3) && m(4)<=1.2*mb(4) && m(5)<=1.2*mb(5) && pp(1)*9.81/pp(3)<=0.25;
  if feasible && m(2)<best
   best=m(2); pm=pp;
  end
 end
end
assert(~isempty(pm),'No feasible grid candidate');
[zm,am,lm,dm,mm]=simulate(pm,t,zeros(4,1),0.5,false);
metricNames={'Peak body mm';'RMS body mm';'RMS acceleration m_s2';'Peak suspension mm';'Peak dynamic tyre load N';'Peak acceleration m_s2';'Peak wheel mm'};
disp(table(metricNames,mb',mm','VariableNames',{'Metric','Baseline','Modified'}));
disp('Modified parameters [Ms Mu ks kt cs ct]'); disp(pm);
figure; tiledlayout(3,1);
nexttile; plot(t,1000*[zb(:,1) zm(:,1)]); ylabel('Body (mm)'); legend('Baseline','Modified'); grid on;
nexttile; plot(t,[ab am]); ylabel('Acceleration (m/s^2)'); grid on;
nexttile; plot(t,[lb lm]); ylabel('Dynamic tyre load (N)'); xlabel('Time (s)'); grid on;
% Test road-duration dependence, retaining the 30 mm height.
for T=[0.25 0.5 1]
 [~,~,~,~,b]=simulate(p,t,zeros(4,1),T,false);
 [~,~,~,~,m]=simulate(pm,t,zeros(4,1),T,false);
 fprintf('T=%.2f s: body RMS %.6f -> %.6f mm; accel RMS %.6f -> %.6f m/s^2\n',T,b(2),m(2),b(3),m(3));
end
% Harmonic displacement transmissibility.
f=linspace(0.1,20,1000); H=zeros(numel(f),2);
for j=1:2
 if j==1, pp=p; else, pp=pm; end
 [M,K,C,~]=matrices(pp);
 for i=1:numel(f)
  w=2*pi*f(i); h=(K-w^2*M+1i*w*C)\[0;pp(4)+1i*w*pp(6)]; H(i,j)=abs(h(1));
 end
end
figure; plot(f,H); xlabel('Frequency (Hz)'); ylabel('|Xs/Xr|'); legend('Baseline','Modified'); grid on;
writetable(table(t,r,zc(:,1),zc(:,2),ac,dc,lc,'VariableNames',{'Time_s','Road_m','Body_m','Wheel_m','BodyAcceleration_m_s2','SuspensionTravel_m','DynamicTyreLoad_N'}),'Quarter_Car_23103907_results.csv');
function [M,K,C,A]=matrices(p)
 ms=p(1); mu=p(2); ks=p(3); kt=p(4); cs=p(5); ct=p(6);
 M=diag([ms mu]); K=[ks -ks;-ks ks+kt]; C=[cs -cs;-cs cs+ct];
 A=[zeros(2) eye(2);-M\K -M\C];
end
function [r,v]=road(t,T)
 r=zeros(size(t)); v=r; mask=t>=1 & t<=1+T;
 r(mask)=0.015*(1-cos(2*pi*(t(mask)-1)/T));
 v(mask)=0.03*pi/T*sin(2*pi*(t(mask)-1)/T);
end
function [z,acc,load,travel,m]=simulate(p,t,z0,T,free)
 [~,~,~,A]=matrices(p);
 opts=odeset('RelTol',1e-8,'AbsTol',1e-10,'MaxStep',0.01);
 [~,z]=ode45(@(q,y) rhs(q,y,A,p,T,free),t,z0,opts);
 if free, r=zeros(size(t)); v=r; else, [r,v]=road(t,T); end
 travel=z(:,1)-z(:,2);
 acc=-(p(3)*travel+p(5)*(z(:,3)-z(:,4)))/p(1);
 load=p(4)*(r-z(:,2))+p(6)*(v-z(:,4));
 a=t>=1; duration=t(end)-1;
 m=[max(abs(z(a,1)))*1000 sqrt(trapz(t(a),z(a,1).^2)/duration)*1000 sqrt(trapz(t(a),acc(a).^2)/duration) max(abs(travel(a)))*1000 max(abs(load(a))) max(abs(acc(a))) max(abs(z(a,2)))*1000];
end
function dy=rhs(t,y,A,p,T,free)
 if free, r=0; v=0; else, [r,v]=road(t,T); end
 dy=A*y+[0;0;0;(p(4)*r+p(6)*v)/p(2)];
end