close all
clear all

gamma = 1.4;

V0 = 0.5*231;  %in^3

V1 = (0.5-.15)*231;  

Aa = pi*3^2;  %6 Inch diameter

lo = V0/Aa;
l1 = (Aa*lo-V1)/Aa;

Ap = (pi/4)*(2.5^2-1.75^2);

z1 = 0;  %Min stroke
z2 = 25;  %Max stroke
z = z1:(z2-z1)/100:z2;

F1 = 100;
P1 = F1/Ap;  %Precharge (psi)


l = Ap*(z-z1)/Aa+l1;
V = Aa*(lo-l);

CR = V1./V;


Piso = P1*(lo-l1)./(lo-l1-Ap*(z-z1)/Aa);
Fiso = Piso*Ap;

Padiabatic= P1*((lo-l1)./(lo-l1-Ap*(z-z1)/Aa)).^gamma;
Fadiabatic = Padiabatic*Ap;









figure
subplot(4,1,1);
plot(z,Piso,z,Padiabatic);
ylabel('Pressure (psi)');
xlabel('Rod position (in)');
legend('isothermal','adiabatic','Location','NorthWest');

subplot(4,1,2);
plot(z,Fiso,z,Fadiabatic);
ylabel('Force (lbs)');
xlabel('Rod position (in)');
legend('isothermal','adiabatic','Location','NorthWest');

subplot(4,1,3);
plot(z,Fiso,z,Fadiabatic);
ylabel('Force (lbs)');
xlabel('Rod position (in)');
legend('isothermal','adiabatic','Location','NorthWest');

subplot(4,1,4);
plot(z(2:end),(4.45*39.4)*diff(Fiso)./diff(z),z(2:end),(4.45*39.4)*diff(Fadiabatic)./diff(z));
ylabel('k (N/m)');
xlabel('Rod position (in)');
legend('isothermal','adiabatic','Location','NorthWest');




figure
subplot(2,1,1);
plot(z,l);
ylabel('Accumulator Piston Position (in)');
xlabel('Rod position (in)');

subplot(2,1,2);
plot(z,CR);
ylabel('Compression Ratio');
xlabel('Rod position (in)');