clear;

// *** 入力パラメータ ***
// 物理定数
r = 8.31		// 気体定数 (J/K/mol)

// 二酸化炭素のファンデルワールス状態方程式
a = 3.65E-1;		// Pa m^6 / mol^2
b = 4.28E-5;		// m^3 / mol
// 温度 (K)
t = 250;

// *** ファンデルワールスの状態方程式 **
function p = pvdw(v)
  p = r .* t ./ (v - b) - a ./ v ./ v
endfunction

// *** 解くべき方程式 ***
function e = func(pk)
  Vx = roots([pk, -1 * (b * pk + r * t), a, - a * b]);
  vg = max(real(Vx));
  vl = min(real(Vx));
  s1 = integrate('pvdw(v)','v',vl,vg);
  s2 = (vg - vl) * pk;
  e = s1 - s2;
endfunction

// *** 非線形方程式の数値解 ***
pk0 = %eps;
pk = fsolve(pk0,func);

Vx = roots([pk, -1 * (b * pk + r * t), a, - a * b]);
vg = max(real(Vx));
vl = min(real(Vx));

// *** 圧力の計算とプロット ***
// モル体積ベクトル
Vl = linspace(b,vl,1000) + %eps;	// (m^3 / mol)
Vm = linspace(vl,vg,1000);	// (m^3 / mol)
Vg = linspace(vg,1e-3,1000);	// (m^3 / mol)
// 温度
plot(Vl,pvdw(Vl),'-r');
plot(Vm,pk,'-r');
plot(Vm,pvdw(Vm),'--r');
plot(Vg,pvdw(Vg),'-r');
// *** グラフの装飾 ***
xlabel("Molar volume (m^3/mol)");
ylabel("Pressure (Pa)");
plot([b,b],[-1E8,2E7],'--k');
plot([0,0.001],[0,0],'--k');
zoom_rect([0,-2E6,1E-3,2E7]);