ESTE PROGRAMA ESCRITO EN EL LENGUAJE DE PROGRAMACIóN WOLFRAM, VERSIóN MATHEMATICA 11.1, SE BASA íNTEGRAMENTE EN LA SIGUIENTE REFERENCIA (Versión 2, programada en Marzo de 2017):

International Association for the Properties of Water and Steam, IAPWS R6-95(2016), Revised Release on the IAPWS Formulation 1995 for the Thermodynamic Properties of Ordinary Water Substance for General and Scientific Use (2016). Dresden, Germany, 11-16 September, 2016.

“ This formulation provides the most accurate representation of the thermodynamic properties of the fluid phases of water substance over a wide range of conditions available at the time this release was prepared. ”

ESTA VERSIóN 2 SE TERMINó DE PROGRAMAR LIBRE DE ERRORES, EL 22 DE MARZO DE 2017. LA VERSIóN 1 PREVIA, FUE HECHA POR JOSE DARIO ROJAS HERNANDEZ EN SU TESIS DE LICENCIATURA: “ Programación de Propiedades Termofísicas del Agua Geotérmica usando el Potencial de Helmholtz (IAPWS-95) en el Rango p [0.1 MPa, 100 GPa] y T [0, 5000] °C ”  PRESENTADA EN 2013 EN MORELIA Y EXHIBIDA EN LA REUNION DE LA International Geothermal Association IGA EN MELBOURNE, AUSTRALIA EN Abril, 2015. LA DIFERENCIA ENTRE AMBAS VERSIONES ES QUE ESTA 2da ES UNA VERSION AMPLIADA Y MEJORADA, DONDE AGREGUé EL CáLCULO DE: COMPRESIBILIDAD ISOTERMICA, EXPANSIVIDAD Y CONDUCTIVIDAD TéRMICA, VELOCIDAD DEL SONIDO Y LA PRESENTACIóN DE RESULTADOS ES EN FORMA DE TABLA. AGREGUE LA CAPACIDAD DE INVERTIR LA DEPENDENCIA DE Φ(δ, τ) a Φ (δ(p, T), τ) USANDO EL ALGORITMO LINEAL CLASICO DE NEWTON-RAPHSON. PUEDE ADEMAS CALCULAR LA TEMPERATURA DE SATURACION DADA CUALQUIER PRESION EN LA LINEA DE SATURACION K o REGION BIFASICA. EL PROGRAMA SE VERIFICó CORRECTAMENTE CON LOS DATOS PUBLICADOS EN LA REFERENCIA ANTERIOR Y CON RESULTADOS QUE ARROJA EL SOFTWARE DE NIST/ASME DE JUNIO 2013.

1) En una sola fase, líquido o vapor, las variables independientes de entrada para el programa son la densidad en kg/(m^3) y la temperatura en K. En la línea de Saturación o Función K, donde ambas fases coexisten, la variable independiente es la la temperatura T en K. En esta versión ampliada se puede calcular (Subscript[T, sat])(p) con un algoritmo de inversión preciso, incluso en la vecindad del punto crítico. En ambos casos la  temperatura se introduce en °C por ser la unidad de uso común en geotermia.

2) En una o dos fases, las propiedades calculadas son: Temperatura, presión, densidad (en dos fases), entalpía específica, energía interna, entropía específica, calores específicos isocórico e isobárico, conductividad térmica, expansividad, viscosidad dinámica, compresibilidad isotérmica, módulo volumétrico, velocidad del sonido y tensión superficial (en dos fases). En total, junto con la temperatura, 16 funciones son presentadas en forma tabular para ambas fases.

3) Debido a que las unidades de campo comúnes en geotermia son la presión y la temperatura del fluido, esta versión contiene un algoritmo de inversión eficiente, calculando 1ro la densidad del líquido aproximándola con la IFC67, en función de presión y temperatura como variables de entrada. Enseguida se corrige la presión variando la densidad con la IAPWS95. La mejor técnica de inversión tanto para el vapor como para el líquido es el algoritmo de linearización de Newton-Raphson, pues utiliza derivadas parciales ya incluidas en la misma IAPWS95. Así el usuario tiene la opción desde el principio de ingresar ya sea (ρ, T) o  (p, T). Por ello incluyo en la misma celda, el cálculo de la densidad ρ(p, T) de ambas fases, usando la IFC-67 como primer valor aproximado para iniciar el algoritmo de inversión (backward formulation).

4) Este es un programa de cálculo muy potente, utiliza varias de las características avanzadas del lenguaje de programación Wolfram contenido en el software Mathematica versión 11.1, sin embargo se mantiene su compatibilidad con versiones anteriores, de la 6 en adelante. Es un programa para uso científico general, en el rango indicado arriba. En geotermia puede usarse certeramente para toda clase de sistemas geotérmicos de baja, media, alta entalpía, así como en sistemas supercríticos. Está ciertamente indicado en la amplia región geotérmica que abarca [0.1, 10,000] bar y [0, 1000] °C. Su aplicabilidad al agua magmática arriba de este rango es por extrapolación y hasta 1 millón de bars y 4727 °C. Fuera de estos rangos, la formulación IAPWS-95 ya no es aplicable.

ClearAll["Global`*"]; (* Memoria global y variables limpias. Versión 23/03/2017 *)
(* ************************************************************************ *)
(* DENSIDAD DEL AGUA EN FUNCIÓN DE (p, T) USANDO EL POTENCIAL REDUCIDO DE GIBBS (Subscript[G, 1])(Θ = Θv,β = βr) TANTO EN LA REGIÓN DE LÍQUIDO COMO EN LA ZONA DE VAPOR (Datos de Schmidt & Grigull,1979. Properties of Water & Steam, formulación IFC 1967)*)
(* ************************************************************************ *)
(* 1.- Densidad ρl(p,T) de la fase líquida, formulación IFC_1967 *)
Liquid[T0_,P0_]:=Module[{T=T0,p=P0},
AL={6824.687741,-542.2063673,-20966.66205,39412.86787,-67332.77739,99023.81028,-109391.1774,85908.41667,-45111.68742,14181.38926,-2017.271113,7.982692717,-0.02616571843,1.52241179 (10^(-3)),2.284279054 (10^(-2)),242.1647003,1.269716088 (10^(-10)),2.074838328 (10^(-7)),2.17402035 (10^(-8)),1.105710498 (10^(-9)),12.93441934,1.308119072 (10^(-5)),6.047626338 (10^(-14))};   sA={0.8438375405,5.362162162 (10^(-4)),1.72000,7.342278489 (10^(-2)),4.975858870 (10^(-2)),0.65371543,1.1500 (10^(-6)),1.51080 (10^(-5)),0.141880,7.002753165,2.995284926 (10^(-4)),0.2040};
tr=((T+273.15)/647.096); (* Temperatura reducida[K ] *) Pr=(p/(220.64 (10^5))); (* Presión reducida [Pa] *)
YL=1.0-sA[[1]] (tr^2)-sA[[2]](tr^(-6));
ZP=sA[[3]] (YL^2)-2.0 sA[[4]] tr+2.0 sA[[5]]Pr;
ZL=YL+Sqrt[ZP]; PAR1=AL[[12]](((sA[[5]])/(ZL^((5.0/17)))));
PAR2=AL[[13]]+AL[[14]]tr+AL[[15]](tr^2)+AL[[16]](((sA[[6]]-tr))^10)+(((AL[[17]])/(sA[[7]]+(tr^19))));
PAR3=(AL[[18]]+2.0 AL[[19]]Pr+3.0 AL[[20]] (Pr^2))/(sA[[8]]+(tr^11));
PAR4=AL[[21]](tr^18)(sA[[9]]+(tr^2))(((-3.0)/(((sA[[10]]+Pr))^4))+sA[[11]]);PAR5=3.0 AL[[22]] (sA[[12]]-tr) (Pr^2)+4.0 (((AL[[23]])/(tr^20))) (Pr^3);
v1=(PAR1+PAR2-PAR3-PAR4+PAR5) 3.17 (10^(-3));
If[Im[v1]≠0,v1=Re[v1]];  ρ1=1.0/v1;
]; (* Termina módulo densidad del Liquido[T,p] IFC_1967 *)
(* ************************************************************************ *)

(* 2.- Densidad ρv(p,T) de la fase vapor, formulación IFC_1967 *)
Steam[T0_,P0_]:=Module[{T=T0,p=P0},
B0=16.835992740; B01=28.560677960; B02=-54.389233290; B03=0.43306628340;B04=-0.6547711697; B05=0.085651820580; B11=0.066703759180; B12=1.3889838010;B21=0.083901043280; B22=0.026146708930; B23=-0.033734394530; B31=0.45209189040;B32=0.10690366140; B41=-0.59753367070; B42=-0.088475358040; B51=0.59580516090;B52=-0.51593033730; B53=0.20750211220; B61=0.11906102710; B62=-0.098671741320;B71=0.16839988030; B72=-0.058094380010; B81=6.552390126 (10^(-3)); B82=5.7102186490 (10^(-4)); B90=193.65875580; B91=-1388.5224250; B92=4126.6072190;B93=-6508.2116770; B94=5745.9840540; B95=-2693.0883650; B96=523.57186230;
SB=0.76333333330; SB61=0.40060739480; SB71=0.086360816270; SB81=-0.85323229210;SB82=0.34602088610;
θv=((T+273.15)/647.096);βr=(p/(22.064 (10^6))); I1=4.260321148;X=Exp[SB(1-θv)];
βL=15.74373327-34.17061978θv+19.31380707((θv)^2);RL=((βr)/(βL)); 
ξ2=I1((θv)/(βr)); SC=(B11*(X^10)+B12)(X^3); ξ2=ξ2-SC; SC=B21*(X^18)+B22*(X^2)+B23*X; ξ2=ξ2-2βr*SC; SC=(B31*(X^8)+B32)(X^10); ξ2=ξ2-3((βr)^2)SC; SC=(B41*(X^11)+B42)(X^14); ξ2=ξ2-4((βr)^3)SC;SC=(B51*(X^8)+B52*(X^4)+B53)(X^24); ξ2=ξ2-5((βr)^4)SC; SD1=((βr)^(-4))+SB61*(X^14); SD2=((βr)^(-5))+SB71*(X^19); SD3=((βr)^(-6))+(SB81*(X^27)+SB82)(X^27); SN=(B61*X+B62)(X^11);ξ2=ξ2-(SN/SD1*4((βr)^(-5)))/SD1; SN=(B71*(X^6)+B72)(X^18); ξ2=ξ2-(SN/SD2*5((βr)^(-6)))/SD2; SN=(B81*(X^10)+B82)(X^14); ξ2=ξ2-(SN/SD3*6((βr)^(-7)))/SD3;
SC=B96; SC=SC*X+B95; SC=SC*X+B94; SC=SC*X+B93; SC=SC*X+B92; SC=SC*X+B91; SC=SC*X+B90; ξ2=ξ2+11(RL^10)SC; 
vv=0.00317ξ2; If[MatchQ[vv,_Complex],vv=Re[vv]]; (* evita el plano complejo *) ρ2=1.0/vv  
]; (* Termina módulo densidad del Vapor[T,p] IFC_1967 *)
(* ************************************************************************ *)

(* _*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_ *)
(* _*_*_*_*_*_*_*_*_*_*_*_*_ IAPWS95_2016 _*_*_*_*_*_*_*_*_*_*_*_*_ *)
(* _*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_ *)
IAPWS95[ρ_,T_]:=Module[{δ=(ρ/322.0)(* Densidad reducida *),τ=(647.096/T)(* Temperatura reducida inversa *)},
Tc=647.096; (* [(Subscript[T, c])]= K, Temperatura crítica.  *)
Pc=22.064 (10^6); (* [(Subscript[P, c])]= Pa, Presión crítica. *)
ρc=322.0; (* [(Subscript[ρ, c])]= kg (m^(-3)), Densidad crítica. *)
μc=1.0 (10^(-6)); (* [(Subscript[μ, c])]= Pa.s, Viscosidad crítica. *)
λc=1.0 (10^(-3)); (* W/m/K, conductividad de referencia. *)
R=0.46151805; (* kJ/kg/K specific gas constant *)

(* Coeficientes y parámetros de la parte de gas ideal (ϕ^0) de la energía libre de Helmholtz adimensional f(ρ,T)/RT = (ϕ^0)[δ,τ] + (ϕ^r)[δ,τ] = ϕ[δ,τ]. *)
n0={-8.3204464837497,6.6832105275932,3.00632,0.012436,0.97315,1.27950,0.96956,0.24873};
γ0={0,0,0,1.28728967,3.53734222,7.74073708,9.24437796,27.5075105};

(* Parte del gas ideal (ϕ^0)[δ,τ] de la energía libre de Helmholtz adimensional. *)
(* ********************************************************** *)
ϕ0[δ_,τ_]:=Log[δ]+n0[[1]]+n0[[2]]τ+n0[[3]]Log[τ]+(∑_(i=4)^8)n0[[i]]Log[1-Exp[-γ0[[i]]τ]];
(* ********************************************************** *)

(* Derivadas de la parte de gas ideal (ϕ^0) de la energía libre de Helmholtz adimensional. *)
ϕ0δ[δ_,τ_]:=(1/δ); (* (Subscript[[((∂(ϕ^0))/(∂δ))], τ]) *)  ϕ0δδ[δ_,τ_]:=-(1/(δ^2)); (* (Subscript[[(((∂^2)(ϕ^0))/(∂(δ^2)))], τ]) *)
ϕ0τ[δ_,τ_]:=n0[[2]]+((n0[[3]])/τ)+(∑_(i=4)^8)n0[[i]]γ0[[i]]((((1-Exp[-γ0[[i]]τ]))^(-1))-1); (* (Subscript[[((∂(ϕ^0))/(∂τ))], δ]) *)
ϕ0ττ[δ_,τ_]:=-((n0[[3]])/(τ^2))-(∑_(i=4)^8)n0[[i]](((γ0[[i]]))^2)Exp[-γ0[[i]]τ](((1-Exp[-γ0[[i]]τ]))^(-2)); (* (Subscript[[(((∂^2)(ϕ^0))/(∂(τ^2)))], δ]) *)
ϕ0δτ=0.0; (* [(((∂^2)(ϕ^0))/(∂δ∂τ))] *)

(* Coeficientes y parámetros para la parte residual (ϕ^r) de la energía libre de Helmholtz adimensional. *)
Array[a,2,55]; a[55]=a[56]=3.5; Array[b,2,55]; b[55]=0.85; b[56]=0.95;
c={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,0,0,0};
d={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};
t={-0.5,0.875,1,0.5,0.75,0.375,1,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};
n={0.12533547935523 (10^(-1)),7.8957634722828,-0.87803203303561 (10^1),0.31802509345418,-0.26145533859358,-0.78199751687981 (10^(-2)),0.88089493102134 (10^(-2)),-0.66856572307965,0.20433810950965,-0.66212605039687 (10^(-4)),-0.19232721156002,-0.25709043003438,0.16074868486251,-0.40092828925807 (10^(-1)),0.39343422603254 (10^(-6)),-0.75941377088144 (10^(-5)),0.56250979351888 (10^(-3)),-0.15608652257135 (10^(-4)),0.11537996422951 (10^(-8)),0.36582165144204 (10^(-6)),-0.13251180074668 (10^(-11)),-0.62639586912454 (10^(-9)),-0.10793600908932,0.17611491008752 (10^(-1)),0.22132295167546,-0.40247669763528,0.58083399985759,0.49969146990806 (10^(-2)),-0.31358700712549 (10^(-1)),-0.74315929710341,0.47807329915480,0.20527940895948 (10^(-1)),-0.13636435110343,0.14180634400617 (10^(-1)),0.83326504880713 (10^(-2)),-0.29052336009585 (10^(-1)),0.38615085574206 (10^(-1)),-0.20393486513704 (10^(-1)),-0.16554050063734 (10^(-2)),0.19955571979541 (10^(-2)),0.15870308324157 (10^(-3)),-0.16388568342530 (10^(-4)),0.43613615723811 (10^(-1)),0.34994005463765 (10^(-1)),-0.76788197844621 (10^(-1)),0.22446277332006 (10^(-1)),-0.62689710414685 (10^(-4)),-0.55711118565645 (10^(-9)),-0.19905718354408,0.31777497330738,-0.11841182425981,-0.31306260323435 (10^2),0.31546140237781 (10^2),-0.25213154341695 (10^4),-0.14874640856724,0.31806110878444};
Array[α,3,52]; α[52]=α[53]=α[54]=20;
Array[β,5,52]; β[52]=β[53]=150;  β[54]=250;  β[55]= β[56]=0.3;
Array[γ,3,52]; γ[52]=γ[53]=1.21; γ[54]=1.25;
Array[A,2,55]; A[55]=A[56]=0.32;
Array[B,2,55]; B[55]=B[56]=0.2;
Array[Ci,2,55]; Ci[55]=28; Ci[56]=32; 
Array[Di,2,55]; Di[55]=700; Di[56]=800; (* Los simbolos C y D están protegidos en Mathematica *)

(* Función exponencial ψ y sus derivadas. *)
ψ[δ_,τ_]:=Exp[-Ci[i](((δ-1))^2)-Di[i](((τ-1))^2)];
(* ********************************************************** *)
ψδ[δ_,τ_]:=-2Ci[i](δ-1)ψ[δ, τ];  (* ((∂ψ)/(∂δ)) *)
ψδδ[δ_,τ_]:=(2Ci[i](((δ-1))^2)-1)2Ci[i]ψ[δ,τ]; (* (((∂^2)ψ)/(∂(δ^2))) *)
ψτ[δ_,τ_]:=-2Di[i](τ-1)ψ[δ,τ]; (* ((∂ψ)/(∂τ)) *)
ψττ[δ_,τ_]:=(2Di[i](((τ-1))^2)-1)2Di[i]ψ[δ,τ]; (* (((∂^2)ψ)/(∂(τ^2))) *)
ψδτ[δ_,τ_]:=4Ci[i]Di[i](δ-1)(τ-1)ψ[δ,τ]; (* (((∂^2)ψ)/(∂δ∂τ)) *)
(* ********************************************************** *)

(* Función Δ y sus derivadas, Theta es griega Mayúscula *)
Θ[δ_,τ_]:=(1-τ)+A[i]((((((δ-1))^2)))^(1/(2β[i]))); Δ[δ_,τ_]:=(((Θ[δ,τ]))^2)+B[i]((((((δ-1))^2)))^(a[i]));
(* ********************************************************** *)
Δδ[δ_,τ_]:=(δ-1)(A[i]((2Θ[δ,τ])/(β[i]))((((((δ-1))^2)))^(((1/(2β[i]))-1)))+2B[i]a[i]((((((δ-1))^2)))^(a[i]-1))); (* ((∂Δ)/(∂δ)) *)
Δδδ[δ_,τ_]:=((Δδ[δ,τ])/((δ-1)))+(((δ-1))^2)(4B[i]a[i](a[i]-1)((((((δ-1))^2)))^(a[i]-2))+2(((((A[i])/(β[i]))))^2)(((((((((δ-1))^2)))^(((1/(2β[i]))-1)))))^2)+4Θ[δ,τ]((A[i])/(β[i]))((1/(2β[i]))-1)((((((δ-1))^2)))^(((1/(2β[i]))-2)))); (* (((∂^2)Δ)/(∂(δ^2))) *)
(* ********************************************************** *)

(* Derivadas de la función distancia (Δ^(Subscript[b, i])). *)
Δbδ[δ_,τ_]:=b[i]((Δ[δ,τ])^(b[i]-1))Δδ[δ,τ]; (* ((∂(Δ^(Subscript[b, i])))/(∂δ)) *)
Δbδδ[δ_,τ_]:=b[i](((Δ[δ,τ])^(b[i]-1))Δδδ[δ,τ]+(b[i]-1)((Δ[δ,τ])^(b[i]-2))((Δδ[δ,τ])^2)); (* (((∂^2)(Δ^(Subscript[b, i])))/(∂(δ^2))) *)
Δbτ[δ_,τ_]:=-2Θ[δ,τ]b[i]((Δ[δ,τ])^(b[i]-1)); (* ((∂(Δ^(Subscript[b, i])))/(∂τ)) *)
Δbττ[δ_,τ_]:=2b[i]((Δ[δ,τ])^(b[i]-1))+4((Θ[δ,τ])^2)b[i](b[i]-1)((Δ[δ,τ])^(b[i]-2)); (* (((∂^2)(Δ^(Subscript[b, i])))/(∂(τ^2))) *)
Δbδτ[δ_,τ_]:=-A[i]b[i]*(2/(β[i]))((Δ[δ,τ])^(b[i]-1))(δ-1)((((((δ-1))^2)))^(((1/(2β[i]))-1)))-2Θ[δ,τ]b[i](b[i]-1)((Δ[δ,τ])^(b[i]-2))Δδ[δ,τ]; (* (((∂^2)(Δ^(Subscript[b, i])))/(∂δ∂τ)) *)
(* ********************************************************** *)
(* Parte residual (ϕ^r) de la energía libre de Helmholtz adimensional. *)
(* ********************************************************** *)
ϕr[δ_,τ_]:=(∑_(i=1)^7)n[[i]](δ^(d[[i]]))(τ^(t[[i]]))+(∑_(i=8)^51)n[[i]](δ^(d[[i]]))(τ^(t[[i]]))Exp[-(δ^(c[[i]]))]+(∑_(i=52)^54)n[[i]](δ^(d[[i]]))(τ^(t[[i]]))Exp[-α[i](((δ-1))^2)-β[i](((τ-γ[i]))^2)]+(∑_(i=55)^56)n[[i]]((Δ[δ,τ])^(b[i]))δ*ψ[δ,τ];
(* ********************************************************** *)

(* Derivadas de la parte residual (ϕ^(r )) de la energía libre de Helmholtz adimensional. *)
(* PRIMERA (Subscript[DERIVADA [((∂(ϕ^r))/(∂δ))], τ]) *)
ϕrδ[δ_,τ_]:=(∑_(i=1)^7)n[[i]]d[[i]](δ^((d[[i]]-1)))(τ^(t[[i]]))+(∑_(i=8)^51)n[[i]]Exp[-(δ^(c[[i]]))](δ^((d[[i]]-1)))(τ^(t[[i]]))(d[[i]]-c[[i]](δ^(c[[i]])))+(∑_(i=52)^54)n[[i]](δ^(d[[i]]))(τ^(t[[i]]))Exp[-α[i](((δ-1))^2)-β[i](((τ-γ[i]))^2)](((d[[i]])/δ)-2α[i](δ-1))+(∑_(i=55)^56)n[[i]](((Δ[δ,τ])^(b[i]))(ψ[δ,τ]+δ ψδ[δ,τ])+δ Δbδ[δ,τ]ψ[δ, τ]); (* (Subscript[[((∂(ϕ^r))/(∂δ))], τ]) *)
(* DERIVADA SEGUNDA ∂(Subscript[δδ [(((∂^2)(ϕ^r))/(∂(δ^2)))], τ]) *)
ϕrδδ[δ_,τ_]:=(∑_(i=1)^7)n[[i]]d[[i]](d[[i]]-1)(δ^((d[[i]]-2)))(τ^(t[[i]]))+(∑_(i=8)^51)n[[i]]Exp[-(δ^(c[[i]]))](δ^((d[[i]]-2)))(τ^(t[[i]]))((d[[i]]-c[[i]](δ^(c[[i]])))(d[[i]]-1-c[[i]](δ^(c[[i]])))-((c[[i]])^2)(δ^(c[[i]])))+(∑_(i=52)^54)n[[i]](τ^(t[[i]]))Exp[-α[i](((δ-1))^2)-β[i](((τ-γ[i]))^2)](-2α[i](δ^(d[[i]]))+4((α[i])^2)(δ^(d[[i]]))(((δ-1))^2)-4d[[i]]α[i](δ^(d[[i]]-1))(δ-1)+d[[i]](d[[i]]-1)(δ^(d[[i]]-2)))+(∑_(i=55)^56)n[[i]](((Δ[δ,τ])^(b[i]))(2ψδ[δ,τ]+δ ψδδ[δ,τ])+2Δbδ[δ,τ](ψ[δ,τ]+δ ψδ[δ,τ])+δ Δbδδ[δ,τ]ψ[δ, τ]); (* (Subscript[[(((∂^2)(ϕ^r))/(∂(δ^2)))], τ]) *)
(* PRIMERA DERIVADA ∂(Subscript[τ [((∂(ϕ^r))/(∂τ))], δ]) *)
ϕrτ[δ_,τ_]:=(∑_(i=1)^7)n[[i]]t[[i]](δ^(d[[i]]))(τ^(t[[i]]-1))+(∑_(i=8)^51)n[[i]]t[[i]](δ^(d[[i]]))(τ^(t[[i]]-1))Exp[-(δ^(c[[i]]))]+(∑_(i=52)^54)n[[i]](δ^(d[[i]]))(τ^(t[[i]]))Exp[-α[i](((δ-1))^2)-β[i](((τ-γ[i]))^2)](((t[[i]])/τ)-2β[i](τ-γ[i]))+(∑_(i=55)^56)n[[i]]δ(Δbτ[δ,τ]ψ[δ,τ]+((Δ[δ,τ])^(b[i]))ψτ[δ,τ]); (* (Subscript[[((∂(ϕ^r))/(∂τ))], δ]) *)
(* DERIVADA SEGUNDA ∂(Subscript[ττ [(((∂^2)(ϕ^r))/(∂(τ^2)))], δ]) *)
ϕrττ[δ_,τ_]:=(∑_(i=1)^7)n[[i]]t[[i]](t[[i]]-1)(δ^(d[[i]]))(τ^(t[[i]]-2))+(∑_(i=8)^51)n[[i]]t[[i]](t[[i]]-1)(δ^(d[[i]]))(τ^(t[[i]]-2))Exp[-(δ^(c[[i]]))]+(∑_(i=52)^54)n[[i]](δ^(d[[i]]))(τ^(t[[i]]))Exp[-α[i](((δ-1))^2)-β[i](((τ-γ[i]))^2)]((((((t[[i]])/τ)-2β[i](τ-γ[i])))^2)-((t[[i]])/(τ^2))-2β[i])+(∑_(i=55)^56)n[[i]]δ(Δbττ[δ,τ]ψ[δ,τ]+2Δbτ[δ,τ]ψτ[δ,τ]+((Δ[δ,τ])^(b[i]))ψττ[δ,τ]); (* (Subscript[[(((∂^2)(ϕ^r))/(∂(τ^2)))], δ]) *)
(* DERIVADA SEGUNDA ∂(Subscript[δτ [(((∂^2)(ϕ^r))/(∂δ ∂τ))], δ]) *)
ϕrδτ[δ_,τ_]:=(∑_(i=1)^7)n[[i]]d[[i]]t[[i]](δ^(d[[i]]-1))(τ^(t[[i]]-1))+(∑_(i=8)^51)n[[i]]t[[i]](δ^(d[[i]]-1))(τ^(t[[i]]-1))(d[[i]]-c[[i]](δ^(c[[i]])))Exp[-(δ^(c[[i]]))]+(∑_(i=52)^54)n[[i]](δ^(d[[i]]))(τ^(t[[i]]))Exp[-α[i](((δ-1))^2)-β[i](((τ-γ[i]))^2)](((d[[i]])/δ)-2α[i](δ-1))(((t[[i]])/τ)-2β[i](τ-γ[i]))+(∑_(i=55)^56)n[[i]](((Δ[δ,τ])^(b[i]))(ψτ[δ,τ]+δ ψδτ[δ,τ])+δ Δbδ[δ,τ]ψτ[δ, τ]+Δbτ[δ,τ](ψ[δ,τ]+δ ψδ[δ,τ])+δ Δbδτ[δ,τ]ψ[δ,τ]); (* [(((∂^2)(ϕ^r))/(∂δ∂τ))] *)

(* **************************************************************** *)
(* DERIVADA PARA CALCULAR CONDUCTIVIDAD, EXPANSIVIDAD, COMPRESIBILIDAD ISOTéRMICA e (Subscript[INVERSIóN [((∂δ)/(∂Π))], θ]) *)
δΠ[δ_,τ_]:=((Pc τ)/(R ρc Tc))(((1+2δ ϕrδ[δ,τ]+(δ^2)ϕrδδ[δ,τ]))^(-1));   (* deducida el 15/03/2017 *)
(* **************************************************************** *)

(* Relaciones entre las propiedades termodinámicas con la parte del gas ideal (ϕ^0) y la parte residual (ϕ^r) de la energía libre de Helmholtz adimensional f y sus derivadas. *)
p[δ_,τ_]:=(R ρ T)(1+δ ϕrδ[δ,τ]); (* Presión [p]=10 MPa. *)
u[δ_,τ_]:=R T τ(ϕ0τ[δ,τ]+ϕrτ[δ,τ]); (* Energía interna específica [u]=kJ (kg^(-1)). *)
h[δ_,τ_]:=R T(1+τ(ϕ0τ[δ,τ]+ϕrτ[δ,τ])+δ ϕrδ[δ,τ]); (* Entalpía específica [h]=kJ (kg^(-1)). *)
s[δ_,τ_]:=R(τ(ϕ0τ[δ,τ]+ϕrτ[δ,τ])-ϕ0[δ,τ]-ϕr[δ,τ]); (* Entropía específica [s]=kJ (kg^(-1)) (K^(-1)). *)
cv[δ_,τ_]:=-R (τ^2)(ϕ0ττ[δ,τ]+ϕrττ[δ,τ]); (* Capacidad calorífica isocórica específica [(Subscript[c, v])]=kJ (kg^(-1)) (K^(-1)). *)
cp[δ_,τ_]:=R(-(τ^2)(ϕ0ττ[δ,τ]+ϕrττ[δ,τ])+((((1+δ ϕrδ[δ,τ]-δ τ ϕrδτ[δ,τ]))^2)/(1+2δ ϕrδ[δ,τ]+(δ^2)ϕrδδ[δ,τ])));
(* Capacidad calorífica isobárica específica [(Subscript[c, p])]=kJ (kg^(-1)) (K^(-1)). *)
Kf[δ_,τ_]:= ((R ρc Tc)/τ)(δ+2(δ^2)ϕrδ[δ,τ]+(δ^3) ϕrδδ[δ,τ]);
Cf[δ_,τ_]:=(1/(Kf[δ,τ])); (* Compresibilidad (Pa^(-1)) y Módulo Volumétrico a θ constante *)
(* velocidad del sonido en m/s : *)
vs[δ_,τ_]:=(((R T(1+2δ ϕrδ[δ,τ]+(δ^2) ϕrδδ[δ,τ]-((((1+δ ϕrδ[δ,τ]-δ τ ϕrδτ[δ,τ]))^2)/((τ^2) (ϕ0ττ[δ,τ]+ϕrττ[δ,τ]))))(10^3)))^0.5);
(* expansividad térmica en 1/K *)
γf[δ_,τ_]:=(τ/Tc)(1+δ ϕrδ[δ,τ]-δ τ ϕrδτ[δ,τ])(((1+2δ ϕrδ[δ,τ]+(δ^2)ϕrδδ[δ,τ]))^(-1));
(* *************************************************************** *)

(* _*_*_*_*_* Conductividad Térmica (versión IAPWS95 09/2011) *_*_*_*_*_*_*_*_ *)
(* Constantes de la Región Crítica. *)
Λ=177.8514; qD=(((0.4 (10^(-9))))^(-1)); ν0=0.630; γl=1.239; ξ0=0.13 (10^(-9));Γ0=0.06; TR=1.5;  
(* Constantes para λ0 y λ1. *)
Array[L,{5,6},0];
(Subscript[L, 0])=2.443221 (10^(-3)); (Subscript[L, 1])=1.323095 (10^(-2)) ; (Subscript[L, 2])=6.770357 (10^(-3)); (Subscript[L, 3])=-3.454586 (10^(-3)); (Subscript[L, 4])=4.096266 (10^(-4)); 
L[0,0]=1.60397357; L[0,1]=−0.646013523;  L[0,2]=0.111443906;  L[0,3]=0.102997357;  L[0,4]=−0.0504123634; L[0,5]=0.00609859258; L[1,0]=2.33771842; L[1,1]=−2.78843778; L[1,2]=1.53616167; L[1,3]=−0.463045512; L[1,4]=0.0832827019; L[1,5]=−0.00719201245; 
L[2,0]=2.19650529; L[2,1]=−4.54580785; L[2,2]=3.55777244; L[2,3]=−1.40944978; L[2,4]=0.275418278; L[2,5]=−0.0205938816; L[3,0]=−1.21051378; L[3,1]=1.60812989; L[3,2]=−0.621178141; L[3,3]=0.0716373224; L[3,4]=0; L[3,5]=0; L[4,0]=−2.7203370; L[4,1]=4.57586331; L[4,2]=−3.18369245; L[4,3]=1.1168348; L[4,4]=−0.19268305; L[4,5]=0.012913842;
(* Primer factor (Subscript[λ, 0]) de la conductividad térmica. *)  λ0[τ_]:=((Sqrt[(τ^(-1))])/((∑_(i=0)^4)((Subscript[L, i])(τ^i)))); 
(* La contribución del 2do factor (Subscript[λ, 1]) es debida a la densidad: 
The second factor (Subscript[λ, 1]) of Kter represents the contribution to thermal conductivity due to finite density . *)
λ1[δ_,τ_]:=Exp[δ(∑_(i=0)^4)(((τ-1))^i)(∑_(j=0)^5)L[i,j](((δ-1))^j)];
(* Tercer factor (Subscript[λ, 2]) de la conductividad térmica (Critical Enhancement). Este factor para la región crítica lo corregí el 19 de Marzo de 2017 con la 2da correlación del NIST/ASME Steam Properties software en FORTRAN versión 3.0 de Junio 2013. Subrutina TRNS_TCX.FOR/function TCCR_2010  escrita por Allan H. Harvey & Eric W. Lemmon. La fórmula de la IAPWS95 está incompleta y no es correcta en toda la región crítica: *)

ΔX[δ_,τ_]:=δ (δΠ[δ,τ]-TR τ δΠ[δ,TR]);  (* critical enhancement is zero for low densities, to avoid complex numbers *)
If[ΔX[δ,τ]≤ 0.0,Λ=0.0,
ξ[δ_,τ_]:=ξ0 ((( ((ΔX[δ,τ])/(Γ0))))^((ν0)/(γl))); ξ1:=ξ[δ,τ]; Y[δ_,τ_]:=qD ξ[δ,τ] ;  y0:=Y[δ,τ] ;
If[ξ1>1.5 (10^(-12)),Z[y_]:=(2/(π y))(((1-(1/(κ[δ,τ])))ArcTan[y]+(y/(κ[δ,τ])))-(1-Exp[((-1)/((y^(-1))+((y^2)/(3(δ^2)))))])),
 Z[y_]:=(2/(π ))((y/2)+(y^2)((1/(3κ[δ,τ]))-(1/2))+(y^3)((1/(3(δ^2)))+(1/24))-(y^4)((1/(3(δ^2)))+(1/(5κ[δ,τ]))-(3/120))+(y^5)((1/(6(δ^2)))+(1/270))) ]
 
]; If[y0<1.2 (10^(-7)),Z[y0]=0.0];  κ[δ_,τ_]:= ((cp[δ,τ])/(cv[δ,τ]));
μf0:=μ[(ρ/322),(T/647.096)](10^(-6)); λ2[δ_,τ_]:=(Λ/τ) ((δ cp[δ,τ] μc)/(R μf0 )) Z[y0]; 

KT[δ_,τ_]:=λ0[τ]λ1[δ,τ]+λ2[δ,τ](* Conductividad Térmica *)

(* Fin de Ecuación para la conductividad térmica *)
(* _*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_ *)
];(*Fin del módulo IAPWS95_2016*)

(* _*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_ *)
(* _*_*_*_*_*_*_*_*_*_*_*_*_* VISCOSIDAD _*_*_*_*_*_*_*_*_*_*_*_*_ *)
(* _*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_ *)

Viscosidad[ρ_,T_]:=Module[{δ=(ρ/322), (* Densidad reducida. *)θ=(T/647.096) (* Temperatura reducida. Theta minúscula, no confundir con Θ mayúscula de la parte residual (ϕ^r) *)},
Tc=647.096;  (* [(Subscript[T, c])]= K, Temperatura crítica. *)
ρc=322.0; (* [(Subscript[ρ, c])]=kg (m^(-3)), Densidad crítica. *)
μc=1.0 (10^(-6)); (* [(Subscript[μ, c])]=Pa s, Viscosidad crítica. *)

(* Viscosidad de gas diluido (Subscript[μ, 0]). *)
(Subscript[H, 0])=1.67752; (Subscript[H, 1])=2.20462 ;(Subscript[H, 2])=0.6366564; (Subscript[H, 3])=-0.241605;
μ0[δ_,θ_]:=((100(Sqrt[θ]))/((∑_(i=0)^3)(((Subscript[H, i])/(θ^i)))));
(* Contribución a la viscosidad debida a la viscosidad finita (Subscript[μ, 1]). *)
Do[(Subscript[H, i,j])=0,{i,0,6},{j,0,6}];
(Subscript[H, 0,0])=5.20094*(10^(-1)); (Subscript[H, 1,0])=8.50895*(10^(-2)); (Subscript[H, 2,0])=-1.08374; (Subscript[H, 3,0])=-2.89555*(10^(-1)); (Subscript[H, 0,1])=2.22531*(10^(-1)); (Subscript[H, 1,1])=9.99115*(10^(-1)); (Subscript[H, 2,1])=1.88797; (Subscript[H, 3,1])=1.26613; (Subscript[H, 5,1])=1.20573*(10^(-1)); (Subscript[H, 0,2])=-2.81378*(10^(-1)); (Subscript[H, 1,2])=-9.06851*(10^(-1)); (Subscript[H, 2,2])=-7.72479*(10^(-1)); (Subscript[H, 3,2])=-4.89837*(10^(-1));(Subscript[H, 4,2])=-2.57040*(10^(-1)); (Subscript[H, 0,3])=1.61913*(10^(-1)); (Subscript[H, 1,3])=2.57399*(10^(-1));(Subscript[H, 0,4])=-3.25372*(10^(-2)); (Subscript[H, 3,4])=6.98452*(10^(-2)); (Subscript[H, 4,5])=8.72102*(10^(-3)); (Subscript[H, 3,6])=-4.35673*(10^(-3));(Subscript[H, 5,6])=-5.93264*(10^(-4));
μ1[δ_,θ_]:=Exp[δ(∑_(i=0)^5)((((1/θ)-1))^i)*(∑_(j=0)^6)(Subscript[H, i,j])(((δ-1))^j)];
(* Ecuación para la viscosidad. *)
μ[δ_,θ_]:=μ0[δ,θ]μ1[δ,θ];
];(*Termina módulo de Viscosidad*)

(* _*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_ *)
(* _*_*_*_*_*_*_*_*_*_*_ Presión de Saturación *_*_*_*_*_*_*_*_*_*_ *)
(* _*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_ *)
Saturacion[Temp_]:=Module[{T=Temp},
Tc=647.096; (* [(Subscript[T, c])]=K, Temperatura crítica. *)
ρc=322; (* [(Subscript[ρ, c])]=kg (m^(-3)), Densidad crítica. *)
pc=22.064; (* [(Subscript[p, c])]=MPa, Presión crítica. *)

(* Coeficientes para el cálculo de la presión de vapor (Subscript[p, σ]) y su derivada ((Subscript[dp, σ])/dT). *)
aσ={-7.85951783,1.84408259,-11.7866497,22.6807411,-15.9618719,1.80122502};
(* Coeficientes para el cálculo de la densidad de saturación líquida (ρ^(,)). *)
bσ={1.99274064,1.09965342,-0.510839303,-1.75493479,-45.5170352,-6.74694450*(10^5)};
(* Coeficientes para el cálculo de la densidad de saturación de vapor (ρ^(,,)). *)
cσ={-2.03150240,-2.68302940,-5.38626492,-17.2991605,-44.7586581,-63.9201063};
(* Coeficientes para el cálculo de la entalpía y entropía específica de saturación. *)
dσ={-5.65134998*(10^(-8)),2690.66631,127.287297,-135.003439,0.981825814,-1135.905627715,2319.5246};

(* Ecuación para la presión de saturación (Subscript[p, σ]) *)
pσ[T_]:=pc*Exp[(Tc/T)(aσ[[1]](1-(T/Tc))+aσ[[2]](((1-(T/Tc)))^1.5)+aσ[[3]](((1-(T/Tc)))^3)+aσ[[4]](((1-(T/Tc)))^3.5)+aσ[[5]](((1-(T/Tc)))^4)+aσ[[6]](((1-(T/Tc)))^7.5))]; (*[(Subscript[p, σ])]=(10^4) Pa.*)

(* Ecuación para la derivada de la presión de saturación ((Subscript[dp, σ])/dT). *)
pσT[T_]:=-((pσ[T])/T)(Log[((pσ[T])/pc)]+aσ[[1]]+1.5aσ[[2]](((1-(T/Tc)))^0.5)+3aσ[[3]](((1-(T/Tc)))^2)+3.5aσ[[4]](((1-(T/Tc)))^2.5)+4aσ[[5]](((1-(T/Tc)))^3)+7.5aσ[[6]]*(((1-(T/Tc)))^6.5)); (*[((Subscript[dp, σ])/dT)]=(10^4) Pa (K^(-1)).*)

(* Ecuación para la densidad del líquido en estado de saturación (ρ^(,)). *)
ρl[T_]:=ρc(1+bσ[[1]](((1-(T/Tc)))^(1/3))+bσ[[2]](((1-(T/Tc)))^(2/3))+bσ[[3]](((1-(T/Tc)))^(5/3))+bσ[[4]](((1-(T/Tc)))^(16/3))+bσ[[5]](((1-(T/Tc)))^(43/3))+bσ[[6]](((1-(T/Tc)))^(110/3))); (*[(δ^(,))]=kg (m^(-3)).*)

(* Ecuación para la densidad del vapor en estado de saturación (ρ^(,,)). *)
ρv[T_]:=ρc*Exp[cσ[[1]](((1-(T/Tc)))^(2/6))+cσ[[2]](((1-(T/Tc)))^(4/6))+cσ[[3]](((1-(T/Tc)))^(8/6))+cσ[[4]](((1-(T/Tc)))^(18/6))+cσ[[5]](((1-(T/Tc)))^(37/6))+cσ[[6]](((1-(T/Tc)))^(71/6))]; (*[(δ^(,,))]=kg (m^(-3)).*)

(* Ecuaciones auxiliares para la entalpía y entropía específicas de saturación. *)
ασ[T_]:=1000(dσ[[6]]+dσ[[1]]((((T/Tc)))^(-19))+dσ[[2]]((T/Tc))+dσ[[3]]((((T/Tc)))^4.5)+dσ[[4]]((((T/Tc)))^5)+dσ[[5]]((((T/Tc)))^54.5)); (*[(Subscript[α, σ])]=J (kg^(-1)).*)
ψσ[T_]:=(1000/Tc)(dσ[[7]]+(19/20)dσ[[1]]((((T/Tc)))^(-20))+dσ[[2]]Log[(T/Tc)]+(9/7)dσ[[3]]((((T/Tc)))^3.5)+(5/4)dσ[[4]]((((T/Tc)))^4)+(109/107)dσ[[5]]((((T/Tc)))^53.5)); (*[(Subscript[ψ, σ])]=J (kg^(-1)) (K^(-1)).*)

(* Ecuaciones para la entalpía específica de líquido y vapor saturado. *)
hl[T_]:=(ασ[T]+((T*(10^6)pσT[T])/(ρl[T]))); (*La derivada de (Subscript[p, σ]) está en (10^4) Pa (K^(-1)).*)
hv[T_]:=(ασ[T]+((T*(10^6)pσT[T])/(ρv[T]))); (*[h]=J (kg^(-1)).*)

(* Ecuaciones para la entropía específica de líquido y vapor saturado. *)
sl[T_]:=(ψσ[T]+(((10^6)pσT[T])/(ρl[T]))); (*La derivada de (Subscript[p, σ]) está en (10^4) Pa (K^(-1)).*)
sv[T_]:=(ψσ[T]+(((10^6)pσT[T])/(ρv[T]))); (*[s]=J (kg^(-1)) (K^(-1)).*)

(* Ecuaciones para la energía interna específica de líquido y vapor saturado. *)
ul[T_]:=(ασ[T]+(10^6)(((T*pσT[T]-pσ[T])/(ρl[T])))); (*La derivada de (Subscript[p, σ]) y (Subscript[p, σ]) está en (10^4) Pa (K^(-1)).*)
uv[T_]:=(ασ[T]+(10^6)(((T*pσT[T]-pσ[T])/(ρv[T])))); (*[u]=J (kg^(-1)).*)
(* Ecuación para la tensión superficial del agua ordinaria bifásica. *)
σ[T_]:=235.8(((1-(T/Tc)))^1.256)(1-0.625(1-(T/Tc))) (* [σ]=mN (m^(-1)). *)
];(*Fin del modulo Presión de Saturación*)
(* **************************************************************** *)

(* _*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_ *)
(* _*_*_*_*_*_*_*_*_*_*_*_*_*_ PropTermo _*_*_*_*_*_*_*_*_*_*_*_*_ *)
(* _*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_ *)

PropTermo[Nedo_]:=Module[{state=Nedo},
(* ********* Inicializo Vectores Termodinámicos ********** *)
Array[Agua,12,0]; Array[Likido,12,0]; Array[Vapor,12,0]; Array[Mezcla,12,0]; 
Agua[0]={" "," "," "," "," "," "," "," "," "," "," "," "," "," "," "," "};
datH2O={Agua[0]}; j=0; (* Inicializo bucle y vector para agua en una sola fase *)
Mezcla={Agua[0]};                (* Inicializo bucle y vector para agua en dos fases *)
(* ************** Termina inicialización **************** *)

(Label[1];
accion=ChoiceDialog[Style["Choice one Option: \n\n-1. - DISPLAY PROGRAM CAPABILITIES \n\n0. - FIND A THERMODYNAMIC ZONE GIVEN (p, T). \n\n1. - COMPUTE LIQUID PHASE THERMODYNAMICS  f (ρ ,T) \n\n2. - COMPUTE VAPOR PHASE THERMODYNAMICS  f (ρ ,T) \n\n3. - BACKWARD LIQUID PHASE THERMODYNAMICS  ρL(p,T)⟶ f (p,T) \n\n4. - BACKWARD STEAM PHASE THERMODYNAMICS  ρV(p,T)⟶ f (p,T) \n\n5. - COMPUTE STATE OF SATURATION FOR A GIVEN T (°C). \n\n6. - COMPUTE STATE OF SATURATION FOR A GIVEN P (bar) \n\n7. - PLOT A GIVEN PROPERTY IN 2D or 3D \n\n8. - EXIT",Blue,Bold,FontSize->14],{-1,0,1,2,3,4,5,6,7,EXIT}];
If[accion==EXIT,Print[DateString[]]Abort[]];
If[accion==-1,
inicio=ChoiceDialog["ESTE PROGRAMA CALCULA LAS SIGUIENTES PROPIEDADES TERMODINáMICAS DEL AGUA EN FASE LíQUIDA o VAPOR : \n\n Monofase en función de T (°C -> K) y densidad (kg/m3) \n o en función de la temperatura (°C -> K) y presión p (bar): \n\n01)  p - Presión (bar) o Densidad (kg/m3) \n02)  e - Energía Interna (kJ/kg) \n03)  h - Entalpía Específica (kJ/kg) \n04)  s - Entropía Específica (kJ/kg/K) \n05)  cp - Calor Específico Isobárico (kJ/kg/K) \n06)  cv - Calor Específico Isocórico (kJ/kg/K) \n07)  Kt - Conductividad Térmica (W/m/K) \n08)   μ - Viscosidad Dinámica (Pa.s) \n09)  Cf - Compresibilidad Isotérmica (1/bar) \n10)  Kf - Módulo Volumétrico (bar) \n11)  γf - Expansividad Térmica (1/K) \n12)  vs - Velocidad del Sonido (m/s)\n\nEn 2-fases curva de saturación en función de T (°C) \no en función de la presión pσ (bar) en zona bifásica: \n\n01)  pσ - Presión de Saturación o Tσ de Saturación \n02)  ρ - Densidades de Líquido y Vapor Saturado \n03)  u - Energías Internas de Líquido y Vapor Saturado \n04)  h - Entalpía de Líquido y Vapor Saturado \n05)  s - Entropía de Líquido y Vapor Saturado \n06)  μ - Viscosidad de Líquido y Vapor Saturado \n07)  cp - Capacidad Calorífica Isobárica del Líquido y Vapor \n08)  cv - Capacidad Calorífica Isocórica del Líquido y Vapor \n09)  Cf - Compresibilidad Isotérmica de Líquido y Vapor (1/bar) \n10)  Kf - Módulo Volumétrico de Líquido y Vapor (bar) \n11)  γf - Expansividad Térmica de Líquido y Vapor (1/K)  \n12)  vs - Velocidad del Sonido en Líquido y en Vapor (m/s) \n13)  σ - Tensión Superficial \n14)  Propiedades de la Mezcla Bifásica dada Xv, Sl o Hf",{Continuar,Terminar}];
If[inicio==Terminar,Print[DateString[]]Abort[]]; Goto[1] ];

If[accion==7,tipoG=InputString["Introduce el tipo de Gráfica <2D> o <3D>:"]; 
If[tipoG=="3D"||tipoG=="3d",var3=ChoiceDialog["Escoge una propiedad a graficar función de (ρ,T):   \n\n01)  p - Presión (bar)  \n02)  e - Energía Interna (kJ/kg) \n03)  h - Entalpía Específica (kJ/kg) \n04)  s - Entropía Específica (kJ/kg/K) \n05)  cp - Calor Específico Isobárico (kJ/kg/K) \n06)  cv - Calor Específico Isocórico (kJ/kg/K) \n07)  Kt - Conductividad Térmica (W/m/K) \n08)   μ - Viscosidad Dinámica (Pa.s) \n09)  Cf - Compresibilidad Isotérmica (1/bar) \n10)  Kf - Módulo Volumétrico (bar) \n11)  γf - Expansividad Térmica (1/K) \n12)  vs - Velocidad del Sonido (m/s) \n13) - Ninguna",{1,2,3,4,5,6,7,8,9,10,11,12,Salir}],var2=ChoiceDialog["Escoge propiedad a graficar en 2D función de ρ, p o T:  \n\n01)  p - Presión (bar)  \n02)  e - Energía Interna (kJ/kg) \n03)  h - Entalpía Específica (kJ/kg) \n04)  s - Entropía Específica (kJ/kg/K) \n05)  cp - Calor Específico Isobárico (kJ/kg/K) \n06)  cv - Calor Específico Isocórico (kJ/kg/K) \n07)  Kt - Conductividad Térmica (W/m/K) \n08)   μ - Viscosidad Dinámica (Pa.s) \n09)  Cf - Compresibilidad Isotérmica (1/bar) \n10)  Kf - Módulo Volumétrico (bar) \n11)  γf - Expansividad Térmica (1/K) \n12)  vs - Velocidad del Sonido (m/s) \n13) - Ninguna",{1,2,3,4,5,6,7,8,9,10,11,12,Salir}],
var1=InputString["Introduce la variable independiente ρ (ro), p o T"]; 

If[var3==Salir||var2==Salir,Print[DateString[]]Abort[]] ]]);

If[accion==5||accion==6,rp=InputString["¿Despliego propiedades de la Mezcla Bifásica <s,n> ?"];If[rp=="s"||rp=="S",
var=ChoiceDialog["La Zona Bifásica requiere una variable extra. Escoge una opción: \n\n1 - CALIDAD del Vapor (%) \n2 - SATURACIóN del Líquido (%). \n3 - ENTALPIA de la Mezcla (kJ/kg) \n4 - Ninguna",{1,2,3,Salir}];
If[var==Salir,Print[DateString[]]Abort[]];
If[var==1,Xv=Input["Introduce la Calidad del Vapor (ad): "]; Print["Xv = ",Xv]];
If[var==2,Sl=Input["Introduce la Saturación del Líquido (ad): "]; Print["Sl = ",Sl]];
If[var==3,Hf=Input["Introduce la Entalpía de la Mezcla (kJ/kg): "]; Print["Hm = ",Hf]]; 
Print[" ----------------------------------------------------------------------------- "];
Print["LOS VALORES DE LA MEZCLA BIFáSICA SE DESPLIEGAN CON: tabla2"] ] ];

(Label[iniciaB]; j=j+1;
If[accion==7,Goto[finB]];
If[accion≤4,
TC=Input[" - Introduce temperatura [°C]:",j];
While[TC==Null,TC=Input["Introduce temperatura [°C]:",j];];
If[TC==$Canceled,Print[DateString[]]Abort[]]; T=TC+273.15];

If[accion==1||accion==2,
ρ=Input[" -Introduce densidad [kg/m^3]: ",j];
While[ρ==Null,ρ=Input["Introduce densidad [kg/m^3]:",j];];
If[ρ==$Canceled,Print[DateString[]]Abort[]];
Viscosidad[ρ,T]; IAPWS95[ρ,T] ];

If[accion==3||accion==4||accion==0, i=0;
pB=Input["- Introduce presión [bar]: ",j];
While[pB==Null,pB=Input["Introduce presión [bar]:",j]];
If[pB==$Canceled,Print[DateString[]]Abort[]]; 

(* ______AQUI SE DETERMINA LA REGIóN TERMODINáMICA SEGUN (p,T)__________ *)
If[accion==0,Saturacion[T]; ps=pσ[T];If[MatchQ[ps,_Complex],ps=Re[ps]];
 If[10ps<pB,indT=1; Print["Estado Termodinámico: LíQUIDO COMPRIMIDO"]]; 
If[10ps>pB,indT=2; Print["Estado Termodinámico: VAPOR SUPERCALIENTE"]]; 
If[10ps==pB,indT=3;Print["Estado Termodinámico: ZONA BIFáSICA"]] ] ; 
(* ________________________________________________________________________ *)
(* Aproxima Densidad del Líquido y del Vapor ρ(p,T) con Newton-Raphson lineal, ρ0 con IFC-67 *)
If[accion==3||indT==1,Liquid[TC,(10^5)pB]; ρ0=ρ1]; If[accion==4||indT==2,Steam[TC,(10^5)pB]; ρ0=ρ2]; 
p0=(10^2)pB; T0=TC+273.15; τ0=(Tc/T0); δn=((ρ0)/(ρc)); (* 1er valor con IFC67 para iniciar iteración *)
IAPWS95[ρ0,T0];  (* Calcula 1er valor aprox con p(δn,τ0) IAPWS95 *)

(* Aproxima δ correspondiente a p0 con el algoritmo de Newton-Raphson lineal *) 
f[δ_]:=R ρc T0(δ+(δ^2)ϕrδ[δ,τ0])-p0; fp[δ_]:=R ρc T0(1+2δ ϕrδ[δ,τ0]+(δ^2)ϕrδδ[δ,τ0]);
(Label[10]; δn1= δn-((f[δn])/(fp[δn])); If[Abs[f[δn1]]≤(10^(-8))||i≥10,Goto[end],δn=δn1;i=i+1; Goto[10]]; Label[end]); ρ=δn1 ρc; 
Print[i," δn1 = ",δn1,", ρ = ",ρ,", pn = ",p[δn1,τ0](10^(-2)),", f[δn1] = ",f[δn1]];
(* Termina linearización de Densidad del Líquido/Vapor ρ(p,T) con Newton-Raphson  *)
(* _____________________________________________________________________ *)

Viscosidad[ρ,T]; IAPWS95[ρ,T]];

If[accion≤4,
If[accion==0||accion==3||accion==4,pf=pB,pf=SetPrecision[p[ρ/322,647.096/T](10^(-2)),5] ]; 
ρf=SetPrecision[ρ,7];
ef=SetPrecision[u[ρ/322,647.096/T],5]; hf=SetPrecision[h[ρ/322,647.096/T],5]; sf=SetPrecision[s[ρ/322,647.096/T],5];  cvf=SetPrecision[cv[ρ/322,647.096/T](10^3),5]; cpf=SetPrecision[cp[ρ/322,647.096/T](10^3),5]; ktf=SetPrecision[KT[ρ/322,647.096/T](10^(-3)),6];
μf=SetPrecision[μ[ρ/322,T/647.096],6](10^(-6)) ;
Cw=(10^2)/Kf[ρ/322,647.096/T]; Kw=Kf[ρ/322.0,647.096/T](10^(-2));
γw=γf[ρ/322,647.096/T]; vw=vs[ρ/322,647.096/T];δf=(ktf/(ρf cpf));
Agua[j]={TC,pf,ρf,hf,ef,sf,cvf,cpf,ktf,μf,Cw,Kw,γw,δf,vw,"indef"}; 
datH2O=Append[datH2O,Agua[j]];If[j<state,Goto[iniciaB]]  ];

If[accion==5,
TC=Input["Introduce temperatura de saturación [°C]:",j];
While[TC==Null,TC=Input["Introduce temperatura de saturación [°C]:"];];
If[TC==$Canceled,Print[DateString[]]Abort[]]; T=TC+273.15;
			While[T<273.16||T>647.096,
T=Input["Temperatura fuera de rango, valores Tϵ[0.1, 374]°C"];
			While[T==Null,T=Input["Introduce temperatura de saturación Tϵ[273.16, 647.096] K"];];
			If[T==$Canceled,Print[DateString[]]Abort[]];
			]];
If[accion==6, (* Temperatura de Saturación (Subscript[T, sat])(p) correspondiente a (Subscript[p, sat]) *)
PB=Input["Introduce presión de saturación en bar:"];  PX=(10^(-1))PB;  (* en MPa *)
 i=0; TX0=(4606.0/(24.02-Log[(10^5)PB]));   (* valor burdo en K para iniciar iteración *)
Ts=TX0; DT=(10^(-9))Ts; TsD=Ts+DT; Saturacion[Ts] ;
(Label[11]; i=i+1; (* Regresa p como presión de saturación con IAPWS95, MatchQ evita el plano complejo cerca del punto crítico *)
Ps=pσ[Ts]; If[MatchQ[Ps,_Complex],Ps=Re[Ps]]; If[Abs[((PX-Ps)/PX)]≤(10^(-12))||i>10,Goto[Fin]];
TsD=Ts+DT; PsD=pσ[TsD]; If[MatchQ[PsD,_Complex],PsD=Re[PsD]]; Ts=Ts+((PX-Ps)/(PsD-Ps))DT;
 Goto[ 11]; Label[Fin]); T=Ts; TC=SetPrecision[Ts-273.15,5];   Print["Tsat(",PB," bar) = ",SetPrecision[Ts-273.15,7],"°C"];
];  (* Termina cálculo de (Subscript[T, sat])(p) *)

If[accion==5||accion==6||indT==3,Saturacion[T];  (* Presento resultados en forma de Tabla. *)

(* PROPIEDADES DE LA FASE LíQUIDA (L) Y DE LA FASE VAPOR (V) *)
pS=SetPrecision[pσ[T]10,6]; σT=SetPrecision[(10^(-3))σ[T],5];
ρL=SetPrecision[ρl[T],7]; ρV=SetPrecision[ρv[T],7]; Viscosidad[ρl[T],T];
eL=SetPrecision[(10^(-3))ul[T],6]; eV=SetPrecision[(10^(-3))uv[T],6]; 
hL=SetPrecision[(10^(-3))hl[T],6]; hV=SetPrecision[(10^(-3))hv[T],6]; 
sL=SetPrecision[(10^(-3))sl[T],6];  sV=SetPrecision[(10^(-3))sv[T],6];
 IAPWS95[ρL,T];  
cvL=SetPrecision[cv[ρL/322,647.096/T](10^3),5]; cpL=SetPrecision[cp[ρL/322,647.096/T](10^3),5]; ktL=SetPrecision[KT[ρL/322,647.096/T](10^(-3)),6]; 
μL=SetPrecision[μ[ρL/322,T/647.096],6](10^(-6)) ;
CL=(10^2)/Kf[ρL/322,647.096/T]; KL=Kf[ρL/322.0,647.096/T](10^(-2)); 
vL=vs[ρL/322,647.096/T]; γL=γf[ρL/322,647.096/T]; δL=(ktL/(ρL cpL));
Viscosidad[ρv[T],T]; IAPWS95[ρV,T]; 
cvV=SetPrecision[cv[ρV/322,647.096/T](10^3),5]; cpV=SetPrecision[cp[ρV/322,647.096/T](10^3),5]; ktV=SetPrecision[KT[ρV/322,647.096/T](10^(-3)),6]; 
μV=SetPrecision[μ[ρV/322,T/647.096],6](10^(-6)) ;
CV=(10^2)/Kf[ρV/322,647.096/T]; KV=Kf[ρV/322.0,647.096/T](10^(-2)); 
vV=vs[ρV/322,647.096/T];  γV=γf[ρV/322,647.096/T]; δV=(ktV/(ρV cpV));
Likido[j]={TC,pS,ρL,hL,eL,sL,cvL,cpL,ktL,μL,CL,KL,γL,δL,vL, σT}; 
datH2O=Append[datH2O,Likido[j]];
Vapor[j]={TC,pS,ρV,hV,eV,sV,cvV,cpV,ktV,μV,CV,KV,γV,δV,vV,σT}; 
datH2O=Append[datH2O,Vapor[j]];

(* AQUI CALCULO PROPIEDADES DE LA MEZCLA BIFASICA LíQUIDO + VAPOR  *)
If[rp=="s"||rp=="S",
If[var==1,Hf=hL+Xv(hV-hL);Sl=((ρV(hV-Hf))/(Hf(ρL-ρV)-ρL hL+ρV hV)); ρF=ρL Sl+(1-Sl)ρV];
If[var==2,ρF=ρL Sl+(1-Sl)ρV; Xv=((ρV)/(ρF))(1-Sl);Hf=hL+Xv(hV-hL)];
If[var==3,Sl=((ρV(hV-Hf))/(Hf(ρL-ρV)-ρL hL+ρV hV)); ρF=ρL Sl+(1-Sl)ρV;  Xv=((ρV)/(ρF))(1-Sl)];
Ef=eL+Xv(eV-eL); Sf=sL+Xv(sV-sL); CVf=cvL Sl+(1-Sl)cvV; CPf=cpL Sl+(1-Sl)cpV;
KTf=ktL Sl+(1-Sl)ktV; μF=μL Sl+(1-Sl)μV; CW=CL Sl+(1-Sl)CV; KW=KL Sl+(1-Sl)KV;
γW=γL Sl+(1-Sl)γV; δF=δL Sl+(1-Sl)δV; VW=vL Sl+(1-Sl)vV;

Agua[j]={TC,pS,ρF,Hf,Ef,Sf,CVf,CPf,KTf,μF,CW,KW,γW,δF,VW,σT}; 
Mezcla=Append[Mezcla,Agua[j]]];
 If[j<state,Goto[iniciaB]] ]; 
Label[finB]); 
Print[" ----------------------------------------------------------------------------- "];

(Label[IniciaG];  If[accion<7,Clear[g2,g3];Goto[finG]];
If[accion==7&&tipoG=="2D"||tipoG=="2d", edo=InputString["Introduce el estado a graficar líquido <liq> o vapor <vap>:"];If[edo=="liq",ρmax=1000.0; ρmin=887.0]; If[edo=="vap",ρmax=0.6; ρmin=0.2],Goto[3]];
If[var2==1&&var1=="p"||var1=="T",Print["LAS VARIABLES A GRAFICAR SON INDEPENDIENTES EN UNA SOLA FASE"];Goto[IniciaG]];
T0=Input["Introduce una temperatura fija en °C"];  IAPWS95[ρ,T0+273.15];
(*If[var2==1&&var1=="ro",];*)
g2=Plot[p[(ρ/322),(647.096/(T0+273.15))](10^(-2)),{ρ,ρmin,ρmax},PlotLabel->"Pressure 1-phase -> STEAM",AxesLabel->{"ρ (kg/m3)","p (bar)"}];
Goto[finG];

Label[3]; ρmax=1000.0; ρmin=887.0; Tmax=800; Tmin=100;
T0=Input["Introduce una temperatura fija en °C"];  IAPWS95[ρ,T0+273.15];g3=Plot3D[p[(ρ/322),(647.096/(T+273.15))](10^(-2)),{ρ,ρmin,ρmax},{T,Tmin,Tmax},PlotLabel->"Presión del Líquido (bar)",AxesLabel->{"ρ (kg/m3)","T (°C)","p"}];
(*If[accion==7&&tipoG="3D"||tipoG="3d", edo=InputString["Introduce el estado a graficar líquido <liq> o vapor <vap>:"];
If[edo=="liq",ρmax=1000.0; ρmin=887.0]; If[edo=="vap",ρmax=0.6; ρmin=0.2];Tmax=800;Tmin=100,Goto[IniciaG]];
 IAPWS95[ρ,T]; g3=Plot3D[p[(ρ/322),(647.096/(T+273.15))](10^(-2)),{ρ,ρmin,ρmax},{T,Tmin,Tmax},PlotLabel->"Presión del Líquido (bar)",AxesLabel->{"ρ (kg/m3)","T (°C)","p"}];*)
Label[finG]); If[accion==7&&tipoG=="2d",Clear[g3]; g2];
If[accion==7&&tipoG=="3d",Clear[g2];g3];
];

Print["  TABLA CON PROPIEDADES TERMODINáMICAS DEL AGUA - IAPWS95, fecha: ",DateString[] ]; 

tabla2=
Text@Grid[Prepend[Mezcla,{"T (°C)","p  (bar)","ρ ((kg/(m^3)))","h ((kJ/kg))","e ((kJ/kg))","s ((kJ/kg))","(Subscript[C, V]) (\!\(\*FractionBox[\(J\), \(kg\\ K\)]\))","(Subscript[C, P]) (\!\(\*FractionBox[\(J\), \(kg\\ K\)]\))","\!\(\*SubscriptBox[\(K\), \(\(T\)\(\\ \\ \)\)]\)(\!\(\*FractionBox[\(W\), \(m\\ °C\)]\))","\!\(\*SubscriptBox[\(μ\), \(\(f\)\(\\ \)\)]\) (Pa·s)","(Subscript[C, f]) ((1/bar))","(Subscript[K, f]) (bar)","γf (1/K)","δf ((m^2)/s)","vs (m/s)","σT (N/m)"} ],Background->{None,{Lighter[Orange,.5],{White,Lighter[Blend[{Green,Orange}],.8]}}},Dividers->{{Darker[Gray,.6],{Lighter[Gray,.5]},Darker[Gray,.6]},{Darker[Gray,.6],Darker[Gray,.6],{False},Darker[Gray,.6]}},Alignment->{{Center}},ItemSize->Automatic,Frame->Darker[Gray,.6],ItemStyle->{12,{Bold},{Blue}},Spacings->{Automatic,.8}];

tabla1=Text@Grid[Prepend[datH2O,{"T (°C)","p  (bar)","ρ ((kg/(m^3)))","h ((kJ/kg))","e ((kJ/kg))","s ((kJ/kg))","(Subscript[C, V]) (\!\(\*FractionBox[\(J\), \(kg\\ K\)]\))","(Subscript[C, P]) (\!\(\*FractionBox[\(J\), \(kg\\ K\)]\))","\!\(\*SubscriptBox[\(K\), \(\(T\)\(\\ \\ \)\)]\)(\!\(\*FractionBox[\(W\), \(m\\ °C\)]\))","\!\(\*SubscriptBox[\(μ\), \(\(f\)\(\\ \)\)]\) (Pa·s)","(Subscript[C, f]) ((1/bar))","(Subscript[K, f]) (bar)","γf (1/K)","δf ((m^2)/s)","vs (m/s)","σT (N/m)"} ],Background->{None,{Lighter[Yellow,.5],{White,Lighter[Blend[{Blue,Green}],.8]}}},Dividers->{{Darker[Gray,.6],{Lighter[Gray,.5]},Darker[Gray,.6]},{Darker[Gray,.6],Darker[Gray,.6],{False},Darker[Gray,.6]}},Alignment->{{Center}},ItemSize->Automatic,Frame->Darker[Gray,.6],ItemStyle->{12,{Bold},{Blue}},Spacings->{Automatic,.8}];

(* Termina módulo PropTermo y construcción de Tablas Termodinámicas *)
(* ************************************************************************* *)

La celda siguiente hace uso del módulo PropTermo, el cual debe ejecutarse primero.
El argumento de PropTermo[Nedo] es el número de estados que se desean calcular.

fecha=Date[];  dm=fecha[[5]]; ds=fecha[[6]];

Print["ATENCIóN --> SI SE VA A GRAFICAR, INTRODUCE Nedo = 1"]
Nedo=Input["Introduce No. de estados a calcular:"];

PropTermo[Nedo]  (* Límites del Programa Exacto: 0.1°C ≤ T ≤ 1000°C, 99.605916952189602398 *) 

Print["No. de procesadores utilizados: ",$ProcessorCount,", Tipo de Procesador: ",$ProcessorType]; 
Print["Termino cálculos termodinámicos en ",Abs[dm-Date[][[5]] ]," minutos, ",Abs[ds-Date[][[6]] ]," segundos"]

g2

g3

tabla1

tabla2

(* Verificación de las funciones ϕ0, ϕr y sus derivadas (correctas) según datos de la IAPWS *)
{Column[{"ϕ0 =","ϕ0δ =","ϕ0δδ =","ϕ0τ =","ϕ0ττ =","ϕ0δτ ="}],SetPrecision[Column[{ϕ0[(838.025/322),(647.096/500)],ϕ0δ[(838.025/322),(647.096/500)],ϕ0δδ[(838.025/322),(647.096/500)],ϕ0τ[(838.025/322),(647.096/500)],ϕ0ττ[(838.025/322),(647.096/500)],ϕ0δτ}],10],Column[{"ϕr =","ϕrδ =","ϕrδδ =","ϕrτ =","ϕrττ =","ϕrδτ ="}],SetPrecision[Column[{ϕr[(838.025/322),(647.096/500)],ϕrδ[(838.025/322),(647.096/500)],ϕrδδ[(838.025/322),(647.096/500)],ϕrτ[(838.025/322),(647.096/500)],ϕrττ[(838.025/322),(647.096/500)],ϕrδτ[(838.025/322),(647.096/500)]}],10]}

tabla1

tabla2

En estas tablas se verifican los algoritmos de inversión (p, T) para cada fase: Líquido, Vapor y 2-Fases.

1.- Estas 1ras Tablas son para comparar el Líquido con f (p, T) y f (ρ,T) respectivamente. Solo los ultimos decimales difieren.

2.- Esta 2da Tabla es para el Vapor y Región Supercrítica para (p, T) y (ρ,T) respectivamente. Solo el último decimal difiere.

Clear[δ,τ,fd]
fd[δ_,τ_]=1+2δ ϕrδ[δ,τ]+(δ^2)ϕrδδ[δ,τ];  ((((1+δ ϕrδ[δ,τ]-δ τ ϕrδτ[δ,τ]))^2)/(1+2δ ϕrδ[δ,τ]+(δ^2)ϕrδδ[δ,τ]));
Plot3D[fd[δ,τ]/(10^15),{δ,0.75,1.3},{τ,1.9,2.0},PlotRange->Full,AxesLabel->Automatic,ColorFunction->"Rainbow"]

Clear[cp] (* Comportamiento extraño del cp en el rango indicado *)
cp[δ_,τ_]:=R(-(τ^2)(ϕ0ττ[δ,τ]+ϕrττ[δ,τ])+((((1+δ ϕrδ[δ,τ]-δ τ ϕrδτ[δ,τ]))^2)/(1+2δ ϕrδ[δ,τ]+(δ^2)ϕrδδ[δ,τ]))); (* δ=(ρ/322.0); [241.5, 418.6] τ=(647.096/T); [340.6,]*)
Plot3D[cp[δ,τ]/(10^18),{δ,0.75,1.3},{τ,1.98,2.1},PlotRange->Full,AxesLabel->Automatic,ColorFunction->"Rainbow"]

En las siguientes 3 celdas desarrollé los algoritmos de inversión (p, T) para las tres fases.

(* 1. Aproxima Densidad del Líquido ρ(p,T) con Newton-Raphson, ρ0 con IFC-67 *)
DateString[]
Tc=647.096; (* [(Subscript[T, c])]= K, Temperatura crítica.  *) ρc=322.0; (* [(Subscript[ρ, c])]= kg (m^(-3)), Densidad crítica. *)
R=0.46151805; (* kJ/kg/K specific gas constant *)
TC=Input["- Introduce temperatura [°C]: "]; pB=Input["- Introduce presión [bar]: "];
p0=(10^2)pB; T0=TC+273.15; 
Liquid[TC,(10^5)pB];ρ0=ρ1;  (* 1er valor con IFC67 para iniciar iteración *)
τ0=(Tc/T0); δn=((ρ0)/(ρc)); j=0;
IAPWS95[ρ0,T0];  (* Calcula 1er valor aprox para p(δ,τ0) *)

(* Aproxima δ correspondiente a p0 con el algoritmo de Newton-Raphson *) 
f[δ_]:=R ρc T0(δ+(δ^2)ϕrδ[δ,τ0])-p0; fp[δ_]:=R ρc T0(1+2δ ϕrδ[δ,τ0]+(δ^2)ϕrδδ[δ,τ0]);
(Label[1]; Print["j = ",j,", δn = ",δn,", ρ = ",δn ρc,", pn = ",p[δn,τ0](10^(-2)),", f[δn] = ",f[δn]];
 δn1= δn-((f[δn])/(fp[δn])); If[Abs[f[δn1]]≤(10^(-8))||j≥10,Goto[end],δn=δn1;j=j+1;Goto[1]]; 
Label[end]); 
Print["δn1 = ",δn1,", ρ = ",δn ρc,", pn = ",p[δn,τ0](10^(-2)),", f[δn1] = ",f[δn1]]

DateString[];

(* Aproxima Densidad del Vapor ρ(p,T) con Newton-Raphson, ρ0 con IFC-67 *)
DateString[]
Tc=647.096; (* [(Subscript[T, c])]= K, Temperatura crítica.  *) ρc=322.0; (* [(Subscript[ρ, c])]= kg (m^(-3)), Densidad crítica. *)
R=0.46151805; (* kJ/kg/K specific gas constant *)
TC=Input["- Introduce temperatura [°C]: "]; pB=Input["- Introduce presión [bar]: "];
p0=(10^2)pB; T0=TC+273.15; 
Steam[TC,p0];ρ0=ρ2;  (* valor IFC67 para iniciar iteración *)
τ0=(Tc/T0); δn=((ρ0)/(ρc)); j=0;
IAPWS95[ρ0,T0];  (* Calcula 1er valor aprox para p(δ,τ0) *)

(* Aproxima δ correspondiente a p0 con el algoritmo de Newton-Raphson *) 
f[δ_]:=R ρc T0(δ+(δ^2)ϕrδ[δ,τ0])-p0; fp[δ_]:=R ρc T0(1+2δ ϕrδ[δ,τ0]+(δ^2)ϕrδδ[δ,τ0]);
(Label[1]; Print["j = ",j,", δn = ",δn,", ρ = ",δn ρc,", pn = ",p[δn,τ0]10,", f[δn] = ",f[δn]];
 δn1= δn-((f[δn])/(fp[δn])); If[Abs[f[δn1]]≤(10^(-10))||j≥20,Goto[end],δn=δn1;j=j+1;Goto[1]]; 
Label[end]); 
Print["δn1 = ",δn1,", ρ = ",δn ρc,", pn = ",p[δn,τ0]10,", f[δn1] = ",f[δn1]]

MatchQ[1,_Complex]

(* 3.- Aproxima Temperatura de Saturación (Subscript[T, sat])(p) correspondiente a (Subscript[p, sat]) *)
(* Válida en la zona bifásica para 0.013 × (10^5) ≤ P0 ≤ 221.2 × (10^5) Pascales. *)
DateString[]

PB=Input["Introduce presión de saturación en bar:"];  PX=(10^(-1))PB;  (* MPa *)
 i=0; 
TX0=(4606.0/(24.02-Log[(10^5)PB]));   (* valor burdo en K para iniciar iteración *)
Ts=TX0; DT=(10^(-9))Ts; TsD=Ts+DT;  Saturacion[Ts] ;
Print["P0 = ",PX,", TX0 = ",Ts,", DT = ",DT,", TsD = ",TsD];
(Label[1]; i=i+1; (* Regresa p como presión de saturación en IAPWS95, MatchQ es para evitar el plano complejo cerca del punto crítico *) 
Ps=pσ[Ts];If[MatchQ[Ps,_Complex],Ps=Re[Ps]]; If[Abs[((PX-Ps)/PX)]≤(10^(-10))||i>10,Goto[Fin]];
TsD=Ts+DT; PsD=pσ[TsD];If[MatchQ[PsD,_Complex],PsD=Re[PsD]]; Ts=Ts+((PX-Ps)/(PsD-Ps))DT;
Print["Psat = ",Ps," MPa, Tsat = ",Ts," = ",Ts-273.15,"°C, DPX = ",((PX-Ps)/PX),", i = ",i];
 Goto[ 1];
Label[Fin]);  Print["Tsat = ",SetPrecision[Ts-273.15,8],"°C"]

??SetDirectory

SetDirectory["D:\\MarioCesar\\Documents\\IAPWS95-R6_2016\\DataH2O"]

Directory[]
FileNames[]

Helm=OpenRead["HelmholtzTD1b.txt"];
HTD=ReadList[Helm,{Number,Number,Number}];
MatrixForm[HTD]

ListPlot3D[HTD,InterpolationOrder->3,ColorFunction->Hue,Mesh->All,PlotLegends->BarLegend[Automatic,LegendMarkerSize->100,LegendLabel->"(Subscript[f, w]) "],PlotLabel->Style["Superficie de Helmholtz  (Subscript[f, w]) (T, ρ) a 1 bar",Bold,14],AxesLabel->{Style["T (°C)",Bold,Blue,11],Style["ρ (kg/(m^3))",Bold,Blue,11],Style["(Subscript[f, w]) (kJ/kg) ",Bold,Blue,11]}]

Close[Helm]

Helm=OpenRead["HelmholtzTD400b.txt"];
HTD=ReadList[Helm,{Number,Number,Number}];
MatrixForm[HTD]

ListPlot3D[HTD,InterpolationOrder->3,ColorFunction->Hue,Mesh->All,PlotLegends->BarLegend[Automatic,LegendMarkerSize->100,LegendLabel->"(Subscript[f, w]) "],PlotLabel->Style["Superficie de Helmholtz  (Subscript[f, w]) (T, ρ), 400 bar, 2700°C ",Bold,14],AxesLabel->{Style["T (°C)",Bold,Blue,11],Style["ρ (kg/(m^3))",Bold,Blue,11],Style["(Subscript[f, w]) (kJ/kg) ",Bold,Blue,11]}]

Close[Helm]

HTDL[[1,2]]

Bell=OpenRead["hpTbell1000.txt"];
HTD=ReadList[Bell,{Number,Number,Number}];
Dimensions[HTD]

sph1=Graphics3D[{Red,Sphere[]},ImageSize->15,Boxed->False];
ListPlot3D[HTD,InterpolationOrder->6,ColorFunction->Hue,Mesh->All,Epilog-> Inset[sph1,{.45,.425}],PlotLegends->BarLegend[Automatic,LegendMarkerSize->100,LegendLabel->"T (°C) "],PlotRange->{{0,4100},Full,{0,800}},ClippingStyle->None,PlotLabel->Style["Temperatura del Agua Geotérmica  T (h, p), 1000 bar, 800°C ",Bold,13],AxesLabel->{Style["h (kJ/kg)",Bold,Blue,12],Style["p (bar)       ",Bold,Blue,12],Style["T (°C) ",Bold,Blue,12]},Ticks->{Automatic,{0,200,400,600,800,1000},{0,200,400,600,800,1000}}]

Close[Bell]

500-273