This is wengine.c in view mode; [Download] [Up]
#include "wstruct.h"
#include "wconst.c"
/* The engine code. Adapted, optimized, fixed, hassled with, and torqued
by: M.L. (Loki) Demsey
University of Arizona
Geo Sciences Department
for and help from: Dr. Denis Norton
Original code, molded in FORTRAN from personal code and code from Haar,
Gallagher, and Kell by James Johnson
This is the commented-to-death version of the code; it is recommended
that this be used only for reference and not compilation.
The user realizes that this code is not supported and may
be superceded at any moment.
*/
/* **********************************************************************
*
* specs - Input unit, triple point, saturation, and option specs:
*
***** it, id, ip, ih, itripl, isat, iopt, useLVS, epseqn, icrit;
*
* note that the returned value of isat may differ from
* its input value and that icrit need not be specified
* prior to invocation.
*
*
* states - State variables:
*
***** temp, pres, dens(1), dens(2);
*
* note that the first three of these must be specified prior
* to invocation and that, in the case of saturation, vapor
* density is returned in dens(1), liquid in dens(2).
*
*
* props - Thermodynamic, transport, electrostatic, and combined
* property values:
*
***** A, G, S, U, H, Cv, Cp, Speed, alpha, beta, diel, visc,
***** tcond, surten, tdiff, Prndtl, visck, albe,
***** ZBorn, YBorn, QBorn, daldT, XBorn
*
*
* error - LOGICAL argument that indicates success ("FALSE") or
* failure ("TRUE") of the call, the latter value in
* response to out-of-bounds specs or states variables.
*
********************************************************************** */
void h2o92(specs,states,props,error)
int *error,specs[10];
double states[4],props[46];
{
extern void unit(),LVSeqn(),HGKeqn(),load(),LVScon(),HGKcon();
extern int valid(),crtreg();
extern double TdegUS();
double Dens[2],tempy;
int useLVS;
HGKcon();
LVScon();
unit(specs[0],specs[1],specs[2],specs[3],specs[4]);
if (! valid(specs[0],specs[1],specs[2],specs[3],specs[4], &specs[5],specs[6],specs[7],specs[8],states[0],states[1],states[2]))
{ *error = 1;
return;
}
else
{ *error = 0; }
if ( crtreg(specs[5],specs[6],specs[0],&states[0],&states[1],&states[2]))
{ specs[9] = 1;
useLVS = (specs[7] == 1);
}
else
{ specs[9] = useLVS = 0; }
if (useLVS)
{ Dens[0] = states[2];
LVSeqn(&specs[5],&specs[6],specs[4],&states[0],&states[1],Dens,specs[8]);
Dens[0] /= 1000.0;
if (specs[5] == 1)
{ Dens[1] /= 1000.0; }
}
else
{ Dens[0] = states[2] / 1000.0;
HGKeqn(&specs[5],specs[6],specs[4],&states[0],&states[1],Dens,specs[8]);
}
load(1,wpvals.wprops,props);
if (specs[5] == 1)
{ tempy = Dens[0];
Dens[0] = Dens[1];
Dens[1] = tempy;
load(2,wpvals.wpliq,props);
}
states[0] = TdegUS(specs[0],states[0]);
states[1] *= units.fp;
states[2] = Dens[0] / units.fd;
if (specs[5] == 1)
{ states[3] = Dens[1] / units.fd; }
return;
}
/* -----------------------------------------------------------------
*** valid - Returns "TRUE" if unit and equation specifications
* are valid and input state conditions fall within
* the HGK equation's region of validity;
* returns "FALSE" otherwise. */
int valid(it,id,ip,ih,itripl,isat,iopt,useLVS,epseqn,Temp,Pres,Dens)
int it,id,ip,ih,itripl,iopt,useLVS,epseqn,*isat;
double Temp,Pres,Dens;
{ /* note switch from usual return declarations to suit ifs */
extern int valspc(),valTD(),valTP();
extern double TdegK();
double T,D,P,Ttripl,Tcrit,Pcrit;
int val;
/* ensure validity of input specifications */
if (! valspc(it,id,ip,ih,itripl,*isat,iopt,useLVS,epseqn))
return 0;
/* convert to degC, bars, g/cm3 */
T = TdegK(it,Temp) - 273.15;
D = Dens * units.fd;
P = (Pres / units.fp) * 10.0;
Ttripl = tpoint.Ttripl - 273.15;
Tcrit = crits.Tc - 273.15;
Pcrit = crits.Pc * 10.0;
if (*isat == 0)
{ if (iopt == 1)
val = valTD(T,D,isat,epseqn);
else
val = valTP(T,P);
}
else
{ if (iopt == 1)
val = (((T + tolers.FPTOL) >= Ttripl) && ((T - tolers.FPTOL) <= Tcrit));
else
val = (((P + tolers.FPTOL) >= tpoint.Ptripl) && ((P - tolers.FPTOL) <= Pcrit));
}
return val;
}
/* -----------------------------------------------------------------
*** valspc - Returns "TRUE" if it, id, ip, ih, itripl, isat, iopt,
* useLVS, and epseqn values all define valid input;
* returns "FALSE" otherwise. */
int valspc(it,id,ip,ih,itripl,isat,iopt,useLVS,epseqn)
int it,id,ip,ih,itripl,isat,iopt,useLVS,epseqn;
{
return ((1 <= it) && (it <= 4) && (1 <= id) && (id <= 4) && (1 <= ip) && (ip <= 5) && (1 <= ih) && (ih <= 6) && (0 <= itripl) && (itripl <= 1) && (0 <= isat) && (isat <= 1) && (1 <= iopt) && (iopt <= 2) && (0 <= useLVS) && (useLVS <= 1) && (1 <= epseqn) && (epseqn <= 5));
}
/* -----------------------------------------------------------------
*** valTD - Returns "TRUE" if T-D defines liquid or vapor H2O
* within validity limits of the HGK equation of state;
* returns "FALSE" otherwise. */
int valTD(T,D,isat,epseqn)
double T,D;
int *isat;
int epseqn;
{
extern void pcorr(),denLVS(),bb(),denHGK();
extern double Pfind(),Psublm();
int val,istemp;
double Tk,Tcrit,Ttripl,Dlimit,P,Ps,Dv,Dl,PMPa,Dguess,Dsublm,dPdD;
if (((T - tolers.FPTOL) > HGKbnd.Ttop) || ((T + tolers.FPTOL) < HGKbnd.Tbtm) || ((D - tolers.FPTOL) > HGKbnd.Dtop) || ((D + tolers.FPTOL) < HGKbnd.Dbtm))
return 0;
Tcrit = crits.Tc - 273.15;
Ttripl = tpoint.Ttripl - 273.15;
if (((T + tolers.FPTOL) >= Tcrit) || ((T >= liqice.TnIB30) && (D >= tpoint.Dltrip)))
{ Dlimit = liqice.sDIB30 * (T - liqice.TnIB30) + HGKbnd.Dtop;
val = ((D - tolers.FPTOL) <= Dlimit);
}
else
{ if ((D - tolers.FPTOL) <= tpoint.Dltrip)
{ if (T >= Ttripl)
{ val = 1;
Tk = T + 273.15;
if (Tk < coefs1.x[0])
{ RTcurr.rt = aconst.gascon * Tk;
pcorr(0,Tk,&Ps,&Dl,&Dv,epseqn);
}
else
{ istemp = 1;
satur.DH2O = 0.0;
P = Pfind(istemp,&Tk,&satur.DH2O);
denLVS(istemp,&Tk,P);
Dv = satur.Dvap / 1000.0;
Dl = satur.Dliq / 1000.0;
}
if ((D >= Dv) && (D <= Dl))
*isat = 1;
}
else
{ P = Psublm(T);
PMPa = P / 10.0;
Tk = T + 273.15;
Dguess = PMPa / Tk / 0.4;
RTcurr.rt = aconst.gascon * Tk;
bb(Tk);
denHGK(&Dsublm,PMPa,Dguess,Tk,&dPdD);
val = ((D - tolers.FPTOL) <= Dsublm);
}
}
else
{ if (D <= liqice.Dli13)
{ Dlimit = liqice.sDli1 * (T - liqice.Tli13) + liqice.Dli13;
val = ((D + tolers.FPTOL) >= Dlimit);
}
else
{ Dlimit = liqice.sDli37 * (T - liqice.Tli13) + liqice.Dli13;
val = ((D - tolers.FPTOL) <= Dlimit);
}
}
}
return val;
}
/* -----------------------------------------------------------------
*** valTP - Returns "TRUE" if T-P defines liquid or vapor H2O
* within validity limits of the HGK equation of state;
* returns "FALSE" otherwise. */
int valTP(T,P)
double T,P;
{
int val;
double Plimit,Psubl;
extern double Psublm();
if (((T - tolers.FPTOL) > HGKbnd.Ttop) || ((T + tolers.FPTOL) < HGKbnd.Tbtm) || ((P - tolers.FPTOL) > HGKbnd.Ptop) || ((P + tolers.FPTOL) < HGKbnd.Pbtm))
return 0;
else
val = 1;
if (P >= liqice.Pli13)
{ Plimit = liqice.sPli37 * (T - liqice.Tli13) + liqice.Pli13;
val = ((P - tolers.FPTOL) <= Plimit);
}
else
{ if (P >= tpoint.Ptripl)
{ Plimit = liqice.sPli1 * (T - liqice.Tli13) + liqice.Pli13;
val = ((P + tolers.FPTOL) >= Plimit);
}
else
{ Psubl = Psublm(T);
val = ((P - tolers.FPTOL) <= Psubl);
}
}
return val;
}
/* -----------------------------------------------------------------
*** Psublm - Returns Psublimation(T) computed from the
* equation given by Washburn (1924): Monthly
* Weather Rev., v.52, pp.488-490. */
double Psublm(Temp)
double Temp;
{
double T,PmmHg;
T = Temp + 273.1;
PmmHg = pow(10.0, (-2445.5646/T + 8.2312*log10(T) - 0.01677006*T + 0.0000120514*T*T - 6.757169));
return (PmmHg * 0.00133322);
/* converted from mmHg to bars */
}
/* -----------------------------------------------------------------
*** unit - Sets internal parameters according to user-specified
* choice of units. Internal program units are degK(T),
* and gm/cm**3(D); all other properties are computed in
* dimensionless form and dimensioned at output time.
* NOTE: conversion factors for j/g ---> cal/(g,mole)
* (ffh [3 & 4]) are consistent with those given in
* Table 1, Helgeson & Kirkham (1974a) for thermal calories,
* and differ slightly with those given by Haar et al (1984)
* for international calories. */
void unit(it,id,ip,ih,itripl)
int it,id,ip,ih,itripl;
{
extern void tpset();
double fft[4] = {1.0,1.0,0.555555556,0.555555556};
double ffd[4] = {0.001,1.0,0.0180152,0.016018};
double ffvd[4] = {1.0,10.0,0.555086816,0.671968969};
double ffvk[4] = {1.0,10000.0,10000.0,10.76391042};
double ffs[4] = {1.0,100.0,100.0,3.280833};
double ffp[5] = {1.0,10.0,9.86923266716,145.038,10.1971};
double ffh[6] = {1.0,1.0,18.0152,0.23901,4.305816,0.4299226};
double ffst[4] = {1.0,1000.0,55.5086816,2.205061};
double ffcd[4] = {1.0,0.01,0.01,0.3048};
double ffch[6] = {0.001,1.0,1.0,0.23901,0.23901,0.000947244};
units.ft = fft[it-1];
units.fd = ffd[id-1];
units.fvd = ffvd[id-1];
units.fvk = ffvk[id-1];
units.fs = ffs[id-1];
units.fp = ffp[ip-1];
units.fh = ffh[ih-1];
units.fst = ffst[id-1];
units.fc = ffcd[id-1] * ffch[ih-1];
if (itripl == 1)
tpset();
return;
}
/* -----------------------------------------------------------------
*** crtreg - Returns "TRUE" if input state conditions fall within
* the critical region of H2O; otherwise returns "FALSE".
* T, P, D, input in user-specified units, are returned in
* degK, MPa, kg/m3. */
int crtreg(isat,iopt,it,T,P,D)
int isat,iopt,it;
double *T,*P,*D;
{
int llim,ulim,val,isat1;
double Pmin,Pmax,ddummy,Pstest;
extern double Pfind();
*T = TdegK(it,*T);
if (isat == 0)
{ if (iopt == 1)
{ *D *= units.fd * 1000.0;
val = ((*T >= coefs1.x[0]) && (*T <= coefs1.x[2]) && (*D >= coefs1.x[3]) && (*D <= coefs1.x[4]));
}
else
{ *P /= units.fp;
if ((*T < coefs1.x[0]) || (*T > coefs1.x[2]))
val = 0;
else
{ Pmin = coefs1.x[5] + coefs1.x[9] * (*T - coefs1.x[0]);
Pmax = coefs1.x[6] + coefs1.x[10] * (*T - coefs1.x[1]);
llim = (*P >= Pmin);
ulim = (*P <= Pmax);
if (llim && ulim)
val = 1;
else
{ if (llim && (*T <= coefs1.x[1]))
{ isat1 = 1;
ddummy = 0.0;
Pstest = Pfind(isat1,T,&ddummy);
val = (*P <= Pstest);
}
else
val = 0;
}
}
}
}
else
{ if (iopt == 1)
val = (*T >= coefs1.x[0]);
else
{ *P /= units.fp;
val = (*P >= coefs1.x[5]);
}
}
return val;
}
/* -----------------------------------------------------------------
*** HGKeqn - Computes thermodynamic and transport properties of
* of H2O from the equation of state given by
* Haar, Gallagher, & Kell (1984). */
void HGKeqn(isat,iopt,itripl,Temp,Pres,Dens,epseqn)
int *isat,iopt,itripl,epseqn;
double *Temp,*Pres,Dens[2];
{
int i;
extern void bb(),calcv3(),thmHGK(),dimHGK();
RTcurr.rt = aconst.gascon * *Temp;
HGKsat(isat,iopt,itripl,Temp,Pres,Dens,epseqn);
if (*isat == 0)
{ bb(*Temp);
calcv3(iopt,itripl,*Temp,Pres,&Dens[0],epseqn);
thmHGK(Dens[0],*Temp);
dimHGK(*isat,itripl,*Temp,*Pres,Dens[0],epseqn);
}
else
{ for (i=0;i < 23;i++)
wpvals.wpliq[i] = wpvals.wprops[i];
dimHGK(2,itripl,*Temp,*Pres,Dens[1],epseqn);
}
return;
}
/* -----------------------------------------------------------------
*** HGKsat - If isat=1, computes Psat(T) or Tsat(P) (iopt=1,2),
* liquid and vapor densities, and associated
* thermodynamic and transport properties.
* If isat=0, checks whether T-D or T-P (iopt=1,2)
* falls on or within TOL of the liquid-vapor
* surface; if so, sets isat <- 1 and computes
* properties. */
void HGKsat(isat,iopt,itripl,Temp,Pres,Dens,epseqn)
int iopt,itripl,epseqn,*isat;
double *Temp,*Pres,Dens[2];
{
double dltemp,dvtemp,Ptemp;
extern double mfabs();
extern void pcorr(),tcorr();
if (*isat == 1)
{ if (iopt == 1)
pcorr(itripl,*Temp,Pres,&Dens[0],&Dens[1],epseqn);
else
tcorr(itripl,Temp,*Pres,&Dens[0],&Dens[1],epseqn);
}
else
{ if ((*Temp > crits.Tc) || (*Temp < tpoint.Ttripl) || ((iopt == 2) && (*Pres > crits.Pc)))
return;
else
{ pcorr(itripl,*Temp,&Ptemp,&dltemp,&dvtemp,epseqn);
if (((iopt == 2) && (mfabs(*Pres-Ptemp) <= tolers.PTOL)) || ((iopt == 1) && ((mfabs(Dens[0] - dltemp) <= tolers.DTOL) || (mfabs(Dens[0] - dvtemp) <= tolers.DTOL))))
{ *isat = 1;
*Pres = Ptemp;
Dens[0] = dltemp;
Dens[1] = dvtemp;
}
}
}
return;
}
/* -----------------------------------------------------------------
*** calcv3 - Compute the dependent state variable. */
void calcv3(iopt,itripl,Temp,Pres,Dens,epseqn)
int iopt,itripl,epseqn;
double Temp,*Pres,*Dens;
{
double ps,dll,dvv,dguess;
extern void pcorr(),resid(),base(),ideal();
if (iopt == 1)
{ resid(Temp,*Dens);
base(*Dens,Temp);
ideal(Temp);
*Pres = RTcurr.rt * *Dens * aconst.zb + qqqq.q;
}
else
{ if (Temp < aconst.tz)
pcorr(itripl,Temp,&ps,&dll,&dvv,epseqn);
else
{ ps = 20000.0;
dll = 0.0;
}
if (*Pres > ps)
dguess = dll;
else
dguess = *Pres / Temp / 0.4;
denHGK(Dens,*Pres,dguess,Temp,&fcts.dpdd);
ideal(Temp);
}
return;
}
/* -----------------------------------------------------------------
*** thmHGK - Computes thermodynamic functions in dimensionless
* units from the HGK equation of state: Helmholtz, Gibbs,
* internal energy, and enthalpy functions (ad, gd, ud, hd) are
* per RT; entropy and heat capacities (sd, cvd, cpd) are per R. */
void thmHGK(d,t)
double d,t;
{
double z;
z = aconst.zb + qqqq.q/RTcurr.rt/d;
fcts.dpdd = RTcurr.rt * (aconst.zb + aconst.yb * aconst.dzb) + qqqq.q5;
fcts.ad = basef.ab + resf.ar + idf.ai - aconst.uref/t + aconst.sref;
fcts.gd = fcts.ad + z;
fcts.ud = basef.ub + resf.ur + idf.ui - aconst.uref/t;
fcts.dpdt = RTcurr.rt * d * basef.dpdtb + resf.dpdtr;
fcts.cvd = basef.cvb + resf.cvr + idf.cvi;
fcts.cpd = fcts.cvd + t * fcts.dpdt * fcts.dpdt / (d * d * fcts.dpdd * aconst.gascon);
fcts.hd = fcts.ud + z;
fcts.sd = basef.sb + resf.sr + idf.si - aconst.sref;
fcts.dvdt = fcts.dpdt / fcts.dpdd / d / d;
fcts.cjtt = 1.0 / d - t * fcts.dvdt;
fcts.cjth = -fcts.cjtt / fcts.cpd / aconst.gascon;
return;
}
/* -----------------------------------------------------------------
*** bb - Computes molecular parameters b, the "excluded volume"
* (eq A.3), and B, the second virial coefficient (eq A.4),
* in cm3/g (b1,b2) and their first and second derivatives
* with respect to temperature (b1t,b1tt,b2t,b2tt). */
void bb(t)
double t;
{
double v[10];
int i;
v[0] = 1.0;
for (i=1; i < 10; i++)
v[i] = v[i-1] * aconst.tz / t;
ellcon.b1 = bconst.bp[0] + bconst.bp[1] * log(1.0 / v[1]);
ellcon.b2 = bconst.bq[0];
ellcon.b1t = bconst.bp[1] * v[1] / aconst.tz;
ellcon.b2t = ellcon.b1tt = ellcon.b2tt = 0.0;
for (i=2; i < 10; i++)
{ ellcon.b1 += bconst.bp[i] * v[i-1];
ellcon.b2 += bconst.bq[i] * v[i-1];
ellcon.b1t -= (i-1) * bconst.bp[i] * v[i-1] / t;
ellcon.b2t -= (i-1) * bconst.bq[i] * v[i-1] / t;
ellcon.b1tt += bconst.bp[i] * (i-1) * (i-1) * v[i-1] / t / t;
ellcon.b2tt += bconst.bq[i] * (i-1) * (i-1) * v[i-1] / t / t;
}
ellcon.b1tt -= ellcon.b1t / t;
ellcon.b2tt -= ellcon.b2t / t;
return;
}
/* -----------------------------------------------------------------
*** base - Computes Abase, Gbase, Sbase, Ubase, Hbase, Cvbase
* -- all per RT (dimensionless) -- as well as Pbase & dP/dT
* -- both per (DRT) -- for the base function (ab, gb, sb, ub,
* hb, cvb, pb, dpdtb). See Haar, Gallagher & Kell (1979), eq(1). */
void base(d,t)
double d,t;
{
double x,z0,dz0,bb2tt;
aconst.yb = 0.25 * ellcon.b1 * d;
x = 1.0 - aconst.yb;
z0 = (1.0 + ellcon.g1 * aconst.yb + ellcon.g2 * aconst.yb * aconst.yb) / (x * x * x);
aconst.zb = z0 + 4.0 * aconst.yb * (ellcon.b2/ellcon.b1 - ellcon.gf);
dz0 = (ellcon.g1 + 2.0 * ellcon.g2 * aconst.yb) / (x * x * x) + 3.0 * (1.0 + ellcon.g1 * aconst.yb + ellcon.g2 * aconst.yb * aconst.yb) / (x * x * x * x);
aconst.dzb = dz0 + 4.0 * (ellcon.b2/ellcon.b1 - ellcon.gf);
basef.pb = aconst.zb;
basef.ab = - log(x) - (ellcon.g2 - 1.0)/x + 28.16666667/x/x + 4.0*aconst.yb* (ellcon.b2 /ellcon.b1 - ellcon.gf) + 15.166666667 +log(d*t*aconst.gascon/0.101325);
basef.gb = basef.ab + aconst.zb;
basef.ub = -t * ellcon.b1t * (aconst.zb - 1.0 - d * ellcon.b2)/ellcon.b1 - d*t*ellcon.b2t;
basef.sb = basef.ub - basef.ab;
basef.hb = aconst.zb + basef.ub;
bb2tt = t * t * ellcon.b2tt;
basef.cvb = 2.0 * basef.ub + (z0 - 1.0)*(((t*ellcon.b1t/ellcon.b1)* (t*ellcon.b1t/ellcon.b1)) - t*t*ellcon.b1tt/ellcon.b1) - d*(bb2tt - ellcon.gf* ellcon.b1tt*t*t) - (t*ellcon.b1t/ellcon.b1)*(t*ellcon.b1t/ ellcon.b1)*aconst.yb*dz0;
basef.dpdtb = basef.pb/t + d*(aconst.dzb*ellcon.b1t/4.0 + ellcon.b2t - ellcon.b2 /ellcon.b1 *ellcon.b1t);
return;
}
/* -----------------------------------------------------------------
*** resid - Computes residual contributions to pressure (q), the
* Helmloltz function (ar) , dP/dD (q5), the Gibbs function
* (gr), entropy (sr), internal energy (ur), enthalpy (hr),
* isochoric heat capacity (cvr), and dP/dT. The first 36
* terms of the residual function represent a global
* least-squares fit to experimental data outside the
* critical region, terms 37-39 affect only the immediate
* vicinity of the critical point, and the last term (40)
* contributes only in the high pressure, low temperature
* region. */
void resid(t,d)
double t,d;
{
int i,k,l,km,zz;
double dadt,e,q10,q20,v,dfdt,d2f,dpt,q2a,ddz,del,ex1,dex,tex;
double att,tx,tau,ex2,qm,qp,fct,q5t,qr[11],qt[10];
qr[0] = qqqq.q5 = qqqq.q = resf.ar = dadt = resf.cvr = resf.dpdtr = 0.0;
e = exp(-aconst.aa * d);
qr[1] = q10 = d * d * e;
q20 = 1.0 - e;
v = aconst.tz / t;
qt[0] = t / aconst.tz;
for (i=1;i<10;i++)
{ qr[i+1] = qr[i] * q20;
qt[i] = qt[i-1] * v;
}
for (i=0;i<nconst.n;i++)
{ k = nconst.ii[i] + 1;
l = nconst.jj[i];
zz = k;
if (k == 1)
qp = nconst.g[i] * aconst.aa * qr[1] * qt[l];
else
qp = nconst.g[i] * aconst.aa * qr[k] * qt[l];
qqqq.q += qp;
qqqq.q5 += aconst.aa*(2.0/d - aconst.aa*(1.0 - e*(k-1)/q20))*qp;
resf.ar += nconst.g[i] * qr[k+1] * qt[l] / q10 / zz / RTcurr.rt;
dfdt = pow(q20,(double) k) * (1-l) * qt[l+1] / aconst.tz / k;
d2f = l * dfdt;
dpt = dfdt * q10 * aconst.aa * k / q20;
dadt += nconst.g[i] * dfdt;
resf.dpdtr += nconst.g[i] * dpt;
resf.cvr += nconst.g[i] * d2f / aconst.gascon;
}
qp = q2a = 0.0;
for (i=36;i<40;i++)
{ if (nconst.g[i] != 0.0)
{ k = nconst.ii[i];
km = nconst.jj[i];
ddz = addcon.adz[i-36];
del = d / ddz - 1.0;
if (mfabs(del) < 0.0000000001)
del = 0.0000000001;
ex1 = -addcon.aad[i-36] * pow(del,(double) k);
if (ex1 < tolers.EXPTOL)
dex = 0.0;
else
dex = exp(ex1) * pow(del,(double) km);
att = addcon.aat[i-36];
tx = addcon.atz[i-36];
tau = t / tx - 1.0;
ex2 = -att * tau * tau;
if (ex2 < tolers.EXPTOL)
tex = 0.0;
else
tex = exp(ex2);
q10 = dex * tex;
qm = km / del - k * addcon.aad[i-36] * pow(del,(double) (k-1));
fct = qm * d * d * q10 / ddz;
q5t = fct * (2.0 / d + qm / ddz) - (d / ddz) * (d / ddz) * q10 * (km / del / del + k*(k-1)*addcon.aad[i-36] * pow(del,(double) (k-2)));
qqqq.q5 += q5t * nconst.g[i];
qp += nconst.g[i] * fct;
dadt -= 2.0 * nconst.g[i] * att * tau * q10 / tx;
resf.dpdtr -= 2.0 * nconst.g[i] * att * tau * fct / tx;
q2a += t * nconst.g[i] * (4.0 * att * ex2 + 2.0 * att) * q10 / tx / tx;
resf.ar += q10 * nconst.g[i] / RTcurr.rt;
}
}
resf.sr = -dadt / aconst.gascon;
resf.ur = resf.sr + resf.ar;
resf.cvr += q2a / aconst.gascon;
qqqq.q += qp;
return;
}
/* -----------------------------------------------------------------
*** ideal - Computes thermodynamic properties for H2O in the
* ideal gas state using equations given by Woolley (1979). */
void ideal(t)
double t;
{
double tt,tl,emult;
double c[18]={19.730271018,20.9662681977,-.483429455355, 6.05743189245,22.56023885,-9.87532442,-4.3135538513,.458155781,-.047754901883,.0041238460633,-.00027929052852,.000014481695261,-.00000056473658748,.000000016200446,-.0000000003303822796,.00000000000451916067368,-.0000000000000370734122708,.000000000000000137546068238};
int i;
tt = t / 100.0;
tl = log(tt);
idf.gi = -(c[0]/tt + c[1]) * tl;
idf.hi = (c[1] + c[0]*(1 - tl)/tt);
idf.cpi = c[1] - c[0]/tt;
for (i=2;i<18;i++)
{ emult = pow(tt,(double)(i-5));
idf.gi -= c[i] * emult;
idf.hi += c[i] * (i-5) * emult;
idf.cpi += c[i] * (i-5) * (i-4) * emult;
}
idf.ai = idf.gi - 1;
idf.ui = idf.hi - 1;
idf.cvi = idf.cpi - 1;
idf.si = idf.ui - idf.ai;
return;
}
/* -----------------------------------------------------------------
*** dalHGK - Computes/returns (d(alpha)/dt)p(d,t,alpha)
* for the Haar et al. (1983) equation of state. */
double dalHGK(d,t,alpha)
double d,t,alpha;
{
double k,l,km,lm,kp,lp,x;
double dydtp,dbdd,db2dd,db2ddt,db2dtp,db3ddt,db3dtt,drdd,dr2dd,dr2ddt;
double e1,e2,tzt,dr2dtp,dr3ddt,dr3dtt,xtzt,ai,bi,di,ti,tau,del,ex1,dex;
double ex2,tex,ex12,qm,xdell,xdelk;
extern double mfabs();
int i;
/* evaluate derivatives for the base function */
aconst.yb = 0.25 * ellcon.b1 * d;
x = 1.0 - aconst.yb;
dydtp = (d/4.0)*(ellcon.b1t - ellcon.b1*alpha);
dbdd = aconst.gascon*t* ((ellcon.b1/4.0/x) * (1.0 - (ellcon.g2-1.0)/x + (ellcon.g1+ellcon.g2+1.0)/x/x) + ellcon.b2 - ellcon.b1*ellcon.gf + 1.0/d);
db2dd = aconst.gascon*t* ((ellcon.b1*ellcon.b1/16.0/x/x) * (1.0 - 2.0*(ellcon.g2-1.0)/x + 3.0 * (ellcon.g1+ellcon.g2+1.0)/x/x) - 1.0/d/d);
db2ddt = aconst.gascon*t * ((ellcon.b1t/4.0/x/x) * (1.0 - (ellcon.g2-1.0)*(1.0+aconst.yb)/x + (ellcon.g1+ellcon.g2+1.0) * (1.0+2.0*aconst.yb)/x/x) + ellcon.b2t - ellcon.gf*ellcon.b1t) + dbdd/t;
db2dtp = dbdd/t + aconst.gascon*t* ( (ellcon.b1*dydtp/4.0/x/x/x) * (1.0 - ellcon.g2 + 2.0 * (ellcon.g1+ellcon.g2+1.0)/x) + ((x*ellcon.b1t + ellcon.b1*dydtp)/4.0/x/x) * (1.0 - (ellcon.g2-1.0)/x + (ellcon.g1+ellcon.g2+1.0) /x/x) + ellcon.b2t - ellcon.gf*ellcon.b1t + alpha/d );
db3ddt = db2dd/t + aconst.gascon*t * ( (ellcon.b1*ellcon.b1*dydtp/8.0/x/x/x/x) * (1.0 - ellcon.g2 + 3.0 * (ellcon.g1+ellcon.g2+1.0)/x) + (ellcon.b1*(x*ellcon.b1t + ellcon.b1*dydtp)/8.0/x/x/x) * (1.0 - 2.0*(ellcon.g2-1.0)/x + 3.0 * (ellcon.g1+ellcon.g2+1.0)/x/x) - 2.0*alpha/d/d );
db3dtt = (db2ddt - dbdd/t)/t + aconst.gascon*t* ((ellcon.b1t*dydtp/2.0/x/x/x/x) * (1.0 - ellcon.g2 + (ellcon.g1+ellcon.g2+1.0)*(2.0+aconst.yb)/x) + ((x*ellcon.b1tt + 2.0*ellcon.b1t*dydtp)/4.0/x/x/x) * (1.0 - (ellcon.g2-1.0) *(1+aconst.yb)/x + (ellcon.g1+ellcon.g2+1.0)*(1.0+2.0*aconst.yb)/x/x) + ellcon.b2tt - ellcon.gf*ellcon.b1tt ) + (t*db2dtp - dbdd)/t/t;
/* evaluate derivatives for the residual function */
e1 = exp(-aconst.aa * d);
e2 = 1.0 - e1;
tzt = aconst.tz / t;
drdd = dr2dd = dr2ddt = dr2dtp = dr3ddt = dr3dtt = 0.0;
for (i=0;i<nconst.n;i++)
{ k = (double)(nconst.ii[i]) + 1.0;
l = (double)(nconst.jj[i]) - 1.0;
km = k - 1.0;
lm = l - 1.0;
kp = k + 1.0;
lp = l + 1.0;
xtzt = pow(tzt,l);
drdd += nconst.g[i] * xtzt*pow(e2,km)*e1;
dr2dd += nconst.g[i] * e1*xtzt*pow(e2,km) * (km*e1/e2 - 1.0);
dr2ddt -= nconst.g[i]*e1*l*pow(e2,km)*pow(tzt,lp)/aconst.tz;
dr2dtp += nconst.g[i]*e1*pow(e2,km)*xtzt * ( d*alpha - l/t - km*e1*d * alpha/e2 );
dr3ddt += nconst.g[i]*( km*d*alpha*e1*e1*xtzt* pow(e2,k-3.0) + e1* xtzt * pow(e2,km)* (km*e1/e2 - 1.0) * (d*alpha - l/t - km*d*alpha*e1/e2) );
dr3dtt += nconst.g[i]*l*e1*pow(e2,km)*pow(tzt,lp)/aconst.tz *( lp/t + d*alpha * km*e1/e2 - d*alpha );
}
for (i=36;i<40;i++)
{ k = (double)(nconst.ii[i]);
l = (double)(nconst.jj[i]);
km = k - 1.0;
lm = l - 1.0;
kp = k + 1.0;
lp = l + 1.0;
ai = addcon.aad[i-36];
bi = addcon.aat[i-36];
di = addcon.adz[i-36];
ti = addcon.atz[i-36];
tau = t/ti - 1.0;
del = d/di - 1.0;
if (mfabs(del) < 0.0000000001)
del = 0.0000000001;
ex1 = -ai * pow(del,k);
if (ex1 < tolers.EXPTOL)
dex = 0.0;
else
dex = exp(ex1);
ex2 = -bi * tau * tau;
if (ex2 <= tolers.EXPTOL)
tex = 0.0;
else
tex = exp(ex2);
ex12 = dex * tex;
qm = l/del - k*ai*pow(del,km);
xdell = pow(del,l);
xdelk = pow(del,k);
drdd += nconst.g[i]*xdell*ex12/di*qm;
dr2dd += nconst.g[i]*xdell*ex12/di/di * (qm*qm - l/di/di - ai*k*km* pow(del,k-2.0));
dr2ddt -= nconst.g[i]*2.0*bi*tau*ex12*xdell/ti/di*qm;
dr2dtp += nconst.g[i]/di*( d*alpha*xdell*ex12/di/del/del * (l + ai*k*km *xdelk) + qm * ( ex12 * ( xdell* (k*ai*d*alpha*pow(del,km)/di - 2.0*bi*tau/ti) - l*d*alpha*pow(del,lm)/di)));
dr3ddt += nconst.g[i]/di/di*( xdell*ex12* (2.0*qm*(l*d*alpha/di/del/del + ai*k*km*d*alpha*pow(del,k-2.0)/di) - 2.0*l*d*alpha/di/del/del/del + ai*k*km* (k-2.0)*pow(del,k-3.0)*d*alpha/di) + (qm*qm - l/del/del - ai*k*km* pow(del,k-2.0)) *(ex12*xdell*( ai*k*pow(del,k-1.0)*d*alpha/di - 2.0*bi*tau/ti ) - ex12*l*pow(del,l-1.0)*d*alpha/di));
dr3dtt -= 2.0*nconst.g[i]*bi/ti/di * ( tau*xdell*ex12*d* alpha/del/del / di * (l + ai*k*km*pow(del,k)) + qm*( xdell*ex12*( ai*k*d*alpha*tau*pow(del,km) /di + (1.0 - 2.0*bi*tau*tau)/ti - tau*l*d*alpha/di/del )));
}
return ((db3dtt + dr3dtt)*(2.0*(dbdd + drdd) + d*(db2dd + dr2dd)) - (db2ddt + dr2ddt)*(2.0*(db2dtp + dr2dtp) + d*(db3ddt + dr3ddt) - d*alpha*(db2dd + dr2dd))) / (2.0*(dbdd + drdd) + d*(db2dd + dr2dd)) / (2.0*(dbdd + drdd) + d*(db2dd + dr2dd));
/* returns (d(alpha)/dT)P */
}
/* -----------------------------------------------------------------
*** denHGK - Computes density (d in g/cm3) and dP/dD (dPdd) as
* f(p(MPa),t(degK)) from an initial density guess (dguess). */
void denHGK(d,p,dguess,t,dpdd)
double *d,*dpdd,p,dguess,t;
{
int i,except1;
double pp,dpdx,x,dp;
extern double mfabs();
i = 0;
*d = dguess;
do
{ i++;
except1 = 0;
if (*d <= 0.0)
*d = 0.00000001;
if (*d > 1.9)
*d = 1.9;
resid(t,*d);
base(*d,t);
pp = RTcurr.rt * (*d) * basef.pb + qqqq.q;
*dpdd = RTcurr.rt * (aconst.zb + aconst.yb * aconst.dzb) + qqqq.q5;
/* if dpdd < 0 assume d in 2-phase region and adjust accordingly */
if (*dpdd <= 0.0)
{ except1 = 1;
if (dguess >= 0.2967)
*d *= 1.02;
else
*d *= 0.98;
}
if ((i <= 10) && except1)
{ ; }
else
{ dpdx = *dpdd * 1.1;
if (dpdx < 0.1)
dpdx = 0.1;
dp = mfabs(1.0 - pp/p);
if ((dp < 0.00000001) || ((dguess > 0.3) && (dp < 0.0000001)) || ((dguess > 0.7) && (dp < 0.000001)))
return;
x = (p - pp) / dpdx;
if (mfabs(x) > 0.1)
x *= 0.1 / mfabs(x);
*d += x;
if (*d <= 0.0)
*d = 0.00000001;
}
} while (i <= 30);
return;
}
/* -----------------------------------------------------------------
*** PsHGK - Returns an approximation to Psaturation(T) that agrees
* to within 0.02% of that predicted by the HGK surface
* for temperatures up to within roughly a degree of
* the critical point. */
double PsHGK(t)
double t;
{
double a[8] = {-7.8889166,2.5514255,-6.716169,33.239495,-105.38479,174.35319, -148.39348,48.631602};
double val,v,w,pl,b,q;
extern double mfabs();
int i;
if (t <= 314.0)
{ pl = 6.3573118 - 8858.843 / t + 607.56335 * pow(t,-0.6);
val = 0.1 * exp(pl);
}
else
{ v = t / 647.25;
w = mfabs(1.0 - v);
b = 0.0;
for (i=0;i<8;i++)
b += a[i] * pow(w,(((double) i) + 2.0) / 2.0);
q = b / v;
val = 22.093 * exp(q);
}
return val;
}
/* -----------------------------------------------------------------
*** TsHGK - Returns Tsaturation(P). */
double TsHGK(p)
double p;
{
double val,pl,tg,pp,dp;
extern double PsHGK(),TdsdT(),mfabs();
int k,okay;
val = 0.0;
if (p > 22.05)
return val;
k = 0;
pl = 2.302585 + log(p);
tg = 372.83 + pl * (27.7589 + pl * (2.3819 + pl * (0.24834 + pl*0.0193855)));
okay = 0;
while (!okay)
{ if (tg < 273.15)
tg = 273.15;
if (tg > 674.0)
tg = 674.0;
if (k >= 8)
{ val = tg;
okay = 1;
}
else
{ k++;
pp = PsHGK(tg);
dp = TdPsdT(tg);
if (mfabs(1.0 - pp/p) < 0.00001)
{ val = tg;
okay = 1;
}
else
tg *= (1.0 + (p - pp)/dp);
}
}
return val;
}
/* -----------------------------------------------------------------
*** TdPsdT - Returns T*(dPsat/dT). */
double TdPsdT(t)
double t;
{
double a[8] = {-7.8889166,2.5514255,-6.716169,33.239495,-105.38479,174.35319,
-148.39348,48.631602};
double v,w,b,c,y,q;
int i;
v = t / 647.25;
w = 1.0 - v;
b = c = 0.0;
for (i=0;i<8;i++)
{ y = a[i] * pow(w,(((double) i) + 2.0) / 2.0);
c += y/w*(0.5 - 0.5*(i+1) - 1.0/v);
b += y;
}
q = b / v;
return 22.093 * exp(q) * c;
}
/* -----------------------------------------------------------------
*** corr - Computes liquid and vapor densities (dliq & dvap)
* and (Gl-Gv)/RT (delg) for T-P conditions on or
* near the saturation surface. */
void corr(itripl,t,p,dl,dv,delg,epseqn)
int itripl,epseqn;
double t,p,*dl,*dv,*delg;
{
double dguess,gl,gv,dpdd;
extern void denHGK(),thmHGK(),dimHGK(),ideal(),bb();
bb(t);
dguess = *dl;
if (*dl <= 0.0)
dguess = 1.11 - 0.0004 * t;
denHGK(dl,p,dguess,t,&dpdd);
ideal(t);
thmHGK(*dl,t);
/* save liquid properties */
dimHGK(1,itripl,t,p,*dl,epseqn);
gl = fcts.gd;
dguess = *dv;
if (*dv <= 0.0)
dguess = p / RTcurr.rt;
denHGK(dv,p,dguess,t,&dpdd);
if (*dv < 0.0000005)
*dv = 0.0000005;
ideal(t);
thmHGK(*dv,t);
/* *** vapor properties will be available
*** in COMMON /fcts/ (dimensionless) after
*** pcorr's final call of corr (delg < 10d-4) */
gv = fcts.gd;
*delg = gl - gv;
return;
}
/* -----------------------------------------------------------------
*** pcorr - Computes Psaturation(T) (p) and liquid and vapor
* densities (dl & dv) from refinement of an initial
* approximation (PsHGK(t)) in accord with Gl = Gv. */
void pcorr(itripl,t,p,dl,dv,epseqn)
int itripl,epseqn;
double t,*p,*dl,*dv;
{
double dp,delg;
extern double PsHGK(),mfabs();
extern void corr();
int okay;
*p = PsHGK(t);
*dl = *dv = 0.0;
okay = 0;
while (! okay)
{ corr(itripl,t,*p,dl,dv,&delg,epseqn);
dp = delg * aconst.gascon * t / (1.0/(*dv) - 1.0/(*dl));
*p += dp;
if (mfabs(delg) <= 0.0001)
okay = 1;
}
return;
}
/* -----------------------------------------------------------------
*** tcorr - Computes Tsaturation(P) (t) and liquid and vapor
* densities (dl & dv) from refinement of an initial
* approximation (TsHGK(p)) in accord with Gl = Gv. */
void tcorr(itripl,t,p,dl,dv,epseqn)
int itripl,epseqn;
double *t,p,*dl,*dv;
{
double dp,delg;
extern double TsHGK(),TdPsdT(),mfabs();
extern void corr();
int okay;
*t = TsHGK(p);
if (*t == 0.0)
return;
*dl = *dv = 0.0;
okay = 0;
while (!okay)
{ RTcurr.rt = *t * aconst.gascon;
corr(itripl,*t,p,dl,dv,&delg,epseqn);
dp = delg * aconst.gascon * (*t) / (1.0/(*dv) - 1.0/(*dl));
*t *= (1.0 - dp/TdPsdT(*t));
if (mfabs(delg) <= 0.0001)
okay = 1;
}
return;
}
/* -----------------------------------------------------------------
*** LVSeqn - Computes thermodynamic and transport properties of
* critical region H2O (369.85-419.85 degC,
* 0.20-0.42 gm/cm3) from the fundamental equation given
* by Levelt Sengers, et al (1983): J.Phys.Chem.Ref.Data,
* V.12, No.1, pp.1-28. */
void LVSeqn(isat,iopt,itripl,T,P,Dens,epseqn)
int *isat,*iopt,itripl,epseqn;
double *T,*P,Dens[2];
{
extern void LVSsat(),thmLVS(),dimLVS(),denLVS(),cpswrap();
double dl,dv,cdens;
int cpoint,okay,ioptsv = 0;
cpoint = okay = 0;
satur.DH2O = Dens[0];
while (!okay)
{ LVSsat(*iopt,isat,T,P,&satur.DH2O);
if ((*isat != 0) || (*iopt != 1))
denLVS(*isat,T,*P);
if (*isat == 0)
Dens[0] = satur.DH2O;
else
{ Dens[0] = satur.Dliq;
Dens[1] = satur.Dvap;
}
if (*isat == 0)
{ thmLVS(*isat,*T,param.r1,param.th1);
dimLVS(*isat,itripl,param.th1,*T,(*P*10.0),&dl,&dv,wpvals.wprops, epseqn);
if (cpoint)
{ cpswrap();
Dens[0] = Dens[1] = cdens;
*isat = 1;
*iopt = ioptsv;
}
okay = 1;
}
else
{ param.th1 = -1.0;
thmLVS(*isat,*T,param.r1,param.th1);
dimLVS(*isat,itripl,param.th1,*T,(*P*10.0),&dl,&dv,wpvals.wprops, epseqn);
param.th1 = 1.0;
thmLVS(*isat,*T,param.r1,param.th1);
dimLVS(*isat,itripl,param.th1,*T,(*P*10.0),&dl,&dv,wpvals.wprops, epseqn);
if (dl == dv)
{ cpoint = 1;
cdens = dl;
*T = 647.0670000003;
*P = 22.0460000008;
ioptsv = *iopt;
*iopt = 2;
*isat = 0;
}
else
okay = 1;
}
}
return;
}
/* -----------------------------------------------------------------
*** cpswrap - Load critical point A, G, U, H, S, Vs, Di, ZB,
* albe values from wpliq into wprops and
* approximations to critical Cv, Cp, alpha, beta,
* visc, tcond, Prndtl, tdiff, visck, YB, QB, XB,
* daldT, st values from wprops into wpliq. */
void cpswrap()
{
wpvals.wprops[0] = wpvals.wpliq[0];
wpvals.wprops[1] = wpvals.wpliq[1];
wpvals.wprops[2] = wpvals.wpliq[2];
wpvals.wprops[3] = wpvals.wpliq[3];
wpvals.wprops[4] = wpvals.wpliq[4];
wpvals.wprops[10] = wpvals.wpliq[10];
wpvals.wprops[18] = wpvals.wpliq[18];
wpvals.wprops[13] = wpvals.wpliq[13];
wpvals.wpliq[5] = wpvals.wprops[5];
wpvals.wpliq[6] = wpvals.wprops[6];
wpvals.wpliq[8] = wpvals.wprops[8];
wpvals.wpliq[9] = wpvals.wprops[9];
wpvals.wpliq[19] = wpvals.wprops[19];
wpvals.wpliq[20] = wpvals.wprops[20];
wpvals.wpliq[22] = wpvals.wprops[22];
wpvals.wpliq[12] = wpvals.wprops[12];
wpvals.wpliq[14] = wpvals.wprops[14];
wpvals.wpliq[15] = wpvals.wprops[15];
wpvals.wpliq[21] = wpvals.wprops[21];
wpvals.wpliq[17] = wpvals.wprops[17];
wpvals.wprops[7] = wpvals.wpliq[7] = 42.9352766443498 * units.fs;
wpvals.wprops[11] = wpvals.wpliq[11] = wpvals.wprops[16] = wpvals.wpliq[16] = 1000000.0;
return;
}
/* -----------------------------------------------------------------
*** LVSsat - If isat=1, computes Psat(T) or Tsat(P) (iopt=1,2).
* If isat=0, checks whether T-D or T-P (iopt=1,2)
* falls on or within TOL of the liq-vap surface; if so,
* isat <- 1 and T <- Tsat. */
void LVSsat(iopt,isat,T,P,D)
int iopt,*isat;
double *T,*P,*D;
{
double ERRTOL = 0.000000000001,TCTOL = 0.01,Tsat;
extern double mfabs(),TsLVS(),Pfind();
extern void backup(),restor();
if (*isat == 1)
{ if (iopt == 1)
*P = Pfind(*isat,T,D);
*T = TsLVS(*isat,*P);
}
else
{ if (iopt == 1)
*P = Pfind(*isat,T,D);
if ((*P - ERRTOL) > crits.Pc)
return;
else
{ backup();
Tsat = TsLVS(*isat,*P);
if (mfabs(Tsat-(*T)) < TCTOL)
{ *T = Tsat;
*isat = 1;
}
else
restor();
}
}
return;
}
/* -----------------------------------------------------------------
*** denLVS - Calculates DH2O(T,P) or Dvap,Dliq(T,P) from the
* Levelt Sengers, et al (1983) critical region
* equation of state. */
void denLVS(isat,T,P)
int isat;
double *T,P;
{
double s[2],sd[2],Pnext,Pdif,delD,Tw,dTw,rho1,rho2;
extern double Pfind(),mfabs();
extern void ss();
int i;
if (isat == 0)
{ satur.DH2O = crits.rhoc;
for (i=0;i<20;i++)
{ Pnext = Pfind(isat,T,&satur.DH2O);
Pdif = Pnext - P;
if (satur.iphase == 2)
{ if (mfabs(Pdif) <= 0.0)
return;
if (Pdif < 0.0)
satur.DH2O = coefs1.x[4];
else
satur.DH2O = coefs1.x[3];
}
else
{ delD = -Pdif/deri2.dPdD;
satur.DH2O += delD;
if (satur.DH2O < coefs1.x[3])
satur.DH2O = coefs1.x[3];
if (satur.DH2O > coefs1.x[4])
satur.DH2O = coefs1.x[4];
if (mfabs(delD/satur.DH2O) < 0.000001)
return;
}
}
}
else
{ Tw = -crits.Tc / *T;
dTw = 1.0 + Tw;
ss(param.r1,param.th1,s,sd);
rho1 = 1.0 + coefs1.q[8]*dTw + coefs1.a[0] * (s[0] + s[1]);
rho2 = coefs1.a[6]*pow(param.r1,coefs1.a[5]) + coefs1.a[11]* pow(param.r1,coefs1.q[15]);
satur.Dvap = crits.rhoc * (rho1 - rho2);
satur.Dliq = crits.rhoc * (rho1 + rho2);
}
return;
}
/* -----------------------------------------------------------------
*** TsLVS - Returns saturation T(P) */
double TsLVS(isat,P)
int isat;
double P;
{
double val,D,Pnext,dT;
extern double Pfind(),mfabs();
int i;
val = crits.Tc - 1.0;
D = crits.rhoc;
for (i=0;i<20;i++)
{ Pnext = Pfind(isat,&val,&D);
dT = (Pnext - P) / deri2.dPdT;
val -= dT;
if (val > crits.Tc)
val = crits.Tc;
else
{ if (mfabs(dT/val) < 0.00000001)
return val;
}
}
return val;
}
/* -----------------------------------------------------------------
*** Pfind - Returns P(T,D). Computes (dP/dD)T when invoked by SUB
* Dfind (isat=0) and (dP/dT)D when invoked by SUB TsLVS
* (isat=1). Also computes 1st & 2nd partial derivatives
* the singular part of the potential (Delta P tilde) that
* are used in SUB thmLVS. */
double Pfind(isat,T,D)
int isat;
double *T,*D;
{
double sd[2],val,tee,rho,tt1,tt2,Pw0,dPw0,da4,dPw;
double Uw,dPdTcd,dPwdTw,rho1,rhodi,Cvcoex,err,Pwmu;
extern double mfabs();
extern void aux(),conver(),ss();
/* EQUIV CRAP */
deriv.xk[0] = coefs1.a[6];
deriv.xk[1] = coefs1.a[11];
if (mfabs(*T-crits.Tc) < tolers.FPTOL)
*T = crits.Tc;
tee = (*T - crits.Tc) / crits.Tc;
deriv.Tw = -crits.Tc / (*T);
deriv.dTw = 1.0 + deriv.Tw;
if (isat == 0)
{ rho = *D / crits.rhoc;
conver(rho,tee,&deriv.amu,¶m.th1,¶m.r1,&rho1,deriv.s,&rhodi,&err);
}
else
{ param.th1 = abc2.th = -1.0;
param.r1 = abc2.r = deriv.dTw/(1.0 - coefs1.a[8]);
ss(param.r1,param.th1,deriv.s,sd);
rho = param.th1 * (coefs1.a[6]*pow(param.r1,coefs1.a[5]) + coefs1.a[11]*pow(param.r1, coefs1.q[15])) + coefs1.a[0]*(deriv.s[0] + deriv.s[1]);
rho += 1.0 + coefs1.q[8]*deriv.dTw;
deriv.amu = 0.0;
*D = rho * crits.rhoc;
}
tt1 = param.th1 * param.th1;
tt2 = tt1 * tt1;
Pw0 = 1.0 + deriv.dTw * (coefs1.a[4] + deriv.dTw*(coefs1.a[3] + deriv.dTw* coefs1.a[1]));
if (isat == 0)
Pwmu = deriv.amu*rhodi;
else
Pwmu = 0.0;
deriv.p0th = coefs1.q[10] + coefs1.q[11]*tt1 + coefs1.q[12]*tt2;
deriv.p1th = coefs1.q[17] + coefs1.q[18]*tt1 + coefs1.q[19]*tt2;
dPw0 = coefs1.a[6]* deriv.p0th*pow(param.r1,2.0-coefs1.q[9]);
da4 = coefs1.a[11] * deriv.p1th * pow(param.r1,2.0-coefs1.q[14]);
dPw = coefs1.a[9]*(dPw0+da4);
deriv.Pw = Pw0 + Pwmu + dPw;
val = deriv.Pw * crits.Pcon * (*T);
if (mfabs(param.th1) < 1.0)
satur.iphase = 1;
else
{ satur.iphase = 2;
deriv.dP0dT = coefs1.a[4] + deriv.dTw*(2.0*coefs1.a[3] + 3.0* coefs1.a[1]*deriv.dTw);
deriv.dM0dT = coefs1.a[13] + deriv.dTw*(2.0*coefs1.a[14] + 3.0* coefs1.a[15]*deriv.dTw);
Uw = deriv.dP0dT - rho*deriv.dM0dT + coefs1.q[8]*deriv.amu + deriv.s[0] + deriv.s[1];
dPdTcd = Uw + rho*deriv.dM0dT;
dPwdTw = deriv.Pw - deriv.Tw*dPdTcd;
deri2.dPdT = crits.Pcon * dPwdTw;
}
aux(param.r1,param.th1,&deriv.d2PdT2,&deriv.d2PdMT,&deriv.d2PdM2,coefs1.a[9], deriv.xk,sd,&Cvcoex);
if (satur.iphase == 1)
deri2.dPdD = crits.dPcon * (*D) * (*T) / deriv.d2PdM2;
return val;
}
/* -----------------------------------------------------------------
*** aux - Calculates some second derivatives of the
* anomalous part of the equation of state. */
void aux(r1,th1,d2PdT2,d2PdMT,d2PdM2,aa,xk,sd,Cvcoex)
double r1,th1,*d2PdT2,*d2PdMT,*d2PdM2,aa,xk[2],sd[2],*Cvcoex;
{
double s[2],w[2],y[2],z[2],coex[2],deli,ww,yy,zz,gamma,tt1,ter,g;
double alhi,beti,gami,a1,a2,a4,f1;
int i;
deli = ww = yy = zz = *Cvcoex = 0.0;
s[0] = coefs1.a[16] + coefs1.a[17]*th1*th1;
s[1] = coefs1.a[18] + coefs1.a[19]*th1*th1;
sd[0] = 2.0 * th1 * coefs1.a[17];
sd[1] = 2.0 * th1 * coefs1.a[19];
gamma = coefs1.a[5] * (coefs1.a[10] - 1.0);
tt1 = th1 * th1;
ter = 2.0 * coefs1.a[5] * coefs1.a[10] - 1.0;
g = (1.0 + (coefs1.a[8] * ter - 3.0) * tt1 - coefs1.a[8] * (ter - 2.0) * tt1 * tt1);
for (i=0;i<2;i++)
{ alhi = coefs1.q[9] - deli;
beti = coefs1.a[5] + deli;
gami = gamma - deli;
if (r1 != 0.0)
{ w[i] = (1.0 - alhi)*(1.0-3.0*tt1)*s[i]-coefs1.a[5]*coefs1.a[10] * (1.0-tt1)*th1*(sd[i]);
w[i] = (w[i] * pow(r1,-alhi)) * xk[i]/g;
ww += w[i];
y[i] = beti*(1.0 - 3.0*tt1)*th1 - coefs1.a[5]*coefs1.a[10]*(1.0 - tt1)*th1;
y[i] = (y[i] * pow(r1,beti-1.0)) * xk[i] / g;
yy += y[i];
z[i] = 1.0-coefs1.a[8]*(1.0-(2.0*beti))*tt1;
z[i] = (z[i] * pow(r1,-gami)) * xk[i] / g;
zz += z[i];
a1 = (coefs1.a[5]*(coefs1.a[10]-3.0)-3.0*deli-coefs1.a[8]*alhi*gami)/ (2.0*coefs1.a[8]*coefs1.a[8]*(2.0- alhi)*(1.0-alhi)*alhi);
a2 = -(1 +((coefs1.a[5]*(coefs1.a[10]-3.0)-3.0*deli-coefs1.a[8]*alhi* ter)/(2.0*coefs1.a[8]*(1.0 - alhi)*alhi)));
a4 = 1.0 + ((ter - 2.0)/(2.0*alhi));
f1 = a1 + a2 + a4;
coex[i]=((2.0-alhi)*(1.0-alhi)*pow(r1,-alhi)*f1*xk[i]);
*Cvcoex += coex[i];
}
deli = 0.5;
}
*d2PdT2 = aa * ww;
*d2PdMT = yy + aa * coefs1.a[0] * ww;
*d2PdM2 = zz/aa + 2.0 * coefs1.a[0] * yy + coefs1.a[0] * coefs1.a[0] * aa * ww;
return;
}
/* -----------------------------------------------------------------
*** conver - Transforms T,D to parametric variables r,theta
* according to the revised and scaled equations. */
void conver(rho,Tee,amu,th1,r1,rho1s,s1,rhodi,error1)
double rho,Tee,*amu,*th1,*r1,*rho1s,s1[2],*rhodi,*error1;
{
double sd[2],Tstar,dtstin,rhodit,drho,rho1co,twofaz,y1,den1,tt,rhoweg;
double den12,den2,slope,hold,error2,rho1s2;
extern void ss(),rtheta();
extern double mfabs(),copysign();
int i;
Tstar = Tee + 1.0;
dtstin = 1.0 - (1.0 / Tstar);
*r1 = dtstin;
if (dtstin <= 0.0)
{ *r1 = dtstin / (1.0 - coefs1.a[8]);
*th1 = 1.0;
}
else
*th1 = 0.0;
ss(*r1,*th1,s1,sd);
*rhodi = 1.0 + coefs1.q[8] * dtstin;
rhodit = *rhodi + coefs1.a[0] * s1[0] + coefs1.a[0] * s1[1];
drho = rho - rhodit;
*amu = 0.0;
if (dtstin <= 0.0)
{ rho1co = coefs1.a[6] * pow(*r1,coefs1.a[5]) + coefs1.a[11] * pow(*r1,coefs1.q[15]);
twofaz = rho1co;
if (mfabs(drho) <= twofaz)
{ *rho1s = copysign(rho1co,drho) + coefs1.a[0] * s1[0];
*th1 = copysign(1.0,drho);
*error1 = 1.0;
abc2.r = *r1;
abc2.th = *th1;
return;
}
}
if (drho == 0.0)
{ *th1 = 0.0;
*r1 = dtstin;
*rho1s = coefs1.a[0] * s1[0];
}
/* rule for first pass */
y1 = dtstin;
den1 = rho - rhodit;
rtheta(r1,th1,den1,y1);
tt = *th1 * *th1;
*amu = coefs1.a[9] * pow(*r1,coefs1.a[5]*coefs1.a[10]) * *th1 * (1.0 - tt);
y1 = dtstin + coefs1.a[0] * *amu;
ss(*r1,*th1,s1,sd);
rhoweg = coefs1.a[11] * pow(*r1,coefs1.q[15]) * *th1 + coefs1.a[0] * s1[1];
*rho1s = den1 + coefs1.a[0] * s1[0] + rhoweg;
*error1 = rho - *rhodi - *rho1s;
abc2.r = *r1;
abc2.th = *th1;
if (mfabs(*error1) < 0.00001)
return;
/* rule for second pass */
den12 = rho - *rhodi - coefs1.a[0] * s1[0] + rhoweg;
if (den12 == den1)
den12 = den1 - 0.000001;
rtheta(r1,th1,den12,y1);
tt = *th1 * *th1;
*amu = coefs1.a[9] * pow(*r1,coefs1.a[5]*coefs1.a[10]) * *th1 * (1.0 - tt);
y1 = dtstin + coefs1.a[0] * *amu;
ss(*r1,*th1,s1,sd);
rhoweg = coefs1.a[11] * pow(*r1,coefs1.q[15]) * *th1 + coefs1.a[0] * s1[1];
rho1s2 = den12 + coefs1.a[0] * s1[0] + rhoweg;
error2 = rho - *rhodi - rho1s2;
if (mfabs(error2) < 0.00001)
{ abc2.r = *r1;
abc2.th = *th1;
*error1 = error2;
*rho1s = rho1s2;
return;
}
/* rule for nth pass */
den2 = den12;
for (i=0;i<10;i++)
{ slope = (error2 - *error1)/(den2-den1);
hold = den2;
den2 = den1 - (*error1/slope);
rtheta(r1,th1,den2,y1);
tt = *th1 * *th1;
*amu = coefs1.a[9] * pow(*r1,(coefs1.a[5]*coefs1.a[10])) * *th1 * (1.0 - tt);
y1 = dtstin + coefs1.a[0] * *amu;
ss(*r1,*th1,s1,sd);
rhoweg = coefs1.a[11] * pow(*r1,coefs1.q[15]) * *th1 + coefs1.a[0] * s1[1];
*rho1s = den2 + coefs1.a[0] * s1[0] + rhoweg;
*error1 = error2;
error2 = rho - *rhodi - rho1s2;
abc2.r = *r1;
abc2.th = *th1;
if (mfabs(error2) < 0.000001)
return;
}
return;
}
/* -----------------------------------------------------------------
*** rtheta - Fits data for 1.0 < theta < 1.000001.
* Solves:
* rho = em*theta*(r**beta)
* Tee = r*(1.0d0-besq*theta*theta)
*
* Routine given by Moldover (1978): Jour. Res. NBS, v. 84, n. 4,
* p. 329 - 334. From Loki: Moldover's code used extensive goto-ing; I
* rewrote this using boolean checks and do-whiles */
void rtheta(r,theta,rho,Tee)
double *r,*theta,rho,Tee;
{
int do600,do496,i;
extern double copysign(),mfabs();
double bee,z,c,z2,z3,dz,absrho;
do600 = do496 = 0;
if ((coefs1.a[6] <= 0.0) || (coefs1.a[8] <= 1.0) || ((absrho=mfabs(rho)) < 0.000000000001))
{ if (mfabs(Tee) >= 0.000000000001)
do600 = 1;
}
else
{ bee = sqrt(coefs1.a[8]);
if (mfabs(Tee) < 0.000000000001)
{ *theta = copysign(1.0,rho)/bee;
*r = pow(rho/(coefs1.a[6]*(*theta)),1.0/coefs1.a[5]);
return;
}
if (Tee < 0.0)
{ z = 1.0 - (1.0 - bee) * Tee / (1.0 - coefs1.a[8]) * pow(coefs1.a[6]/absrho ,1.0/coefs1.a[5]); }
else
{ z = pow(1.0+Tee*pow(coefs1.a[6]/bee/absrho,1.0/coefs1.a[5]), -coefs1.a[5]); }
if (z > (1.00234 * bee))
{ do496 = do600 = 1; }
else
{ c = -rho * bee/coefs1.a[6]/pow(mfabs(Tee),coefs1.a[5]);
z = copysign(z,rho);
for (i=0;i<16;i++)
{ z2 = z * z;
z3 = 1.0 - z2;
dz = z3 * (z + c * pow(mfabs(z3), coefs1.a[5]))/(z3 + 2.0 * coefs1.a[5] * z2);
z -= dz;
if (mfabs(dz/z) < 0.000000000001)
{ *theta = z/bee;
*r = mfabs(Tee/(1.0 - z*z));
return;
}
}
}
}
if (!do600)
{ if (mfabs(*theta) > 1.0001)
*theta /= mfabs(*theta);
return;
}
if ((Tee < 0.0) || do496)
{ *theta = copysign(1.0,rho);
*r = mfabs(Tee/(1.0- coefs1.a[8]));
return;
}
*theta = 0.000000000001;
*r = Tee;
return;
}
/* -----------------------------------------------------------------
*** ss - Computes terms of the summation that defines dPotl/dT
* and the 1st derivative of the theta (s) square polynomial. */
void ss(r,th,s,sd)
double r,th,s[2],sd[2];
{
double sx[2],tt;
tt = th * th;
sx[0] = coefs1.a[16] + coefs1.a[17] * tt;
sd[0] = 2.0 * coefs1.a[17] * th;
sx[1] = coefs1.a[18] + coefs1.a[19] * tt;
sd[1] = 2.0 * coefs1.a[19] * th;
s[0] = sx[0] * coefs1.a[9] * coefs1.a[6] * pow(r,1.0-coefs1.q[9]);
s[1] = sx[1] * coefs1.a[9] * coefs1.a[11] * pow(r,1.0-coefs1.q[14]);
abc1.dPdM = pow(r,coefs1.a[5])*coefs1.a[6]*th + coefs1.a[0]* pow(r,1.0-coefs1.q[9])*coefs1.a[9]*coefs1.a[6]*sx[0] + pow(r,coefs1.q[15]) *coefs1.a[11]*th + coefs1.a[0]*pow(r,1.0-coefs1.q[14])*coefs1.a[9]* coefs1.a[11]*sx[1];
return;
}
/* -----------------------------------------------------------------
*** thmLVS - Calculates thermodynamic and transport properties
* of critical region H2O using the Levelt Sengers, et al
* (1983) equation of state. */
void thmLVS(isat,T,r1,th1)
int isat;
double T,r1,th1;
{
double sd[2],d2P0dT,d2M0dT,rho,Uw,dPdT2,dPwdTw,dPdTal,CviTw2;
double Cvw,Cpw,Hw,Sw,Scond,Cvcoex;
extern void aux(),ss();
d2P0dT = 2.0*coefs1.a[3] + 6.0*coefs1.a[1]*deriv.dTw;
d2M0dT = 2.0*coefs1.a[14] + 6.0*coefs1.a[15]*deriv.dTw;
deriv.dP0dT = coefs1.a[4] + deriv.dTw*(2.0*coefs1.a[3] + 3.0*coefs1.a[1]* deriv.dTw);
deriv.dM0dT = coefs1.a[13] + deriv.dTw*(2.0*coefs1.a[14] + 3.0*coefs1.a[15]* deriv.dTw);
if (isat == 0)
{ rho = satur.DH2O / crits.rhoc;
Uw = deriv.dP0dT - rho*deriv.dM0dT + coefs1.q[8]*deriv.amu + deriv.s[0] + deriv.s[1];
}
else
{ rho = th1 * (coefs1.a[6] * pow(r1,coefs1.a[5]) + coefs1.a[11] * pow(r1,coefs1.q[15])) + coefs1.a[0]*(deriv.s[0]+deriv.s[1]);
rho += 1.0 + coefs1.q[8]*deriv.dTw;
Uw = deriv.dP0dT - rho*deriv.dM0dT + coefs1.q[8] * deriv.amu + deriv.s[0] + deriv.s[1];
satur.DH2O = rho * crits.rhoc;
dPdT2 = deriv.Pw - deriv.Tw*(Uw + rho*deriv.dM0dT);
therm.heat = 1000.0*T*(crits.Pcon*dPdT2)*(1.0/satur.Dvap - 1.0/satur.Dliq);
ss(r1,th1,deriv.s,sd);
aux(r1,th1,&deriv.d2PdT2,&deriv.d2PdMT,&deriv.d2PdM2,coefs1.a[9], deriv.xk,sd,&Cvcoex);
if (r1 != 0.0)
deri2.dPdD = crits.dPcon * satur.DH2O * T / deriv.d2PdM2;
}
if (r1 != 0.0)
{ abc3.dPdTcd = deriv.dP0dT + coefs1.q[8]*(deriv.amu - rho/deriv.d2PdM2) + deriv.s[0] + deriv.s[1] - deriv.d2PdMT*rho/deriv.d2PdM2;
dPwdTw = deriv.Pw - deriv.Tw*abc3.dPdTcd;
dPdTal = crits.Pcon * dPwdTw;
CviTw2 = d2P0dT - rho*d2M0dT + deriv.d2PdT2 - (coefs1.q[8]+deriv.d2PdMT)* (coefs1.q[8] + deriv.d2PdMT)/deriv.d2PdM2;
Cvw = CviTw2 * deriv.Tw * deriv.Tw;
Cpw = Cvw + deriv.d2PdM2*dPwdTw*dPwdTw / (rho*rho);
therm.betaw = 1.0 / (satur.DH2O*deri2.dPdD);
therm.alphw = therm.betaw * dPdTal;
therm.Speed = 1000.0 * sqrt(Cpw/Cvw*deri2.dPdD);
}
else
{ Cvw = Cpw = therm.betaw = therm.alphw = 1.0;
therm.Speed = 0.0;
}
Hw = deriv.Pw - deriv.Tw*Uw;
Sw = Hw - rho*(deriv.amu+coefs1.a[12]+deriv.dTw*(coefs1.a[13]+deriv.dTw* (coefs1.a[14]+deriv.dTw*coefs1.a[15])));
Scond = crits.Scon / satur.DH2O;
therm.U = Uw * crits.Ucon / satur.DH2O;
therm.H = Hw * Scond * T;
therm.Entrop = Sw * Scond;
therm.AE = therm.U - T * therm.Entrop;
therm.GE = therm.H - T * therm.Entrop;
therm.Cv = Cvw * Scond;
therm.Cp = Cpw * Scond;
return;
}
/* -----------------------------------------------------------------
*** dalLVS - Computes/returns (d(alpha)/dt)p(D,T,alpha)
* for the Levelt Sengers et al. (1983)
* equation of state. Note that D (kg/m**3),
* T (degK), P (MPa), alpha (degK**-1). */
double dalLVS(D,T,P,alpha)
double D,T,P,alpha;
{
double dsdT[2],sp[2],dspdT[2],k[2],calpha[2],cbeta[2],delT,b1,b2;
double cgamma[2],u[2],v[2],w[2],dudT[2],dvdT[2],dwdT[2],ar1,a01,ar2,a02;
double amult,d0dT,drdT,q,dqdT,dP0dTT,ddelMT,dPdTT,dPdMMT,dPdMTT,dPPTT,pterm;
extern double mfabs();
int i;
if (abc2.r == 0.0)
return 1000000.0;
k[0] = coefs1.a[6];
k[1] = coefs1.a[11];
calpha[0] = coefs1.q[9];
calpha[1] = coefs1.q[14];
cbeta[0] = coefs1.a[5];
cbeta[1] = coefs1.q[15];
cgamma[0] = cbeta[0] * (coefs1.a[10] - 1.0);
cgamma[1] = cgamma[0] - coefs1.q[13];
delT = (T - crits.Tc) / T;
deriv.s[0] = coefs1.a[16] + coefs1.a[17] * abc2.th * abc2.th;
deriv.s[1] = coefs1.a[18] + coefs1.a[19] * abc2.th * abc2.th;
sp[0] = 2.0 * coefs1.a[17] * abc2.th;
sp[1] = 2.0 * coefs1.a[19] * abc2.th;
/* *********************************************************************
***
*** Compute drdT and d0dT from solution of the linear system
***
*** ax = b
***
*** d(dPdM)/dT = -D/Dc*alpha - P11*Tc/T**2 = ar1*drdT + a01*d0dT = b1
*** d(delT)/dT = Tc/T**2 = ar2*drdT + a02*d0dT = b2
*** */
b1 = -D / crits.rhoc * alpha - coefs1.q[8]*crits.Tc/T/T;
b2 = crits.Tc / (T * T);
ar1 = a01 = 0.0;
for (i=0;i<2;i++)
{ ar1 += k[i]*(cbeta[i]*abc2.th*pow(abc2.r,cbeta[i]-1.0) + coefs1.a[9]*coefs1.a[0]*(1.0 - calpha[i])* pow(abc2.r,-calpha[i])* deriv.s[i]);
a01 += k[i]*(pow(abc2.r,cbeta[i]) + coefs1.a[9]*coefs1.a[0]*sp[i]*pow(abc2.r,1.0 - calpha[i]));
}
ar2 = 1.0 - coefs1.a[8]*abc2.th*abc2.th - coefs1.a[9]*coefs1.a[0]*cbeta[0]*coefs1.a[10]*(1.0-abc2.th*abc2.th)*abc2.th *pow(abc2.r,(cbeta[0]*coefs1.a[10] - 1.0));
a02 = 3.0*coefs1.a[9]*coefs1.a[0]*abc2.th*abc2.th*pow(abc2.r,cbeta[0]*coefs1.a[10]) - 2.0*coefs1.a[8]*abc2.r*abc2.th - coefs1.a[9]*coefs1.a[0]*pow(abc2.r,cbeta[0] * coefs1.a[10]);
/* *********************************************************************
*** solve the linear system with simplistic GE w/ partial pivoting
********************************************************************* */
if (mfabs(ar1) > mfabs(ar2))
{ amult = -ar2 / ar1;
d0dT = (b2 + amult*b1) / (a02 + amult*a01);
drdT = (b1 - a01*d0dT) / ar1;
}
else
{
amult = -ar1 / ar2;
d0dT = (b1 + amult*b2) / (a01 + amult*a02);
drdT = (b2 - a02*d0dT) / ar2;
}
/* *********************************************************************
***
*** Compute theta polynomials and their tempertaure derivatives
*** */
dspdT[0] = 2.0*coefs1.a[17]*d0dT;
dspdT[1] = 2.0*coefs1.a[19]*d0dT;
dsdT[0] = dspdT[0]*abc2.th;
dsdT[1] = dspdT[1]*abc2.th;
q = 1.0 + (coefs1.a[8]*(2.0*cbeta[0]*coefs1.a[10]-1.0)-3.0)*abc2.th*abc2.th - coefs1.a[8]*(2.0*cbeta[0]*coefs1.a[10] - 3.0)*abc2.th*abc2.th*abc2.th*abc2.th;
dqdT = 2.0*(coefs1.a[8]*(2.0*cbeta[0]*coefs1.a[10]-1.0)-3.0)*abc2.th*d0dT - 4.0*coefs1.a[8]*(2.0*cbeta[0]* coefs1.a[10] -3.0)*abc2.th*abc2.th*abc2.th*d0dT;
for (i=0;i<2;i++)
{ u[i] = (1.0-coefs1.a[8]*(1.0-2.0*cbeta[i])*abc2.th*abc2.th)/q;
dudT[i] = (-2.0*coefs1.a[8]*(1.0-2.0*cbeta[i])*abc2.th*d0dT - u[i]*dqdT)/q;
v[i]= ((cbeta[i]-cbeta[0]*coefs1.a[10])*abc2.th+(cbeta[0]*coefs1.a[10]-3.0*cbeta[i])*abc2.th*abc2.th*abc2.th)/q;
dvdT[i] = ((cbeta[i]-cbeta[0]*coefs1.a[10])*d0dT+3.0*(cbeta[0]*coefs1.a[10]-3.0*cbeta[i])*abc2.th*abc2.th*d0dT - v[i]*dqdT) / q;
w[i] = ((1.0 - calpha[i])*(1.0-3.0*abc2.th*abc2.th)*deriv.s[i] - cbeta[0]*coefs1.a[10]*(abc2.th-abc2.th*abc2.th*abc2.th)* sp[i]) / q;
dwdT[i] = ((1.0-calpha[i])*((1.0-3.0*abc2.th*abc2.th)*dsdT[i]-6.0*abc2.th*deriv.s[i]*d0dT)-cbeta[0] *coefs1.a[10]*((abc2.th-abc2.th*abc2.th*abc2.th)*dspdT[i] + sp[i]*(d0dT - 3.0*abc2.th*abc2.th*d0dT)) - w[i]*dqdT) / q;
}
/* *********************************************************************
***
*** Compute dP0dTT, ddelMT, dPdTT, dPdMMT, dPdMTT, dPPTT
*** */
dP0dTT = crits.Tc/(T*T) * (2.0*coefs1.a[3] + 6.0*coefs1.a[1]*delT);
ddelMT = coefs1.a[9]*pow(abc2.r,cbeta[0]*coefs1.a[10])*(cbeta[0]*coefs1.a[10]*abc2.th/abc2.r*(1.0-abc2.th*abc2.th)*drdT + (1.0-3.0*abc2.th*abc2.th)*d0dT);
dPdTT = dPdMMT = dPdMTT = 0.0;
for (i=0;i<2;i++)
{ dPdTT += coefs1.a[9]*k[i] * (pow(abc2.r,1.0-calpha[i])*dsdT[i] + deriv.s[i]*(1.0 - calpha[i])*pow(abc2.r,-calpha[i])*drdT);
dPdMMT += k[i] * ((pow(abc2.r,-cgamma[i])*dudT[i] - u[i]*cgamma[i]*pow(abc2.r, -1.0-cgamma[i])*drdT) / coefs1.a[9] + 2.0*coefs1.a[0]*(pow(abc2.r,cbeta[i]-1.0)*dvdT[i] + v[i]*(cbeta[i] - 1.0)*pow(abc2.r,cbeta[i]-2.0)*drdT) + coefs1.a[9]*coefs1.a[0]*coefs1.a[0]*(pow(abc2.r,-calpha[i])*dwdT[i] - calpha[i]*w[i] *pow(abc2.r,-1.0-calpha[i])*drdT));
dPdMTT += k[i] * (pow(abc2.r,cbeta[i]-1.0)*dvdT[i] +v[i]*(cbeta[i] - 1.0)* pow(abc2.r,cbeta[i]-2.0)*drdT + coefs1.a[9]*coefs1.a[0]*(pow(abc2.r,-calpha[i])*dwdT[i] - calpha[i]*pow(abc2.r,-1.0 - calpha[i])*drdT*w[i]));
}
dPPTT = dP0dTT + dPdTT + coefs1.q[8]*ddelMT - D/crits.rhoc*dPdMTT/deriv.d2PdM2 +(coefs1.q[8] + deriv.d2PdMT)*(D/crits.rhoc*alpha/deriv.d2PdM2+D/crits.rhoc*dPdMMT /(deriv.d2PdM2*deriv.d2PdM2));
pterm = P/crits.Pc + abc3.dPdTcd;
/* return (d(alpha)/dT)P */
return (crits.Tc*crits.rhoc*crits.rhoc/(D*D*T*T) * (-2.0/T*deriv.d2PdM2*pterm + 2.0*alpha*deriv.d2PdM2*pterm + pterm*dPdMMT +deriv.d2PdM2*dPPTT));
}
/* -----------------------------------------------------------------
*** backup - Save Pfind COMMON values during saturation check. */
void backup()
{
store.isav1 = satur.iphase;
store.sav2 = param.r1;
store.sav3 = param.th1;
store.sav4 = deriv.amu;
store.sav5 = deriv.s[0];
store.sav6 = deriv.s[1];
store.sav7 = deriv.Pw;
store.sav8 = deriv.Tw;
store.sav9 = deriv.dTw;
store.sav10 = deriv.dM0dT;
store.sav11 = deriv.dP0dT;
store.sav12 = deriv.d2PdM2;
store.sav13 = deriv.d2PdMT;
store.sav14 = deriv.d2PdT2;
store.sav15 = deriv.p0th;
store.sav16 = deriv.p1th;
store.sav17 = deriv.xk[0];
store.sav18 = deriv.xk[1];
store.sav19 = deri2.dPdD;
return;
}
/* -----------------------------------------------------------------
*** restor - Restore Pfind COMMON values after saturation check. */
void restor()
{
satur.iphase = store.isav1;
param.r1 = store.sav2;
param.th1 = store.sav3;
deriv.amu = store.sav4;
deriv.s[0] = store.sav5;
deriv.s[1] = store.sav6;
deriv.Pw = store.sav7;
deriv.Tw = store.sav8;
deriv.dTw = store.sav9;
deriv.dM0dT = store.sav10;
deriv.dP0dT = store.sav11;
deriv.d2PdM2 = store.sav12;
deriv.d2PdMT = store.sav13;
deriv.d2PdT2 = store.sav14;
deriv.p0th = store.sav15;
deriv.p1th = store.sav16;
deriv.xk[0] = store.sav17;
deriv.xk[1] = store.sav18;
deri2.dPdD = store.sav19;
return;
}
/* -----------------------------------------------------------------
*** load - Load thermodynamic and transport property values from
* ptemp into props. */
void load(phase,ptemp,props)
int phase;
double ptemp[23],props[46];
/* The indexing system was switched over from a two-dimensional array to a one-dimensional array. */
{
int i,key[46] = {1,3,5,7,9,11,13,15,17,19,21,23,25,27,29,31,33,35,37,39,41,43,45,2,4,6,8,10,12,14,16,18,20,22,24,26,28,30,32,34,36,38,40,42,44,46};
for (i=0;i<23;i++)
props[key[(phase-1)*23 + i] - 1] = ptemp[i];
return;
}
/* -----------------------------------------------------------------
*** tpset - Dimension triple point U, S, H, A, G values (in J/g from
* Table 2, Helgeson & Kirkham, 1974a) into user-specified units. */
void tpset()
{
double Utr= -15766.0,Str= 3.5144,Htr= -15971.0,Atr= -12870.0,Gtr= -13073.0;
tpoint.Utripl = Utr*units.fh;
tpoint.Stripl = Str*units.fh;
tpoint.Htripl = Htr*units.fh;
tpoint.Atripl = Atr*units.fh;
tpoint.Gtripl = Gtr*units.fh;
return;
}
/* -----------------------------------------------------------------
*** triple - Convert U, S, H, A, G values computed with reference to
* zero triple point properties (Haar et al., 1984;
* Levelt Sengers et al., 1983) into values referenced to
* triple point properties given by Helgeson and Kirkham, 1974a. */
void triple(T,wpzero)
double T,wpzero[23];
{
int A=0,G=1,S=2,U=3,H=4;
double TS;
wpzero[S] += tpoint.Stripl;
TS = T * wpzero[S] - tpoint.Ttripl*tpoint.Stripl;
wpzero[G] = wpzero[H] - TS + tpoint.Gtripl;
wpzero[A] = wpzero[U] - TS + tpoint.Atripl;
wpzero[H] += tpoint.Htripl;
wpzero[U] += tpoint.Utripl;
return;
}
/* -----------------------------------------------------------------
*** TdegK - Returns input temperature t converted from
* user-specified units to degK. */
double TdegK(it,t)
int it;
double t;
{
switch (it)
{ case 1 : return t;
case 2 : return t + 273.15;
case 3 : return t / 1.8;
default: return (t+459.67)/1.8;
}
}
/* -----------------------------------------------------------------
*** TdegUS - Returns input temperature t converted
* from degK to user-specified units. */
double TdegUS(it,t)
int it;
double t;
{
switch (it)
{ case 1 : return t;
case 2 : return t - 273.15;
case 3 : return t * 1.8;
default: return t * 1.8 - 459.67;
}
}
/* -----------------------------------------------------------------
*** dimHGK - Dimensions thermodynamic and transport property values
* computed from the HGK equation of state per user-specified
* choice of units. */
void dimHGK(isat,itripl,t,p,d,epseqn)
int isat,itripl,epseqn;
double t,p,d;
{
double pbars,dkgm3,betaPa,betab,CpJKkg;
extern double viscos(),thcond(),surten(),dalHGK(),mfabs();
extern void Born92(),triple();
int aw=0,gw=1,sw=2,uw=3,hw=4,cvw=5,cpw=6,vsw=7,alw=8,bew=9,diw=10,viw=11,tcw=12, stw=13,tdw=14,Prw=15,vikw=16,albew=17,ZBw=18,YBw=19,QBw=20,dalwdT=21,XBw=22;
wpvals.wprops[aw] = fcts.ad*RTcurr.rt*units.fh;
wpvals.wprops[gw] = fcts.gd*RTcurr.rt*units.fh;
wpvals.wprops[sw] = fcts.sd*aconst.gascon*units.fh*units.ft;
wpvals.wprops[uw] = fcts.ud*RTcurr.rt*units.fh;
wpvals.wprops[hw] = fcts.hd*RTcurr.rt*units.fh;
wpvals.wprops[cvw] = fcts.cvd*aconst.gascon*units.fh*units.ft;
wpvals.wprops[cpw] = fcts.cpd*aconst.gascon*units.fh*units.ft;
wpvals.wprops[vsw] = sqrt(mfabs(fcts.cpd*fcts.dpdd*1000.0/fcts.cvd)) * units.fs;
wpvals.wprops[bew] = 1.0/(d*fcts.dpdd*units.fp);
wpvals.wprops[alw] = d * fcts.dvdt;
wpvals.wprops[dalwdT] = dalHGK(d,t,wpvals.wprops[alw]);
pbars = p*10.0;
dkgm3 = d * 1000.0;
betaPa = wpvals.wprops[bew]*units.fp / 1000000.0;
betab = wpvals.wprops[bew]*units.fp / 10.0;
CpJKkg = wpvals.wprops[cpw]/units.fh/units.ft * 1000.0;
wpvals.wprops[viw] = viscos(t,pbars,dkgm3,betaPa) * units.fvd;
wpvals.wprops[tcw] = thcond(t,pbars,dkgm3,wpvals.wprops[alw],betaPa) * units.fc * units.ft;
if ((isat == 0) || (isat == 2))
wpvals.wprops[stw] = 0.0;
else
wpvals.wprops[stw] = surten(t) * units.fst;
Born92(t,pbars,(dkgm3/1000.0),betab,wpvals.wprops[alw],wpvals.wprops[dalwdT],&wpvals.wprops[diw],&wpvals.wprops[ZBw],&wpvals.wprops[QBw],&wpvals.wprops[YBw],&wpvals.wprops[XBw],epseqn);
wpvals.wprops[tdw] = wpvals.wprops[tcw]/units.fc/units.ft / (dkgm3 * CpJKkg) * units.fvk;
if (wpvals.wprops[tcw] != 0.0)
wpvals.wprops[Prw] = wpvals.wprops[viw]/units.fvd * CpJKkg / (wpvals.wprops[tcw]/units.fc/units.ft);
else
wpvals.wprops[Prw] = 0.0;
wpvals.wprops[vikw] = wpvals.wprops[viw]/ units.fvd / dkgm3 * units.fvk;
wpvals.wprops[albew] = wpvals.wprops[alw] / wpvals.wprops[bew];
if (itripl == 1)
triple(t,wpvals.wprops);
return;
}
/* -----------------------------------------------------------------
*** dimLVS - Dimension critical region properties per user-specs
* and load into tprops. */
void dimLVS(isat,itripl,theta,T,Pbars,dl,dv,tprops,epseqn)
int isat,itripl,epseqn;
double theta,T,Pbars,*dl,*dv,tprops[23];
{
double CpJKkg,betaPa,betab,dkgm3;
extern double mfabs(),dalLVS(),visoc(),thcond(),surten();
extern void triple(),Born92();
int aw=0,gw=1,sw=2,uw=3,hw=4,cvw=5,cpw=6,vsw=7,alw=8,bew=9,diw=10,viw=11,tcw=12,stw=13,tdw=14,Prw=15,vikw=16,albew=17,ZBw=18,YBw=19,QBw=20,dalwdT=21,XBw=22;
if (isat == 1)
{ *dv = satur.Dvap;
*dl = satur.Dliq;
}
tprops[aw] = therm.AE * units.fh;
tprops[gw] = therm.GE * units.fh;
tprops[sw] = therm.Entrop * units.fh * units.ft;
tprops[uw] = therm.U * units.fh;
tprops[hw] = therm.H * units.fh;
tprops[cvw] = therm.Cv * units.fh * units.ft;
tprops[cpw] = therm.Cp * units.fh * units.ft;
tprops[vsw] = therm.Speed * units.fs;
tprops[bew] = therm.betaw / units.fp;
tprops[alw] = therm.alphw;
abc2.th = theta;
tprops[dalwdT] = dalLVS(satur.DH2O,T,Pbars/10.0,tprops[alw]);
CpJKkg = therm.Cp * 1000.0;
betaPa = therm.betaw / 1000000.0;
betab = therm.betaw / 10.0;
if (mfabs(theta) != 1.0)
{ dkgm3 = satur.DH2O;
tprops[stw] = 0.0;
}
else
{ if (theta < 0.0)
{ dkgm3 = satur.Dvap;
tprops[stw] = 0.0;
}
else
{ dkgm3 = satur.Dliq;
tprops[stw] = surten(T) * units.fst;
}
}
Born92(T,Pbars,(dkgm3/1000.0),betab,tprops[alw],tprops[dalwdT],&tprops[diw], &tprops[ZBw],&tprops[QBw],&tprops[YBw],&tprops[XBw],epseqn);
tprops[viw] = viscos(T,Pbars,dkgm3,betaPa) * units.fvd;
tprops[tcw] = thcond(T,Pbars,dkgm3,tprops[alw],betaPa) * units.fc * units.ft;
tprops[tdw] = tprops[tcw]/units.fc/units.ft / (dkgm3 * CpJKkg) * units.fvk;
tprops[Prw] = tprops[viw]/units.fvd * CpJKkg / (tprops[tcw]/units.fc/units.ft);
tprops[vikw] = tprops[viw]/units.fvd / dkgm3 * units.fvk;
tprops[albew] = tprops[alw] / tprops[bew];
if (itripl == 1)
triple(T,tprops);
return;
}
/* -----------------------------------------------------------------
*** viscos - Returns dynamic viscosity of H2O in kg/m*s (= Pa*s)
* if Tk, Pbars falls within the validity region (specified
* by the initial IF statement) of the Watson et al. (1980)
* equation; otherwise returns zero. See equations 3.1-2 and
* 4.1-5 and Tables 1, 6, and 8 from Sengers and
* Kamgar-Parsi (1984). */
double viscos(Tk,Pbars,Dkgm3,betaPa)
double Tk,Pbars,Dkgm3,betaPa;
{
double a[4] = {0.0181583,0.0177624,0.0105287,-0.0036744};
double b[42]= {0.5132047, 0.3205656, 0.0,0.0,
-0.7782567, 0.1885447, 0.2151778, 0.7317883,
1.2410440, 1.4767830, 0.0,0.0,
-0.2818107, -1.0707860, -1.2631840, 0.0,
0.0,0.0,0.1778064, 0.4605040,
0.2340379, -0.4924179, 0.0,0.0,
-0.0417661, 0.0,0.0,0.1600435,
0.0,0.0,0.0, -0.01578386,
0.0,0.0,0.0,0.0,0.0,0.0,0.0,-0.003629481,0.0,0.0};
double tol = 0.01,TdegC,T,D,sum,u0,u1,xt,u2;
int i,j;
TdegC = Tk - 273.15;
if ((Pbars > (5000.0+tol)) || ((Pbars > (3500.0+tol)) && (TdegC > (150.0+tol))) || ((Pbars > (3000.0+tol)) && (TdegC > (600.0+tol))) || (TdegC > (900.0+tol)))
return 0.0;
T = Tk / 647.27;
D = Dkgm3 / 317.763;
sum = 0.0;
for (i=0;i<4;i++)
sum += a[i]/pow(T,(double) i);
u0 = 0.000001 * sqrt(T) / sum;
sum = 0.0;
for (i=0;i<6;i++)
for (j=0;j<7;j++)
sum += b[j*6+i] * pow((1.0/T-1),(double) i) * pow((D-1),(double) j);
u1 = exp(D*sum);
if ((0.997 <= T) && (T <= 1.0082) && (0.755 <= D) && (D <= 1.2900))
{ xt = 22115000.0 / (100973.324169) * betaPa * Dkgm3 * Dkgm3;
if (xt < 22.0)
u2 = 1.0;
else
u2 = 0.922 * pow(xt,0.0263);
}
else
u2 = 1.0;
return u0*u1*u2;
}
/* -----------------------------------------------------------------
*** thcond - Returns thermal conductivity of H2O in J/m*deg*s (=W/m*deg)
* if Tk, Pbars falls within the validity region (specified
* by the initial IF statement) of the Sengers et al. (1984)
* equation; returns zero otherwise. See equations 3.2-14
* and tables 2-5 and I.5-6 from the above reference. */
double thcond(Tk,Pbars,Dkgm3,alph,betaPa)
double Tk,Pbars,Dkgm3,alph,betaPa;
{
int i,j;
double T,D,sum,L0,L1,L2,TdegC,u0,u1,xt,dPdT;
double aL[4] = { 2.02223,14.11166,5.25597,-2.0187 };
double au[4] = { .0181583,.0177624,.0105287,-.0036744 };
double bL[30] = { 1.329304600, -0.404524370, 0.244094900,
0.018660751, -0.129610680, 0.044809953,
1.701836300, -2.215684500, 1.651105700,
-0.767360020, 0.372833440, -0.112031600,
5.224615800, -10.12411100, 4.987468700,
-0.272976940, -0.430833930, 0.133338490,
8.712767500, -9.500061100, 4.378660600,
-0.917837820, 0.0, 0.0,
-1.852599900, 0.934046900, 0.0,
0.0, 0.0, 0.0};
double bu[30] = { 0.5019380, 0.2356220, -0.2746370, 0.1458310,
-0.0270448, 0.1628880, 0.7893930, -0.7435390,
0.2631290, -0.0253093, -0.1303560, 0.6736650,
-0.9594560, 0.3472470, -0.0267758, 0.9079190,
1.2075520, -0.6873430, 0.2134860, -0.0822904,
-0.5511190, 0.0670665, -0.4970890, 0.1007540,
0.0602253, 0.1465430, -0.0843370, 0.1952860,
-0.0329320, -0.0202595};
double tol = .01;
TdegC = Tk - 273.15;
if ((Pbars > (4000.0+tol)) || ((Pbars > (2000.0+tol)) && (TdegC > (125.0+tol))) || ((Pbars > (1500.0+tol)) && (TdegC > (400.0+tol))) || (TdegC > (800.0+tol)))
return 0.0;
T = Tk / 647.27;
D = Dkgm3 / 317.763;
sum = 0.0;
for (i=0;i<4;i++)
sum += aL[i]/pow(T,(double) i);
L0 = sqrt(T) / sum;
sum = 0.0;
for (i=0;i<5;i++)
for (j=0;j<6;j++)
sum += bL[i*6+j] * pow((1.0/T-1),(double) i) * pow((D-1),(double) j);
L1 = exp(D*sum);
sum = 0.0;
for (i=0;i<4;i++)
sum += au[i]/pow(T,(double) i);
u0 = 0.000001 * sqrt(T) / sum;
sum = 0.0;
for (i=0;i<6;i++)
for (j=0;j<5;j++)
sum += bu[i*5+j] * pow((1.0/T-1),(double) i) * pow((D-1),(double) j);
u1 = exp(D*sum);
xt = 22115000.0 / (100973.324169) * betaPa * Dkgm3 * Dkgm3;
dPdT = (647.27/22115000.0) * alph/betaPa;
L2 = 0.000000037711 / (u0*u1) * T*T/(D*D) * dPdT * dPdT * pow(xt,0.4678) * sqrt(D) * exp(-18.66*(T-1)*(T-1) - pow((D-1),4.0));
return L0*L1+L2;
}
/* -----------------------------------------------------------------
*** surten - Returns the surface tension of vapor-saturated liquid
* H2O in MPa*cm (converted from N/m) as computed from
* the Vargaftik et al. (1983) equation. See equations
* 10.1-2, Kestin et al. (1984); compare also equation
* C.5 and table 11, Haar et al. (1984). */
double surten(Tsatur)
double Tsatur;
{
double Tnorm;
if ((Tsatur < 273.16) || (Tsatur > 647.067))
return 0.0;
if (Tsatur >= (647.067-0.0000000001))
Tnorm = 0.0;
else
Tnorm = (0.999686 - Tsatur/647.27) / 0.999686;
return 0.2358*pow(Tnorm,1.256)*(1.0 - 0.625*Tnorm);
}
/* -----------------------------------------------------------------
*** Born92 - Computes the Z, Q, Y, and X Born functions at TK, Pbars.
*** epseqn = 4 ...... use Johnson-Norton (1991) equation */
void Born92(TK,Pbars,dgcm3,betab,alphaK,daldT,eps,Z,Q,Y,X,epseqn)
int epseqn;
double TK,Pbars,dgcm3,betab,alphaK,daldT,*eps,*Z,*Q,*Y,*X;
{
double TdegC,dedP,dedT,d2edT2;
extern void JN91(),epsBrn();
*eps = *Z = *Y = *Q = *X = 0.0;
TdegC = TK - 273.15;
if ((TdegC > (1000.0+0.001)) || (Pbars > (5000.0+0.001)))
return;
if (epseqn == 4)
{ JN91(TK,dgcm3,betab,alphaK,daldT,eps,&dedP,&dedT,&d2edT2);
epsBrn(*eps,dedP,dedT,d2edT2,Z,Q,Y,X);
return;
}
/* more axed code per .f */
return;
}
/* -----------------------------------------------------------------
*** JN91 - Compute (eps, dedP, dedT, d2edT2)(T,D) using equations
*** given by Johnson and Norton (1991); fit parameters
*** regressed from least squares fit to dielectric data
*** consistent with the HK74 equation and low temperatures,
*** and with the Pitz83 equation at high temperatures.
***
*** Units: T ............... K
*** D ............... g/cm**3
*** beta, dedP ...... bar**(-1)
*** alpha, dedT ..... K**(-1)
*** daldT, d2edT2 ... K**(-2) */
void JN91(T,D,beta,alpha,daldT,eps,dedP,dedT,d2edT2)
double T,D,beta,alpha,daldT,*eps,*dedP,*dedT,*d2edT2;
{ /*OK OK*/
int i;
double Tn,T2,T3,Tr2;
double c[5],dcdT[5],dc2dTT[5];
double Tref = 298.15;
double a[10] = { 14.70333593,212.8462733,-115.4445173,19.55210915, -83.30347980,32.13240048,-6.694098645,-37.86202045,68.87359646, -27.29401652 };
Tn = T / Tref;
T2 = T * T;
T3 = T2 * T;
Tr2 = Tref * Tref;
c[0] = 1.0;
dcdT[0] = dc2dTT[0] = 0.0;
c[1] = a[0]/Tn;
dcdT[1] = -a[0]*Tref/T2;
dc2dTT[1] = 2.0*a[0]*Tref/T3;
c[2] = a[1]/Tn + a[2] + a[3]*Tn;
dcdT[2] = -a[1]*Tref/T2 + a[3]/Tref;
dc2dTT[2] = 2.0*a[1]*Tref/T3;
c[3] = a[4]/Tn + a[5]*Tn + a[6]*Tn*Tn;
dcdT[3] = -a[4]*Tref/T2 + a[5]/Tref + 2.0*a[6]*T/Tr2;
dc2dTT[3] = 2.0*a[4]*Tref/T3 + 2.0*a[6]/(Tref*Tref);
c[4] = a[7]/(Tn*Tn) + a[8]/Tn + a[9];
dcdT[4] = -2.0*a[7]*Tr2/T3 - a[8]*Tref/T2;
dc2dTT[4] = 6.0*a[7]*Tr2/(T2*T2) + 2.0*a[8]*Tref/T3;
*eps = *dedP = *dedT = *d2edT2 = 0.0;
for (i=0;i<5;i++)
*eps += c[i]*pow(D,(double) i);
for (i=0;i<5;i++)
*dedP += i*c[i]*pow(D,(double) i);
*dedP *= beta;
for (i=0;i<5;i++)
*dedT += pow(D,(double) i)*(dcdT[i] - i*alpha*c[i]);
for (i=0;i<5;i++)
*d2edT2 += pow(D,(double) i)*(dc2dTT[i] - i*(alpha*dcdT[i] + c[i]*daldT) - i*alpha* (dcdT[i] - i*alpha*c[i]));
return;
}
/* -----------------------------------------------------------------
*** epsBrn - Compute the Z, Q, Y, and X Born functions from their
*** eps, dedP, dedT, and d2edT2 counterparts. */
void epsBrn(eps,dedP,dedT,d2edT2,Z,Q,Y,X)
double eps,dedP,dedT,d2edT2,*Z,*Q,*Y,*X;
{
*Z = -1.0/eps;
*Q = 1.0/(eps*eps) * dedP;
*Y = 1.0/(eps*eps) * dedT;
*X = 1.0/(eps*eps) * d2edT2 - 2.0*eps * (*Y) * (*Y);
return;
}
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.