/*
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 <strstream>
#include <fstream.h>
#include <iomanip.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <string>
#include <donlp.h>
#include "bacuy.h"

/*
//int MAIN__(void);
double xBa = 0.3, xCu = 0.3, xY = 0.4;
double T = 1000;
*/

ofstream of, lf;
ifstream df;
int print;
int interactive;
int command_line;
int init_LP = 0;
double hBaY = 8.596e+04;
double sBaY = 0.;

static double Tfinish, Tstep;

static char* LICENSE = "\n\
Evgenii Rudnyi 1995, (c) All rights reserved \n\n\
Chemistry Department, Moscow State University, 119899 Moscow, Russia \n\
tel. 7(095)939-5452, fax 7(095)939-1205, 7(095)932-8846 \n\
http://www.chem.msu.su/people/rudnyi/welcome.html \n\
rudnyi@mch.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 something like this to another system - we can have a deal. \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.\n";

static const char* HELP ="\n\
Evgenii Rudnyi 1995, (c) All rights reserved \n\
beta-release, version 0.5 (20 April 1995) \n\
computing the equilibrium composition in the ternary Ba-Cu-Y mixture \n\n\
    bacuy_eq [-options] datafile \n\
      options \n\
      -a  about,  -l  license,  -h  help  \n\
      -k  use Konetzki's Ba-Y interaction parameter \n\
      -p  no DONLP output (default) \n\
      -p1 small DONLP output (only results) \n\
      -p2 medium DONLP output (follow each iteration) \n\
      -p3 large DONLP output (free space at the disk) \n\
      -g N set number of grids for LP (3 < N < 50, default 50) \n\
      -e  uses LP estimates just once (default) \n\
      -e1 uses LP estimates for each run \n\n\
    bacuy_eq -i [output] \n\
      interactive mode (results to output.log and output.lst) \n\n\
    bacuy_eq -c xCu xY Tstart [Tfinish Tstep output] \n\
      command-line mode (results to output.log and output.lst) \n";

static const char* ABOUT ="\n\
Evgenii Rudnyi 1995, (c) All rights reserved \n\n\
Chemistry Department, Moscow State University, 119899 Moscow, Russia \n\
tel. 7(095)939-5452, fax  7(095)939-1205 7(095)932-8846 \n\
http://www.chem.msu.su/people/rudnyi/welcome.html \n\
rudnyi@mch.chem.msu.su \n\n\
The program takes Gibbs energies of the phases in the Ba-Cu-Y system \n\
estimated in \n\
E.B. Rudnyi, G.F. Voronin, 1995, to be submitted to CALPHAD \n\
and computes equilibrium composition in the temperature range 400-4000 K. \n\n\
It employs DONLP (DO NonLinear Programming) subroutine \n\
(available at netlib/opt/donlp.shar) \n\
developed by P. Spellucci (spellucci@mathematik.th-darmstadt.de) \n\
with initial estimates obtained by linear programing. \n\n\
See file README for more info. \n\n\
International Science Foundation is acknowledged for generous financial \n\
support (grant MRL-000). \n";

static const char* WRONGOPTION = "\nNo such option \n\n";
static const char* WRONGARG = "\nWrong set of arguments \n\n";
static const char* NOTEXIST = " - such file does not exist \n";
static const char* CANTCREATE = "\nCan not create a file - ";

void open_streams(int argc, char* argv[])
{
  string datafile, logfile, outfile;
  string spoint(".");
  if (argc == 1)
  {
    cout << HELP;
    exit(1);
  }
  for (int i = 1; i < argc; i++)
  {
    if (argv[i][0] == '-' || argv[i][0] == '/')
    {
      switch (tolower(argv[i][1]))
      {
        case 'h' : case '?' :
          cout << HELP;
          exit(1);
        case 'l' :
          cout << LICENSE;
          exit(1);
        case 'a' :
          cout << ABOUT;
          exit(1);
        case 'p' :
          if(argv[i][2])
            print = atoi(argv[i]+2);
          else
            print = 0;
          if (print)
            o6stpa_1.silent = 0;
          else
            o6stpa_1.silent = 1;
          break;
        case 'k' :
          hBaY = 200700.;
          sBaY = 22.3;
          break;
        case 'i' :
          interactive = 1;
          command_line = 0;
          break;
        case 'c' :
          command_line = 1;
          interactive = 0;
          if (i + 3 >= argc)
          {
            cout << WRONGARG;
            exit(1);
          }
          else
          {
            i++;
            xCu = atof(argv[i]);
            i++;
            xY = atof(argv[i]);
            i++;
            T = atof(argv[i]);
            if (i + 1 < argc)
            {
              i++;
              Tfinish = atof(argv[i]);
            }
            else
              Tfinish = T;
            if (i + 1 < argc)
            {
              i++;
              Tstep = atof(argv[i]);
            }
            else
              Tstep = Tfinish - T;
            if (!Tstep) Tstep = 1.;
          }
          break;
        case 'g' :
          if (i + 1 >= argc)
          {
            cout << WRONGARG;
            exit(1);
          }
          else
          {
            i++;
            grid = atol(argv[i]);
            if (grid < 3) grid = 3;
            if (grid > 50) grid = 50;
          }
          break;
        case 'e' :
          if(argv[i][2] == '1')
            init_LP = 1;
          else
            init_LP = 0;
          break;
        default:
          cout << WRONGOPTION;
          exit(1);
      }
    }
    else if (!datafile.length())
      datafile = argv[i];
    else
    {
      cout << WRONGARG;
      exit(1);
    }
  }
  if (!datafile.length() && !(command_line || interactive))
  {
    cout << WRONGARG;
    exit(1);
  }
  if (datafile.length())
  {
		int pos = datafile.find(spoint);
		string base;
    if (pos == string::npos)
		{
			base = datafile;
      datafile = datafile + ".dat";
		}
		else
			base = datafile.substr(0, pos);
    if (!(command_line || interactive))
    {
      df.open(datafile.c_str(), ios::in);
      if (!df)
      {
        cout << endl << datafile << NOTEXIST;
        exit(1);
      }
    }
    logfile = base + ".log";
    lf.open(logfile.c_str());
    if (!lf)
    {
      cout << CANTCREATE << logfile << endl;
      exit(1);
    }
    outfile = base + ".lst";
    of.open(outfile.c_str());
    if (!of)
    {
      cout << CANTCREATE << outfile << endl;
      exit(1);
    }
  }
  cout.precision(4);
  of.precision(4);
  lf.precision(4);
}

void donlp_estimates(void)
{
  double xCu1, xY1;
  cout << "     ";
  donlp_();
  if (singular)
  {
    o6itin_1.optite = 9;
    cout << setw(2) << o6itin_1.optite;
  }
  else
  {
    cout << setw(2) << o6itin_1.optite;
    int iter = 0;
    while (o6itin_1.optite && iter++ < 10)
    {
      cout << "     ";
      donlp_();
      if (singular)
      {
        o6itin_1.optite = 9;
        break;
      }
      cout << setw(2) << o6itin_1.optite;
    }
  }
  if (o6itin_1.optite)
  {
    cout << endl;
    long oldgrid = grid;
    grid = 50;
    lp_estimates();
    grid = oldgrid;
  }
  cout << endl;
  cout << "DONLP results,      T=" << T << ", ier="
       << o6itin_1.optite << ", G=" << o6xdat_1.fx << endl;
  lf   << "DONLP results,      T=" << T << ", ier="
       << o6itin_1.optite << ", G=" << o6xdat_1.fx << endl;
  cout << "           n(mole)           xBa            xCu            xY" << endl;
  lf   << "           n(mole)           xBa            xCu            xY" << endl;
  if (o6xdat_1.x[0] > 1e-10) cout << "Ba(s)    " << setw(10) << o6xdat_1.x[0] << endl;
  if (o6xdat_1.x[1] > 1e-10) cout << "Cu(s)    " << setw(10) << o6xdat_1.x[1] << endl;
  if (o6xdat_1.x[2] > 1e-10) cout << "Y(s)     " << setw(10) << o6xdat_1.x[2] << endl;
  if (o6xdat_1.x[3] > 1e-10) cout << "BaCu     " << setw(10) << o6xdat_1.x[3] << endl;
  if (o6xdat_1.x[4] > 1e-10) cout << "BaCu13   " << setw(10) << o6xdat_1.x[4] << endl;
  if (o6xdat_1.x[5] > 1e-10) cout << "Cu6Y     " << setw(10) << o6xdat_1.x[5] << endl;
  if (o6xdat_1.x[6] > 1e-10) cout << "Cu4Y     " << setw(10) << o6xdat_1.x[6] << endl;
  if (o6xdat_1.x[7] > 1e-10) cout << "Cu7Y2    " << setw(10) << o6xdat_1.x[7] << endl;
  if (o6xdat_1.x[8] > 1e-10) cout << "Cu2Y     " << setw(10) << o6xdat_1.x[8] << endl;
  if (o6xdat_1.x[9] > 1e-10) cout << "CuY      " << setw(10) << o6xdat_1.x[9] << endl;
  if (o6xdat_1.x[0] > 1e-10) lf   << "Ba(s)    " << setw(10) << o6xdat_1.x[0] << endl;
  if (o6xdat_1.x[1] > 1e-10) lf   << "Cu(s)    " << setw(10) << o6xdat_1.x[1] << endl;
  if (o6xdat_1.x[2] > 1e-10) lf   << "Y(s)     " << setw(10) << o6xdat_1.x[2] << endl;
  if (o6xdat_1.x[3] > 1e-10) lf   << "BaCu     " << setw(10) << o6xdat_1.x[3] << endl;
  if (o6xdat_1.x[4] > 1e-10) lf   << "BaCu13   " << setw(10) << o6xdat_1.x[4] << endl;
  if (o6xdat_1.x[5] > 1e-10) lf   << "Cu6Y     " << setw(10) << o6xdat_1.x[5] << endl;
  if (o6xdat_1.x[6] > 1e-10) lf   << "Cu4Y     " << setw(10) << o6xdat_1.x[6] << endl;
  if (o6xdat_1.x[7] > 1e-10) lf   << "Cu7Y2    " << setw(10) << o6xdat_1.x[7] << endl;
  if (o6xdat_1.x[8] > 1e-10) lf   << "Cu2Y     " << setw(10) << o6xdat_1.x[8] << endl;
  if (o6xdat_1.x[9] > 1e-10) lf   << "CuY      " << setw(10) << o6xdat_1.x[9] << endl;
  if (o6xdat_1.x[10] > 1e-10)
  {
    cout << "L1       " << setw(10) << o6xdat_1.x[10];
    lf   << "L1       " << setw(10) << o6xdat_1.x[10];
    ztox(xCu1, xY1, o6xdat_1.x[12], o6xdat_1.x[13]);
    cout << setw(15) << 1. - xCu1 - xY1 << setw(15) << xCu1
         << setw(15) << xY1 << endl;
    lf   << setw(15) << 1. - xCu1 - xY1 << setw(15) << xCu1
         << setw(15) << xY1 << endl;
  }
  if (o6xdat_1.x[11] > 1e-10)
  {
    cout << "L2       " << setw(10) << o6xdat_1.x[11];
    lf   << "L2       " << setw(10) << o6xdat_1.x[11];
    ztox(xCu1, xY1, o6xdat_1.x[14], o6xdat_1.x[15]);
    cout << setw(15) << 1. - xCu1 - xY1 << setw(15) << xCu1
         << setw(15) << xY1 << endl;
    lf   << setw(15) << 1. - xCu1 - xY1 << setw(15) << xCu1
         << setw(15) << xY1 << endl;
  }
  cout << endl;
  lf   << endl;
  if (of)
  {
    of << setw(2) << o6itin_1.optite
       << setw(10) << xBa
       << setw(10) << xCu
       << setw(10) << xY
       << setw(10) << T;
		int i;
    for (i = 0; i < 11; i++)
      of << setw(12) << o6xdat_1.x[i];
    ztox(xCu1, xY1, o6xdat_1.x[12], o6xdat_1.x[13]);
    of << setw(12) << 1. - xCu1 - xY1 << setw(12) << xCu1
       << setw(12) << xY1;
    of << setw(12) << o6xdat_1.x[i];
    ztox(xCu1, xY1, o6xdat_1.x[14], o6xdat_1.x[15]);
    of << setw(12) << 1. - xCu1 - xY1 << setw(12) << xCu1
       << setw(12) << xY1
       << endl;
  }
}

int main(int argc, char *argv[])
{
  o6stpa_1.silent = 1;
  df.open("fort.10", ios::in);
  if (df)
  {
    df.close();
    of.open("fort.10");
    of.close();
  }
  df.open("fort.20", ios::in);
  if (df)
  {
    df.close();
    of.open("fort.20");
    of.close();
  }
  open_streams(argc, argv);
  of << "ier  xBa       xCu        xY        T     "
     << "    Ba_s        Cu_s         Y_s        BaCu       BaCu13   "
     << "    Cu6Y        Cu4Y        Cu7Y2       Cu2Y         CuY    "
     << "     L1         xBa1        xCu1        xY1          L2     "
     << "    xBa2        xCu2        xY2"
     << endl;
  if (command_line)
  {
    xBa = 1. - xCu - xY;
    cout << "xBa=" << xBa << ", xCu=" << xCu << ", xY=" << xY << endl;
    lf   << "xBa=" << xBa << ", xCu=" << xCu << ", xY=" << xY << endl;
    if (xBa < 0.005 || xCu < 0.005 || xY < 0.005)
    {
      cout << "composition is out of range" << endl;
      lf   << "composition is out of range" << endl;
      exit(1);
    }
    int first = 0;
    for (; (Tstep > 0.) ? T < Tfinish + Tstep/2.: T > Tfinish + Tstep/2.;
           T += Tstep)
    {
      if (T > 4000. || T < 300.)
      {
        cout << "temperature is out of range" << endl;
        lf << "temperature is out of range" << endl;
        of << "temperature is out of range" << endl;
        exit(1);
      }
      if (init_LP || !first)
      {
        lp_estimates();
        first++;
      }
      donlp_estimates();
    }
  }
  else if (interactive)
  {
    int first = 0;
    while (1)
    {
      cout << "Cu, Y? ";
      cin >> xCu >> xY;
      xBa = 1. - xCu - xY;
      cout << "xBa=" << xBa << ", xCu=" << xCu << ", xY=" << xY << endl;
      lf   << "xBa=" << xBa << ", xCu=" << xCu << ", xY=" << xY << endl;
      if (xBa < 0.005 || xCu < 0.005 || xY < 0.005)
      {
        cout << "composition is out of range" << endl;
        lf   << "composition is out of range" << endl;
        exit(1);
      }
      while(1)
      {
        cout << "T? ";
        cin >> T;
        if (T > 4000. || T < 300.)
        {
          cout << "temperature is out of range" << endl;
          lf << "temperature is out of range" << endl;
          break;
        }
        if (init_LP || !first)
        {
          lp_estimates();
          first++;
        }
        donlp_estimates();
      }
    }
  }
  else
  {
    char line[81];
    int first = 0;
    while (df.peek() != EOF)
    {
      df.getline(line, 80);
      if (line[0] != '#')
      {
        istrstream in(line);
        in >> xCu >> xY >> T;
        xBa = 1. - xCu - xY;
        cout << "xBa=" << xBa << ", xCu=" << xCu << ", xY=" << xY << endl;
        lf   << "xBa=" << xBa << ", xCu=" << xCu << ", xY=" << xY << endl;
        if (xBa < 0.005 || xCu < 0.005 || xY < 0.005)
        {
          cout << "composition is out of range" << endl;
          lf   << "composition is out of range" << endl;
          of   << "composition is out of range" << endl;
          continue;
        }
        if (T > 4000. || T < 300.)
        {
          cout << "temperature is out of range" << endl;
          lf << "temperature is out of range" << endl;
          of << "temperature is out of range" << endl;
          continue;
        }
        if (init_LP || !first)
        {
          lp_estimates();
          first++;
        }
        donlp_estimates();
      }
    }
  }
  return 0;
}