Implement code for multi range of temperature

Hi

I need to execute and plot the below code for range of temperatures (T) instead of one value (T=333)

exampel range of T value between 310 and 350 with 5 incremental steps

clc; clear all; close all;
 
T =333;
Tc = 369.8;
Tr = T/Tc;
Vr = linspace(0.5,5,500);
funcPr = @(Vr) 8/3 * Tr./(Vr - 1/3) - 3./(Vr.^2);
Pr = funcPr(Vr);
 
if Tr < 1
    Pr_l = 1.0;                                         
    vdW_Pr_l = [1 -1/3*(1+8*Tr/Pr_l) 3/Pr_l -1/Pr_l];
    v=sort(roots(vdW_Pr_l));
    A1 = (v(2)-v(1))*Pr_l - integral(funcPr,v(1),v(2));
    A2 = integral(funcPr,v(2),v(3)) - (v(3)-v(2))*Pr_l;
    AreaDiff = abs (A1-A2);
  while AreaDiff > 0.0001
    vdW_Pr_l = [1 -1/3*(1+8*Tr/Pr_l) 3/Pr_l -1/Pr_l];
    v=sort(roots(vdW_Pr_l));
    funcPr = @(Vr) 8/3 * Tr./(Vr - 1/3) - 3./(Vr.^2);
    A1 = (v(2)-v(1))*Pr_l - integral(funcPr,v(1),v(2));
    A2 = integral(funcPr,v(2),v(3)) - (v(3)-v(2))*Pr_l;
    AreaDiff = abs (A1-A2);
    Pr_l = Pr_l - 0.00001;
  end
end
 
% Plot
plot(Vr,Pr); hold all
xlim([0.25 4])
ylim([-0.5 2])
xlabel('Vr')
ylabel('Pr')
plot([0.25 4],[Pr_l Pr_l],'k--')
plot(v(1),Pr_l,'bo',v(2),Pr_l,'bo',v(3),Pr_l,'bo')

Answer 1

[Sajeer Modavan]

I don't know what you are looking exactly. I did something here, please let me know is it correct or you need further help

clc; clear all; close all;
 
T   = 310:5:350;
Tc  = 369.8;
Tr1 = T./Tc;
Vr  = linspace(0.5,5,500);

figure
for ii = 1:length(Tr1)
    Tr = Tr1(ii);
funcPr = @(Vr) 8/3 * Tr./(Vr - 1/3) - 3./(Vr.^2);
Pr = funcPr(Vr);

if Tr < 1
    Pr_l = 1.0;                                         
    vdW_Pr_l = [1 -1/3*(1+8*Tr/Pr_l) 3/Pr_l -1/Pr_l];
    v=sort(roots(vdW_Pr_l));
    A1 = (v(2)-v(1))*Pr_l - integral(funcPr,v(1),v(2));
    A2 = integral(funcPr,v(2),v(3)) - (v(3)-v(2))*Pr_l;

    AreaDiff = abs (A1-A2);
  while AreaDiff > 0.0001
    vdW_Pr_l = [1 -1/3*(1+8*Tr/Pr_l) 3/Pr_l -1/Pr_l];
    v=sort(roots(vdW_Pr_l));
    funcPr = @(Vr) 8/3 * Tr./(Vr - 1/3) - 3./(Vr.^2);
    A1 = (v(2)-v(1))*Pr_l - integral(funcPr,v(1),v(2));
    A2 = integral(funcPr,v(2),v(3)) - (v(3)-v(2))*Pr_l;

    AreaDiff = abs (A1-A2);
    Pr_l = Pr_l - 0.00001;
  end
end
 
plot(Vr,Pr); hold all
xlim([0.25 4])
ylim([-0.5 2])
xlabel('Vr')
ylabel('Pr')
plot([0.25 4],[Pr_l Pr_l],'k--')
plot(v(1),Pr_l,'bo',v(2),Pr_l,'bo',v(3),Pr_l,'bo')
end

[Sajeer Modavan]

Do you need any further help, accept the answer and close this issue if you don't need any further help.