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

Program bacuy_eq: Computing Phase Diagram of Ba-Cu-Y

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

#include <iostream.h>
#include <iomanip.h>
#include <math.h>
#include <stdlib.h>
#include <donlp.h>
#include "bacuy.h"

extern "C" int zx3lp_(double *a, long *ia, double *b,
        double *c, long *n, long *m1, long *m2, double *s,
        double *psol, double *dsol, double *rw, long *iw,
        long *ier);

double T;
double xBa, xCu, xY;
static double liq_small = 1e-15;
long grid = 50l;

void init_solution(long n, double *psol, double (*a)[5], int nophase)
{
  int ii;
  for (ii = 10; ii < n; ii++)
    if (psol[ii] > 1e-10)
    {
      if (!o6xdat_1.x[10])
      {
        if (!nophase)
          o6xdat_1.x[10] = psol[ii];
        else
          o6xdat_1.x[10] = liq_small;
        o6xdat_1.x[12] = a[ii][1l];
        o6xdat_1.x[13] = a[ii][2l];
      }
      else if (!o6xdat_1.x[11])
      {
        if (!nophase)
          o6xdat_1.x[11] = psol[ii];
        else
          o6xdat_1.x[11] = liq_small;
        o6xdat_1.x[14] = a[ii][1l];
        o6xdat_1.x[15] = a[ii][2l];
      }
      else
      {
        if (fabs(o6xdat_1.x[15] - o6xdat_1.x[13]) <
            fabs(a[ii][2l] - o6xdat_1.x[13]))
        {
          if (fabs(o6xdat_1.x[15] - o6xdat_1.x[13]) <
              fabs(o6xdat_1.x[15] - a[ii][2l]))
          {
            o6xdat_1.x[10] += o6xdat_1.x[11];
            o6xdat_1.x[11] = 0.;
          }
          if (!nophase)
            o6xdat_1.x[11] += psol[ii];
          else
            o6xdat_1.x[11] = liq_small;
          o6xdat_1.x[14] = a[ii][1l];
          o6xdat_1.x[15] = a[ii][2l];
        }
        else if (!nophase)
        {
          if (fabs(a[ii][2l] - o6xdat_1.x[13]) <
              fabs(a[ii][2l] - o6xdat_1.x[15]))
            o6xdat_1.x[10] += psol[ii];
          else
            o6xdat_1.x[11] += psol[ii];
        }
      }
/*
        cout << setw(10) << a[ii][0l]
             << setw(10) << a[ii][1l]
             << setw(10) << a[ii][2l]
             << setw(10) << psol[ii]
             << endl;
*/
    }
}

int lp_estimates()
{
  double xBai, xCui, xYi;
  double step = 1./grid;
  long m1 = 0l;
  long m2 = 3l;
  long n = 10l + (grid+2)*(grid+1l)/2l;
  //                       Ba  Cu   Y
  double (*a)[5] = new double[n][5];
  double tmp[10][5] =
                         { 1.,  0., 0., 0., 0.,    // Ba(s)
                           0.,  1., 0., 0., 0.,    // Cu(s)
                           0.,  0., 1., 0., 0.,    // Y(s)
                           1.,  1., 0., 0., 0.,    // BaCu(s)
                           1., 13., 0., 0., 0.,    // BaCu13(s)
                           0.,  6., 1., 0., 0.,    // Cu6Y(s)
                           0.,  4., 1., 0., 0.,    // Cu4Y(s)
                           0.,  7., 2., 0., 0.,    // Cu72Y(s)
                           0.,  2., 1., 0., 0.,    // Cu2Y(s)
                           0.,  1., 1., 0., 0.};    // CuY(s)
  memcpy(a, tmp, sizeof(tmp));
  double* b = new double[m1+m2+2];
  double* c = new double[n];
  double* psol = new double[n];
  double* dsol = new double[m1+m2+2];
  double* rw = new double[(m1+m2+2)*(m1+m2+2)+3*m1+2*m2+4];
  long* iw = new long[2*m2+3*m1+4];
  if (!a || !b || !c || !psol || !dsol || !rw || !iw)
  {
    cout << "not enough memory for bacuy_lp" << endl;
    lf   << "not enough memory for bacuy_lp" << endl;
    exit(1);
  }
  long dim = m1+m2+2;
  double s;
  long ier;
  long i,j;
//  memory t;
/*
  for (i = 0; i < dim; i++)
  {
    b[i] = 0.;
    dsol[i] = 0.;
  }
  for (i = 0; i < (m1+m2+2)*(m1+m2+2)+3*m1+2*m2+4; i++)
    rw[i] = 0.;
  for (i = 0; i < 2*m2+3*m1+4; i++)
    iw[i] = 0;
*/
  b[0] = xBa;
  b[1] = xCu;
  b[2] = xY;
  c[0l] = -gBas(T);
  c[1l] = -gCus(T);
  c[2l] = -gYs(T);
  c[3l] = -gBaCu(T);
  c[4l] = -gBaCu13(T);
  c[5l] = -gCu6Y(T);
  c[6l] = -gCu4Y(T);
  c[7l] = -gCu7Y2(T);
  c[8l] = -gCu2Y(T);
  c[9l] = -gCuY(T);
/*
  for (i = 0; i < 10; i++)
    c[i] = -1e30;
*/
  long ii = 10l;
  double small = 0.0001;
  for (i = 0; i <= grid; i++)
  {
    if (i == 0)
      xCui = small;
    else if (i == grid)
      xCui = 1. - small;
    else
      xCui = i*step;
    for (j = 0; j < grid - i + 1; j++)
    {
      if (i == grid)
        xYi = small/2.;
      else if (j == 0)
        xYi = small;
      else if (j == (grid - i))
        xYi = 1. - xCui - small;
      else
        xYi = j*step;
      xBai = 1. - xCui - xYi;
      c[ii] = -(xBai*gBal(T, xCui, xYi) + xCui*gCul(T, xCui, xYi)
                                  + xYi*gYl(T, xCui, xYi));
/*
      cout << setw(3) << i
           << setw(3) << j
           << setw(5) << ii
           << setw(9) << xCui
           << setw(9) << xYi
           << setw(9) << xBai
           << setw(9) << c[ii]
           << endl;
*/
      a[ii][0l] = xBai;
      a[ii][1l] = xCui;
      a[ii][2l] = xYi;
//      a[ii][3l] = 0.;
//      a[ii][4l] = 0.;
//      psol[ii] = 0.;
      ii++;
      if (i == grid)
        break;
    }
//    cout << setw(2) << i;
  }
//  cout << N << endl;
//  cout.flush();
  if (ii != n)
  {
    cout << ii << " and " << n << " - something wrong in bacuy_lp" << endl;
    lf   << ii << " and " << n << " - something wrong in bacuy_lp" << endl;
    exit(1);
  }
  zx3lp_((double*)a, &dim, b, c, &n, &m1, &m2, &s,
         psol, dsol, rw, iw, &ier);
/*
  for (j = 0; j < m1+m2; j++)
  {
    for (i = 0; i < n; i++)
      cout << setw(5) << a[i][j];
    cout << setw(5) << b[j] << endl;
  }
*/
  cout << "initial LP guesses, T=" << T
       << ", ier=" << ier << ", G=" << -s << ", grids = " << grid << endl;
  lf   << "initial LP guesses, T=" << T
       << ", ier=" << ier << ", G=" << -s << ", grids = " << grid << endl;
  for (i = 0; i < 10; i++)
  {
    o6xdat_1.x[i] = psol[i];
  }
  if (psol[0]) cout << "Ba(s)    " << setw(10) << psol[0] << endl;
  if (psol[1]) cout << "Cu(s)    " << setw(10) << psol[1] << endl;
  if (psol[2]) cout << "Y(s)     " << setw(10) << psol[2] << endl;
  if (psol[3]) cout << "BaCu     " << setw(10) << psol[3] << endl;
  if (psol[4]) cout << "BaCu13   " << setw(10) << psol[4] << endl;
  if (psol[5]) cout << "Cu6Y     " << setw(10) << psol[5] << endl;
  if (psol[6]) cout << "Cu4Y     " << setw(10) << psol[6] << endl;
  if (psol[7]) cout << "Cu7Y2    " << setw(10) << psol[7] << endl;
  if (psol[8]) cout << "Cu2Y     " << setw(10) << psol[8] << endl;
  if (psol[9]) cout << "CuY      " << setw(10) << psol[9] << endl;
  if (psol[0]) lf   << "Ba(s)    " << setw(10) << psol[0] << endl;
  if (psol[1]) lf   << "Cu(s)    " << setw(10) << psol[1] << endl;
  if (psol[2]) lf   << "Y(s)     " << setw(10) << psol[2] << endl;
  if (psol[3]) lf   << "BaCu     " << setw(10) << psol[3] << endl;
  if (psol[4]) lf   << "BaCu13   " << setw(10) << psol[4] << endl;
  if (psol[5]) lf   << "Cu6Y     " << setw(10) << psol[5] << endl;
  if (psol[6]) lf   << "Cu4Y     " << setw(10) << psol[6] << endl;
  if (psol[7]) lf   << "Cu7Y2    " << setw(10) << psol[7] << endl;
  if (psol[8]) lf   << "Cu2Y     " << setw(10) << psol[8] << endl;
  if (psol[9]) lf   << "CuY      " << setw(10) << psol[9] << endl;
  o6xdat_1.x[10] = 0.;
  o6xdat_1.x[11] = 0.;
  init_solution(n, psol, a, 0);
  if (!o6xdat_1.x[11] && o6xdat_1.x[10])
  {
    o6xdat_1.x[11] = liq_small;
    o6xdat_1.x[14] = o6xdat_1.x[12];
    o6xdat_1.x[15] = o6xdat_1.x[13];
  }
  if (!o6xdat_1.x[10])
  {
    long nnew = n - 10;
    zx3lp_((double*)a + 5*10, &dim, b, c + 10, &nnew, &m1, &m2, &s,
           psol + 10, dsol, rw, iw, &ier);
//    cout << "without solids" << endl;
    init_solution(n, psol, a, 1);
//    cout << "s   " << s << endl;
  }
  if (o6xdat_1.x[10] > 1e-10)
  {
    cout << "L1       " << setw(10) << o6xdat_1.x[10];
    cout << setw(15) << 1. - o6xdat_1.x[12] - o6xdat_1.x[13]
         << setw(15) << o6xdat_1.x[12]
         << setw(15) << o6xdat_1.x[13]
         << endl;
    lf   << "L1       " << setw(10) << o6xdat_1.x[10];
    lf   << setw(15) << 1. - o6xdat_1.x[12] - o6xdat_1.x[13]
         << setw(15) << o6xdat_1.x[12]
         << setw(15) << o6xdat_1.x[13]
         << endl;
  }
  if (o6xdat_1.x[11] > 1e-10)
  {
    cout << "L2       " << setw(10) << o6xdat_1.x[11];
    cout << setw(15) << 1. - o6xdat_1.x[14] - o6xdat_1.x[15]
         << setw(15) << o6xdat_1.x[14]
         << setw(15) << o6xdat_1.x[15]
         << endl;
    lf   << "L2       " << setw(10) << o6xdat_1.x[11];
    lf   << setw(15) << 1. - o6xdat_1.x[14] - o6xdat_1.x[15]
         << setw(15) << o6xdat_1.x[14]
         << setw(15) << o6xdat_1.x[15]
         << endl;
  }
  xtoz(o6xdat_1.x[12], o6xdat_1.x[13], o6xdat_1.x[12], o6xdat_1.x[13]);
  xtoz(o6xdat_1.x[14], o6xdat_1.x[15], o6xdat_1.x[14], o6xdat_1.x[15]);
  delete a;
  delete b;
  delete c;
  delete psol;
  delete dsol;
  delete rw;
  delete iw;
  return ier;
}

/*
int main()
{
//  memory test;
  cout << "Ba, Cu, Y? ";
  cin >> xBa >> xCu >> xY;
  if (!xBa || !xCu || !xY)
    return 1;
  while (T != 1005.)
  {
    cout << "T? ";
    cin >> T;
//    test();
    lp_estimates();
  }
  return 0;
}
*/