/*
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 <math.h>

extern double hBaY;
extern double sBaY;

double R = 8.31441;
double A0BAY;
double A1BAY;
double A0BACUY;
double A0BACU;
double A1BACU;
double A2BACU;
double A3BACU;
double A0CUY;
double A1CUY;
double A2CUY;
double A3CUY;
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 Dmixg1(double T, double XCU, double XY)
{
  if (XCU <= 0. || XY <= 0. || (1. - XCU - XY) <= 0.)
    return 0.;
  // Cu
  double ANS1=-3.*A2CUY*pow(XCU,3)*XY+6.*A2CUY*pow(XCU,2)*pow(XY,2)
     +3.*A2CUY*pow(XCU,2)*XY-3.*A2CUY*XCU*pow(XY,3)-4.*A2CUY*XCU*pow(XY,2)
     +A2CUY*pow(XY,3)+32.*A3BACU*pow(XCU,5)+80.*A3BACU*pow(XCU,4)*XY
     -100.*A3BACU*pow(XCU,4)+72.*A3BACU*pow(XCU,3)*pow(XY,2)
     -188.*A3BACU*pow(XCU,3)*XY+116.*A3BACU*pow(XCU,3)
     +28.*A3BACU*pow(XCU,2)*pow(XY,3)-117.*A3BACU*pow(XCU,2)*pow(XY,2)
     +150.*A3BACU*pow(XCU,2)*XY-61.*A3BACU*pow(XCU,2)
     +4.*A3BACU*XCU*pow(XY,4)-26.*A3BACU*XCU*pow(XY,3)
     +54.*A3BACU*XCU*pow(XY,2)-46.*A3BACU*XCU*XY+14.*A3BACU*XCU
     -A3BACU*pow(XY,4)+4.*A3BACU*pow(XY,3)-6.*A3BACU*pow(XY,2)
     +4.*A3BACU*XY-A3BACU+4.*A3CUY*pow(XCU,4)*XY
     -12.*A3CUY*pow(XCU,3)*pow(XY,2)-4.*A3CUY*pow(XCU,3)*XY
	  +12.*A3CUY*pow(XCU,2)*pow(XY,3)+9.*A3CUY*pow(XCU,2)*pow(XY,2)
     -4.*A3CUY*XCU*pow(XY,4)-6.*A3CUY*XCU*pow(XY,3)+A3CUY*pow(XY,4);

  return log(XCU)*R*T+2.*A0BACUY*pow(XCU,2)*XY+2.*A0BACUY*XCU*pow(XY,2)
    -3.*A0BACUY*XCU*XY-A0BACUY*pow(XY,2)+A0BACUY*XY+A0BACU*pow(XCU,2)
    +A0BACU*XCU*XY-2.*A0BACU*XCU-A0BACU*XY+A0BACU+A0BAY*XCU*XY
    +A0BAY*pow(XY,2)-A0BAY*XY-A0CUY*XCU*XY+A0CUY*XY+4.*A1BACU*pow(XCU,3)
    +6.*A1BACU*pow(XCU,2)*XY-9.*A1BACU*pow(XCU,2)+2.*A1BACU*XCU*pow(XY,2)
    -8.*A1BACU*XCU*XY+6.*A1BACU*XCU-A1BACU*pow(XY,2)+2.*A1BACU*XY-A1BACU
    +2.*A1BAY*pow(XCU,2)*XY+6.*A1BAY*XCU*pow(XY,2)-4.*A1BAY*XCU*XY
    +4.*A1BAY*pow(XY,3)-6.*A1BAY*pow(XY,2)+2.*A1BAY*XY
    +2.*A1CUY*pow(XCU,2)*XY-2.*A1CUY*XCU*pow(XY,2)-2.*A1CUY*XCU*XY
    +A1CUY*pow(XY,2)+12.*A2BACU*pow(XCU,4)+24.*A2BACU*pow(XCU,3)*XY
    -32.*A2BACU*pow(XCU,3)+15.*A2BACU*pow(XCU,2)*pow(XY,2)
    -44.*A2BACU*pow(XCU,2)*XY+29.*A2BACU*pow(XCU,2)
	 +3.*A2BACU*XCU*pow(XY,3)-16.*A2BACU*XCU*pow(XY,2)+23.*A2BACU*XCU*XY
    -10.*A2BACU*XCU-A2BACU*pow(XY,3)+3.*A2BACU*pow(XY,2)
    -3.*A2BACU*XY+A2BACU+ANS1;
}

double Dmixg2(double T, double XCU, double XY)
{
  if (XCU <= 0. || XY <= 0. || (1. - XCU - XY) <= 0.)
    return 0.;
  // Y
  double ANS1=-4.*A2CUY*pow(XCU,2)*XY-3.*A2CUY*XCU*pow(XY,3)+3.*A2CUY*XCU*pow(XY,2)
     +32.*A3BACU*pow(XCU,5)+80.*A3BACU*pow(XCU,4)*XY-80.*A3BACU*pow(XCU,4)
     +72.*A3BACU*pow(XCU,3)*pow(XY,2)-144.*A3BACU*pow(XCU,3)*XY
     +72.*A3BACU*pow(XCU,3)+28.*A3BACU*pow(XCU,2)*pow(XY,3)
     -84.*A3BACU*pow(XCU,2)*pow(XY,2)+84.*A3BACU*pow(XCU,2)*XY
	  -28.*A3BACU*pow(XCU,2)+4.*A3BACU*XCU*pow(XY,4)
     -16.*A3BACU*XCU*pow(XY,3)+24.*A3BACU*XCU*pow(XY,2)
     -16.*A3BACU*XCU*XY+4.*A3BACU*XCU+4.*A3CUY*pow(XCU,4)*XY-A3CUY*pow(XCU,4)
     -12.*A3CUY*pow(XCU,3)*pow(XY,2)+6.*A3CUY*pow(XCU,3)*XY
     +12.*A3CUY*pow(XCU,2)*pow(XY,3)-9.*A3CUY*pow(XCU,2)*pow(XY,2)
     -4.*A3CUY*XCU*pow(XY,4)+4.*A3CUY*XCU*pow(XY,3);

  return log(XY)*R*T+2.*A0BACUY*pow(XCU,2)*XY-A0BACUY*pow(XCU,2)
    +2.*A0BACUY*XCU*pow(XY,2)-3.*A0BACUY*XCU*XY+A0BACUY*XCU
    +A0BACU*pow(XCU,2)+A0BACU*XCU*XY-A0BACU*XCU+A0BAY*XCU*XY-A0BAY*XCU
    +A0BAY*pow(XY,2)-2.*A0BAY*XY+A0BAY-A0CUY*XCU*XY+A0CUY*XCU
    +4.*A1BACU*pow(XCU,3)+6.*A1BACU*pow(XCU,2)*XY-6.*A1BACU*pow(XCU,2)
    +2.*A1BACU*XCU*pow(XY,2)-4.*A1BACU*XCU*XY+2.*A1BACU*XCU
    +2.*A1BAY*pow(XCU,2)*XY-A1BAY*pow(XCU,2)+6.*A1BAY*XCU*pow(XY,2)
    -8.*A1BAY*XCU*XY+2.*A1BAY*XCU
	 +4.*A1BAY*pow(XY,3)-9.*A1BAY*pow(XY,2)+6.*A1BAY*XY-A1BAY
    +2.*A1CUY*pow(XCU,2)*XY-A1CUY*pow(XCU,2)-2.*A1CUY*XCU*pow(XY,2)
    +2.*A1CUY*XCU*XY+12.*A2BACU*pow(XCU,4)+24.*A2BACU*pow(XCU,3)*XY
    -24.*A2BACU*pow(XCU,3)+15.*A2BACU*pow(XCU,2)*pow(XY,2)
    -30.*A2BACU*pow(XCU,2)*XY+15.*A2BACU*pow(XCU,2)
    +3.*A2BACU*XCU*pow(XY,3)-9.*A2BACU*XCU*pow(XY,2)+9.*A2BACU*XCU*XY
    -3.*A2BACU*XCU-3.*A2CUY*pow(XCU,3)*XY+A2CUY*pow(XCU,3)
    +6.*A2CUY*pow(XCU,2)*pow(XY,2)+ANS1;
}

double Dmixg3(double T, double XCU, double XY)
{
  if (XCU <= 0. || XY <= 0. || (1. - XCU - XY) <= 0.)
    return 0.;
  // Ba
  double ANS1=42.*A3BACU*pow(XCU,2)*XY-7.*A3BACU*pow(XCU,2)
     +4.*A3BACU*XCU*pow(XY,4)-12.*A3BACU*XCU*pow(XY,3)
     +12.*A3BACU*XCU*pow(XY,2)-4.*A3BACU*XCU*XY+4.*A3CUY*pow(XCU,4)*XY
     -12.*A3CUY*pow(XCU,3)*pow(XY,2)+12.*A3CUY*pow(XCU,2)*pow(XY,3)
     -4.*A3CUY*XCU*pow(XY,4);

  return log(-XCU-XY+1.)*R*T+2.*A0BACUY*pow(XCU,2)*XY
    +2.*A0BACUY*XCU*pow(XY,2)-A0BACUY*XCU*XY+A0BACU*pow(XCU,2)
    +A0BACU*XCU*XY+A0BAY*XCU*XY+A0BAY*pow(XY,2)-A0CUY*XCU*XY
    +4.*A1BACU*pow(XCU,3)+6.*A1BACU*pow(XCU,2)*XY-3.*A1BACU*pow(XCU,2)
    +2.*A1BACU*XCU*pow(XY,2)-2.*A1BACU*XCU*XY+2.*A1BAY*pow(XCU,2)*XY
    +6.*A1BAY*XCU*pow(XY,2)-2.*A1BAY*XCU*XY+4.*A1BAY*pow(XY,3)
    -3.*A1BAY*pow(XY,2)+2.*A1CUY*pow(XCU,2)*XY-2.*A1CUY*XCU*pow(XY,2)
    +12.*A2BACU*pow(XCU,4)+24.*A2BACU*pow(XCU,3)*XY
    -16.*A2BACU*pow(XCU,3)+15.*A2BACU*pow(XCU,2)*pow(XY,2)
	 -20.*A2BACU*pow(XCU,2)*XY+5.*A2BACU*pow(XCU,2)
    +3.*A2BACU*XCU*pow(XY,3)-6.*A2BACU*XCU*pow(XY,2)+3.*A2BACU*XCU*XY
    -3.*A2CUY*pow(XCU,3)*XY+6.*A2CUY*pow(XCU,2)*pow(XY,2)
    -3.*A2CUY*XCU*pow(XY,3)+32.*A3BACU*pow(XCU,5)+80.*A3BACU*pow(XCU,4)*XY
    -60.*A3BACU*pow(XCU,4)+72.*A3BACU*pow(XCU,3)*pow(XY,2)
    -108.*A3BACU*pow(XCU,3)*XY+36.*A3BACU*pow(XCU,3)
    +28.*A3BACU*pow(XCU,2)*pow(XY,3)-63.*A3BACU*pow(XCU,2)*pow(XY,2)+ANS1;
}

void coeff(double T)
{
  A0BAY = hBaY-sBaY*T;
//  A0BAY = 8.596e+04;
//  A0BAY = 200700-22.3*T;
  A1BAY = 0.;
  A0BACUY = 0.;
  A0BACU = -7191. + T*3.363;
  A1BACU = 0.;
  A2BACU = 0.;
  A3BACU = 0.;
  A0CUY = A0 + T*B0 + C0*T*log(T);
  A1CUY = A1 + T*B1 + C1*T*log(T);
  A2CUY = 0.;
  A3CUY = 0.;
}

double gCus(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 gCul(double T, double xCu = 0., double xY = 0.)
{
// 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
  if (xCu && xY)
  {
    coeff(T);
    g += Dmixg1(T, xCu, xY);
  }
  return g;
}

double gYs(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 gYl(double T, double xCu = 0., double xY = 0.)
{
// 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
  if (xCu && xY)
  {
    coeff(T);
	 g += Dmixg2(T, xCu, xY);
  }
  return g;
}

double gBas(double T)
{
// Ba - 91Din (based on TPIS)
  return (T < 1000.) ?
       -17685.226+233.78606*T-42.889*T*log(T)-1.8314e-3*pow(T,2)
       -0.000095e-6*pow(T,3)+705880/T:    //  298.15,1000
       -64873.614+608.188389*T-94.2824199*T*log(T)+19.504772e-3*pow(T,2)
       -1.051353e-6*pow(T,3)+8220192/T;   //  1000,2995
}

double gBal(double T, double xCu = 0., double xY = 0.)
{
// Ba - 91Din (based on TPIS)
  double g;
  if (T < 1000.)
    g = -9738.988+229.540143*T-43.4961089*T*log(T)-2.346416e-3*pow(T,2)
        +0.991223e-6*pow(T,3)+723016/T;   //  298.15,1000
  else
    g = -7381.093+235.49642*T-45.103*T*log(T)+2.154e-3*pow(T,2)
        +0.000027e-6*pow(T,3)-365/T;      //  1000,2995
  if (xCu && xY)
  {
    coeff(T);
    g += Dmixg3(T, xCu, xY);
  }
  return g;
}

double gCu6Y(double T)
{
  return -1.261e+04*7 + gCus(T)*6 + gYs(T);
}

double gCu4Y(double T)
{
  return -1.656e+04*5 + gCus(T)*4 + gYs(T);
}

double gCu7Y2(double T)
{
  return -1.728e+04*9 + gCus(T)*7 + gYs(T)*2;
}

double gCu2Y(double T)
{
  return -2.01e+04*3 + gCus(T)*2 + gYs(T);
}

double gCuY(double T)
{
  return -2.018e+04*2 + gCus(T) + gYs(T);
}

double gBaCu(double T)
{
  return -2594*2 + gBas(T) + gCus(T);
}

double gBaCu13(double T)
{
  return -924.6*14 + gBas(T) + gCus(T)*13;
}

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);
}

void dxdz(double& dx1dz1, double& dx1dz2,
          double& dx2dz1, double& dx2dz2,
          double z1, double z2)
{
  long double x11 = 1./(1. + exp((long double)-z1));
  long double x21 = 1./(1. + exp((long double)-z2));
  long double dx11dz1 = exp((long double)-z1)*x11*x11;
  long double dx21dz2 = exp((long double)-z2)*x21*x21;
  long double a, dadz1, dadz2;
  if (x21 >= x11)
  {
    a = x11/x21;
    dadz1 = dx11dz1/x21;
    dadz2 = - x11/x21/x21*dx21dz2;
  }
  else
  {
    a = x21/x11;
    dadz1 = - x21/x11/x11*dx11dz1;
    dadz2 = dx21dz2/x11;
  }
//  x1 = x11/(1. + a);
  dx1dz1 = dx11dz1/(1. + a) - x11*dadz1/(1. + a)/(1. + a);
  dx1dz2 = - x11*dadz2/(1. + a)/(1. + a);
//  x2 = x21/(1. + a);
  dx2dz1 = - x21*dadz1/(1. + a)/(1. + a);
  dx2dz2 = dx21dz2/(1. + a) - x21*dadz2/(1. + a)/(1. + a);
}