#pragma TextEncoding = "UTF-8" #pragma rtGlobals=3 // Use modern global access method and strict wave access #pragma DefaultTab={3,20,4} // Set default tab width in Igor Pro 9 and later // PARAMETERS from NPROP.CMN: Static Constant NPROP = 43 Static Constant NRIMAX = 4 // PARAMETER from NIST.c: Static Constant DELTAT = 0.0001 // used to bias the temperaure at he saturation line // Next, the PARAMETERS from WCONST.CMN: Static Constant TCW = 647.096 Static Constant PCW = 22.064 Static Constant RHOCW = 322.0 Static Constant RW = 0.46151805 Static Constant WMMASS = 18.015268 Static Constant TTRIPW = 273.16 Static Constant PTRIPW = 611.657e-6 // PARAMETERS from WLIMIT.CMN: Static Constant TMELTX = 251.165 Static Constant TLOWER = 90.0 Static Constant TUPPER = 5000.0 Static Constant PLOWER = 1.0e-20 Static Constant PUPPER = 10000.0 Static Constant TSMIN = 253.15 Static Constant PSMIN = 125.45844e-6 // // Functions from NIST.c // // // Functions from auxpk.c // Static Function DLSAT(double TR) double E1 = 1.0/3.0; double E2 = 2.0/3.0; double E3 = 5.0/3.0; double E4 = 16.0/3.0; double E5 = 43.0/3.0; double E6 = 110.0/3.0; Make/D/FREE B = { 0.0, 1.99274064, 1.09965342, -0.510839303, -1.75493479, -45.5170352, -6.7469445e5 }; double TAU; double dlsat; TAU = 1.0 - TR; dlsat = 1.0 + B[1]*TAU^E1 + B[2]*TAU^E2 + B[3]*TAU^E3 + B[4]*TAU^E4 + B[5]*TAU^E5 + B[6]*TAU^E6; return dlsat End Static Function DVSAT(double TR) double E1 = 1.0/3.0; double E2 = 2.0/3.0; double E3 = 4.0/3.0; double E4 = 3.0; double E5 = 37.0/6.0; double E6 = 71.0/6.0; Make/D/FREE C = { 0.0, -2.03150240, -2.68302940, -5.38626492, -17.2991605, -44.7586581, -63.9201063 }; double TAU; double dvsat; TAU = 1.0 - TR; dvsat = C[1]*TAU^E1 + C[2]*TAU^E2 + C[3]*TAU^E3 + C[4]*TAU^E4 + C[5]*TAU^E5 + C[6]*TAU^E6; return exp(dvsat) End Static Function DLSATC(double TR) double A1 = 1.281714; double B = 0.326; double TAU; double dlsatc; TAU = 1.0 - TR; dlsatc = 1.0 + A1*TAU^B; return dlsatc End Static Function DVSATC(double TR) double A2 = 1.293771; double B = 0.326; double TAU; double dvsatc; TAU = 1.0 - TR; dvsatc = 1.0 - A2*TAU^B; return dvsatc End Static Function PVAP(double TR, double &PS, double &DPSDT) Make/D/FREE A = { 0.0, -7.85951783, 1.84408259, -11.7866497, 22.6807411, -15.9618719, 1.80122502 }; double TAU, SUMM, PSLN, DSDTAU, DLNPDT; TAU = 1.0 - TR; SUMM = TAU * (A[1] + TAU*TAU*(A[3] + A[5]*TAU) + sqrt(TAU)*(A[2] + TAU*TAU*(A[4] + A[6]*TAU^4))); PSLN = SUMM / TR; PS = exp(PSLN); DSDTAU = A[1] + TAU*TAU*(3.0*A[3] + 4.0*A[5]*TAU) + sqrt(TAU)*(1.50*A[2] + 3.5*A[4]*TAU*TAU + 7.5*A[6]*TAU^6); DLNPDT = -(PSLN + DSDTAU) / TR; DPSDT = PS * DLNPDT; End Static Function TVAP(double PR) double TOL = 1.0e-7; double PRLOG, PL, TG, PSNEW, DPSDT, PERRLN, TGINV, DLNPTI, TGINV2; PRLOG = ln(PR); // FOR INITIAL GUESS, USE GUESS FROM HGK ROUTINE TSAT PL = 7.1226152 + 2.302585*PRLOG; TG = 0.57602 + PL*(0.042887 + PL*(0.00368 + PL*(3.837e-4 + PL*3.0e-5))); if (TG < 0.422) TG = 0.422; endif if (TG > 0.9996) TG = 0.9996; endif do PVAP(TG, PSNEW, DPSDT); // get new iterate value for PSNEW (and hence TG). PERRLN = PRLOG - ln(PSNEW); TGINV = 1.0 / TG; DLNPTI = -TG*TG*DPSDT/PSNEW; TGINV2 = TGINV + PERRLN/DLNPTI; TG = 1.0 / TGINV2; if (TG >= 1.0) TG = 0.999999; endif while(abs((PSNEW - PR)/PR) > TOL) return TG End Static Function PMELT1(double TK, double &PMPA) double TN = 273.16; double TMIN = 251.165; double PN = 611.657e-6; int IERR; double THETA, TERM2, TERM3, PRAT; PMPA = 0.0; IERR = 0; if (TK < TMIN) IERR = 1; elseif (TK > TN) IERR = 2; else THETA = TK / TN; TERM2 = -0.626e6 * (1.0 - THETA^(-3)); TERM3 = 0.197135e6 * (1.0 - THETA^21.2); PRAT = 1.0 + TERM2 + TERM3; PMPA = PRAT * PN; endif return IERR End Static Function PMELT2(double TK, double &PMPA, int &IFORM) double TMIN3 = 251.165; double PN3 = 209.9; double TMIN5 = 256.164; double PN5 = 350.1; double TMIN6 = 273.31; double PN6 = 632.4; double TMIN7 = 355.0; double PN7 = 2216.0; double TMAX7 = 715.0; int IERR; double THETA, PRAT, TERM1, TERM2, TERM3; PMPA = 0.0; IERR = 0; IFORM = 0; if (TK < TMIN3) IERR = 1; elseif (TK > TMAX7) IERR = 2; elseif (TK <= TMIN5) THETA = TK / TMIN3; PRAT = 1.0 - 0.295252*(1.0 - THETA^60); PMPA = PRAT * PN3; IFORM = 3; elseif (TK <= TMIN6) THETA = TK / TMIN5; PRAT = 1.0 - 1.18721*(1.0 - THETA^8); PMPA = PRAT * PN5; IFORM = 5; elseif (TK <= TMIN7) THETA = TK / TMIN6; PRAT = 1.0 - 1.07476*(1.0 - THETA^4.6); PMPA = PRAT * PN6; IFORM = 6; else THETA = TK / TMIN7; TERM1 = 1.73683 * (1.0 - 1.0/THETA); TERM2 = -0.0544606 * (1.0 - THETA^5); TERM3 = 0.806106e-7 * (1.0 - THETA^22); PRAT = exp(TERM1 + TERM2 + TERM3); PMPA = PRAT * PN7; IFORM = 7; endif return IERR End // // Functions from eospk.c // Static Function AIDEAL(int LPHI, int LPHID, int LPHIT, int LPHIDD, int LPHITT, int LPHIDT, int LPHDDD, int LPHDDT, int LPHDTT, double &PHI, double &PHID, double &PHIT, double &PHIDD, double &PHITT, double &PHIDT, double &PHIDDD, double &PHIDDT, double &PHIDTT, double TAU, double DEL) Make/D/FREE WA0 = { 0.0, -8.32044648375, 6.683210527593, 3.00632, 0.012436, 0.97315, 1.27950, 0.96956, 0.24873 }; Make/D/FREE WGAM0 = { 0.0, 0.0, 0.0, 0.0, 1.28728967, 3.53734222, 7.74073708, 9.24437796, 27.5075105 }; double EGAM, DIFF; PHI = 0.0; PHID = 0.0; PHIDD = 0.0; PHIT = 0.0; PHITT = 0.0; PHIDT = 0.0; PHIDDD = 0.0; PHIDDT = 0.0; PHIDTT = 0.0; if (LPHI) PHI = ln(DEL) + WA0[1] + WA0[2]*TAU + WA0[3]*ln(TAU); endif if (LPHID) PHID = 1.0/DEL; endif if (LPHIDD) PHIDD = -1.0/(DEL*DEL); endif if (LPHIT) PHIT = WA0[2] + WA0[3]/TAU; endif if (LPHITT) PHITT = -WA0[3]/(TAU*TAU); endif if (LPHDDD) PHIDDD = 2.0 / DEL^3; endif if (LPHI || LPHIT || LPHITT) int i for (i = 4; i<9; i++) EGAM = exp(-WGAM0[i]*TAU); DIFF = 1.0 - EGAM; if (LPHI) PHI = PHI + WA0[i]*ln(DIFF); endif if (LPHIT) PHIT = PHIT + WA0[i]*WGAM0[i]*(1.0/DIFF - 1.0); endif if (LPHITT) PHITT = PHITT - WA0[i]*WGAM0[i]*WGAM0[i]*EGAM/(DIFF*DIFF); endif endfor endif End Static Function ARESID(int LPHI, int LPHID, int LPHIT, int LPHIDD, int LPHITT, int LPHIDT, int LPHDDD, int LPHDDT, int LPHDTT, double &PHI, double &PHID, double &PHIT, double &PHIDD, double &PHITT, double &PHIDT, double &PHIDDD, double &PHIDDT, double &PHIDTT, double TAU, double DEL) Make/I/FREE IWC = { 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 4, 6, 6, 6, 6 }; Make/I/FREE IWD = { 0, 1, 1, 1, 2, 2, 3, 4, 1, 1, 1, 2, 2, 3, 4, 4, 5, 7, 9, 10, 11, 13, 15, 1, 2, 2, 2, 3, 4, 4, 4, 5, 6, 6, 7, 9, 9, 9, 9, 9, 10, 10, 12, 3, 4, 4, 5, 14, 3, 6, 6, 6, 3, 3, 3 }; Make/D/FREE WT = { 0.0, -0.5, 0.875, 1.0, 0.5, 0.75, 0.375, 1.0 }; Make/I/FREE IWT = { 0, 0, 0, 0, 0, 0, 0, 0, 4, 6, 12, 1, 5, 4, 2, 13, 9, 3, 4, 11, 4, 13, 1, 7, 1, 9, 10, 10, 3, 7, 10, 10, 6, 10, 10, 1, 2, 3, 4, 8, 6, 9, 8, 16, 22, 23, 23, 10, 50, 44, 46, 50, 0, 1, 4 }; Make/D/FREE WN = { 0.0, 1.2533547935523e-02, 7.8957634722828, -8.7803203303561, 3.1802509345418e-01, -2.6145533859358e-01, -7.8199751687981e-03, 8.8089493102134e-03, -6.6856572307965e-01, 2.0433810950965e-01, -6.6212605039687e-05, -1.9232721156002e-01, -2.5709043003438e-01, 1.6074868486251e-01, -4.0092828925807e-02, 3.9343422603254e-07, -7.5941377088144e-06, 5.6250979351888e-04, -1.5608652257135e-05, 1.1537996422951e-09, 3.6582165144204e-07, -1.3251180074668e-12, -6.2639586912454e-10, -1.0793600908932e-01, 1.7611491008752e-02, 2.2132295167546e-01, -4.0247669763528e-01, 5.8083399985759e-01, 4.9969146990806e-03, -3.1358700712549e-02, -7.4315929710341e-01, 4.7807329915480e-01, 2.0527940895948e-02, -1.3636435110343e-01, 1.4180634400617e-02, 8.3326504880713e-03, -2.9052336009585e-02, 3.8615085574206e-02, -2.0393486513704e-02, -1.6554050063734e-03, 1.9955571979541e-03, 1.5870308324157e-04, -1.6388568342530e-05, 4.3613615723811e-02, 3.4994005463765e-02, -7.6788197844621e-02, 2.2446277332006e-02, -6.2689710414685e-05, -5.5711118565645e-10, -1.9905718354408e-01, 3.1777497330738e-01, -1.1841182425981e-01, -3.1306260323435e+01, 3.1546140237781e+01, -2.5213154341695e+03, -1.4874640856724e-01, 3.1806110878444e-01 }; Make/I/FREE/N=55 IWALF = 0 IWALF[52] = 20 IWALF[53] = 20 IWALF[54] = 20 Make/I/FREE/N=55 IWBET = 0 IWBET[52] = 150 IWBET[53] = 150 IWBET[54] = 250 Make/D/FREE/N=55 WGAM = 0.0 WGAM[52] = 1.21 WGAM[53] = 1.21 WGAM[54] = 1.25 Make/I/FREE/N=55 IWEPS = 0 IWEPS[52] = 1 IWEPS[53] = 1 IWEPS[54] = 1 Make/D/FREE/N=57 WA = 0.0 WA[55] = 3.5 WA[56] = 3.5 Make/D/FREE/N=57 WB = 0.0 WB[55] = 0.85 WB[56] = 0.95 Make/D/FREE/N=57 WBB = 0.0 WBB[55] = 0.2 WBB[56] = 0.2 Make/I/FREE/N=57 IWCC = 0 IWCC[55] = 28 IWCC[56] = 32 Make/I/FREE/N=57 IWDD = 0 IWDD[55] = 700 IWDD[56] = 800 Make/D/FREE/N=57 WAA = 0.0 WAA[55] = 0.32 WAA[56] = 0.32 Make/D/FREE/N=57 WBET = 0.0 WBET[55] = 0.3 WBET[56] = 0.3 int IC2, IC3, IS1, IS2, IS3, IS4, IA2E, IAE2; double DINV, DINV2, DINV3, TINV, TINV2, DELPOW, TAUPOW, WNDT, DELC, WNDTE, PROD1; double DEPS, TGAM, ALFEXP, ALFBET, BETEXP, Q, DM1, TM1, PEXP, PSI, THETA, DELTA; double BM1, BM2, BM3, AM1, AM2, AM3, DB, DBM1, DBM2, DBM3, DDDT, DDBDT, DPDT; double DQDD, OVERB, OBM1, DTHDD, DDDD, DDBDD, DPDD, D2DBTD, D2PDTD, OBM2; double D2THDD, D2DDD, D2DBDD, D2PDD, D2DBDT, D2PDT, D3THDD, D3DDD, D3DBDD, D3PDD; double D3DBDT, D3PDDT, D3DBTD, D3PDTT, EPOW, S2; PHI = 0.0; PHID = 0.0; PHIDD = 0.0; PHIT = 0.0; PHITT = 0.0; PHIDT = 0.0; PHIDDD = 0.0; PHIDDT = 0.0; PHIDTT = 0.0; DINV = 1.0 / DEL; DINV2 = DINV*DINV; DINV3 = DINV2*DINV; TINV = 1.0 / TAU; TINV2 = TINV*TINV; int i for (i=1; i<8; i++) DELPOW = DEL^IWD[i]; TAUPOW = TAU^WT[i]; WNDT = WN[i]*DELPOW*TAUPOW; if (LPHI) PHI = PHI + WNDT; endif if (LPHID) PHID = PHID + WNDT*IWD[i]*DINV; endif if (LPHIDD) PHIDD = PHIDD + WNDT*IWD[i]*(IWD[i]-1)*DINV2; endif if (LPHIT) PHIT = PHIT + WNDT*WT[i]*TINV; endif if (LPHITT) PHITT = PHITT + WNDT*WT[i]*(WT[i]-1.0)*TINV2; endif if (LPHIDT) PHIDT = PHIDT + WNDT*IWD[i]*DINV*WT[i]*TINV; endif if (LPHDDD) PHIDDD = PHIDDD + WNDT*IWD[i]*(IWD[i]-1)*(IWD[i]-2)*DINV3; endif if (LPHDDT) PHIDDT = PHIDDT + WNDT*IWD[i]*(IWD[i]-1)*WT[i]*TINV*DINV2; endif if (LPHDTT) PHIDTT = PHIDTT + WNDT*IWD[i]*WT[i]*(WT[i]-1.0)*TINV2*DINV; endif endfor for (i=8; i<52; i++) DELPOW = DEL^IWD[i]; TAUPOW = TAU^IWT[i]; DELC = DEL^IWC[i]; if (DELC > 700.) EPOW = 0.0; else EPOW = exp(-DELC); endif WNDTE = WN[i]*DELPOW*TAUPOW*EPOW; if (LPHI) PHI = PHI + WNDTE; endif if (LPHID) PHID = PHID + WNDTE*DINV*(IWD[i]-IWC[i]*DELC); endif if (LPHIDD) PROD1 = IWC[i]*DELC; PHIDD = PHIDD + WNDTE*DINV2*((IWD[i]-PROD1)*(IWD[i]-1.0-PROD1) - IWC[i]*PROD1); endif if (LPHIT) PHIT = PHIT + WNDTE*IWT[i]*TINV; endif if (LPHITT) PHITT = PHITT + WNDTE*IWT[i]*(IWT[i]-1)*TINV2; endif if (LPHIDT) PHIDT = PHIDT + WNDTE*DINV*IWT[i]*TINV*(IWD[i] - IWC[i]*DELC); endif if (LPHDDD) IC2 = IWC[i]*IWC[i]; IC3 = IC2*IWC[i]; IS1 = -2*IWC[i] + 3*IC2 - IC3 + 6*IWC[i]*IWD[i] - 3*IWD[i]*(IC2 + IWC[i]*IWD[i]); IS2 = 3 * (IC3 - IC2 + IC2*IWD[i]); PHIDDD = PHIDDD + WNDTE*DINV3*(IWD[i]*(2-3*IWD[i]+IWD[i]*IWD[i]) + IS1*DELC+ IS2*DELC*DELC - IC3*DELC^3); endif if (LPHDDT) IC2 = IWC[i]*IWC[i]; IS1 = IWC[i] - IC2 - 2*IWC[i]*IWD[i]; PHIDDT = PHIDDT + WNDTE*IWT[i]*TINV*DINV2*(IWD[i]*IWD[i] - IWD[i] + IS1*DELC + IC2*DELC*DELC); endif if (LPHDTT) PHIDTT = PHIDTT + WNDTE*IWT[i]*(IWT[i]-1)*TINV2*DINV*(IWD[i] - IWC[i]*DELC); endif endfor for (i=52; i<55; i++) DELPOW = DEL^IWD[i]; TAUPOW = TAU^IWT[i]; DEPS = DEL - IWEPS[i]; TGAM = TAU - WGAM[i]; ALFEXP = IWALF[i]*DEPS*DEPS; BETEXP = IWBET[i]*TGAM*TGAM; ALFBET = ALFEXP + BETEXP; if (ALFBET > 700.0) EPOW = 0.0; else EPOW = exp(-ALFBET); endif WNDTE = WN[i]*DELPOW*TAUPOW*EPOW; if (LPHI) PHI = PHI + WNDTE; endif if (LPHID) PHID = PHID + WNDTE*(IWD[i]*DINV - 2.0*IWALF[i]*DEPS); endif if (LPHIDD) PHIDD = PHIDD + WNDTE*(-2.0*IWALF[i] + 4.0*IWALF[i]*IWALF[i]*DEPS*DEPS - 4.0*IWD[i]*IWALF[i]*DINV*DEPS + IWD[i]*(IWD[i] - 1)*DINV2); endif if (LPHIT) PHIT = PHIT + WNDTE * (IWT[i]*TINV - 2.0*IWBET[i]*TGAM); endif if (LPHITT) PHITT = PHITT + WNDTE*((IWT[i]*TINV - 2.0*IWBET[i]*TGAM)^2 - IWT[i]*TINV2 - 2.0*IWBET[i]); endif if (LPHIDT) PHIDT = PHIDT + WNDTE*(IWD[i]*DINV - 2.0*IWALF[i]*DEPS)*(IWT[i]*TINV - 2.0*IWBET[i]*TGAM); endif if (LPHDDD) IA2E = IWALF[i]*IWALF[i]*IWEPS[i]; IAE2 = IWALF[i]*IWEPS[i]*IWEPS[i]; IS1 = 6*IWALF[i]*IWD[i]*IWEPS[i]*(IWD[i] - 1); IS2 = 6*IWALF[i]*IWD[i]*(2*IAE2 - IWD[i]); IS3 = 4*IA2E*(2*IAE2 - 6*IWD[i] - 3); IS4 = 12*IWALF[i]*IWALF[i]*(1 + IWD[i] - 2*IAE2); PHIDDD = PHIDDD + WNDTE*DINV3*(IWD[i]*(2 - 3*IWD[i] + IWD[i]*IWD[i]) + DEL*(IS1 \ + DEL*(IS2 + DEL*(IS3 + DEL*(IS4 + DEL*(24*IWALF[i]*IA2E - 8*DEL*IWALF[i]^3)))))); endif if (LPHDDT) IS2 = 2*IWALF[i]*(2*IWALF[i]*IWEPS[i]*IWEPS[i] - 2*IWD[i] - 1); PHIDDT = PHIDDT + WNDTE*DINV2*TINV*(IWD[i]*(IWD[i] - 1) + DEL*(4*IWALF[i]*IWD[i]*IWEPS[i] \ + DEL*(IS2 + DEL*IWALF[i]*IWALF[i]*(-8*IWEPS[i] + 4*DEL)))) * (IWT[i] + 2*IWBET[i]*TAU*(WGAM[i] - TAU)); endif if (LPHDTT) S2 = 2*IWBET[i]*(2*IWBET[i]*WGAM[i]*WGAM[i] - 2*IWT[i] - 1); PHIDTT = PHIDTT + WNDTE*DINV*TINV2*(IWT[i]*(IWT[i] - 1) \ + TAU*(4*IWBET[i]*WGAM[i]*IWT[i] + TAU*(S2 + TAU*(-8*IWBET[i]*IWBET[i]*WGAM[i] \ + 4*IWBET[i]*IWBET[i]*TAU))))*(IWD[i] + 2*IWALF[i]*DEL*(IWEPS[i] - DEL)); endif endfor Q = (1.0 - DEL)^2; DM1 = DEL - 1.0; TM1 = TAU - 1.0; for (i=55; i<57; i++) PEXP = -IWCC[i]*Q - IWDD[i]*TM1*TM1; if (PEXP < -700.0) PSI = 0.0; else PSI = exp(PEXP); endif THETA = 1.0 - TAU + WAA[i]*Q^(0.50/WBET[i]); DELTA = THETA*THETA + WBB[i]*Q^WA[i]; BM1 = WB[i] - 1.0; BM2 = BM1 - 1.0; BM3 = BM2 - 1.0; AM1 = WA[i] - 1.0; AM2 = AM1 - 1.0; AM3 = AM2 - 1.0; DB = DELTA^WB[i]; if (DELTA == 0.0) DBM1 = 0.0; DBM2 = 0.0; DBM3 = 0.0; else DBM1 = DELTA^BM1; DBM2 = DELTA^BM2; DBM3 = DELTA^BM3; endif if (LPHI) PHI = PHI + WN[i]*DB*DEL*PSI; endif if (LPHIT || LPHITT || LPHIDT || LPHDDT || LPHDTT) DDDT = -2.0*THETA; DDBDT = WB[i]*DBM1*DDDT; DPDT = -2.0*IWDD[i]*TM1*PSI; endif if (LPHID || LPHIDD || LPHIDT || LPHDDD || LPHDDT || LPHDTT) DQDD = 2.0*DM1; OVERB = 0.50/WBET[i]; OBM1 = OVERB - 1.0; DTHDD = WAA[i]*OVERB * (Q^OBM1)*DQDD; DDDD = 2.0*THETA*DTHDD + WBB[i]*WA[i]*(Q^AM1)*DQDD; DDBDD = WB[i]*DBM1*DDDD; DPDD = -2.0*IWCC[i]*DM1*PSI; endif if (LPHIDT || LPHDDT || LPHDTT) D2DBTD = -2.0*WB[i]*(THETA*BM1*DBM2*DDDD + DTHDD*DBM1); D2PDTD = 4.0*IWCC[i]*IWDD[i]*DM1*TM1*PSI; endif if (LPHIDD || LPHDDD || LPHDDT) OBM2 = OBM1 - 1.0; if (Q == 0.0) D2THDD = 0.0; D2DDD = 2.0*(THETA*D2THDD + DTHDD*DTHDD); else D2THDD = WAA[i]*OVERB*(OBM1*(Q^OBM2)*DQDD*DQDD + 2.0*Q^OBM1); D2DDD = 2.0*(THETA*D2THDD + DTHDD*DTHDD) + WBB[i]*WA[i]*(AM1*(Q^AM2)*DQDD*DQDD + 2.0*Q^AM1); endif D2DBDD = WB[i]*(BM1*DBM2*DDDD*DDDD + DBM1*D2DDD); D2PDD = 2.0*IWCC[i]*PSI*(2.0*IWCC[i]*Q - 1.0); endif if (LPHITT || LPHDTT) D2DBDT = 2.0*WB[i]*(2.0*THETA*THETA*BM1*DBM2 + DBM1); D2PDT = 2.0*IWDD[i]*PSI*(2.0*IWDD[i]*TM1*TM1 - 1.0); endif if (LPHID) PHID = PHID + WN[i]*(DB*(PSI+DEL*DPDD) + DDBDD*DEL*PSI); endif if (LPHIDD) PHIDD = PHIDD + WN[i]*(DB*(2.0*DPDD + DEL*D2PDD) + 2.0*DDBDD*(PSI + DEL*DPDD) + D2DBDD*DEL*PSI); endif if (LPHIT) PHIT = PHIT + WN[i]*DEL*(DDBDT*PSI + DB*DPDT); endif if (LPHITT) PHITT = PHITT + WN[i]*DEL*(D2DBDT*PSI + 2.0*DDBDT*DPDT + DB*D2PDT); endif if (LPHIDT) PHIDT = PHIDT + WN[i]*(DB*(DPDT + DEL*D2PDTD) + DEL*DDBDD*DPDT + DDBDT*(PSI + DEL*DPDD) + D2DBTD*DEL*PSI); endif if (LPHDDD) if (Q == 0.0) D3THDD = 0.0; D3DDD = 0.0; else D3THDD = WAA[i]*OVERB*OBM1 * (OBM2*(Q^(OVERB-3.0))*DQDD^3 + 6.0*(Q^OBM2)*DQDD); D3DDD = 6.0*DTHDD*D2THDD + 2.0*THETA*D3THDD + WBB[i]*WA[i]*AM1 \ * (AM2*(Q^AM3)*(DQDD^3) + 6.0*(Q^AM2)*DQDD); endif D3DBDD = WB[i] * (BM1*(BM2*DBM3*DDDD^3 + 3.0*DBM2*DDDD*D2DDD) + DBM1*D3DDD); D3PDD = 4*IWCC[i]*IWCC[i]*DM1*PSI * (3 - 2*IWCC[i]*DM1*DM1); PHIDDD = PHIDDD + WN[i]*(DB*(3.0*D2PDD + DEL*D3PDD) + 3.0*DDBDD*(2.0*DPDD + DEL*D2PDD) \ + 3.0*D2DBDD*(PSI + DEL*DPDD) + D3DBDD*DEL*PSI); endif if (LPHDDT) D3DBDT = -2.0*WB[i] * (BM1*(BM2*THETA*DBM3*DDDD*DDDD + 2.0*DBM2*DDDD*DTHDD \ + THETA*DBM2*D2DDD) + DBM1*D2THDD); D3PDDT = -4*IWCC[i]*IWDD[i]*TM1*PSI * (2*IWCC[i]*DM1*DM1 - 1.0); PHIDDT = PHIDDT + WN[i]*(DB*(2.0*D2PDTD + DEL*D3PDDT) + DDBDT*(2.0*DPDD + DEL*D2PDD) \ + 2.0*(DDBDD*(DPDT + DEL*D2PDTD) + D2DBTD*(PSI + DEL*DPDD)) \ + D2DBDD*DEL*DPDT + D3DBDT*DEL*PSI); endif if (LPHDTT) D3DBTD = 2.0*WB[i]*BM1 * (4.0*THETA*DTHDD*DBM2 + 2.0*THETA*THETA*BM2*DBM3*DDDD + DBM2*DDDD); D3PDTT = -4*IWCC[i]*IWDD[i]*DM1*PSI * (2*IWDD[i]*TM1*TM1 -1.0); PHIDTT = PHIDTT + WN[i]*(DB*(D2PDT + DEL*D3PDTT) + 2.0*DDBDT*(DPDT + DEL*D2PDTD) \ + DDBDD*DEL*D2PDT + D2DBDT*(PSI + DEL*DPDD) + 2.0*D2DBTD*DEL*DPDT + D3DBTD*DEL*PSI); endif endfor End // // Functions from intpk.c // Static Function DENS(double &DOUT, double PMPA, double D, double TK, double &DPD, wave IWORK, wave PROPR) int IERR, IERS, IFORM; double TSOLID, DR, PR, TR, DOUTR, DPDR; if ((PMPA < PLOWER) || (PMPA > PUPPER)) IERR = 5; DOUT = D; DPD = 0.0; return IERR; endif if ((TK < TLOWER) || (TK > TUPPER)) IERR = 6; DOUT = D; DPD = 0.0; return IERR; endif if (PMPA <= PTRIPW) IERS = TSUB(TSOLID, PMPA); if ((IERS != 1) && (TK < (TSOLID - 35.0))) IERR = 7; DOUT = 1.e-20; DPD = 0.0; return IERR; endif else IERS = TMELT(TSOLID, PMPA, IFORM); if ((IERS != 2) && (TK < (TSOLID - 35.0))) IERR = 7; DOUT = 1.0e-20; DPD = 0.0; return IERR; endif endif DR = D / RHOCW; PR = PMPA / PCW; TR = TK / TCW; IERR = DFIND(DOUTR, PR, DR, TR, &DPDR, IWORK, PROPR); DOUT = DOUTR * RHOCW; DPD = DPDR * PCW / RHOCW; return IERR End Static Function DENS0(double &DOUT, double PMPA, double TK, double &DPD, wave IWORK, wave PROPR) int IERR, IERS, IFORM; double TSOLID, TR, PR, DOUTR, DPDR; if ((PMPA < PLOWER) || (PMPA > PUPPER)) IERR = 5; DOUT = 1.0e-20; DPD = 0.0; return IERR; endif if ((TK < TLOWER) || (TK > TUPPER)) IERR = 6; DOUT = 1.0e-20; DPD = 0.0; return IERR; endif if (PMPA <= PTRIPW) IERS = TSUB(TSOLID, PMPA); if ((IERS != 1) && (TK < (TSOLID - 35.0))) IERR = 7; DOUT = 1.0e-20; DPD = 0.0; return IERR; else IERS = TMELT(&TSOLID, PMPA, &IFORM); if ((IERS != 2) && (TK < (TSOLID - 35.0))) IERR = 7; DOUT = 1.0e-20; DPD = 0.0; return IERR; endif endif endif TR = TK / TCW; PR = PMPA / PCW; IERR = DFIND0(DOUTR, PR, TR, DPDR, IWORK, PROPR); DOUT = DOUTR * RHOCW; DPD = DPDR * PCW / RHOCW; return IERR End Static Function PROPS(wave IWANT, double TK, double RHO, wave PROPSI, wave PROPR, int I2PHCK, int &I2PH, int ISCHK, int &ISFLG, int ICCHK, int &ICFLG, wave IPCHK, wave IPFLG, int IGFLG, int NRI, wave WAVRI, wave RI, wave IRIFLG) Make/D/FREE/N=(NPROP + 1) TPROP; Make/I/FREE/N=(NPROP + 1) IWTMP; Make/I/FREE/N=6 IPDUM; int ICLIQ, ICVAP, IERRPS, IPMAX; double Q, RHOV, RHOL, Q1, VVAP, VLIQ, DPSDT, DVVDT, DVLDT, DUVDP, DULDP, DUVDT, DULDT; double DUVDTS, DULDTS, DQDT, DELU, CV2PH, PDUM, PMPA; I2PH = 0; if ((I2PHCK != 0) && (IGFLG != 1)) if (TK < TCW) IERRPS = PSAT(TK, PMPA, RHOL, RHOV, IWTMP, PROPR); if (IERRPS == 1) I2PH = 3; else if ((RHO < RHOL) && (RHO > RHOV)) if (IERRPS == 4) I2PH = 4; else I2PH = 2; endif Q = (1.0 / RHO - 1.0 / RHOL) / (1.0 / RHOV - 1.0 / RHOL); int i for (i = 1; i <= NPROP; i++) IWTMP[i] = IWANT[i]; endfor if (IWANT[8] != 0) IWTMP[2] = 1; IWTMP[6] = 1; IWTMP[9] = 1; IWTMP[10] = 1; IWTMP[14] = 1; IWTMP[15] = 1; endif PROPS1(IWTMP, IGFLG, TK, RHOV, TPROP, PROPR); PROPS1(IWTMP, IGFLG, TK, RHOL, PROPSI, PROPR); Q1 = 1.0 - Q; if (IWANT[8] != 0) VVAP = 1.0 / RHOV; VLIQ = 1.0 / RHOL; DPSDT = (TPROP[6] - PROPSI[6]) / (TK * (VVAP - VLIQ)); DPSDT = DPSDT / 1.0e3; DVVDT = VVAP * (TPROP[15] - TPROP[14] * DPSDT); DVLDT = VLIQ * (PROPSI[15] - PROPSI[14] * DPSDT); DUVDP = 1.0e3 * VVAP * (TPROP[14] * TPROP[2] - TPROP[15] * TK); DULDP = 1.0e3 * VLIQ * (PROPSI[14] * PROPSI[2] - PROPSI[15] * TK); DUVDT = TPROP[9] - 1.0e3 * TPROP[2] * VVAP * TPROP[15]; DULDT = PROPSI[9] - 1.0e3 * PROPSI[2] * VLIQ * PROPSI[15]; DUVDTS = DUVDT + DUVDP * DPSDT; DULDTS = DULDT + DULDP * DPSDT; DQDT = -(Q * DVVDT + Q1 * DVLDT) / (VVAP - VLIQ); DELU = TPROP[10] - PROPSI[10]; CV2PH = Q * DUVDTS + Q1 * DULDTS + DQDT * DELU; PROPSI[8] = CV2PH; endif if (IWANT[3] != 0) PROPSI[3] = RHO; if (IWANT[4] != 0) PROPSI[4] = 1.0 / RHO; endif if (IWANT[6] != 0) PROPSI[6] = Q * TPROP[6] + Q1 * PROPSI[6]; endif if (IWANT[7] != 0) PROPSI[7] = Q * TPROP[7] + Q1 * PROPSI[7]; endif if (IWANT[10] != 0) PROPSI[10] = Q * TPROP[10] + Q1 * PROPSI[10]; endif if (IWANT[11] != 0) PROPSI[11] = Q * TPROP[11] + Q1 * PROPSI[11]; endif PROPSI[9] = 0.0; for (i = 14; i <= 43; i++) PROPSI[i] = 0.0; endfor if (NRI > 0) int i for (i = 1; i <= NRI; i++) RI[i] = 0.0; endfor endif if (I2PH == 4) BNDCHK(TK, PDUM, RHOV, 2, 0, ISFLG, 0, ICFLG, IPCHK, IPFLG, IWTMP, TPROP, PROPR); BNDCHK(TK, PDUM, RHOL, 2, 0, ISFLG, 0, ICFLG, IPCHK, IPDUM, IWTMP, TPROP, PROPR); int i for (i = 1; i <= 5; i++) IPFLG[i] = (IPFLG[i] > IPDUM[i]) ? IPFLG[i] : IPDUM[i]; endfor else if (ICCHK == 1) int i for (i = 1; i <= 5; i++) IPDUM[i] = 0; endfor BNDCHK(TK, PDUM, RHOV, 2, 0, ISFLG, ICCHK, &ICVAP, IPDUM, IPFLG, IWTMP, TPROP, PROPR); BNDCHK(TK, PDUM, RHOL, 2, 0, ISFLG, ICCHK, &ICLIQ, IPDUM, IPFLG, IWTMP, TPROP, PROPR); ICFLG = (ICVAP > ICLIQ) ? ICVAP : ICLIQ; endif if ((ISCHK == 1) && (I2PH == 4)) ISFLG = 1; endif endif else I2PH = 1; endif endif if ((I2PH != 2) && (I2PH != 4)) Q = 0.0; BNDCHK(TK, PDUM, RHO, 2, ISCHK, ISFLG, ICCHK, ICFLG, IPCHK, IPFLG, IWTMP, TPROP, PROPR); IPMAX = 0; int i for (i = 1; i <= 5; i++) IPMAX = (IPMAX > IPFLG[i]) ? IPMAX : IPFLG[i]; endfor IPMAX = (IPMAX > ISFLG) ? IPMAX : *ISFLG; if (IPMAX < 2) PROPS1(IWANT, IGFLG, TK, RHO, PROPSI, PROPR); if (NRI > 0) RIND(TK, RHO, WAVRI, RI, IRIFLG, NRI); endif endif endif if ((I2PHCK != 0) && (I2PH == 0)) if ((PROPR[2] < 1.0) && (RHO <= RHOV)) I2PH = 1; else I2PH = -1; endif endif PROPSI[5] = Q; End Static Function BNDCHK(double TK, double PMPA, double RHO, int MODE, int ISCHK, int &ISFLG, int ICCHK, int &ICFLG, wave IPCHK, wave IPFLG, wave IWORK, wave PROPSI, wave RWORK) double RCMNS = 320.39; double RCPLS = 323.61; double TCMNS = 647.076; double TCPLS = 652.096; int IHAVEP, IHAVER, INEEDP, INEEDR, IERRD, IERR, IT, IP, IFORM; double TC, DPD, TSOLID; ISFLG = 0; ICFLG = 0; int i for (i = 1; i <= 5; i++) IPFLG[i] = 0; endfor TC = TK - 273.15; if (MODE == 0) IHAVEP = 1; IHAVER = 1; else if (MODE == 1) IHAVEP = 1; IHAVER = 0; else IHAVEP = 0; IHAVER = 1; endif if ((IHAVEP == 0) && ((ISCHK == 1) || (IPCHK[1] == 1) || (IPCHK[2] == 1) || (IPCHK[3] == 1) || (IPCHK[4] == 1) || (IPCHK[5] == 1))) INEEDP = 1; else INEEDP = 0; endif if ((IHAVER == 0) && ((ICCHK == 1) || (IPCHK[5] == 1))) INEEDR = 1; else INEEDR = 0; endif if (TK < 0.0) if (ISCHK == 1) ISFLG = 2; endif for (int i = 1; i <= 5; i++) if (IPCHK[i] == 1) IPFLG[i] = 2; endif endfor return; endif if (IHAVER == 1) if (RHO <= 0.0) for (int i = 1; i <= 5; i++) if (IPCHK[i] == 1) IPFLG[i] = 2; endif endfor return; endif endif if (INEEDP == 1) IVZERO(IWORK, NPROP); IWORK[2] = 1; PROPS1(IWORK, 0, TK, RHO, PROPSI, RWORK); PMPA = PROPSI[2]; endif if (INEEDR == 1) IERRD = DENS0(RHO, PMPA, TK, DPD, IWORK, RWORK); endif if (ISCHK == 1) if (PMPA <= PTRIPW) IERR = TSUB(TSOLID, PMPA); if (IERR != 1) if (TK < (TSOLID - 35.0)) ISFLG = 2; elseif (TK < TSOLID) ISFLG = 1; endif else IERR = TMELT(&TSOLID, PMPA, IFORM); if (IERR != 2) if (TK < (TSOLID - 35.0)) ISFLG = 2; elseif (TK < TSOLID) ISFLG = 1; endif endif endif endif endif if ((ICCHK == 1) && ((IERRD == 0) || (abs(IERRD) == 3))) if ((TK > TCMNS) && (TK < TCPLS) && (RHO > RCMNS) && (RHO < RCPLS)) ICFLG = 1; endif endif if (IPCHK[1] == 1) IT = 0; IP = 0; if ((TK > TUPPER) || (TK < TLOWER)) IT = 2; elseif ((TK > 1273.15) || (TK < 273.15)) IT = 1; endif if ((PMPA > PUPPER) || (PMPA < PLOWER)) IP = 2; elseif ((PMPA > 1.0e3) || (PMPA <= 0.0)) IP = 1; endif IPFLG[1] = (IP > IT) ? IP : IT; endif if (IPCHK[2] == 1) IT = 0; IP = 0; if ((TK > TUPPER) || (TK < TLOWER)) IT = 2; else if ((TC > 900.0) || ((TC > 600.0) && (PMPA > 300.0)) || ((TC > 150.0) && (PMPA > 350.0)) || (PMPA > 500.0) || (TC < 0.0)) IT = 1; endif if ((PMPA > PUPPER) || (PMPA < PLOWER)) IP = 2; else if (PMPA <= 0.0) IP = 1; endif IPFLG[2] = (IP > IT) ? IP : IT; endif if (IPCHK[3] == 1) IT = 0; IP = 0; if ((TK > TUPPER) || (TK < TLOWER)) IT = 2; else if ((TC > 800.0) || ((TC > 400.0) && (PMPA > 100.0)) || ((TC > 250.0) && (PMPA > 150.0)) || ((TC > 125.0) && (PMPA > 200.0)) || (PMPA > 400.0) || (TC < 0.0)) IT = 1; endif if ((PMPA > PUPPER) || (PMPA < PLOWER)) IP = 2; else if (PMPA <= 0.0) IP = 1; endif IPFLG[3] = (IP > IT) ? IP : IT; endif if (IPCHK[4] == 1) IT = 0; IP = 0; if ((TK > TUPPER) || (TK < TLOWER)) IT = 2; else if ((TK > 873.0) || (TK < 238.0)) IT = 1; endif if ((PMPA > PUPPER) || (PMPA < PLOWER)) IP = 2; else if ((PMPA > 1.0e3) || (PMPA <= 0.0)) IP = 1; endif IPFLG[4] = (IP > IT) ? IP : IT; endif if (IPCHK[5] == 1) IT = 0; IP = 0; if ((TK > TUPPER) || (TK < TLOWER)) IT = 2; else if ((TK > 773.15) || (TK < 261.15)) IT = 1; endif if ((PMPA > PUPPER) || (PMPA < PLOWER)) IP = 2; else if (RHO > 1060.0) IP = 1; endif IPFLG[5] = (IP > IT) ? IP : IT; endif End Static Function PSAT(double TK, double &PMPA, double &RHOL, double &RHOV, wave IWORK, wave PROPR) int IERR; double TR, PR, RHOLR, RHOVR; TR = TK / TCW; IERR = PCOEX(TR, PR, RHOLR, RHOVR, IWORK, PROPR); PMPA = PR * PCW; RHOV = RHOVR * RHOCW; RHOL = RHOLR * RHOCW; return IERR End Static Function TSAT(double PMPA, double &TK, double &RHOL, double &RHOV, wave IWORK, wave PROPR) int IERR; double PR, TR, RHOLR, RHOVR; PR = PMPA / PCW; IERR = TCOEX(PR, TR, RHOLR, RHOVR, IWORK, PROPR); TK = TR * TCW; RHOL = RHOLR * RHOCW; RHOV = RHOVR * RHOCW; return IERR End Static Function TSUB(double &TK, double PMPA) double TMAX = 273.16; double TMIN = 190.0; double TOL = 1.0e-8; int IERR; double PMIN, PMAX, FHIGH, FLOW, TAUHI, TAULOW, FRDIST, TNEW, TAUNEW, FNEW, PNEW; IERR = PSUB(TMIN, &PMIN); IERR = PSUB(TMAX, &PMAX); if (PMPA < PMIN) TK = TMIN; IERR = 1; else if (PMPA > PMAX) TK = TMAX; IERR = 2; else FHIGH = ln(PMAX / PMPA); FLOW = ln(PMIN / PMPA); TAULOW = 1.0 / TMIN; TAUHI = 1.0 / TMAX; int i for (i = 1; i <= 100; i++) FRDIST = -FLOW / (FHIGH - FLOW); if (FRDIST < 0.01) FRDIST = 0.01; endif if (FRDIST > 0.99) FRDIST = 0.99; endif TAUNEW = TAULOW + FRDIST * (TAUHI - TAULOW); TNEW = 1.0 / TAUNEW; IERR = PSUB(TNEW, &PNEW); FNEW = ln(PNEW / PMPA); if (fabs(FNEW) < TOL) IERR = 0; TK = TNEW; return IERR; endif if (FNEW > 0.0) FHIGH = FNEW; TAUHI = TAUNEW; else FLOW = FNEW; TAULOW = TAUNEW; endif endfor IERR = 3; TK = TNEW; return IERR; endif return IERR End Static Function PSUB(double TK, double &PMPA) double TN = 273.16; double TMIN = 190.0; double PN = 611.657e-6; double A1 = -13.928169; double A2 = 34.7078238; int IERR; double THETA, FSUB; IERR = 0; if (TK < TMIN) IERR = 1; else if (TK > TN) IERR = 2; endif THETA = TK / TN; FSUB = A1 * (1.0 - THETA^(-1.5)) + A2 * (1.0 - THETA^(-1.25)); PMPA = exp(FSUB) * PN; return IERR End Static Function TMELT(double &TK, double PMPA, int &IFORM) double PMAX = 20618.0; double PMIN = 611.657e-6; double TTOP = 715.0; double TBOT = 273.16; double T13 = 251.165; double P13 = 209.9; double TOL = 1.0e-8; int LHIGHP; int IERR, IDUM; double THIGH, TLOW, PTOP, FHIGH, FLOW, PLOW, FRDIST, TNEW, FNEW, PNEW; i f (PMPA < PMIN) TK = TBOT; IERR = 1; else if (PMPA > PMAX) TK = TTOP; IERR = 2; else if (PMPA > P13) LHIGHP = 1; THIGH = TTOP; TLOW = T13; IERR = PMELT2(TTOP, PTOP, IDUM); FHIGH = ln(PTOP / PMPA); FLOW = ln(P13 / PMPA); else LHIGHP = 0; IFORM = 1; THIGH = TBOT; TLOW = T13; IERR = PMELT1(TLOW, PLOW); FHIGH = -ln(PMIN / PMPA); FLOW = -ln(PLOW / PMPA); endif int i for (i = 1; i <= 100; i++) FRDIST = -FLOW / (FHIGH - FLOW); if (FRDIST < 0.01) FRDIST = 0.01; endif if (FRDIST > 0.99) FRDIST = 0.99; endif TNEW = TLOW + FRDIST * (THIGH - TLOW); if (LHIGHP) IERR = PMELT2(TNEW, PNEW, IFORM); FNEW = ln(PNEW / PMPA); else IERR = PMELT1(TNEW, PNEW); FNEW = -ln(PNEW / PMPA); endif if (fabs(FNEW) < TOL) IERR = 0; TK = TNEW; return IERR; endif if (FNEW > 0.0) FHIGH = FNEW; THIGH = TNEW; else FLOW = FNEW; TLOW = TNEW; endif endfor IERR = 3; TK = TNEW; return IERR; endif return IERR End Static Function PMELT(double TK, int &NROOTS, double &PMPA1, double &PMPA2, int &IFORM) double TMAX = 715.0; double TMIN = 251.165; double T1MAX = 273.16; int IERR, IDUM; PMPA1 = 0.0; PMPA2 = 0.0; IFORM = 0; NROOTS = 0; IERR = 0; if (TK < TMIN) IERR = 1; else if (TK > TMAX) IERR = 2; else if (TK <= T1MAX) IDUM = PMELT1(TK, PMPA1); IDUM = PMELT2(TK, PMPA2, IFORM); NROOTS = 2; else IDUM = PMELT2(TK, PMPA1, IFORM); NROOTS = 1; endif return IERR End Static Function SURF(double TK, double &SURFT) int IERR; double TAU; if (TK < TSMIN) IERR = 1; SURFT = 0.0; else if (TK > TCW) IERR = 2; SURFT = 0.0; else if (TK < TTRIPW) IERR = 3; else IERR = 0; TAU = 1.0 - TK / TCW; SURFT = 0.2358 * TAU^1.256 * (1.0 - 0.625 * TAU); endif return IERR End // // Functions from proppk.c // Static Function KFACT(double TR, double RHOL, double RHOV, double &FLIQ, double &FVAP, double &XK, wave IWANT, wave PROPR) IWANT = 0 IWANT[13] = 1; PROPS2(IWANT, 0, TR, RHOL, PROPR); FLIQ = PROPR[13]; PROPS2(IWANT, 0, TR, RHOV, PROPR); FVAP = PROPR[13]; XK = FLIQ / FVAP; End // // Functions from solvpk.c // Static Function DFIND(double &DOUT, double P, double D, double TR, double &DPD, wave IWANT, wave PROPR) int IERR; if (TR >= 1.0) IERR = DFIND1(DOUT, P, D, TR, DPD, IWANT, PROPR); else IERR = DFIND2(DOUT, P, D, TR, DPD, IWANT, PROPR); return IERR End Static Function DFIND0(double &DOUT, double PR, double TR, double &DPD, wave IWANT, wave PROPR) int IERR; double PBOTR, P0AUX, DP0, P0, RHOL, RHOV, RGUESS; PBOTR = PSMIN / PCW; if (TR < 1.0) { PVAP(TR, &P0AUX, &DP0); if (fabs(1.0 - P0AUX/PR) < 0.05) IERR = PCOEX(TR, &P0, &RHOL, &RHOV, IWANT, PROPR); if (IERR == 1) P0 = P0AUX; RHOL = DLSAT(TR); if (PR > 0.1) RHOV = DVSAT(TR); endif endif else P0 = P0AUX; RHOL = DLSAT(TR); if (PR > 0.1) RHOV = DVSAT(TR); endif endif if ((PR < P0) || (PR < PBOTR)) RGUESS = PR*PCW*1.0e3 / (RW*TR*TCW); RGUESS = RGUESS / RHOCW; if (PR > 0.1) RGUESS = RGUESS + (PR/P0)*(RHOV-RGUESS); endif else RGUESS = RHOL; endif else RGUESS = PR*PCW*1.0e3 / (RW*TR*TCW); RGUESS = RGUESS / RHOCW; endif IERR = DFIND(DOUT, PR, RGUESS, TR, DPD, IWANT, PROPR); return IERR End Static Function DFIND1(double &DOUT, double P, double D, double TR, double &DPD, wave IWANT, wave PROPR) double TCLOSE = 1.1; double TOL = 1.0e-8; int IERR; double DD, PNEW, DOLD, ERROLD, STEP, DPNEW, ERRNEW, DNEW, FLOW, FHIGH, FRDIST; double POLD, PHIGH, DHIGH, PLOW, DLOW, FNEW; DD = D; if (DD <= 1.0e-8) DD = 1.0e-8; endif if (DD > 5.5) DD = 5.5; endif PDP(DD, TR, PNEW, DPD, IWANT, PROPR); if (TR > TCLOSE) { DOLD = DD; ERROLD = PNEW - P; STNO_10: STEP = -ERROLD / *DPD; if (fabs(STEP) > 0.5*DOLD) STEP = copysign(0.5*DOLD, STEP); DNEW = DOLD + STEP; STNO_20: PDP(DNEW, TR, &PNEW, &DPNEW, IWANT, PROPR); /* C CHECK FOR CONVERGENCE */ if (fabs(1.0 - PNEW/P) < TOL) IERR = 0; DOUT = DNEW; DPD = DPNEW; return IERR endif /* C IF NOT CONVERGED, SEE IF ERROR HAS DECREASED. IF NOT, CUT STEP IN HALF C IF IT HAS, TAKE THE NEXT NEWTON STEP */ ERRNEW = PNEW - P; if (fabs(ERRNEW) > fabs(ERROLD)) { DNEW = 0.5*(DNEW + DOLD); goto STNO_20; } ERROLD = ERRNEW; DOLD = DNEW; *DPD = DPNEW; goto STNO_10; } else { /* C IF CLOSE TO TC, WORK BY BOUNDING THE SOLUTION */ if (PNEW > P) { /* C IF P TOO HIGH, GO DOWN IN DENSITY UNTIL SOLUTION BOUNDED */ DOLD = DD; POLD = PNEW; for (int i=1; i<=9; i++) { DNEW = (10-i)*0.1*DD; PDP(DNEW, TR, &PNEW, &DPNEW, IWANT, PROPR); if (PNEW <= P) { PHIGH = POLD; DHIGH = DOLD; PLOW = PNEW; DLOW = DNEW; goto STNO_222; } DOLD = DNEW; POLD = PNEW; } /* C IF YOU GET HERE, TRY SUCCESSIVE HALVING OF THE DENSITY */ for (int i=1; i<=20; i++) { DNEW = 0.5*DOLD; PDP(DNEW, TR, &PNEW, &DPNEW, IWANT, PROPR); if (PNEW <= P) { PHIGH = POLD; DHIGH = DOLD; PLOW = PNEW; DLOW = DNEW; goto STNO_222; } DOLD = DNEW; POLD = PNEW; } /* C IF YOU STILL HAVEN'T BRACKETED IT, GIVE UP */ IERR = -2; *DOUT = DNEW; *DPD = DPNEW; return IERR; } else { /* C IF P TOO LOW, GO UP IN DENSITY UNTIL SOLUTION BOUNDED */ DOLD = DD; POLD = PNEW; for (int i=1; i<=100; i++) { DNEW = (10 + i)*0.1*DD; PDP(DNEW, TR, &PNEW, &DPNEW, IWANT, PROPR); if (PNEW >= P) { PLOW = POLD; DLOW = DOLD; PHIGH = PNEW; DHIGH = DNEW; goto STNO_222; } DOLD = DNEW; POLD = PNEW; } /* C IF YOU STILL HAVEN'T BRACKETED IT, GIVE UP */ IERR = -2; *DOUT = DNEW; *DPD = DPNEW; return IERR; } STNO_222: /* C SOLUTION IS NOW BRACKETED, USE BOUNDED SECANT TO SOLVE */ FLOW = PLOW - P; FHIGH = PHIGH - P; int i for (i=1; i<=100; i++) FRDIST = -FLOW / (FHIGH-FLOW); if (FRDIST < 0.01) FRDIST = 0.01; endif if (FRDIST > 0.99) FRDIST = 0.99; endif DNEW = DLOW + FRDIST*(DHIGH - DLOW); PDP(DNEW, TR, PNEW, DPNEW, IWANT, PROPR); if (fabs(1.0 - PNEW/P) < TOL) IERR = 0; DOUT = DNEW; DPD = DPNEW; return IERR endif FNEW = PNEW - P; if (FNEW >= 0.0) FHIGH = FNEW; DHIGH = DNEW; else FLOW = FNEW; DLOW = DNEW; endif endfor IERR = -3; DOUT = DNEW; DPD = DPNEW; return IERR; endif End