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

Software for paper 
E. B. Rudnyi, G. F. Voronin. Thermodynamic assessment of the ternary Ba-Cu-Y system. Direct optimization of the miscibility gap. CALPHAD, 1996, v. 20, N 3, p. 297-305. 

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 <imsl.h>
#include <sumsqr.h>

double dGdxCu(double& T, double& XCU, double& XY);
double dGdxY(double& T, double& XCU, double& XY);
double d2GdxCu2(double& T, double& XCU, double& XY);
double d2GdxY2(double& T, double& XCU, double& XY);
double d2GdxCudxY(double& T, double& XCU, double& XY);
double d3GdxCu3(double& T, double& XCU, double& XY);
double d3GdxCu2dxY(double& T, double& XCU, double& XY);
double d3GdxY3(double& T, double& XCU, double& XY);
double d3GdxCudxY2(double& T, double& XCU, double& XY);
double DmixG(double& T, double& XCU, double& XY);
double Dmixg1(double& T, double& XCU, double& XY);
double Dmixg2(double& T, double& XCU, double& XY);
double Dmixg3(double& T, double& XCU, double& XY);

char* LICENSE = "\
Evgenii Rudnyi 1994, (c) All rights reserved \n\n\
Chemistry Department \n\
Moscow State University \n\
119899 Moscow, Russia \n\
tel. (095)939-5452 \n\
rudnyi@mch.chem.msu.su \n\n\
This program is freeware (public domain). \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 = 3;
param par[3];

double& A0BAY = par[0].x;
double& A1BAY = par[1].x;
double& A0BACUY = par[2].x;

double Tsys = 1273.;
// /*   my work
double A0BACU = -7191 + Tsys*3.363;
double A1BACU = 0.;
double A2BACU = 0.;
double A3BACU = 0.;
double A0  = -8.246e+04;
double B0  =      14.74;
double C0  =        0.0;
double A1  =  3.848e+04;
double B1  =     -6.607;
double C1  =        0.0;
double A0CUY = A0 + Tsys*B0 + C0*Tsys*log(Tsys);
double A1CUY = A1 + Tsys*B1 + C1*Tsys*log(Tsys);
double A2CUY = 0.;
double A3CUY = 0.;
// */

/* // 90Mey
double A0BACU = -764.8 - Tsys*18.96614;
double A1BACU = -2595.8 + Tsys*16.93798;
double A2BACU = -2879.6;
double A3BACU = -1970.9;
double A0  = -124182.5;
double A1  = +116962.9;
double A2  =  -80027.4;
double B0  =  311.58059;
double B1  = -432.28393;
double B2  =  518.93980;
double C0  = -33.929066;
double C1  =  48.942898;
double C2  = -63.802135;
double A0CUY = A0 + Tsys*B0 + C0*Tsys*log(Tsys);
double A1CUY = A1 + Tsys*B1 + C1*Tsys*log(Tsys);
double A2CUY = A2 + Tsys*B2 + C2*Tsys*log(Tsys);
double A3CUY = 0.;
*/

double GCu = -824.39;
double GY = -4110.3;

double R = 8.31441;

double partmp[5];

int x_switch = 1;
double xCu1calc;
double xY1calc;
double xCu2calc;
double xY2calc;
double xCuctr;
double xYctr;
double Ttmp;

double xCu_cryt;
double xY_cryt;
boolean cryt;
boolean init_cryt;

double mu1(double* reg)
{
  double& T = reg[0];
  double& xCu1 = reg[1];
  double& xY1  = reg[2];
  double& xCu2 = reg[3];
  double& xY2  = reg[4];
  return Dmixg1(T, xCu1, xY1) - Dmixg1(T, xCu2, xY2);
}

double mu2(double* reg)
{
  double& T = reg[0];
  double& xCu1 = reg[1];
  double& xY1  = reg[2];
  double& xCu2 = reg[3];
  double& xY2  = reg[4];
  return Dmixg2(T, xCu1, xY1) - Dmixg2(T, xCu2, xY2);
}

double mu3(double* reg)
{
  double& T = reg[0];
  double& xCu1 = reg[1];
  double& xY1  = reg[2];
  double& xCu2 = reg[3];
  double& xY2  = reg[4];
  return Dmixg3(T, xCu1, xY1) - Dmixg3(T, xCu2, xY2);
}

char* EQ = "\
There are equations as follows \n\
\teq_ID              equation \n\n\
\tmu1   mu1(T, xCu1, xY1, teta_vec) = mu1(T, xCu2, xY2) + error\n\
\tmu2   mu2(T, xCu1, xY1, teta_vec) = mu2(T, xCu2, xY2) + error\n\
\tmu3   mu3(T, xCu1, xY1, teta_vec) = mu3(T, xCu2, xY2) + error\n\
\tdelx  xi = xi(teta_vec) + error\n\
availale in this release \n\
and three unknown parameters \n\
\t\t A0BaY    term for xBa*xY \n\
\t\t A1BaY    term for xBa*xY*(xY-xBa)  \n\
\t\t A0BaCuY  term for xBa*xY*xCu \n";

void init_par()
{
  par[0].name = "A0BaY";
  par[0].atr = true;
  par[1].name = "A1BaY";
  par[1].atr = false;
  par[2].name = "A0BaCuY";
  par[2].atr = false;
  init_cryt = false;
}

void init_est()
{
  cryt = false;
}

void xtoz1(double x, double& z, double scale)
{
  z = log(x/(scale - x));
}

void ztox1(double& x, double z, double scale)
{
  x = scale/(1 + exp(-z));
}

//     0 < x1 < 1,      0 < x2 < 1,     x1 + x2 < 1
//  -inf < z1 < inf, -inf < z1 < inf

void xtoz(double x1, double x2, double& z1, double &z2)
{
  double a;
  if (x2 >= x1)
    a = x1/x2;
  else
    a = x2/x1;
  double x11 = x1*(1. + a);
  double x21 = x2*(1. + a);
  z1 = log(x11/(1. - x11));
  z2 = log(x21/(1. - x21));
}

void ztox(double& x1, double& x2, double z1, double z2)
{
  double x11 = 1./(1 + exp(-z1));
  double x21 = 1./(1 + exp(-z2));
  double a;
  if (x21 >= x11)
    a = x11/x21;
  else
    a = x21/x11;
  x1 = x11/(1. + a);
  x2 = x21/(1. + a);
}

extern "C" {
void difcryt(double* x, int *m, int *n, double* f)
{
  double xCu, xY;
  ztox(xCu, xY, x[0], x[1]);
  f[0] = d2GdxCu2(Ttmp, xCu, xY)/d2GdxCudxY(Ttmp, xCu, xY)
         *d2GdxY2(Ttmp, xCu, xY)/d2GdxCudxY(Ttmp, xCu, xY)
         - 1.;
/*
  f[0] = (d2GdxCu2(Ttmp, xCu, xY)
         *d2GdxY2(Ttmp, xCu, xY)
         - pow(d2GdxCudxY(Ttmp, xCu, xY), 2));
*/
  double tmp1 = d3GdxCu3(Ttmp, xCu, xY)
                *d2GdxY2(Ttmp, xCu, xY)
                + d2GdxCu2(Ttmp, xCu, xY)
                *d3GdxCudxY2(Ttmp, xCu, xY)
                - 2.*d2GdxCudxY(Ttmp, xCu, xY)
                *d3GdxCu2dxY(Ttmp, xCu, xY);
  double tmp2 = d3GdxCu2dxY(Ttmp, xCu, xY)
                *d2GdxY2(Ttmp, xCu, xY)
                + d2GdxCu2(Ttmp, xCu, xY)
                *d3GdxY3(Ttmp, xCu, xY)
                - 2.*d2GdxCudxY(Ttmp, xCu, xY)
                *d3GdxCudxY2(Ttmp, xCu, xY);
  f[1] = tmp1/tmp2*d2GdxY2(Ttmp, xCu, xY)/d2GdxCudxY(Ttmp, xCu, xY)
         - 1.;
/*
  f[1] = (tmp1*d2GdxY2(Ttmp, xCu, xY) -
         tmp2*d2GdxCudxY(Ttmp, xCu, xY));
*/
  if (print == all_iter)
    of << "cryt "
       << setw(wdth) << xCu
       << setw(wdth) << xY
       << setw(wdth) << f[0]
       << setw(wdth) << f[1]
       << setw(wdth) << sqrt(pow(f[0],2) + pow(f[1],2))
       << endl;
}
}

void cryt_solve(void)
{
  int n = 2;
  int m = n;
  double x[2];
  double ff[2];
  double xCu;
  double xY;
  double min = HUGE_VAL;
  double xCumin;
  double xYmin;
  double tmp;
  if (!init_cryt)
  {
    for (xCu = 0.65; xCu < 0.901 ; xCu += 0.02)
      for (xY = 0.01; xY < 1. - xCu ; xY += 0.02)
      {
        xtoz(xCu, xY, x[0], x[1]);
        difcryt(x, &m, &n, ff);
        if ((tmp = pow(ff[0],2) + pow(ff[1],2)) < min)
        {
          min = tmp;
          xCumin = xCu;
          xYmin = xY;
        }
      }
    init_cryt = true;
    xtoz(xCumin, xYmin, x[0], x[1]);
  }
  else
    xtoz(xCu_cryt, xY_cryt, x[0], x[1]);
  int nsig = 10;
  double eps = 0;
  double delta = 0;
  int maxfn = 300;
  int iopt = 1;
  int ixjac;
  int ier;
  int infer;
  double parm[4];
  double ssq;
  ixjac = m;
  double* f = new double[m];
  double* xjac = new double[n*m];
  double* xjtj = new double[n*(n+1)/2];
  double* work = new double[5*n + 2*m + (n+1)*n/2];
  double zxssq_machine_eps = 12;
  zxssq1_(difcryt, &m, &n, &nsig, &eps, &delta, &maxfn, &iopt,
        parm, x, &ssq, f, xjac, &ixjac, xjtj, work, &infer, &ier, &zxssq_machine_eps);
  delete f;
  delete xjac;
  delete xjtj;
  delete work;
  ztox(xCu_cryt, xY_cryt, x[0], x[1]);
  if (sqrt(ssq) > 1e-10)
  {
    of << "critical point - warning" << endl;
    print_zxssq_1(infer, ier, work);
    of << setw(wdth) << xCu_cryt << setw(wdth) << xY_cryt;
    of << setw(wdth) << sqrt(ssq) << endl << endl;
  }
}

extern "C" {
void difternary(double* x, int *m, int *n, double* f)
{
  ztox(partmp[3], partmp[4], x[1], x[2]);
/*
  // y = (y2-y1)/(x2-x1)*(x-x1) + y1
  partmp[1] = (partmp[3] - xCuctr)
             /(partmp[4] - xYctr)
             *partmp[2] + xCuctr;
*/
  double a = (partmp[3] - xCuctr)
            /(partmp[4] - xYctr);
  double b = xCuctr - a*xYctr;
  double d = 1./(1. + exp(-x[0]));
  partmp[2] = (d - b)/(1. + a);
  partmp[1] = a*partmp[2] + b;
  if (partmp[1] < 0.1e-8) partmp[1] = 0.1e-8;
  if (partmp[1] > partmp[2])
  {
    partmp[2] = (-(b - d*a) + sqrt(pow((b - d*a),2) + 4.*(a + 1.)*d*b))
               /2./(a + 1.);
    partmp[1] = a*partmp[2] + b;
    if (partmp[1] < partmp[2])
    {
      of << "something is wrong in difternary - one" << endl;
      exit(1);
    }
  }
  if (partmp[1] + partmp[2] >= 1.)
  {
    partmp[2] = (0.999999 - b)/(a + 1.);
    partmp[1] = a*partmp[2] + b;
  }
/*
  double tmp1, tmp2;
  xtoz(partmp[1], partmp[2], tmp1, tmp2);
  if (fabs(tmp2 - x[0]) > 1e-10*fabs(x[0]))
  {
    of << "something is wrong in difternary - two" << endl;
    exit(1);
  }
*/
  f[0] = mu1(partmp);
  f[1] = mu2(partmp);
  f[2] = mu3(partmp);
  if (print == all_iter)
    of << setw(wdth) << partmp[1]
       << setw(wdth) << partmp[2]
       << setw(wdth) << partmp[3]
       << setw(wdth) << partmp[4]
       << setw(wdth) << 1. - partmp[4] - partmp[3]
       << setw(wdth) << sqrt(pow(f[0],2) + pow(f[1],2) + pow(f[2],2))
       << endl;
}
}

void ternary_solve(double* reg)
{
  Ttmp = reg[0];
  xCu1calc = reg[1];
  xY1calc = reg[2];
  xCu2calc = reg[3];
  xY2calc = reg[4];
  if (!cryt)
  {
    cryt_solve();
    cryt = true;
  }
  double x[3];
  double tmp;
  partmp[0] = Ttmp;
  xtoz(xCu1calc, xY1calc, tmp, x[0]);
  xtoz(xCu2calc, xY2calc, x[1], x[2]);
  int n = 3;
  int m = n;
  int nsig = 10;
  double eps = 0;
  double delta = 0;
  int maxfn = 300;
  int iopt = 1;
  int ixjac;
  int ier;
  int infer;
  double parm[4];
  double ssq;
  ixjac = m;
  double* f = new double[m];
  double* xjac = new double[n*m];
  double* xjtj = new double[n*(n+1)/2];
  double* work = new double[5*n + 2*m + (n+1)*n/2];
  double zxssq_machine_eps = 12;
  zxssq1_(difternary, &m, &n, &nsig, &eps, &delta, &maxfn, &iopt,
        parm, x, &ssq, f, xjac, &ixjac, xjtj, work, &infer, &ier, &zxssq_machine_eps);
  difternary(x, &m, &n, f);
  delete f;
  delete xjac;
  delete xjtj;
  delete work;
  if (ssq > 1e-4)
  {
    if (print == all_iter)
    {
      of << "ternary solve - ";
      print_zxssq_1(infer, ier, work);
      of << setw(wdth) << partmp[1]
         << setw(wdth) << partmp[2]
         << setw(wdth) << partmp[3]
         << setw(wdth) << partmp[4]
         << setw(wdth) << 1. - partmp[4] - partmp[3]
         << setw(wdth) << sqrt(ssq)
         << endl;
    }
    xCu1calc = xCu2calc = xCu_cryt;
    xY1calc = xY2calc = xY_cryt;
  }
  else
  {
    xCu1calc = partmp[1];
    xY1calc = partmp[2];
    xCu2calc = partmp[3];
    xY2calc = partmp[4];
    if (fabs(xCu1calc - xCu2calc) < 0.005)
    {
      xCu1calc = xCu2calc = xCu_cryt;
      xY1calc = xY2calc = xY_cryt;
    }
  }
}

double delx(double* reg)
{
  double& xCu1 = reg[1];
  double& xY1  = reg[2];
  double& xCu2 = reg[3];
  double& xY2  = reg[4];
  xCuctr = (xCu2 + xCu1)/2.;
  xYctr = (xY2 + xY1)/2.;
  switch (x_switch)
  {
    case 1:
      ternary_solve(reg);
      x_switch = 2;
      return xCu1 - xCu1calc;
    case 2:
      x_switch = 3;
      return xY1 - xY1calc;
    case 3:
      x_switch = 4;
      return xCu2 - xCu2calc;
    case 4:
      x_switch = 1;
      return xY2 - xY2calc;
    default:
      of << "something is wrong with DELX - exit" << endl;
      exit(1);
      return 0.;
  }
}

void init_res()
{
  res.insert(residual("mu1",        mu1,       5, 0x7));
  res.insert(residual("mu2",        mu2,       5, 0x7));
  res.insert(residual("mu3",        mu3,       5, 0x7));
  res.insert(residual("delx",      delx,       5, 0x7));
}

void post_analysis()
{
  int N = 8;
  double data[8][4] =
    {0.154,   0.017 ,  0.583,   0.413,
     0.342,   0.0083,  0.683,   0.308,
     0.433,   0.017 ,  0.733,   0.25,
     0.508,   0.0083,  0.775,   0.217,
     0.55 ,   0.017 ,  0.775,   0.208,
     0.592,   0.025 ,  0.825,   0.158,
     0.717,   0.042 ,  0.883,   0.1,
     0.767,   0.05  ,  0.9  ,   0.083};
  double par[5];
  par[0] = Tsys;
  int i, j;
  of << endl << "******     post analysis    *****" << endl;
  of << "literature - normal" << endl;
  of << "  xCu1    xY1     xCu2    xY2    xCu3    xY3    xCu4    xY4 \
   xCu5    xY5  xCu6    xY6   xCu7    xY7   xCu8    xY8" << endl;
  for (j = 0; j < 4; j +=2)
  {
    for (i = 0; i < N; i++)
    {
      of << setw(wdth) << data[i][j]
         << setw(wdth) << data[i][j+1];
    }
    of << endl;
  }
  of << endl;
  of << "literature " << endl;
  of << "  xCu1    xY1     xCu2    xY2    xCu3    xY3    xCu4    xY4 \
   xCu5    xY5  xCu6    xY6   xCu7    xY7   xCu8    xY8" << endl;
  for (j = 0; j < 4; j +=2)
  {
    for (i = 0; i < N; i++)
    {
      of << setw(wdth) << sqrt(3.)/2.*data[i][j]
         << setw(wdth) << data[i][j+1] + 1./2.*data[i][j];
    }
    of << endl;
  }
  of << endl;
  for (i = 0; i < N; i++)
  {
    par[1] = data[i][0];
    par[2] = data[i][1];
    par[3] = data[i][2];
    par[4] = data[i][3];
    xCuctr = (par[1] + par[3])/2.;
    xYctr = (par[2] + par[4])/2.;
    ternary_solve(par);
    data[i][0] = xCu1calc;
    data[i][1] = xY1calc;
    data[i][2] = xCu2calc;
    data[i][3] = xY2calc;
  }
  of << "calc - connodes - normal" << endl;
  of << "  xCu1    xY1     xCu2    xY2    xCu3    xY3    xCu4    xY4 \
   xCu5    xY5  xCu6    xY6   xCu7    xY7   xCu8    xY8" << endl;
  for (j = 0; j < 4; j +=2)
  {
    for (i = 0; i < N; i++)
    {
      of << setw(wdth) << data[i][j]
         << setw(wdth) << data[i][j+1];
    }
    of << endl;
  }
  of << endl;
  of << "calc - connodes" << endl;
  of << "  xCu1    xY1     xCu2    xY2    xCu3    xY3    xCu4    xY4 \
   xCu5    xY5  xCu6    xY6   xCu7    xY7   xCu8    xY8" << endl;
  for (j = 0; j < 4; j +=2)
  {
    for (i = 0; i < N; i++)
    {
      of << setw(wdth) << sqrt(3.)/2.*data[i][j]
         << setw(wdth) << data[i][j+1] + 1./2.*data[i][j];
    }
    of << endl;
  }
  of << endl;
  of << "calc - lines" << endl;
  of << "   xCu1          xY1           xCu2           xY2 " << endl;
  int ncalc = 100;
  xCu1calc = 0.01;
  xY1calc = 1e-5;
  xCu2calc = 0.09;
  xY2calc = 1 - 0.09 - 1e-5;
  of << setw(wdth) << 0.
     << setw(wdth) << 0.
     << setw(wdth) << 0.
     << setw(wdth) << 1.
     << endl;
  double Cu1_Y;
  double Cu2_Y;
  double Y1_Y;
  double Y2_Y;
  for (i = 0; i < ncalc; i++)
  {
    xCuctr = 0.05 + (xCu_cryt - 0.05)/ncalc*i;
    xYctr = 0.5 + (xY_cryt - 0.5)/ncalc*i;
    par[1] = xCu1calc;
    par[2] = xY1calc;
    par[3] = xCu2calc;
    par[4] = xY2calc;
    ternary_solve(par);
    of << setw(wdth) << sqrt(3.)/2.*xCu1calc
       << setw(wdth) << xY1calc + 1./2.*xCu1calc
       << setw(wdth) << sqrt(3.)/2.*xCu2calc
       << setw(wdth) << xY2calc + 1./2.*xCu2calc
       << endl;
    if (fabs(Dmixg2(Tsys, xCu2calc, xY2calc) - GY) < 250.)
    {
      cout << "Y - yes" << endl;
      Cu1_Y = xCu1calc;
      Cu2_Y = xCu2calc;
      Y1_Y = xY1calc;
      Y2_Y = xY2calc;
    }
  }
  of << setw(wdth) << sqrt(3.)/2.*xCu_cryt
     << setw(wdth) << xY_cryt + 1./2.*xCu_cryt
     << setw(wdth) << sqrt(3.)/2.*xCu_cryt
     << setw(wdth) << xY_cryt + 1./2.*xCu_cryt
     << endl;
  of << "Y" << endl;
  of << "   xCu1          xY1           xCu2           xY2 " << endl;
  of << setw(wdth) << sqrt(3.)/2.*Cu1_Y
     << setw(wdth) << Y1_Y + 1./2.*Cu1_Y
     << setw(wdth) << sqrt(3.)/2.*Cu2_Y
     << setw(wdth) << Y2_Y + 1./2.*Cu2_Y
     << endl;
}