/*
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 = 27;
param par[27];

double& A0 = par[10].x;
double& A1 = par[11].x;
double& A2 = par[12].x;
double& A3 = par[13].x;
double& B0 = par[14].x;
double& B1 = par[15].x;
double& B2 = par[16].x;
double& B3 = par[17].x;
double& C0 = par[23].x;
double& C1 = par[24].x;
double& C2 = par[25].x;
double& C3 = par[26].x;

double g1s(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 g1l(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 + Dmixg1(T, x) : g;
}

double g2s(double& T)
{
// Y - 91Din (based on Hultgren)
  double g;
  if (T < 1500.)
    g = -7347.055+117.532124*T-23.8685*T*log(T)-3.845475e-3*pow(T,2)
        +0.011125e-6*pow(T,3)-16486./T;
                                      // 298.15, 1500 hcp
  else if (T < 1752.)
    g = -15802.62+229.831717*T-40.2851*T*log(T)+6.8095e-3*pow(T,2)
        -1.14182e-6*pow(T,3);
                                      // 1500,   1799 hcp
  else if (T < 1799.)
    g = -10207.724+195.741984*T-35.0201*T*log(T);
                                      // 1752,   1799 bcc
  else
    g = 104813.954-386.167564*T+39.8075986*T*log(T)-19.918739e-3*pow(T,2)
        +0.841308e-6*pow(T,3)-31549963./T;
                                      // 1799,   3700 bcc
  return g;
}

double g2l(double& T, double x)
{
// Y - 91Din (based on Hultgren)
  double g;
  if (T < 1799.)
    g = 3934.121+59.921688*T-14.8146562*T*log(T)-15.623487e-3*pow(T,2)
        +1.442946e-6*pow(T,3)-140695./T;
                                      // 298.15, 1799
  else
    g = -13337.609+258.004539*T-43.0952*T*log(T);
                                      // 1799,   3700
  return (x) ? g + Dmixg2(T, x) : g;
}

/*                                 2        3       7        -1        -9
G := A + B*TE + LOG(TE)*C*TE + D*TE  + EE*TE  + F*TE  + GG*TE   + HH*TE

                    2          3         7          -1           -9
H := A - C*TE - D*TE  - 2*EE*TE  - 6*F*TE  + 2*GG*TE   + 10*HH*TE
*/

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

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

double h2s(double& T)
{
// Y - 91Din (based on Hultgren)
  double h;
  if (T < 1500.)
    h = -7347.055+23.8685*T+3.845475e-3*pow(T,2)
        -2.*0.011125e-6*pow(T,3)-2.*16486./T;
                                      // 298.15, 1500 hcp
  else if (T < 1752.)
    h = -15802.62+40.2851*T-6.8095e-3*pow(T,2)
        +2.*1.14182e-6*pow(T,3);
                                      // 1500,   1799 hcp
  else if (T < 1799.)
    h = -10207.724+35.0201*T;
                                      // 1752,   1799 bcc
  else
    h = 104813.954-39.8075986*T+19.918739e-3*pow(T,2)
        -2.*0.841308e-6*pow(T,3)-2.*31549963./T;
                                      // 1799,   3700 bcc
  return h;
}

double h2l(double& T, double x = 0.)
{
// Y - 91Din (based on Hultgren)
  double h;
  if (T < 1799.)
    h = 3934.121+14.8146562*T+15.623487e-3*pow(T,2)
        -2.*1.442946e-6*pow(T,3)-2.*140695./T;
                                      // 298.15, 1799
  else
    h = -13337.609+43.0952*T;
                                      // 1799,   3700
  return (x) ? h + Dmixh2(T, x) : h;
}

elem e1("Cu", 0., 1358.);
elem e2("Y", 1., 1799.);

int nlp = 5;
linephase lp[5] =
  {linephase("Cu6Y", 1./7., par[0].x, par[1].x, par[18].x),
   linephase("Cu4Y", 0.2, par[2].x, par[3].x, par[19].x),
   linephase("Cu7Y2", 2./9., par[4].x, par[5].x, par[20].x),
   linephase("Cu2Y", 1./3., par[6].x, par[7].x, par[21].x),
   linephase("CuY", 0.5, par[8].x, par[9].x, par[22].x)};

int nnon = 6;
nonvar non[6] = {{1133., 0.093},
                 {1183., 0.12},
                 {1193., 0.25},
                 {1113., 0.28},
                 {1103., 0.42},
                 {1043., 0.67}};

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 lp2_L(double *reg)
{
  subst1 = &lp[2];
  return T_solve(reg);
}

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

double lp4_L(double *reg)
{
  subst1 = &lp[4];
  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 e1_L_lp0_x(double* reg)
{
  e1_L_lp0(reg);
  return reg[0] - non[0].x;
}

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

double lp0_L_lp1_x(double* reg)
{
  lp0_L_lp1(reg);
  return reg[0] - non[1].x;
}

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

double lp2_L_lp3(double* reg)
{
  subst1 = &lp[2];
  subst2 = &lp[3];
  return Tx_solve(reg, 3);
}

double lp2_L_lp3_x(double* reg)
{
  lp2_L_lp3(reg);
  return reg[0] - non[3].x;
}

double lp3_L_lp4(double* reg)
{
  subst1 = &lp[3];
  subst2 = &lp[4];
  return Tx_solve(reg, 4);
}

double lp3_L_lp4_x(double* reg)
{
  lp3_L_lp4(reg);
  return reg[0] - non[4].x;
}

double lp4_L_e2(double* reg)
{
  subst1 = &lp[4];
  subst2 = &e2;
  return Tx_solve(reg, 5);
}

double lp4_L_e2_x(double* reg)
{
  lp4_L_e2(reg);
  return reg[0] - non[5].x;
}

double H(double* reg)
{
  double& H = reg[0];
  double& x = reg[1];
  double& T = reg[2];
  double Hc = DmixH(T, x);
//  return (H - Hc)/fabs(H);
//  return (H - Hc);
  return (H - Hc)/x/(1.-x);
}

double Hsl(double* reg)
{
  double& H = reg[0];
  double& x = reg[1];
  double& T = reg[2];
  double Hc = DmixH(T, x);
  double Hm = h2l(T) - h2s(T);
//  return (H - Hm*x - Hc)/fabs(H);
//  return (H - Hm*x - Hc);
  return (H - Hm*x - Hc)/x/(1.-x);
}

double HTst(double& x)
{
  if (x < 0. || x > 1.)
  {
    warning("value of x is out of range in HTst");
    exit(1);
    return 0.;
  }
  subst1 = &e1;
  subst2 = &e2;
  for (int i = 0; i < nlp; i++)
  {
    if (x <= lp[i].xs)
    {
      subst2 = &lp[i];
      break;
    }
    subst1 = &lp[i];
  }
  double Tst = 298.15;
  return ((x - subst1->xs)*subst2->h(Tst)
          + (subst2->xs - x)*subst1->h(Tst))
         /(subst2->xs - subst1->xs);
}

double HTs(double* reg)
{
  double& H = reg[0];
  double& T = reg[1];
  double& x = reg[2];
  double Hst = HTst(x);
  double Hc = ((x - subst1->xs)*subst2->h(T)
               + (subst2->xs - x)*subst1->h(T))
              /(subst2->xs - subst1->xs) - Hst;
  return H - Hc;
}

double HTl(double* reg)
{
  double& H = reg[0];
  double& T = reg[1];
  double& x = reg[2];
  double Hst = HTst(x);
  double Hc = (1. - x)*h1l(T, x) + x*h2l(T, x) - Hst;
  return H - Hc;
}

double Hlp4(double* reg)
{
  double& H = reg[0];
  double& T = reg[1];
  return H - (lp[4].H + lp[4].C*T);
}

void init_res()
{
// C3 C2 C1 C0 Ce Cd Cc Cb Ca B3 B2 B1 B0 A3 A2 A1 A0 Be Ae Bd Ad Bc Ac Bb Ab Ba Aa
// 27 26 25 24 23 22 21 20 19 18 17 16 15 14 13 12 11 10  9  8  7  6  5  4  3  2  1
  res.insert(residual("Cu_L",          e1_L,        2, 0x783FC00));
  res.insert(residual("Cu6Y_L",        lp0_L,       2, 0x7C3FC03));
  res.insert(residual("Cu4Y_L",        lp1_L,       2, 0x7A3FC0C));
  res.insert(residual("Cu7Y2_L",       lp2_L,       2, 0x793FC30));
  res.insert(residual("Cu2Y_L",        lp3_L,       2, 0x78BFCC0));
  res.insert(residual("CuY_L",         lp4_L,       2, 0x787FF00));
  res.insert(residual("Y_L",           e2_L,        2, 0x783FC00));
  res.insert(residual("Cu_L_Cu6Y",     e1_L_lp0,    1, 0x7C3FC03));
  res.insert(residual("Cu_L_Cu6Y_x",   e1_L_lp0_x,  1, 0x7C3FC03));
  res.insert(residual("Cu6Y_L_Cu4Y",   lp0_L_lp1,   1, 0x7E3FC0F));
  res.insert(residual("Cu6Y_L_Cu4Y_x", lp0_L_lp1_x, 1, 0x7E3FC0F));
  res.insert(residual("Cu4Y_L_Cu7Y2",  lp1_L_lp2,   1, 0x7B3FC3C));
  res.insert(residual("Cu7Y2_L_Cu2Y",  lp2_L_lp3,   1, 0x79BFCF0));
  res.insert(residual("Cu7Y2_L_Cu2Y_x",lp2_L_lp3_x, 1, 0x79BFCF0));
  res.insert(residual("Cu2Y_L_CuY",    lp3_L_lp4,   1, 0x78FFFC0));
  res.insert(residual("Cu2Y_L_CuY_x",  lp3_L_lp4_x, 1, 0x78FFFC0));
  res.insert(residual("CuY_L_Y",       lp4_L_e2,    1, 0x787FF00));
  res.insert(residual("CuY_L_Y_x",     lp4_L_e2_x,  1, 0x787FF00));
  res.insert(residual("H",             H,           3, 0x7803C00));
  res.insert(residual("Hsl",           Hsl,         3, 0x7803C00));
  res.insert(residual("HTs",           HTs,         3, 0x07C0000));
  res.insert(residual("HTl",           HTl,         3, 0x7FC3D55));
  res.insert(residual("HCuY",          Hlp4,        2, 0x0040100));
}

char* EQ = "\
There are equations as follows \n\
\t\t    eq_ID                    equation \n\n\
\t\tCu_L                T = T(x, teta_vec) + error\n\
\t\tCu6Y_L              T = T(x, teta_vec) + error\n\
\t\tCu4Y_L              T = T(x, teta_vec) + error\n\
\t\tCu7Y2_L             T = T(x, teta_vec) + error\n\
\t\tCu2Y_L              T = T(x, teta_vec) + error\n\
\t\tCuY_L               T = T(x, teta_vec) + error\n\
\t\tY_L                 T = T(x, teta_vec) + error\n\
\t\tCu_L_Cu6Y           T = T(teta_vec) + error\n\
\t\tCu_L_Cu6Y_x         x = x(teta_vec) + error\n\
\t\tCu6Y_L_Cu4Y         T = T(teta_vec) + error\n\
\t\tCu6Y_L_Cu4Y_x       x = x(teta_vec) + error\n\
\t\tCu4Y_L_Cu7Y2        T = T(teta_vec) + error\n\
\t\tCu7Y2_L_Cu2Y        T = T(teta_vec) + error\n\
\t\tCu7Y2_L_Cu2Y_x      x = x(teta_vec) + error\n\
\t\tCu2Y_L_CuY          T = T(teta_vec) + error\n\
\t\tCu2Y_L_CuY_x        x = x(teta_vec) + error\n\
\t\tCuY_L_Cu            T = T(teta_vec) + error\n\
\t\tCuY_L_Cu_x          x = x(teta_vec) + error\n\
\t\tH                   H = H(x, T, teta_vec) + error\n\
\t\tHsl                 Hsl = Hsl(x, T, teta_vec) + error\n\
\t\tHTs                 HTs - H298 = H(T, x, teta_vec) + error\n\
\t\tHTl                 HTl - H298 = H(T, x, teta_vec) + error\n\
\t\tHCuy                H = Ae + error\n\
availale in this release \n\
and twelve unknown parameters \n\
\t\tDelH  DelS   DelCp\n\
\t\t Aa    Ba     Ca   6/7 Cu + 1/7 Y = 1/7 Cu6Y \n\
\t\t Ab    Bb     Cb   4/5 Cu + 1/5 Y = 1/5 Cu4Y \n\
\t\t Ac    Bd     Cd   7/9 Cu + 2/9 Y = 1/9 Cu7Y2 \n\
\t\t Ad    Bd     Cd   2/3 Cu + 1/3 Y = 1/3 Cu2Y \n\
\t\t Ae    Be     Ce   1/2 Cu + 1/2 Y = 1/2 CuY \n\
\t\t A0    B0     C0   x(1-x)          terms for CuY(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 = "Ac";
  par[4].atr = true;
  par[5].name = "Bc";
  par[5].atr = false;
  par[6].name = "Ad";
  par[6].atr = true;
  par[7].name = "Bd";
  par[7].atr = false;
  par[8].name = "Ae";
  par[8].atr = true;
  par[9].name = "Be";
  par[9].atr = false;
  par[10].name = "A0";
  par[10].atr = false;
  par[11].name = "A1";
  par[11].atr = false;
  par[12].name = "A2";
  par[12].atr = false;
  par[13].name = "A3";
  par[13].atr = false;
  par[14].name = "B0";
  par[14].atr = false;
  par[15].name = "B1";
  par[15].atr = false;
  par[16].name = "B2";
  par[16].atr = false;
  par[17].name = "B3";
  par[17].atr = false;
  par[18].name = "Ca";
  par[18].atr = false;
  par[19].name = "Cb";
  par[19].atr = false;
  par[20].name = "Cc";
  par[20].atr = false;
  par[21].name = "Cd";
  par[21].atr = false;
  par[22].name = "Ce";
  par[22].atr = false;
  par[23].name = "C0";
  par[23].atr = false;
  par[24].name = "C1";
  par[24].atr = false;
  par[25].name = "C2";
  par[25].atr = false;
  par[26].name = "C3";
  par[26].atr = false;
  log_name = file_name + "log";
  if (strlen(file_name.c_str()) < 9)
  {
    file_name[strlen(file_name.c_str()) - 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 = 800.;
      x = lp[0].xs;
      break;
    case 4:
      T = non[1].T;
      x = lp[0].xs;
      break;
    case 5:
      T = non[1].T;
      x = non[1].x;
      break;
    case 6:
      T = non[1].T;
      x = lp[1].xs;
      break;
    case 7:
      T = 800.;
      x = lp[1].xs;
      break;
    case 8:
      T = lp[1].Tm;
      x = lp[1].xs;
      break;
    case 9:
      T = non[2].T;
      x = lp[1].xs;
      break;
    case 10:
      T = non[2].T;
      x = non[2].x;
      break;
    case 11:
      T = non[2].T;
      x = lp[2].xs;
      break;
    case 12:
      T = 800.;
      x = lp[2].xs;
      break;
    case 13:
      T = non[3].T;
      x = lp[2].xs;
      break;
    case 14:
      T = non[3].T;
      x = lp[3].xs;
      break;
    case 15:
      T = 800.;
      x = lp[3].xs;
      break;
    case 16:
      T = lp[3].Tm;
      x = lp[3].xs;
      break;
    case 17:
      T = non[4].T;
      x = lp[3].xs;
      break;
    case 18:
      T = non[4].T;
      x = lp[4].xs;
      break;
    case 19:
      T = 800.;
      x = lp[4].xs;
      break;
    case 20:
      T = lp[4].Tm;
      x = lp[4].xs;
      break;
    case 21:
      T = non[5].T;
      x = lp[4].xs;
      break;
    case 22:
      T = non[5].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 = 3000.;
  double T3 = 400.;
  ofpd << setw(wdth) << txt
     << setw(wdth) << T
     << setw(wdth) << x
     << setw(wdth) << DmixG(T1, x)
     << setw(wdth) << DmixG(T2, x)
     << setw(wdth) << DmixG(T3, x);
  print_non();
}

void post_analysis_H()
{
  ofstream of((file_name + "ax3").c_str());
  double T1 = 1373.;
  double Hm1 = h2l(T1) - h2s(T1);
  double T2 = 1410.;
  double Hm2 = h2l(T2) - h2s(T2);
  double T3 = 1963.;
  double x;
  ::of << "T  " << setw(wdth) << T1
     << "  Hm  " << setw(wdth) << Hm1 << endl;
  ::of << "T  " << setw(wdth) << T2
     << "  Hm  " << setw(wdth) << Hm2 << endl;
  of << "   x        H1373       H1410        H1963    "
     << "   Hn1373       Hn1410        Hn1963"
     << endl;
  for (x = 0.001; x < 1.01; x+= 0.05)
  {
    if (x > 1.) x = 0.999;
    of << setw(wdth) << x
       << setw(wdth) << DmixH(T1, x)
       << setw(wdth) << DmixH(T2, x)
       << setw(wdth) << DmixH(T3, x);
    of << setw(wdth) << DmixH(T1, x)/x/(1.-x)
       << setw(wdth) << DmixH(T2, x)/x/(1.-x)
       << setw(wdth) << DmixH(T3, x)/x/(1.-x);
    of << endl;
  }
}

void post_analysis_HT()
{
  ofstream of((file_name + "ax4").c_str());
  double reg[3];
  reg[0] = 0.;
  double& Temp = reg[1];
  double& xc = reg[2];
  double reg1[2];
  double& T1 = reg1[0];
  double& x1 = reg1[1];
  double Tm[8][2] = {1155., 1155.,
                     1240., 1240.,
                     1127., 1127.,
                     1170., 1170.,
                     1110., 1165.,
                     1110., 1165.,
                     1220., 1220.,
                     1070., 1070.};
  double x[8] = {0.093,
                 0.2,
                 0.278,
                 0.333,
                 0.341,
                 0.434,
                 0.5,
                 0.67};

  Tm[0][0] = non[0].T;
  T1 = Tm[0][1];
  x1 = x[0];
  if (non[0].x <= x1)
    subst1 = &lp[0];
  else
    subst1 = &e1;
  Tm[0][1] = T1 - T_solve(reg1);
  Tm[1][0] = Tm[1][1] = lp[1].Tm;
  Tm[2][0] = non[3].T;
  T1 = Tm[2][1];
  x1 = x[2];
  if (non[3].x <= x1)
    subst1 = &lp[3];
  else
    subst1 = &lp[2];
  Tm[2][1] = T1 - T_solve(reg1);
  Tm[3][0] = Tm[3][1] = lp[3].Tm;
  Tm[4][0] = non[4].T;
  T1 = Tm[4][1];
  x1 = x[4];
  if (non[4].x <= x1)
    subst1 = &lp[4];
  else
    subst1 = &lp[3];
  Tm[4][1] = T1 - T_solve(reg1);
  Tm[5][0] = non[4].T;
  T1 = Tm[5][1];
  x1 = x[5];
  if (non[4].x <= x1)
    subst1 = &lp[4];
  else
    subst1 = &lp[3];
  Tm[5][1] = T1 - T_solve(reg1);
  Tm[6][0] = Tm[6][1] = lp[4].Tm;
  Tm[7][0] = non[5].T;
  T1 = Tm[7][1];
  x1 = x[7];
  if (non[5].x <= x1)
    subst1 = &e2;
  else
    subst1 = &lp[4];
  Tm[7][1] = T1 - T_solve(reg1);

  of << "     T1          HT1         T2          HT2         T3     \
       HT3         T4          HT4         T5          HT5         T6      \
     HT6         T7          HT7         T8          HT8" << endl;
  for (double T = 950.; T < 1400.1; T +=25.)
  {
    for (int i = 0; i < 8; i++)
    {
      xc = x[i];
      if (T + 25. <= Tm[i][0])
      {
        Temp = T;
        of << setw(wdth) << Temp
           << setw(wdth) << -HTs(reg);
      }
      else if (T <= Tm[i][0])
      {
        Temp = Tm[i][0];
        of << setw(wdth) << Temp
           << setw(wdth) << -HTs(reg);
      }
      else if (T <= Tm[i][1])
      {
        of << setw(wdth) << "mis"
           << setw(wdth) << "mis";
      }
      else if (T - 25. > Tm[i][1])
      {
        Temp = T;
        of << setw(wdth) << Temp
           << setw(wdth) << -HTl(reg);
      }
      else if (T >= Tm[i][1])
      {
        Temp = Tm[i][1];
        of << setw(wdth) << Temp
           << setw(wdth) << -HTl(reg);
      }
    }
    of << endl;
  }
}

void post_analysis1()
{
  ofpd.open((file_name + "axm").c_str());
  cString txt;
  double reg[2];
  double& T = reg[0];
  double& x = reg[1];
  iii = 1;
  T = 975. + 273;
  x = lp[1].xs;
  subst1 = &lp[1];
  lp[1].Tm = T - T_solve(reg);
  T = 935. + 273;
  x = lp[3].xs;
  subst1 = &lp[3];
  lp[3].Tm = T - T_solve(reg);
  T = 935. + 273;
  x = lp[4].xs;
  subst1 = &lp[4];
  lp[4].Tm = T - T_solve(reg);

  T = 1000.;
  if (!non[0].est) e1_L_lp0(reg);
  if (!non[1].est) lp0_L_lp1(reg);
  if (!non[2].est) lp1_L_lp2(reg);
  if (!non[3].est) lp2_L_lp3(reg);
  if (!non[4].est) lp3_L_lp4(reg);
  if (!non[5].est) lp4_L_e2(reg);
  ofpd << "equilibria        T          x        "
     << "      DmixG1       DmixG2        DmixG3         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;
    if (i == -1)
      subst1 = &e1;
    else if (i == nlp)
      subst1 = &e2;
    else
      subst1 = &lp[i];
    txt = cString(subst1->name) + cString("_L");
    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)
    {
      T -= T_solve(reg);
      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);
//      if (!non[i+1].conv)
//        of << "WARNING - there was no solution "
//           << setw(wdth) << non[i+1].T1
//           << setw(wdth) << non[i+1].T2  << endl;
      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) &&
//      (par[4].use || par[5].use) && (par[6].use || par[7].use) &&
//      (par[8].use || par[9].use))

//  if (par[0].use && par[2].use && par[4].use && par[6].use && par[8].use)

void post_analysis()
{
  int i;
  file_name[strlen(file_name.c_str()) - 2] = '0' + irun++;
  of << endl << "******     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.c_str(), 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 <= 4; ++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("cuy_x.dat");
  read_errors("cuy_e.dat");
  deviates();
  eb_vs_ea();
  post_analysis_H();
  post_analysis_HT();
}

/*   test
int main()
{
  double i;
  cout.precision(4);
  cout.setf(ios::showpoint);
  cout << "T    Cul    Cus     dif" << endl;
  for (i = 1350.; i <= 1360.; i += 1.)
  {
    cout << setw(19) << i
         << setw(19) << h1l(i)
         << setw(19) << h1s(i)
         << setw(19) << h1l(i) - h1s(i) << endl;
  }
  cout << "T    Yl    Ys     dif" << endl;
  for (i = 1375.; i <= 1380.; i += 1.)
  {
    cout << setw(19) << i
         << setw(19) << h2l(i)
         << setw(19) << h2s(i)
         << setw(19) << h2l(i) - h2s(i) << endl;
  }
  for (i = 1497.; i <= 1503.; i += 1.)
  {
    cout << setw(19) << i
         << setw(19) << h2l(i)
         << setw(19) << h2s(i)
         << setw(19) << h2l(i) - h2s(i) << endl;
  }
  for (i = 1750.; i <= 1755.; i += 1.)
  {
    cout << setw(19) << i
         << setw(19) << h2l(i)
         << setw(19) << h2s(i)
         << setw(19) << h2l(i) - h2s(i) << endl;
  }
  for (i = 1794.; i <= 1805.; i += 1.)
  {
    cout << setw(19) << i
         << setw(19) << h2l(i)
         << setw(19) << h2s(i)
         << setw(19) << h2l(i) - h2s(i) << endl;
  }
  cout << "T    Cul    Cus     dif" << endl;
  for (i = 1350.; i <= 1360.; i += 1.)
  {
    cout << setw(19) << i
         << setw(19) << g1l(i)
         << setw(19) << g1s(i)
         << setw(19) << g1l(i) - g1s(i) << endl;
  }
  cout << "T    Yl    Ys     dif" << endl;
  for (i = 1750.; i <= 1800.; i += 1.)
  {
    cout << setw(19) << i
         << setw(19) << g2l(i)
         << setw(19) << g2s(i)
         << setw(19) << g2l(i) - g2s(i) << endl;
  }
  i = 1357.77;
  cout << setw(19) << i
       << setw(19) << g2l(i)
       << setw(19) << g2s(i)
       << setw(19) << g2l(i) - g2s(i) << endl;
  A0 = 100000.;
  A1 = 100.;
  A2 = .1;
  A3 = 1e-4;
  B0 = 100.;
  B1 = .1;
  B2 = 1e-4;
  B3 = 1e-7;
  C0 = .01;
  C1 = 1e-3;
  C2 = 1e-5;
  C3 = 1e-7;
  double T = 1500.;
  double x;
  cout << "DmixG    DmixH-T*DmixS   (1-x)Dmixg1+xDmixg2" << endl;
  for (x = 0.1; x <= 0.9; x += 0.1)
  {
    cout << setw(19) << DmixG(T, x)
         << setw(19) << DmixH(T, x) - T*DmixS(T, x)
         << setw(19) << (1.-x)*Dmixg1(T, x) + x*Dmixg2(T, x) << endl;
  }
  cout << "Dmixg1    Dmixh1-T*Dmixs1   Dmixg2    Dmixh2-T*Dmixs2" << endl;
  for (x = 0.1; x <= 0.9; x += 0.1)
  {
    cout << setw(19) << Dmixg1(T, x)
         << setw(19) << Dmixh1(T, x) - T*Dmixs1(T, x);
    cout << setw(19) << Dmixg2(T, x)
         << setw(19) << Dmixh2(T, x) - T*Dmixs2(T, x);
    cout << endl;
  }
  return 0;
}

 end of test */