/*program cohsmix.c
 *
 *calculates equilibrium fractions of gas species in CO-CO2-SO2-H2S-H2O
 *(or subset) gas mixtures at temperature.
 *
 *Victor Kress Jan. 1994
 *
 *CAUTION:  Most arrays and vector are zero-offset (indexed from 0-n).
 *Arrays passed to Numerical Recipes routines are unit ofset to conform
 *to their conventions.
 *
 *on DOS machines, compile with large memory model
 *link with nrfiles.obj
 */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define NCOMP 4
#define NSPEC 31
#define NINP 6
#define MAXITER 5000
#define R 8.31468
#define SMALL 1.e-25
#define CONVERGE_TOL 1.e-10

/**********function prototypes*************************/
int speciate(int ncomp,double *spec,
                double *Gf,double tk,int nspec,double **spstoich);
double Gs(double tk);
double Gg(double tk);
double O2g(double tk);
double S2g(double tk);
double H2g(double tk);
double Sg(double tk);
double S3g(double tk);
double COg(double tk);
double CO2g(double tk);
double SO2g(double tk);
double SOg(double tk);
double S2Og(double tk);
double COSg(double tk);
double CSg(double tk);
double CS2g(double tk);
double Hg(double tk);
double OHg(double tk);
double HO2g(double tk);
double H2Og(double tk);
double H2O2g(double tk);
double H2SO4g(double tk);
double H2Sg(double tk);
double CHg(double tk);
double COHg(double tk);
double CH2g(double tk);
double COH2g(double tk);
double CH3g(double tk);
double CH4g(double tk);
double C2Hg(double tk);
double C2H2g(double tk);
double C2H4g(double tk);
double C2OH4g(double tk);
void  svbksb(double **u, double *w, double **v, int m, int n, double *b,double *x);
void  svdcmp(double **a, int m, int n, double *w, double **v);
double *dvector(int,int);
double **dmatrix(int,int,int,int);
int *ivector(int,int);
int **imatrix(int,int,int,int);
void free_dvector(double *,int,int);
void free_ivector(int *,int,int);
void free_dmatrix(double **,int,int,int,int);
void free_imatrix(int **,int,int,int,int);
void nrerror(char *);

/*global veriables and initialized arrays*/
char *splab[]={"C"   ,"O2"  ,"S2"   ,"H2",
               "S"   ,"S3"  ,"CO"   ,"CO2",
               "SO2" ,"SO"  ,"S2O"  ,"COS",
               "CS"  ,"CS2" ,"H"    ,"OH",
               "HO2" ,"H2O" ,"H2O2" ,"H2SO4",
               "H2S" ,"CH"  ,"COH"  ,"CH2",
               "COH2","CH3" ,"CH4"  ,"C2H",
               "C2H2","C2H4","C2OH4"};
int inputdex[NINP]={1,3,6,7,8,17},
    nsup=1,suppress[]={0};/*list of indices of species to suppress*/
double
         spvec[]={
         1.0,0.0,0.0,0.0, 0.0,1.0,0.0,0.0, 0.0,0.0,1.0,0.0, 0.0,0.0,0.0,1.0,
         0.0,0.0,0.5,0.0, 0.0,0.0,1.5,0.0, 1.0,0.5,0.0,0.0, 1.0,1.0,0.0,0.0,
         0.0,1.0,0.5,0.0, 0.0,0.5,0.5,0.0, 0.0,0.5,1.0,0.0, 1.0,0.5,0.5,0.0,
         1.0,0.0,0.5,0.0, 1.0,0.0,1.0,0.0, 0.0,0.0,0.0,0.5, 0.0,0.5,0.0,0.5,
         0.0,1.0,0.0,0.5, 0.0,0.5,0.0,1.0, 0.0,1.0,0.0,1.0, 0.0,2.0,0.5,1.0,
         0.0,0.0,0.5,1.0, 1.0,0.0,0.0,0.5, 1.0,0.5,0.0,0.5, 1.0,0.0,0.0,1.0,
         1.0,0.5,0.0,1.0, 1.0,0.0,0.0,1.5, 1.0,0.0,0.0,2.0, 2.0,0.0,0.0,0.5,
         2.0,0.0,0.0,1.0, 2.0,0.0,0.0,2.0, 2.0,0.5,0.0,2.0
         },
         (*energies[])(double)={Gg,   O2g,  S2g,   H2g,
                                Sg,   S3g,  COg,   CO2g,
                                SO2g, SOg,  S2Og,  COSg,
                                CSg,  CS2g, Hg,    OHg,
                                HO2g, H2Og, H2O2g, H2SO4g,
                                H2Sg, CHg,  COHg,  CH2g,
                                COH2g,CH3g, CH4g,  C2Hg,
                                C2H2g,C2H4g,C2OH4g};

main(){
char  buffer[273];
int i,j,k,nzcomp,compdex[NCOMP],nzspec,specdex[NSPEC],iO2,iCO2;
double **spstoich,**master_spstoich,tk,sum,bigspec,
         ingas[NINP],comp[NCOMP],spec[NSPEC],Gf[NSPEC];

printf(
  "\nProgram COHSmix\nby Victor Kress\nUniversity of Washington\nJan. 1994\n\n");

/*arange pointers and allocate space*/
master_spstoich=(double **)calloc((size_t) NSPEC,sizeof(double *));
for (i=0;i<NSPEC;i++) {
  master_spstoich[i]=spvec+i*NCOMP;
  specdex[i]=i;
  }

/*input temperature*/
printf("\nInput temperature in C : ");
gets(buffer);
sscanf(buffer,"%lf",&tk);
tk += 273.15;

/*input gas mix*/
for (i=0,sum=0.0;i<NINP;i++){
  ingas[i]=0.0;
  printf("\nInput %5s :",splab[inputdex[i]]);
  gets(buffer);
  sscanf(buffer,"%lf",ingas+i);
  sum += ingas[i];
  }

/*ideal gas assumed so volume fraction equal to mole fraction*/
for (i=0;i<NINP;i++) ingas[i] /= sum;

/*print mole fractions*/
printf("\n\nInput mole fractions are:");
for (i=0;i<NINP;i++){
  printf("\n%5s %f:",splab[inputdex[i]],ingas[i]);
  }

/*convert to components*/
for(i=0;i<NCOMP;i++){
  for(j=0,comp[i]=0.;j<NINP;j++){
    comp[i] += ingas[j]*master_spstoich[inputdex[j]][i];
    }
  }

/*normalize components*/
for(i=0,sum=0.0;i<NCOMP;i++)sum+=comp[i];
for(i=0;i<NCOMP;i++)comp[i]/=sum;

/*put into species vector*/
for(i=0;i<NSPEC;i++)spec[i]=100.*SMALL;
for(i=0,j=NSPEC;i<NINP;i++) if (ingas[i]>0.0) j--;  /*j is # zero spec*/
for(i=0;i<NINP;i++){
  if(ingas[i]>0.0){
    spec[inputdex[i]]=(1.0-j*100.*SMALL)*ingas[i];
    }
  }
/*print mole fractions*/
printf("\n\nComponent mole fractions are:");
for (i=0;i<NCOMP;i++){
  printf("\n%5s %f:",splab[i],comp[i]);
  }

/*Cut out suppressed species*/
for (i=nsup-1,nzspec=NSPEC;i>=0;i--){
  for(j=suppress[i];j<nzspec;j++)specdex[j]=specdex[j+1];
  nzspec--;
  }

/*Cut out zero components and species*/
  /*first setup index lists*/
for(i=0,nzcomp=0;i<NCOMP;i++){
  if (comp[i]==0.0) {
    for(j=0;j<nzspec;j++){
      if(master_spstoich[specdex[j]][i]>0.0){ /*species includes zero component*/
        for(k=j;k<nzspec;k++)specdex[k]=specdex[k+1];/*delete species from list*/
        nzspec--;
        j--;                           /*to null j++ at end of loop*/
        }
      }
    }
  else{
    compdex[nzcomp++]=i;/*add to component index*/
    }
  }

  /*shuffle components and species*/
for(i=0;i<nzcomp;i++)comp[i]=comp[compdex[i]];
for(i=0;i<nzspec;i++)spec[i]=spec[specdex[i]];

  /*fill spstoich*/
spstoich=dmatrix(0,nzspec-1,0,nzcomp-1);
for (i=0;i<nzspec;i++){
  for (j=0;j<nzcomp;j++){
    spstoich[i][j]=master_spstoich[specdex[i]][compdex[j]];
    }
  }

  /*fill Gf*/
for(i=0;i<nzspec;i++){
  Gf[i]=(*energies[specdex[i]])(tk);
  }

printf("\n");

speciate(nzcomp,spec,Gf,tk,nzspec,spstoich);


/*see if graphite is supersaturated*/
iO2=iCO2= -1;
for(i=0;i<nzspec;i++){
  if(specdex[i]==1)iO2=i;
  if(specdex[i]==7)iCO2=i;
  }
if (iO2>=0&&iCO2>=0){
  if ((Gs(tk)+O2g(tk)-CO2g(tk)+R*tk*log(spec[iO2]/spec[iCO2]))<0.0){
    printf("\n\nGraphite is supersaturated.  Results metastable.");
    }
  }

/*print results*/
for(i=0,bigspec=0.0;i<nzspec;i++)if(bigspec<spec[i])bigspec=spec[i];
bigspec /= 60.; /*bigspec now is 1/60th largest spec*/

printf("\n\nEquilibrium species (fractions and log10):");
for(i=0;i<nzspec;i++){
  printf("\n%5s = %4.2lf %4.2f ",splab[specdex[i]],spec[i],log10(spec[i]));
  for (sum=bigspec;sum<=spec[i];sum+=bigspec)printf("*");
  }

exit(0);
}

/****************************************************************/
#define T0 1000.0
/*all values from JANAF 1985 in J at 1000.0K*/

double Gs(double tk){
/*graphite solid*/
  double Gf;
  Gf= 0.0 - 24.457*(tk-T0) + 21.610*(tk-T0-tk*log(tk/T0));
  return(Gf);
  }
double Gg(double tk){
/*carbon gas*/
  double Gf;
  Gf= 560654.0 - 183.278*(tk-T0) + 20.791*(tk-T0-tk*log(tk/T0));
  return(Gf);
  }
double O2g(double tk){
/*O2 gas*/
  double Gf;
  Gf= 0.0 - 243.578*(tk-T0) + 34.870*(tk-T0-tk*log(tk/T0));
  return(Gf);
  }
double S2g(double tk){
/*S2 gas*/
  double Gf;
  Gf= 0.0 - 270.807*(tk-T0) + 37.277*(tk-T0-tk*log(tk/T0));
  return(Gf);
  }
double H2g(double tk){
/*H2 gas*/
  double Gf;
  Gf= 0.0 - 166.216*(tk-T0) + 30.205*(tk-T0-tk*log(tk/T0));
  return(Gf);
  }
double Sg(double tk){
/*S gas*/
  double Gf;
  Gf= 156096. - 195.142*(tk-T0) + 21.489*(tk-T0-tk*log(tk/T0));
  return(Gf);
  }
double S3g(double tk){
/*S3 gas*/
  double Gf;
  Gf= 21078.0 - 334.323*(tk-T0) + 57.482*(tk-T0-tk*log(tk/T0));
  return(Gf);
  }
double COg(double tk){
/*carbon monoxide gas*/
  double Gf;
  Gf= -200275.0 - 234.538*(tk-T0) + 33.183*(tk-T0-tk*log(tk/T0));
  return(Gf);
  }
double CO2g(double tk){
/*carbon dioxide gas*/
  double Gf;
  Gf= -395886.0 - 269.299*(tk-T0) + 54.308*(tk-T0-tk*log(tk/T0));
  return(Gf);
  }
double SO2g(double tk){
/*sulfur dioxide gas*/
  double Gf;
  Gf= -288725.0 - 305.767*(tk-T0) + 54.484*(tk-T0-tk*log(tk/T0));
  return(Gf);
  }
double SOg(double tk){
/*sulfur monoxide gas*/
  double Gf;
  Gf= -64382.0 - 262.155*(tk-T0) + 36.053*(tk-T0-tk*log(tk/T0));
  return(Gf);
  }
double S2Og(double tk){
/*disulfur monoxide gas*/
  double Gf;
  Gf= -179507.0 - 328.491*(tk-T0) + 55.656*(tk-T0-tk*log(tk/T0));
  return(Gf);
  }
double COSg(double tk){
/*carbon oxide sulfide gas*/
  double Gf;
  Gf= -212544.0 - 291.828*(tk-T0) + 56.990*(tk-T0-tk*log(tk/T0));
  return(Gf);
  }
double CSg(double tk){
/*carbon sulfide gas*/
  double Gf;
  Gf= 124839.0 - 250.094*(tk-T0) + 35.612*(tk-T0-tk*log(tk/T0));
  return(Gf);
  }
double CS2g(double tk){
/*carbon disulfide gas*/
  double Gf;
  Gf= -17306.0 - 302.345*(tk-T0) + 59.283*(tk-T0-tk*log(tk/T0));
  return(Gf);
  }
double Hg(double tk){
/*monatomic H gas*/
  double Gf;
  Gf= 165485.0 - 139.871*(tk-T0) + 20.786*(tk-T0-tk*log(tk/T0));
  return(Gf);
  }
double OHg(double tk){
/*OH gas*/
  double Gf;
  Gf= 23391.0 - 219.736*(tk-T0) + 30.676*(tk-T0-tk*log(tk/T0));
  return(Gf);
  }
double HO2g(double tk){
/*HO2 gas*/
  double Gf;
  Gf= 46783.0 - 278.535*(tk-T0) + 47.604*(tk-T0-tk*log(tk/T0));
  return(Gf);
  }
double H2Og(double tk){
/*H2O gas*/
  double Gf;
  Gf= -192590.0 - 232.738*(tk-T0) + 41.268*(tk-T0-tk*log(tk/T0));
  return(Gf);
  }
double H2O2g(double tk){
/*H2O2 gas*/
  double Gf;
  Gf= -28578.0 - 297.897*(tk-T0) + 62.844*(tk-T0-tk*log(tk/T0));
  return(Gf);
  }
double H2SO4g(double tk){
/*H2SO4 gas*/
  double Gf;
  Gf= -440854.0 - 432.545*(tk-T0) + 132.690*(tk-T0-tk*log(tk/T0));
  return(Gf);
  }
double H2Sg(double tk){
/*H2S gas*/
  double Gf;
  Gf= -40984.0 - 252.579*(tk-T0) + 45.786*(tk-T0-tk*log(tk/T0));
  return(Gf);
  }
double CHg(double tk){
/*CH gas*/
  double Gf;
  Gf= 481496.0 - 219.334*(tk-T0) + 32.565*(tk-T0-tk*log(tk/T0));
  return(Gf);
  }
double COHg(double tk){
/*COH gas*/
  double Gf;
  Gf= -4783.0 - 273.454*(tk-T0) + 48.056*(tk-T0-tk*log(tk/T0));
  return(Gf);
  }
double CH2g(double tk){
/*CH2 gas*/
  double Gf;
  Gf= 331620.0 - 241.075*(tk-T0) + 45.571*(tk-T0-tk*log(tk/T0));
  return(Gf);
  }
double COH2g(double tk){
/*COH2 gas*/
  double Gf;
  Gf= -88068.0 - 275.624*(tk-T0) + 61.995*(tk-T0-tk*log(tk/T0));
  return(Gf);
  }
double CH3g(double tk){
/*CH3 gas*/
  double Gf;
  Gf= 159817.0 - 251.521*(tk-T0) + 58.954*(tk-T0-tk*log(tk/T0));
  return(Gf);
  }
double CH4g(double tk){
/*CH4 gas*/
  double Gf;
  Gf= 19492.0 - 247.549*(tk-T0) + 71.795*(tk-T0-tk*log(tk/T0));
  return(Gf);
  }
double C2Hg(double tk){
/*C2H gas*/
  double Gf;
  Gf= 346268.0 - 260.550*(tk-T0) + 50.943*(tk-T0-tk*log(tk/T0));
  return(Gf);
  }
double C2H2g(double tk){
/*C2H2 gas*/
  double Gf;
  Gf= 315144.0 - 269.192*(tk-T0) + 68.275*(tk-T0-tk*log(tk/T0));
  return(Gf);
  }
double C2H4g(double tk){
/*C2H4 gas*/
  double Gf;
  Gf= 119122.0 - 300.408*(tk-T0) + 93.899*(tk-T0-tk*log(tk/T0));
  return(Gf);
  }
double C2OH4g(double tk){
/*C2OH4 gas*/
  double Gf;
  Gf= 94907.0 - 340.800*(tk-T0) + 114.963*(tk-T0-tk*log(tk/T0));
  return(Gf);
  }

/****************************************************************/
int speciate(int ncomp,double *spec,
              double *Gf,double tk,int nspec,double **spstoich){
/*Implements RAND speciation algorithm as presented in:
 *
 *Smith and Missen (1982) Chemical Reaction Equilibrium Analysis:
 *Theory and Algorithms. Krieger.
 *
 *Gf is vector of standard state chemical potentials at T
 *spstoich[0:nspec-1][0:ncomp-1]
 */


int i,j,k,max=MAXITER,niter,converged;
double sum,*rhs,*psi,**A,*b,*comp,omega,*delta,
       *mu,dGdw[2],**v,*w,mdel;
A=dmatrix(1,ncomp+1,1,ncomp+1);
v=dmatrix(1,ncomp+1,1,ncomp+1);
w=dvector(1,ncomp+1);
rhs=dvector(1,ncomp+1);
psi=dvector(0,ncomp);
b=dvector(0,ncomp-1);
comp=dvector(0,ncomp-1);
delta=dvector(0,nspec-1);
mu=dvector(0,nspec-1);

/*********fill b and comp vector*************************/
for(k=0;k<ncomp;k++){
  for(j=0,b[k]=0.0;j<nspec;j++){
    b[k]+=spstoich[j][k]*spec[j];
    }
  comp[k]=b[k]; /*debug*/
  }

/*****************main iteration loop*************************/
/*printf("\n");*/
for (niter=0;niter<=max;niter++) {
  /*********calculate potentials****************/
  for(i=0,sum=0.0;i<nspec;i++)sum+=spec[i];
  for(i=0;i<nspec;i++) mu[i]=Gf[i]+R*tk*log(spec[i]/sum);
  if(errno==EDOM)printf("line #%d\n",__LINE__);

  /***********fill A and rhs***(1 offset for N.R.)*****************/
  /*Equation 6.3-26*/
    /*lhs*/
  for (j=0;j<ncomp;j++){
    for (i=0;i<ncomp;i++){
      for (k=0,A[j+1][i+1]=0.0;k<nspec;k++){
        A[j+1][i+1] += spstoich[k][i]*spstoich[k][j]*spec[k];
        }
      }
    A[j+1][ncomp+1]=b[j]; /*this is the 'u' variable in Smith and Missen*/
    }
    /*rhs*/
  for (j=0;j<ncomp;j++){
    for (k=0,rhs[j+1]=0.0;k<nspec;k++)
      rhs[j+1]+=spstoich[k][j]*spec[k]*mu[k]/(R*tk);
    rhs[j+1] += comp[j]-b[j];
    }

  /*Equation 6.3-27*/
  for (i=0;i<ncomp;i++) A[ncomp+1][i+1]=b[i];
  A[ncomp+1][ncomp+1]=0.0;/*no innert species (nz=0)*/
  for (k=0,rhs[ncomp+1]=0.0;k<nspec;k++){
    rhs[ncomp+1] += spec[k]*mu[k]/(R*tk);
    }

  /*************solve for Psi*******************************/
  svdcmp(A,ncomp+1,ncomp+1,w,v);
  svbksb(A,w,v,ncomp+1,ncomp+1,rhs,psi-1);

  /*************calculate increment*(Eq. 6.3-24)*****************/
  for (j=0;j<nspec;j++){
    for (k=0,delta[j]=0.0;k<ncomp;k++)delta[j]+=spstoich[j][k]*psi[k];
    delta[j] += psi[ncomp]-mu[j]/(R*tk);
    delta[j] *= spec[j];
    }

  /*calculate factor necessary to prevent negative fractions*/
  for (j=0,omega=1.0;j<nspec;j++){
    if (delta[j]<0.&&spec[j]+delta[j]<SMALL) {
      double temp;
      temp=spec[j]*(1.e-7 - 1.)/delta[j];
      if(temp<0.){
        printf("negative value for temp at line #%d\n",__LINE__);
	}
      if (temp<omega)omega=temp;
      }
    }
  if (omega<1.0){
    for(j=0;j<nspec;j++){
      delta[j]*=omega; /*to prevent iteration to negative fractions*/
      }
    }

  /*************calculate omega**********************************/
  /*see Smith and Missen p. 115*/
  for (i=0,sum=0.0;i<nspec;i++)sum += spec[i]+delta[i];
  for (i=0,dGdw[0]=dGdw[1]=0.0;i<nspec;i++){
    dGdw[0]+=mu[i]*delta[i];
    dGdw[1]+=(Gf[i]
           +R*tk*log((spec[i]+delta[i])/sum))*delta[i];
    if(errno==EDOM)printf("line #%d\n",__LINE__);
    }
  if (dGdw[0]<0.0) {
    if (dGdw[1]<=0.0) omega=1.0;
    else omega=dGdw[0]/(dGdw[0]-dGdw[1]);
    }
  else {
    /*printf("\nwarning: speciation function not incrementing 'downhill'");*/
    omega=0.5;
    }

  /*************increment species********************************/
  for (i=0,mdel=0.0;i<nspec;i++) {
    double temp;
    if (spec[i]>10.*SMALL){  /*no tiny species in convergence test*/
      temp=fabs(delta[i])/spec[i];
      if (mdel<temp) mdel=temp;
      }
    spec[i] += omega*delta[i];
    if (spec[i]<=0.0){
      printf("\nAttempt to drive species %d negative at iteration %d",
	 i,niter);
      printf("\nSpec[%d]=%10.1e  delta=%10.1e",spec[i],delta[i]);
      spec[i]=SMALL;
      }
    }

  /*************test for convergence*****************************/
  if (mdel<CONVERGE_TOL) break;
  } /*end of main loop*/

/*****************wrap up******************************/
for (j=0,sum=0.0;j<nspec;j++) sum+=spec[j];
for (j=0;j<nspec;j++) spec[j]/=sum;
converged = (sum>0&&niter<=max);
printf("\n%s in %d iterations",(converged)? "Converged":"Not converged",niter);
printf("\nFinal iteration max change %10.1g percent.",mdel*100.);
free_dmatrix(A,1,ncomp+1,1,ncomp+1);
free_dmatrix(v,1,ncomp+1,1,ncomp+1);
free_dvector(w,1,ncomp+1);
free_dvector(rhs,1,ncomp+1);
free_dvector(psi,0,ncomp);
free_dvector(b,0,ncomp-1);
free_dvector(comp,0,ncomp-1);
free_dvector(delta,0,nspec-1);
free_dvector(mu,0,nspec-1);
return(converged);
}