N = [2000:2000:12000]*(1/60);
y = cr+ct-(sqrt(cr.^2 - (ct.^2).*sin(CADrad).^2)+ct.*cos(CADrad));
V = V2 + 0.25.*pi.*(b.^2).*y;
P = zeros(length(N),length(Vc));
T = zeros(length(N),length(Vc));
Wcycle = zeros(1,length(N));
mdot = zeros(1,length(N));
Wnet = zeros(1,length(N));
T1s = Tatm*(PR)^((g-1)/g);
T1 = Tatm - ((Tatm-T1s)/etac);
M(count) = P(count,1)*Vc(1)/(R*T(count,1));
for count1 = 2:length(CADc(1:181))
P(count,count1) = P(count,count1-1)*(Vc(count1-1)/Vc(count1))^g;
T(count,count1) = P(count,count1)*Vc(count1)/(M(count)*R);
Mf(count) = M(count)/(1+FA);
Qin(count) = Mf(count)*LHV;
T(count,182) = T(count,181) + Qin(count)/(M(count)*cv);
P(count,182) = M(count)*R*T(count,182)/Vc(182);
for count2 = 183:length(CADc)
P(count,count2) = P(count,count2-1)*(Vc(count2-1)/Vc(count2))^g;
T(count,count2) = P(count,count2)*Vc(count2)/(M(count)*R);
Wcycle(count) = 2.*N(count).*trapz(Vc,P(count,:));
mdot(count) = 2.*M(count).*N(count);
Ws(count) = -mdot(count).*cps.*(T1-Tatm);
Wnet(count) = Wcycle(count) + Ws(count);
eta(count) = Wnet(count)./(2.*Qin(count)*N(count));
legend('No supercharger','Supercharger')
plot(N,Wnet*2/N,'-o','LineWidth',2)
xlabel('Supercharger Pressure Ratio')
plot(N,eta,'-o','LineWidth',2)
xlabel('Supercharger Pressure Ratio')
ylabel('Thermal Efficiency')