Source: prq.js

/**
 * Funciones, métodos y propiedades empleadas en el diseño y dimensionamiento de equipos.
 * 
 * @autor Hernan Peñaranda Villavicencio <hernanpv@gmail.com>
 * @version 1.0.0
 * 
 */
"use strict";


/**
 * Clase predefinida Array
 * @external Array 
 * @see https://developer.mozilla.org/es/docs/Web/JavaScript/Reference/Global_Objects/Array
 */


/**
 * Calcula, con la ecuación de Antoine, la presión de vapor (en mmHg) si se manda la Temperatura (T) o la temperatura de saturación (en Kelvin) si se manda la presión (P). Por defecto el compuesto es el agua, para otros compuestos se deben enviar las constantes de Antoine: antA, antB y antC, obtenidas estando la temperatura en Kelvin y la presión en mmHg.
 * @param {Object} datos - Temperatura/presión y constantes de Antoine. 
 * @param {number} datos.T - Temperatura de operación en Kelvin.
 * @param {number} datos.P - Presión de operación en mmHg.
 * @param {number} [datos.antA=18.3036] - Primera constante de la ecuación de Antoine
 * @param {number} [datos.antB=3816.44] - Segunda constante de la ecuación de Antoine
 * @param {number} [datos.antC=-46.13] - Tercera constante de la ecuación de Antoine.
 * @returns {number} Presión de vapor en mmHg o temperatura de saturación en Kelvin.
 */
function antoine({T, P, antA=18.3036, antB=3816.44, antC=-46.13}={}) {
    return (P==undefined ? Math.exp(antA-antB/(T+antC)) : antB/(antA-Math.log(P))-antC);
} 





/**
 * Calcula, con la ecuación de Harlacher, la presión de vapor (en mmHg) si se manda la Temperatura (T) o la temperatura de saturación (en Kelvin) si se manda la presión (P). Por defecto el compuesto es el agua, para otros compuestos se deben enviar las constantes de Harlacher: harA, harB, harC y harD, obtenidas estando la temperatura en Kelvin y la presión en mmHg.
 * @param {Object} datos - Temperatura/presión y constantes de Harlacher.
 * @param {number} datos.T - Temperatura de operación en Kelvin.
 * @param {number} datos.P - Presión de operación en mmHg.
 * @param {number} [datos.harA=55.336] - Primera constante de la ecuación de Harlacher
 * @param {number} [datos.harB=-6869.50] - Segunda constante de la ecuación de Harlacher
 * @param {number} [datos.harC=-5.115] - Tercera constante de la ecuación de Harlacher
 * @param {number} [datos.harD=1.05] - Cuarta constante de la ecuación de Harlacher
 * @returns {number} Presión de vapor en mmHg o Temperatura de saturación en Kelvin.
 */
function harlacher({T, P, harA=55.336, harB=-6869.50, harC=-5.115, harD=1.05}={}) {
    if (T===undefined) {
        let xi = harB/(Math.log(P)-harA);
        const fp = T => harA+harB/T+harC*Math.log(T)+harD*P/(T*T)-Math.log(P);
        return fp.newton({xi});
    } else  {
        const xi = Math.exp(harA+harB/T+harC*Math.log(T));
        let ft = P => harA+harB/T+harC*Math.log(T)+harD*P/(T*T)-Math.log(P);
        return ft.newton({xi});
    }
};




/**
 * Calcula, con la ecuación de Thek-Stiel, la presión de vapor (en mmHg) si se manda la Temperatura (T) o la temperatura de saturación (en Kelvin) si se manda la presión (P). Por defecto el compuesto es el agua, para otros compuestos se deben enviar las siguientes propiedades del compuesto: Tb, temperatura de ebullición en Kelvin; Tc, temperatura crítica en Kelvin; Pc, presión crítica en atm y Hvb, entalpía de vaporización a la temperatura normal de ebullición en kj/kg.
 * @param {Object} datos - Temperatura/Presión y propiedades del compuesto.
 * @param {number} [datos.Tb=373.15] - Temperatura normal de ebullición en Kelvin.
 * @param {number} [datos.Tc=647.3] - Temperatura crítica en Kelvin.
 * @param {number} [datos.Pc=217.6] - Presión crítica en atm.
 * @param {number} [datos.Hvb=9717] - Entalpía de vaporización a la temperatura normal de ebullición en kJ/Kg.
 * @returns {number} Presión de vapor en mmHg o Temperatura de saturación en Kelvin.
 */
function thekStiel({T, P, Tb=373.15, Tc=647.3, Pc=217.6, Hvb=9717}={}) { 
    const R=1.987207;
    const Tbr = Tb/Tc;
    const A = Hvb/(R*Tc*(1-Tbr)**0.375);
    let B = 1.14893-1/Tbr-0.11719*Tbr-0.03174*Tbr*Tbr-0.375*Math.log(Tbr);
    const h = Tbr*Math.log(Pc)/(1-Tbr);
    const C = 5.2691+2.0753*A-3.1738*h;
    const ac = (0.46284*A+(Math.log(1/Pc)-A*B)/((Tbr**C-1)/C+0.040*(1/Tbr-1)))/1.042;
    const D = 1.042*ac-0.46284*A;
    if (T!==undefined) {
        const Tr = T/Tc;
        const B = 1.14893-1/Tr-0.11719*Tr-0.03174*Tr*Tr-0.375*Math.log(Tr);
        return Math.exp(A*B+D*((Tr**C-1)/C+0.040*(1/Tr-1)))*Pc*760;
    } else {
        const Pr = P/Pc/760;
        const ft = T => {
            const Tr = T/Tc;
            const B = 1.14893-1/Tr-0.11719*Tr-0.03174*Tr*Tr-0.375*Math.log(Tr);
            return A*B+D*((Tr**C-1)/C+0.040*(1/Tr-1))-Math.log(Pr);
        }
        return ft.newton({xi:Tb}); 
    }
}



/**
 * Calcula la humedad de un vapor en un gas, dadas la presión parcial (Pp) del vapor, la presión total (P) del gas, el peso molecular del vapor y el peso molecular del gas. Por defecto el peso molecular del vapor es el del agua y el peso molecular del gas es el del aire.
 * @param {Object} datos - Presiones y pesos moleculares.
 * @param {number} datos.Pp - Presiones parcial del líquido (mismas unidades que P)
 * @param {number} datos.P - Presión total del gas (mismas unidades que Pp)
 * @param {number} [datos.pmv=18.02] - Peso molecular del vapor (por defecto agua).
 * @param {number} [datos.pmg=28.97] - Peso molecular del gas seco (por defecto aire).
 * @returns {number} Humedad del gas
 */
function humedad({Pp, P, pmv=18.02, pmg=28.97}={}) {
    return pmv/pmg*Pp/(P-Pp);
}




/**
 * Calcula la humedad de un gas saturado, dadas la temperatura (T) en Kelvin, la presión (P) en mmHg y las constantes de la ecuación de Antoine: antA, antB y antC, obtenidas estando la temperatura en Kelvin y la presión en mmHg. Los valores por defecto para las constante de Antoine, corresponden al agua.
 * @param {Object} datos - Temperatura, presión y constantes de la ecuación de Antoine.
 * @param {datos.T} - Temperatura del gas en Kelvin.
 * @param {datos.P} - Presión del gas en mmHg.
 * @param {number} [datos.antA=18.3036] - Primera constante de la ecuación de Antoine
 * @param {number} [datos.antB=3816.44] - Segunda constante de la ecuación de Antoine
 * @param {number} [datos.antC=-46.13] - Tercera constante de la ecuación de Antoine.
 * @returns {number} humedad del gas saturado.
 */
function humSat({T, P, antA=18.3036, antB=3816.44, antC=-46.13}={}) {
    return humedad({Pp:antoine({T, antA, antB, antC}), P});
}




/**
 * Calcula el porcentaje (%) de humedad del gas, dadas la humedad (h) y la humedad de saturación (hs) del gas.
 * @param {Object} datos - Humedad del gas y humedad del gas saturado.
 * @param {number} datos.h - Humedad del gas
 * @param {number} datos.hs - Humedad del gas saturado
 * @returns {number} Porcentaje de humedad del gas
 */
function porHum({h, hs}={}) {
    return h/hs*100;
}




/**
 * Calcula la humedad relativa (%) de un compuesto en un gas, dadas la presión parcial del compuesto (Pp) y la presión de vapor del compuesto (Pv) a la temperatura del gas.
 * @param {Object} datos - Presión parcial y presión de vapor.
 * @param {number} datos.Pp - Presión parcial del compuesto en el gas (misma unidad que Pv)
 * @param {number} datos.Pv - Presión del vapor (misma unidad que Pp)
 * @returns {number} Humedad relativa del gas (%)
 */
function humRel({Pp, Pv}={}) {
    return Pp/Pv*100;
}




/**
 * Calcula la presión parcial de un compuesto (vapor) en un gas, dadas la humedad (h) del gas, la presión total del gas (P) y los pesos moleculares del compuesto y del gas. Pesos moleculares por defecto: agua (vapor) y aire (gas).
 * @param {Object} datos - Humedad, presión y pesos moleculares.
 * @param {number} datos.h - Humedad del gas.
 * @param {number} datos.P - Presión total del gas (unidad de presión).
 * @param {number} [datos.pmv=18.02] - Peso molecular del vapor
 * @param {number} [datos.pmg=28.97] - Peso molecular del gas seco
 * @returns {number} Presión parcial del gas (misma unidad que P)
 */
function prePar({h, P, pmv=18.02, pmg=28.97}={}) {
    return P/(1+pmv/(pmg*h));
}




/**
 * Calcula el volumen húmedo de un gas, dadas la temperatura (T), presión (P) y humedad (h) del gas, así como los pesos moleculares del vapor y del gas. Pesos moleculares por defecto: agua (vapor) y aire (gas).
 * @param {Object} datos - Temperatura, presión, humedad y pesos moleculares.
 * @param {number} datos.T - Temperatura del gas húmedo (K)
 * @param {number} datos.P - Presión del gas húmedo (mmHg)
 * @param {number} datos.h - Humedad del gas.
 * @param {number} [datos.pmv=18.02] - Peso molecular del vapor en el gas
 * @param {number} [datos.pmg=28.97] - Peso molecular del gas seco
 * @returns {number} Volumen del gas húmedo (m^3/kmol)
 */
function volHum({T, P, h, pmv=18.02, pmg=28.97}={}) {
    return (1/pmg+h/pmv)*62.3635984185*T/P;
}




/**
 * Calcula el calor sensible de un gas húmedo, para la humedad (h), capacidad calorífica promedio del gas y capacidad calorífica promedio del vapor. Capacidades caloríficas promedio por defecto: aire (gas), agua (vapor).
 * @param {Object} datos - humedad y capacidades caloríficas promedio.
 * @param {number} datos.h - Humedad del gas.
 * @param {number} [datos.cpg=1.005] - Capacidad calorífica promedio del gas en kJ/kg.
 * @param {number} [datos.cpv=1.88] - Capacidad calorífica promedio del vapor en kJ/kg.
 * @returns {number} Calor sensible del gas húmedo (kJ/kg*K)
 */
function csGas({h, cpg=1.005, cpv=1.88}={}) {
    return cpg+cpv*h;
}




/**
 * Calcula la entalpía de un gas húmedo, dadas la temperatura (T), humedad (h), capacidades caloríficas promedio del gas (cpg), capacidad calorífica promedio del vapor (cpv) y entalpía de vaporización (Hv) a la temperatura de solidificación (Ts) del vapor (estado de referencia temperatura de solidificación del vapor). Vapores por defecto: aire (gas), agua (vapor) 
 * @param {Object} datos - Temperatura, humedad, capacidades caloríficas promedio, entalpía de vaporización y temperatura de solidificación.
 * @param {number} datos.T - Temperatura del gas en Kelvin.
 * @param {number} datos.h - Humedad del aire.
 * @param {number} [datos.cpg=1.005] - Capacidad calorífica promedio del gas en kJ/kg.
 * @param {number} [datos.cpv=1.88] - Capacidad calorífica promedio del vapor en kJ/kg.
* @param {number} [datos.Hv=2501.4] - Entalpía de vaporización del vapor, a la temperatura de solidificación, en kJ/kg.
 * @returns {number} Entalpía del aire húmedo en kJ/kg de aire seco
 */
function entGas({T, h, cpg=1.005, cpv=1.88, Hv=2501.4, Ts=273.15}={}) {
    return csGas({h, cpg, cpv})*(T-Ts)+Hv*h;
}




/**
 * Calcula la entalpía del líquido, dadas la temperatura actual del líquido (T), la capacidad calorífica media del líquido y la temperatura de solidificación (Ts) del líquido. Líquido por defecto: agua.
 * @param {Object} datos - Temperatura del líquido, capacidad calorífica media y temperatura de solidificación del líquido.
 * @param {number} datos.T - Temperatura del líquido en Kelvin.
 * @param {number} [datos.Ts=273.15] - Temperatura de solidificación del líquido en Kelvin.
 * @param {number} [datos.cpl=4.1868] - Capacidad calorífica promedio del líquido en kJ/kg.
 * @returns {number} Entalpía del líquido (kJ/kg) 
 */
function entLiq({T, Ts=273.15, cpl=4.1868}={}) {
    return cpl*(T-Ts);
}




/**
 * Calcula la temperatura (T2) y/o la humedad (h2) a la salida de un saturador adiabático. Calcula la temperatura y la humedad si se mandan la presión (P) y el porcentaje de saturación (ps), calcula la temperatura si se manda la humedad final (h2) y calcula la temperatura si se manda la temperatura final. Por defecto, las propiedades del gas son las del aire y las del líquido las del agua.
 * @param {Object} datos - Temperatura, humedad, presión, porcentaje de saturación y propiedades del gas y vapor. 
 * @param {number} datos.T1 - Temperatura inicial del aire en Kelvin.
 * @param {number} datos.h1 - Humedad inicial del aire.
 * @param {number} [datos.P] - Presión inicial del aire en mmHg. 
 * @param {number} [datos.ps] - Porcentaje de humedad del aire a la salida del saturador.
 * @param {number} [datos.T2] - Temperatura del aire, en Kelvin, a la salida del saturador.
 * @param {number} [datos.h2] - Humedad del aire, en Kelvin, a la salida del saturador.
 * @param {number} [datos.antA=18.3036] - Primera constante de la ecuación de Antoine.
 * @param {number} [datos.antB=3816.44] - Segunda constante de la ecuación de Antoine.
 * @param {number} [datos.antC=-46.13] - Tercera constante de la ecuación de Antoine.
 * @param {number} [datos.cpg=1.005] - Capacidad calorífica promedio del gas en kJ/kg.
 * @param {number} [datos.cpv=1.88] - Capacidad calorífica promedio del vapor en kJ/kg.
 * @param {number} [datos.cpl=4.1868] - Capacidad calorífica promedio del líquido en kJ/kg.
 * @param {number} [datos.Ts=273.15] - Temperatura de solidificación del líquido en Kelvin.
* @param {number} [datos.Hv=2501.4] - Entalpía de vaporización del vapor, a la temperatura de solidificación, en kJ/kg.
 * @returns {Object} Temperatura (T2) en K y humedad (h2) del aire a la salida del saturador.
*/
function temSatAdiabática({T1, h1, P, ps, h2, T2, antA=18.3036, antB=3816.44, antC=-46.13, cpg=1.005, cpv=1.88, cpl=4.1868, Ts=273.15, Hv=2501.4}={}) {
    if (P!==undefined && ps!==undefined) {
        const Hy1 = entGas({T:T1, h:h1, cpg, cpv, Hv, Ts});
        let Pp, hs, h2, Hx2, Hy2;
        function ft2(T2) {
            Pp = antoine({T:T2, antA, antB, antC});
            hs = humedad({Pp, P});
            h2 = ps*hs/100;
            Hx2 = entLiq({T:T2, Ts, cpl});
            Hy2 = entGas({T:T2, h:h2, cpg, cpv, Hv, Ts});
            return Hy1+Hx2*(h2-h1)-Hy2;
        };
        const T2 = ft2.newton({xi:T1*0.90, err:12}).round(2);
        h2 = h2.round(6);
        return {T2, h2};
    } else if (h2!==undefined) {
        const Hy1 = entGas({T:T1, h:h1, cpg, cpv, Hv, Ts});
        let Hx2, Hy2;
        function ft2(T2) {
            Hx2 = entLiq({T:T2, Ts, cpl});
            Hy2 = entGas({T:T2, h:h2, cpg, cpv, Hv, Ts});
            return Hy1+Hx2*(h2-h1)-Hy2;
        };
        const T2 = ft2.newton({xi:T1*0.90, err:12}).round(2);
        return {T2};    
    } else if (T2!==undefined) {
        const Hy1 = entGas({T:T1, h:h1, cpg, cpv, Hv, Ts});
        const Hx2 = entLiq({T:T2, Ts, cpl});
        let Hy2;
        function ft2(h2) {
            Hy2 = entGas({T:T2, h:h2, cpg, cpv, Hv, Ts});
            return Hy1+Hx2*(h2-h1)-Hy2;
        }
        const h2 = ft2.newton({xi:h1/0.90, err:12}).round(6);
        return {h2};
    } else {
        throw new TypeError("Se debe enviar la Presión y porcentaje de saturación o la humedad final o la temperatura final.");
    }
}




/**
 * Calcula el calor latente de vaporización (entalpía de vaporización) a la temperatura T, conociendo la temperatura crítica (Tc) del compuesto y su entalpía de vaporización (Hv0) a la temperatura T0. El compuesto por defecto es el agua.
 * @param {Object} datos - Temperatura, temperatura crítica y entalpía de vaporización a la temperatura T0.
 * @param {number} datos.T - Temperatura del compuesto en Kelvin.
 * @param {number} [datos.Hv0=2501.4] - Entalpía de vaporización (calor latente de vaporización) del compuesto, a la temperatura T0.
 * @param {number} [datos.T0=173.15] - Temperatura en Kelvin, a la que se conoce la entalpía de vaporización.
 * @param {number} [datos.Tc=647.3] - Temperatura crítica del compuesto en K
 * @returns {number} Calor latente de vaporización del compuesto a la temperatura dada (en las mismas unidades que Hv0).
 */
function entVap({T, Hv0=2501.4, T0=273.15, Tc=647.3}={}) {
    return Hv0*((1-T/Tc)/(1-T0/Tc))**0.375;
}




/**
 * Calcula la temperatura (T2) y/o la humedad (h2) a la salida de un saturador adiabático. Calcula la temperatura y la humedad si se mandan la presión (P) y el porcentaje de saturación (ps), calcula la temperatura si se manda la humedad final (h2) y calcula la temperatura si se manda la temperatura final. Por defecto, las propiedades del gas son las del aire y las del líquido las del agua.
 * @param {Object} datos - Temperatura, humedad, presión, porcentaje de saturación y propiedades del gas y vapor. 
 * @param {number} datos.T1 - Temperatura inicial del aire en Kelvin.
 * @param {number} datos.h1 - Humedad inicial del aire.
 * @param {number} [datos.P] - Presión inicial del aire en mmHg. 
 * @param {number} [datos.ps] - Porcentaje de humedad del aire a la salida del saturador.
 * @param {number} [datos.T2] - Temperatura del aire, en Kelvin, a la salida del saturador.
 * @param {number} [datos.h2] - Humedad del aire, en Kelvin, a la salida del saturador.
 * @param {number} [datos.antA=18.3036] - Primera constante de la ecuación de Antoine.
 * @param {number} [datos.antB=3816.44] - Segunda constante de la ecuación de Antoine.
 * @param {number} [datos.antC=-46.13] - Tercera constante de la ecuación de Antoine.
 * @param {number} [datos.cpg=1.005] - Capacidad calorífica promedio del gas en kJ/kg.
 * @param {number} [datos.cpv=1.88] - Capacidad calorífica promedio del vapor en kJ/kg.
 * @param {number} [datos.Hv0=2501.4] - Entalpía de vaporización (calor latente de vaporización) del compuesto, a la temperatura T0.
 * @param {number} [datos.T0=173.15] - Temperatura en Kelvin, a la que se conoce la entalpía de vaporización.
 * @param {number} [datos.Tc=647.3] - Temperatura crítica del compuesto en K
 * @returns {Object} Temperatura (T2) en K y humedad (h2) del aire a la salida del saturador.
*/
function temSatAdiabática2({T1, h1, P, ps, h2, T2, antA=18.3036, antB=3816.44, antC=-46.13, cpg=1.005, cpv=1.88, Hv0=2501.4, T0=273.15, Tc=647.3}={}) {
    if (P!==undefined && ps!==undefined) {
        const cs1 = csGas({h:h1, cpg, cpv});
        let Pp, hs, h2;
        function ft2(T2) {
            Pp = antoine({T:T2, antA, antB, antC});
            hs = humedad({Pp, P});
            h2 = ps*hs/100;
            return (h1-h2)/(T1-T2)+cs1/entVap({T:T2, Hv0, T0, Tc});
        };
        const T2 = ft2.secante({xi:T1*0.98, xf:T1*1.01, err:12}).round(2);
        h2 = h2.round(6);
        return {T2, h2};
    } else if (h2!==undefined) {
        const cs1 = csGas({h:h1, cpg, cpv});
        function ft2(T2) {
            return (h1-h2)/(T1-T2)+cs1/entVap({T:T2, Hv0, T0, Tc});
        };
        const T2 = ft2.secante({xi:T1*0.98, xf:T1*1.01, err:12}).round(2);
        return {T2};
    } else if (T2!==undefined) {
        const cshv = csGas({h:h1, cpg, cpv})/entVap({T:T2, Hv0, T0, Tc});
        const dt = T1-T2;
        function fh2(h2) {
            return (h1-h2)/dt+cshv;
        };
        const h2 = fh2.secante({xi:T1*0.98, xf:T1*1.01, err:12}).round(6);
        return {h2};
    } else {
        throw new TypeError("Se debe enviar la Presión y porcentaje de saturación o la humedad final o la temperatura final.");
    }
}




/**
 * Calcula la humedad del aire en base a la temperatura del aire (T1), la temperatura de bulbo húmedo (T2) y la presión (P) a la que se encuentra.
 * @param {Object} datos - Parámetros nombrados, valores por defecto para el aire húmedo
 * @param {number} datos.T - Temperatura del aire en Kelvin.
 * @param {number} datos.Tbh - Temperatura de bulbo húmedo en Kelvin. 
 * @param {number} datos.P - Presión del aire en mmHg.
  * @param {number} [datos.antA=18.3036] - Primera constante de la ecuación de Antoine.
 * @param {number} [datos.antB=3816.44] - Segunda constante de la ecuación de Antoine.
 * @param {number} [datos.antC=-46.13] - Tercera constante de la ecuación de Antoine.
 * @param {number} [datos.cpg=1.005] - Capacidad calorífica promedio del gas en kJ/kg.
 * @param {number} [datos.cpv=1.88] - Capacidad calorífica promedio del vapor en kJ/kg.
 * @param {number} [datos.cpl=4.1868] - Capacidad calorífica promedio del líquido en kJ/kg.
 * @param {number} [datos.Ts=273.15] - Temperatura de solidificación del líquido en Kelvin.
* @param {number} [datos.Hv=2501.4] - Entalpía de vaporización del vapor, a la temperatura de solidificación, en kJ/kg.
 * @returns {Number} Humedad (h1) en kg de vapor de agua/kg de aire seco (humedad del aire)
*/
function humBulbo({T, Tbh, P, antA=18.3036, antB=3816.44, antC=-46.13, cpg=1.005, cpv=1.88, cpl=4.1868, Ts=273.15, Hv=2501.4}={}) {
    const Pp = antoine({T:Tbh, antA, antB, antC});
    const h2 = humedad({Pp, P});
    const Hx2 = entLiq({T:Tbh, Ts, cpl});
    const Hy2 = entGas({T:Tbh, h:h2, cpg, cpv, Hv, Ts});
    let Hy1;
    const fh1 = h => {
        Hy1 = entGas({T, h, cpg, cpv, Hv, Ts});
        return Hy1+Hx2*(h2-h)-Hy2;
    };
    return fh1.newton({xi:h2/0.90, err:12}).round(6);
}




/**
 * Calcula la entalpía de la interfaz gas líquido en una torre de enfriamiento, dadas la entalpía del gas (Hy), la temperatura del líquido (TL), la presión de operación (P) y el valor de la relación hl*a/MB*KG*a*P (hMKP).
 * @param {Object} datos - Entalpía, temperatura y presión del aire y valor de la relación hLa/MB*Kga*P.
 * @param {number} datos.Hy - Entalpía del gas en kJ/kg.
 * @param {number} datos.TL -Temperatura del agua en K.
 * @param {number} datos.P - Presión de operación en mmHg.
 * @param {number} datos.hMKP - Valor de la relación hl*a/MB*KG*a*P en kj/kg*K.
 */
function entIntGasLiq({Hy, TL, P, hMKP}={}) {
    let Hyi;
    const fti = T => {
        Hyi = entGas({T, h:humSat({T, P})});
        return Hy-hMKP*(T-TL)-Hyi; 
    }
    const Ti = fti.newton({xi:TL*0.99, err:9});
    return {Hyi: Hyi.round(6), Ti: Ti.round(2)};
}




/**
 * Calcula la altura de una Torre de Enfriamiento de agua con Coeficientes Individuales. La función recibe la velocidad superficial del líquido (L), la temperatura inicial y final del líquido (TL1 y TL2), la velocidad superficial del aire (G), la temperatura inicial del aire (TG1), la humedad inicial del aire (h1), la presión de operación (P), la relación hla/MB.KGa.P (hMKP) y el coeficiente volumétrico de transferencia de masa (KGa).<br>
 * Devuelve un objeto con la altura de la torre (z), la altura equivalente a una unidad de transferencia (HOG) y el número de unidades de transferencia (NOG).
 * @param {Object} datos - Datos necesarios para el cálculo de la altura de la torre.
 * @param {number} datos.L - Velocidad superficial del agua en m/s.m2.
 * @param {number} datos.TL1 -Temperatura del agua a la salida de la torre en K.
 * @param {number} datos.TL2 - Temperatura del agua a la entrada a la torre en K.
 * @param {number} datos.G - Velocidad superficial del aire en m/s.m2.
 * @param {number} datos.TG1 - Temperatura del aire a la entrada a la torre en K.
 * @param {number} datos.h1 - Humedad del aire a la entrada a la torre.
 * @param {number} datos.P - Presión de operación en mmHg.
 * @param {number} datos.hMKP - Relación hla/M_B.KGa.P) en kJ/kg.K.
 * @param {number} datos.KGa - Coeficiente volumétrico de transferencia de masa en kmol/s.m^3.mmHg.
 * @return {Object} Atura de la torre (z en m), altura equivalente a una unidad de transferencia (HOG en m) y número de unidades de transferencia (NOG).
 */
function alturaTECI({L, TL1, TL2, G, TG1, h1, P, hMKP, KGa}={}) {
    const Hyi = (Hy, TL) => {
        let Hi, Pp, h;
        const fti = T => {
            Pp = antoine({T});
            h = humedad({Pp, P});
            Hi = entGas({T, h});
            return Hy-hMKP*(T-TL)-Hi; 
        }
        fti.newton({xi:TL*0.99, err:9});
        return Hi.precision(6);
    };
    const fhy = Hy => {
        const TL = G/(L*CL)*(Hy-Hy1)+TL1;
        const Hi = Hyi(Hy, TL);
        return 1/(Hi-Hy);
    };
    const MB = 28.97, CL = 4.187;
    const Hy1 = entGas({T:TG1, h:h1});
    const Hy2 = L*CL/G*(TL2-TL1)+Hy1;
    let NOG = fhy.simpson({a:Hy1, b:Hy2, n:10});
    let HOG = G/(MB*KGa*P);
    let z = NOG*HOG;
    return {z:z.round(2), HOG:HOG.round(2), NOG:NOG.round(2)};
}




/**
 * Calcula la altura de una torre de enfriamiento de agua empleando Coeficientes Globales de transferencia. La función recibe la velocidad superficial del líquido (L), la temperatura inicial y final del líquido (TL1 y TL2), la velocidad superficial del aire (G), la temperatura inicial del aire (TG1), la humedad inicial del aire (h1), la presión de operación (P) y el coeficiente volumétrico de transferencia de masa (KGa).<br>
 * Devuelve un objeto con la altura de la torre (z), la altura equivalente a una unidad de transferencia (HOG) y el número de unidades de transferencia (NOG).
* @param {Object} datos - Datos necesarios para el cálculo de la altura de la torre.
 * @param {number} datos.L - Velocidad superficial del agua en m/s.m2.
 * @param {number} datos.TL1 -Temperatura del agua a la salida de la torre en K.
 * @param {number} datos.TL2 - Temperatura del agua a la entrada a la torre en K.
 * @param {number} datos.G - Velocidad superficial del aire en m/s.m2.
 * @param {number} datos.TG1 - Temperatura del aire a la entrada a la torre en K.
 * @param {number} datos.h1 - Humedad del aire a la entrada a la torre.
 * @param {number} datos.P - Presión de operación en mmHg.
 * @param {number} datos.hMKP - Relación hla/M_B.KGa.P) en kJ/kg.K.
 * @param {number} datos.KGa - Coeficiente volumétrico de transferencia de masa en kmol/s.m^3.mmHg.
 * @return {Object} Atura de la torre (z en m), altura equivalente a una unidad de transferencia (HOG en m) y número de unidades de transferencia (NOG).
*/
function alturaTECG({L, TL1, TL2, G, TG1, h1, P, KGa}={}) {
    const fhy = Hy => {
        const T = G/(L*CL)*(Hy-Hy1)+TL1;
        const h = humSat({T, P});
        const Hyi = entGas({T, h});
        return 1/(Hyi-Hy);
    };
    const MB = 28.97, CL = 4.187;
    const Hy1 = entGas({T:TG1, h:h1});
    const Hy2 = L*CL/G*(TL2-TL1)+Hy1;
    let NOG = fhy.simpson({a:Hy1, b:Hy2, n:10});
    let HOG = G/(MB*KGa*P);
    let z = NOG*HOG;
    return {z:z.round(2), HOG:HOG.round(2), NOG:NOG.round(2)};
};



/**
 * Calcula la altura de una torre de enfriamiento de agua empleando Coeficientes Globales de transferencia de masa. La función recibe la velocidad superficial del líquido (L), la temperatura inicial del líquido (TL2), la velocidad superficial del aire (G), la temperatura inicial del aire (TG1), la humedad inicial y final del aire (h1 y h2), la presión de operación (P) y el coeficiente volumétrico de transferencia de masa (KGa).<br>
 * Devuelve un objeto con la altura de la torre (z), la altura equivalente a una unidad de transferencia (HOG) y el número de unidades de transferencia (NOG).
 * @param {Object} datos - Datos necesarios para el cálculo de la altura.
 * @param {number} datos.L - Velocidad superficial del agua (m/(s.m2))
 * @param {number} datos.TL2 - Temperatura del agua a la entrada de la torre (K)
 * @param {number} datos.G - Velocidad superficial del aire (m/(s.m2))
 * @param {number} datos.TG1 - Temperatura del aire a la entrada de la torre (K)
 * @param {number} datos.h1 - Humedad del aire a la entrada a la torre
 * @param {number} datos.h2 - Humedad del aire a la salida de la torre
 * @param {number} datos.P - Presión de operación (mmHg)
 * @param {number} datos.KGa - Coeficiente volumétrico global de transferencia de masa (kmol/s.m^3.mmHg)
 * @return {Object} Atura de la torre (z en m), altura equivalente a una unidad de transferencia (HOG en m) y número de unidades de transferencia (NOG).
*/
function alturaTECGh2({L, TL2, G, TG1, h1, h2, P, KGa}={}) {
    const fhy = Hy => {
        const T = G/(L*CL)*(Hy-Hy2)+TL2;
        const h = humSat({T, P});
        const Hyi = entGas({T, h});
        return 1/(Hy-Hyi);
    };
    const MB = 28.97, CL = 4.187;
    const Hy1 = entgas({T:TG1, h:h1});
    const Pp = prePar({h:h2, P});
    const TG2 = antoine({P:Pp});
    const Hy2 = entGas({T:TG2, h:h2})
    let NOG = fhy.simpson({a:Hy2, b:Hy1, n:10});
    let HOG = G/(MB*KGa*P);
    let z = NOG*HOG;
    return {z:z.round(2), HOG:HOG.round(2), NOG:NOG.round(2)};
};




/**
 * Para una torre de enfriamiento dada, calcula el coeficiente volumétrico global de transferencia de masa (KGa, en kmol/s.m3.mmHg) y la altura equivalente a una unidad de transferencia (HOG, en m). La función recibe la velocidad superficial del líquido (L), la temperatura inicial y final del líquido (TL1 y TL2), la velocidad superficial del aire (G), la temperatura inicial del aire (TG1), la humedad inicial del aire (h1), la presión de operación (P) y la altura de la torre (z) en m.<br>
 * Devuelve un objeto con el coeficiente volumétrico de transferencia de masa (KGa) en kmol/s.m3.mmHg y la altura equivalente a una unidad de transferencia (HOG) en m.
 * @param {Object} datos - Parámetros nombrados
 * @param {number} datos.L - Velocidad superficial del agua (m/(s.m2))
 * @param {number} datos.TL1 - Temperatura del agua a la salida de la torre (K)
 * @param {number} datos.TL2 - Temperatura del agua a la entrada de la torre (K)
 * @param {number} datos.G - Velocidad superficial del aire (m/(s.m2))
 * @param {number} datos.TG1 - Temperatura del aire a la entrada de la torre (K)
 * @param {number} datos.h1 - Humedad del aire a la entrada de la torre
 * @param {number} datos.P - Presión de operación (mmHg)
 * @param {number} datos.z - Altura de la torre de enfriamiento (m)
 * @return {Object} Coeficiente global de transferencia de masa (KGa) en kmol/s.m3.mmHg y altura equivalente a una unidad de transferencia (HOG) en m.
*/
function coeficienteGlobalTE({L, TL1, TL2, G, TG1, h1, P, z}={}) {
    const fhy = Hy => {
        const T = G/(L*CL)*(Hy-Hy1)+TL1;
        const h = humSat({T, P});
        const Hyi = entGas({T, h});
        return 1/(Hyi-Hy);
    };
    const MB = 28.97, CL = 4.187;
    const Hy1 = entGas({T:TG1, h:h1});
    const Hy2 = L*CL/G*(TL2-TL1)+Hy1;
    let NOG = fhy.simpson({a:Hy1, b:Hy2, n:10});
    let HOG = z/NOG;
    let KGa = G/(MB*HOG*P);
    return {KGa:KGa.precision(4), HOG:HOG.round(2)};
};




/**
 * Dada la altura de una torre de enfriamiento, calcula el coeficiente volumétrico de transferencia de mas (KGa en kmol/s*m3*mmHg) y la altura equivalente a una unidad de transferencia (HOG en m). La función recibe la velocidad superficial del líquido (L en kg/s*m2), la temperatura de entrada del líquido (TL2 en K), la velocidad superficial del aire (G en kg/s*m2), la temperatura inicial del gas (TG! en K), las humedades inicial y final del gas (h1 y h2), la presión de operación (P en mmHg) y la altura de la torre (z en m).<br>
 * Devuelve un objeto con el coeficiente volumétrico de transferencia de masa (KGa) y la altura equivalente a una unidad de transferencia (HOG).
 * @param {Object} datos - Datos de operación para la torre de enfriamiento.
 * @param {number} datos.L - Velocidad superficial del agua en kg/s.m2.
 * @param {number} datos.TL2 - Temperatura del agua a la entrada de la torre en K.
 * @param {number} datos.G - Velocidad superficial del aire en m/s.m2.
 * @param {number} datos.TG1 - Temperatura del aire a la entrada de la torre en K.
 * @param {number} datos.h1 - Humedad del aire a la entrada de la torre.
 * @param {number} datos.h2 - Humedad del aire a la salida de la torre.
 * @param {number} datos.P - Presión de operación en mmHg.
 * @param {number} datos.z - Altura de la torre de enfriamiento en m.
 * @returns {Object} Objeto con el coeficiente volumétrico de transferencia de masa (KGa en kmol/s*m3*mmHg) y la altura equivalente a una unidad de transferencia (HOG en m).
 */
function coeficienteGlobalTEh2({L, TL2, G, TG1, h1, h2, P, z}={}) {
    const fhy = Hy => {
        const T = G/(L*CL)*(Hy-Hy2)+TL2;
        const h = humSat({T, P})
        const Hyi = entGas({T, h});
        return 1/(Hyi-Hy);
    };
    const MB = 28.97, CL = 4.187;
    const Hy1 = entGas({T:TG1, h:h1});
    const Pp = prePar({h:h2, P});
    const TG2 = antoine({P:Pp});
    const Hy2 = entGas({T:TG2, h:h2});
    let NOG = fhy.simpson({a:Hy2, b:Hy1, n:10});
    let HOG = z/NOG;
    let KGa = G/(MB*HOG*P);
    return {KGa:KGa.precision(4), HOG:HOG.round(2)};
};




/**
 * Calcula el flujo superficial mínimo de aire en kg/s.m2, en una torre de enfriamiento que enfría agua de TL2 a TL1.
 * @param {Object} datos - Parámetros nombrados
 * @param {number} datos.L - Velocidad superficial del agua (m/(s.m2))
 * @param {number} datos.TL1 - Temperatura del agua a la salida de la torre (K)
 * @param {number} datos.TL2 - Temperatura del agua a la entrada de la torre (K)
 * @param {number} datos.TG1 - Temperatura del aire a la entrada de la torre (K)
 * @param {number} datos.h1 - Humedad del aire a la entrada de la torre
 * @param {number} datos.P - Presión de operación (mmHg)
 * @returns {number} Flujo superficial mínimo del aire, en kg/s.m2, para lograr el enfriamiento requerido.
 */
var Gmin = ({L, TL1, TL2, TG1, h1, P}) => {
    const CL = 4.187;
    const Hy1 = entGas({T:TG1, h:h1});
    const h2 = humSat({T:TL2, P});
    const Hy2 = entGas({T:TL2, h:h2});
    return L*CL*(TL2-TL1)/(Hy2-Hy1);
};




/**
 * Calcula el flujo superficial mínimo de agua (L en kg/s.m2), para secar aire desde una humedad inicial h1, hasta que sale saturado a la temperatura de entrada del agua (TL2). La función recibe la temperatura de entrada del líquido (TL2 en K), la velocidad superficial del aire (G en kg/s*m2), la temperatura inicial del aire (TG1 en K), la humedad inicial del aire (h1) y la presión de operación (P en mmHG).<br>
 * Devuelve el flujo mínimo de agua necesario para lograr el secado requerido.
 * @param {Object} datos - Datos de operación para el secado de aire.
 * @param {number} datos.TL2 - Temperatura del agua a la entrada de la torre en K.
 * @param {number} datos.G - Velocidad superficial del aire en m/s.m2.
 * @param {number} datos.TG1 - Temperatura del aire a la entrada de la torre en K.
 * @param {number} datos.h1 - Humedad del aire a la entrada de la torre.
 * @param {number} datos.P - Presión de operación en mmHg.
 * @returns {number} Velocidad superficial del líquido en kg/s*m2.
 */
function Lmin({TL2, G, TG1, h1, P}={}) {
    const CL = 4.187;
    const Hy1 = entGas({T:TG1, h:h1});
    const Hy2 = entGas({T:TL2, h:humSat({T:TL2,P})});
    return G/CL*(Hy1-Hy2)/(TL1-TL2);
};




/**
 * Calcula, para una etapa de absorción, los flujos y fracciones molares finales de las fases gaseosa y líquida (L1, x1, V1 y y1).<br>La función recibe los flujos y fracciones molares iniciales de las fases gaseosa y líquida (L0, x0, V0, y0), así como las funciones de equilibrio (fx y fy).
 * @param {Object} datos - L0, x0, V0, y0, fx y fy.
 * @param {number} datos.L0 - Flujo molar inicial de líquido en kmol/h.
 * @param {number} datos.x0 - Fracción molar inicial de soluto en el líquido.
 * @param {number} datos.V0 - Flujo molar inicial de gas en kmol/h.
 * @param {number} datos.y0 - Fracción molar inicial de soluto en el gas.
 * @param {number} datos.fx - Función que devuelve "y" para un valor dado de "x".
 * @param {number} datos.fy - Función que devuelve "x" para un valor dado de "y".
 * @returns {Object} Objeto con el flujo molar final del líquido (L1 en kmol/h), la fracción molar final de soluto en el líquido (x1), el flujo molar final de gas (V1 en kmol/h) y la fracción molar final de soluto en el gas (y1).
 */
function absor1l0x0v0y0({L0,x0,V0,y0,fx,fy}={}) {
    const M=L0+V0;
    const B=L0*(1-x0);
    const A=L0*x0+V0*y0;
    let y1, L1, V1;
    function fx1(x1) {
      y1=fx(x1);
      L1=B/(1-x1);
      V1=M-L1;
      return L1*x1+V1*y1-A;
    }
    const x1=fx1.refa({xi:0, xf:fy(y0)});
    return {L1, x1, V1, y1};
}




/**
 * Calcula, para una etapa de absorción, los flujos molares inicial y final de la fase líquida, la fracción molar final de la fase líquida y el flujo y fracción  molar finales de la fase gaseosa. <br>La función recibe la fracción molar inicial en la fase líquida (x0), el flujo y fracción molar iniciales en la fase gaseosa (V0, y0), el porcentaje de absorción (pa) y la función de equilibrio dependiente "y" (fy).
 * @param {Object} datos - x0, V0, y0, pa y fy.
 * @param {number} datos.x0 - Fracción molar inicial de soluto en el líquido.
 * @param {number} datos.V0 - Flujo molar inicial de gas en kmol/h.
 * @param {number} datos.y0 - Fracción molar inicial de soluto en el gas.
 * @param {number} datos.pa - Porcentaje de soluto que se absorbe.
 * @param {number} datos.fy - Función que devuelve "x" para un valor dado de "y".
 * @returns {Object} Objeto con los flujos molares inicial y final del líquido (L0 y L1 en kmol/h), la fracción molar final de soluto en el líquido (x1), el flujo molar final de gas (V1 en kmol/h) y la fracción molar final de soluto en el gas (y1).
 */
function absor1x0v0y0pa({x0,V0,y0,pa,fy}={}) {
    const V1y1 = (1-pa/100)*V0*y0;
    const C = V0*(1-y0);
    const V1 = C+V1y1;
    const y1 = V1y1/V1;
    const x1 = fy(y1);
    const D = V0-V1;
    const L0 = (1-x1)/(x1-x0)*D;
    const L1 = L0*(1-x0)/(1-x1);
    return {L0, L1, x1, V1, y1};
}




/**
 * Calcula, para una etapa de absorción, los flujos molares inicial y final de la fase líquida (L0, L1), el flujo molar y la fracción molar final en de la fase gaseosa (Vi, y1). <br>Recibe las fracciones molares inicial y final en la fase líquida (x0, x1), el flujo molar y la fracción molar iniciales en la fase gaseosa (V0, y0) y la función de equilibrio dependiente de "x" (fx).
 * @param {Object} datos - x0, x1, V0, y0, fx.
 * @param {number} datos.L0 - Flujo molar inicial de líquido en kmol/h.
 * @param {number} datos.x0 - Fracción molar inicial de soluto en el líquido.
 * @param {number} datos.V0 - Flujo molar inicial de gas en kmol/h.
 * @param {number} datos.y0 - Fracción molar inicial de soluto en el gas.
 * @param {number} datos.fx - Función que devuelve "y" para un valor dado de "x".
 * @returns {Object} Objeto con los flujos molares inicial y final del líquido (L0 y L1 en kmol/h), el flujo molar final de gas (V1 en kmol/h) y la fracción molar final de soluto en el gas (y1).
 */
function absor1x0x1v0y0({x0, x1, V0, y0, fx}={}) {
    const y1 = fx(x1);
    const C = V0*(1-y0);
    const V1 = C/(1-y1);
    const DV = V1-V0;
    const L0 = DV*(1-x1)/(x0-x1);
    const L1 = L0-DV;
    return {L0, L1, V1, y1};
}




/**
 * Calcula, para una torre con múltiples etapas de absorción, los flujos y fracciones molares finales de las fases gaseosa y líquida (Ln, xn, V1 y y1). <br>La función recibe los flujos y fracciones molares iniciales de las fases gaseosa y líquida (L0, x0, V0, y0) y el porcentaje de absorción (pa).
 * @param {Object} datos - L0, x0, V0, y0 y pa.
 * @param {number} datos.L0 - Flujo molar inicial de líquido en kmol/h.
 * @param {number} datos.x0 - Fracción molar inicial de soluto en el líquido.
 * @param {number} datos.V0 - Flujo molar inicial de gas en kmol/h.
 * @param {number} datos.y0 - Fracción molar inicial de soluto en el gas.
 * @param {number} datos.pa - Porcentaje de soluto que se absorbe.
 * @returns {Object} Objeto con el flujo molar final del líquido (Ln en kmol/h), la fracción molar final de soluto en el líquido (xn), el flujo molar final de gas (V1 en kmol/h) y la fracción molar final de soluto en el gas (y1).
 */
function absorNl0x0v0y0pa({L0, x0, V0, y0, pa}={}) {
    const M = L0+V0;
    const C = V0*(1-y0);
    const A = L0*x0+V0*y0;
    const V1y1 = (1-pa/100)*V0*y0;
    const V1 = C+V1y1;
    const y1 = V1y1/V1;
    const Ln = M-V1;
    const xn = (A-V1y1)/Ln;
    return {Ln, xn, V1, y1};
}




/**
 * Calcula, para una torre con múltiples etapas de absorción, el número de etapas de equilibrio necesarias para lograr la separación requerida (n). <br>La función recibe el flujo molar inicial de líquido (L0), las fracciones molares iniciales y final de la fase líquida (x0, xn), el flujo y fracción molar finales de la fase gaseosa (V1, y1) y la función de equilibrio dependiente de "y" (fy).
 * @param {Object} datos - L0, x0, xn, V1, y1 y fy.
 * @param {number} datos.L0 - Flujo molar inicial de líquido en kmol/h.
 * @param {number} datos.x0 - Fracción molar inicial de soluto en el líquido.
 * @param {number} datos.xn - Fracción molar final de soluto en el líquido.
 * @param {number} datos.V1 - Flujo molar final de gas en kmol/h.
 * @param {number} datos.y1 - Fracción molar final de soluto en el gas.
 * @param {number} datos.fy - Función que devuelve "x" para un valor dado de "y".
 * @returns {number} Número de etapas teóricas de equilibrio necesarias para lograr la separación requerida.
 */
function absorNl0x0xnv1y1({L0, x0, xn, V1, y1, fy}={}) {
    const DM = L0-V1;
    const DA = L0*x0-V1*y1;
    const B = L0*(1-x0);
    let i = 0, xi, Li, Vj, yi;
    for(;;) {
      xi = fy(y1);
      if (xi>=xn) break;
      x0 = xi;
      i++;
      if (i>500) return (i.round(2));
      Li = B/(1-xi);
      Vj = Li-DM;
      y1 = (Li*xi-DA)/Vj;
    }
    const n = i+(xn-x0)/(xi-x0);
    return n.round(2);
}




/**
 * Calcula, para una torre de absorción, el flujo mínimo de líquido para lograr la separación requerida (número infinito de etapas). <br>La función recibe la fracción inicial en la fase líquida (x0), el flujo y fracción molar iniciales de la fase gaseosa (V0, y0), el porcentaje de absorción (pa) y  la función de equilibrio dependiente de "y" (fy).
 * @param {Object} datos - x0, V0, y0, pa y fy.
 * @param {number} datos.x0 - Fracción molar inicial de soluto en el líquido.
 * @param {number} datos.V0 - Flujo molar inicial de gas en kmol/h.
 * @param {number} datos.y0 - Fracción molar inicial de soluto en el gas.
 * @param {number} datos.pa - Porcentaje de soluto que se absorbe.
 * @param {number} datos.fy - Función que devuelve "x" para un valor dado de "y".
 * @returns {number} Flujo mínimo de líquido, en kmol/h, con el que se puede lograr la separación requerida.
 */
function absorLmin({x0, V0, y0, pa, fy}={}) {
    const xmax = fy(y0);
    const V1y1 = (1-pa/100)*V0*y0;
    const Bmin = (V0*y0-V1y1)/(xmax/(1-xmax)-x0/(1-x0));
    return Bmin*(1-x0);
}




/**
 * Calcula, empleando el método de Kremser, el número de etapas para lograr la separación requerida en una torre de absorción. <br>La función recibe el flujo y  fracción molar iniciales del líquido (L0, x0), el flujo molar final del líquido (Ln), el flujo y fracción molar iniciales de la fase gaseosa (V0, y0), el flujo y fracción molar finales de la fase gaseosa (V1, y1) y la constante de la ley de Henry (H).
 * @param {Object} datos - Lo, x0, Ln, V0, y0, V1, y1 y H.
 * @param {number} datos.L0 - Flujo molar inicial de líquido en kmol/h.
 * @param {number} datos.x0 - Fracción molar inicial de soluto en el líquido.
 * @param {number} datos.Ln - Flujo molar final de líquido en kmol/h.
 * @param {number} datos.V0 - Flujo molar inicial de gas en kmol/h.
 * @param {number} datos.y0 - Fracción molar inicial de soluto en el gas.
 * @param {number} datos.V1 - Flujo molar final de gas en kmol/h.
 * @param {number} datos.y1 - Fracción molar final de soluto en el gas.
 * @param {number} datos.H - Constante de la ley de Henry.
 * @returns {number} Número de etapas teóricas de equilibrio necesarias para lograr la separación requerida.
 */
function absorKremser({L0, x0, Ln, V0, y0, V1, y1, H}={}) {
    const A = Math.sqrt((L0/(H*V1)*(Ln/(H*V0))));
    return Math.log((y0-H*x0)/(y1-H*x0)*(1-1/A)+1/A)/Math.log(A);
}





/**
 * Calcula, para una torre de absorción rellena, los flujos y fracciones molares finales de las fases líquida (L1, x1) y gaseosa (V2, y2).<br>La función recibe el flujo y la fracción molar iniciales en la fase líquida (L2, x2), el flujo y fracción molar iniciales de la fase gaseosa (V1, y1) y el porcentaje de absorción (pa) de soluto.
 * @param {Object} datos - L2, x2, V1, y1, pa.
 * @param {number} datos.L2 - Flujo molar inicial de líquido en kmol/h.
 * @param {number} datos.x2 - Fracción molar inicial de soluto en el líquido.
 * @param {number} datos.V1 - Flujo molar inicial de gas en kmol/h.
 * @param {number} datos.y1 - Fracción molar inicial de soluto en el gas.
 * @param {number} datos.pa - Porcentaje de absorción del soluto.
 * @returns {Object} Objecto con los flujos (en kmol/h) y fracciones molares finales de las fases líquida (L1, x1) y gaseosa (V2, y2).
 */
function abrc1({L2, x2, V1, y1, pa}={}) {
    const B = L2*(1-x2);
    const C = V1*(1-y1);
    const A = B*(x2/(1-x2))+C*(y1/(1-y1));
    const V2y2 = V1*y1*(1-pa/100);
    const V2 = C+V2y2;
    const y2 = V2y2/V2;
    const D = (A-C*(y2/(1-y2)))/B;
    const x1 = D/(1+D);
    const L1 = B/(1-x1);
    return {L1, x1, V2, y2};
}





/**
 * Calcula, para una torre de absorción rellena, las fracciones molares (xi, yi) en la interfaz gas líquido y los flujos molares superficiales de soluto en las fases líquida y gaseosa (Nax, Nay).<br>La función recibe las fracciones molares de soluto en las fases gaseosa y líquida (x, y), los coeficientes individuales ideales de las fases líquida y gaseosa (kx y ky) y las funciones de equilibrio dependientes de las fracciones molares en las fases líquida y gaseosa (fx, fy).
 * @param {Object} datos - x, y, kx, ky, fx, fy.
 * @param {number} datos.x - Fracción molar de soluto en la fase líquida.
 * @param {number} datos.y - Fracción molar de soluto en la fase gaseosa.
 * @param {number} datos.kx - Coeficiente individual ideal en la fase líquida en kmol/s*m2.
 * @param {number} datos.ky - Coeficiente individual ideal en la fase gaseosa en kmol/s*m2.
 * @param {number} datos.fx - Función de equilibrio dependiente de la fracción molar en la fase líquida.
 * @param {number} datos.fy - Función de equilibrio dependiente de la fracción molar en la fase gaseosa.
 * @returns {Object} Objecto con las fracciones molares en la interfaz gas-líquido (xi, yi) y los flujos molares superficiales de soluto (en kmol/s*m2) en las fases líquida y gaseosas (Nax, Nay).
 */
function xiyi1({x, y, kpx, kpy, fx, fy}={}) {
    let yi;
    function f(xi) {
      yi = fx(xi);
      return kpx/kpy+(y-yi)/(x-xi);
    }
    const xi = f.refa({xi:x+1e-5, xf: fy(y)-1e-5});
    const Nax = kpx*(xi-x);
    const Nay = kpy*(y-yi);
    return {xi, yi, Nax, Nay}
}




/**
 * Calcula, para una torre de absorción rellena, las fracciones molares (xi, yi) en la interfaz gas líquido, los coeficientes individuales reales (kx, ky) las medias logarítmicas de xi y yi (xim, yim) y los flujos molares superficiales de soluto en las fases líquida y gaseosa (Nax, Nay).<br>La función recibe las fracciones molares de soluto en las fases gaseosa y líquida (x, y), los coeficientes individuales ideales de las fases líquida y gaseosa (kpx y kpy) y las funciones de equilibrio dependientes de las fracciones molares en las fases líquida y gaseosa (fx, fy).
 * @param {Object} datos - x, y, kpx, kpy, fx, fy.
 * @param {number} datos.x - Fracción molar de soluto en la fase líquida.
 * @param {number} datos.y - Fracción molar de soluto en la fase gaseosa.
 * @param {number} datos.kpx - Coeficiente individual ideal en la fase líquida en kmol/s*m2.
 * @param {number} datos.kpy - Coeficiente individual ideal en la fase gaseosa en kmol/s*m2.
 * @param {number} datos.fx - Función de equilibrio dependiente de la fracción molar en la fase líquida.
 * @param {number} datos.fy - Función de equilibrio dependiente de la fracción molar en la fase gaseosa.
 * @returns {Object} Objecto con las fracciones molares en la interfaz gas-líquido (xi, yi), las medias logarítmicas de x-xi (xim) y yi-y (yim), los coeficientes individuales reales (kx, ky) y los flujos molares superficiales de soluto (en kmol/s*m2) en las fases líquida y gaseosas (Nax, Nay).
 */
function xiyi2({x, y, kpx, kpy, fx, fy}={}) {

    // Función de la fracción molar en el líquido en la interfaz gas-líquido
    let yi, xim, yim, kx, ky;
    function fxi(xi) {
        yi = fx(xi);
        xim = ((1-x)-(1-xi))/Math.log((1-x)/(1-xi));
        yim = ((1-yi)-(1-y))/Math.log((1-yi)/(1-y));
        kx = kpx/xim;
        ky = kpy/yim;
        return kx/ky+(y-yi)/(x-xi);
    }
    
    // Resolución de la función "fxi"
    const xi = fxi.refa({xi:x+1e-5, xf:fy(y)-1e-5});
    
    // Flujos de soluto
    const Nax = kx*(xi-x);
    const Nay = ky*(y-yi);
    
    return {xi, yi, xim, yim, kx, ky, Nax, Nay};
    
}




/**
 * Calcula el coeficiente global ideal de transferencia de masa, para la absorción en torres rellenas, basado en la fase líquida.<br>La función recibe la fracción molar en la fase gaseosa (y), las fracciones molares en la interfaz gas-líquido (xi, yi), los coeficientes individuales ideales de transferencia de masa (kpx, kpy) y la fracción molar en la fase líquida (xa) en equilibrio con la fracción molar en la fase gaseosa.
 * @param {Object} datos - y, xi, yi, kpx, kpy, xa
 * @param {number} datos.y - Fracción molar de soluto en la fase gaseosa.
 * @param {number} datos.xi - Fracción molar de soluto en el líquido, en la interfaz gas-líquido.
 * @param {number} datos.yi - Fracción molar de soluto en el gas, en la interfaz gas-líquido.
 * @param {number} datos.kpx - Coeficiente individual ideal en la fase líquida en kmol/s*m2.
 * @param {number} datos.kpy - Coeficiente individual ideal en la fase gaseosa en kmol/s*m2.
 * @param {number} datos.xa - Fracción molar de soluto en la fase líquida en equilibrio con la fracción molar en la fase gaseosa (y).
 * @returns {number} Coeficiente global ideal de transferencia de masa, basado en la fase líquida.
 */
function Kx1({y, xi, yi, kpx, kpy, xa}={}) {
    const m = (y-yi)/(xa-xi);
    return 1/(1/(m*kpy)+1/kpx);
}




/**
 * Calcula el coeficiente global ideal de transferencia de masa, para la absorción en torres rellenas, basado en la fase gaseosa.<br>La función recibe la fracción molar en la fase líquida (x), las fracciones molares en la interfaz gas-líquido (xi, yi), los coeficientes individuales ideales de transferencia de masa (kpx, kpy) y la fracción molar en la fase gaseosa (ya) en equilibrio con la fracción molar en la fase líquida.
 * @param {Object} datos - x, xi, yi, kpx, kpy, ya
 * @param {number} datos.x - Fracción molar de soluto en la fase líquida.
 * @param {number} datos.xi - Fracción molar de soluto en el líquido, en la interfaz gas-líquido.
 * @param {number} datos.yi - Fracción molar de soluto en el gas, en la interfaz gas-líquido.
 * @param {number} datos.kpx - Coeficiente individual ideal en la fase líquida en kmol/s*m2.
 * @param {number} datos.kpy - Coeficiente individual ideal en la fase gaseosa en kmol/s*m2.
 * @param {number} datos.ya - Fracción molar de soluto en la fase gaseosa en equilibrio con la fracción molar en la fase líquida (x).
 * @returns {number} Coeficiente global ideal de transferencia de masa, basado en la fase gaseosa.
 */
function Ky1({x, xi, yi, kpx, kpy, ya}={}) {
    const m = (yi-ya)/(xi-x);
    return 1/(1/kpy+m/kpx);
}




/**
 * Calcula el coeficiente global real de transferencia de masa, para la absorción en torres rellenas, basado en la fase líquida.<br>La función recibe la fracción molar en la fase líquida (x), el coeficiente global ideal basado en la fase líquida (Kpx) y la fracción molar en la fase líquida (xa) en equilibrio con la fracción molar en la fase gaseosa.
 * @param {Object} datos - x, Kpx, xa
 * @param {number} datos.x - Fracción molar de soluto en la fase líquida.
 * @param {number} datos.Kpx - Coeficiente global ideal en la fase líquida en kmol/s*m2.
 * @param {number} datos.xa - Fracción molar de soluto en la fase líquida en equilibrio con la fracción molar en la fase gaseosa (y).
 * @returns {number} Coeficiente global real de transferencia de masa, basado en la fase líquida.
 */
function Kx2({x, Kpx, xa}={}) {
    const ml = ((1-x)-(1-xa))/Math.log((1-x)/(1-xa));
    return Kpx/ml
}




/**
 * Calcula el coeficiente global real de transferencia de masa, para la absorción en torres rellenas, basado en la fase gaseosa.<br>La función recibe la fracción molar en la fase gaseosa (y), el coeficiente global ideal basado en la fase gaseosa (Kpy) y la fracción molar en la fase gaseosa (ya) en equilibrio con la fracción molar en la fase líquida.
 * @param {Object} datos - y, Kpy, ya
 * @param {number} datos.y - Fracción molar de soluto en la fase gaseosa.
 * @param {number} datos.Kpy - Coeficiente global ideal en la fase gaseosa en kmol/s*m2.
 * @param {number} datos.ya - Fracción molar de soluto en la fase gaseosa en equilibrio con la fracción molar en la fase gaseosa (y).
 * @returns {number} Coeficiente global real de transferencia de masa, basado en la fase gaseosa.
 */
function Ky2({y, Kpy, ya}={}) {
    const ml = ((1-ya)-(1-y))/Math.log((1-ya)/(1-y));
    return Kpy/ml;
}



/**
 * Calcula la altura de una torre de absorción rellena (mezclas gaseosas diluidas) basado en la fase líquida, con coeficientes individuales ideales.<br>La función recibe los flujos y fracciones molares iniciales y finales del líquido (L2, x2, L1, x1), las fracciones molares inicial y final en el gas (y1, y2), los coeficientes individuales ideales de las fases líquida y gaseosa (kpx y kpy), las funciones de equilibrio dependientes de las fracciones molares en las fases líquida y gaseosa (fx, fy) y el área transversal de la columna (S).
 * @param {Object} datos - L2, x2, L1, x1, y1, y2, kpx, kpy, fx, fy, S.
 * @param {number} datos.L2 - Flujo molar inicial de la fase líquida.
 * @param {number} datos.x2 - Fracción molar inicial de soluto en la fase líquida.
 * @param {number} datos.L1 - Flujo molar final de la fase líquida.
 * @param {number} datos.x1 - Fracción molar final de soluto en la fase líquida.
 * @param {number} datos.y1 - Fracción molar inicial de soluto en la fase líquida.
 * @param {number} datos.y2 - Fracción molar final de soluto en la fase líquida.
 * @param {number} datos.kpx - Coeficiente individual ideal en la fase líquida en kmol/s*m2.
 * @param {number} datos.kpy - Coeficiente individual ideal en la fase gaseosa en kmol/s*m2.
 * @param {number} datos.fx - Función de equilibrio dependiente de la fracción molar en la fase líquida.
 * @param {number} datos.fy - Función de equilibrio dependiente de la fracción molar en la fase gaseosa.
 * @param {number} datos.S - Area transversal de la columna (m2).
 * @returns {Object} Objecto con la altura equivalente a una unidad de transferencia basada en la fase líquida (HL), el número de unidades de transferencia basado en la fase líquida (NL) y la altura de la torre de absorción (z) en m.
 */
function zabkxi({L2, x2, L1, x1, y1, y2, kpx, kpy, fx, fy, S}={}) {
    const L = (L1+L2)/2;
    const HL = L/(kpx*S);
    const {xi:xi1, xim:xim1} = xiyi2({x:x1, y:y1, kpx, kpy, fx, fy});
    const {xi:xi2, xim:xim2} = xiyi2({x:x2, y:y2, kpx, kpy, fx, fy});
    const xmi = ((xi1-x1)-(xi2-x2))/Math.log((xi1-x1)/(xi2-x2));
    const xm = (xim1/(1-x1)+xim2/(1-x2))/2;
    const NL = xm*(x1-x2)/xmi;
    const z = HL*NL;
    return {HL, NL, z};
}



/**
 * Calcula la altura de una torre de absorción rellena (mezclas gaseosas diluidas), basado en la fase líquida, con coeficientes individuales reales.<br>La función recibe los flujos y fracciones molares iniciales y finales del líquido (L2, x2, L1, x1), las fracciones molares inicial y final en el gas (y1, y2), los coeficientes individuales ideales de las fases líquida y gaseosa (kpx y kpy), las funciones de equilibrio dependientes de las fracciones molares en las fases líquida y gaseosa (fx, fy) y el área transversal de la columna (S).
 * @param {Object} datos - L2, x2, L1, x1, y1, y2, kpx, kpy, fx, fy, S.
 * @param {number} datos.L2 - Flujo molar inicial de la fase líquida.
 * @param {number} datos.x2 - Fracción molar inicial de soluto en la fase líquida.
 * @param {number} datos.L1 - Flujo molar final de la fase líquida.
 * @param {number} datos.x1 - Fracción molar final de soluto en la fase líquida.
 * @param {number} datos.y1 - Fracción molar inicial de soluto en la fase gaseosa.
 * @param {number} datos.y2 - Fracción molar final de soluto en la fase gaseosa.
 * @param {number} datos.kpx - Coeficiente individual ideal en la fase líquida en kmol/s*m2.
 * @param {number} datos.kpy - Coeficiente individual ideal en la fase gaseosa en kmol/s*m2.
 * @param {number} datos.fx - Función de equilibrio dependiente de la fracción molar en la fase líquida.
 * @param {number} datos.fy - Función de equilibrio dependiente de la fracción molar en la fase gaseosa.
 * @param {number} datos.S - Area transversal de la columna (m2).
 * @returns {Object} Objecto con la altura equivalente a una unidad de transferencia basada en la fase líquida (HL), el número de unidades de transferencia basado en la fase líquida (NL) y la altura de la torre de absorción (z) en m.
 */
function zabkx({L2, x2, L1, x1, y1, y2, kx, ky, fx, fy, S}={}) {
    const {xi:xi1} = xiyi1({x:x1, y:y1, kpx:kx, kpy:ky, fx, fy});
    const {xi:xi2} = xiyi1({x:x2, y:y2, kpx:kx, kpy:ky, fx, fy});
    const xim1 = ((1-x1)-(1-xi1))/Math.log((1-x1)/(1-xi1));
    const xim2 = ((1-x2)-(1-xi2))/Math.log((1-x2)/(1-xi2));
    const Lm = (L1/xim1+L2/xim2)/2;
    const HL = Lm/(kx*S);
  	const xmi = ((xi1-x1)-(xi2-x2))/Math.log((xi1-x1)/(xi2-x2));
    const xm = (xim1/(1-x1)+xim2/(1-x2))/2;
    const NL = xm*(x1-x2)/xmi;
    const z = HL*NL;
    return {HL, NL, z};
}



/**
 * Calcula la altura de una torre de absorción rellena (mezclas gaseosas diluidas), basado en la fase gaseosa, con coeficientes individuales ideales.<br>La función recibe los flujos y fracciones molares iniciales y finales del gas (V1, y1, V2, y2), las fracciones molares inicial y final en el líquido (x2, x1), los coeficientes individuales ideales de las fases líquida y gaseosa (kpx y kpy), las funciones de equilibrio dependientes de las fracciones molares en las fases líquida y gaseosa (fx, fy) y el área transversal de la columna (S).
 * @param {Object} datos - V1, y1, V2, y2, x2, x1, kpx, kpy, fx, fy, S.
 * @param {number} datos.V1 - Flujo molar inicial de la fase gaseosa.
 * @param {number} datos.y1 - Fracción molar inicial de soluto en la fase gaseosa.
 * @param {number} datos.V2 - Flujo molar final de la fase gaseosa.
 * @param {number} datos.y2 - Fracción molar final de soluto en la fase gaseosa.
 * @param {number} datos.x2 - Fracción molar inicial de soluto en la fase líquida.
 * @param {number} datos.y2 - Fracción molar final de soluto en la fase gaseosa.
 * @param {number} datos.kpx - Coeficiente individual ideal en la fase líquida en kmol/s*m2.
 * @param {number} datos.kpy - Coeficiente individual ideal en la fase gaseosa en kmol/s*m2.
 * @param {number} datos.fx - Función de equilibrio dependiente de la fracción molar en la fase líquida.
 * @param {number} datos.fy - Función de equilibrio dependiente de la fracción molar en la fase gaseosa.
 * @param {number} datos.S - Area transversal de la columna (m2).
 * @returns {Object} Objecto con la altura equivalente a una unidad de transferencia basada en la fase gaseosa (HG), el número de unidades de transferencia basado en la fase gaseosa (NG) y la altura de la torre de absorción (z) en m.
 */
function zabkyi({V1, y1, V2, y2, x2, x1, kpx, kpy, fx, fy, S}={}) {
    const V = (V1+V2)/2;
    const HG = V/(kpy*S);
    const {yi:yi1, yim:yim1} = xiyi2({x:x1, y:y1, kpx, kpy, fx, fy});
    const {yi:yi2, yim:yim2} = xiyi2({x:x2, y:y2, kpx, kpy, fx, fy});
    const ymi = ((y1-yi1)-(y2-yi2))/Math.log((y1-yi1)/(y2-yi2));
    const ym = (yim1/(1-y1)+yim2/(1-y2))/2;
    const NG = ym*(y1-y2)/ymi;
    const z = HG*NG;
    return {HG, NG, z};
}



/**
 * Calcula, la altura de una torre de absorción rellena (mezclas gaseosas diluidas), basado en la fase líquida, con coeficientes individuales reales.<br>La función recibe los flujos y fracciones molares iniciales y finales del gas (V2, y1, V2, y2), las fracciones molares inicial y final en el líquido (x2, x1), los coeficientes individuales ideales de las fases líquida y gaseosa (kpx y kpy), las funciones de equilibrio dependientes de las fracciones molares en las fases líquida y gaseosa (fx, fy) y el área transversal de la columna (S).
 * @param {Object} datos - V1, y1, V2, y2, x2, x1, kpx, kpy, fx, fy, S.
 * @param {number} datos.V1 - Flujo molar inicial de la fase gaseosa.
 * @param {number} datos.y1 - Fracción molar inicial de soluto en la fase gaseosa.
 * @param {number} datos.V2 - Flujo molar final de la fase gaseosa.
 * @param {number} datos.y2 - Fracción molar final de soluto en la fase gaseosa.
 * @param {number} datos.x2 - Fracción molar inicial de soluto en la fase líquida.
 * @param {number} datos.x1 - Fracción molar final de soluto en la fase líquida.
 * @param {number} datos.kpx - Coeficiente individual ideal en la fase líquida en kmol/s*m2.
 * @param {number} datos.kpy - Coeficiente individual ideal en la fase gaseosa en kmol/s*m2.
 * @param {number} datos.fx - Función de equilibrio dependiente de la fracción molar en la fase líquida.
 * @param {number} datos.fy - Función de equilibrio dependiente de la fracción molar en la fase gaseosa.
 * @param {number} datos.S - Area transversal de la columna (m2).
 * @returns {Object} Objecto con la altura equivalente a una unidad de transferencia basada en la fase gaseosa (HG), el número de unidades de transferencia basado en la fase gaseosa (NG) y la altura de la torre de absorción (z) en m.
 */
function zabky({V1, y1, V2, y2, x2, x1, kx, ky, fx, fy, S}={}) {
    const {yi:yi1} = xiyi1({x:x1, y:y1, kpx:kx, kpy:ky, fx, fy});
    const {yi:yi2} = xiyi1({x:x2, y:y2, kpx:kx, kpy:ky, fx, fy});
    const yim1 = ((1-yi1)-(1-y1))/Math.log((1-yi1)/(1-y1));
    const yim2 = ((1-yi2)-(1-y2))/Math.log((1-yi2)/(1-y2));
  	const ymi = ((y1-yi1)-(y2-yi2))/Math.log((y1-yi1)/(y2-yi2));
    const Vm = (V1/yim1+V2/yim2)/2;
    const HG = Vm/(ky*S);
    const ym = (yim1/(1-y1)+yim2/(1-y2))/2;
    const NG = ym*(y1-y2)/ymi;
    const z = HG*NG;
    return {HG, NG, z};
}



/**
 * Calcula la altura de una torre de absorción rellena (mezclas gaseosas diluidas), basado en la fase líquida, con el coeficiente global ideal.<br>La función recibe los flujos y fracciones molares iniciales y finales del líquido (L2, x2, L1, x1), las fracciones molares inicial y final en el gas (y1, y2), el coeficiente global ideal de la fase líquida (Kpx), la función de equilibrio dependiente de la fracción molar en las fases gaseosa (fy) y el área transversal de la columna (S).
 * @param {Object} datos - L2, x2, L1, x1, y1, y2, Kpx, fy, S.
 * @param {number} datos.L2 - Flujo molar inicial de la fase líquida.
 * @param {number} datos.x2 - Fracción molar inicial de soluto en la fase líquida.
 * @param {number} datos.L1 - Flujo molar final de la fase líquida.
 * @param {number} datos.x1 - Fracción molar final de soluto en la fase líquida.
 * @param {number} datos.y1 - Fracción molar inicial de soluto en la fase gaseosa.
 * @param {number} datos.y2 - Fracción molar final de soluto en la fase gaseosa.
 * @param {number} datos.Kpx - Coeficiente global ideal en la fase líquida en kmol/s*m2.
 * @param {number} datos.fy - Función de equilibrio dependiente de la fracción molar en la fase gaseosa.
 * @param {number} datos.S - Area transversal de la columna (m2).
 * @returns {Object} Objecto con la altura equivalente a una unidad de transferencia basada en la fase líquida (HOL), el número de unidades de transferencia basado en la fase líquida (NOL) y la altura de la torre de absorción (z) en m.
 */
function zabKxi({L2, x2, L1, x1, y1, y2, Kpx, fy, S}={}) {
    const L = (L1+L2)/2;
    const HOL = L/(Kpx*S);
  	const xa1 = fy(y1);
  	const xa2 = fy(y2);
    const xam1 = ((1-x1)-(1-xa1))/Math.log((1-x1)/(1-xa1));
    const xam2 = ((1-x2)-(1-xa2))/Math.log((1-x2)/(1-xa2));
    const xma = ((xa1-x1)-(xa2-x2))/Math.log((xa1-x1)/(xa2-x2));
    const xm = (xam1/(1-x1)+xam2/(1-x2))/2;
    const NOL = xm*(x1-x2)/xma;
    const z = HOL*NOL;
    return {HOL, NOL, z};
}




/**
 * Calcula la altura de una torre de absorción rellena (mezclas gaseosas diluidas), basado en la fase líquida, con el coeficiente global real.<br>La función recibe los flujos y fracciones molares iniciales y finales del líquido (L2, x2, L1, x1), las fracciones molares inicial y final en el gas (y1, y2), el coeficiente global ideal de la fase líquida (Kpx), la función de equilibrio dependiente de la fracción molar en las fases gaseosa (fy) y el área transversal de la columna (S).
 * @param {Object} datos - L2, x2, L1, x1, y1, y2, Kx, fy, S.
 * @param {number} datos.L2 - Flujo molar inicial de la fase líquida.
 * @param {number} datos.x2 - Fracción molar inicial de soluto en la fase líquida.
 * @param {number} datos.L1 - Flujo molar final de la fase líquida.
 * @param {number} datos.x1 - Fracción molar final de soluto en la fase líquida.
 * @param {number} datos.y1 - Fracción molar inicial de soluto en la fase gaseosa.
 * @param {number} datos.y2 - Fracción molar final de soluto en la fase gaseosa.
 * @param {number} datos.Kx - Coeficiente global real en la fase líquida en kmol/s*m2.
 * @param {number} datos.fy - Función de equilibrio dependiente de la fracción molar en la fase gaseosa.
 * @param {number} datos.S - Area transversal de la columna (m2).
 * @returns {Object} Objecto con la altura equivalente a una unidad de transferencia basada en la fase líquida (HOL), el número de unidades de transferencia basado en la fase líquida (NOL) y la altura de la torre de absorción (z) en m.
 */
function zabKx({L2, x2, L1, x1, y1, y2, Kx, fy, S}={}) {
    const xa1 = fy(y1);
    const xa2 = fy(y2);
    const xam1 = ((1-x1)-(1-xa1))/Math.log((1-x1)/(1-xa1));
    const xam2 = ((1-x2)-(1-xa2))/Math.log((1-x2)/(1-xa2));
    const xma = ((xa1-x1)-(xa2-x2))/Math.log((xa1-x1)/(xa2-x2));
    const Lm= (L1/xam1+L2/xam2)/2;
    const HOL = Lm/(Kx*S);
    const xm = (xam1/(1-x1)+xam2/(1-x2))/2;
    const NOL = xm*(x1-x2)/xma;
    const z = HOL*NOL;
    return {HOL, NOL, z};
}




/**
 * Calcula la altura de una torre de absorción rellena (mezclas gaseosas diluidas), basado en la fase gaseosa, con el coeficiente global ideal.<br>La función recibe los flujos y fracciones molares iniciales y finales del gas (V1, y1, V2, y2), las fracciones molares inicial y final en el líquido (x2, x1), el coeficiente global ideal de la fase gaseosa (kpy), la función de equilibrio dependientes de la fracción molar en las fase gaseosa (fy) y el área transversal de la columna (S).
 * @param {Object} datos - V1, y1, V2, y2, x2, x1, Kpy, fx, S.
 * @param {number} datos.V1 - Flujo molar inicial de la fase gaseosa.
 * @param {number} datos.y1 - Fracción molar inicial de soluto en la fase gaseosa.
 * @param {number} datos.V2 - Flujo molar final de la fase gaseosa.
 * @param {number} datos.y2 - Fracción molar final de soluto en la fase gaseosa.
 * @param {number} datos.x2 - Fracción molar inicial de soluto en la fase líquida.
 * @param {number} datos.x1 - Fracción molar final de soluto en la fase líquida.
 * @param {number} datos.Kpy - Coeficiente individual ideal en la fase gaseosa en kmol/s*m2.
 * @param {number} datos.fx - Función de equilibrio dependiente de la fracción molar en la fase líquida.
 * @param {number} datos.fy - Función de equilibrio dependiente de la fracción molar en la fase gaseosa.
 * @param {number} datos.S - Area transversal de la columna (m2).
 * @returns {Object} Objecto con la altura equivalente a una unidad de transferencia basada en la fase gaseosa (HG), el número de unidades de transferencia basado en la fase gaseosa (NG) y la altura de la torre de absorción (z) en m.
 */
function zabKyi({V1, y1, V2, y2, x2, x1, Kpy, fx, S}={}) {
    const V = (V1+V2)/2;
    const HOG = V/(Kpy*S);
  	const ya1 = fx(x1);
  	const ya2 = fx(x2);
    const yam1 = ((1-ya1)-(1-y1))/Math.log((1-ya1)/(1-y1));
    const yam2 = ((1-ya2)-(1-y2))/Math.log((1-ya2)/(1-y2));
    const yma = ((y1-ya1)-(y2-ya2))/Math.log((y1-ya1)/(y2-ya2));
    const ym = (yam1/(1-y1)+yam2/(1-y2))/2;
    const NOG = ym*(y1-y2)/yma;
    const z = HOG*NOG;
    return {HOG, NOG, z};
}




/**
 * Calcula, la altura de una torre de absorción rellena (mezclas gaseosas diluidas), basado en la fase gaseosa, con el coeficiente global real.<br>La función recibe los flujos y fracciones molares iniciales y finales del gas (V1, y1, V2, y2), las fracciones molares inicial y final en el líquido (x2, x1), el coeficiente global real de la fase gaseosa (Ky), las función de equilibrio dependientes de la fracción molar en la fase líquida (fx) y el área transversal de la columna (S).
 * @param {Object} datos - V1, y1, V2, y2, x2, x1, Ky, fx, S.
 * @param {number} datos.V1 - Flujo molar inicial de la fase gaseosa.
 * @param {number} datos.y1 - Fracción molar inicial de soluto en la fase gaseosa.
 * @param {number} datos.V2 - Flujo molar final de la fase gaseosa.
 * @param {number} datos.y2 - Fracción molar final de soluto en la fase gaseosa.
 * @param {number} datos.x2 - Fracción molar inicial de soluto en la fase líquida.
 * @param {number} datos.x1 - Fracción molar final de soluto en la fase líquida.
 * @param {number} datos.Ky - Coeficiente individual ideal en la fase gaseosa en kmol/s*m2.
 * @param {number} datos.fx - Función de equilibrio dependiente de la fracción molar en la fase líquida.
 * @param {number} datos.S - Area transversal de la columna (m2).
 * @returns {Object} Objecto con la altura equivalente a una unidad de transferencia basada en la fase gaseosa (HOG), el número de unidades de transferencia basado en la fase gaseosa (NOG) y la altura de la torre de absorción (z) en m.
 */
function zabKy({V1, y1, V2, y2, x2, x1, Ky, fx, S}={}) {
    const ya1 = fx(x1);
    const ya2 = fx(x2);
    const yam1 = ((1-ya1)-(1-y1))/Math.log((1-ya1)/(1-y1));
    const yam2 = ((1-ya2)-(1-y2))/Math.log((1-ya2)/(1-y2));
    const yma = ((y1-ya1)-(y2-ya2))/Math.log((y1-ya1)/(y2-ya2));
    const Vm = (V1/yam1+V2/yam2)/2;
    const HOG = Vm/(Ky*S);
    const ym = (yam1/(1-y1)+yam2/(1-y2))/2;
    const NOG = ym*(y1-y2)/yma;
    const z = HOG*NOG;
    return {HOG, NOG, z};
}



/**
 * Calcula la altura de una torre de absorción rellena (mezclas gaseosas concentradas), basado en la fase líquida, con los coeficientes individuales ideales.<br>La función recibe la fracción molar inicial en la fase líquida (x2), el flujo y fracción molar finales del líquido (L1, x1), el flujo y fracción molar iniciales en el gas (V1, y1), la función que devuelve los coeficientes individuales ideales de las fases líquida y gaseosa (fkpxkpy), las funciones de equilibrio dependientes de las fracciones molares en las fases gaseosa y líquida (fx, fy) y el área transversal de la columna (S).
 * @param {Object} datos - x2, L1, x1, V1, y1, fkpxkpy, fx, fy, S.
 * @param {number} datos.x2 - Fracción molar inicial de soluto en la fase líquida.
 * @param {number} datos.L1 - Flujo molar final de la fase líquida.
 * @param {number} datos.x1 - Fracción molar final de soluto en la fase líquida.
 * @param {number} datos.V1 - Flujo molar inicial de la fase gaseosa.
 * @param {number} datos.y1 - Fracción molar inicial de soluto en la fase gaseosa.
 * @param {number} datos.fkpxkpy - Función que devuelve los coeficientes individuales ideales en las fases líquida y gaseosa, en kmol/s*m2.
 * @param {number} datos.fx - Función de equilibrio dependiente de la fracción molar en la fase líquida.
 * @param {number} datos.fy - Función de equilibrio dependiente de la fracción molar en la fase gaseosa.
 * @param {number} datos.S - Area transversal de la columna (m2).
 * @returns {number} Altura de la torre de absorción en m.
 */
function abkxi({x2, L1, x1, V1, y1, fkpxkpy, fx, fy, S}={}) {
	const B = L1*(1-x1);
    const C = V1*(1-y1);
    const D1 = B*x1/(1-x1)-C*y1/(1-y1);
    const fkxi = x => {
        const D2 = (B*x/(1-x)-D1)/C;
        const y = D2/(1+D2);
        const L = B/(1-x);
        const V = C/(1-y);
        const {kpx, kpy} = fkpxkpy({L, x, V, y, S});
        const {xi, xim} = xiyi2({x, y, kpx, kpy, fx, fy});
        return xim*L/(kpx*S*(1-x)*(xi-x));
    };
    return fkxi.simpson({a:x2, b:x1, n:10});
}



/**
 * Calcula la altura de una torre de absorción rellena (mezclas gaseosas concentradas), basado en la fase líquida, con los coeficientes individuales reales.<br>La función recibe la fracción molar inicial en la fase líquida (x2), el flujo y fracción molar finales del líquido (L1, x1), el flujo y fracción molar iniciales en el gas (V1, y1), la función que devuelve los coeficientes individuales reales de las fases líquida y gaseosa (fkxky), las funciones de equilibrio dependientes de las fracciones molares en las fases gaseosa y líquida (fx, fy) y el área transversal de la columna (S).
 * @param {Object} datos - x2, L1, x1, V1, y1, fkxky, fx, fy, S.
 * @param {number} datos.x2 - Fracción molar inicial de soluto en la fase líquida.
 * @param {number} datos.L1 - Flujo molar final de la fase líquida.
 * @param {number} datos.x1 - Fracción molar final de soluto en la fase líquida.
 * @param {number} datos.V1 - Flujo molar inicial de la fase gaseosa.
 * @param {number} datos.y1 - Fracción molar inicial de soluto en la fase gaseosa.
 * @param {number} datos.fkxky - Función que devuelve los coeficientes individuales reales en las fases líquida y gaseosa, en kmol/s*m2.
 * @param {number} datos.fx - Función de equilibrio dependiente de la fracción molar en la fase líquida.
 * @param {number} datos.fy - Función de equilibrio dependiente de la fracción molar en la fase gaseosa.
 * @param {number} datos.S - Area transversal de la columna (m2).
 * @returns {number} Altura de la torre de absorción en m.
 */
function abkx({x2, L1, x1, V1, y1, fkxky, fx, fy, S}={}) {
    const B = L1*(1-x1);
    const C = V1*(1-y1);
    const D1 = B*x1/(1-x1)-C*y1/(1-y1);
    const fkxi = x => {
        const D2 = (B*x/(1-x)-D1)/C;
        const y = D2/(1+D2);
        const L = B/(1-x);
        const V = C/(1-y);
        const {kx, ky} = fkxky({L, x, V, y, S});
        const {xi} = xiyi1({x, y, kpx:kx, kpy:ky, fx, fy});
        return L/(kx*S*(1-x)*(xi-x));
    };
    return fkxi.simpson({a:x2, b:x1, n:10});
}



/**
 * Calcula la altura de una torre de absorción rellena (mezclas gaseosas concentradas), basado en la fase gaseosa, con los coeficientes individuales ideales.<br>La función recibe el flujo y fracción molar finales del líquido (L1, x1), el flujo y fracción molar iniciales en el gas (V1, y1), la fracción molar final del gas (y2), la función que devuelve los coeficientes individuales ideales de las fases líquida y gaseosa (fkpxkpy), las funciones de equilibrio dependientes de las fracciones molares en las fases gaseosa y líquida (fx, fy) y el área transversal de la columna (S).
 * @param {Object} datos - L1, x1, V1, y1, y2, fkpxkpy, fx, fy, S.
 * @param {number} datos.L1 - Flujo molar final de la fase líquida.
 * @param {number} datos.x1 - Fracción molar final de soluto en la fase líquida.
 * @param {number} datos.V1 - Flujo molar inicial de la fase gaseosa.
 * @param {number} datos.y1 - Fracción molar inicial de soluto en la fase gaseosa.
 * @param {number} datos.y2 - Fracción molar final de soluto en la fase gaseosa.
 * @param {number} datos.fkpxkpy - Función que devuelve los coeficientes individuales ideales en las fases líquida y gaseosa, en kmol/s*m2.
 * @param {number} datos.fx - Función de equilibrio dependiente de la fracción molar en la fase líquida.
 * @param {number} datos.fy - Función de equilibrio dependiente de la fracción molar en la fase gaseosa.
 * @param {number} datos.S - Area transversal de la columna (m2).
 * @returns {number} Altura de la torre de absorción en m.
 */
function abkyi({L1, x1, V1, y1, y2, fkpxkpy, fx, fy, S}={}) {
    const B = L1*(1-x1);
    const C = V1*(1-y1);
    const D1 = B*x1/(1-x1)-C*y1/(1-y1);
    const fkxi = y => {
        const D3 = (C*y/(1-y)+D1)/B;
        const x = D3/(1+D3);
        const L = B/(1-x);
        const V = C/(1-y);
        const {kpx, kpy} = fkpxkpy({L, x, V, y, S});
        const {yi, yim} = xiyi2({x, y, kpx, kpy, fx, fy});
        return yim*V/(kpy*S*(1-y)*(y-yi));
    };
    return fkxi.simpson({a:y2, b:y1, n:10});
}



/**
 * Calcula la altura de una torre de absorción rellena (mezclas gaseosas concentradas), basado en la fase gaseosa, con los coeficientes individuales reales.<br>La función recibe el flujo y fracción molar finales del líquido (L1, x1), el flujo y fracción molar iniciales en el gas (V1, y1), la fracción molar final del gas (y2), la función que devuelve los coeficientes individuales reales de las fases líquida y gaseosa (fkxky), las funciones de equilibrio dependientes de las fracciones molares en las fases gaseosa y líquida (fx, fy) y el área transversal de la columna (S).
 * @param {Object} datos - L1, x1, V1, y1, y2, fkxky, fx, fy, S.
 * @param {number} datos.L1 - Flujo molar final de la fase líquida.
 * @param {number} datos.x1 - Fracción molar final de soluto en la fase líquida.
 * @param {number} datos.V1 - Flujo molar inicial de la fase gaseosa.
 * @param {number} datos.y1 - Fracción molar inicial de soluto en la fase gaseosa.
 * @param {number} datos.y2 - Fracción molar final de soluto en la fase gaseosa.
 * @param {number} datos.fpxky - Función que devuelve los coeficientes individuales reales en las fases líquida y gaseosa, en kmol/s*m2.
 * @param {number} datos.fx - Función de equilibrio dependiente de la fracción molar en la fase líquida.
 * @param {number} datos.fy - Función de equilibrio dependiente de la fracción molar en la fase gaseosa.
 * @param {number} datos.S - Area transversal de la columna (m2).
 * @returns {number} Altura de la torre de absorción en m.
 */
function abky({L1, x1, V1, y1, y2, fkxky, fx, fy, S}={}) {
    const B = L1*(1-x1);
    const C = V1*(1-y1);
    const D1 = B*x1/(1-x1)-C*y1/(1-y1);
    const fkxi = y => {
        const D3 = (C*y/(1-y)+D1)/B;
        const x = D3/(1+D3);
        const L = B/(1-x);
        const V = C/(1-y);
        const {kx, ky} = fkxky({L, x, V, y, S});
        const {yi} = xiyi1({x, y, kpx:kx, kpy:ky, fx, fy});
        return V/(ky*S*(1-y)*(y-yi));
    };
    return fkxi.simpson({a:y2, b:y1, n:10});
}



/**
 * Calcula la altura de una torre de absorción rellena (mezclas gaseosas concentradas), basado en la fase líquida, con el coeficiente global ideal de la fase líquida.<br>La función recibe la fracción molar inicial en la fase líquida (x2), el flujo y fracción molar finales del líquido (L1, x1), el flujo y fracción molar iniciales en el gas (V1, y1), la función que devuelve el coeficiente global ideal de la fase líquida (fKpx), la función de equilibrio dependiente de la fracción molar en la fase gaseosa (fy) y el área transversal de la columna (S).
 * @param {Object} datos - x2, L1, x1, V1, y1, fKpx, fy, S.
 * @param {number} datos.x2 - Fracción molar inicial de soluto en la fase líquida.
 * @param {number} datos.L1 - Flujo molar final de la fase líquida.
 * @param {number} datos.x1 - Fracción molar final de soluto en la fase líquida.
 * @param {number} datos.V1 - Flujo molar inicial de la fase gaseosa.
 * @param {number} datos.y1 - Fracción molar inicial de soluto en la fase gaseosa.
 * @param {number} datos.fKpx - Función que devuelve el coeficiente global ideal en la fase líquida, en kmol/s*m2.
 * @param {number} datos.fy - Función de equilibrio dependiente de la fracción molar en la fase gaseosa.
 * @param {number} datos.S - Area transversal de la columna (m2).
 * @returns {number} Altura de la torre de absorción en m.
 */
function abKxi({x2, L1, x1, V1, y1, fKpx, fy, S}={}) {
    const B = L1*(1-x1);
    const C = V1*(1-y1);
    const D1 = B*x1/(1-x1)-C*y1/(1-y1);
    const fKxi = x => {
        const D2 = (B*x/(1-x)-D1)/C;
        const y = D2/(1+D2);
        const L = B/(1-x);
        const V = C/(1-y);
        const Kpx = fKpx();
      	const xa = fy(y);
        const xam = ((1-x)-(1-xa))/Math.log((1-x)/(1-xa));
        return xam*L/(Kpx*S*(1-x)*(xa-x));
    };
    const z = fKxi.simpson({a:x2, b:x1, n:10});
    return z;
}



/**
 * Calcula la altura de una torre de absorción rellena (mezclas gaseosas concentradas), basado en la fase líquida, con el coeficiente global real.<br>La función recibe la fracción molar inicial en la fase líquida (x2), el flujo y fracción molar finales del líquido (L1, x1), el flujo y fracción molar iniciales en el gas (V1, y1), la función que devuelve el coeficiente global real de la fase líquida (fKx), las función de equilibrio dependiente de la fracción molar en la fase gaseosa (fy) y el área transversal de la columna (S).
 * @param {Object} datos - x2, L1, x1, V1, y1, fKx, fy, S.
 * @param {number} datos.x2 - Fracción molar inicial de soluto en la fase líquida.
 * @param {number} datos.L1 - Flujo molar final de la fase líquida.
 * @param {number} datos.x1 - Fracción molar final de soluto en la fase líquida.
 * @param {number} datos.V1 - Flujo molar inicial de la fase gaseosa.
 * @param {number} datos.y1 - Fracción molar inicial de soluto en la fase gaseosa.
 * @param {number} datos.fKx - Función que devuelve el coeficiente global real en la fase líquida, en kmol/s*m2.
 * @param {number} datos.fy - Función de equilibrio dependiente de la fracción molar en la fase gaseosa.
 * @param {number} datos.S - Area transversal de la columna (m2).
 * @returns {number} Altura de la torre de absorción en m.
 */
function abKx({x2, L1, x1, V1, y1, fKx, fy, S}={}) {
    const B = L1*(1-x1);
    const C = V1*(1-y1);
    const D1 = B*x1/(1-x1)-C*y1/(1-y1);
    const fK = x => {
        const D2 = (B*x/(1-x)-D1)/C;
        const y = D2/(1+D2);
        const L = B/(1-x);
        const V = C/(1-y);
        const Kx = fKx();
      	const xa = fy(y);
        return L/(Kx*S*(1-x)*(xa-x));
    };
    const z = fK.simpson({a:x2, b:x1, n:10});
    return z;
}



/**
 * Calcula la altura de una torre de absorción rellena (mezclas gaseosas concentradas), basado en la fase gaseosa, con el coeficiente global ideal de la fase gaseosa.<br>La función recibe el flujo y fracción molar finales del líquido (L1, x1), el flujo y fracción molar iniciales en el gas (V1, y1), la fracción molar final en la fase gaseosa (y2), la función que devuelve el coeficiente global ideal de la fase gaseosa (fKpy), la función de equilibrio dependiente de la fracción molar en la fase líquida (fx) y el área transversal de la columna (S).
 * @param {Object} datos - L1, x1, V1, y1, y2, fKpy, fy, S.
 * @param {number} datos.L1 - Flujo molar final de la fase líquida.
 * @param {number} datos.x1 - Fracción molar final de soluto en la fase líquida.
 * @param {number} datos.V1 - Flujo molar inicial de la fase gaseosa.
 * @param {number} datos.y1 - Fracción molar inicial de soluto en la fase gaseosa.
 * @param {number} datos.y2 - Fracción molar final de soluto en la fase gaseosa.
 * @param {number} datos.fKpy - Función que devuelve el coeficiente global ideal en la fase gaseosa, en kmol/s*m2.
 * @param {number} datos.fx - Función de equilibrio dependiente de la fracción molar en la fase líquida.
 * @param {number} datos.S - Area transversal de la columna (m2).
 * @returns {number} Altura de la torre de absorción en m.
 */
function abKyi({L1, x1, V1, y1, y2, fKpy, fx, S}={}) {
    const B = L1*(1-x1);
    const C = V1*(1-y1);
    const D1 = B*x1/(1-x1)-C*y1/(1-y1);
    const fKyi = y => {
        const D3 = (C*y/(1-y)+D1)/B;
        const x = D3/(1+D3);
        const L = B/(1-x);
        const V = C/(1-y);
        const Kpy = fKpy();
      	const ya = fx(x);
        const yam = ((1-ya)-(1-y))/Math.log((1-ya)/(1-y));
        return yam*V/(Kpy*S*(1-y)*(y-ya));
    };
    const z = fKyi.simpson({a:y2, b:y1, n:10});
    return z;
}



/**
 * Calcula la altura de una torre de absorción rellena (mezclas gaseosas concentradas), basado en la fase gaseosa, con el coeficiente global real de la fase gaseosa.<br>La función recibe el flujo y fracción molar finales del líquido (L1, x1), el flujo y fracción molar iniciales en el gas (V1, y1), la fracción molar final en la fase gaseosa (y2), la función que devuelve el coeficiente global real de la fase gaseosa (fKy), la función de equilibrio dependiente de la fracción molar en la fase líquida (fx) y el área transversal de la columna (S).
 * @param {Object} datos - L1, x1, V1, y1, y2, fKy, fy, S.
 * @param {number} datos.L1 - Flujo molar final de la fase líquida.
 * @param {number} datos.x1 - Fracción molar final de soluto en la fase líquida.
 * @param {number} datos.V1 - Flujo molar inicial de la fase gaseosa.
 * @param {number} datos.y1 - Fracción molar inicial de soluto en la fase gaseosa.
 * @param {number} datos.y2 - Fracción molar final de soluto en la fase gaseosa.
 * @param {number} datos.fKy - Función que devuelve el coeficiente global real en la fase gaseosa, en kmol/s*m2.
 * @param {number} datos.fx - Función de equilibrio dependiente de la fracción molar en la fase líquida.
 * @param {number} datos.S - Area transversal de la columna (m2).
 * @returns {number} Altura de la torre de absorción en m.
 */
function abKy({L1, x1, V1, y1, y2, fKy, fx, S}={}) {
    const B = L1*(1-x1);
    const C = V1*(1-y1);
    const D1 = B*x1/(1-x1)-C*y1/(1-y1);
    const fK = y => {
        const D3 = (C*y/(1-y)+D1)/B;
        const x = D3/(1+D3);
        const L = B/(1-x);
        const V = C/(1-y);
        const Ky = fKy();
      	const ya = fx(x);
        return V/(Ky*S*(1-y)*(y-ya));
    };
    const z = fK.simpson({a:y2, b:y1, n:10});
    return z;
}




/**
 * Calcula los flujos y composiciones finales, en la extracción líquido-líquido, en una etapa de equilibrio conformada por el soluto (A), el solvente de la fase refinado (B) y el solvente de la fase extracto (C).
 * @param {Object} datos - Flujos, composiciones iniciales de las fases refinado y extracto y funciones de equilibrio.
 * @param {number} datos.L0 - Flujo inicial de la fase refinado.
 * @param {number} datos.xA0 - Fracción inicial de soluto en la fase refinado.
 * @param {number} datos.xB0 - Fracción inicial del solvente B, en la fase refinado.
 * @param {number} datos.V0 - Flujo inicial de la fase extracto.
 * @param {number} datos.yA0 - Fracción inicial de soluto en la fase extracto.
 * @param {number} datos.yB0 - Fracción inicial del solvente B, en la fase extracto.
 * @param {Function} datos.FXAXB - Función de equilibrio que, para un valor conocido de xA, devuelve el correspondiente valor de XB.
 * @param {Function} datos.FXAYA - Función de equilibrio que, para un valor conocido de xA, devuelve el correspondiente valor de yA.
 * @param {Function} datos.FXAYB - Función de equilibrio que, para un valor conocido de xA, devuelve el correspondiente valor de yB.
 * @returns {Object} Objeto con los flujos y composiciones finales de la fase refinado (L1, xA1, XB1) y de la fase extracto (V1, yA1, yB1).
 */
function ELL1({L0, xA0, xB0, V0, yA0, yB0, FXAXB, FXAYA, FXAYB}={}) {
    const M = L0+V0;
    const CA0 = L0*xA0+V0*yA0;
    const CB0 = L0*xB0+V0*yB0;
    let xB1, yA1, yB1, L1, V1;
    const fXA1 = xA1 => {
        xB1 = FXAXB(xA1);
        yA1 = FXAYA(xA1);
        yB1 = FXAYB(xA1);
        L1 = (CA0*yB1-CB0*yA1)/(xA1*yB1-xB1*yA1);
        V1 = (CA0-L1*xA1)/yA1;
        return L1+V1-M;
    };
    const xA1 = fXA1.newton({xi:CA0/M/2, err:9});
    return {L1, xA1, xB1, V1, yA1, yB1};
}



/**
 * Calcula los flujos y composiciones finales faltantes, en la extracción líquido-líquido, en múltiples etapas de equilibrio. El sistema está conformado por el soluto (A), el solvente de la fase refinado (B) y el solvente de la fase extracto (C).
 * @param {Object} datos - Flujos, composiciones iniciales de las fases refinado y extracto y funciones de equilibrio.
 * @param {number} datos.L0 - Flujo inicial de la fase refinado.
 * @param {number} datos.xA0 - Fracción inicial de soluto en la fase refinado.
 * @param {number} datos.xB0 - Fracción inicial del solvente B, en la fase refinado.
 * @param {number} datos.xAN - Fracción final de soluto en la fase refinado.
 * @param {number} datos.V0 - Flujo inicial de la fase extracto.
 * @param {number} datos.yA0 - Fracción inicial de soluto en la fase extracto.
 * @param {number} datos.yB0 - Fracción inicial del solvente B, en la fase extracto.
 * @param {Function} datos.FXAXB - Función de equilibrio que, para un valor conocido de xA, devuelve el correspondiente valor de XB.
 * @param {Function} datos.FYAYB - Función de equilibrio que, para un valor conocido de yA, devuelve el correspondiente valor de yB.
 * @returns {Object} Objeto con los flujos y composiciones finales faltantes de la fase refinado (L1, XBN) y de la fase extracto (V1, yA1, yB1).
 */
function ELLN1({L0, xA0, xB0, xAN, V0, yA0, yB0, FXAXB, FYAYB}={}) {
    const M = L0+V0;
    const CA0 = L0*xA0+V0*yA0;
    const CB0 = L0*xB0+V0*yB0;
    const xBN = FXAXB(xAN);
    let yB1, LN, V1;
    const fYA1 = yA1 => {
        yB1 = FYAYB(yA1);
        LN = (CA0*yB1-CB0*yA1)/(xAN*yB1-xBN*yA1);
        V1 = (CA0-LN*xAN)/yA1;
        return LN+V1-M;
    };
    const yA1 = fYA1.newton({xi:CA0/M});
    return {LN, xBN, V1, yA1, yB1};
}



/**
 * Calcula los flujos y composiciones finales faltantes, en la extracción líquido-líquido, en múltiples etapas de equilibrio. El sistema está conformado por el soluto (A), el solvente de la fase refinado (B) y el solvente de la fase extracto (C).
 * @param {Object} datos - Flujos, composiciones iniciales de las fases refinado y extracto y funciones de equilibrio.
 * @param {number} datos.L0 - Flujo inicial de la fase refinado.
 * @param {number} datos.xA0 - Fracción inicial de soluto en la fase refinado.
 * @param {number} datos.xB0 - Fracción inicial del solvente B, en la fase refinado.
 * @param {number} datos.LNxAN - Flujo final de soluto en la fase refinado.
 * @param {number} datos.V0 - Flujo inicial de la fase extracto.
 * @param {number} datos.yA0 - Fracción inicial de soluto en la fase extracto.
 * @param {number} datos.yB0 - Fracción inicial del solvente B, en la fase extracto.
 * @param {Function} datos.FXAXB - Función de equilibrio que, para un valor conocido de xA, devuelve el correspondiente valor de XB.
 * @param {Function} datos.FYAYB - Función de equilibrio que, para un valor conocido de yA, devuelve el correspondiente valor de yB.
 * @returns {Object} Objeto con los flujos y composiciones faltantes de la fase refinado (LN, xAN, xBN,) y de la fase extracto (V1, yA1, yB1).
 */
function ELLN2({L0, xA0, xB0, LNxAN, V0, yA0, yB0, FXAXB, FYAYB}={}) {
    const M = L0+V0;
    const CA0 = L0*xA0+V0*yA0;
    const CB0 = L0*xB0+V0*yB0;
    let yB1, LN, V1, xAN, xBN;
    const fYA1 = yA1 => {
        yB1 = FYAYB(yA1);
        V1 = (CA0-LNxAN)/yA1;
        LN = M-V1;
      	xAN = LNxAN/LN;
      	xBN = FXAXB(xAN);
        return LN*xBN+V1*yB1-CB0;
    };
    const yA1 = fYA1.newton({xi:CA0/M});
    return {LN, xAN, xBN, V1, yA1, yB1};
}



/**
 * Calcula el número de etapas requeridas, en la extracción líquido-líquido, en múltiples etapas de equilibrio. El sistema está conformado por el soluto (A), el solvente de la fase refinado (B) y el solvente de la fase extracto (C).
 * @param {Object} datos - Flujo y composición inicial y final de las fases refinado; flujo y composición final de la fase extracto y funciones de equilibrio.
 * @param {number} datos.L0 - Flujo inicial de la fase refinado.
 * @param {number} datos.xA0 - Fracción inicial de soluto en la fase refinado.
 * @param {number} datos.xB0 - Fracción inicial del solvente B, en la fase refinado.
 * @param {number} datos.xAN - Fracción final de soluto en la fase refinado.
 * @param {number} datos.V1 - Flujo final de la fase extracto.
 * @param {number} datos.yA1 - Fracción final de soluto en la fase extracto.
 * @param {number} datos.yB1 - Fracción final del solvente B, en la fase extracto.
 * @param {Function} datos.FYAYB - Función de equilibrio que, para un valor conocido de yA, devuelve el correspondiente valor de yB.
 * @param {Function} datos.FYAXA - Función de equilibrio que, para un valor conocido de yA, devuelve el correspondiente valor de xA.
 * @param {Function} datos.FYAXA - Función de equilibrio que, para un valor conocido de yA, devuelve el correspondiente valor de xB.
 * @returns {number} Número de etapas de equilibrio.
 */
function ELLN({L0, xA0, xB0, xAN, V1, yA1, yB1, FYAYB, FYAXA, FYAXB}={}) {
    const DM = L0-V1;
    const DA = L0*xA0-V1*yA1;
    const DB = L0*xB0-V1*yB1;
    let yAI = yA1, xAI, xBI, N=0;
    var fYAJ = yAJ => {
        const yBJ = FYAYB(yAJ);
        const LI = (DA*yBJ-DB*yAJ)/(xAI*yBJ-xBI*yAJ);
        const VJ = (LI*xAI-DA)/yAJ;
        return LI-VJ-DM
    }
    for(;;) {
        xAI = FYAXA(yAI);
        xBI = FYAXB(yAI);
        if (xAI<=xAN) break;
        yAI = fYAJ.newton({xi:yAI*0.5});
        xA0 = xAI;
        N++
    }
    return N+(xAN-xA0)/(xAI-xA0);
}




/**
 * Calcula el flujo mínimo de refinado, en la extracción líquido-líquido, en múltiples etapas de equilibrio. El sistema está conformado por el soluto (A), el solvente de la fase refinado (B) y el solvente de la fase extracto (C).
 * @param {Object} datos - Flujo y composición inicial y final de las fases refinado; composición inicial de la fase extracto y funciones de equilibrio.
 * @param {number} datos.L0 - Flujo inicial de la fase refinado.
 * @param {number} datos.xA0 - Fracción inicial de soluto en la fase refinado.
 * @param {number} datos.xB0 - Fracción inicial del solvente B, en la fase refinado.
 * @param {number} datos.xAN - Fracción final de soluto en la fase refinado.
 * @param {number} datos.yA0 - Fracción inicial de soluto en la fase extracto.
 * @param {number} datos.yB0 - Fracción inicial del solvente B, en la fase extracto.
 * @param {Function} datos.FXAXB - Función de equilibrio que, para un valor conocido de xA, devuelve el correspondiente valor de xB.
 * @param {Function} datos.FXAYA - Función de equilibrio que, para un valor conocido de xA, devuelve el correspondiente valor de yA.
 * @param {Function} datos.FXAYB - Función de equilibrio que, para un valor conocido de xA, devuelve el correspondiente valor de yB.
 * @returns {number} Flujo mínimo de refinado.
 */
function ELMIN({L0, xA0, xB0, xAN, yA0, yB0, FXAXB, FXAYA, FXAYB}={}) {
    const xBN = FXAXB(xAN);
    const yA1 = FXAYA(xA0);
    const yB1 = FXAYB(xA0);
    const Vmin = ((yA1-xAN)*L0*(xB0-xBN)-(yB1-xBN)*L0*(xA0-xAN))/
          ((yA1-xAN)*(xBN-yB0)-(yB1-xBN)*(xAN-yA0));
    return {Vmin};
}



/**
 * Calcula, en base a las fracciones en peso en base a la solución sin sólido inerte, las fracciones en peso estándar, es decir, en base a la solución total (incluído el sólido inerte).
 * @param {Object} datos - Fracciones en peso con relación a la solución sin sólido inerte.
 * @param {Array} datos.Ny - Fracciones en peso de sólido inerte con relación a la solución sin sólido inerte (fase sólida).
 * @param {Array} datos.Nx - Fracciones en peso de solido inerte con relación a la solución sin sólido inerte (fase líquida).
 * @param {Array} datos.yAA - Fracciones en peso de soluto con relación a la solución sin sólido inerte (fase sólida).
 * @param {Array} datos.xAA - Fracciones en peso de soluto con relación a la solución sin sólod inerte (fase líquida).
 * @returns {Object} Objeto con las fracciones en peso de sólido inerte en las fases sólida y líquida (vxb, vyb) y las fracciones en peso de soluto en las fases sólida y líquida (vxa, vya).
 */
function LIX1({Ny, yAA, Nx, xAA}={}) {
    const xB = Ny.map(n=>n/(1+n));
    const xA = yAA.map((y,i)=>y/(1+Ny[i]));
    const yB = Nx.map(n=>n/(1+n));
    const yA = xAA.map((x,i)=>x/(1+Nx[i]));
    return {xB, xA, yB, yA}
}



/**
 * Calcula, para un equipo de adsorción discontinuo (por lotes), las concentraciones finales de soluto en la solución y en el adsorbente. 
 * @param {Object} datos - 
 * @param {Array} datos.cf - Concentración inicial de soluto en la solución.
 * @param {Array} datos.qf - Concentración inicial de soluto en el adsorbente.
 * @param {Array} datos.S - Volumen de solución a tratar.
 * @param {Array} datos.M - Masa de adsorbente empleada.
 * @param {Array} datos.fc - Función de equilibrio. Devuelve la concentración de soluto en el adsorbente para un valor conocido de la concentración de soluto en la solución.
 * @returns {Object} Objeto con la concentración final de soluto en la solución (c) y la concentración final de soluto en el adsrobente (q).
 */
function Adb1({cf, qf, S, M, fc}={}) {
    const C0 = qf*M+cf*S;
    let q;
    function f(c) {
      q = fc(c);
      return q*M+c*S-C0;
    }
    const c = f.newton({xi:cf/2, err:9});
    return {c, q};
}




/**
 * Calcula la temperatura del punto de burbuja en K, para una mezcla binaria. La función recibe la presión de operación (P en mmHg), la fracción molar del compuesto 1 (x1), las constantes de Antoine para los compuestos 1 y 2 (antA1, antB1, antC1, antA2, antB2, antC2), las temperaturas de ebullición de los compuestos 1 y 2 (Tb1 y Tb2 en K), los parámetros de la ecuación de antoine (alfa=alfa12=alfa21, deltag12 y deltag21). Los valores por defecto corresponden a la mezcla binaria etanol - agua (compuestos 1 y 2).<br>
 * La función devuelve la temperatura del punto de burbuja en K.
 * @param {Object} datos - Presión, fracción molar y parámetros de las ecuaciones de Antoine y NRTL.
 * @param {number} datos.P - Presión de operación en mmHg.
 * @param {number} datos.x1 - Fracción molar del compuesto 1.
 * @param {number} [datos.antA1=18.9119] - Primera constante de la ecuación de Antoine para el compuesto 1.
 * @param {number} [datos.antB1=3803.98] - Segunda constante de la ecuación de Antoine para el compuesto 1.
 * @param {number} [datos.antC1=-41.68] - Tercera constante de la ecuación de Antoine para el compuesto 1.
 * @param {number} [datos.Tb1=351.4] - Temperatura de ebullición normal del compuesto 1 en K.
 * @param {number} [datos.antA2=18.3036] - Primera constante de la ecuación de Antoine para el compuesto 2.
 * @param {number} [datos.antB2=3816.44] - Segunda constante de la ecuación de Antoine para el compuesto 2.
 * @param {number} [datos.antC2=-46.13] - Tercera constante de la ecuación de Antoine para el compuesto 2.
 * @param {number} [datos.Tb2=373.15] - Temperatura de ebullición normal del compuesto 1 en K.
 * @param {number} [datos.alfa=0.2974] - Parámetro no aleatorio entre los compuestos 1 y 2 y 2 y 1 de la ecuación NRTL.
 * @param {number} [datos.deltag12=-121.27] - Parámetro de energía de interacción entre los compuestos 1 y 2 de la ecuación NRTL.
 * @param {number} [datos.deltag21=1337.86] - Parámetro de energía de interacción entre los compuestos 2 y 1 de la ecuación NRTL.
 * @returns {number} Temperatura del punto de burbuja en K.
 * @example <caption>Ejemplos de uso de la función.</caption>
 * // Cálculo de la temperatura de burbuja para la una mezcla de etanol-agua a 543 mmHg, con una fracción molar de etanol igual a 0.78.
 * tBurbujaNRTL({P:542, x1:0.78})
 * 343.31107809691065  // resultado obtenido
 * 
 * // Temperatura de burbuja para la mezcla etanol-agua a 760 mmHg, con una fracción molar de etanol igual a 0.46. 
 * tBurbujaNRTL({P:760, x1:0.46})
 * 353.35790057312795  // resultado obtenido
 */
function tBurbujaNRTL({P, x1,
                      antA1=18.9119, antB1=3803.98, antC1=-41.68, Tb1=351.4,
                      antA2=18.3036, antB2=3816.44, antC2=-46.13, Tb2=373.15,
                      alfa=0.2974, deltag12=-121.27, deltag21=1337.86}={}) {
    const R = 1.987207;
  	const x2 = 1-x1;
    let y1, y2, k1, k2;

    function ft(T) {
        const tau12 = deltag12/(R*T);
            const tau21 = deltag21/(R*T);
        const G12 = Math.exp(-alfa*tau12);
        const G21 = Math.exp(-alfa*tau21);
        const gamma1 = Math.exp(x2**2*(tau21*(G21/(x1+x2*G21))**2+tau12*G12/(x2+x1*G12)**2));
        const gamma2 = Math.exp(x1**2*(tau12*(G12/(x2+x1*G12))**2+tau21*G21/(x1+x2*G21)**2));
        const Pv1 = antoine({T, antA:antA1, antB:antB1, antC:antC1});
        const Pv2 = antoine({T});
        K1 = gamma1*Pv1/P;
        K2 = gamma2*Pv2/P;
        y1 = K1*x1;
        y2 = K2*x2;
        return y1+y2-1;
    }

    const Tburbuja = ft.newton({xi:(Tb1+Tb2)/2});
    return {Tburbuja, y1, y2, K1, K2};

}



/**
 * Calcula la constante de equilibrio líquido-vapor, a una temperatura (T) y presión (P) dadas, para el compuesto especificado (compuesto), empleando la ecuación de McWilliams.
 * @function kMcWilliams
 * @param {Object} datos - T, P y nombre del compuesto.
 * @param {number} datos.T - Temperatura de operación en R.
 * @param {number} datos.P - Presión de operación en psi.
 * @param {string} datos.compuesto - Nombre del compuesto para el cual se calcula la constante de equilibrio.
 * @returns {number} Constante de equilibrio calculada.
 * @example <caption>Ejemplos de uso de la función</caption>
 * // Cálculo de las constantes de equilibrio para el isopentano, n-pentano y n-hexano a 581.58 R y 14.7 psi:
 * kMcWilliams({T:581.58, P:14.7, compuesto:"isopentano"})
 * 2.0065279392461237 // Resultado devuelto
 * 
 * kMcWilliams({T:581.58, P:14.7, compuesto:"n-pentano"})
 * 1.5324821601252245 // Resultado devuelto
 * 
 * kMcWilliams({T:581.58, P:14.7, compuesto:"n-hexano"})
 * 0.5675864021894845 // Resultado devuelto
 * 
 * _
 */
const kMcWilliams = (function() {

    // constantes conocidas de la ecuación de McWilliams
    const cMcW = [
      ["metano",-292.860,0,8.2445,-.8951,59.8465,0],
      ["etileno",-600076.875,0,7.90595,-.84677,42.94594,0],
      ["etano",-687248.25,0,7.90694,-.8860,49.02654,0],
      ["propileno",-923484.6875,0,7.71725,-.87871,47.67624,0],
      ["propano",-970688.5625,0,7.15059,-.76984,0,6.90224],
      ["isobutano",-1166846,0,7.72668,-.92213,0,0],
      ["n-butano",-1280557,0,7.94986,-.96455,0,0],
      ["isopentano",-1481583,0,7.58071,-.93159,0,0],
      ["n-pentano",-1524891,0,7.33129,-.89143,0,0],
      ["n-hexano",-1778901,0,6.96783,-.84634,0,0],
      ["n-heptano",-2013803,0,6.52914,-.79543,0,0],
      ["n-octano",0,-7646.81641,12.48457,-.73152,0,0],
      ["n-nonano",-2551040,0,5.69313,-.67818,0,0],
      ["n-decano",0,-9760.45703,13.80354,-.71470,0,0]
    ];
  
    // Ecuación de McWilliams
    return function({T, P, compuesto}={}) {
        let i = cMcW.findIndex(e=>e[0]==compuesto.toLowerCase());
        if (i<0) throw new RangeError(`No existe la constante de McWilliams para el ${compuesto}.`);
        const [aT1, aT2, aT6, ap1, ap2, ap3] = cMcW[i].slice(1);
        return Math.exp(aT1/T**2+aT2/T+aT6+ap1*Math.log(P)+ap2/P**2+ap3/P);  
    };
    
})();




/**
 * Calcula la temperatura de burbuja, a una presión (P) dada, de los compuestos con una composición en la fase líquida conocida.
 * @param {Object} datos - Temperatura asumida (Ti), presión conocida (P), compuestos y fracciones molares en la fase líquida.
 * @param {number} datos.ti - Temperatura inicial asumida (en R) para la temperatura de burbuja.
 * @param {number} datos.P - Presión de operación conocida (en psia).
 * @param {Array} datos.compuestos - Array con los nombres de los compuestos que conforman la mezcla.
 * @param {Array} datos.x - Array con las fracciones molares en la fase líquida.
 * @returns {number} Temperatura de burbuja calculada (en R).
 * 
 */
function tBurbujaMW({Ti, P, compuestos, x}={}) {

    // Sumatoria de las fracciones molares en la fase gaseosa, menos 1. a la temperatura T.
    function sumyi(T) {
        let s = 0;
        for (let i=0; i<compuestos.length; i++) {
            Ki = kMcWilliams({T, P, compuesto:compuestos[i]});
            s += Ki*x[i];
        }
        return s-1;
    }
    
    // Resolución la función sumyi, con el método de Newton
    return sumyi.newton({xi:Ti, err:9});
    
}




/**
 * Calcula la presión de burbuja, a una temperatura (T) dada, de los compuestos con una composición en la fase líquida conocida.
 * @param {Object} datos - Presión asumida (Pi), temperatura conocida (T), compuestos y fracciones molares en la fase líquida.
 * @param {number} datos.Pi - Presión inicial asumida (en psi) para la presión de burbuja.
 * @param {number} datos.T - Temperatura de operación conocida (en R).
 * @param {Array} datos.compuestos - Array con los nombres de los compuestos que conforman la mezcla.
 * @param {Array} datos.x - Array con las fracciones molares en la fase líquida.
 * @returns {number} Presión de burbuja calculada (en psi).
 * 
 */
function pBurbujaMW({Pi, T, compuestos, x}={}) {

    // Sumatoria de las fracciones molares en la fase gaseosa, menos 1, a la presión P.
    function sumyi(P) {
      let s = 0;
      for (let i=0; i<compuestos.length; i++) {
        Ki = kMcWilliams({T, P, compuesto:compuestos[i]});
        s += Ki*x[i];
      }
      return s-1;
    }

    // Resolución la función sumyi, con el método de Newton
    return sumyi.newton({xi:Pi, err:9});
  
}




/**
 * Calcula, para una presión (P) dada, la temperatura de rocío de los compuestos con una composición en la fase gaseosa conocida.
 * @param {Object} datos - Temperatura asumida (Ti), presión conocida (P), compuestos y fracciones molares en la fase gaseosa.
 * @param {number} datos.ti - Temperatura inicial asumida (en R) para la temperatura de rocío.
 * @param {number} datos.P - Presión de operación conocida (en psi).
 * @param {Array} datos.compuestos - Array con los nombres de los compuestos que conforman la mezcla.
 * @param {Array} datos.y - Array con las fracciones molares en la fase gaseosa.
 * @returns {number} Temperatura de rocío calculada (en R).
 * 
 */
function tRocíoMW({Ti, P, compuestos, y}={}) {

    // Sumatoria de las fracciones molares en la fase líquida, menos 1. a la temperatura T.
    function sumxi(T) {
      let s = 0;
      for (let i=0; i<compuestos.length; i++) {
        Ki = kMcWilliams({T, P, compuesto:compuestos[i]});
        s += y[i]/Ki;
      }
      return s-1;
    }

    // Resolución la función sumxi, con el método de Newton
    return sumxi.newton({xi:Ti, err:9});

}




/**
 * Calcula la presión de rocío, a una temperatura (T) dada, de los compuestos con una composición en la fase gaseosa conocida.
 * @param {Object} datos - Presión asumida (Pi), temperatura conocida (T), compuestos y fracciones molares en la fase gaseosa.
 * @param {number} datos.Pi - Presión inicial asumida (en psi) para la presión de rocío.
 * @param {number} datos.T - Temperatura de operación conocida (en R).
 * @param {Array} datos.compuestos - Array con los nombres de los compuestos que conforman la mezcla.
 * @param {Array} datos.x - Array con las fracciones molares en la fase gaseosa.
 * @returns {number} Presión de rocío calculada (en psi).
 * 
 */
function pRocíoMW({Pi, T, compuestos, y}={}) {

    // Sumatoria de las fracciones molares en la fase líquida, menos 1, a la presión P.
    function sumxi(P) {
      let s = 0;
      for (let i=0; i<compuestos.length; i++) {
        Ki = kMcWilliams({T, P, compuesto:compuestos[i]});
        s += y[i]/Ki;
      }
      return s-1;
    }
    
    // Resolución la función sumxi, con el método de Newton
    return sumxi.newton({xi:Pi, err:9});
    
}



/**
 * Calcula la composición líquido-vapor (x-y) en un flash adiabático.<br>La función recibe la temperatura de operación (T), la presión de operación (P), la composición (compuestos) y las fracciones molares iniciales del fluido que ingresa al equipo (z). 
 * @param {Object} datos - Temperatura (T), Presión (p), compuestos y fracciones molares (z).
 * @param {number} datos.T - Temperatura de operación del equipo en R.
 * @param {number} datos.P - Presión de operación del equipo en psi.
 * @param {Array} datos.compuestos - Array con los nombres de los compuestos que conforman la mezcla.
 * @param {Array} datos.z - Array con las fracciones molares iniciales de los compuestos en la mezcla.
 * @returns {Object} Objeto con la fracción molar de gas en la mezcla (phi), las fracciones molares de los compuestos en la fase líquida (x) y las fracciones molares de los compuestos en la fase gaseosa (y).
 */
function compLVMWR({T, P, compuestos, z}={}) {
    // Cálculo de las constantes de equilibrio
    const K = [];
    for (compuesto of compuestos) {
        K.push(kMcWilliams({T, P, compuesto}))
    }
    // Ecuación de Rachford-Rice
    function rachford(phi) {
        let s = 0;
        for (let i=0; i<z.length; i++) {
            s += z[i]*(1-K[i])/(1+phi*(K[i]-1));
        }
        return s;
    }
    // Cálculo del parámetro phi
    const phi = rachford.newton({xi:0.5, err:9}).round(9);
    // Cálculo de las composiciones líquido-vapor
    const x = [], y = [];
    for (let i=0; i<z.length; i++) {
        x.push(z[i]/(1+phi*(K[i]-1)));
        y.push(K[i]*x[i]);
    }
    // La función devuelve el parámetro phi y las composiciones en las fases líquido - vapor
    return {phi: phi.round(5), x: x.round(5), y: y.round(5)};
}





/**
 * Redondea él o los números del objeto que recibe, al número de dígitos después del punto recibido como segundo parámetro.
 * @param {Object} objeto - Objeto con los núemros a redondear.
 * @param {number} [n=0] - Número de dígitos después del punto.  
 * @returns {Object|number} objeto con sus propiedades numéricas redondeadas o número redondeado.
 */
function round(objeto, n=0) {
    switch (objeto.constructor.name) {
        case "Number": return objeto.round(n);
        case "Array": return objeto.round(n);
        case "Object": 
            const keys = Object.keys(objeto);
            const values = Object.values(objeto);
            let newValues = [], nv;
            for (let v of values) {
                switch (v.constructor.name) {
                    case "Number":
                    case "Array": nv = v.round(n); break;
                    case "objeto": nv = round(v,n); break;
                    default: nv = v;
                }
                newValues.push(nv);
            }
            let r = {};
            for (let i=0; i<keys.length; i++) {
                r[keys[i]] = newValues[i];
            }
            return r;
    }
}




/**
 * Redondea él o los números del objeto que recibe, al número de dígitos de precisión recibido como segundo parámetro.
 * @param {Object} objeto - Objeto con los núemros a redondear.
 * @param {number} [n=0] - Número de dígitos de precisión.
 * @returns {Object|number} Objeto con sus propiedades numéricas redondeadas o número redondeado.
 */
function precision(objeto, n=6) {
    switch (objeto.constructor.name) {
        case "Number": return objeto.precision(n);
        case "Array": return objeto.precision(n);
        case "Object": 
            const keys = Object.keys(objeto);
            const values = Object.values(objeto);
            let newValues = [], nv;
            for (v of values) {
                switch (v.constructor.name) {
                    case "Number":
                    case "Array": nv = nv.precision(n); break;
                    case "objeto": nv = precision(v,n); break;
                    default: nv = v;
                }
                newValues.push(nv);
            }
            let r = {};
            for (let i=0; i<keys.length; i++) {
                r[keys[i]] = newValues[i];
            }
            return r;
    }
}




/**
 * Redondea él o los números del objeto que recibe, al número de dígitos de precisión recibido como segundo parámetro, devolviendo los números como strings.
 * @param {Object} objeto - Objeto con los números a redondear.
 * @param {number} [n=0] - Número de dígitos de precisión.
 * @returns {Object|number} objeto con sus propiedades numéricas redondeadas como strings o número redondeado como string.
 */
function toPrecision(objeto, n=6) {
    switch (objeto.constructor.name) {
        case "Number": return objeto.toPrecision(n);
        case "Array": return objeto.toPrecision(n);
        case "Object": 
            const keys = Object.keys(objeto);
            const values = Object.values(objeto);
            let newValues = [], nv;
            for (v of values) {
                switch (v.constructor.name) {
                    case "Number":
                    case "Array": nv = v.toPrecision(n); break;
                    case "objeto": nv = toPrecision(nv,n); break;
                    default: nv = v;
                }
                newValues.push(nv);
            }
            let r = {};
            for (let i=0; i<keys.length; i++) {
                r[keys[i]] = newValues[i];
            }
            return r;
    }
}




/**
 * Redondea él o los números del objeto que recibe, al número de dígitos después del punto recibido como segundo parámetro, devolviendo los números redondeados como strings.
 * @param {Object|number} objeto - Objeto o número a redondear.
 * @param {number} [n=0] - Número de dígitos después del punto.
 * @returns {Object|Number} Objeto con sus propiedades numéricas redondeadas como strings o número redondeado como string.
 */
function toFixed(objeto, n=6) {
    switch (objeto.constructor.name) {
        case "Number": return objeto.toFixed(n);
        case "Array": return objeto.toFixed(n);
        case "Object": 
            const keys = Object.keys(objeto);
            const values = Object.values(objeto);
            let newValues = [], nv;
            for (v of values) {
                switch (v.constructor.name) {
                    case "Number":
                    case "Array": nv = v.toFixed(n); break;
                    case "objeto": nv = toFixed(nv,n); break;
                    default: nv = v;
                }
                newValues.push(nv);
            }
            let r = {};
            for (let i=0; i<keys.length; i++) {
                r[keys[i]] = newValues[i];
            }
            return r;
    }
}




/**
 * Redondea él o los números del objeto que recibe, al número de dígitos de precisión recibido como segundo parámetro, devolviendo los números redondeados como strings, en formato exponencial.
 * @param {Object|number} objeto - Objeto o número a redondear.
 * @param {number} [n=0] - Número de dígitos de precisión.
 * @returns {Object|Number} Objeto con sus propiedades numéricas redondeadas como strings o número redondeado como string en formato exponencial. 
 */
function toExponential(objeto, n=6) {
    switch (objeto.constructor.name) {
        case "Number": return objeto.toExponential(n);
        case "Array": return objeto.toExponential(n);
        case "Object": 
            const keys = Object.keys(objeto);
            const values = Object.values(objeto);
            let newValues = [], nv;
            for (v of values) {
                switch (v.constructor.name) {
                    case "Number":
                    case "Array": nv = v.toExponential(n); break;
                    case "objeto": nv = toExponential(nv,n); break;
                    default: nv = v;
                }
                newValues.push(nv);
            }
            let r = {};
            for (let i=0; i<keys.length; i++) {
                r[keys[i]] = newValues[i];
            }
            return r;
    }
}





/**
 * Calcula, empleando el método del trapecio, con los datos de la tabla desde la cual se llama al método [[x<sub>0</sub>, x<sub>1</sub>, ...], [y<sub>0</sub>, y<sub>1</sub>, ...]], la integral entre los límites de la tabla (entre x<sub>0</sub> y x<sub>n-1</sub>) o entre los límites inicial (<i>a</i>) y final (<i>b</i>) proporcionados.<br>
 * Los datos correspondientes a la variable independiente (x) deben estar en orden ascendente.
 * @method external:Array#trapecio
 * @param {Object} [limites] - Objecto con los límites de integración
 * @param {number} [limites.a] - Límite inferior
 * @param {number} [limites.b] - Límite superior
 * @returns {number} Integral calculada con los datos tabulados
 * @example <caption>Ejemplos de uso del método</caption>
 * // Con los siguientes datos para las variables independiente (vx) y dependiente (vy)
 * var vx = [2.0, 2.5, 2.7, 3.0, 3.3, 3.8, 4.2, 4.5, 5.1, 5.7, 6.0, 6.4, 6.9, 7.2, 7.5];
 * var vy = [4.800, 6.625, 7.509, 9.00, 10.689, 13.944, 16.944, 19.425, 24.981, 31.329, 34.800, 39.736, 46.401, 50.664, 55.125];
 * var td = [vx, vy]; // Tabla con los datos x-y
 * 
 * // Llamando al método del trapecio, para los siguientes casos, se obtienen los resultados mostrados
 * 
 * td.trapecio()
 * 134.49425 // Integral entre 2.0 y 7.5
 * 
 * td.trapecio({a:2.5, b:7.2})
 * 115.76965000000001 // Integral entre 2.5 y 7.2
 * 
 * td.trapecio({a:2.6, b:6.7})
 * 91.5117 // Integral entre 2.6 y 6.7 (límites interpolados)
 * 
 * td.trapecio({a:3.1})
 * 126.8201 // Integral entre 3.1 (límite inicial interpolado) y 7.5
 * 
 * td.trapecio({b:7.1})
 * 113.63054999999999 // Integral entre 2.0 y 7.1 (límite final interpolado)
 * 
 * td.trapecio({a:1.5, b:6.0})
 * 69.56845 // Integral entre 1.5 (límite inicial extrapolado) y 6.0
 * 
 * td.trapecio({a:3.0, b:7.8})
 * 144.95489999999998 // Integral entre 3.0 y 7.8 (límite final extrapolado)
 * 
 * td.trapecio({a:1.3, b:8.1})
 * 172.71159999999998 // Integral entre 1.3 y 8.1 (límites inferior y superior extrapolados)
 * 
 * _
 */
Array.prototype.trapecio = function({a, b}={}) {

    if (this.length<2) throw "trapecio: la tabla debe tener dos filas: datos 'x' y datos 'y'";
    let [x, y] = [this[0].slice(), this[1].slice()];
    if (x.length != y.length) throw "trapecio: las dos filas deben tener el mismo número de elementos."
    if (x.length<2) throw "trapecio: La tabla debe tener por lo menos 2 puntos";

    // Función local de interpolación lineal
    const inLin = (xc,i,j) => {
      const b = (y[i]-y[j])/(x[i]-x[j]);
      const a = y[i]-b*x[i];
      return a+b*xc;
    };
  
    //  Si se ha mandado el límite inferior de integración
    if (a!=undefined) {
        // Y el límite es menor al menor valor tabulado, se extrapola el punto y se añade a los datos
        if (a<x[0]) {
            let yc = inLin(a,0,1);	
            x.unshift(a);
            y.unshift(yc);
        // Si no corresponde a uno de los puntos tabulados se interpola el punto, se eliminan los puntos previos y se añade al principio el punto interpolado
        } else if (x.indexOf(a)<0) {
            let j = x.findIndex(e => e>a);
            if (j<0 || j>x.length-1) throw "trapecio: el límite inferior debe ser menor al último punto tabulado"
            let yc = inLin(a,j-1,j);	
            x.splice(0,j); x.unshift(a);
            y.splice(0,j); y.unshift(yc);
        // Si es uno de los puntos tabulados, se eliminan los puntos previos
        } else {
            let i = x.indexOf(a);
            x.splice(0,i);
            y.splice(0,i);
        }
    }
  
    //  Si se ha mandado el límite superior de integración
    if (b!=undefined) {
        // Y el límite es mayor al mayor valor tabulado, se extrapola el punto y se añade a los datos
        if (b>x[x.length-1]) {
            let j = x.length-1;
            let yc = inLin(b,j-1,j)
            x.push(b);
            y.push(yc);
        // Si no corresponde a uno de los puntos tabulados se interpola el punto, se eliminan los puntos posteriores y se añade al filan el punto interpolado
        } else if (x.indexOf(b)<0) {
            let j = x.findIndex(e => e>b);
            if (j<0 || j<1) throw "trapecio: el límite superior debe ser mayor al primer punto tabulado"
            let yc = inLin(b,j-1,j);	
            x.splice(j); x.push(b);
            y.splice(j); y.push(yc);
            // Si es uno de los puntos tabulados, se eliminan los puntos posteriores
        } else {
            let k = x.indexOf(b);
            x.splice(k+1);
            y.splice(k+1);
        }
    }
  
    // Se calcula la integral con el método del trapecio
    let s = 0;
    for (let i=0; i<x.length-1; i++) {
      s += (x[i+1]-x[i])/2*(y[i]+y[i+1]);
    }
    return s;  

};




/**
 * Calcula, empleando el método del Simpson, con los datos de la tabla desde la cual se llama al método [[x<sub>0</sub>, x<sub>1</sub>, ...], [y<sub>0</sub>, y<sub>1</sub>, ...]], la integral entre los límites de la tabla (entre x<sub>0</sub> y x<sub>n-1</sub>) o entre los límites inicial (<i>a</i>) y final (<i>b</i>) proporcionados.<br>
 * Los datos correspondientes a la variable independiente (x) deben estar en orden ascendente.
 * @method external:Array#simpson
 * @param {Object} [limites] - Objecto con los límites de integración
 * @param {number} [limites.a] - Límite inferior
 * @param {number} [limites.b] - Límite superior
 * @returns {number} Integral calculada con los datos tabulados
 * @example <caption>Ejemplos de uso del método</caption>
 * // Con los siguientes datos para las variables independiente (vx) y dependiente (vy)
 * var vx = [2.0, 2.5, 2.7, 3.0, 3.3, 3.8, 4.2, 4.5, 5.1, 5.7, 6.0, 6.4, 6.9, 7.2, 7.5];
 * var vy = [4.800, 6.625, 7.509, 9.00, 10.689, 13.944, 16.944, 19.425, 24.981, 31.329, 34.800, 39.736, 46.401, 50.664, 55.125];
 * var td = [vx, vy]; // Tabla con los datos x-y
 * 
 * // Llamando al método de Simpson, para los siguientes casos, se obtienen los resultados mostrados
 * 
 * td.simpson()
 * 134.29166666666669 // Integral entre 2.0 y 7.5
 * 
 * td.simpson({a:2.5, b:7.2})
 * 115.59493333333334 // Integral entre 2.5 y 7.2
 * 
 * td.simpson({a:2.6, b:6.7})
 * 91.35073333333335 // Integral entre 2.6 y 6.7 (límites interpolados)
 * 
 * td.simpson({a:3.1})
 * 126.64813333333333 // Integral entre 3.1 (límite inicial interpolado) y 7.5
 * 
 * td.simpson({b:7.1})
 * 113.43419999999998 // Integral entre 2.0 y 7.1 (límite final interpolado)
 * 
 * td.simpson({a:1.5, b:6.0})
 * 69.525 // Integral entre 1.5 (límite inicial extrapolado) y 6.0
 * 
 * td.simpson({a:3.0, b:7.8})
 * 144.8064 // Integral entre 3.0 y 7.8 (límite final extrapolado)
 * 
 * td.simpson({a:1.3, b:8.1})
 * 172.9081333333333 // Integral entre 1.3 y 8.1 (límites inferior y superior extrapolados)
 * 
 * _
 */
Array.prototype.simpson = function({a, b}={}) {
  
    if (this.length<2) throw "simpson: la tabla debe tener dos filas: datos 'x' y datos 'y'";
    let [x, y] = [this[0].slice(), this[1].slice()];
    if (x.length != y.length) throw "simpson: las dos filas deben tener el mismo número de elementos."
    if (x.length<3) throw "simpson: La tabla debe tener por lo menos 3 puntos";
  
    // Cálculo de los coeficientes de la ecuación cuadrática
    const coeficientes = (i,j,k) => {
        const a = ((y[i]-y[k])*(x[i]-x[j])-(y[i]-y[j])*(x[i]-x[k]))/
        ((x[i]**2-x[k]**2)*(x[i]-x[j])-(x[i]**2-x[j]**2)*(x[i]-x[k]));
        const b = ((y[i]-y[j])-a*(x[i]**2-x[j]**2))/(x[i]-x[j]);
        const c = y[i]-a*x[i]**2-b*x[i];
        return [a,b,c];
    };
  
    // Cálculo del área debajo de la curva cuadrática
    const area = (a,b,c,i,k) => a/3*(x[k]**3-x[i]**3)+b/2*(x[k]**2-x[i]**2)+c*(x[k]-x[i]);
  
    //  Si se ha mandado el límite inferior de integración
    if (a!=undefined) {
        // Y el límite es menor al menor valor tabulado, se extrapola el punto y se añade a los datos
        if (a<x[0]) {
            const [ca,cb,cc] = coeficientes(0,1,2);
            x.unshift(a);
            y.unshift(ca*a**2+cb*a+cc);
        // Si no corresponde a uno de los puntos tabulados se interpola el punto, se eliminan los puntos previos y se añade al principio el punto interpolado
        } else if (x.indexOf(a)<0) {
            let j = x.findIndex(e => e>a);
            if (j<0 || j>x.length-2) throw "simpson: el límite inferior debe ser menor al penúltimo punto tabulado"
            const [ca,cb,cc] = coeficientes(j-1,j,j+1);
            x.splice(0,j); x.unshift(a);
            y.splice(0,j); y.unshift(ca*a**2+cb*a+cc);
        // Si es uno de los puntos tabulados, se eliminan los puntos previos
        } else {
            let i = x.indexOf(a);
            x.splice(0,i);
            y.splice(0,i);
        }
    }
  
    //  Si se ha mandado el límite superior de integración
    if (b!=undefined) {
        // Y el límite es mayor al mayor valor tabulado, se extrapola el punto y se añade a los datos
        if (b>x[x.length-1]) {
            let k = x.length-1;
            const [ca,cb,cc] = coeficientes(k-2,k-1,k);
            x.push(b);
            y.push(ca*b**2+cb*b+cc);
        // Si no corresponde a uno de los puntos tabulados se interpola el punto, se eliminan los puntos posteriores y se añade al filan el punto interpolado
        } else if (x.indexOf(b)<0) {
            let j = x.findIndex(e => e>b);
            if (j<0 || j<2) throw "simpson: el límite superior debe ser mayor al segundo punto tabulado"
            const [ca,cb,cc] = coeficientes(j-2,j-1,j);
            x.splice(j); x.push(b);
            y.splice(j); y.push(ca*b**2+cb*b+cc);
        // Si es uno de los puntos tabulados, se eliminan los puntos posteriores
        } else {
            let k = x.indexOf(b);
            x.splice(k+1);
            y.splice(k+1);
        }
    }

    // Se calcula la integral con el método de Simpson
    let s=0, i=0, n=(x.length-1).quot(2);
    if (x.length%2===0) {
      const [a,b,c] = coeficientes(0,1,2);
      s += area(a,b,c,0,1);
      i++;
    }
    for (let l=0; l<n; l++) {
      const [a,b,c] = coeficientes(i,i+1,i+2);
      s += area(a,b,c,i,i+2);
      i+=2;
    }
    return s;

};