//###
//*****************************************************************************
//                                                                            *
//     Function Name: matrix                                                  *
//                                                                            *
//     File Name: matrix.cpp                                                  *
//                                                                            *
//     Purpose:  This routine computes the elements of a symetrical matrix    *
//               and stores them in their proper location.  Only the upper    *
//               diagonal matix is set up in this routine and reflected to the*
//               lower part in subroutine gauss.                              *
//                                                                            *
//     Author:    Robert L. Reese                                             *
//                                                                            *
//     Date: June 01, 1973                                                    *
//                                                                            *
//     Language: Borland C++ Compiler, Version 4.5                            *
//                                                                            *
//     Hardware: IBM Compatible PC                                            *
//                                                                            *
//     Operating System: Microsoft Windows 98                                 *
//                                                                            *
//     Update History:                                                        *
//                                                                            *
//     Name        Date      Revision          Changes                        *
//     ----        ----      ---------         -------                        *
//                                                                            *
//    R.L. Reese   06/01/73      0             Initial Release                *
//    R.L. Reese   11/23/03      A             Converted to C++               *
//                                                                            *
//----------------------------------------------------------------------------*
//                                                                            *
//  Usage:                                                                    *
//                                                                            *
//  matrix(output_unit) where:                                                *
//                                                                            *
//     output_unit - FILE pointer to logical output listing unit              *
//                                                                            *
//                                                                            *
//----------------------------------------------------------------------------*
//                                                                            *
//     Definitions:                                                           *
//                                                                            *
//     Local Variables:                                                       *
//     ----------------                                                       *
//                                                                            *
//     f         - Double Precision variable used as Gibbs free energy for    *
//                 gaseous species                                            *
//     i         - Integer variable used as do loop index                     *
//     iq2       - Integer variable used as number of columns in matrix       *
//     j         - Integer variable used as do loop index                     *
//     k         - Integer variable used as do loop index                     *
//     kk        - Integer variable used as number of condensed species being *
//                 considered in the present matrix                           *
//     kmat      - Integer variable used to store number of columns in matrix *
//     loop_end  - Logical variable used to go to end of loop                 *
//     term      - Real variable used to sum mole fraction times stoiciometric*
//                 coefficients of gaseous species                            *
//                                                                            *
//     Pointers                                                               *
//     -----------------                                                      *
//                                                                            *
//     None                                                                   *
//                                                                            *
//     Global Variables:                                                      *
//     -----------------                                                      *
//                                                                            *
//     double.h                                                               *
//     g            - Double precision array used as elements of the matrix   *
//     indx.h                                                                 *
//     imat         - Integer variable used as number of rows in matrix       *
//     iq1          - Integer variable used as number of elements plus the    *
//                    number of condensed species included in matrix          *
//                    equations plus one                                      *
//     nlm          - Integer variable used as number of elements in the      *
//                    system                                                  *
//     npt          - Integer variable used as index for data saved for output*
//     ns           - Integer variable used as number of species being        *
//                    considered in the system                                *
//     misc.h                                                                 *
//     bo           - Real array used as assigned Kg-atoms of element i per   *
//                    Kg of total reactant                                    *
//     enn          - Real variable used as moles of mixture, KG-mole/KG      *
//     gamma        - Real array used as fugacity coefficients for gases      *
//     pp           - Real variable used as assigned P for current point      *
//     sumn         - Real variable used as moles of mixture                  *
//     tm           - Real variable used as log P/n                           *
//     tt           - Real variable used as assigned T in kelvin for current  *
//                    point                                                   *
//     vol          - Real variable used as molar volume of species           *
//     speces.h                                                               *
//     a            - Real variable used as stoichiometric coefficient of ith *
//                    element in jth species                                  *
//     en           - Real array used as Kg-mole of jth species/ Kg of mixture*
//     enln         - Real array used as log of en                            *
//     ho           - Real array used as dimensionless enthalpy for species   *
//     iuse         - Integer array used to identify status of ith specie     *
//     s            - Real array used as dimensionless entropy calculated from*
//                    coefficients for reaction species                       *
//                                                                            *
//     Constants:                                                             *
//     ---------                                                              *
//                                                                            *
//     local                                                                  *
//     FALSE       - Defines value of logical false = 0                       *
//     r           - Real value of universal gas constant                     *
//     TRUE        - Defines value of logical true = 1                        *
//                                                                            *
//     Include Files:                                                         *
//     --------------                                                         *
//                                                                            *
//     arrays.h          - Contains sizes for include files                   *
//     double.h          - Contains double precision variables                *
//     indx.h            - Contains index variables                           *
//     misc.h            - Contains miscellaneous variables                   *
//     speces.h          - Contains species variables                         *
//     stdio.h           - Contains C++ I/O functions                         *
//                                                                            *
//     External Software:                                                     *
//     -----------------                                                      *
//                                                                            *
//     None                                                                   *
//                                                                            *
//     Files:                                                                 *
//     -----                                                                  *
//                                                                            *
//       Input Files:                                                         *
//                                                                            *
//       None                                                                 *
//                                                                            *
//       Intermediate Files:                                                  *
//                                                                            *
//       None                                                                 *
//                                                                            *
//       Output Files:                                                        *
//                                                                            *
//       None                                                                 *
//                                                                            *
//*****************************************************************************
//###
 
#include "arrays.h"
#include "double.h"
#include "indx.h"
#include "misc.h"
#include "speces.h"
#include <stdio.h>
 
#define FALSE 0
#define TRUE 1
 
void matrix(FILE *output_unit)
{
 
#include "constant.h"
 
//  const float r = 8314.3/4184. ;
 
  int loop_end ;
 
  double f ;
 
  float term ;
 
  int i ;
  int iq2 ;
  int j ;
  int k ;
  int kk ;
  int kmat ;
 
  iq2  = iq1 + 1 ;
  kmat = iq2 ;
  imat = iq1 ;
 
//     Clear matrix storages to zero
 
  for(i=0;i<imat;i++)
  {
    for(j=0;j<kmat;j++)
    {
      g[i][j] = 0.0 ;
    }
  }
  g[iq2-1][iq1-1] = 0. ;
 
//     Begin set up of iteration matrix
 
  kk = nlm ;
 
//start matrix_1
  for(j=0;j<ns;j++)
  {
    loop_end = FALSE ;
//  start matrix_2
    if(iuse[j] >= 0)
    {
//    start matrix_3
      if(iuse[j] <= 0)
      {
        f = (ho[j]-s[j]+enln[j] + tm + gamma[j])*en[j][npt] ;
 
//     Calculate the elements g(i,k)
 
//      start matrix_4
        for(i=0;i<nlm;i++)
        {
//        start matrix_5
          if(a[i][j] != 0.0)
          {
            term = a[i][j]*en[j][npt] ;
//          start matrix_6
            for(k=i;k<nlm;k++)
            {
              g[i][k] = g[i][k] + a[k][j]*term ;
            } // end matrix_6
 
            g[i][iq1-1] = g[i][iq1-1]+term ;
            g[i][iq2-1] = g[i][iq2-1] + a[i][j] * f ;
          } // end matrix_5
        } // end matrix_4
        g[iq1-1][iq2-1] = g[iq1-1][iq2-1] + f ;
        loop_end = TRUE ;
      } // end matrix_3
 
//     Condensed species.
 
//    start matrix_7
      if( ! loop_end)
      {
        kk = kk + 1 ;
//      start matrix_8
        for(i=0;i<nlm;i++)
        {
          g[i][kk-1] = a[i][j] ;
          g[i][kmat-1] = g[i][kmat-1] - a[i][j]*en[j][npt] ;
        } // end matrix_8
        g[kk-1][kmat-1] = ho[j]-s[j] ;
//      start matrix_9
        if(pp > 1.0)
        {
          g[kk-1][kmat-1] = ho[j]-s[j]+((vol[j]*(pp-1.0))/(r*tt)) ;
        } // end matrix_9
      } // matrix_7
    } // end matrix_2
  } // end matrix_1
  g[iq1-1][iq1-1] = sumn - enn ;
 
//     Reflect symmetric portions of the matrix
 
//start matrix_10
  for(i=0;i<iq1;i++)
  {
//  start matrix_11
    for(j=i;j<iq1;j++)
    {
      g[j][i] = g[i][j] ;
    } // end matrix_11
  } // end matrix_10
 
//     Complete the right hand side
 
//start matrix_12
  for(i=0;i<nlm;i++)
  {
    g[i][kmat-1] = g[i][kmat-1]+(bo[i]-g[i][iq1-1]) ;
  } // end matrix_12
  g[iq1-1][kmat-1] = g[iq1-1][kmat-1] + enn - sumn ;
 
//start debug
//  for(i=0;i<4;i++)
//  {
//      fprintf(output_unit,"%15.8f %15.8f %15.8f %15.8f\n",
//                            g[0][i],g[1][i],g[2][i],g[3][i]);
//  } // end debug
//  fprintf(output_unit,"\n");
 
} // end function matrix
//###
//*****************************************************************************
//                                                                            *
//     Function Name: molewt                                                  *
//                                                                            *
//     File Name: molewt.cpp                                                  *
//                                                                            *
//     Purpose:  Molewt computes new base values for the reactants.  This     *
//               routine will handle four types of problems.  It will alter   *
//               the values between the reactants without adding or removing  *
//               any reactants.  It will consider a ternary system by         *
//               considering first two, then all three components             *
//               simultaneously.  This modification involved mking changes to *
//               react to recompute the mole fractions every time a new base  *
//               value was set, and to search to either store species to be   *
//               considered later or remove species from storage.  Molewt will*
//               also consider a ternary system where all three reactants are *
//               in the system at the same time at the start.  The fourth type*
//               of problem involves starting out with n reactants and adding *
//               one reactant at a time until all the reactants have been     *
//               added.  Molewt calls react and init every time a new base    *
//               value is computed.  Search is called only when a reactant is *
//               added or removed from the system.                            *
//                                                                            *
//     Author:    Robert L. Reese                                             *
//                                                                            *
//     Date: June 01, 1973                                                    *
//                                                                            *
//     Language: Borland C++ Compiler, Version 4.5                            *
//                                                                            *
//     Hardware: IBM Compatible PC                                            *
//                                                                            *
//     Operating System: Microsoft Windows 98                                 *
//                                                                            *
//     Update History:                                                        *
//                                                                            *
//     Name        Date      Revision          Changes                        *
//     ----        ----      ---------         -------                        *
//                                                                            *
//    R.L. Reese   06/01/73      0             Initial Release                *
//    R.L. Reese   11/23/03      A             Converted to C++               *
//                                                                            *
//----------------------------------------------------------------------------*
//                                                                            *
//  Usage:                                                                    *
//                                                                            *
//  molewt(eof,                                                               *
//         output_unit,                                                       *
//         input_unit,                                                        *
//         scratch_unit,                                                      *
//         therm_unit) where:                                                 *
//                                                                            *
//     eof         - integer pointer to error flag                            *
//                  = 0, read ok                                              *
//                 != 0, eof or read error occurred                           *
//     output_unit - FILE pointer to output listing logical unit              *
//     input_unit  - FILE pointer to input file logical unit                  *
//     scratch_unit- FILE pointer to temporary output unit for thermo data    *
//     therm_unit  - FILE pointer to thermodynamic data input logical unit    *
//                                                                            *
//----------------------------------------------------------------------------*
//                                                                            *
//     Definitions:                                                           *
//                                                                            *
//     Local Variables:                                                       *
//     ----------------                                                       *
//                                                                            *
//     add_new            - Logical variable used to test for adding a new    *
//                          reactant to the system                            *
//     add_or_increment   - Logical variable used to test if a new reactant   *
//                          should be added or should the constant reactant   *
//                          be incremented                                    *
//     add_reactant       - Logical variable used to decide if to exit the add*
//                          new reactant loop                                 *
//     bypass_search      - Logical variable used to test for calling         *
//                          function search                                   *
//     i                  - Integer variable used as do loop index            *
//     reset_initial      - Logical variable used to determine if the intial  *
//                          estimates need to be reset                        *
//     store_values       - Logical variable used to determine if values need *
//                          to be stored                                      *
//     tim                - Real variable used as multiplier of tri value in  *
//                          plus or systm problems                            *
//                                                                            *
//     Pointers                                                               *
//     -----------------                                                      *
//                                                                            *
//     eof_pt       - integer pointer to eof values returned from routines    *
//                                                                            *
//     Global Variables:                                                      *
//     -----------------                                                      *
//                                                                            *
//     local                                                                  *
//     add_all      - Logical variable used to test for adding all the        *
//                    reactants to the system, or increment each              *
//                    reactant                                                *
//     add_remove   - Logical variable used to test for adding or             *
//                    removing a reactant                                     *
//     call_react   - Logical variable used to test for calling               *
//                    function react                                          *
//     loup_exit    - Logical variable used to determine if to exit loop      *
//     increment    - Logical variable used to determine if to increment      *
//                    moles between included reactants                        *
//     process_loop - Logical variable used to control the outer              *
//                    process loop                                            *
//     rtest        - Real variableused as absolute value of increment        *
//                    of ith reactant mole fraction                           *
//     zerot        - Real variable used as test for rtest test               *
//     therm.cpp                                                              *
//     aa           - Double precision array used to store increment          *
//                    value of reactant mole fraction                         *
//     bb           - Real array used to store base value of reactant         *
//                    mole fraction                                           *
//     ll           - Integer variable used to determine what reactants       *
//                    are to be added or removed                              *
//     removed      - Integer variable used to determine if the first         *
//                    reactant in a ternary system is to be removed           *
//     to           - Real array used to store final value of reactant        *
//                    mole fraction                                           *
//     indx.h                                                                 *
//     alter        - Logical variable used to indicate reactant mole         *
//                    fractions with unequal top and base values are to be    *
//                    altered                                                 *
//     newmol       - Integer variable used to determine what species are to  *
//                    be added or deleted from the current system             *
//     nmin         - Integer variable used as number of elements that are put*
//                    into temporary storage for system and plus problems     *
//     nmix         - Integer variable used to increment multiplier for mole  *
//                    fraction computations in plus and systm problems        *
//     npt          - Integer variable used as index for data saved for       *
//                    output                                                  *
//     nrat         - Integer variable used to determine if the reactants     *
//                    should be added or incremented in systm and plus        *
//                    problems                                                *
//     nreac        - Integer variable used as number of reactant records     *
//     plus         - Logical variable used to indicate if reactants are to be*
//                    added to the system or if a ternary system is being     *
//                    considered                                              *
//     systm        - Logical variable used to denote that a ternary system is*
//                    being considered                                        *
//     misc.h                                                                 *
//     add          - Real array used as increment of ith reactant mole       *
//                    fraction                                                *
//     base         - Real array containing first value of mole fraction of   *
//                    ith reactant                                            *
//     top          - Real variable containing final value of mole fraction of*
//                    ith reactant                                            *
//     tri          - Real variable containing increment value of mole        *
//                    fraction for all reactants of ternary system            *
//                                                                            *
//     Constants:                                                             *
//     ---------                                                              *
//                                                                            *
//     local                                                                  *
//     FALSE           - Defines value of logical false = 0                   *
//     TRUE            - Defines value of logical true = 1                    *
//     arrays.h                                                               *
//     element_max     - Integer value of maximum number of elements allowed  *
//     constant.h                                                             *
//     newmol_test     - Integer value to use for test to determine if new    *
//                       search is to be bypassed                             *
//     switch_index    - Integer value to use for index in switching values   *
//     third_add       - Integer value to use for third index of add array    *
//     third_base      - Integer value to use for third index of base array   *
//     third_top       - Integer value to use for third index of top array    *
//                                                                            *
//     Include Files:                                                         *
//     -------------                                                          *
//                                                                            *
//     arrays.h          - Contains sizes for include files                   *
//     constants.h       - Contains various constants                         *
//     indx.h            - Contains index variables                           *
//     iostream.h        - Contains C++ output routines                       *
//     math.h            - Contains C++ math routines                         *
//     misc,h            - Contains miscellaneous variables                   *
//     stdio.h           - C++ file containing I/O routines                   *
//     stdlib.h          - C++ file containing standard routines              *
//                                                                            *
//     External Software:                                                     *
//     -----------------                                                      *
//                                                                            *
//     increment_add      - Routine that decides if reactant is to be removed *
//                          added or incremented                              *
//     increment_only     - Routine to increment the base values only         *
//     init               - Routine to set initial values                     *
//     react              - Routine to read reactant records and compute mole *
//                          fractions                                         *
//     remove_and_add     - Routine to add or remove elements from system     *
//     search             - Routine to search file for thermo data for species*
//                          to be considered                                  *
//                                                                            *
//     Files:                                                                 *
//     -----                                                                  *
//                                                                            *
//       Input Files:                                                         *
//                                                                            *
//       None                                                                 *
//                                                                            *
//       Intermediate Files:                                                  *
//                                                                            *
//       None                                                                 *
//                                                                            *
//       Output Files:                                                        *
//                                                                            *
//       Outputs to screen                                                    *
//                                                                            *
//*****************************************************************************
//###
 
#include "arrays.h"
#include "indx.h"
#include <iostream.h>
#include <math.h>
#include "misc.h"
#include <stdio.h>
#include <stdlib.h>
 
#define FALSE 0
#define TRUE 1
 
// used by molewt and remove_and_add
 
extern double aa[element_max] ;
 
extern float bb[element_max] ;
extern float to[element_max] ;
 
extern int ll;
extern int removed;
 
int add_all;
int call_react ;
 
// used by molewt and increment_add
 
int add_remove ;
int loup_exit ;
int increment ;
int process_loop ;
 
// used by molewt and increment_only
 
float rtest ;
float zerot ;
 
void increment_add(void);
void increment_only(void);
void init();
void react(int *eof_pt,
           FILE *output_unit,
           FILE *input_unit);
void remove_and_add();
void search(int *eof_pt,
            FILE *scratch_unit,
            FILE *therm_unit,
            FILE *output_unit);
 
void molewt(int *eof_pt,
            FILE *output_unit,
            FILE *input_unit,
            FILE *scratch_unit,
            FILE *therm_unit)
{
 
#include "constant.h"
 
  int add_new ;
  int add_or_increment ;
  int add_reactant ;
  int bypass_search ;
  int i ;
  int reset_initial ;
  int store_values ;
 
  float tim ;
 
  add_all          = FALSE ;
  add_new          = FALSE ;
  add_or_increment = FALSE ;
  add_remove       = FALSE ;
  call_react       = FALSE ;
  increment        = FALSE ;
  reset_initial    = FALSE ;
  store_values     = FALSE ;
  process_loop     = TRUE ;
 
  nmin = 0 ;
 
//start outer_loop
  do
  {
 
//  start finish_up
    if( ! reset_initial &&
        ! add_all       &&
        ! add_remove    &&
        ! increment     &&
        ! call_react    &&
        ! reset_initial)
    {
 
//  Determine which system is being considered.
 
//  start molewt_1
      if(alter)
      {
        increment = TRUE ;
      }
      else if(plus && tri != 0.0 && newmol < newmol_test)
      {
        store_values = TRUE ;
      }
      else if(plus && tri != 0.0 && newmol >= newmol_test)
      {
        add_or_increment = TRUE ;
      }
      else if(plus)
      {
        add_new = TRUE ;
      }
      else if(systm && newmol >= newmol_test)
      {
        add_or_increment = TRUE ;
      }
      else
      {
        removed = 1 ;
        ll = 1 ;
        store_values = TRUE ;
      } // end molewt_1 ;
    } // end finish_up
 
//  start process_flags
 
    if(add_or_increment)
    {
 
//     Should a new reactant be added or should the constant reactant
//     be incremented.
 
//    fprintf(output_unit," at add_or_increment\n");
      add_or_increment = FALSE ;
      increment_add();
 
    }
    else if(increment)
    {
 
//     Increment moles between included reactants.
 
//    fprintf(output_unit," at increment\n");
      increment = FALSE ;
      increment_only();
      react(eof_pt,output_unit,input_unit);
 
//    start molewt_21
      if(*eof_pt !=0)
      {
//      exit outer_loop
        break;
      } // end molewt_21 ;
      increment = FALSE ;
      reset_initial = TRUE ;
 
    }
    else if(add_remove)
    {
 
//     Add and remove one reactant.
 
//    fprintf(output_unit," at add_remove\n");
      add_remove = FALSE ;
      remove_and_add();
 
//     Store base(i), top(i), and add(i) values.
 
    }
    else if(store_values)
    {
//    fprintf(output_unit," at store_values\n");
//    start molewt_35
      for(i = 0;i<nreac;i++)
      {
        bb [i] = base [i] ;
        to [i] = top [i] ;
        aa [i] = add [i] ;
      } // molewt_35 ;
      store_values = FALSE ;
//    start molewt_36
      if(systm)
      {
        increment = TRUE ;
      }
      else
      {
        add_all = TRUE ;
      } // end molewt_36 ;
 
    }
    else if(add_all)
    {
 
//     Add all the reactants to the system, or increment each reactant.
 
//    fprintf(output_unit," at add_all\n");
      add_all = FALSE ;
      nmix = nmix+1 ;
      tim = nmix ;
//    start molewt_37
      if(plus)
      {
        nrat = 1 ;
      } // end molewt_37 ;
      base [0] = bb [0]-(tri*tim) ;
      top [0] = to [0] ;
      add [0] = -tri ;
      base [switch_index] = bb [switch_index] ;
      top [switch_index] = to [switch_index]-(tri*tim) ;
      add [switch_index] = tri ;
      base [third_base] = aa [third_add]+(tri*tim) ;
      top [third_top] = base [third_base] ;
      add [third_add] = 0. ;
      call_react = TRUE ;
    }
    else if(add_new)
    {
 
//     Add a new reactant to the system.
 
//    fprintf(output_unit," at add_new\n");
      add_reactant = FALSE ;
//    start molewt_38
      for(i = 0;i<nreac;i++)
      {
//      start molewt_39
        for(;;)
        {
//        start molewt_40
          if(base [i] != 0.0)
          {
//          exit molewt_39
            break;
          } // end molewt_40 ;
          base [i] = add [i] ;
          add [i] = 0.0 ;
          add_reactant = TRUE ;
//        exit molewt_39
          break;
        } // molewt_39 ;
//      start molewt_41
        if(add_reactant)
        {
//        exit molewt_38
          break;
        } // end molewt_41 ;
      } // molewt_38 ;
      add_new = FALSE ;
      call_react = TRUE ;
 
    }
    else if(call_react)
    {
 
//     Call function react to compute new mole fraction.
 
//    fprintf(output_unit," at call_react\n");
      react(eof_pt,output_unit,input_unit);
//    start molewt_42
      if(*eof_pt != 0)
      {
//      exit outer_loop
        break;
      } // end molewt_42 ;
//    start molewt_43
      if(npt < 0)
      {
        npt = 0 ;
      } // end molewt_43 ;
      bypass_search = FALSE ;
//    start molewt_44
      if(plus && tri != 0.0 && newmol >= newmol_test)
      {
        bypass_search = TRUE ;
      }
      else
      {
//      start molewt_45
        if(systm && nmix > 1)
        {
          bypass_search = TRUE ;
        } // end molewt_45 ;
      } // end molewt_44 ;
 
//     Call function search to add or subtract species.
 
//    start molewt_46
      if( ! bypass_search)
      {
        search(eof_pt,scratch_unit,therm_unit,output_unit);
//      start molewt_47
        if(*eof_pt !=0)
        {
//        exit outer_loop
          break;
        } // end molewt_47 ;
      } // end molewt_46 ;
      call_react = FALSE ;
      reset_initial = TRUE ;
 
    }
    else if(reset_initial)
    {
 
//     Reset the initial estimates.
 
//    fprintf(output_unit," at reset_initial\n");
      npt = 0 ;
 
//      fprintf(output_unit," Alter, Systm, Plus = ");
//    start alter_conversion
//      if(alter)
//      {
//        fprintf(output_unit,"T ");
//      }
//      else
//      {
//        fprintf(output_unit,"F ");
//      } // end alter_conversion
//    start systm_conversion
//      if(systm)
//      {
//        fprintf(output_unit,"T ");
//      }
//      else
//      {
//        fprintf(output_unit,"F ");
//      } // end systm_conversion
//    start plus_conversion
//      if(plus)
//      {
//        fprintf(output_unit,"T \n");
//      }
//      else
//      {
//        fprintf(output_unit,"F \n");
//      } // end plus_conversion
 
      init() ;
      reset_initial = FALSE ;
      process_loop = FALSE ;
    } // end process_flags
  } while (process_loop) ;
//end outer_loop
} // end function molewt
 
//###
//*****************************************************************************
//                                                                            *
//     Function Name: omit_insert                                             *
//                                                                            *
//     File Name: omit_insert.cpp                                             *
//                                                                            *
//     Purpose:  Function called from search_coeff.  It processes the         *
//               omitted and inserted species and selects the species to be   *
//               included in the current system or put into temporary         *
//               storage.                                                     *
//                                                                            *
//     Author:    Robert L. Reese                                             *
//                                                                            *
//     Date: January 07, 2004                                                 *
//                                                                            *
//     Language: Borland C++ Compiler, Version 4.5                            *
//                                                                            *
//     Hardware: IBM Compatible PC                                            *
//                                                                            *
//     Operating System: Microsoft Windows 98                                 *
//                                                                            *
//     Update History:                                                        *
//                                                                            *
//     Name        Date      Revision          Changes                        *
//     ----        ----      ---------         -------                        *
//                                                                            *
//    R.L. Reese   01/07/04      0             Initial Release                *
//                                                                            *
//----------------------------------------------------------------------------*
//                                                                            *
//  Usage:                                                                    *
//                                                                            *
//  omit_insert((break_flag_pt)                                               *
//     where:                                                                 *
//                                                                            *
//     break_flag_pt - Logical pointer to flag indicating if calling          *
//                     function needs to exit loop                            *
//                                                                            *
//----------------------------------------------------------------------------*
//                                                                            *
//     Definitions:                                                           *
//                                                                            *
//     Local Variables:                                                       *
//     ----------------                                                       *
//                                                                            *
//     exit_loop      - Logical variable used to decide when to exit a loop   *
//     i              - Integer value used as do loop index                   *
//     j              - Integer value used as do loop index                   *
//     k              - Integer value used as do loop index                   *
//     kount          - Integer value used to compute maxns                   *
//     kount1         - Integer value used as index for omits and inserts     *
//     l              - Integer value used as do loop index                   *
//     ni             - Integer value used as first index of array a for      *
//                      specie being considered                               *
//     nmin_test      - Logical variable used to test for minimum n           *
//                                                                            *
//     Pointers                                                               *
//     -----------------                                                      *
//                                                                            *
//     None                                                                   *
//                                                                            *
//     Global Variables:                                                      *
//     -----------------                                                      *
//                                                                            *
//     Used by search_coeff and omit_insert                                   *
//     b            - Real array containing stoichiometric coefficients of    *
//                    ith element read from the coefficient records           *
//     begin_loop   - Logical variable used as loop control                   *
//     insert_flag  - Integer array used to flag legal insert species         *
//     jj           - Integer value used as counter for the number of         *
//                    elements to be put into temporary storage               *
//     lmt          - Character*4 array used as alphanumeric storage for the  *
//                    name of elements to be put into temporary storage in    *
//                    systm and plus problems                                 *
//     maxns        - Integer value used as maimum number of species allowed  *
//                    to be considered in the present system                  *
//     mt           - Character*4 array used as alphanumeric symbol of the    *
//                    elemen read from the coefficient records                *
//     nam          - Character*4 array used as alphanumeric symbol of        *
//                    species name read from the coefficient records          *
//     nm           - Integer value used as number of species put into        *
//                    temporary storage                                       *
//     omit_flag    - Integer array used to flag legal omit species           *
//     return_to_top- Logical variable used to test for returning to top of   *
//                    loop                                                  *
//     skip_to      - Logical variable to test for skipping code              *
//     test         - Logical variable used to test if species is to be put   *
//                    into temporary storage                                  *
//     indx.h                                                                 *
//     newmol       - Integer value used to determine if species should be    *
//                    input from the scratch file or the thermo file.         *
//     nlm          - Integer value used as number of elements in the system  *
//     nmin         - Integer value used as number of elements that are put   *
//                    into temporary storage for systm and plus problems      *
//     nomit        - Integer value containing number of species to be omitted*
//     ns           - Integer value of number of species being considered in  *
//                    the system                                              *
//     nsert        - Integer number of species listed in insert records      *
//     misc.h                                                                 *
//     llmt         - Character*4 array containing alphameric symbols for the *
//                    elements in the system being studied                    *
//     vol          - Real array containing molar volume of species           *
//     speces.h                                                               *
//     a            - Real array containing stoichiometric coefficients of ith*
//                    element in jth species                                  *
//     ensert       - Character*4 array of ith specie of n species which are  *
//                    specified to be considered or input                     *
//     juse         - Integer array containing variables used to identify     *
//                    which phase of a condensed specie should be used or if  *
//                    the spcie fits the T range being considered             *
//     omit         - Character*4 array containing names of specie to be      *
//                    omitted from current system                             *
//                                                                            *
//     Constants:                                                             *
//     ---------                                                              *
//                                                                            *
//     local                                                                  *
//     FALSE            - Defines value of logical false = 0                  *
//     TRUE             - Defines value of logical true = 1                   *
//     arrays.h                                                               *
//     b_size           - Integer variable used as size of b array            *
//     element_max      - Integer value of maximum number of elements allowed *
//     mt_size          - Integer variable used as size of mt array           *
//     size_char        - Integer value used as len of char to replace default*
//                        length of 4 characters in one word for FORTRAN 4 and*
//                        FORTRAN 77 conversion                               *
//     species_max      - Integer value of maximum number of species allowed  *
//     sub_size         - Integer variable used as maximum number of          *
//                        characters allowed in a name                        *
//     constant.h                                                             *
//     iblank           - Character*1 containing a blank                      *
//     first_start      - Integer value of start of first name                *
//     jj_bug_fix       - Integer value of fix for bug in temporary storage   *
//                        counter                                             *
//     kount_reset      - Integer value used to reset kount variable          *
//                                                                            *
//     Include Files:                                                         *
//     -------------                                                          *
//                                                                            *
//     arrays.h          - Contains sizes for include files                   *
//     constant.h        - Contains constants and sizes                       *
//     indx.h            - Contains index variables                           *
//     misc.h            - Contains miscellaneous variables                   *
//     speces.h          - Contains species variables                         *
//                                                                            *
//                                                                            *
//     External Software:                                                     *
//     -----------------                                                      *
//                                                                            *
//     None                                                                   *
//                                                                            *
//     Files:                                                                 *
//     -----                                                                  *
//                                                                            *
//       Input Files:                                                         *
//                                                                            *
//       None                                                                 *
//                                                                            *
//       Intermediate Files:                                                  *
//                                                                            *
//       None                                                                 *
//                                                                            *
//       Output Files:                                                        *
//                                                                            *
//       None                                                                 *
//                                                                            *
//*****************************************************************************
//###
 
#include "arrays.h"
#include "indx.h"
#include "misc.h"
#include "speces.h"
 
#define FALSE 0
#define TRUE 1
 
 
extern  int begin_loop ;
extern  int jj;
extern  int maxns ;
extern  int nm ;
extern  int return_to_top ;
extern  int skip_to ;
extern  int test ;
 
extern float b[b_size] ;
 
extern int   insert_flag[species_max];
extern int   omit_flag[species_max];
 
extern char  lmt[size_char][element_max] ;
extern char  mt[size_char][mt_size] ;
extern char  nam[size_char][sub_size] ;
 
 
void omit_insert(int *break_flag_pt)
{
 
  int exit_loop ;
  int i;
  int j;
  int k;
  int kount ;
  int kount1;
  int l ;
  int ni ;
  int nmin_test ;
 
#include "constant.h"
 
  *break_flag_pt = FALSE;
 
//start omit_insert_12
  for(;;)
  {
//  start omit_insert_19
    if( ! begin_loop)
    {
      *break_flag_pt = TRUE;
//    exit omit_insert_12
      break;
    } // end omit_insert_19
 
//  set jj to fix bug
    jj = jj_bug_fix ;
 
//  start omit_insert_20
    if(newmol < 1)
    {
 
//    start omit_insert_21
      if(nomit != 0)
      {
 
//  Is this species to be ommited?
 
        return_to_top = FALSE ;
//      start omit_insert_22
        for(i=0;i<nomit;i++)
        {
          kount = 0;
          kount1 = 0;
          exit_loop = FALSE ;
//        start omit_insert_23
          for(j=0;j<first_start;j++)
          {
//          start omit_insert_24
            if(omit[j][i] != iblank)
            {
//            start omit_insert_24a
              if(omit[j][i] != nam[kount][kount1])
              {
                exit_loop = TRUE ;
//              exit omit_insert_23
                break;
              } // end omit_insert_24a
            }
            else
            {
//            exit omit_insert_23
              break;
            } // end omit_insert_24
            kount++;
//          start omit_insert_24b
            if(kount == kount_reset)
            {
              kount = 0;
              kount1++;
            } // end omit_insert_24b
          } // end omit_insert_23
//        start omit_insert_25
          if( ! exit_loop)
          {
            return_to_top = TRUE ;
            omit_flag[i] = 1;
//          exit omit_insert_22
            break;
          } // end omit_insert_25
        } // end omit_insert_22
//      start omit_insert_26
        if(return_to_top)
        {
          *break_flag_pt = TRUE;
//        exit omit_insert_12
          break;
        } // end omit_insert_26
      } // end omit_insert_21
 
//    start omit_insert_27
      if(nsert != 0)
      {
 
//  Is this species to be included?
 
//      start omit_insert_28
        for(i=0;i<nsert;i++)
        {
          kount = 0;
          kount1 = 0;
          skip_to = FALSE ;
          exit_loop = FALSE ;
//        start omit_insert_29
          for(j=0;j<first_start;j++)
          {
//          start omit_insert_30
            if(ensert[j][i] != iblank)
            {
//            start omit_insert_30a
              if(ensert[j][i] != nam[kount][kount1])
              {
                exit_loop = TRUE ;
//              exit omit_insert_29
                break;
              } // end omit_insert_30
            }
            else
            {
//            exit omit_insert_29
              break;
            } // end omit_insert_30
            kount++;
//          start omit_insert_30b
            if(kount == kount_reset)
            {
              kount = 0;
              kount1++;
            } //  end omit_insert_30b
          } // end omit_insert_29
//        start omit_insert_31
          if( ! exit_loop)
          {
            skip_to = TRUE ;
            insert_flag[i] = 1;
//          exit omit_insert_28
            break;
          } // end omit_insert_31
        } // end omit_insert_28
//      start omit_insert_32
        if( ! skip_to)
        {
          *break_flag_pt = TRUE;
//        exit omit_insert_12
          break;
        } // end omit_insert_32
      } // end omit_insert_27
    } // end omit_insert_20
    test = FALSE ;
//  start omit_insert_33
    for(k=0;k<mt_size;k++)
    {
      return_to_top = FALSE ;
      nmin_test = FALSE ;
//    start omit_insert_34
      if(b[k] == 0.0)
      {
//      exit omit_insert_33
        break;
      } // end omit_insert_34
//    start omit_insert_35 :
      for(i=0;i<nlm;i++)
      {
        ni = i ;
//      start omit_insert_36
        if(llmt[0][i] == mt[0][k] &&
           llmt[1][i] == mt[1][k])
        {
          nmin_test = TRUE ;
//        exit omit_insert_35
          break;
        } // end omit_insert_36
      } // end omit_insert_35
//    start omit_insert_37
      if( ! nmin_test)
      {
//      start omit_insert_38
        if(jj != 0)
        {
//        start omit_insert_39
          for(i=0;i<jj;i++)
          {
            ni = i ;
//          start omit_insert_40
            if(lmt[0][i] == mt[0][k] &&
               lmt[1][i] == mt[1][k])
            {
              nmin_test = TRUE ;
//            exit omit_insert_39
              break;
            } // end omit_insert_40
          } // end omit_insert_39
        } // end omit_insert_38
 
//  This species will not be included.
 
//      start omit_insert_41
        if( ! nmin_test)
        {
//        start omit_insert_42
          for(j=0;j<nlm;j++)
          {
            a[j][ns] = 0. ;
          } //  end omit_insert_42
// Special case for systm problem after first set of mole fractions
//        start omit_insert_42a
          if(!systm)
          {
            return_to_top = TRUE ;
//          exit omit_insert_33
            break;
          }
          else if(newmol <=0)
          {
            return_to_top = TRUE;
//          exit omit_insert_33
            break;
          } // end search 42a
        } // end omit_insert_41
      } // end omit_insert_37
//    start omit_insert_43
      if(nmin >= 1)
      {
        skip_to = FALSE ;
//      start omit_insert_44
        for(l=0;l<jj;l++)
        {
//        start omit_insert_45
          if(llmt[0][l] == mt[0][k] &&
             llmt[1][l] == mt[1][k])
          {
            skip_to = TRUE ;
//           exit omit_insert_44
            break;
          } // end omit_insert_45
        } // end omit_insert_44
//      start omit_insert_46
        if(! skip_to)
        {
          test = TRUE ;
        } // end omit_insert_46
      } // end omit_insert_43
//    start omit_insert_47
      if((ns+nm) < maxns)
      {
        a[ni][ns] = b[k] ;
      } // end omit_insert_47
    } // end omit_insert_33
//  exit omit_insert_12
    break;
 
  } // end omit_insert_12
 
} // end function omit_insert
//###
//*****************************************************************************
//                                                                            *
//     Function Name: out1                                                    *
//                                                                            *
//     File Name: out1.cpp                                                    *
//                                                                            *
//     Purpose:  Prints out the results of up to nine P-T points.  Also prints*
//               out the mole fractions of species with values less than      *
//               5 x 10-6.  Out1 computes the values nd then prints out the   *
//               log 10 values                                                *
//                                                                            *
//     Author:    Robert L. Reese                                             *
//                                                                            *
//     Date: June 01, 1973                                                    *
//                                                                            *
//     Language: Borland C++ Compiler, Version 4.5                            *
//                                                                            *
//     Hardware: IBM Compatible PC                                            *
//                                                                            *
//     Operating System: Microsoft Windows 98                                 *
//                                                                            *
//     Update History:                                                        *
//                                                                            *
//     Name        Date      Revision          Changes                        *
//     ----        ----      ---------         -------                        *
//                                                                            *
//    R.L. Reese   06/01/73      0             Initial Release                *
//    J. Holloway  07/23/89      A             Added IW                       *
//    R.L. Reese   11/23/03      B             Converted to C++               *
//                                                                            *
//----------------------------------------------------------------------------*
//                                                                            *
//  Usage:                                                                    *
//                                                                            *
//  out1(output_unit)    where:                                               *
//                                                                            *
//     output_unit - FILE pointer to output listing logical unit              *
//                                                                            *
//----------------------------------------------------------------------------*
//                                                                            *
//     Definitions:                                                           *
//                                                                            *
//     Local Variables:                                                       *
//     ----------------                                                       *
//                                                                            *
//     fug       - Real array of fugacity values                              *
//     i         - Integer variable used as do loop index                     *
//     j         - Integer variable used as do loop index                     *
//     jj        - Integer variable used as do loop index                     *
//     jt        - Integer variable used to save value of j                   *
//     k         - Integer variable used as do loop index                     *
//     n         - Integer variable used as do loop index                     *
//     not_empty - Logical variable used to determine if log10 mole values    *
//                 are present for output                                     *
//     qfmrel    - Real array for printout of qfi fo2 values                  *
//     tc        - Real array for current temperatures in Celsius             *
//     xiw       - Real variable used to calculate qfmrel                     *
//                                                                            *
//     Pointers                                                               *
//     -----------------                                                      *
//                                                                            *
//     None                                                                   *
//                                                                            *
//     Global Variables:                                                      *
//     -----------------                                                      *
//     dltra     - Double precision variable used as dlog10 of dtra           *
//     dtra      - Double precision variable used as limit for trace species  *
//     loop_exit - Logical variable used to exit do loop                      *
//     ncond_out - Integer variable used as counter for mole fractions of gas *
//     td        - Double precision variable for totn and gtotn calculations  *
//     v         - Double precision array used to store mole fraction for     *
//                 output                                                     *
//     indx.h                                                                 *
//     kase         - Integer value of identification number of current       *
//                    problem being considered                                *
//     npt          - Integer value of index for data saved for output        *
//     nreac        - Integer number of reactant records                      *
//     ns           - Integer value of number of species being considered in  *
//                    the system                                              *
//     misc.h                                                                 *
//     anum         - Real value of stochiometric coefficients of elements on *
//                    reactant records                                        *
//     base         - Real value of first value of mole fraction of ith       *
//                    reactant                                                *
//     name         - Character*4 symbols for the elements in the system      *
//                    being studied                                           *
//     points.h                                                               *
//     fo2c         - Real array of fugacity of O2C                           *
//     fo2h         - Real array of fugacity of O2H                           *
//     free_energy  - Real array of free energy                               *
//     gas_input    - Logical value of Gas Input flag                         *
//     gtotn        - Real array of values used in gas mole fractions         *
//     newest       - Logical value for first time in routine                 *
//     oxfug        - Real array of oxygen fugacity values                    *
//     ppp          - Real value of static P stored  for output               *
//     totn         - Real value used in calculation of mole fractions        *
//     ttt          - Real value of assigned T in degrees Kelvin for current  *
//                    point                                                   *
//     speces.h                                                               *
//     en           - Real value of Kg-moles of jth species/Kg of mixture.    *
//     sub          - Character*4 names of species read from thermo data file.*
//                                                                            *
//     Constants:                                                             *
//     ---------                                                              *
//                                                                            *
//     local                                                                  *
//     FALSE           - Defines value of logical false = 0                   *
//     TRUE            - Defines value of logical true = 1                    *
//     arrays.h                                                               *
//     npt_out         - Integer value of total points output on a page       *
//     reacts_max      - Integer value of maximum mole values per record      *
//     constant.h                                                             *
//     ib              - Character contining blank                            *
//     jt_value        - Integer value of maximum size of output line values  *
//     kelvin          - Real factor for Kelvin - Centigrade conversion       *
//     oxfug_val       - Real value of maximum size of oxygen fugacity        *
//     qfi_factor_1    - Real value of first value in calculating qfi         *
//     qfi_factor_2    - Real value of second value in calculating qfi        *
//     qfi_factor_3    - Real value of third value in calculating qfi         *
//                                                                            *
//     Include Files:                                                         *
//     -------------                                                          *
//                                                                            *
//     arrays.h          - Contains sizes for include files                   *
//     constant.h        - Contains out1 constants                            *
//     indx.h            - Contains index variables                           *
//     math.h            - Contains C++ math routines                         *
//     misc.h            - Contains miscellaneous variables                   *
//     points.h          - Contains temperature and pressure variables        *
//     speces.h          - Contains species variables                         *
//     stdio.h           - Contains C++ I/O routines                          *
//                                                                            *
//                                                                            *
//     External Software:                                                     *
//     -----------------                                                      *
//                                                                            *
//     log10         - C++ math library routine to find log value             *
//     out1_moles    - Function to output mole fractions total and gas        *
//                                                                            *
//     Files:                                                                 *
//     -----                                                                  *
//                                                                            *
//       Input Files:                                                         *
//                                                                            *
//       None                                                                 *
//                                                                            *
//       Intermediate Files:                                                  *
//                                                                            *
//       None                                                                 *
//                                                                            *
//       Output Files:                                                        *
//                                                                            *
//       output_unit - Integer value of logical unit for output listing       *
//                                                                            *
//*****************************************************************************
//###
 
#include "arrays.h"
#include "indx.h"
#include <math.h>
#include "misc.h"
#include "points.h"
#include "speces.h"
#include <stdio.h>
 
#define FALSE 0
#define TRUE 1
 
void out1_moles(FILE *output_unit);
 
double dltra ;
double dtra ;
int    loop_exit ;
int    ncond_out ;
double td ;
double v[npt_out] ;
 
void out1(FILE *output_unit)
{
 
#include "constant.h"
 
  float fug[npt_out] ;
  float qfmrel[npt_out] ;
  float tc[npt_out] ;
 
  float xiw ;
 
  int i ;
  int j ;
  int jj ;
  int jt ;
  int k ;
  int n ;
 
// logical variables
 
  int not_empty ;
 
// Zero out fug
 
//start zero_fug
  for(i=0;i<=npt;i++)
  {
    fug[i] = 0.0;
  } // end zero_fug
 
//start out1_1
  if(newest)
  {
    fprintf(output_unit,"\t\t      THERMODYNAMIC EQUILIBRIUM PROPERTIES AT\n\n");
    fprintf(output_unit,"\t\t         ASSIGNED TEMPERATURE AND PRESSURE\n\n");
    fprintf(output_unit,(" CASE NO. %8d\n"),kase);
    fprintf(output_unit," INPUT COMPOSITION\t\t\t\t\t\t MOLES \n");
    ncond_out  =  0 ;
//  start out1_2
    for( n = 0; n<nreac;n++)
    {
      loop_exit = FALSE ;
//    start out1_3
      for( j = 0; j<reacts_max;j++)
      {
        jt = j ;
//      start out1_4
        if(name[0][n][j] == ib)
        {
          loop_exit = TRUE ;
//        exit out1_3
          break;
        } // end out1_4
      } // end out1_3
//    start out1_5
      if( ! loop_exit)
      {
        jt = jt_value ;
      } // end out1_5
      j = jt-1 ;
 
      fprintf(output_unit," ");
//    start fprintf_1
      for(jj=0;jj<5;jj++)
      {
        fprintf(output_unit,"%c%c ",name[0][n][jj],
                                       name[1][n][jj]);
        fprintf(output_unit,"%7.4f ",anum[n][jj]);
      } // end fprintf_1
      fprintf(output_unit,"     %9.3f\n",base[n]);
 
    } // end out1_2
    newest = FALSE;
  } // end out1_1
  fprintf(output_unit,"\n\n THERMODYNAMIC PROPERTIES\n");
 
//     Pressure.
 
//start fprintf_2
  for(j=0;j<=npt;j++)
  {
//  start fprintf_3
    if(j == 0)
    {
      fprintf(output_unit,"\n P,BARS\t\t         %8.0f.",ppp[j]);
    }
    else
    {
      fprintf(output_unit,"%8.0f.",ppp[j]);
    } // end fprintf_3
  } // end fprintf_2
  fprintf(output_unit,"\n");
 
//     Temperature.
 
//start out1_6
  for( i = 0;i<=npt;i++)
  {
    tc[i] = ttt[i] - kelvin ;
  } // end out1_6
//start fprintf_4
  for(j=0;j<=npt;j++)
  {
//  start fprintf_5
    if(j == 0)
    {
      fprintf(output_unit," T, DEG. C\t\t\t %8.0f.",tc[j]);
    }
    else
    {
      fprintf(output_unit,"%8.0f.",tc[j]);
    } // end fprintf_5
  } // end fprintf_4
  fprintf(output_unit,"\n");
 
//     Free energy.
 
//start fprintf_6
  for(j=0;j<=npt;j++)
  {
//  start fprintf_7
    if(j == 0)
    {
      fprintf(output_unit," TOT. DELTA G\t\t\t %9.3f",free_energy[j]);
    }
    else
    {
      fprintf(output_unit,"%9.3f",free_energy[j]);
    } // end fprintf_7
  } // end fprintf_6
  fprintf(output_unit,"\n");
 
//     Mole fractions - equilibrium.
 
  out1_moles(output_unit);
 
  not_empty = FALSE ;
//start out1_20
  for( k = 0;k < ns;k++)
  {
    loop_exit = FALSE ;
//  start out1_21
    for( i = 0;i<=npt;i++)
    {
//    start out1_22
      if( ! gas_input)
      {
        td = totn[i] ;
        v[i] = en[k][i]/td ;
        loop_exit = TRUE ;
      } // end out1_22
//    start out1_23
      if( ! loop_exit)
      {
        td = gtotn[i] ;
        v[i] = en[k][i]/td ;
      } // end out1_23
    } // end out1_21
//  start out1_24
    for( i = 0;i<=npt;i++)
    {
//    start out1_25
      if(v[i] != 0.0)
      {
        v[i] = log10(v[i]) ;
//      start out1_26
        if(sub[0][0][k] == 'O' && sub[1][0][k] == '2')
        {
          fug[i] = v[i]+ log10(ppp[i]) ;
        } // end out1_26
//
      } // end out1_25
    } // end out1_24
    loop_exit = FALSE ;
//  start out1_27
    for( j = 0;j<=npt;j++)
    {
//    start out1_28
      if(v[j] < dltra)
      {
        loop_exit = TRUE ;
//      exit out1_27
        break;
      } // end out1_28
    } // end out1_27
//  start out1_29
    if(loop_exit)
    {
//    start out1_30
      if( ! not_empty)
      {
        fprintf(output_unit,"\n LOG(10) (MOLE FRACTION)\n");
      } // end out1_30
      fprintf(output_unit," ");
      not_empty = TRUE;
//    start fprintf_12
      for(i=0;i<8;i++)
      {
        fprintf(output_unit,"%c%c%c%c",sub[0][i][k],
                                       sub[1][i][k],
                                       sub[2][i][k],
                                       sub[3][i][k]);
      } // end fprintf_12
//    start fprintf_13
      for(i=0;i<=npt;i++)
      {
        fprintf(output_unit,"%9.3f",v[i]);
      } // end fprintf_13
      fprintf(output_unit,"\n");
    } // end out1_29
  } // end out1_20
 
//       Oxygen fugacity calculations
 
//start out1_31
  for( n = 0;n<=npt;n++)
  {
 
//      QFM = 9.0-25738./ttt(n)+0.092*(ppp(n)-1.)/ttt(n)
//       This version uses qfi (HUEBNER) as reference value
 
    xiw = qfi_factor_1 - qfi_factor_2/ttt[n]+qfi_factor_3*
          (ppp[n]-1.0)/ttt[n] ;
    qfmrel[n] = xiw ;
//  start out1_32
    if(ppp[n] >= 1.0 && v[n] >= v_limit)
    {
//    start out1_33
      if(oxfug[n] > oxfug_val)
      {
        oxfug[n] = oxfug_val ;
      } // end out1_33
      fug[n] = fug[n]+oxfug[n] ;
    } // end out1_32
  } // end out1_31
  fprintf(output_unit,"\n OXYGEN FUGACITY CALCULATIONS:\n");
  fprintf(output_unit,
  " FO2 = XO2*FGO2\t\t   ");
// start fprintf_14
  for(n=0;n<=npt;n++)
  {
    fprintf(output_unit,"%7.3f  ",fug[n]);
  } // end fprintf_14
  fprintf(output_unit,"\n");
 
  fprintf(output_unit,
  " FO2 = CO/CO2  \t\t   ");
// start fprintf_15
  for(n=0;n<=npt;n++)
  {
    fprintf(output_unit,"%7.3f  ",fo2c[n]);
  } // end fprintf_15
  fprintf(output_unit,"\n");
 
  fprintf(output_unit,
  " FO2 = H2/H2O  \t\t   ");
// start fprintf_16
  for(n=0;n<=npt;n++)
  {
    fprintf(output_unit,"%7.3f  ",fo2h[n]);
  } // end fprintf_16
  fprintf(output_unit,"\n");
 
  fprintf(output_unit,
  " QFI FO2       \t\t   ");
// start fprintf_17
  for(n=0;n<=npt;n++)
  {
    fprintf(output_unit,"%7.3f  ",qfmrel[n]);
  } // end fprintf_17
  fprintf(output_unit,"\n");
 
} // end function out1
//###
//*****************************************************************************
//                                                                            *
//     Function Name: out1_moles                                              *
//                                                                            *
//     File Name: out1_moles.cpp                                              *
//                                                                            *
//     Purpose:  Used by out1 to output mole fractions and mole fractions for *
//               gases.                                                       *
//                                                                            *
//     Author:    Robert L. Reese                                             *
//                                                                            *
//     Date: January 29, 2004                                                 *
//                                                                            *
//     Language: Borland C++ Compiler, Version 4.5                            *
//                                                                            *
//     Hardware: IBM Compatible PC                                            *
//                                                                            *
//     Operating System: Microsoft Windows 98                                 *
//                                                                            *
//     Update History:                                                        *
//                                                                            *
//     Name        Date      Revision          Changes                        *
//     ----        ----      ---------         -------                        *
//                                                                            *
//    R.L. Reese   01/29/04      0             Initial Release                *
//                                                                            *
//----------------------------------------------------------------------------*
//                                                                            *
//  Usage:                                                                    *
//                                                                            *
//  out1_moles(output_unit)    where:                                         *
//                                                                            *
//     output_unit - FILE pointer to output listing logical unit              *
//                                                                            *
//----------------------------------------------------------------------------*
//                                                                            *
//     Definitions:                                                           *
//                                                                            *
//     Local Variables:                                                       *
//     ----------------                                                       *
//                                                                            *
//     i         - Integer variable used as do loop index                     *
//     k         - Integer variable used as do loop index                     *
//                                                                            *
//     Pointers                                                               *
//     -----------------                                                      *
//                                                                            *
//     None                                                                   *
//                                                                            *
//     Global Variables:                                                      *
//     -----------------                                                      *
//     local - used by out1 and out1_moles                                    *
//     dltra     - Double precision variable used as dlog10 of dtra           *
//     dtra      - Double precision variable used as limit for trace species  *
//     loop_exit - Logical variable used to exit do loop                      *
//     ncond_out - Integer variable used as counter for mole fractions of gas *
//     td        - Double precision variable for totn and gtotn calculations  *
//     v         - Double precision array used to store mole fraction for     *
//                 output                                                     *
//     indx.h                                                                 *
//     npt          - Integer value of index for data saved for output        *
//     ns           - Integer value of number of species being considered in  *
//                    the system                                              *
//     points.h                                                               *
//     gas_input    - Logical value of Gas Input flag                         *
//     gtotn        - Real array of values used in gas mole fractions         *
//     totn         - Real value used in calculation of mole fractions        *
//     speces.h                                                               *
//     en           - Real value of Kg-moles of jth species/Kg of mixture.    *
//     iuse         - Integer value used to identify status of ith species    *
//     sub          - Character*4 names of species read from thermo data file.*
//                                                                            *
//     Constants:                                                             *
//     ---------                                                              *
//                                                                            *
//     local                                                                  *
//     FALSE           - Defines value of logical false = 0                   *
//     TRUE            - Defines value of logical true = 1                    *
//     arrays.h                                                               *
//     npt_out         - Integer value of total points output on a page       *
//     constant.h                                                             *
//     trace_limit     - Real value of limit of trace                         *
//                                                                            *
//     Include Files:                                                         *
//     -------------                                                          *
//                                                                            *
//     arrays.h          - Contains sizes for include files                   *
//     constant.h        - Contains out1 constants                            *
//     indx.h            - Contains index variables                           *
//     math.h            - Contains C++ math routines                         *
//     points.h          - Contains temperature and pressure variables        *
//     speces.h          - Contains species variables                         *
//     stdio.h           - Contains C++ I/O routines                          *
//                                                                            *
//                                                                            *
//     External Software:                                                     *
//     -----------------                                                      *
//                                                                            *
//     log10         - C++ math library routine to find log value             *
//                                                                            *
//     Files:                                                                 *
//     -----                                                                  *
//                                                                            *
//       Input Files:                                                         *
//                                                                            *
//       None                                                                 *
//                                                                            *
//       Intermediate Files:                                                  *
//                                                                            *
//       None                                                                 *
//                                                                            *
//       Output Files:                                                        *
//                                                                            *
//       output_unit - Integer value of logical unit for output listing       *
//                                                                            *
//*****************************************************************************
//###
 
#include "arrays.h"
#include "indx.h"
#include <math.h>
#include "points.h"
#include "speces.h"
#include <stdio.h>
 
#define FALSE 0
#define TRUE 1
 
extern double dltra ;
extern double dtra ;
extern int    loop_exit ;
extern int    ncond_out ;
extern double td ;
extern double v[npt_out] ;
 
void out1_moles(FILE *output_unit)
{
#include "constant.h"
 
  int i;
  int k;
 
//    Define limit for trace species
 
  dtra = trace_limit ;
  dltra = log10(dtra) ;
 
//     Mole fractions - equilibrium.
  fprintf(output_unit,"\n MOLE FRACTIONS(TOTAL)\n");
//start out1_7
  for( k = 0;k<ns;k++)
  {
//  start out1_8
    for( i = 0;i<=npt;i++)
    {
      td = totn[i] ;
      v[i] = en[k][i]/td ;
    } // end out1_8
    loop_exit = FALSE ;
//  start out1_9
    for( i = 0;i<=npt;i++)
    {
//    start out1_10
      if(v[i] >= dtra)
      {
        loop_exit = TRUE ;
//      exit out1_9
        break;
      } // end out1_10
    } // end out1_9
//  start out1_11
    if(loop_exit)
    {
      fprintf(output_unit," ");
//    start fprintf_8
      for(i=0;i<8;i++)
      {
        fprintf(output_unit,"%c%c%c%c",sub[0][i][k],
                                       sub[1][i][k],
                                       sub[2][i][k],
                                       sub[3][i][k]);
      } // end fprintf_8
//    start fprintf_9
      for(i=0;i<=npt;i++)
      {
        fprintf(output_unit,"%9.5f",v[i]);
      } // end fprintf_9
      fprintf(output_unit,"\n");
//    start out1_12
      if(iuse[k] > 0)
      {
        ncond_out = ncond_out + 1 ;
      } // end out1_12
    } // end out1_11
  } // end out1_7
//start out1_13
  if(gas_input && ncond_out != 0)
  {
    fprintf(output_unit,"\n MOLE FRACTIONS(GASES)\n");
//  start out1_14
    for( k = 0;k<ns;k++)
    {
//    start out1_15
      if(iuse[k] == 0)
      {
//      start out1_16
        for( i = 0;i<=npt;i++)
        {
          td = gtotn[i] ;
          v[i] = en[k][i]/td ;
        } // end out1_16
        loop_exit = FALSE ;
//      start out1_17
        for( i = 0;i<=npt;i++)
        {
//        start out1_18
          if(v[i] >= dtra)
          {
            loop_exit = TRUE ;
//          exit out1_17
            break;
          } // end out1_18
        } // end out1_17
//      start out1_19
        if(loop_exit)
        {
          fprintf(output_unit," ");
//        start fprintf_10
          for(i=0;i<8;i++)
          {
            fprintf(output_unit,"%c%c%c%c",sub[0][i][k],
                                           sub[1][i][k],
                                           sub[2][i][k],
                                           sub[3][i][k]);
          } // end fprintf_10
//        start fprintf_11
          for(i=0;i<=npt;i++)
          {
            fprintf(output_unit,"%9.5f",v[i]);
          } // end fprintf_11
          fprintf(output_unit,"\n");
        } // end out1_19
      } // end out1_15
    } // end out1_14
  } // end out1_13
 
} // end function out1_moles
//###
//*****************************************************************************
//                                                                            *
//     Function Name: process_pressure                                        *
//                                                                            *
//     File Name: process_pressure                                            *
//                                                                            *
//     Purpose:  Reads the pressures from the input file and outputs them to  *
//               the listing file. It checks for valid values. Pressures must *
//               be greater than zero bars.  Temperatures are converted to    *
//               Kelvin and must range from zero degrees to 5000 degrees      *
//               Kelvin. Blanks are not interpreted as zero.  If there are    *
//               nine values in a record, there must be a blank record        *
//               following that record.                                       *
//                                                                            *
//     Author:    Robert L. Reese                                             *
//                                                                            *
//     Date: January 08, 2004                                                 *
//                                                                            *
//     Language: Borland C++ Compiler, Version 4.5                            *
//                                                                            *
//     Hardware: IBM Compatible PC                                            *
//                                                                            *
//     Operating System: Microsoft Windows 98                                 *
//                                                                            *
//     Update History:                                                        *
//                                                                            *
//     Name        Date      Revision          Changes                        *
//     ----        ----      ---------         -------                        *
//                                                                            *
//    R.L. Reese   01/07/04      0             Initial Release                *
//                                                                            *
//----------------------------------------------------------------------------*
//                                                                            *
//  Usage:                                                                    *
//                                                                            *
//     process_pressure(output_unit,                                          *
//                      input_unit,                                           *
//                      break_value_pt,                                       *
//                      read_error_pt)                                        *
//      where:                                                                *
//                                                                            *
//     output_unit    - FILE pointer to logical output listing unit           *
//     input_unit     - FILE pointer to logical input unit                    *
//     break_flag_pt  - Integer logical pointer to value indicating if loop   *
//                      needs to be exited upon return                        *
//     read_error_pt  - Integer logical pointer indicating if error occurred  *
//                      on read                                               *
//                                                                            *
//----------------------------------------------------------------------------*
//                                                                            *
//     Definitions:                                                           *
//                                                                            *
//     Local Variables:                                                       *
//     ----------------                                                       *
//                                                                            *
//     adata          - Character buffer containing data returned from        *
//                      init_read                                             *
//     blank          - Character buffer containing blanks                    *
//     empty_file     - Integer value of status returned from init_read       *
//     eof            - Integer value of status returned from init_read       *
//     i              - Integer for loop index                                *
//     input_end      - Logical flag for end of temperature or pressure       *
//                      values                                                *
//     j              - Integer for loop index                                *
//     kount          - Integer value of index for temp and pressures         *
//     m              - Integer for loop index                                *
//     nbytes         - Integer value of number of bytes from init_read       *
//     nbytes_count   - Integer value used to determine if current p or       *
//                      t has been defaulted                                  *
//     nmove          - Integer value of number of digits to move from        *
//                      input record for reactant values                      *
//     pressure_index - Integer do loop index for processing pressure         *
//                      inputs                                                *
//     real_buf       - Real buffer used to contain current number to be      *
//                      converted from characters to float                    *
//     temp_index     - Integer do loop index for processing temperature      *
//                      inputs                                                *
//     tp_index       - Integer value of index for outputting pressure        *
//                      temperature values                                    *
//     tp_kount       - Integer value for number of temps or pressures        *
//                                                                            *
//                                                                            *
//     Pointers                                                               *
//     -----------------                                                      *
//                                                                            *
//     adata_pt      - Character pointer to data from init_read               *
//     empty_file_pt - Integer pointer to empty file value from               *
//                     init_read                                              *
//     end_ptr       - Character pointer to error in strtod                   *
//     eof_pt        - Integer pointer to eof value from init_read            *
//     nbytes_pt     - Integer pointer to number of bytes read from           *
//                     init_read                                              *
//     real_buf_pt   - Character pointer to buffer used for conversions       *
//                                                                            *
//     Global Variables:                                                      *
//     -----------------                                                      *
//                                                                            *
//     points.h                                                               *
//     p               - Real values of P in bars that have been input        *
//     pressure_number - Integer number of input pressure values              *
//     t               - Real values of T in Kelvin that have been input      *
//     temp_number     - Integer number of input temperature values           *
//                                                                            *
//     Constants:                                                             *
//     ---------                                                              *
//                                                                            *
//     local                                                                  *
//     FALSE           - Defines value of logical false = 0                   *
//     TRUE            - Defines value of logical true = 1                    *
//     arrays.h                                                               *
//     adata_len       - Integer variable used as length of adata array       *
//     blank_length    - Integer variable used as length of blank array       *
//     real_buf_len    - Integer variable used as size of real_buf            *
//     constants.h                                                            *
//     blank_index     - Integer value pointing to start of blanked field     *
//     count_decrement - Integer value to subtract from nbytes_count for      *
//                       pressure default test                                *
//     field_width     - Integer value of output field width for pressures    *
//     iblank          - Character*1 containing a blank                       *
//     initial_nbytes  - Integer value to ending position of first input      *
//                       pressure                                             *
//     kelvin          - Integer value of Kelvin conversion factor            *
//     max_out_values  - Integer value to number of output values per line    *
//     max_t_p         - Integer value of maximum number of values allowed    *
//     max_temp        - Real value of maximum temperature allowed            *
//     round_up        - Real value of value used to round up for below zero  *
//                       Kelvin test                                          *
//     string_end      - Character containing end of string control character *
//     t_p_default     - Real value of input temperature or pressure input    *
//     t_p_records     - Integer value of number of temp or press records     *
//     values_record   - Integer value of maximum values per input record     *
//                                                                            *
//     Include Files:                                                         *
//     -------------                                                          *
//                                                                            *
//     arrays.h          - Contains sizes for include files                   *
//     constant.h        - Contains program constants                         *
//     iostream.h        - Contains I/O routintes                             *
//     points.h          - Contains temperature and pressure variables        *
//     stdlib.h          - Contains C++ standard library routines             *
//     string.h          - Contains C++ string routines                       *
//                                                                            *
//                                                                            *
//     External Software:                                                     *
//     -----------------                                                      *
//                                                                            *
//     init_read        - Routine to read input records                       *
//     stdio.h           - Contains FILE declaration                          *
//     strncmp          - C++ routine to compare strings                      *
//     strncpy          - C++ routine to copy strings                         *
//     strtod           - C++ routine to convert from ascii to float          *
//                                                                            *
//     Files:                                                                 *
//     -----                                                                  *
//                                                                            *
//       Input Files:                                                         *
//                                                                            *
//       FILE input_unit  - Logical device for input reading                  *
//                                                                            *
//       Intermediate Files:                                                  *
//                                                                            *
//       None                                                                 *
//                                                                            *
//       Output Files:                                                        *
//                                                                            *
//       FILE output_unit - logical device for output listings                *
//                                                                            *
//*****************************************************************************
//###
 
#include "arrays.h"
#include <iostream.h>
#include "points.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
 
#define FALSE 0
#define TRUE  1
 
void init_read(char *idata,
               int *eof,
               FILE *opu,
               FILE *inp_unit,
               int *empty_file,
               int *nbytes);
 
 
void process_pressure(FILE *output_unit,
                      FILE *input_unit,
                      int *break_value_pt,
                      int *read_error_pt)
{
 
#include "constant.h"
 
  int empty_file;
  int eof;
  int i;
  int input_end;
  int j;
  int kount;
  int m;
  int nbytes;
  int nbytes_count;
  int nmove;
  int pressure_index ;
  int temp_index ;
  int tp_index;
  int tp_kount;
 
  char  blank[blank_length];
  char  adata[adata_len];
  char  real_buf[real_buf_len];
 
  char *adata_pt;
  int  *empty_file_pt;
  char *end_ptr = "abcde";
  int  *eof_pt;
  int  *nbytes_pt;
  char *real_buf_pt;
 
  adata_pt       = &adata[0];
  empty_file_pt  = &empty_file;
  eof_pt         = &eof;
  nbytes_pt      = &nbytes;
  real_buf_pt    = &real_buf[0];
 
//start initialize blank
  for(i=0;i<blank_length;i++)
  {
    blank[i] = iblank;
  } // end initialize blank
 
// process the input pressure values
 
//start do_while_loop
  for(;;)
  {
    tp_kount = 0;
    input_end = FALSE;
//  start read_pressure_values
    for(i=0;i<t_p_records;i++)
    {
      *eof_pt = 0;
      init_read(adata_pt,eof_pt,output_unit,input_unit,
                empty_file_pt,nbytes_pt);
//    start pressure_read_error
      if(*eof_pt != 0)
      {
//      exit read_pressure_values
        break ;
      } // end pressure_error
      kount = 0;
//    start convert_pressures
      nbytes_count = initial_nbytes;
      for(j=0;j<values_record;j++)
      {
//      start pressure_default_test
        if(*nbytes_pt < (nbytes_count-count_decrement))
        {
          p[tp_kount] = t_p_default;
          input_end = TRUE;
//        exit convert_pressures
          break;
        }
        else
        {
          strncpy(real_buf_pt,&blank[0],blank_length);
          *(real_buf_pt+blank_length) = string_end;
//        start compute_pressure_moves
          if(*nbytes_pt > (nbytes_count-1))
          {
            nmove = field_width;
          }
          else
          {
            nmove = *nbytes_pt - (nbytes_count-field_width);
          } // end compute_pressure_moves
          strncpy(real_buf_pt,adata_pt+kount,nmove);
//        start press_all_blanks
          if(strncmp(real_buf_pt,&blank[0],blank_length)==0)
          {
            p[tp_kount] = t_p_default;
            input_end = TRUE;
//           exit convert_pressures
            break;
          }
          else
          {
            p[tp_kount] = strtod(real_buf_pt,&end_ptr);
//            start pressure_find_blanks
            for(m=blank_index;m>=0;m=m-1)
            {
//             start pressure_found_nonblank
               if(*(real_buf_pt+m) != iblank)
               {
//               exit pressure_find_blanks
                 break;
               } // end pressure_found_nonblank
            } // end pressure_find_blanks
//          start pressure_error
            if(end_ptr != (real_buf_pt+(m+1)))
            {
              fprintf(output_unit,
                " Error in PRESSURE Value number %d\n",(tp_kount+1));
              cout << " Error in PRESSURE Value number "
                   << (tp_kount+1)<<endl;
              *read_error_pt = TRUE;
//            exit convert_pressures
              break;
            } // end pressure_error
          } // end press_all_blanks
        } // end pressure_default_test
        tp_kount++;
        kount = kount + field_width;
        nbytes_count = nbytes_count + field_width;
      } // end convert_pressures
//    start pressure_input_end
      if(input_end)
      {
//      exit read_pressure values
        break;
      } // end pressure_input_end
//    start exit_read_pressure_values
      if(*read_error_pt)
      {
//      exit read_pressure values
        break;
      } // end exit_read_pressure_values
    } // end read_pressure_values
//      start exit_do_read
       if(*read_error_pt)
    {
//    exit do_while_loop
      *break_value_pt = TRUE;
      break;
    } // end exit_do_read
 
// output the pressure values
 
    pressure_number = tp_kount;
 
//  start first_p_index
    if(tp_kount < max_out_values)
    {
      tp_index = tp_kount;
    }
    else
    {
      tp_index = max_out_values;
    } // end first_p_index
//  start pressure_output_first
    for (i=0;i<tp_index;i++)
    {
//    start pressure_output_label
      if(i == 0)
      {
        fprintf(output_unit,"\n P     = %10.1f",p[i]);
      }
      else
      {
        fprintf(output_unit,"%10.1f",p[i]);
      } // end pressure_output_label
    } // end pressure_output_first
//  start pressure_output_rest
    for (i=max_out_values;i<max_t_p;i = i + max_out_values)
    {
//    start pressure_output_others
      for (j=i;j<i+max_out_values;j++)
      {
//      start pressure_output_others_first
        if(j  >= tp_kount)
        {
//        exit pressure_output_others
          break;
        } // end pressure_output_others_first
//      start pressure_output_others_second
        if(j == i)
        {
          fprintf(output_unit,"\n         %10.1f",p[j]);
        }
        else
        {
//        start output_pressure_partial
          if(j < max_t_p)
          {
            fprintf(output_unit,"%10.1f",p[j]);
          }// end output_pressure_partial
        } // end pressure_output_others_second
      } // end pressure_output_others
//    start pressure_output_others_third
      if(j  >= tp_kount)
      {
//      exit pressure_output_rest
        break;
      } // end pressure_output_others_third
 
    } // end pressure_output_rest
//  start pressure_check
    for(pressure_index=0;pressure_index<pressure_number;pressure_index++)
    {
//    start negative_pressure
      if(p[pressure_index] <= 0.0)
      {
        fprintf(output_unit,
         "\n Input Pressure %6.1f below or equal zero\n",p[pressure_index]);
        cout <<
         " Input Pressure " << p[pressure_index] << " below or equal zero"<<endl;
        *read_error_pt = TRUE;
//      exit pressure_check
        break;
      } // end negative_pressure
    } // end pressure_check
//  start pressure_exit
    if(*read_error_pt)
    {
//    exit do_while_loop
      *break_value_pt = TRUE;
      break;
    }// end pressure_exit
 
    fprintf(output_unit,"\n");
    fprintf(output_unit,"\n");
 
//  start kelvin_convert
    for( temp_index = 0;temp_index < temp_number ;temp_index++)
    {
//    start temps_present
      if(t[temp_index] <  -(kelvin + round_up))
      {
        fprintf(output_unit,
         "\n Error - Input Temperature %6.1f below zero Kelvin\n",t[temp_index]);
        cout <<
         " Error - Input Temperature " << t[temp_index]
         <<" below zero Kelvin"<<endl;
        *read_error_pt = TRUE;
//      exit kelvin_convert
        break;
      }
      else if (t[temp_index] > (max_temp - kelvin))
      {
        fprintf(output_unit,
         "\n Error - Input Temperature %6.1f above %7.1f degrees Kelvin\n",
              t[temp_index],max_temp);
        cout <<
         " Error - Input Temperature "<<t[temp_index]
        <<" above << max_temp << degrees Kelvin"<<endl;
        *read_error_pt = TRUE;
//      exit kelvin_convert
        break;
      }
      else if(t[temp_index] == t_p_default)
      {
//      exit kelvin_convert
        break;
      }
      else
      {
        t[temp_index] = t[temp_index] + kelvin ;
      } // end temps_present
    } // end kelvin_convert
//  start temperature_exit
    if(*read_error_pt)
    {
//      exit do_while_loop
      *break_value_pt = TRUE;
      break;
    }// end temperature_exit
 
//  exit do_while_loop
    break;
 
  } // end do_while_loop
 
} // end function process_pressure
//###
//*****************************************************************************
//                                                                            *
//     Function Name: react                                                   *
//                                                                            *
//     File Name: react.cpp                                                   *
//                                                                            *
//     Purpose:  This function is called from therm and molewt.  When it is   *
//               called from therm, it reads the reactant records and stores  *
//               the element symbols in the llmt array.  This function also   *
//               stores the number of elements in nlm and the number of       *
//               reactant records in nreac.  It then compares the mole        *
//               fraction of each element.  It also determines what type of   *
//               mixture problems is being considered (see molewt).  React    *
//               will store the number of elements to be stored for later use *
//               in nmin if a plus or systm problem is being considered.  It  *
//               is called from molewt every time a new base value is computed*
//               and only computes the new mole fractions.  React will reset  *
//               nmin if necessary.                                           *
//                                                                            *
//     Author:    Robert L. Reese                                             *
//                                                                            *
//     Date: June 01, 1973                                                    *
//                                                                            *
//     Language: Borland C++ Compiler, Version 4.5                            *
//                                                                            *
//     Hardware: IBM Compatible PC                                            *
//                                                                            *
//     Operating System: Microsoft Windows 98                                 *
//                                                                            *
//     Update History:                                                        *
//                                                                            *
//     Name        Date      Revision          Changes                        *
//     ----        ----      ---------         -------                        *
//                                                                            *
//    R.L. Reese   06/01/73      0             Initial Release                *
//    R.L. Reese   11/23/03      A             Converted to C++               *
//                                                                            *
//----------------------------------------------------------------------------*
//                                                                            *
//  Usage:                                                                    *
//                                                                            *
//     react(read_status,                                                     *
//           output_file,                                                     *
//           input_file) where:                                               *
//                                                                            *
//     read_status - pointer to integer value of status of read               *
//     output_file - FILE assignment for output listings                      *
//     input_file  - FILE assignment for input data file                      *
//                                                                            *
//                                                                            *
//----------------------------------------------------------------------------*
//                                                                            *
//     Definitions:                                                           *
//                                                                            *
//     Local Variables:                                                       *
//     ----------------                                                       *
//                                                                            *
//     adata        - Real array used as summation of molecular weight of each*
//                    element in the reactant records                         *
//     atom_temp    - Character containing temporay storage for current atom  *
//     break_flag   - Logical variable used to determine if function needs to *
//                    break from outer loop on return from function           *
//     empty_file   - Integer file status flag returned from init_read        *
//     exit_loop    - Logical variable used to exit do loop                   *
//     go_back      - Logical variable used to loop back in the code          *
//     go_to_rm     - Logical variable used to skip to RM calculation         *
//     rdata        - Character buffer containg data returned from init_read  *
//     j            - Integer variable used as do loop index                  *
//     jj           - Integer variable used as number of elements being       *
//                    considered in a reactant record                         *
//     kk           - Integer variable used as do loop index                  *
//     kl           - Integer variable used as do loop index                  *
//     kt           - Integer variable used to save kk do loop index          *
//     loop_again   - Logical variable used to loop back again                *
//     loop_end     - Logical variable used to go to end of loop              *
//     n            - Integer variable used as number of reactant records     *
//                    being considered                                        *
//     nbytes       - Integer variable containg number of bytes from init_read*
//     nj           - Integer variable used as counter for number of elements *
//                    in reactant records                                     *
//     pos          - Fpos_t variable used to backspace input file one record *
//     read_error   - Integer logical variable used to detect conversion error*
//     rm           - Real variable used as summation of total molecular      *
//                    weight of all the elements in a reactant record         *
//     skip_code    - Logical variable used to skip code                      *
//     wp           - Real variable used as summation of total molecular      *
//                    weight of all elements being considered                 *
//                                                                            *
//     Pointers                                                               *
//     -----------------                                                      *
//                                                                            *
//     adata_pt       - Real pointer to adata array                           *
//     atom_temp_pt   - Character pointer to temporay storage for current atom*
//     break_flag_pt  - Integer logical pointer to break flag returned from   *
//                      function                                              *
//     empty_file_pt  - Integer pointer to flag returned from init_read       *
//     rdata_pt       - Character pointer to data returned from init_read     *
//     loop_again_pt  - Integer logical pointer to loop_again flag            *
//     n_pt           - Integer pointer to variable used as number of reactant*
//                      records being considered                              *
//     nbytes_pt      - Integer pointer to number of bytes read in init_read  *
//     read_error_pt  - Integer pointer to a read error                       *
//     read_status_pt - Integer pointer to  status of a read                  *
//                                                                            *
//     Global Variables:                                                      *
//     -----------------                                                      *
//                                                                            *
//     indx.h                                                                 *
//     alter        - Logical variable indicating reactant mole factions with *
//                    unequal top and base values are to be altered           *
//     newmol       - Integer variable used to determine if reactant records  *
//                    have been read, and if type of mixture problem has been *
//                    determined.                                             *
//     nlm          - Integer variable used as number of elements in the      *
//                    system                                                  *
//     nmin         - Integer variable used as number of elements that are put*
//                    into temporary storage for systm and plus problems      *
//     nreac        - Integer variable used number of reactant records read   *
//     plus         - Logical variable used to indicate if reactants are to be*
//                    added to the syste or if a ternary system is being      *
//                    considered                                              *
//     systm        - Logical variable used to denote that a ternary system is*
//                    being considered                                        *
//     misc.h                                                                 *
//     anum         - Real array containing stochiometric coefficients of     *
//                    elements on reactant records                            *
//     base         - Real array containing first value of mole fraction of   *
//                    ith reactant                                            *
//     bo           - Real array containing assigned Kg-atoms of element i per*
//                    Kg of total reactant                                    *
//     llmt         - Character*4 containing alphameric symbols for the       *
//                    elements in the system being studied                    *
//     name         - Character*4 containing alphameric symbols for the       *
//                    elements as given in the reactant records               *
//                                                                            *
//     Constants:                                                             *
//     ---------                                                              *
//                                                                            *
//     local                                                                  *
//     FALSE           - Defined value for false logical test = 0             *
//     TRUE            - Defined value for true logical test = 1              *
//     arrays.h                                                               *
//     atom_temp_len   - Integer value of length of atom_temp array           *
//     element_max     - Integer containing maximum number of elements allowed*
//     max_atoms       - Integer containing maximum number of atoms           *
//     reacts_max      - Integer containing maximum number of reactants in a  *
//                       record                                               *
//     rdata_len       - Integer value of length of rdata_array               *
//     constant.h                                                             *
//     iblank          - Character*1 containing a blank                       *
//     atoms.h                                                                *
//     atom            - Real array containing atomic weight of the elements  *
//     iatom           - Character*2 array containing atomic symbol of        *
//                       elements                                             *
//                                                                            *
//     Include Files:                                                         *
//     -------------                                                          *
//                                                                            *
//     arrays.h          - Contains sizes for include files                   *
//     atoms.h           - Contains atomic symbols and weights                *
//     contant.h         - Contains various constants                         *
//     indx.h            - Contains index variables                           *
//     iostream.h        - Contains I/O routines                              *
//     misc.h            - Contains miscellaneous variables                   *
//     stdio.h           - Contains I/O routines                              *
//     string.h          - Contains string routines                           *
//     stdlib.h          - Contains standard library routines                 *
//                                                                            *
//     External Software:                                                     *
//     -----------------                                                      *
//                                                                            *
//     clear_arrays- Function to initialize various arrays                    *
//     init_read   - Function to read input records and check for errors      *
//     read_react  - Reads the reactant records and processes them            *
//     stdio.h     - C++ function for I/O                                     *
//     strncmp     - C++ function to compare strings                          *
//     strncpy     - C++ function to copy strings                             *
//     strtod      - C++ routine to convert from ascii to float               *
//     system_type - Function to determine type of system being used, alter,  *
//                   systm or plus                                            *
//                                                                            *
//     Files:                                                                 *
//     -----                                                                  *
//                                                                            *
//       Input Files:                                                         *
//                                                                            *
//       input_unit - Logical unit number of file containing input records    *
//                                                                            *
//       Intermediate Files:                                                  *
//                                                                            *
//       None                                                                 *
//                                                                            *
//       Output Files:                                                        *
//                                                                            *
//       output_unit - Logical unit number of file containing output listing  *
//                                                                            *
//*****************************************************************************
//###
 
#include "arrays.h"
#include "indx.h"
#include <iostream.h>
#include "misc.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
 
#define FALSE 0
#define TRUE 1
 
void clear_arrays(float *adata);
 
void init_read(char *rdata,
               int *eof,
               FILE *opu,
               FILE *inp_unit,
               int *empty_file,
               int *nbytes);
 
void read_react(FILE *output_unit,
                FILE *input_unit,
                int  *read_status_pt,
                int  *read_error_pt,
                int  *n_pt,
                int  *loop_again_pt,
                int  *break_flag_pt);
 
void system_type(FILE *output_unit,
                 int *read_status_pt);
 
void react(int *read_status_pt,
           FILE *output_unit,
           FILE *input_unit)
{
 
#include "constant.h"
 
#include "atoms.h"
 
  float adata[element_max] ;
 
  float rm ;
  float wp ;
 
  int break_flag;
  int empty_file;
  int j ;
  int jj ;
  int kk ;
  int kl ;
  int kt ;
  int n ;
  int nbytes;
  int nj ;
  int read_error;
 
//  pointers
 
  float *adata_pt;
  char  *atom_temp_pt;
  int   *break_flag_pt;
  int   *empty_file_pt;
  char  *rdata_pt;
  int   *loop_again_pt;
  int   *n_pt;
  int   *nbytes_pt;
  int   *read_error_pt;
 
//  logical variables
 
  int exit_loop ;
  int go_back ;
  int go_to_rm ;
  int loop_again ;
  int loop_end ;
  int skip_code ;
 
  char rdata[rdata_len];
  char atom_temp[atom_temp_len];
 
  fpos_t pos;
 
  strncpy(atom_temp,ablank,atom_temp_len);
 
// set pointers
 
  adata_pt       = &adata[0];
  break_flag_pt  = &break_flag;
  empty_file_pt  = &empty_file;
  rdata_pt       = &rdata[0];
  loop_again_pt  = &loop_again;
  n_pt           = &n;
  nbytes_pt      = &nbytes;
  atom_temp_pt   = &atom_temp[0];
  read_error_pt  = &read_error;
  read_error = FALSE;
  wp = 0.0 ;
 
// Clear arrays
 
  clear_arrays(adata_pt);
 
  n = 0 ;
  nlm = 1 ;
 
//    read and write reactant records
 
//  start outer_loop
  for (;;)
  {
    skip_code = FALSE ;
    loop_again = FALSE ;
 
//  start react_1
    if(newmol != 0)
    {
      skip_code = TRUE ;
//    start blank_test
      if(name[0][n][0] == iblank)
      {
//      exit outer_loop
        break;
      } // end blank_test
 
    }
    else
    {
      break_flag = FALSE;
 
//  read and process each reactant record
      read_react(output_unit,
                 input_unit,
                 read_status_pt,
                 read_error_pt,
                 n_pt,
                 loop_again_pt,
                 break_flag_pt);
//    start break_test
      if(br