#include <iostream>
#include <Vector>
#include <math.h>

#define OUTPUT
using namespace std;

const static double R1 = 0.008314472;
const static double R = 8314.472;
const static double Rgas = 0.0831434;
const static double cal2J = 4.1868;
const static double dp = 1.0;
const static double dv = 1.0e-7;
const static double Rcal = R/cal2J;
const static double pbar = 1.01325e+5;

class properties
{
    double Mw_;

    double Pc_;
    double rhoc_;
    double Tc_;
  
    double Pref_;
    double Tref_;
    double h0_;
    double s0_;

    const double *a_;
    const double *b_;
    const double *c_;
    const double *d_;
    const double *theta_;

    const double *Nk_;
    const double *ik_;
    const double *jk_;
    const double *lk_;
    
    const double *eta0_;
    const double *etab_;
    const double *etat_;

public:

  properties
  (
      const double& Mw,
      const double& Pc,
      const double& rhoc,
      const double& Tc,
      const double& Pref,
      const double& Tref,
      const double& h0,
      const double& s0,
      const double *a,
      const double *b,
      const double *c,
      const double *d,
      const double *theta,
      const double *Nk,
      const double *ik,
      const double *jk,
      const double *lk,
      const double *eta0,
      const double *etab,
      const double *etat
  )
  :
      Mw_(Mw),
      Pc_(Pc),
      rhoc_(rhoc),
      Tc_(Tc),
      Pref_(Pref),
      Tref_(Tref),
      h0_(h0),
      s0_(s0),
      a_(a),
      b_(b),
      c_(c),
      d_(d),
      theta_(theta),
      Nk_(Nk),
      ik_(ik),
      jk_(jk),
      lk_(lk),
      eta0_(eta0),
      etab_(etab),
      etat_(etat)
  {
    /*
    for(int i=0;i<23;i++)
      {
	cout << "i = " << i+1
	     << ", Nk = " << Nk_[i]
	     << ", ik = " << ik_[i]
	     << endl;
      }
    */
  }


  // access
    const double* a() const
    {
        return a_;
    }

    const double* b() const
    {
        return b_;
    }

    const double* c() const
    {
        return c_;
    }

    const double* d() const
    {
        return d_;
    }

    const double Mw() const
    {
        return Mw_;
    }

    const double Pc() const
    {
        return Pc_;
    }

    const double Tc() const
    {
        return Tc_;
    }

    const double rhoc() const
    {
        return rhoc_;
    }

    double mu
    (
        const double& d,
        const double& t
    ) const
    {
        const double ek = 362.6;
        const double sigma = 0.453e-9;
        const double Na = 6.0221415e+23;
	const double c1 = 23.7222995;
	const double c2 = -3.38264462;
	const double c3 = 12.7568864;

        double etaa_[4][3];
	etaa_[2][0] = 0.131194057;
        etaa_[2][1] = -0.382240694;
	etaa_[2][2] = 0.0;
	etaa_[3][0] = -0.0805700894;
	etaa_[3][1] = 0.153811778;
	etaa_[3][2] = -0.110578307;

        double tau = 1.0/t;
	double T = Tc_/t;
	double rho = d*rhoc_;
	double Ts = T/ek;
	double Bs = 0.0;
	for(int i=0; i<9; i++)
	{
            Bs += etab_[i]*pow(Ts, etat_[i]);
	}
	double B = Na*pow(sigma, 3.0)*Bs;
	double eta0 = eta0_[0] + eta0_[1]*T + eta0_[2]*T*T;
	double d0 = c2 + c3*sqrt(tau);
	double s = 0.0;
	for(int j=2; j<4; j++)
        {
	    for(int k=0; k<3; k++)
            {
	        s += etaa_[j][k]*pow(d, j)*pow(t, k);  
	    }
	}

	s += c1*d*(1.0/(d0 - d) - 1.0/d0);
	s *= 1.0e+3;
	return eta0*(1.0 + B*rho) + s;
    }

    double Pv
    (
        const double& T
    ) const
    {
        double p = 1.0 - T/Tc_;
	p = max(p, 0.0);
	double sum = a_[0]*pow(p, 0.5) + a_[1]*p + a_[2]*pow(p, 3.0) + a_[3]*pow(p, 11.0/2.0);

	return Pc_*exp(Tc_*sum/T);
    }

    double rholSat
    (
        const double& T
    ) const
    {
        double p = 1.0 - T/Tc_;
	p = max(p, 0.0);
	double sum = 1.0 + b_[0]*pow(p, 1.0/3.0) + b_[1]*p + b_[2]*pow(p, 5.0/3.0) + b_[3]*pow(p, 2.0) + b_[4]*pow(p, 8.0/3.0);
	
	return rhoc_*sum;

    }

    double rhovSat
    (
        const double& T
    ) const
    {
        double p = 1.0 - T/Tc_;

	if (p > 0.0)
	{
	    double sum = c_[0]*pow(p, 2.0/3.0) + c_[1]*pow(p, 5.0/3.0) + c_[2]*pow(p, 10.0/3.0) + c_[3]*pow(p, 16.0/3.0);	
	    return rhoc_*exp(sum);
	}
	else
	{
	    return rholSat(T);
	}
    }


    double cp0
    (
        const double& t
    ) const
    {

        if (Tc_/t > 1500.0)
	{
	  cout << "Warning!!! Temperature is out of valid range 1500 K, T = " << Tc_/t << endl;
	}

        double cpt = d_[0];
	
	for(int i=1; i<5; i++)
	{
	    double q = theta_[i]*t/Tc_;
	    double e = exp(q);
	    double nom = q*q*e;
	    double den = (e - 1.0)*(e - 1.0);
	    cpt += d_[i]*nom/den;
	}
	
	return cpt*R;
    }

    double integrateCp
    (
        const double& t,
	const double& t0,
	double& cp1,
	double& cp2
	
    ) const
    {
        double N = 1.0e+5;
	double dt = 0.5*(t-t0)/N;
	double cpa = cp0(t0);
	double cpb = cp0(t);
	cp1 = (cpa/t0 + cpb/t);
	cp2 = (cpa/(t0*t0) + cpb/(t*t));
	for(int i=2; i<2*N; i += 2)
	{
	    double tau = t0 + i*dt;
            double cpt = cp0(tau);
	    cp1 += 4.0*cpt/tau;
	    cp2 += 4.0*cpt/(tau*tau);
	}

	for(int i=1; i<2*N; i += 2)
	{
	    double tau = t0 + i*dt;
            double cpt = cp0(tau);
	    cp1 += 2.0*cpt/tau;
	    cp2 += 2.0*cpt/(tau*tau);
	}

	cp1 *= dt/3.0;
	cp2 *= dt/3.0;
	return 0.0;
    }

    double a0
    (
        const double& d,
	const double& t
    ) const
    {
	 
	double rho0 = Pref_/(R*Tref_);
	double d0 = rho0/rhoc_;
	double t0 = Tc_/Tref_;

	double a1 = h0_*t/Tc_;
	double a2 = 0.0;
	double a3 = 0.0;
	double dummy = integrateCp(t, t0, a2, a3);

	return (a1 - s0_ + a2 - t*a3)/R - 1.0 + log(d*t0/(d0*t));
    } 

    double aRes
    (
        const double& d,
	const double& t
    ) const
    {
	double sum = 0.0;
	for (int i=0; i<23; i++)
	{
	    double gamma = 1.0;
	    if (fabs(lk_[i]) < 1.0e-8)
	    {
	        gamma = 0.0;
	    }
            sum += Nk_[i]*pow(d, ik_[i])*pow(t, jk_[i])*exp(-gamma*pow(d, lk_[i]));
	}

	return sum;
    }

    double daResdd
    (
        const double& d,
	const double& t
    ) const
    {
	double sum = 0.0;
	for (int i=0; i<23; i++)
	{
	    double gamma = 1.0;
	    if (fabs(lk_[i]) < 1.0e-8)
	    {
	        gamma = 0.0;
	    }

	    double a = Nk_[i]*pow(t, jk_[i])*exp(-gamma*pow(d, lk_[i]));
	    double da = -a*gamma*lk_[i]*pow(d, lk_[i] - 1.0);
	    double b = pow(d, ik_[i]);
	    double db = ik_[i]*pow(d, ik_[i] - 1.0);
	    sum += a*db + b*da;
	}

	return sum;
    }

    double daResdt
    (
        const double& d,
	const double& t
    ) const
    {
	double sum = 0.0;
	for (int i=0; i<23; i++)
	{
	    double gamma = 1.0;
	    if (fabs(lk_[i]) < 1.0e-8)
	    {
	        gamma = 0.0;
	    }
            sum += jk_[i]*Nk_[i]*pow(d, ik_[i])*pow(t, jk_[i] - 1.0)*exp(-gamma*pow(d, lk_[i]));
	}

	return sum;
    }

    double da0dt
    (
        const double& d,
	const double& t
    ) const
    {
        double dt = 1.0e-4*t;
	double t1 = t + dt;
	double t0 = t - dt;
	double a1 = a0(d, t1);
	double a00 = a0(d, t0);
	
	double der = 0.5*(a1 - a00)/dt;
	//cout << "da0dt = " << der << endl;
	return der;
    }

    double d2a0dt2
    (
        const double& d,
	const double& t
    ) const
    {
        double dt = 1.0e-4*t;

	double t1 = t + dt;
	double t0 = t - dt;
	
	double a2 = a0(d, t1);
	double a1 = a0(d, t);
	double a00 = a0(d, t0);

	double result = (a2 - 2.0*a1 + a00)/(dt*dt);
	/*
	cout << "d2a0dt2 = " << result 
	     << ", dt = " << dt
	     << endl;
	*/
	return result;
    }

    double d2aResdt2
    (
        const double& d,
	const double& t
    ) const
    {
	double sum = 0.0;
	for (int i=0; i<23; i++)
	{
	    double gamma = 1.0;
	    if (fabs(lk_[i]) < 1.0e-8)
	    {
	        gamma = 0.0;
	    }

	    double term =jk_[i]*(jk_[i] - 1.0)*Nk_[i]*pow(d, ik_[i])*pow(t, jk_[i] - 2.0)*exp(-gamma*pow(d, lk_[i]));
	    /*
	    cout << "i = " << i 
		 << ", jk = " << jk_[i] 
		 << ", term = " << term
		 << endl;
	    */
            sum += term;
	}

	return sum;

    }

    double d2aResdd2
    (
        const double& d,
	const double& t
    ) const
    {
	double sum = 0.0;
	for (int i=0; i<23; i++)
	{
	    double gamma = 1.0;
	    if (fabs(lk_[i]) < 1.0e-8)
	    {
	        gamma = 0.0;
	    }

	    double a = Nk_[i]*pow(t, jk_[i])*exp(-gamma*pow(d, lk_[i]));
	    double da = -a*gamma*lk_[i]*pow(d, lk_[i]-1.0);
	    double dda = -da*gamma*lk_[i]*pow(d, lk_[i]-1.0) - a*gamma*lk_[i]*(lk_[i]-1.0)*pow(d, lk_[i]-2.0);
	    double b = pow(d, ik_[i]);
	    double db = ik_[i]*pow(d, ik_[i]-1.0);
	    double ddb = ik_[i]*(ik_[i]-1.0)*pow(d, ik_[i]-2.0);

	    sum += a*ddb + da*db + db*da + b*dda;

	}

	return sum;

	// ------------------- old ***************************
        double dd = 1.0e-5*d;

	double d1 = d + dd;
	double d0 = d - dd;
	
	double a2 = aRes(d1, t);
	double a1 = aRes(d, t);
	double a00 = aRes(d0, t);

	double result = (a2 - 2.0*a1 + a00)/(dd*dd);

	cout << "d2aResdd2 = " << result 
	     << ", = " << sum
	     << endl;

	return result;
    }

    double d2aResdtdd
    (
        const double& d,
	const double& t
    ) const
    {
	double sum = 0.0;
	for (int i=0; i<23; i++)
	{
	    double gamma = 1.0;
	    if (fabs(lk_[i]) < 1.0e-8)
	    {
	        gamma = 0.0;
	    }

	    double a = Nk_[i]*pow(t, jk_[i])*exp(-gamma*pow(d, lk_[i]));
	    //double da = -a*gamma*lk_[i]*pow(d, lk_[i]-1.0);
	    double b = pow(d, ik_[i]);
	    double db = ik_[i]*pow(d, ik_[i]-1.0);
	    double dat = Nk_[i]*jk_[i]*pow(t, jk_[i]-1.0)*exp(-gamma*pow(d, lk_[i]));
	    double dadt = -dat*gamma*lk_[i]*pow(d, lk_[i]-1.0);
	    sum += dat*db + b*dadt;

            //sum += Nk_[i]*jk_[i]*pow(t, jk_[i]-1.0)*exp(-gamma*pow(d, lk_[i]))*(ik_[i]*pow(d, ik_[i]-1.0) - gamma*lk_[i]*pow(d, ik_[i]+lk_[i]-1.0));
	}

	return sum;
    }

    double u
    (
        const double& d,
	const double& t
    ) const
    {
	return R*Tc_*(da0dt(d, t) + daResdt(d, t));
    }

    double h
    (
        const double& d,
	const double& t
    ) const
    {
	double T = Tc_/t;
        double dt = 1.0e-7*t;

	return u(d, t) + R*T*(d*daResdd(d, t) + 1.0);
    }

    double rho
    (
        const double& P,
	const double& T
    ) const
    {

	double t = Tc_/T;
        double r = 0.0;
	double pv = Pv(T);

	// initial guess for rho

	if (P > pv)
	{
	    r = rholSat(T);
	}
	else
	{
	    r = rhovSat(T);
	}

	double pq = P/(R*T);
	double err = 1.0;
	double tol = 1.0e-12;
	int i = 0;
	int Nmax = 1000;
	double ur = 1.0e-4;

	while ((err > tol) && (i<Nmax))
	{
	    double delta = r/rhoc_;
	    double A = daResdd(delta, t)/rhoc_;
	    double B = d2aResdd2(delta, t)/rhoc_;

	    double d = (pq - A*r*r - r)/(r*r*B + 2.0*A*r + 1.0);
	    err = fabs(d)/r;
	    r = r + 0.5*d;
	    i++;
	    //cout << "i = " << i << ", B = " << B << endl;
	}

	/*
	while ((err > tol) && (i<Nmax))
	{
	    double delta = r/rhoc_;
	    double den = 1.0 + delta*daResdd(delta, t);
	    double rNew = pq/den;
	    err = fabs(r-rNew)/r;
	    if (err < 1.0e+4*tol)
	    {
	      //ur *= 1.1;
	      ur = min(ur, 0.1);
	    }
	    r = r + ur*(rNew - r);
	    i++;
	}
	*/
	
	if (i > Nmax-2)
	{
	  cout << "Warning!!! rho calculation did not converge. err = " << err << endl;
	}
	return r;

    }

    double cv
    (
        const double& d,
	const double& t
    ) const
    {
	return -R*t*t*(d2a0dt2(d, t) + d2aResdt2(d, t));
    }

    double cp
    (
        const double& d,
	const double& t
    ) const
    {
	double help1 = daResdd(d, t);
	double nom = 1.0 + d*help1 - d*t*d2aResdtdd(d, t);
	double den = 1.0 + 2.0*d*help1 + d*d*d2aResdd2(d, t);

	return cv(d, t) + R*nom*nom/den;
    }

    double c
    (
        const double& d,
	const double& t
    ) const
    {
	double T = Tc_/t;
	double den = 1.0 + 2.0*d*daResdd(d, t) + d*d*d2aResdd2(d, t);
	double w2 = R*T*(cp(d, t)/cv(d, t))*den/Mw_;

	return sqrt(w2);
    }
};

int main()
{

    // ethanol property coefficients - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    const double a[] = {-0.0514770863682, -8.27075196981, -5.49245388577, 5.64828891406};
    const double b[] = {2.22194636037, -0.0469267553094, 10.3035666311, -17.2304684516, 8.23564165285};
    const double c[] = {-8.35647816638, -2.38721859682, -39.6946441837, -9.99133502692};
    const double d[] = { 6.41129104405, 1.95988750679, 7.60084166080, 3.89583440622, 4.23238091363 };
    const double theta[] = { 0.0, 694.0, 1549.0, 2911.0, 4659.0};

    const double Nk[] = { 1.14008942201e+1, -3.95227128302e+1, 4.13063408370e+1, -1.88892923721e+1,
			  4.72310314140, -7.78322827052e-3, 1.71707850032e-1, -1.53758307602,
			  1.42405508571, 1.32732097050e-1, -1.14231649761e-1, 3.27686088736e-6,
			  4.95699527725e-4, -7.01090149558e-5, -2.25019381648e-6, -2.55406026981e-1,
			  -6.32036870646e-2, -3.14882729522e-2, 2.56187828185e-2, -3.08694499382e-2,
			  7.22046283076e-3, 2.99286406225e-3, 9.72795913095e-4};

    const double ik[] = { 1.0, 1.0, 1.0, 1.0,
			  1.0, 1.0, 2.0, 2.0,
			  2.0, 3.0, 3.0, 6.0,
			  7.0, 8.0, 8.0, 1.0,
			  3.0, 3.0, 6.0, 7.0,
			  8.0, 2.0, 7.0};

    const double jk[] = { -0.5, 0.0, 0.5, 1.5,
			  2.0, 5.0, -0.5, 1.0,
			  2.0, 0.0, 2.5, 6.0,
			  2.0, 2.0, 4.0, 5.0,
			  3.0, 7.0, 5.5, 4.0,
			  1.0, 22.0, 23.0};

    const double lk[] = { 0.0, 0.0, 0.0, 0.0,
			  0.0, 0.0, 0.0, 0.0,
			  0.0, 0.0, 0.0, 0.0,
			  0.0, 0.0, 0.0, 2.0,
			  2.0, 2.0, 2.0, 2.0,
			  2.0, 4.0, 4.0};

    const double eta0[] = {-1.03116, 3.48379e-2, -6.50264e-6};
    const double etab[] = {-19.572881, 219.73999, -1015.3226, 2471.0125, -3375.1717, 2491.6597, -787.26086, 14.085455, -0.34664158};
    const double etat[] = {0.0, -0.25, -0.50, -0.75, -1.0, -1.25, -1.50, -2.50, -5.50 };

    const properties ethanol(46.06844, 6.148e+6, 5.991, 513.9, 0.001e+6, 273.15, 45800.0e+3, 180.0e+3, a, b, c, d, theta, Nk, ik, jk, lk, eta0, etab, etat);

    //const int Np = 22;
    const int Np = 2;
    //bar
    //const double p[] = { 0.1, 0.3, 0.5, 0.7, 0.9, 1.1, 1.3, 1.5, 1.7, 1.9, 2.1, 2.3, 2.5, 2.7, 2.9, 3.1, 3.3, 3.5, 3.7, 3.9, 4.1, 4.3 };
    //const double p[] = { 0.0, 200.0, 400.0, 600.0, 800.0, 1000.0, 1200.0, 1400.0, 1600.0, 1800.0, 2000.0 };
    const double p[] = { 1.0, 10.0 };

    const int NT = 3;
    // deg C
    const double T[] = { 80.0, 100.0, 120.0 };
    //const double T[] = { 300, 350, 400 };

    double rho[Np][NT];
    double soundSpeed[Np][NT];
    double mu[Np][NT];

    cout.precision(8);
    
    for(int i=0; i<Np; i++)
    {
#ifdef OUTPUT
        cout << "Calculating for p = " << p[i] 
	     << " bar. " 
	     << i+1 
	     << "/" 
	     << Np
	     << ". "
	     << flush;
#endif
        for(int j=0; j<NT; j++)
	{
#ifdef OUTPUT
	    cout << " " << T[j] << flush;
#endif
	    double pressure = (p[i] + 0.0)*1.0e+5;
	    double temp = T[j] + 273.15;
	    //double temp = T[j];
	    rho[i][j] = ethanol.rho(pressure, temp);
	    
	    double d = rho[i][j]/ethanol.rhoc();
	    double t = ethanol.Tc()/temp;
	    
	    soundSpeed[i][j] = ethanol.c(d, t);
	    mu[i][j] = ethanol.mu(d, t);
	  
	}
#ifdef OUTPUT
	cout << endl;
#endif
    }
#ifdef OUTPUT
    cout << endl;
#endif

    // -------------------------------------------------------
    cout.width(15);
    cout << "rho (kg/m3)";
    for(int i=0; i<Np; i++)
    {
        cout.width(15);
	cout << p[i];	
    }
    cout << endl;

    for(int j=0; j<NT; j++)
    {
        cout.width(15);
        cout << T[j];
        for(int i=0; i<Np; i++)
	{
	    cout.width(15);
	    cout << rho[i][j]*ethanol.Mw();
	}
	cout << endl;
    }
    cout << endl;

    // -------------------------------------------------------
    cout.width(15);
    cout << "vel (m/s)";
    for(int i=0; i<Np; i++)
    {
        cout.width(15);
	cout << p[i];	
    }
    cout << endl;

    for(int j=0; j<NT; j++)
    {
        cout.width(15);
        cout << T[j];
        for(int i=0; i<Np; i++)
	{
	    cout.width(15);
	    cout << soundSpeed[i][j];
	}
	cout << endl;
    }
    cout << endl;

    // -------------------------------------------------------
    cout.width(15);
    cout << "bulk N/mm2";	
    for(int i=0; i<Np; i++)
    {
        cout.width(15);
	cout << p[i];	
    }
    cout << endl;

    for(int j=0; j<NT; j++)
    {
        cout.width(15);
        cout << T[j];
        for(int i=0; i<Np; i++)
	{
	  double K = rho[i][j]*ethanol.Mw()*pow(soundSpeed[i][j], 2)/1.0e+6;
	    cout.width(15);
	    cout << K;
	}
	cout << endl;
    }
    cout << endl;

    // -------------------------------------------------------
    cout.width(15);
    cout << "nu mm2/s";	
    for(int i=0; i<Np; i++)
    {
        cout.width(15);
	cout << p[i];	
    }
    cout << endl;

    for(int j=0; j<NT; j++)
    {
        cout.width(15);
        cout << T[j];
        for(int i=0; i<Np; i++)
	{
	  double nu = mu[i][j]/(rho[i][j]*ethanol.Mw());
	    cout.width(15);
	    cout << nu;
	}
	cout << endl;
    }
    cout << endl;

    //    cout << "Done!" << endl;
    return 0;
}