/*
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 <string.h>
#include <stdlib.h>

#include "bacuy.h"
#include "donlp.h"

static double c[10];

extern "C"
{

int setup0_(void)
{
  o6stpa_1.analyt = 1;
  o6der_1.epsdif = 0.;
  o6io_1.prou = 10;
  o6io_1.meu = 20;
//  o6stpa_1.silent = 1;
  o6dim_1.n = 16;
  o6dim_1.nh = 3;
//  o6dim_1.ng = 12;
  o6dim_1.ng = 20;
  o6rst_1.nreset = o6dim_1.n;
  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);
  int i, j;
  strcpy(o6id_1.name, "BaCuY_phase_diagram");
  o6par_1.del0 = .01;
  o6par_1.tau0 = 1.;
//  o6par_1.tau = 1.;
/*
  o6xdat_1.x[0] = 0.01;
  o6xdat_1.x[1] = 0.01;
  o6xdat_1.x[2] = 0.01;
  o6xdat_1.x[3] = 0.01;
  o6xdat_1.x[4] = 0.01;
  o6xdat_1.x[5] = 0.01;
  o6xdat_1.x[6] = 0.01;
  o6xdat_1.x[7] = 0.01;
  o6xdat_1.x[8] = 0.01;
  o6xdat_1.x[9] = 0.01;
  o6xdat_1.x[10] = 0.4;
  o6xdat_1.x[11] = 0.4;
  o6xdat_1.x[12] = 0.1;
  o6xdat_1.x[13] = 0.01;
  o6xdat_1.x[14] = 0.1;
  o6xdat_1.x[15] = 0.8;
*/
//  for (i = 0; i <= o6dim_1.nh; i++)
  for (i = 0; i < o6dim_1.nh + o6dim_1.ng + 1; i++)
  {
    o6gri_1.gunit[3*i+0] =  -1;
    o6gri_1.gunit[3*i+1] =  0;
    o6gri_1.gunit[3*i+2] =  0;
  }
  for (i = o6dim_1.nh + 1; i <= o6dim_1.nh + 16; i++)
  {
    o6gri_1.gunit[3*i+0] =  1;
    o6gri_1.gunit[3*i+1] =  i - o6dim_1.nh;
    o6gri_1.gunit[3*i+2] =  1;
  }
  for (i = o6dim_1.nh + 17; i <= o6dim_1.nh + 20; i++)
  {
    o6gri_1.gunit[3*i+0] =  1;
    o6gri_1.gunit[3*i+1] =  i - o6dim_1.nh - 4;
    o6gri_1.gunit[3*i+2] =  -1;
  }
  return 0;
}

int setup_(void)
{
  o6stpa_1.te0 = 0;
  o6stpa_1.te1 = 0;
  switch (print)
  {
    case 3:
      o6stpa_1.te3 = 1;
      o6stpa_1.te2 = 1;
      break;
    case 2:
      o6stpa_1.te3 = 0;
      o6stpa_1.te2 = 1;
      break;
    case 1:
      o6stpa_1.te3 = 0;
      o6stpa_1.te2 = 0;
      break;
  }
  return 0;
}

int ef_(double *x, double *fx)
{
  ++o6cnt_1.cf;
  long double sum = 0.;
  double xCu1, xY1;
  for (int i = 0; i < 10; i++)
    sum += x[i]*c[i];
  ztox(xCu1, xY1, x[12], x[13]);
  sum += x[10]*((1. - xCu1 - xY1)*gBal(T, xCu1, xY1)
                           + xCu1*gCul(T, xCu1, xY1)
                            + xY1*gYl(T, xCu1, xY1));
  ztox(xCu1, xY1, x[14], x[15]);
  sum += x[11]*((1. - xCu1 - xY1)*gBal(T, xCu1, xY1)
                           + xCu1*gCul(T, xCu1, xY1)
                            + xY1*gYl(T, xCu1, xY1));
  *fx = sum;
  cout << "\b\b\b\b\b" << setw(5) << o6cnt_1.cf;
  return 0;
}

int egradf_(double *x, double *gradf)
{
//  cout << "egradf" << endl;
  ++o6cnt_1.cgf;
  double xCu1, xY1;
  double dxCu1, dxY1;
  double dx1dz1, dx1dz2, dx2dz1, dx2dz2;
  for (int i = 0; i < 10; i++)
    gradf[i] = c[i];
  ztox(xCu1, xY1, x[12], x[13]);
  dxdz(dx1dz1, dx1dz2, dx2dz1, dx2dz2, x[12], x[13]);
  gradf[10] = ((1. - xCu1 - xY1)*gBal(T, xCu1, xY1)
                          + xCu1*gCul(T, xCu1, xY1)
                           + xY1*gYl(T, xCu1, xY1));
  dxCu1 = (gCul(T, xCu1, xY1) - gBal(T, xCu1, xY1))*x[10];
  dxY1 = (gYl(T, xCu1, xY1) - gBal(T, xCu1, xY1))*x[10];
  gradf[12] = dxCu1*dx1dz1 + dxY1*dx2dz1;
  gradf[13] = dxCu1*dx1dz2 + dxY1*dx2dz2;
  ztox(xCu1, xY1, x[14], x[15]);
  dxdz(dx1dz1, dx1dz2, dx2dz1, dx2dz2, x[14], x[15]);
  gradf[11] = ((1. - xCu1 - xY1)*gBal(T, xCu1, xY1)
                          + xCu1*gCul(T, xCu1, xY1)
                           + xY1*gYl(T, xCu1, xY1));
  dxCu1 = (gCul(T, xCu1, xY1) - gBal(T, xCu1, xY1))*x[11];
  dxY1 = (gYl(T, xCu1, xY1) - gBal(T, xCu1, xY1))*x[11];
  gradf[14] = dxCu1*dx1dz1 + dxY1*dx2dz1;
  gradf[15] = dxCu1*dx1dz2 + dxY1*dx2dz2;
  return 0;
}

//Ba  Cu   Y
//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)

int eh_(long *i, double *x, double *hxi)
{
//  cout << "eh " << *i << endl;
  ++o6cnt_1.cres[*i - 1];
  double xCu1, xY1, xCu2, xY2;
  ztox(xCu1, xY1, x[12], x[13]);
  ztox(xCu2, xY2, x[14], x[15]);
  switch (*i)
  {
    case 1:  // Ba
      *hxi = x[0] + x[3] + x[4] + x[10]*(1. - xCu1 - xY1)
             + x[11]*(1. - xCu2 - xY2) - xBa;
      break;
    case 2:  // Cu
      *hxi = x[1] + x[3] + 13.*x[4] + 6.*x[5] + 4.*x[6] + 7.*x[7]
             + 2.*x[8] + x[9] + x[10]*xCu1 + x[11]*xCu2 - xCu;
      break;
    case 3:  // Y
      *hxi = x[2] + x[5] + x[6] + 2.*x[7] + x[8] + x[9]
             + x[10]*xY1 + x[11]*xY2 - xY;
      break;
    default:
      cout << "something wrong - eh_" << endl;
      exit(1);
  }
  return 0;
}

int egradh_(long *i, double *x, double *gradhi)
{
//  cout << "egradh " << *i << endl;
  ++o6cnt_1.cgres[*i - 1];
  for (int j = 0; j < o6dim_1.n; j++)
    gradhi[j] = 0.;
  double xCu1, xY1, xCu2, xY2;
  ztox(xCu1, xY1, x[12], x[13]);
  ztox(xCu2, xY2, x[14], x[15]);
  double dx1dz1a, dx1dz2a, dx2dz1a, dx2dz2a;
  dxdz(dx1dz1a, dx1dz2a, dx2dz1a, dx2dz2a, x[12], x[13]);
  double dx1dz1b, dx1dz2b, dx2dz1b, dx2dz2b;
  dxdz(dx1dz1b, dx1dz2b, dx2dz1b, dx2dz2b, x[14], x[15]);
  switch (*i)
  {
    case 1:  // Ba
//      *hxi = x[0] + x[3] + x[4] + x[10]*(1. - xCu1 - xY1)
//             + x[11]*(1. - xCu2 - xY2) - xBa;
      gradhi[0] = 1.;
      gradhi[3] = 1.;
      gradhi[4] = 1.;
      gradhi[10] = 1. - xCu1 - xY1;
      gradhi[11] = 1. - xCu2 - xY2;
      gradhi[12] = -x[10]*(dx1dz1a + dx2dz1a);
      gradhi[13] = -x[10]*(dx1dz2a + dx2dz2a);
      gradhi[14] = -x[11]*(dx1dz1b + dx2dz1b);
      gradhi[15] = -x[11]*(dx1dz2b + dx2dz2b);
      break;
    case 2:  // Cu
//      *hxi = x[1] + x[3] + 13.*x[4] + 6.*x[5] + 4.*x[6] + 7.*x[7]
//             + 2.*x[8] + x[9] + x[10]*xCu1 + x[11]*xCu2 - xCu;
      gradhi[1] = 1.;
      gradhi[3] = 1.;
      gradhi[4] = 13.;
      gradhi[5] = 6.;
      gradhi[6] = 4.;
      gradhi[7] = 7.;
      gradhi[8] = 2.;
      gradhi[9] = 1.;
      gradhi[10] = xCu1;
      gradhi[11] = xCu2;
      gradhi[12] = x[10]*dx1dz1a;
      gradhi[13] = x[10]*dx1dz2a;
      gradhi[14] = x[11]*dx1dz1b;
      gradhi[15] = x[11]*dx1dz2b;
      break;
    case 3:  // Y
//      *hxi = x[2] + x[5] + x[6] + 2.*x[7] + x[8] + x[9]
//             + x[10]*xY1 + x[11]*xY2 - xY;
      gradhi[2] = 1.;
      gradhi[5] = 1.;
      gradhi[6] = 1.;
      gradhi[7] = 2.;
      gradhi[8] = 1.;
      gradhi[9] = 1.;
      gradhi[10] = xY1;
      gradhi[11] = xY2;
      gradhi[12] = x[10]*dx2dz1a;
      gradhi[13] = x[10]*dx2dz2a;
      gradhi[14] = x[11]*dx2dz1b;
      gradhi[15] = x[11]*dx2dz2b;
      break;
    default:
      cout << "something wrong - egradh_" << endl;
      exit(1);
  }
  return 0;
}

int eg_(long *i, double *x, double *gxi)
{
//  cout << "eg " << *i << endl;
  ++o6cnt_1.cres[*i + o6dim_1.nh - 1];
  if (*i <= 12)
    *gxi = x[*i-1];
  else if (*i <= 16)
    *gxi = x[*i-1] + 20.;
  else if (*i <= 20)
    *gxi = 20. - x[*i-1-4];
  else
  {
    cout << "something wrong - eg_" << endl;
    exit(1);
  }
  return 0;
}

int egradg_(long *i, double *x, double *gradgi)
{
//  cout << "egradg " << *i << endl;
  ++o6cnt_1.cgres[*i + o6dim_1.nh - 1];
  for (int j = 0; j < o6dim_1.n; j++)
    gradgi[j] = 0.;
  if (*i <= 16)
    gradgi[*i-1] = 1.;
  else if (*i <= 20)
    gradgi[*i-1-4] = -1.;
  else
  {
    cout << "something wrong - egradg_" << endl;
    exit(1);
  }
  return 0;
}

}