/*
Copyright (C) 1996-1998 Evgenii Rudnyi, http://Evgenii.Rudnyi.Ru/

Software for paper 
E. B. Rudnyi, V. V. Kuzmenko, G. V. Voronin. Simultaneous assessment of the YBa2Cu6O6+z thermodynamics under the linear error model. J. Phys. Chem. Ref. Data, 1998, v. 27, N 5, p. 855-888.

This software is a copyrighted work licensed under the terms, described
in the file "FREE_LICENSE". 
*/

#include <sumsqr.h>
#include <iomanip.h>
#include <imsl.h>
#include <float.h>
#include "cuox_pl.h"

char* LICENSE = "\
Evgenii Rudnyi 1996, (c) All rights reserved \n\n\
Chemistry Department \n\
Moscow State University \n\
119899 Moscow, Russia \n\
tel. (095)939-5452 \n\
rudnyi@comp.chem.msu.su \n\n\
This program is freeware (public domain). Feel free to use and \n\
distribute it, provided no charge is taken. \n\n\
I will be glad if you like this program. Let me know if you find any \n\
bugs. I would also appreciate your comments. \n\n\
If you need to change the model - try to contact me. \n\n\
Disclaimer of warranty: \n\
This program is supplied as is. I disclaim all warranties, \n\
express or implied, including, without limitation, the warranties of \n\
merchantability and of fitness of this program for any purpose. I assume \n\
no liability for damages direct or consequential, which may result from \n\
the use of this program.";

int nparam;
param par[50];
CuOx_plane Y123;

double R = 8.31441;
//double F = 96520;
double m123_met = 554.204;
double m_O = 15.9994;
double ln10 = log(10.);
func_Tp_ptr Y2O3, BaO, CuO, O2;
double mbar2atm = log10(1./1.01325e3);
double Pa2atm = log(101325.);
double kPa2atm = log(101.325);
double R_l_atm = R/101.325;

double scaleO2 = 10.;
double scaleT = 1000.;

double T_z(double *reg)
{
  double& T = reg[0];
  double& z = reg[1];
  return T - Y123.T_tr(z);
}

double T_lnpO2(double *reg)
{
  double &T = reg[0];
  double lnpO2 = reg[1]*scaleO2;
  return T - Y123.T_tr1(lnpO2, T);
}

double X_(double *reg)
{
  double &x = reg[0];
  double T = reg[1]*scaleT;
  double &lnpO2 = reg[2];
  return x - Y123.x_eq(T, Y123.z(T, lnpO2));
}

double Z_T_lnpO2(double *reg)
{
  double &z = reg[0];
  double T = reg[1]*scaleT;
  double &lnpO2 = reg[2];
  return z - Y123.z(T, lnpO2, z);
}

double O_lnpO2_T(double *reg)
{
  double &lnpO2 = reg[0];
  double &z = reg[1];
  double &T = reg[2];
  return (lnpO2 - Y123.lnpO2(T, z));
}

static double z0_tmp;
static double m_tmp;
static double V_tmp;
static double T_tmp;

extern "C" double solve_V(double *z)
{
  return (z0_tmp - *z) -
    2.*exp(Y123.lnpO2(T_tmp, *z))*V_tmp*(m123_met + m_O*(6. + z0_tmp))
    /R_l_atm/T_tmp/m_tmp;
}

double V(double *reg)
{
  double& lnpO2 = reg[0]; // ln atm
  double& Ti = reg[1];    // 1e3/T
  z0_tmp = reg[2];
  m_tmp = reg[3];  // g
  V_tmp = reg[4];  // l
  T_tmp = 1e3/Ti;
/* old
  double z = z0 - 2.*exp(lnpO2)*V*(m123_met + m_O*(6. + z0))/R_l_atm/T/m;
  return lnpO2 - Y123.lnpO2(T, z);
*/
  double eps = 0;
  int nsig = 11, itmax = 100, ier;
  double zinit = z0_tmp - 2.*exp(Y123.lnpO2(T_tmp, z0_tmp))*
    V_tmp*(m123_met + m_O*(6. + z0_tmp))/R_l_atm/T_tmp/m_tmp;
  if (zinit < 0.) zinit = DBL_EPSILON*100.;
  double a = zinit;
  double b = z0_tmp;
  zbrent_(solve_V, &eps, &nsig, &a, &b, &itmax, &ier);
  if (ier == 130)
  {
    cout << "V - can not find z, ier = " << ier << endl;
    exit(1);
  }
  return lnpO2 - Y123.lnpO2(T_tmp, b);
}

double N(double *reg)
{
  double &lnpO2 = reg[0];
  double T = 1e3/reg[1];
  double &z = reg[2];
  return (lnpO2 - Y123.lnpO2(T, z));
}

double P_T(double *reg)
{
  double &H = reg[0];
  double &z = reg[1];
  double &T = reg[2];
  return H - Y123.HO(T, z)*R*1e-3;
}

double P_P_(double *reg)
{
  double &z1 = reg[1];
  double z2 = 0.36;
  double &H = reg[0];
  double T1 = 298.15;
  double T2 = 1057.15;
  double H1 = Y123.H(T1, z1)*R +
     0.5*Y2O3->H(T1) + 2*BaO->H(T1) + 3*CuO->H(T1) + (z1-0.5)/2.*O2->H(T1);
  double H2 = Y123.H(T2, z2)*R +
     0.5*Y2O3->H(T2) + 2*BaO->H(T2) + 3*CuO->H(T2) + (z1-0.5)/2.*O2->H(T2);
  return H - (H2 - H1)*1e-3;
}

double C(double *reg)
{
  double& Cp = reg[0];
  double T = reg[1]*scaleT;
  double& z = reg[2];
  return (Cp - (Y123.Cpz(T, z)*R +
      0.5*Y2O3->Cp(T) + 2*BaO->Cp(T) + 3*CuO->Cp(T) + (z-0.5)/2.*O2->Cp(T)));
}

double Cz(double *reg)
{
  double& Cp = reg[0];
  double& z = reg[1];
  double& T = reg[2];
  return (Cp - (Y123.Cpz(T, z)*R +
      0.5*Y2O3->Cp(T) + 2*BaO->Cp(T) + 3*CuO->Cp(T) + (z-0.5)/2.*O2->Cp(T)));
}

double S(double *reg)
{
  double& S = reg[0];
  double& z = reg[1];
  double& T = reg[2];
  return (S - (Y123.S(T, z)*R +
      0.5*Y2O3->S(T) + 2*BaO->S(T) + 3*CuO->S(T) + (z-0.5)/2.*O2->S(T)));
}

double H(double *reg)
{
  double& H = reg[0];
  double& z = reg[1];
  double& T = reg[2];
  return H - Y123.H(T, z)*R*1e-3;
}

double G_ox(double *reg)
{
  double& G = reg[0];
  double& z = reg[1];
  double& T = reg[2];
  return G - Y123.G(T, z)*R*1e-3;
}

double G_ox_O(double *reg)
{
  double& G = reg[0];
  double T = reg[1]*scaleT;
  double& lnpO2 = reg[2];
  double z = Y123.z(T, lnpO2);
  return G - Y123.G(T, z)*R*1e-3;
}

double dummy(double *reg)
{
  cout << endl << reg[0] << " " << reg[1] <<
          " this series should be redetermined" << endl;
  exit(1);
  return reg[0];
}

void fix_series()
{
//  ofstream o("tt");
  size_t i;
  double *reg;
  iser = ser.begin();
  while(iser != ser.end())
  {
    series& t = *iser;
    if (t.eq == "T_z" ||
        t.eq == "T_lnpO2" ||
        t.eq == "X_" ||
        t.eq == "Z_T_lnpO2" ||
        t.eq == "O_lnpO2_T" ||
        t.eq == "V" ||
        t.eq == "N" ||
        t.eq == "P_T" ||
        t.eq == "P_P_" ||
        t.eq == "C" ||
        t.eq == "Cz" ||
        t.eq == "S" ||
        t.eq == "H" ||
        t.eq == "G_ox" ||
        t.eq == "G_ox_O" ||
        t.eq == "Z2_T_lnpO2" ||
        t.eq == "C2" ||
        t.eq == "S2" ||
        t.eq == "H2" ||
        t.eq == "equil")
    {
    }
    else if (t.eq == "T_O")
    {
      t.eq = "T_lnpO2";
      t.f = T_lnpO2;
      for (i = 0; i < t.Ntot(); ++i)
      {
        reg = t(i);
        reg[1] = log(reg[1]);
      }
    }
    else if (t.eq == "T_O_t")
    {
      t.eq = "T_lnpO2";
      t.f = T_lnpO2;
      for (i = 0; i < t.Ntot(); ++i)
      {
        reg = t(i);
        reg[0] += 273.15;
        reg[1] = log(reg[1]);
      }
    }
    else if (t.eq == "T_O_t_kPa")
    {
      t.eq = "T_lnpO2";
      t.f = T_lnpO2;
      for (i = 0; i < t.Ntot(); ++i)
      {
        reg = t(i);
        reg[0] += 273.15;
        reg[1] = log(reg[1]) - kPa2atm;
      }
    }
    else if (t.eq == "X")
    {
      t.eq = "X_";
      t.f = X_;
      for (i = 0; i < t.Ntot(); ++i)
      {
        reg = t(i);
        reg[1] += 273.15;
        reg[2] = log(reg[2]);
      }
    }
    else if (t.eq == "Z_t_pO2")
    {
      t.eq = "Z_T_lnpO2";
      t.f = Z_T_lnpO2;
      for (i = 0; i < t.Ntot(); ++i)
      {
        reg = t(i);
        reg[1] += 273.15;
        reg[2] = log(reg[2]);
      }
    }
    else if (t.eq == "Z_t_pO2_bar")
    {
      t.eq = "Z_T_lnpO2";
      t.f = Z_T_lnpO2;
      for (i = 0; i < t.Ntot(); ++i)
      {
        reg = t(i);
        reg[1] += 273.15;
        reg[2] = log(reg[2]*1.01325);
      }
    }
    else if (t.eq == "Z_t_lgpO2")
    {
      t.eq = "Z_T_lnpO2";
      t.f = Z_T_lnpO2;
      for (i = 0; i < t.Ntot(); ++i)
      {
        reg = t(i);
        reg[1] += 273.15;
        reg[2] *= ln10;
      }
    }
    else if (t.eq == "Z_t_lnpO2")
    {
      t.eq = "Z_T_lnpO2";
      t.f = Z_T_lnpO2;
      for (i = 0; i < t.Ntot(); ++i)
      {
        reg = t(i);
        reg[1] += 273.15;
      }
    }
    else if (t.eq == "Z_T_lgpO2")
    {
      t.eq = "Z_T_lnpO2";
      t.f = Z_T_lnpO2;
      for (i = 0; i < t.Ntot(); ++i)
      {
        reg = t(i);
        reg[2] *= ln10;
      }
    }
    else if (t.eq == "Z_Ti_pO2_Pa")
    {
      t.eq = "Z_T_lnpO2";
      t.f = Z_T_lnpO2;
      for (i = 0; i < t.Ntot(); ++i)
      {
        reg = t(i);
        reg[1] = 1000./reg[1];
        reg[2] = log(reg[2]) - Pa2atm;
      }
    }
    else if (t.eq == "Z_T_pO2_Pa")
    {
      t.eq = "Z_T_lnpO2";
      t.f = Z_T_lnpO2;
      for (i = 0; i < t.Ntot(); ++i)
      {
        reg = t(i);
        reg[2] = log(reg[2]) - Pa2atm;
      }
    }
    else if (t.eq == "O_lgpO2_mbar_t")
    {
      t.eq = "O_lnpO2_T";
      t.f = O_lnpO2_T;
      for (i = 0; i < t.Ntot(); ++i)
      {
        reg = t(i);
        reg[0] = (reg[0] + mbar2atm)*ln10;
        reg[2] += 273.15;
      }
    }
    else if (t.eq == "O_lgpO2_t")
    {
      t.eq = "O_lnpO2_T";
      t.f = O_lnpO2_T;
      for (i = 0; i < t.Ntot(); ++i)
      {
        reg = t(i);
        reg[0] *= ln10;
        reg[2] += 273.15;
      }
    }
    else if (t.eq == "O_lgpO2_bar_T")
    {
      t.eq = "O_lnpO2_T";
      t.f = O_lnpO2_T;
      for (i = 0; i < t.Ntot(); ++i)
      {
        reg = t(i);
        reg[0] = reg[0]*ln10 - log10(1.01325);
      }
    }
    else if (t.eq == "O_lgpO2_T_zi")
    {
      t.eq = "O_lnpO2_T";
      t.f = O_lnpO2_T;
      for (i = 0; i < t.Ntot(); ++i)
      {
        reg = t(i);
        reg[0] *= ln10;
        reg[1] = 1. - reg[1];
      }
    }
    else if (t.eq == "O_lnpO2_t")
    {
      t.eq = "O_lnpO2_T";
      t.f = O_lnpO2_T;
      for (i = 0; i < t.Ntot(); ++i)
      {
        reg = t(i);
        reg[2] += 273.15;
      }
    }
    else if (t.eq == "O_pO2_t")
    {
      t.eq = "O_lnpO2_T";
      t.f = O_lnpO2_T;
      for (i = 0; i < t.Ntot(); ++i)
      {
        reg = t(i);
        reg[0] = log(reg[0]);
        reg[2] += 273.15;
      }
    }
    else if (t.eq == "O_pO2_torr_t")
    {
      t.eq = "O_lnpO2_T";
      t.f = O_lnpO2_T;
      for (i = 0; i < t.Ntot(); ++i)
      {
        reg = t(i);
        reg[0] = log(reg[0]/760.);
        reg[2] += 273.15;
      }
    }
    else if (t.eq == "V_tor_t_cm")
    {
      t.eq = "V";
      t.f = V;
      for (i = 0; i < t.Ntot(); ++i)
      {
        reg = t(i);
        reg[0] = log(reg[0]/760.);
        reg[1] = 1e3/(reg[1] + 273.15);
        reg[4] *= 1e-3;
      }
    }
    else if (t.eq == "P_t")
    {
      t.eq = "P_T";
      t.f = P_T;
      for (i = 0; i < t.Ntot(); ++i)
      {
        reg = t(i);
        reg[0] *= 1e-3;
        reg[2] += 273.15;
      }
    }
    else if (t.eq == "P_P")
    {
      t.eq = "P_P_";
      t.f = P_P_;
      for (i = 0; i < t.Ntot(); ++i)
      {
        reg = t(i);
        reg[0] *= (m123_met + m_O*(6. + reg[1]))*1e-3;
      }
    }
    else
    {
      cout << "Not known series: " << t.eq << endl;
      exit(1);
    }
    t.set_av();
//    o << t << endl << endl;
    iser++;
  }
}

void scale_x()
{
  size_t i;
  double *reg;
  iser = ser.begin();
  while(iser != ser.end())
  {
    series& t = *iser;
    if (t.eq == "T_z" ||
        t.eq == "O_lnpO2_T" ||
        t.eq == "P_T" ||
        t.eq == "P_P_" ||
        t.eq == "Cz" ||
        t.eq == "S" ||
        t.eq == "H" ||
        t.eq == "S2" ||
        t.eq == "H2" ||
        t.eq == "V" ||   //??
        t.eq == "N" ||   //??
        t.eq == "G_ox")  //??
    {
    }
    else if (t.eq == "T_lnpO2" ||
        t.eq == "equil")
    {
      for (i = 0; i < t.Ntot(); ++i)
      {
        reg = t(i);
        reg[1] /= scaleO2;
      }
    }
    else if (t.eq == "X_" ||
        t.eq == "Z_T_lnpO2" ||
        t.eq == "C" ||
        t.eq == "G_ox_O" ||
        t.eq == "Z2_T_lnpO2" ||
        t.eq == "C2")
    {
      for (i = 0; i < t.Ntot(); ++i)
      {
        reg = t(i);
        reg[1] /= scaleT;
      }
    }
    else
    {
      cout << "scale_x: Not known series: " << t.eq << endl;
      exit(1);
    }
    t.set_av();
    iser++;
  }
}

double scale(const cString &eq)
{
  if (eq == "T_z" ||
      eq == "O_lnpO2_T" ||
      eq == "P_T" ||
      eq == "P_P_" ||
      eq == "Cz" ||
      eq == "S" ||
      eq == "H" ||
      eq == "S2" ||
      eq == "H2" ||
      eq == "V" ||   //??
      eq == "N" ||   //??
      eq == "G_ox")  //??
  {
    return 1.;
  }
  else if (eq == "T_lnpO2" ||
      eq == "equil")
  {
    return scaleO2;
  }
  else if (eq == "X_" ||
      eq == "Z_T_lnpO2" ||
      eq == "C" ||
      eq == "G_ox_O" ||
      eq == "Z2_T_lnpO2" ||
      eq == "C2")
  {
    return scaleT;
  }
  return 1.;
}

void init_res()
{
  res.insert(residual("T_z", T_z, 2, 0xFFFFFFFF));
  res.insert(residual("T_lnpO2", T_lnpO2, 2, 0xFFFFFFFF));
  res.insert(residual("X_", X_, 3, 0xFFFFFFFF));
  res.insert(residual("Z_T_lnpO2", Z_T_lnpO2, 3, 0xFFFFFFFF));
  res.insert(residual("O_lnpO2_T", O_lnpO2_T, 3, 0xFFFFFFFF));
  res.insert(residual("V", V, 5, 0xFFFFFFFF));
  res.insert(residual("N", N, 3, 0xFFFFFFFF));
  res.insert(residual("P_T", P_T, 3, 0xFFFFFFFF));
  res.insert(residual("P_P_", P_P_, 2, 0xFFFFFFFF));
  res.insert(residual("C", C, 3, 0xFFFFFFFF));
  res.insert(residual("Cz", Cz, 3, 0xFFFFFFFF));
  res.insert(residual("S", S, 3, 0xFFFFFFFF));
  res.insert(residual("H", H, 3, 0xFFFFFFFF));
  res.insert(residual("G_ox", G_ox, 3, 0xFFFFFFFF));
  res.insert(residual("G_ox_O", G_ox_O, 3, 0xFFFFFFFF));
/*
  res.insert(residual("Z2_T_lnpO2", Z2_T_lnpO2, 3, 0xFFFFFFFF));
  res.insert(residual("C2", C2, 3, 0xFFFFFFFF));
  res.insert(residual("S2", S2, 3, 0xFFFFFFFF));
  res.insert(residual("H2", H2, 3, 0xFFFFFFFF));
  res.insert(residual("equil", equil, 2, 0xFFFFFFFF));
  res.insert(residual("equil_", equil_, 2, 0xFFFFFFFF));
*/

  res.insert(residual("T_O", dummy, 2, 0xFFFFFFFF));
  res.insert(residual("T_O_t", dummy, 2, 0xFFFFFFFF));
  res.insert(residual("T_O_t_kPa", dummy, 2, 0xFFFFFFFF));
  res.insert(residual("X", dummy, 3, 0xFFFFFFFF));
  res.insert(residual("Z_t_pO2", dummy, 3, 0xFFFFFFFF));
  res.insert(residual("Z_t_pO2_bar", dummy, 3, 0xFFFFFFFF));
  res.insert(residual("Z_t_lgpO2", dummy, 3, 0xFFFFFFFF));
  res.insert(residual("Z_t_lnpO2", dummy, 3, 0xFFFFFFFF));
  res.insert(residual("Z_T_lgpO2", dummy, 3, 0xFFFFFFFF));
  res.insert(residual("Z_Ti_pO2_Pa", dummy, 3, 0xFFFFFFFF));
  res.insert(residual("Z_T_pO2_Pa", dummy, 3, 0xFFFFFFFF));
  res.insert(residual("O_lgpO2_mbar_t", dummy, 3, 0xFFFFFFFF));
  res.insert(residual("O_lgpO2_t", dummy, 3, 0xFFFFFFFF));
  res.insert(residual("O_lgpO2_bar_T", dummy, 3, 0xFFFFFFFF));
  res.insert(residual("O_lgpO2_T_zi", dummy, 3, 0xFFFFFFFF));
  res.insert(residual("O_lnpO2_t", dummy, 3, 0xFFFFFFFF));
  res.insert(residual("O_pO2_t", dummy, 3, 0xFFFFFFFF));
  res.insert(residual("O_pO2_torr_t", dummy, 3, 0xFFFFFFFF));
  res.insert(residual("V_tor_t_cm", dummy, 5, 0xFFFFFFFF));
  res.insert(residual("P_t", dummy, 3, 0xFFFFFFFF));
  res.insert(residual("P_P", dummy, 2, 0xFFFFFFFF));
  res.insert(residual("Z2_t_pO2", dummy, 3, 0xFFFFFFFF));
  res.insert(residual("Z2_t_pO2_bar", dummy, 3, 0xFFFFFFFF));
  res.insert(residual("equil_t_atm", dummy, 2, 0xFFFFFFFF));
  res.insert(residual("equil_lgp", dummy, 2, 0xFFFFFFFF));
  res.insert(residual("equil_t_lgp", dummy, 2, 0xFFFFFFFF));
  res.insert(residual("equil_t_atm_", dummy, 2, 0xFFFFFFFF));
  res.insert(residual("equil_lgp_", dummy, 2, 0xFFFFFFFF));
  res.insert(residual("equil_t_lgp_", dummy, 2, 0xFFFFFFFF));
}

char* EQ = "";

void init_par()
{
  fix_series();
  scale_x();
  ifstream in("y123.mod");
  in >> Y123 >> Y2O3 >> BaO >> CuO >> O2;
  CuOx_plane::iterator i;
  nparam = 0;
  for (i = Y123.begin(); !(i == Y123.end()); ++i)
  {
    par[nparam].x = (*i).x;
    par[nparam].name = (*i).name.c_str();
    par[nparam].atr = (boolean)(*i).unknown;
    ++nparam;
  }

}

void init_est()
{
  CuOx_plane::iterator i;
  int n = 0;
  for (i = Y123.begin(); !(i == Y123.end()); ++i)
    (*i).x = par[n++].x;
}

void test_convex_Y123()
{
  double p = 1.;
  of << "testing Gibbs energy of the Y123" << endl;
  double old, next;
  for (double T = 200.; T <= 1500.01; T+=25.)
  {
    int okay = 1;
    old = -HUGE_VAL;
    for (double z = 0.01; z <= 0.99; z+= 0.005)
    {
      next = Y123.GO(T, z);
      if ( next <= old)
      {
        if (okay)
        {
          of << setw(wdth) << T
             << setw(wdth) << z;
          okay = 0;
        }
        else if (z > 0.98)
          of << setw(wdth) << z << endl;
      }
      else
      {
        if (!okay)
        {
          of << setw(wdth) << z << endl;
          okay = 1;
        }
      }
      old = next;
    }
  }
  of << endl;
}

void post_analysis()
{
  ofstream out((file_name + "mod").c_str());
  out << Y123 << Y2O3 << BaO << CuO << O2;

  test_convex_Y123();
  of << endl;
}