% vdwvol.m                                       by:  Edward T Peltzer, MBARI
%                                                  revised:  01 Dec 98
%
% CALCULATE THE VOLUME OF A GAS AT A GIVEN T & P
%     using the Van der Waals Equation of State.
%
%   This equation works best at P = 1 atm or less.
%
% INPUT:       	Gas selector:	 1 == methane
%				 2 == carbon dioxide
%				 3 == nitrogen
%				 4 == oxygen
%				 5 == argon
%				 6 == ethane
%				 7 == propane
%				 8 == n-butane
%				 9 == iso-butane
%				10 == n-pentane
%				11 == iso-pentane
%				12 == n-hexane
%
%		Temperature (T) in degrees C.
%		Pressure in kPa  ==  P(mbar) / 10.
%                                    P(dbar) * 10.
%
% OUTPUT:	Volume of 1 mole of gas [vpm] in liters.
%
%			vpm = vdwvol(n,T,P).

function [vpm]=vdwvol(n,T,P)


% DEFINE CONSTANTS FOR EQUATION OF STATE

    A = [2.253 3.592 1.390 1.360 1.345 5.489 8.664 14.47 12.87 19.01 18.05 24.39];
    B = [4.278 4.267 3.913 3.183 3.219 6.380 8.445 12.26 11.42 14.60 14.17 17.35];
    Bc = B ./ 100;
    R = 0.082057;


%  CONVERT INPUT DATA TO APPROPRIATE UNITS

    Tc = T + 273.15;
    Pc = P ./ 101.325;


%  CHECK TO SEE IF PRESSURE 'IN-RANGE'

    if Pc > 5
	disp(' ')
	disp(' ')
	disp('Warning: input pressure out-of-range,')
	disp('           answer may be incorrect!')
    end


% CALCULATE COEFFICIENTS FOR CUBIC EQUATION

    C2 = 0 - ((Pc .* Bc(n)) + (R .* Tc));
    C4 = 0 - A(n) .* Bc(n);


% SOLVE CUBIC

    rts = roots([Pc C2 A(n) C4]);


% FIND REAL ROOT == VPM

    L = length(rts);
    nreal = 0;

    for I = 1:L
	if rts(I) == real(rts(I))
	    nreal = nreal + 1;
	    idx = I;
	end
    end

    if nreal == 1
	vpm = rts(idx);
    else
	vdif = abs(rts - idgvol(T,P));
	[MinDif,idx] = min(vdif);
	vpm = rts(idx);
    end