clear;

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

// 二酸化炭素のファンデルワールス状態方程式
a = 3.65E-1;		// Pa m^6 / mol^2
b = 4.28E-5;		// m^3 / mol
// 理想気体の状態方程式
//a = 0;		// Pa m^6 / mol^2
//b = 0;		// m^3 / mol

// *** ファンデルワールスの状態方程式 **
deff("p = f(v,t)","p = r .* t ./ (v - b) - a ./ v ./ v");

// 体積・温度ベクトルと計算グリッド
v = [5E-5:1E-5:1E-3];	// 体積ベクトル (m^3 / mol)
t = [200:10:500];	// 温度ベクトル (K)
[V,T] = ndgrid(v,t);	// 体積温度空間グリッド

// *** 圧力の計算と2次元グラフのプロット ***
scf(0);
plot(v,f(v,250),'-r');
plot(v,f(v,275),'-g');
plot(v,f(v,300),'-b');
plot(v,f(v,325),'-m');
plot(v,f(v,350),'-c');
// グラフの装飾
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]);
legend(["T = 250 (K)","T = 275 (K)","T = 300 (K)","T = 325 (K)","T = 350 (K)"],1);

// *** 圧力の計算と3次元グラフのプロット ***
scf(1);
P = f(V,T);
mesh(V,T,P);
// グラフの装飾
a = gca();
a.data_bounds = [0,200,0;1E-3,500,2E7];
xlabel("Volume (m^3/mol)");
ylabel("Temperature (K)");
zlabel("Pressure (Pa)");