rkg_model.c

Go to the documentation of this file.

/* src/rkg_model.c
 *
 * Copyright (C) 2011-2018 Dongdong Li
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundataion; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABLITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

//Redlich-Kister-Guggenheim formula
static double rkg_ai (SOLIDSOLUTION_PHASE ss, int index, double *n, double T, double P)
{
  double lnri = 0.0;
  double A0[ss.species_num][ss.species_num];
  double A1[ss.species_num][ss.species_num];
  double x[ss.species_num];
  int i, j, k, l;
  double sumn;

  for (i = 0; i < ss.species_num; i ++)
  {
    for (j = 0; j < ss.species_num; j ++)
    {
      A0[i][j] = 0.0;
      A1[i][j] = 0.0;
    }
  }      

  sumn = 0.0;
  for (k = 0; k < ss.species_num; k ++)
  {
    sumn += n[k];
  }

  for (k = 0; k < ss.species_num; k ++)
  {
    x[k] = n[k] / sumn;
  }

  for (l = 0; l < ss.param_num; l ++)
  {
    for (k = 0; k < ss.species_num; k ++)
    {
      if (strcmp (ss.params[l].species_a, ss.solidsolution_species[k].species_symbol) == 0)
      {
        i = k;
      }
    }

    for (k = 0; k < ss.species_num; k ++)
    {
      if (strcmp (ss.params[l].species_b, ss.solidsolution_species[k].species_symbol) == 0)
      {
        j = k;
      }
    }

    A0[i][j] = ss.params[l].A0;
    A0[j][i] = A0[i][j];
    A1[i][j] = ss.params[l].A1;
    A1[j][i] = A1[i][j];
  }

  for (j = 0; j < ss.species_num; j++)
  {
    lnri += ((1 - x[index]) * x[j] * A0[index][j]/(R * T) + 
                  x[index] * x[j] * (1.0 - 2.0 * x[j]) * A1[index][j]/(R * T));
  }

    if (sumn == 0 || x[index] == 0 
                || x[index] * exp(lnri) < 1.0e-15 
                || x[index] * exp(lnri) > 1.0e+15) {return 1.0E-10;}

  return x[index] * exp(lnri);
}

extern double rkg_a (SOLIDSOLUTION_PHASE ss, int index, double *n, double T, double P)
{
  rkg_ai (ss, index, n, T, P);
}