close all
clear all

%gamma = 7/5;  %Adiabatic compression of a Diatomic ideal gas
gamma = 1;  %Isothermal compression of a Diatomic ideal gas

T0 = 15+273;  %Degrees Kelvin

Amp = 40;

z = -Amp:Amp;

D = 5;  %Bore Diameter
d = 2; %Rod Diameter

Au = 0.25*pi*D^2;
Al = 0.25*pi*(D^2-d^2);

k = 1600;  %pounds per full stroke (40 inches).

Pu0 = 70.0;
Pl0 = (k+Pu0*Au)/Al;


Vu0 = 2432.0;
Vl0 = 3384.0;


Vu = Vu0+Au*z;
Vl = Vl0-Al*z;


Cu = Pu0*Vu0^gamma;
Cl = Pl0*Vl0^gamma;


Pu = Cu./(Vu.^gamma);
Pl = Cl./(Vl.^gamma);

PuIso = Pu0*Vu0./Vu;
PlIso = Pl0*Vl0./Vl;

Fu = Pu*Au;
Fl = Pl*Al;
F = Fl-Fu;

FuIso = PuIso*Au;
FlIso = PlIso*Al;
FIso = FlIso-FuIso


TCu = T0*Vu(ceil(length(z)/2))^(gamma-1);
TCl = T0*Vl(ceil(length(z)/2))^(gamma-1);

Tu = TCu./(Vu.^(gamma-1));
Tl = TCl./(Vl.^(gamma-1));


figure('Position', [55 100 673 844],'PaperOrientation','portrait','PaperPosition',[0 0 8.5 11])

subplot(4,1,1)
plot(z,Vu,z,Vl);

ylabel('in^3');
xlabel('z (inches)')
legend('Upper Chamber Volume','Lower Chamber Volume');
%title('Adiabatic Compression of Diatomic Gas:  k = 1600 lbs/m');
title('Isothermal Compression:  k = 1600 lbs/m');

subplot(4,1,2)
plot(z,Pu,z,Pl);
legend('Upper Chamber Pressure','Lower Chamber Pressure','Location','NorthWest');
ylabel('psi');
xlabel('z (inches)')

subplot(4,1,3)
plot(z,Fu,z,Fl,z,F);
legend('Upper Chamber Force','Lower Chamber Force','Net Force','Location','NorthWest');
ylabel('lbs');
xlabel('z (inches)')

subplot(4,1,4)
plot(z,Tu-273,z,Tl-273);
legend('Upper Chamber Temperature','Lower Chamber Temperature');
ylabel('Degrees C');
xlabel('z (inches)')

%print('-dpdf','SpringProfiles-Adiabatic-k1600');
%print('-dpdf','SpringProfiles-Isothermal-k1600');
