#define _CRT_SECURE_NO_WARNINGS
#ifndef __GNUC__
#define _CRTDBG_MAP_ALLOC
#include <stdlib.h>
#include <crtdbg.h>
#endif
#include "math.h"
#include "stdio.h"
#include <string.h>
#include "CoolProp.h"
#include "FloodProp.h"

	//  ****************************
	//  ***** Gas Variables ********
	//  ****************************
	static int I_N2=0, I_He=1, I_Ne=2, I_Ar=3, I_Kr=4, I_Xe=5, I_CO2=6;
	static int I_Methanol=0, I_Ethanol=1, I_Propanol=2, I_Butanol=3, I_Water=4, I_NH3=5, I_Zerol=6,I_POE=7;
	
	static double cp_A[6], cp_B[6], cp_C[6], cp_D[6], cp_E[6];
	static double kg_A[7], kg_B[7], kg_C[7];
	static double w[6], Pc[6], Tc[6], MM_g[6];
	static double mug_A[7], mug_B[7], mug_C[7];
    //When adding gas, make sure to increase
    // the length of coefficient vectors
	
	//  ****************************
	//  ***** Liq Variables ********
	//  ****************************
	
	#define NL 8

	double rhol_A[NL], rhol_B[NL], rhol_n[NL], rhol_Tc[NL];
	double kl_A[NL], kl_B[NL], kl_C[NL], kl_D[NL];
	double cl_A[NL], cl_B[NL], cl_C[NL], cl_D[NL];
	double mul_A[NL], mul_B[NL], mul_C[NL], mul_D[NL];
	double MM_l[NL];

int isNAN_FP(double x)
{
	// recommendation from http://www.devx.com/tips/Tip/42853
	return x != x;
}

int isINFINITY_FP(double x)
{
	// recommendation from http://www.devx.com/tips/Tip/42853
	if ((x == x) && ((x - x) != 0.0)) 
		return 1;//return (x < 0.0 ? -1 : 1); // This will tell you whether positive or negative infinity
	else 
		return 0;
}

int getIndex(char *Fluid)
{
    int I;
    if (!strcmp(Fluid,"N2")) {I=I_N2;}
    if (!strcmp(Fluid,"He")) {I=I_He;}
    if (!strcmp(Fluid,"Ne")) {I=I_Ne;}
    if (!strcmp(Fluid,"Ar")) {I=I_Ar;}
    if (!strcmp(Fluid,"Kr")) {I=I_Kr;}
    if (!strcmp(Fluid,"Xe")) {I=I_Xe;}
    if (!strcmp(Fluid,"CO2")) {I=I_CO2;}
    
    if (!strcmp(Fluid,"Methanol")) { I=I_Methanol;}
    if (!strcmp(Fluid,"Ethanol")) { I=I_Ethanol;}
    if (!strcmp(Fluid,"Propanol")) {I=I_Propanol;}
    if (!strcmp(Fluid,"Butanol")) {I=I_Butanol;}
    if (!strcmp(Fluid,"Water")) {I=I_Water;}
    if (!strcmp(Fluid,"NH3")) {I=I_NH3;}
    if (!strcmp(Fluid,"Zerol")) {I=I_Zerol;}
	if (!strcmp(Fluid,"POE")) {I=I_POE;}
    
    return I;
}

// // ********************************************************************
// // ************************  CONSTANTS  *******************************
// // ********************************************************************
void setGas()
{

    //Constants for Molar specific heat
    //cp in [kJ/kg-K]
        cp_A[I_N2]=29.342;
        cp_B[I_N2]=-0.0035395;
        cp_C[I_N2]=0.000010076;	
        cp_D[I_N2]=-4.3116E-09;	
        cp_E[I_N2]=2.5935E-13;
        
        cp_A[I_He]=20.786;
        cp_B[I_He]=0.0;
        cp_C[I_He]=0.0;
        cp_D[I_He]=0.0;
        cp_E[I_He]=0.0;

        cp_A[I_Ne]=20.786;
        cp_B[I_Ne]=0.0;
        cp_C[I_Ne]=0.0;
        cp_D[I_Ne]=0.0;
        cp_E[I_Ne]=0.0;

        cp_A[I_Ar]=20.786;
        cp_B[I_Ar]=0.0;
        cp_C[I_Ar]=0.0;
        cp_D[I_Ar]=0.0;
        cp_E[I_Ar]=0.0;

        cp_A[I_Kr]=20.786;
        cp_B[I_Kr]=0.0;
        cp_C[I_Kr]=0.0;
        cp_D[I_Kr]=0.0;
        cp_E[I_Kr]=0.0;

        cp_A[I_Xe]=20.786;
        cp_B[I_Xe]=0.0;
        cp_C[I_Xe]=0.0;
        cp_D[I_Xe]=0.0;
        cp_E[I_Xe]=0.0;
    //Constants for Molar specific heat
    //cp in [kJ/kg-K]
    
    // Thermal Conductivity
    // Units of [W/m-K]
        kg_A[I_N2]=0.00309;     kg_B[I_N2]=7.593e-5;        kg_C[I_N2]=-1.1014e-8;
        kg_A[I_He]=0.05516;     kg_B[I_He]=0.0003254;       kg_C[I_He]=-2.2723E-08;
        kg_A[I_Ne]=0.01379;     kg_B[I_Ne]=0.00012156;      kg_C[I_Ne]=-2.359E-08;
        kg_A[I_Ar]=0.00548;     kg_B[I_Ar]=0.000043869;     kg_C[I_Ar]=-6.8141E-09;
        kg_A[I_Kr]=0.00168;     kg_B[I_Kr]=0.000027493;     kg_C[I_Kr]=-4.7254E-09;
        kg_A[I_Xe]=0.00034;     kg_B[I_Xe]=0.000018809;     kg_C[I_Xe]=-3.0072E-09;
        kg_A[I_CO2]=-0.012;      kg_B[I_CO2]=0.00010208;     kg_C[I_CO2]=-2.2403E-08;
    //Thermal Conductivity
    //Units of [W/m-K]
    
    //Accentric Factor For Gas
    //Dimensionless
        w[I_N2]=0.04;
        w[I_He]=-0.39;
        w[I_Ne]=-0.04;
        w[I_Ar]=0.0;
        w[I_Kr]=0.0;
        w[I_Xe]=0.0;
    //Accentric Factor For Gas
    //Dimensionless

    //Critical Temp For Gas
    //Units of [K]
        Tc[I_N2]=126.1;
        Tc[I_He]=5.2;
        Tc[I_Ne]=44.4;
        Tc[I_Ar]=150.86;
        Tc[I_Kr]=209.35;
        Tc[I_Xe]=289.74;
    //Critical Temp For Gas
    //Units of [K]        

    //Critical Pressure For Gas
    //Units of [kPa]
        Pc[I_N2]=3394.0;
        Pc[I_He]=228.0;
        Pc[I_Ne]=2653.0;
        Pc[I_Ar]=4898.0;
        Pc[I_Kr]=5502.0;
        Pc[I_Xe]=5840.0;
    //Critical Pressure For Gas
    //Units of [kPa]
        
    //Viscosity Constants of Gas
    //Viscosity has units of micro-Poise
        mug_A[I_Ar]=44.997;     mug_B[I_Ar]=0.63892;	mug_C[I_Ar]=-0.00012455;
        mug_A[I_He]=71.094;     mug_B[I_He]=0.443;      mug_C[I_He]=-0.0000518;
        mug_A[I_Kr]=31.096;     mug_B[I_Kr]=0.798;      mug_C[I_Kr]=-0.000179;
        mug_A[I_Ne]=102.964;	mug_B[I_Ne]=0.746;      mug_C[I_Ne]=-0.000136;
        mug_A[I_N2]=42.606;     mug_B[I_N2]=0.475;      mug_C[I_N2]=-0.0000988;
        mug_A[I_Xe]=7.386;      mug_B[I_Xe]=0.787;      mug_C[I_Xe]=-0.000151;
        mug_A[I_CO2]=11.811;    mug_B[I_CO2]=0.49838;   mug_C[I_CO2]=-0.00010851;
    //Viscosity Constants of Gas
    //Viscosity has units of micro-Poise
 
    //Mole Mass of Gas
    //Units of kg/kmol
        MM_g[I_Ar]=39.948;
        MM_g[I_He]=4.003;
        MM_g[I_Kr]=83.8;
        MM_g[I_Ne]=20.18;
        MM_g[I_N2]=28.013;
        MM_g[I_Xe]=131.29;
    //Mole Mass of Gas
    //Units of kg/kmol

}

void setLiq()
{
    //Constants for Molar specific heat
    //cp in [kJ/kg-K]
        cl_A[I_Methanol]=40.152;
        cl_B[I_Methanol]=0.31046;
        cl_C[I_Methanol]=-0.0010291;
        cl_D[I_Methanol]=1.4598E-06;

        cl_A[I_Ethanol]=59.342;	
        cl_B[I_Ethanol]=0.36358;	
        cl_C[I_Ethanol]=-0.0012164;	
        cl_D[I_Ethanol]=0.000001803;

        cl_A[I_Propanol]=72.525;	
        cl_B[I_Propanol]=0.79553;	
        cl_C[I_Propanol]=-0.002633;	
        cl_D[I_Propanol]=3.6498E-06;

        cl_A[I_Butanol]=95.037;	
        cl_B[I_Butanol]=0.56593;
        cl_C[I_Butanol]=-0.0018256;	
        cl_D[I_Butanol]=2.6675E-06;

        cl_A[I_Water]=92.053;
        cl_B[I_Water]=-0.039953;
        cl_C[I_Water]=-0.00021103;
        cl_D[I_Water]=5.3469E-07;
        
        cl_A[I_NH3]=-182.157;
        cl_B[I_NH3]=3.3618;
        cl_C[I_NH3]=-0.014398;
        cl_D[I_NH3]=0.000020371;
        
        cl_A[I_Zerol]=337.116/1000;
        cl_B[I_Zerol]=5.186/1000;
        cl_C[I_Zerol]=0.0;
        cl_D[I_Zerol]=0.0;

		cl_A[I_POE]=1.8; // From Booser page 16
        cl_B[I_POE]=0.0;
        cl_C[I_POE]=0.0;
        cl_D[I_POE]=0.0;

    //Constants for Molar specific heat
    //cp in [kJ/kg-K]
        
    // Thermal Conductivity
    // Units of [W/m-K]        
        kl_A[I_Methanol]=-1.1793;	kl_B[I_Methanol]=0.6191;	kl_C[I_Methanol]=512.58;
        kl_A[I_Ethanol]=-1.3172;	kl_B[I_Ethanol]=0.6987;     kl_C[I_Ethanol]=516.25;
        kl_A[I_Propanol]=-1.3721;	kl_B[I_Propanol]=0.658;     kl_C[I_Propanol]=508.31;
        kl_A[I_Butanol]=-1.4633;	kl_B[I_Butanol]=0.7473;     kl_C[I_Butanol]=536.01;
        kl_A[I_Water]=-0.2758;      kl_B[I_Water]=0.004612;     kl_C[I_Water]=-5.5391E-06;
        kl_A[I_NH3]=1.1606;         kl_B[I_NH3]=-0.002284;      kl_C[I_NH3]=3.1245E-18;
        kl_A[I_Zerol]=0.1700;       kl_B[I_Zerol]=0;            kl_C[I_Zerol]=0;
		kl_A[I_POE]=0.147;          kl_B[I_POE]=0;              kl_C[I_POE]=0;
    // Thermal Conductivity
    // Units of [W/m-K]
        
    //Viscosity Constants of Liquid
    //Viscosity has units of centi-Poise
        mul_A[I_Methanol]=-9.0562;	mul_B[I_Methanol]=1254.2;	mul_C[I_Methanol]=0.022383;     mul_D[I_Methanol]=-0.000023538;
        mul_A[I_Ethanol]=-6.4406;	mul_B[I_Ethanol]=1117.6;	mul_C[I_Ethanol]=0.013721;      mul_D[I_Ethanol]=-0.000015465;
        mul_A[I_Propanol]=-0.7009;	mul_B[I_Propanol]=841.5;	mul_C[I_Propanol]=-0.0086068;	mul_D[I_Propanol]=8.2964E-06;
        mul_A[I_Butanol]=-20.6736;	mul_B[I_Butanol]=3549.3;	mul_C[I_Butanol]=0.040352;      mul_D[I_Butanol]=-0.000030937;
        mul_A[I_Water]=-10.2158;	mul_B[I_Water]=1792.5;      mul_C[I_Water]=0.01773;         mul_D[I_Water]=-0.000012631;
        mul_A[I_NH3]=-8.591;        mul_B[I_NH3]=876.4;         mul_C[I_NH3]=0.02681;           mul_D[I_NH3]=-0.00003612;        
        mul_A[I_Zerol]=0.0102;      mul_B[I_Zerol]=0;           mul_C[I_Zerol]=0;               mul_D[I_Zerol]=0;      
		mul_A[I_POE]=0.0102;        mul_B[I_POE]=0;             mul_C[I_POE]=0;                 mul_D[I_POE]=0;      
    //Viscosity Constants of Liquid
    //Viscosity has units of centi-Poise        
        
    //Density Constants of Liquid
    //Density has units of kg/m^3
        rhol_A[I_Methanol]=0.27197;	rhol_B[I_Methanol]=0.27192;	rhol_n[I_Methanol]=0.2331;	rhol_Tc[I_Methanol]=512.58;
        rhol_A[I_Ethanol]=0.2657;	rhol_B[I_Ethanol]=0.26395;	rhol_n[I_Ethanol]=0.2367;	rhol_Tc[I_Ethanol]=516.25;
        rhol_A[I_Propanol]=0.26785;	rhol_B[I_Propanol]=0.26475;	rhol_n[I_Propanol]=0.243;	rhol_Tc[I_Propanol]=508.31;
        rhol_A[I_Butanol]=0.27343;	rhol_B[I_Butanol]=0.2635;	rhol_n[I_Butanol]=0.2604;	rhol_Tc[I_Butanol]=536.01;
        rhol_A[I_Water]=0.3471;     rhol_B[I_Water]=0.274;      rhol_n[I_Water]=0.28571;	rhol_Tc[I_Water]=647.13;
        rhol_A[I_NH3]=0.23689;      rhol_B[I_NH3]=0.25471;      rhol_n[I_NH3]=0.2887;       rhol_Tc[I_NH3]=405.65;
    //Density Constants of Liquid
    //Density has units of kg/m^3
        
    //Mole Mass of Liquid
    //Units of kg/kmol
        MM_l[I_Methanol]=32.04;
        MM_l[I_Ethanol]=46.06844;
        MM_l[I_Propanol]=60.10;
        MM_l[I_Butanol]=74.1216;
        MM_l[I_Water]=18.0153;
        MM_l[I_NH3]=17.0306;
        MM_l[I_Zerol]=1;  //Dummy value since values for Zerol are given
                          // in terms of kJ/kg and dont need conversion
		MM_l[I_POE]=1;    //Dummy value since values for POE are given
                          // in terms of kJ/kg and dont need conversion
    //Mole Mass of Liquid
    //Units of kg/kmol        
}


// // ********************************************************************
// // ************************  FUNCTIONS  *******************************
// // ********************************************************************

double cK_e(double v_l, double v_g, double x, double w, double flag)
{
    // Equation taken from page 43, equation 4.51 from Chisholm for
    // liquid entrainment in gas.  Value of w 0.4 is recommended from text
    double KE;
    if (x==0 || x==1)
    {
        return 1;
    }
    
    if (flag>0.9 && flag<1.1)
    {
        KE=w+(1.0-w)*sqrt((v_g/v_l+w*(1-x)/x)/(1+w*(1-x)/x));   
    }
    if (flag>1.9 && flag<2.1)
    {    
        // Chisholms sqrt(vh/vl)
        KE=sqrt(1.0+x*(v_g/v_l-1.0));
    }
    if (flag>2.9 && flag<3.1)
    {
        KE=pow(v_g/v_l,0.25*.28);
    }
    return KE;
}

double cv_e(double v_l, double v_g, double K_e, double x, double w, double flag)
{
    double ve;
    double Kc;
    if (flag>0.9 && flag<1.1)
    {
        // using 5.48 and 5.49 from Chisholm
        // if w=0, separated flow results ( flag==2 )
        // if w=1, homogeneous flow results ( flag == 5)
        // So basically this form is general and captures all possibilities, sep, hom, or entrained
        // should use this one

        //     (    [(1-w)^2] ) -1
        // Kc= ( w+ [-------] )
        //     (    [ K_e-w ] )
        
        Kc=1.0/(w+((1.0-w)*(1.0-w))/(K_e-w));
        ve=(x*v_g+K_e*(1.0-x)*v_l)*(x+(1.0-x)/Kc);
    }
    if (flag>1.9 && flag<2.1)
    {
        // Equation 5.13 from Chisholm for separated flow
        ve=(x*v_g+K_e*(1.0-x)*v_l)*(x+(1.0-x)/K_e);
    }
    if (flag>2.9 && flag<3.1)
    {
        // Equation 2.48 from Chisholm
        ve=(1+w*(1-x)/x*v_l/v_g)/(1+w*(1-x)/x)*v_g;
    }
    if (flag>3.9 && flag<4.1)
    {
        // using 15 from Morris
        // If going to use this formula, must change pow() to multiply
        ve=(x*v_g+K_e*(1.0-x)*v_l)*(x+(1.0-x)/K_e*(1+pow(K_e-1,2.0)/(pow(v_g/v_l,0.5)-1)));
    }
    if (flag>4.9 && flag<5.1)
    {
        // homogenous
        ve=(x*v_g+(1.0-x)*v_l);
    }
    return ve;
}

double R(char *Gas)
{
    // output in kJ/kg-K
	int ii;

	return 8.31447215/Props('M','T',0,'P',0,Gas);
    
    setGas();
    ii=getIndex(Gas);
    return 8.31447215/MM_g[ii];
}

double rho_l(char *Liq, double T)
{
    // T in K
    // rho in kg/m^3
    double rhoL;
    int ii;
    setLiq();
    ii=getIndex(Liq);
    if (!strcmp(Liq,"Zerol"))
        rhoL=-.6670*T +1050.865;
	else if (!strcmp(Liq,"POE"))
	{
		// Based on 3MAF POE oil data provided by Emerson Climate
		rhoL=(-0.00074351165052847*(T-273.15)+0.992439584365375)*1000;
	}
    else
    {
        //Given by the form rho=A/B^((1-T/Tc)^n)
        rhoL=rhol_A[ii]/pow(rhol_B[ii],pow(1-T/rhol_Tc[ii],rhol_n[ii]))*1000;
    }
    return rhoL;
}

double rho_g(char *Gas, double T, double P)
{ 
    // input in K, [-]
    // output in kg/m^3
	double Value;

	// Turn on the LUT
	UseSinglePhaseLUT(1);
	Value=Props('D','T',T,'P',P,Gas);
	if (!ValidNumber(Value))
	{
		//ERROR
		printf("rho_g error");
	}
	else
		return Value;

    return P/(R(Gas)*T);
}

double h_g(char *Gas, double T, double P)
{
    // input in K,kPa
    // output in kJ/kg
    int ii;
	double Value;

	// Turn on the LUT
	UseSinglePhaseLUT(1);
    Value=Props('H','T',T,'P',P,Gas);
    if (!ValidNumber(Value))
    {
    	//ERROR
    	printf("h_g error");
    }
    else
    	return Value;

    setGas();
	ii=getIndex(Gas);
	return (cp_A[ii]*(T-298.15) + cp_B[ii]/2.0*(T*T-298.15*298.15) + cp_C[ii]/3.0*(T*T*T-298.15*298.15*298.15) + cp_D[ii]/4.0*(T*T*T*T-298.15*298.15*298.15*298.15) + cp_E[ii]/5.0*(T*T*T*T*T-298.15*298.15*298.15*298.15*298.15))/MM_g[ii];
}

double u_g(char *Gas, double T, double P)
{
    // input in K,kPa
    // output in kJ/kg
	double Value;

	// Turn on the LUT
	UseSinglePhaseLUT(1);
	Value=Props('U','T',T,'P',P,Gas);
	if (!ValidNumber(Value))
	{
		//ERROR
		printf("u_g error");
	}
	else
		return Value;
    
	return h_g(Gas,T,P)-R(Gas)*T;
}

double s_g(char *Gas, double T, double P)
{
    // input in K
    // output in kJ/kg-K
	double Value;
	int ii;

	// Turn on the LUT
	UseSinglePhaseLUT(1);
	Value=Props('S','T',T,'P',P,Gas);
	if (!ValidNumber(Value))
	{
		//ERROR
		printf("s_g error");
	}
	else
		return Value;

    setGas();
    ii=getIndex(Gas); 
    return (cp_A[ii]*log(T/298.15) + cp_B[ii]*(T-298.15) + cp_C[ii]/2.0*(T*T-298.15*298.15) + cp_D[ii]/3.0*(T*T*T-298.15*298.15*298.15))/MM_g[ii]-R(Gas)*log(P/101.325);
}

double c_v(char *Gas, double T, double P)
{ 
    // input in K, [-]
    // output in kJ/kg-K
    int ii;
    double Value;
    
    // Turn on the LUT
	UseSinglePhaseLUT(1);
	Value=Props('O','T',T,'P',P,Gas);
	if (!ValidNumber(Value))
	{
		//ERROR
		printf("c_v error");
	}
	else
		return Value;

    setGas();
    ii=getIndex(Gas);
    return (cp_A[ii] + cp_B[ii]*T + cp_C[ii]*T*T + cp_D[ii]*T*T*T + cp_E[ii]*T*T*T*T)/MM_g[ii]-R(Gas);
}

double c_p(char *Gas, double T, double P)
{ 
    // input in K, [-]
    // output in kJ/kg-K
    int ii;
    double Value;

	// Turn on the LUT
	UseSinglePhaseLUT(1);
	Value=Props('C','T',T,'P',P,Gas);
	if (!ValidNumber(Value))
	{
		//ERROR
		printf("c_p error");
	}
	else
		return Value;
    
    setGas();
    ii=getIndex(Gas);
    return (cp_A[ii] + cp_B[ii]*T + cp_C[ii]*T*T + cp_D[ii]*T*T*T + cp_E[ii]*T*T*T*T)/MM_g[ii];
}

double k_g(char *Gas, double T, double p)
{
    // input in K
    // output in kW/m-K
    int ii;
    double Value;

	// Turn on the LUT
	UseSinglePhaseLUT(1);
	Value=Props('L','T',T,'P',p,Gas);
	if (!ValidNumber(Value))
	{
		//ERROR
		printf("k_g error");
	}
	else
		return Value;

    setGas();
    ii=getIndex(Gas);
    return (kg_A[ii] + kg_B[ii]*T + kg_C[ii]*T*T)/1000.0;
}

double mu_g(char *Gas, double T, double p)
{
    // input in K
    // output in microP --> Pa-s
    int ii;
    double Value;

	// Turn on the LUT
	UseSinglePhaseLUT(1);
	Value=Props('V','T',T,'P',p,Gas);
	if (!ValidNumber(Value))
	{
		//ERROR
		printf("mu_g error");
	}
	else
		return Value;

    setGas();
    ii=getIndex(Gas);
    return (mug_A[ii] + mug_B[ii]*T + mug_C[ii]*T*T)/1.0e7;
}


double rho_m(char *Gas, char *Liq, double T, double P, double xL)
{
    // input in K, kPa, [-]
    // output in kg/m^3
    double vG, vL, x,rhom;

    if (xL==0)
    {
        return rho_g(Gas,T,P);
    }
    if (xL==1)
    {
        return rho_l(Liq,T);
    }
    vG=1/rho_g(Gas,T,P);
    vL=1/rho_l(Liq,T);
    x=1-xL;
	rhom=1/(vG*x+vL*(1-x));
	if (isNAN_FP(rhom))
		printf("rhom is NaN");
	if (isINFINITY_FP(rhom))
		printf("rhom is Infinite");
    return rhom;
}

double u_m(char *Gas, char *Liq, double T, double P, double xL)
{
    // input in K, [-]
    // output in kJ/kg
    double uG, uL,um;
        uL=u_l(Liq,T);
        uG=u_g(Gas,T,P);
		um=xL*uL+(1-xL)*uG;
	if (isNAN_FP(um))
		printf("um is NaN");
	if (isINFINITY_FP(um))
		printf("um is Infinite");
    return um;
}

double h_m(char *Gas, char *Liq, double T, double P, double xL)
{
    // input in K, [-]
    // output in kJ/kg
    double hG,hL,hm;
    hG=h_g(Gas,T,P);
    hL=u_l(Liq,T)+(P-101.325)/rho_l(Liq,T);
	hm=xL*hL+(1-xL)*hG;
	if (isNAN_FP(hm))
		printf("hm is NaN");
	if (isINFINITY_FP(hm))
		printf("hm is Infinite");
    return hm;
}



double c_l(char *Liq, double T)
{
    // input in K, [-]
    // output in kJ/kg
    int ii;
    setLiq();
    ii=getIndex(Liq);
    return (cl_A[ii] + cl_B[ii]*T + cl_C[ii]*T*T + cl_D[ii]*T*T*T)/MM_l[ii];
}

double u_l(char *Liq, double T)
{
    // input in K
    // output in kJ/kg
    int ii;
    setLiq();
    ii=getIndex(Liq);
    return (cl_A[ii]*(T-298.15) + cl_B[ii]/2.0*(T*T-298.15*298.15) + cl_C[ii]/3.0*(T*T*T-298.15*298.15*298.15) + cl_D[ii]/4.0*(T*T*T*T-298.15*298.15*298.15*298.15))/MM_l[ii];
}

double s_l(char *Liq, double T)
{
    // input in K
    // output in kJ/kg-K
    int ii;
    setLiq();
    ii=getIndex(Liq); 
    return (cl_A[ii]*log(T/298.15) + cl_B[ii]*(T-298.15) + cl_C[ii]/2.0*(T*T-298.15*298.15) + cl_D[ii]/3.0*(T*T*T-298.15*298.15*298.15))/MM_l[ii];
}

double k_l(char *Liq, double T)
{
    // input in K
    // output in kW/m-K
    int ii;
    double kL;
    setLiq();
    ii=getIndex(Liq); 
    if (!strcmp(Liq,"NH3") || !strcmp(Liq,"Water") || !strcmp(Liq,"Zerol") || !strcmp(Liq,"POE"))
    {
        kL=(kl_A[ii] + kl_B[ii]*T + kl_C[ii]*T*T)/1000.0;
    }
    else
    {
        kL=(pow(10.0,kl_A[ii] + kl_B[ii]*pow(1-T/kl_C[ii],2.0/7.0)))/1000.0;
    }
       
    return kL;
}

double mu_l(char *Liq, double T)
{
    // input in K
    // output in cP --> Pa-s
    int ii;
    double muL,mu_cSt;
    
    if (!strcmp(Liq,"Zerol"))
    {
        muL=-0.000122996*T+0.048002276;
    }
	else if (!strcmp(Liq,"POE"))
	{
		// Based on 3MAF POE oil data provided by Emerson Climate
		mu_cSt= 0.0002389593*log(T)*log(T) - 0.1927238779*log(T) + 40.3718884485;
		// From cSt to m^s, multiply by 1e-6, then multiply by density
		muL=mu_cSt*1e-6*rho_l(Liq,T);
	}
    else
    {
		setLiq();
		ii=getIndex(Liq); 
        muL=pow(10.0,mul_A[ii] + mul_B[ii]/T + mul_C[ii]*T + mul_D[ii]*T*T)/1.0e3;
    }
    return muL;
}



double dudT_m(char *Gas, char *Liq, double T, double P, double xL)
{
    double delta=.001;
    return (u_m(Gas,Liq,T+delta,P,xL)-u_m(Gas,Liq,T,P,xL))/delta;
}

double dudP_m(char *Gas, char *Liq, double T, double P, double xL)
{
    double delta=.001;
    return (u_m(Gas,Liq,T,P+delta,xL)-u_m(Gas,Liq,T,P,xL))/delta;
}

double dudxL_m(char *Gas, char *Liq, double T, double P, double xL)
{
    double delta=.001;
    return (u_m(Gas,Liq,T,P,xL+delta)-u_m(Gas,Liq,T,P,xL))/delta;
}

double drhodP_m(char *Gas, char *Liq, double T, double P, double xL)
{
    double delta=0.001;
    return (rho_m(Gas, Liq,T,P+delta,xL)-rho_m(Gas, Liq,T,P,xL))/delta;
}

double drhodT_m(char *Gas, char *Liq, double T, double P, double xL)
{
    double delta=0.001;
    return (rho_m(Gas, Liq,T+delta,P,xL)-rho_m(Gas, Liq,T,P,xL))/delta;
}

double drhodxL_m(char *Gas, char *Liq, double T, double P, double xL)
{
    double delta=.001;
    return (rho_m(Gas, Liq,T,P,xL+delta)-rho_m(Gas, Liq,T,P,xL))/delta;
}

double dvdT_m(char *Gas, char *Liq, double T, double P, double xL)
{
    double delta=.001;
    return (1/rho_m(Gas, Liq,T+delta,P,xL)-1/rho_m(Gas, Liq,T,P,xL))/delta;
}

double dvdP_m(char *Gas, char *Liq, double T, double P, double xL)
{
    double delta=.001;
    return (1/rho_m(Gas, Liq,T,P+delta,xL)-1/rho_m(Gas, Liq,T,P,xL))/delta;
}

double cp_mix(char *Gas, char *Liq, double T, double P, double xL)
{
    // input in K,[-]
    // output in kJ/kg
	double cpm;
	cpm=xL*c_l(Liq,T)+(1-xL)*c_p(Gas,T,P); 
	if (isNAN_FP(cpm))
		printf("cpm is NaN");
	if (isINFINITY_FP(cpm))
		printf("cpm is Infinite");
    return cpm;
}

double mu_mix(char *Gas, char *Liq, double T, double p, double xL)
{
    double muL, muG,mum;
    muL=mu_l(Liq,T); //kg/m-s
    muG=mu_g(Gas,T,p);
	mum=1/(xL/muL+(1-xL)/muG);
	if (isNAN_FP(mum))
		printf("mum is NaN");
	if (isINFINITY_FP(mum))
		printf("mum is Infinite");
    return mum;
}

double k_mix(char *Gas, char *Liq, double T, double P, double xL)
{
    //input in K,kPa, [-]
    //output in kW/m-K
    double kL, kG, VF,km;

    kL=k_l(Liq,T); 
    kG=k_g(Gas,T,P);
    VF=VoidFrac(Gas,Liq,T,P,xL);
    if (VF<1e-12)
		return kL;
	else if (VF>1-1e-12)
		return kG;
	km=(1-VF)*kL+VF*kG;
	if (isNAN_FP(km))
		printf("km is NaN");
	if (isINFINITY_FP(km))
		printf("km is Infinite");
    return km;
}

double VoidFrac(char *Gas, char *Liq, double T, double P, double xL)
{
    //input in K, kPa, [-]
    // output in [-]
    double rhoL, rhoG,VF;
    rhoL=rho_l(Liq,T);
    rhoG=rho_g(Gas,T,P);
	VF= ((1-xL)/rhoG)/((1-xL)/rhoG+xL/rhoL);
	if (isNAN_FP(VF))
		printf("VF is NaN");
	if (isINFINITY_FP(VF))
		printf("VF is Infinite");
    return VF;
}

double Pr_mix(char *Gas, char *Liq, double T, double P, double xL)
{
	double Prm;
//	double cp,mu,k;
//	cp=cp_mix(Gas,Liq,T,P,xL);
//	mu=mu_mix(Gas,Liq,T,P,xL);
//	k=k_mix(Gas,Liq,T,P,xL);
	Prm=cp_mix(Gas,Liq,T,P,xL)*mu_mix(Gas,Liq,T,P,xL)/k_mix(Gas,Liq,T,P,xL);
	if (isNAN_FP(Prm))
		printf("Prm is NaN");
	if (isINFINITY_FP(Prm))
		printf("Prm is Infinite");

    return Prm;
}

double kstar_m(char *Gas, char *Liq, double T, double P,double xL)
{
	double kstarm;
	//double cp,cv,cl;
	//cp=c_p(Gas,T,P);
	//cv=c_v(Gas,T,P);
	//cl=c_l(Liq,T);
	kstarm=((1-xL)*c_p(Gas,T,P)+xL*c_l(Liq,T))/((1-xL)*c_v(Gas,T,P)+xL*c_l(Liq,T));
	if (isNAN_FP(kstarm))
		printf("kstarm is NaN");
	if (isINFINITY_FP(kstarm))
		printf("kstarm is Infinite");
    return kstarm;
}

double s_m(char *Gas, char *Liq, double T, double P, double xL)
{
//     input in K, kPa, [-]
//     output in kJ/kg-K
    double sG,sL,sm;
    sG=s_g(Gas,T,P);
    sL=s_l(Liq,T);
	sm=xL*sL+(1-xL)*sG;

	if (isNAN_FP(sm))
		printf("sm is NaN");
	if (isINFINITY_FP(sm))
		printf("sm is Infinite");
    return sm;
}

double e_m(char *Gas, char *Liq, double T, double P, double xL)
{
//     input in K, kPa, [-]
//     output in kJ/kg
    double T0=300,P0=100;
    return (h_m(Gas,Liq,T,P,xL)-h_m(Gas,Liq,T0,P0,xL))-T0*(s_m(Gas,Liq,T,P,xL)-s_m(Gas,Liq,T0,P0,xL));
}

double T_hp_FP(char *Gas, char *Liq, double h, double p, double xL, double T_guess)
{
	double x1,x2,x3,y1,y2,eps,f,T;
	int iter;

	eps=1e-8;
	iter=1;
	f=100.0;
	while ((iter<=3 || fabs(f)>eps) && iter<100)
	{
		if (iter==1){x1=T_guess; T=x1;}
		if (iter==2){x2=T_guess+0.1; T=x2;}
		if (iter>2) {T=x2;}

			f=h_m(Gas,Liq,T,p,xL)-h;

		if (iter==1){y1=f;}
		if (iter==2){y2=f;}
		if (iter>2)
		{
			y2=f;
			x3=x2-y2/(y2-y1)*(x2-x1);
			y1=y2; x1=x2; x2=x3;
		}
		iter=iter+1;
	}
	if (isNAN_FP(T))
		printf("uhoh");
	if (isINFINITY_FP(T))
		printf("uhoh");
	//printf("THP: %d %g %g %g %g\n",iter,f,T,p,h_m(Gas,Liq,T,p+20,xL));
	return T;
	
}

double T_Up(char *Gas, char *Liq, double U, double p, double xL, double V, double T_guess)
{
	double x1,x2,x3,y1,y2,eps,f=999.,T;
	int iter;

	eps=1e-8;
	iter=1;
	while ((iter<=3 || fabs(f)>eps) && iter<100)
	{
		if (iter==1){x1=T_guess; T=x1;}
		if (iter==2){x2=T_guess+0.1; T=x2;}
		if (iter>2) {T=x2;}

			f=U-u_m(Gas,Liq,T,p,xL)*rho_m(Gas,Liq,T,p,xL)*V;

		if (iter==1){y1=f;}
		if (iter>1)
		{
			y2=f;
			x3=x2-y2/(y2-y1)*(x2-x1);
			y1=y2; x1=x2; x2=x3;
		}
		iter=iter+1;
		if (iter>50)
		{
			printf("T_Up not converging with inputs\n U: %g p: %g xL: %g V: %g T_guess: %g\n",U,p,xL,V,T_guess);
		}
	}
	if (isNAN_FP(T))
		printf("uhoh");
	if (isINFINITY_FP(T))
		printf("uhoh");
	return T;
}

double p_Trho(char *Gas, char *Liq, double rho, double T, double xL, double p_guess)
{
	double x1,x2,x3,y1,y2,eps,f,p;
	int iter;

	eps=1e-8;
	iter=1;
    f=999;
	while ((iter<=3 || fabs(f)>eps) && iter<100)
	{
		if (iter==1){x1=p_guess; p=x1;}
		if (iter==2){x2=p_guess+0.1; p=x2;}
		if (iter>2) {p=x2;}

			// Find the pressure which gives the same density
			f=rho_m(Gas,Liq,T,p,xL)-rho;

		if (iter==1){y1=f;}
		if (iter==2){y2=f;}
		if (iter>2)
		{
			y2=f;
			x3=x2-y2/(y2-y1)*(x2-x1);
			y1=y2; x1=x2; x2=x3;
		}
		iter=iter+1;
		if (iter>50)
		{
			printf("p_Trho not converging with inputs\n rho: %g T: %g xL: %g p_guess: %g\n",rho,T,xL,p_guess);
		}
	}
	if (isNAN_FP(x3))
		printf("uhoh");
	if (isINFINITY_FP(x3))
		printf("uhoh");
	return x3;
}

double h_sp_FP(char *Gas, char *Liq, double s, double p, double xL, double T_guess)
{
	double x1,x2,x3,y1,y2,eps,f,T;
	int iter;

	eps=1e-8;
	f=999;
	iter=1;
	while ((iter<=3 || fabs(f)>eps) && iter<100)
	{
		if (iter==1){x1=T_guess; T=x1;}
		if (iter==2){x2=T_guess+0.1; T=x2;}
		if (iter>2) {T=x2;}

			// Find the temperature which gives the same entropy
			f=s_m(Gas,Liq,T,p,xL)-s;

		if (iter==1){y1=f;}
		if (iter==2){y2=f;}
		if (iter>2)
		{
			y2=f;
			x3=x2-y2/(y2-y1)*(x2-x1);
			y1=y2; x1=x2; x2=x3;
		}
		iter=iter+1;
		if (iter>50)
		{
			printf("h_sp not converging with inputs\n s: %g p: %g xL: %gT_guess: %g\n",s,p,xL,T_guess);
		}
	}
	if (isNAN_FP(h_m(Gas,Liq,T,p,xL)))
		printf("uhoh");
	if (isINFINITY_FP(h_m(Gas,Liq,T,p,xL)))
		printf("uhoh");

	//Evaluate the enthalpy at the constant-entropy temp
	return h_m(Gas,Liq,T,p,xL);
}

double T_sp(char *Gas, char *Liq, double s, double p, double xL, double T_guess)
{
	double x1,x2,x3,y1,y2,eps,f,T;
	int iter;

	eps=1e-8;
	f=999;
	iter=1;
	while ((iter<=3 || fabs(f)>eps) && iter<100)
	{
		if (iter==1){x1=T_guess; T=x1;}
		if (iter==2){x2=T_guess+0.1; T=x2;}
		if (iter>2) {T=x2;}

			// Find the temperature which gives the same entropy
			f=s_m(Gas,Liq,T,p,xL)-s;

		if (iter==1){y1=f;}
		if (iter==2){y2=f;}
		if (iter>2)
		{
			y2=f;
			x3=x2-y2/(y2-y1)*(x2-x1);
			y1=y2; x1=x2; x2=x3;
		}
		iter=iter+1;
		if (iter>50)
		{
			printf("h_sp not converging with inputs\n s: %g p: %g xL: %gT_guess: %g\n",s,p,xL,T_guess);
		}
	}
	if (isNAN_FP(h_m(Gas,Liq,T,p,xL)))
		printf("uhoh");
	if (isINFINITY_FP(h_m(Gas,Liq,T,p,xL)))
		printf("uhoh");

	//Return the temperature
	return T;
}

double dpdT_const_v(char *Gas, char *Liq, double T, double p1, double xL)
{
	double x1,x2,x3,y1,y2,eps,f,v1,delta,p2;
	int iter;

	delta=1e-5;
	eps=1e-6;
	iter=1;

	v1=1/rho_m(Gas,Liq,T,p1,xL);
	while ((iter<=3 || fabs(f)>eps) && iter<100)
	{
		if (iter==1){x1=p1; p2=x1;}
		if (iter==2){x2=p1+0.001; p2=x2;}
		if (iter>2) {p2=x2;}

			// Find the pressure which gives the same specific volume
			f=1.0/rho_m(Gas,Liq,T+delta,p2,xL)-v1;

		if (iter==1){y1=f;}
		if (iter>1)
		{
			y2=f;
			x3=x2-y2/(y2-y1)*(x2-x1);
			y1=y2; x1=x2; x2=x3;
		}
		iter=iter+1;
		if (iter>50)
		{
			printf("dpdT_const_v not converging with inputs\n T: %g p1: %g xL: %g\n",T,p1,xL);
		}
	}

	if (isNAN_FP((p2-p1)/delta))
		printf("uhoh");
	if (isINFINITY_FP((p2-p1)/delta))
		printf("uhoh");

	//Evaluate the enthalpy at the constant-entropy temp
	return (p2-p1)/delta;

}

//double c_g(char *Gas, double T, double p)
//{
//	double dp_dv,k,v;
//
//	v=1/rho_g(Gas,T,p);
//	k=c_p(Gas,T,p)/c_v(Gas,T,p);
//	dp_dv=(
//}
