%The following data is extracted from Table A-4 of Cengel and Turner,
%"Fundamentals of Thermal-Fluid Sciences", McGraw-Hill 2001.
%I want to make a phase diagram in the PV plane for water with it.
%Ted Golfinopoulos, 5.16/2008


%The data is organized as temperature (degrees Celsius), pressure (kPa), saturated liquid
%specific volume (m^3/kg), saturated vapor specific volume (m^3/kg).
T=[0.01 5:5:320 330:10:370 374.14]';

vg= [ 206.14
147.12
       106.38
       77.93
       57.79
       43.36
       32.89
       25.22
      19.52
       15.26
       12.03
       9.568
       7.671
       6.197
       5.042
          4.131
    3.407
    2.828
    2.361
    1.982
    1.6729
    1.4194
    1.2102
    1.0366
    0.8919
    0.7706
    0.6685
    0.5822
    0.5089
    0.4463
    0.3928
    0.3468
    0.3071
    0.2727
    0.2428
    0.2168
    0.19405
    0.17409
    0.15654
    0.14105
    0.12436
    0.11521
    0.10441
    0.09479
    0.08619
    0.07849
    0.07158
    0.06537
    0.05976
    0.05471
    0.05013
    0.04598
    0.04221
    0.03877
    0.03564
    0.03279
    0.03017
    0.02777
    0.02557
    0.02354
    0.02167
    0.019948
    0.018350
    0.016867
    0.015488
    0.012996
    0.010797
    0.008813
    0.006945
    0.004925
    0.003155];
    
    
    
    
    
    
    
    
    
    
    
    
    
vf=[    0.001
    0.001
    0.001
    0.001001
    0.001002
    0.001003
    0.001004
    0.001006
0.001008
    0.001010
    0.001012
    0.001015
    0.001017
    0.001020
    0.001023
0.001026
0.001029
0.001033
0.001036
0.001040
0.001044
0.001048
0.001052
0.001056
0.001060
0.001065
0.001070
0.001075
0.001080
0.001085
0.001091
0.001096
0.001102
0.001108
0.001114
0.001121
0.001127
0.001134
0.001141
0.001149
0.001157
0.001164
0.001173
0.001181
0.001190
0.001199
0.001209
0.001219
0.001229
0.001240
0.001251
0.001263
0.001276
0.001289
0.001302
0.001317
0.001332
0.001348
0.001366
0.001384
0.001404
0.001425
0.001447
0.001472
0.001499
0.001561
0.001638
0.001740
0.001893
0.002213
0.003155];

P=[0.6113
0.8721
1.2276
1.7051
2.339
3.169
4.246    
5.628    
7.384
9.593
12.349
15.758
19.940 
25.03
31.19    
38.58   
47.39
57.83
70.14
84.55
101.33
120.82
143.27
169.03
198.53
232.1
270.1
313.0
361.3
415.4
475.8
543.1
617.8
700.5
791.7
892.0
1002.1
1122.7
1254.4
1397.8
1553.9
1723.0
1906.2
2104
2318
2548
2795
3060
3344
3648
3973
4319
4688
5081
5499
5942
6412
6909
7436
7993
8581
9202
9856
10547
11274
12845
14586
16513
18651
21030
22090];

%Superheated water isotherms:
Ts=[350; 
    ];

Piso=[500, 600, 800, 1000, 1200, 1400 1600 1800 2000 2500 3000 3500 4000 4500 5000 6000 7000 8000 9000 10000 12500 15000 20000];% 30000 50000];


%For T=200: [21.825,  4.356,, 2.172, 1.0803, 0.7163, 0.5342, 0.4249, 0.3520, 0.2608,
%0.2060, 0.1693, 14302;...
v350=[0.5701 0.4742 0.3544 0.2825 0.2345 0.2003 0.17456 0.15457 0.13857 0.10976 0.09053 0.07678 0.06645 0.05840 0.05194 0.04223 0.03524 0.02995 0.02580 0.02242 0.016123 0.011470 (0.0015684+0.0018226)*0.5];% (0.0014920+0.0016265)*0.5 (0.0014032+0.0014838)*0.5];
v400=[0.6173 0.5137 0.3843 0.3066 0.2548 0.2178 0.19005 0.16847 0.15120 0.12010 0.09936 0.08453 0.07341 0.06475 0.05781 0.04739 0.03993 0.03432 0.02993 0.02641 0.02000 0.015649 0.009942];% 0.002790 0.0017309 ];



Pcrit=[Piso 25000 30000 35000 40000 50000 60000];
vcrit=[(v350+v400)/2  0.0019731 0.0017892 0.0017003 0.0016407 0.0015594 0.0015028];

%Psuper=linspace(min(Piso), max(Piso), 50);
vsuper=v350(1:length(v350)-1);%logspace(min(vcrit), max(vcrit), 50);
%van der Waals isotherm parameters for water from Wikipedia,
%http://en.wikipedia.org/wiki/Van_der_Waals_constants_%28data_page%29
kg2mol=1000/18; mol2kg=1/kg2mol;
a=5.536/(mol2kg^2)/10; b=0.03049/1000; kB=1.3806504E-23; Na=6.02E23; Ts=450+273;

p=Na*kB*Ts./(vsuper*mol2kg-b)-a.*(vsuper).^-2;
p=p/1000; %Convert to kpa.
%Psuper=fsolve(@(p) ( p+a*(vsuper.^(-2))).*(vsuper-b)-kB*Ts, Pcrit);

xxf=logspace(log10(min(vf)), log10(max(vf)), 100)';
xxg=logspace(log10(min(vg)), log10(max(vg)), 100)';
cs=spline(vf(3:length(vf)), [P(3); P(3:length(P)); P(length(P))]);
yyf=ppval(cs,xxf);
cs=spline(vg, [P(1); P; P(length(P))]);
yyg=ppval(cs,xxg);
figure(1)
clf(1)
hold on
%semilogx(xxf, yyf, 'k', xxg, yyg, 'k')
semilogx(vf(3:length(vf)), P(3:length(P)), 'k', vg, P, 'k')
semilogx(vcrit, Pcrit, 'r.')