/*
Copyright (C) 1995 Evgenii Rudnyi, http://Evgenii.Rudnyi.Ru/

Software for paper 
E. B. Rudnyi. Statistical model of systematic errors: An assessment of the Ba-Cu and Cu-Y phase diagram. Chemometrics and Intelligent Laboratory Systems, 1997, v. 36, p. 213-227.

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

#include <math.h>
#include <float.h>
#include <iomanip.h>
#include <sumsqr.h>
#include "bin_td.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 = 18;
param par[18];

double& A0 = par[4].x;
double& A1 = par[5].x;
double& A2 = par[6].x;
double& A3 = par[7].x;
double& B0 = par[8].x;
double& B1 = par[9].x;
double& B2 = par[10].x;
double& B3 = par[11].x;
double& C0 = par[14].x;
double& C1 = par[15].x;
double& C2 = par[16].x;
double& C3 = par[17].x;

double h1s(double& T)
{
  cerr << "function h1s is not defined" << endl;
  exit(1);
  return T*0.;
}

double h2s(double& T)
{
  cerr << "function h2s is not defined" << endl;
  exit(1);
  return T*0.;
}

double g1s(double& T)
{
// Ba - 91Din (based on TPIS)
  return (T < 1000.) ?
       -17685.226+233.78606*T-42.889*T*log(T)-1.8314e-3*pow(T,2)
       -0.000095e-6*pow(T,3)+705880/T:    //  298.15,1000
       -64873.614+608.188389*T-94.2824199*T*log(T)+19.504772e-3*pow(T,2)
       -1.051353e-6*pow(T,3)+8220192/T;   //  1000,2995
}

double g1l(double& T, double x)
{
// Ba - 91Din (based on TPIS)
  double g;
  if (T < 1000.)
    g = -9738.988+229.540143*T-43.4961089*T*log(T)-2.346416e-3*pow(T,2)
        +0.991223e-6*pow(T,3)+723016/T;   //  298.15,1000
  else
    g = -7381.093+235.49642*T-45.103*T*log(T)+2.154e-3*pow(T,2)
        +0.000027e-6*pow(T,3)-365/T;      //  1000,2995
  return (x) ? g + Dmixg1(T, x) : g;
}

double g2s(double& T)
{
// Cu - 91Din (based on JANAF)
  return (T < 1357.77) ?
    -7770.458+130.485235*T-24.112392*T*log(T)-2.65684e-3*pow(T,2)
    +0.129223e-6*pow(T,3)+52478/T : // 298.15,1357.77
    -13542.026+183.803828*T-31.38*T*log(T)+3.642e29*pow(T,-9);
                                   // 1357.77,3200
}

double g2l(double& T, double x)
{
// Cu - 91Din (based on JANAF)
  double g;
  if (T < 1357.77)
    g = 5194.277+120.973331*T-24.112392*T*log(T)-2.65684e-3*pow(T,2)
        +0.129223e-6*pow(T,3)+52478/T-5.849e-21*pow(T,7);
                                    // 298.15,1357.77
  else
    g = -46.545+173.881484*T-31.38*T*log(T);
                                    // 1357.77,3200
  return (x) ? g + Dmixg2(T, x) : g;
}

elem e1("Ba", 0., 1000.);
elem e2("Cu", 1., 1358.);

int nlp = 2;
double tttt;
linephase lp[5] =
  {linephase("BaCu", 0.5, par[0].x, par[1].x, par[12].x),
   linephase("BaCu13", 13./14., par[2].x, par[3].x, par[13].x),
   linephase("", 0., tttt, tttt, tttt),
   linephase("", 0., tttt, tttt, tttt),
   linephase("", 0., tttt, tttt, tttt)
  };

int nnon = 3;
nonvar non[6] = {{731., 0.239}, {823., 0.50}, {1031., 0.704}};

double e1_L(double *reg)
{
  subst1 = &e1;
  return T_solve(reg);
}

double lp0_L(double *reg)
{
  subst1 = &lp[0];
  return T_solve(reg);
}

double lp1_L(double *reg)
{
  subst1 = &lp[1];
  return T_solve(reg);
}

double e2_L(double *reg)
{
  subst1 = &e2;
  return T_solve(reg);
}

double e1_L_lp0(double* reg)
{
  subst1 = &e1;
  subst2 = &lp[0];
  return Tx_solve(reg, 0);
}

double lp0_L_lp1(double* reg)
{
  subst1 = &lp[0];
  subst2 = &lp[1];
  return Tx_solve(reg, 1);
}

double lp1_L_e2(double* reg)
{
  subst1 = &lp[1];
  subst2 = &e2;
  return Tx_solve(reg, 2);
}

double Hlp0(double* reg)
{
  double& H = reg[0];
  return H - (lp[0].H + lp[0].C*298.15);
//  return H - lp[0].H;
}

double h1(double* reg)
{
  double& H = reg[0];
  double& x = reg[1];
  double& T = reg[2];
  return H - Dmixh1(T, x);
}

double h2(double* reg)
{
  double& H = reg[0];
  double& x = reg[1];
  double& T = reg[2];
  return H - Dmixh2(T, x);
}

void init_res()
{
// C3 C2  C1 C0 Cb Ca  B3 B2 B1 B0  A3 A2 A1 A0  Bb Ab Ba Aa
// 18 17  16 15 14 13  12 11 10  9   8  7  6  5   4  3  2  1
  res.insert(residual("Ba_L", e1_L, 2, 0x3CFF0));
                                          // A0 - B3 C0-C3
  res.insert(residual("BaCu_L", lp0_L, 2, 0x3DFF3));
                                          // Aa Ba A0 - B3 Ca C0-C3
  res.insert(residual("BaCu13_L", lp1_L, 2, 0x3EFFC));
                                          // Ab Bb A0 - B3 Cb C0-C3
  res.insert(residual("Cu_L", e2_L, 2, 0x3CFF0));
                                          // A0 - B3 C0-C3
  res.insert(residual("Ba_L_BaCu", e1_L_lp0, 1, 0x3DFF3));
                                          // Aa Ba A0 - B3 Ca C0-C3
  res.insert(residual("BaCu_L_BaCu13", lp0_L_lp1, 1, 0x3FFFF));
                                          // Aa Ba Ab Bb A0 - B3 Ca Cb C0-C3
  res.insert(residual("BaCu13_L_Cu", lp1_L_e2, 1, 0x3EFFC));
                                          // Ab Bb A0 - B3 Cb C0-C3
  res.insert(residual("HBaCu", Hlp0, 1, 0x1001)); // Aa Ca
  res.insert(residual("HBa", h1, 3, 0x3C0F0));     // A0-A3 C0-C3
  res.insert(residual("HCu", h2, 3, 0x3C0F0));     // A0-A3 C0-C3
}

char* EQ = "\
There are equations as follows \n\
\t\t\t    eq_ID                    equation \n\n\
\t\t\tBa_L                T = T(x, teta_vec) + error\n\
\t\t\tBaCu_L              T = T(x, teta_vec) + error\n\
\t\t\tBaCu13_L            T = T(x, teta_vec) + error\n\
\t\t\tCu_L                T = T(x, teta_vec) + error\n\
\t\t\tBa_L_BaCu           T = T(x, teta_vec) + error\n\
\t\t\tBaCu_L_BaCu13       T = T(x, teta_vec) + error\n\
\t\t\tBaCu13_L_Cu         T = T(x, teta_vec) + error\n\
\t\t\tHBaCu               H = Aa + error\n\
\t\t\tHBa                 HBa = HBa(x, T, teta_vec) + error\n\
\t\t\tHCu                 HCu = HCu(x, T, teta_vec) + error\n\
availale in this release \n\
and twelve unknown parameters \n\
\t\tDelH     DelS \n\
\t\t Aa       Bb     Ca    0.5 Ba + 0.5 Cu = 0.5 BaCu \n\
\t\t Aa       Bb     Cb    1/14 Ba + 13/14 Cu = 1/14 BaCu13 \n\
\t\t A0       B0     C0    x(1-x)          terms for BaCu(l) \n\
\t\t A1       B1     C1    x(1-x)(2x-1)  \n\
\t\t A2       B2     C2    x(1-x)(2x-1)^2  \n\
\t\t A3       B3     C3    x(1-x)(2x-1)^3  \n";

void init_par()
{
  par[0].name = "Aa";
  par[0].atr = true;
  par[1].name = "Ba";
  par[1].atr = false;
  par[2].name = "Ab";
  par[2].atr = true;
  par[3].name = "Bb";
  par[3].atr = false;
  par[4].name = "A0";
  par[4].atr = false;
  par[5].name = "A1";
  par[5].atr = false;
  par[6].name = "A2";
  par[6].atr = false;
  par[7].name = "A3";
  par[7].atr = false;
  par[8].name = "B0";
  par[8].atr = false;
  par[9].name = "B1";
  par[9].atr = false;
  par[10].name = "B2";
  par[10].atr = false;
  par[11].name = "B3";
  par[11].atr = false;
  par[12].name = "Ca";
  par[12].atr = false;
  par[13].name = "Cb";
  par[13].atr = false;
  par[14].name = "C0";
  par[14].atr = false;
  par[15].name = "C1";
  par[15].atr = false;
  par[16].name = "C2";
  par[16].atr = false;
  par[17].name = "C3";
  par[17].atr = false;
  log_name = file_name + "log";
  if (strlen(file_name.getstr()) < 9)
  {
    file_name[strlen(file_name.getstr()) - 1] = '0';
    file_name += ".";
  }
}

void print_non()
{
  double T, x;
  switch (iii)
  {
    case 1:
      T = non[0].T;
      x = 0.;
      break;
    case 2:
      T = non[0].T;
      x = lp[0].xs;
      break;
    case 3:
      T = 600.;
      x = lp[0].xs;
      break;
    case 4:
      T = lp[0].Tm;
      x = lp[0].xs;
      break;
    case 5:
      T = non[1].T;
      x = lp[0].xs;
      break;
    case 6:
      T = non[1].T;
      x = lp[1].xs;
      break;
    case 7:
      T = 600.;
      x = lp[1].xs;
      break;
    case 8:
      T = non[2].T;
      x = lp[1].xs;
      break;
    case 9:
      T = non[2].T;
      x = non[2].x;
      break;
    case 10:
      T = non[2].T;
      x = 1.;
      break;
    default:
      ofpd << endl;
      return;
  }
  iii++;
  ofpd << setw(wdth) << T
     << setw(wdth) << x  << endl;
}

void print_Tx(cString& txt, double&T, double& x)
{
  double T1 = 1273.;
  double T2 = 2000.;
  double h;
  ofpd << setw(wdth) << txt
     << setw(wdth) << T
     << setw(wdth) << x;
  h = Dmixh1(T1, x);
  if (x == 0. || x >= .599)
    ofpd << setw(wdth) << h;
  else
    ofpd << setw(wdth) << "miss";
  h = Dmixh2(T1, x);
  if (x < .599)
    ofpd << setw(wdth) << h;
  else
    ofpd << setw(wdth) << "miss";
  ofpd << setw(wdth) << DmixG(T1, x)
     << setw(wdth) << DmixG(T2, x);
  print_non();
}

void post_analysis1()
{
  ofpd.open((file_name + "axm").c_str());
  cString txt;
  double reg[2];
  double& T = reg[0];
  double& x = reg[1];
  T = 1000;
  if (!non[0].est) e1_L_lp0(reg);
  if (!non[1].est) lp0_L_lp1(reg);
  if (!non[2].est) lp1_L_e2(reg);
  iii = 1;
  T = 830.;
  x = 0.5;
  subst1 = &lp[0];
  lp[0].Tm = T - T_solve(reg);
  ofpd << "equilibria        T          x          H"
     << e1.name << "          H"
     << e2.name << "      DmixG1       DmixG2       nonT       nonx" << endl;
  T = e1.Tm;
  x = 0.;
  txt = cString(e1.name) + cString("_L");
  print_Tx(txt, T, x);
  for (int i = -1; i <= nlp; i++)
  {
    T = (x) ? non[i].T : e1.Tm;
    for (x = (x) ? ((int(non[i].x*100.)/2 + 1)*2)/100. : 0.02;
         (i < nlp) ? x < non[i + 1].x : x < 0.9801;
         x+= 0.02)
    {
      if (i == -1)
        subst1 = &e1;
      else if (i == nlp)
        subst1 = &e2;
      else
        subst1 = &lp[i];
      T -= T_solve(reg);
      txt = cString(subst1->name) + cString("_L");
      print_Tx(txt, T, x);
    }
    if (i < nlp)
    {
      if (i < nlp - 1)
        txt += cString("_") + cString(lp[i+1].name);
      else
        txt += cString("_") + cString(e2.name);
      print_Tx(txt, non[i+1].T, non[i+1].x);
    }
  }
  T = e2.Tm;
  x = 1.;
  txt = cString(e2.name) + cString("_L");
  print_Tx(txt, T, x);
  ofpd.close();
}

//  if ((par[0].use || par[1].use) && (par[2].use || par[3].use))

void post_analysis()
{
  int i;
  file_name[strlen(file_name.getstr()) - 2] = '0' + irun++;
  of << "******************** post analysis *************************"
     << endl;
  for (i = 0; i < nnon; ++i)
    if (!non[i].conv)
      of << "WARNING - no solution, non_point number "
         << i
         << setw(wdth) << non[i].T1
         << setw(wdth) << non[i].T2  << endl;
  test_convex_L();
  test_convex_solids();
  ofstream of1(log_name.getstr(), ios::out | ios::app);
  of1.precision(prcsn);
  of1.setf(ios::showpoint);
  double L, SS;
  eval_L(L, SS);
  of1 << setw(wdth) << L << setw(wdth) << SS;
  for (i = 1; i <= 3; ++i)
    of1 << setw(wdth) << sqrt(sr2t[i])
        << setw(wdth) << sqrt(gat[i])
        << setw(wdth) << sqrt(gbt[i]);
  for (i = 0; i < nparam; ++i)
    of1 << setw(wdth) << par[i].x;
  of1 << endl;
  post_analysis1();
  set_hide();
  read_extra_x("bacu_x.dat");
  read_errors("bacu_e.dat");
  deviates();
  eb_vs_ea();
}