clear all

% Plot and label an accurate Carnot cycle.
% 
% Executes properly on MATLAB version 9.4.0.813654 (R2018a).
% Executes mostly properly on OCTAVE version 4.2.1.
% Joshua C. Vaughan 1/4/2021


T1 = 250; % K
T2 = 300; % K
n = 1;    % mol
R = 8.31; % J/(K mol)
gam = 5/3; % unitless

V = linspace(0.001, 0.5, 2^10);  % m^3

% isotherms here
Piso1 = n*R*T1./V;
Piso2 = n*R*T2./V;

figure(1); clf
plot(V, Piso1, V, Piso2)
ylim([0 1e5])
xlim([0 max(V)])

% find point where adiabat intersects
V_val1 = 0.31;
V_val2 = 0.32;
[jnk, ind1] = min( abs( V - V_val1 ) );
[jnk, ind2] = min( abs( V - V_val2 ) );

% adiabats here
Padi1 = Piso1(ind1)*(V(ind1)^gam)./(V.^gam);
Padi2 = Piso2(ind2)*(V(ind2)^gam)./(V.^gam);

% Get ranges to plot. Need four ranges.

% Adiabatic 1 range
[jnk, adi1_ind1] = min( abs( Padi1-Piso1 ) );
[jnk, adi1_ind2] = min( abs( Padi1-Piso2 ) );
adi1_inds = adi1_ind2:1:adi1_ind1;

% Adiabatic 2 range
[jnk, adi2_ind1] = min( abs( Padi2-Piso1 ) );
[jnk, adi2_ind2] = min( abs( Padi2-Piso2 ) );
adi2_inds = adi2_ind2:1:adi2_ind1;

% Isothermal 1 range
iso1_inds = adi1_ind1:1:adi2_ind1;

% Isothermal 2 range
iso2_inds = adi1_ind2:1:adi2_ind2;

% Plot Carnot cycle
ax_lw = 2;
pl_lw = 3;
ax_fs = 24;
la_fs = 24;
ti_fs = 24;
mk_sz = 10;
fontname = 'Arial';

figure(1); clf
hold on
% plot carnot cycle
plot(V(iso2_inds), Piso2(iso2_inds)/1e4, 'r', 'linewidth', pl_lw)
plot(V(adi2_inds), Padi2(adi2_inds)/1e4, 'r--', 'linewidth', pl_lw)
plot(V(iso1_inds), Piso1(iso1_inds)/1e4, 'b', 'linewidth', pl_lw)
plot(V(adi1_inds), Padi1(adi1_inds)/1e4, 'b--', 'linewidth', pl_lw)


% add vertices
plot(V(adi1_inds([1 end])), Padi1(adi1_inds([1 end]))/1e4, 'ko','markersize',mk_sz,'markerfacecolor','k')
plot(V(adi2_inds([1 end])), Padi2(adi2_inds([1 end]))/1e4, 'ko','markersize',mk_sz,'markerfacecolor','k')

% add text labels
vt_fs = 32;
text(V(adi1_inds(1))*0.94,Padi1(adi1_inds(1))/1e4,'A','fontsize',vt_fs,'horizontalalignment','center','fontweight','bold','color','k')
text(V(adi1_inds(end))*0.94,Padi1(adi1_inds(end))/1e4,'D','fontsize',vt_fs,'horizontalalignment','center','fontweight','bold','color','k')
text(V(adi2_inds(1))*1.06,Padi2(adi2_inds(1))/1e4,'B','fontsize',vt_fs,'horizontalalignment','center','fontweight','bold','color','k')
text(V(adi2_inds(end))*1.04,Padi2(adi2_inds(end))/1e4,'C','fontsize',vt_fs,'horizontalalignment','center','fontweight','bold','color','k')

% add legend
%h = legend('adiabatic, warming','2','3','4');
h = legend('isothermal, {\itT}_{hot}', 'adiabatic, cooling','isothermal, {\itT}_{cool}','adiabatic, warming');
set(h,'fontsize',16);
set(h,'linewidth',1.5);
% 

hold off

xlabel('Volume [m^3]','fontsize',la_fs)
ylabel('Pressure [\times10^4 Pa]','fontsize',la_fs)
%title('Carnot Cycle','fontsize',ti_fs)
set(gca,'fontsize',la_fs)
set(gca,'linewidth',ax_lw)
set(gca,'tickdir','out')
set(gcf,'color',[1 1 1])
set(gca,'fontname',fontname)
xlim([0.21 0.45])
ylim([0.4e4 1.15e4]/1e4)
set(gca,'LooseInset',get(gca,'TightInset'));

