clear all;
clc;

%ADIABATIC ANALYSIS

%Nomenclature

% Ve(j,i)=Maximum Hot end volume
% Vhe=Hot end heat exchanger volume
% Vhd=hot end dead volume
% Vc(j,i)=Cold end volume
% Vc(j,i)e=cold end heat exchanger volume
% Vc(j,i)d=cold end dead volume
% Vd=cold end maximum displacer volume
% Vp=cold end maximum piston volume
% Vr=regenerator volume

% The=hot end heat exchanger temperature
% Tce=cold end heat exchanger temperature
% Tr=regenerator temperature
% Te=expansion space temperature
% Tc=compression space temperature

% me=expansion space mass
% meh=hot end heaat exchanger mass
% mr=regenerator mass
% mc=compression space mass
% mec=cold end heat exchanger mass

% p=instantaneous pressure
% r=Co/Cv for the working gas
% R=gas constant;
% 
% f=crank angle
% AL=phase lag

%Define constant values

R=2.0679*1000;      %Set the value of gas constant for your working fluid
r=1.666666672;      %Set the value of ratio of specific heats for your working fluid
Cp=5.1926*1000;     %Set the value of specific heat of working fluid at constant pressure. Enter the value of specific heat at the heater temperature
Cv=3.1156*1000;     %Set the value of specific heat of working fluid at constant volume. Enter the value of specific heat at the heater temperature


m=0.1;              %This crank angle resolution is optimum when solving explicit relation for pressure
k=10;               %It is known that the convergence should take arounf 5 iterations. However we set the array size of the variables so that it can store upto 10 iteration data. It is not advisable to keep increasing the
                    %array size with each iteration because this will slow
                    %the computational speed.
 
AL=pi/180*90;       %Set the desired/design phase angle between displacer and power piston

flag=0;             %flag will be set to 1 when the error criterion is satiisfied

%Set engine parameters
Vhe=input('Enter hot end heat exchanger volume in cm3: ');
Vce=input('Enter cold end heat exchanger volume in cm3: ');
Vr=input('Enter regenerator section volume in cm3: ');

%Set initial conditions
Th=input('Enter hot end heat exchanger temperature in K: ');
Tk=input('Enter cold end heat exchanger temperature in K: ');
Tr=(Th-Tk)/log(Th/Tk);

n=(360/m)+1;

%Initialize arrays to store variables

p=zeros(k,n);
p(1,1)=input('Enter the initial charging pressure in MPa: ');
Ve=zeros(k,n);
Vc=zeros(k,n);
v=zeros(k,n);
Te=zeros(k,n);
Tc=zeros(k,n);
               

W1=zeros(k,n);

%Volume arrays
Vem=input('Enter expansion end working volume in cm3: ');
Vp=input('Enter compression end piston working volume in cm3: ');
Vcm=Vem;
%Vp=Vp-(Vem/2);
Vd=Vem+Vp;
Ved=30.52;
Vcd=28.68;

Ve(1,1)=(Vd/2*(1-cos(m*pi/180)));       %Start with displacer at topmost position
Vc(1,1)=(Vd/2*(1+cos(m*pi/180))+Vp/2*(1-cos((m*pi/180)-(AL))));
Vpiston=zeros(k,n);
Vcompression=zeros(k,n);
v=zeros(k,n);
v(1,1)=Vc(1,1)+Ve(1,1)+Vr+Vhe+Vce;
v(2,1)=v(1,1);


%Mass arrays
Dmc=zeros(k,n);

%The following variables are nothing but the values of instantaneous mass flowing across
%control volumes
gAck=zeros(k,n);
gAkr=zeros(k,n);
gArh=zeros(k,n);
gAhe=zeros(k,n);
Dmce=zeros(k,n);
Dmr=zeros(k,n);
Dmhe=zeros(k,n);

%These are the instantaneous mass contents of each control volume
mce=zeros(k,n);
mhe=zeros(k,n);
mr=zeros(k,n);
me=zeros(k,n);
mc=zeros(k,n);
meact=zeros(k,n);

%Work and Heat arrays
DQk=zeros(k,n);
Qk=zeros(k,n);
DQr=zeros(k,n);
Qr=zeros(k,n);
DQh=zeros(k,n);
Qh=zeros(k,n);
DW=zeros(k,n);

%Other arrays
f=zeros(k,n);
fd=zeros(k,n);
dexp=zeros(k,n);
dcomp=zeros(k,n);
num=zeros(k,n);
den=zeros(k,n);
Dp=zeros(k,n);
carnot_efficiency=zeros(k,n);
instant_efficiency=zeros(k,n);

Tck=Tc(1,1);
The=Th;

term=r*((Vce/Tck)+(Vr/Tr)+(Vhe/The));


M=(p(1,1)*(Ve(1,1)+Ved)/R/Th)+(p(1,1)*(Vc(1,1)+Vcd)/R/Tk)+(p(1,1)*Vr/R/Tr)+(p(1,1)*Vhe/R/Th)+(p(1,1)*Vce/R/Tk);     %constant mass of gas in engine
mr(1,1)=(p(1,1)*Vr/R/Tr);
mhe(1,1)=(p(1,1)*Vhe/R/Th);
mce(1,1)=(p(1,1)*Vce/R/Tk);
mc(1,1)=(p(1,1)*(Vc(1,1)+Vcd)/R/Tk);
% me(1,1)=M-(mr(1,1)+mhe(1,1)+mce(1,1)+mc(1,1));
Te(1,1)=Th;
Tc(1,1)=Tk;

j=0;

while flag<1
    j=j+1;
    Ve(j,1)=(Vd/2*(1-cos(m*pi/180)));      									 %Start with displacer at topmost position

%Initialize volumes, masses and temperature values 

    Vc(j,1)=(Vd/2*(1+cos(m*pi/180))+Vp/2*(1-cos((m* pi/180)-(AL))));
    Vpiston(j,1)=Vp/2*(1-cos((m* pi/180)-(AL)));
    Vcompression(j,1)=Vd/2*(1+cos(m*pi/180));
    v(j,1)=Vc(j,1)+Ve(j,1)+Vce+Vhe+Vr+Ved+Vcd;
    mr(j,1)=(p(j,1)*Vr/R/Tr);
    mhe(j,1)=(p(j,1)*Vhe/R/Th);
    mce(j,1)=(p(j,1)*Vce/R/Tk);
    mc(j,1)=(p(j,1)*(Vc(j,1)+Vcd)/R/Tc(j,1));
    me(j,1)=M-(mr(j,1)+mc(j,1)+mhe(j,1)+mce(j,1));
    Te(j,1)=p(j,1)*(Ve(j,1)+Ved)/R/me(j,1);
    Tck=Tc(j,1);
    The=Th;

for i=2:n
    fd(j,i)=i*m;
    f(j,i)=i*m*pi/180;
    Ve(j,i)=(Vem/2*(1-cos(f(j,i))));
    dexp(j,i-1)=Ve(j,i)-Ve(j,i-1);
    Vc(j,i)=(Vd/2*(1+cos(f(j,i)))+Vp/2*(1-cos(f(j,i)-AL)));
    dcomp(j,i-1)=Vc(j,i)-Vc(j,i-1);
    Vcompression(j,i)=Vd/2*(1+cos(f(j,i)));
    Vpiston(j,i)=Vp/2*(1-cos(f(j,i)-AL));
    v(j,i)=Vc(j,i)+Ve(j,i)+Vce+Vhe+Vr+Vcd+Ved;
    num(j,i)=-r*p(j,i-1)*(((dcomp(j,i-1))/Tck)+(((dexp(j,i-1))/The)));
    term=r*((Vce/Tck)+(Vr/Tr)+(Vhe/The));
    den(j,i)=((Vc(j,i)+Vcd)/Tck)+term+((Ve(j,i)+Ved)/The);
    Dp(j,i)=num(j,i)/den(j,i);
 
    p(j,i)=p(j,i-1)+Dp(j,i);
    
    Dmc(j,i)=(((p(j,i)*dcomp(j,i-1))+(Vc(j,i)+Vcd)*Dp(j,i)/r))/(R*Tck);
    mc(j,i)=mc(j,i-1)+Dmc(j,i);
    mr(j,i)=(p(j,i)*Vr/R/Tr);
    mhe(j,i)=(p(j,i)*Vhe/R/Th);
    mce(j,i)=(p(j,i)*Vce/R/Tk);
    me(j,i)=M-(mc(j,i)+mr(j,i)+mhe(j,i)+mce(j,i));
    meact(j,i)=p(j,i)*Ve(j,i)/R/Te(j,i-1);
    error(j,1)=abs(me(j,i)-meact(j,i));
    
    Tc(j,i)=p(j,i)*(Vc(j,i)+Vcd)/R/mc(j,i);
    Te(j,i)=p(j,i)*(Ve(j,i)+Ved)/R/me(j,i);
    
    Dmce(j,i)=mce(j,i)*Dp(j,i)/p(j,i);
    Dmr(j,i)=mr(j,i)*Dp(j,i)/p(j,i);
    Dmhe(j,i)=mhe(j,i)*Dp(j,i)/p(j,i);
    
    gAck(j,i)=-Dmc(j,i);
    gAkr(j,i)=gAck(j,i)-Dmce(j,i);
    gArh(j,i)=gAkr(j,i)-Dmr(j,i);
    gAhe(j,i)=gArh(j,i)-Dmhe(j,i);
    
    if gAck(j,i)>0
        Tck=Tc(j,i);
        Tkr=Tk;
    else
        Tck=Tk;
        Tkr=Tk;
    end
    if gAhe(j,i)>0
        The=Th;
        Trh=Th;
    else
        The=Te(j,i);
        Trh=Th;
    end
    DW(j,i)=(p(j,i)+p(j,i-1))/2*(dcomp(j,i-1)+dexp(j,i-1));
    W1(j,i)=W1(j,i-1)+DW(j,i);
    DQk(j,i)=((Vce)*Dp(j,i)*Cv/R)-Cp*((Tck*gAck(j,i))-(Tkr*gAkr(j,i)));
    Qk(j,i)=Qk(j,i-1)+DQk(j,i);
    DQr(j,i)=(Vr*Dp(j,i)*Cv/R)-Cp*((Tkr*gAkr(j,i))-(Trh*gArh(j,i)));
    Qr(j,i)=Qr(j,i-1)+DQr(j,i);
    Qh(j,i)=W1(k,n)-Qk(k,n);
    carnot_efficiency(j,i)=1-(Tc(j,i)/Te(j,i));
end
    p(j+1,1)=p(j,n);
    Tc(j+1,1)=Tc(j,n);
    if j>2
        error(j,1:n)=p(j,1:n)-p(j-1,1:n);
        
        if max(abs(error(j,1:n)))<0.001
            flag=1;
        end
    end
end

    
% fink=50*area/10000/15/M/4/5.19/1000;
% disp('Finkelstein number for the engine = ');
% disp(fink);
% clearance(1,temp1)=h;
% Power(1,temp1)=W1(k,n);
% efficiency(1,temp1)=(W1(k,n)/Qh(k,n))*100;
% temp1=temp1+1;





% figure
% plot(fd(k,1:n),Ve(k,1:n),'-m','linewidth',2);
% hold on;
% plot(fd(k,1:n),Vc(k,1:n),'-r','linewidth',2);
% hold on;
% plot(fd(k,1:n),Vpiston(k,1:n),'-k','linewidth',2);
% hold on;
% plot(fd(k,1:n),v(k,1:n),'-b','linewidth',2);
% hold on;
% legend('expansion working volume','compression working volume','power piston working volume','total engine volume');
% xlabel('Crank angle in degree');
% ylabel('Volume in cm3');
% %  
% % 
% figure;
% plot(fd(k,1:n),me(k,1:n),'-m','linewidth',2);
% hold on;
% plot(fd(k,1:n),mhe(k,1:n),'-r','linewidth',2);
% hold on;
% plot(fd(k,1:n),mr(k,1:n),'-k','linewidth',2);
% hold on;
% plot(fd(k,1:n),mce(k,1:n),'-b','linewidth',2);
% hold on;
% plot(fd(k,1:n),mc(k,1:n),'-c','linewidth',2);
% hold on;
% plot(fd(k,1:n),mc(k,1:n)+me(k,1:n)+mhe(k,1:n)+mce(k,1:n)+mr(k,1:n),'-g');
% legend('hot working volume','hot heat exchanger','regenerator','cold heat exchanger','cold working volume','Total gas mass in engine');
% xlabel('Crank angle in degrees');
% ylabel('Mass in kg');
% 
% figure;
% plot(fd(k,1:n),p(k,1:n),'linewidth',2);
% xlabel('Crank angle in degrees');
% ylabel('Pressure in MPa');
% 
% 
% figure;
% plot(v(k,1:n),p(k,1:n),'linewidth',2);
% xlabel('Volume in cm3');
% ylabel('Pressure in MPa');
% 
% 
% 
% figure;
% plot(fd(k,1:n),gAck(k,1:n),'-r','linewidth',2);
% hold on;
% plot(fd(k,1:n),gAkr(k,1:n),'-b','linewidth',2);
% hold on;
% plot(fd(k,1:n),gArh(k,1:n),'-c','linewidth',2);
% hold on;
% plot(fd(k,1:n),gAhe(k,1:n),'-g','linewidth',2);
% xlabel('Crank angle in degrees');
% ylabel('mass flow in kg');
% legend('Mass flow rate between cold end working volume and heat exchanger','Mass flow rate between cold end heat exchanger and regenerator','Mass flow rate between regenerator and hot end heat exchanger','Mass flow rate between hot end heat exchanger and working volume');
% 
% % 
disp('Compression ratio for the engine = ');
disp(max(v(j,1:n))/min(v(j,1:n)));
Qh(j,n)=W1(j,n)-Qk(j,n);
disp('Heat given per cycle in J = ');
disp(Qh(j,n));
disp('Heat rejected per cycle in J = ');
disp(Qk(j,n));
disp('Regenerator net heat per cycle in J = ');
disp(Qr(j,n));
temp=(W1(j,n)+Qh(j,n)+Qk(j,n))/2;
disp('Work output per cycle in J = ');
disp((W1(j,n)));
disp('Efficiency of the engine in percentage = ');
disp((W1(j,n)/Qh(j,n))*100);
disp('Power output for a 41.72 Hz operating frequency in W = ');
disp(W1(j,n)*16.67);
disp('Heat given per cycle in W = ');
disp(Qh(j,n)*16.67);
disp('Heat rejected per cycle in W = ');
disp(Qk(j,n)*16.67);
% % 
% figure;
% plot(fd(k,1:n),Te(k,1:n),'-r','linewidth',2);
% hold on;
% plot(fd(k,1:n),Tc(k,1:n),'-c','linewidth',2);
% xlabel('Crank angle in degrees');
% ylabel('Temperature in K');
% legend('Expansion space temperature','Compression space temperature');
% 
disp('Mean operating pressure for the engine= ');
disp(mean(p(j,1:n)));

plot(fd(j,2:n),carnot_efficiency(j,2:n)*100);

% 
% figure;
% plot(fd(k,1:n),DW(k,1:n),'-k','linewidth',2);
% xlabel('Crank angle in degrees');
% ylabel('Instantaneous work in J');
% figure;
% plot(fd(k,1:n),W1(k,1:n),'-g');
% hold on;
% plot(fd(k,1:n),Qh(k,1:n),'-r');
% hold on;
% plot(fd(k,1:n),Qr(k,1:n),'-c');
% hold on;
% plot(fd(k,1:n),Qk(k,1:n),'-b');
% xlabel('Crank angle in degrees');
% ylabel('Heat/Work in J');
% legend('Work output','Heat input','Regenerator heat','Heat rejected');