import java.io.*;
import java.util.*;
import java.text.*;
public class combustion1
{  //Combustion of C,H,S,N base fuels
public String fuel_name[];
public String air_name[];
public String product_name[];
public String reactant_name[];
public double Nfuel[]; // fuel amount
public double Nair[];  // air amount
public double Nreactant[];
public double Nproduct[];//product gases amount
boolean isO2=false;
boolean isN2=false;
Gas fuel[];
Gas air[];
Gas product[];
Gas reactant[];
String name[];
double N[][];
int natom,nfatom,ngas,nfuel,nproduct;
double nair,nair_stoichiometric;
int nrt,npt;
String reacName;
Atom atomList[];
Atom fatomList[];
double Csay;
double Nsay;
double Osay;
double Hsay;
double Ssay;
public combustion(String crname,String fueli[],double Nfueli[],String airi[],double Nairi[],double a_f_ratio) 
{reacName=crname;

 setfuel(fueli,Nfueli);
 setcombustion(a_f_ratio);
} 
public void setfuel(String fueli[],double Nfueli[])
{  nfuel=fueli.length;
  fuel_name=new String[nfuel];
  Nfuel=new double[nfuel];
  fuel=new Gas[nfuel];
  natom=0;
  nfatom=0;
  nrt=0;
  npt=0;
for(int i=0;i<nfuel;i++)
{  	fuel_name[i]=fueli[i];
    fuel[i]=new Gas(fuel_name[i]);
}
for(int i=0;i<nfuel;i++)
{Nfuel[i]=Nfueli[i];nrt+=Nfuel[i];
}
arrange_fuel_atoms();
//arrange fuel output
Csay=0;
Nsay=0;
Osay=0;
Hsay=0;
Ssay=0;
boolean sulphur_exist=false;

for(int i=0;i<fatomList.length;i++)
{ if(fatomList[i].symbol.equals("C")) Csay=fatomList[i].N;
  else if(fatomList[i].symbol.equals("H")) Hsay=fatomList[i].N;
  else if(fatomList[i].symbol.equals("S")) {Ssay=fatomList[i].N;sulphur_exist=true;}
  else if(fatomList[i].symbol.equals("O")) Osay=fatomList[i].N;
  else if(fatomList[i].symbol.equals("N")) Nsay=fatomList[i].N;
}

}
public void setair(String airi[],double Nairi[])
{ String air_name[]=airi;
 Nair=Nairi;
}
public void setcombustion(double a_f_ratio) 
{ 
 double xO2=0;
 double xN2=0;
  xO2=Csay+Ssay+Hsay/4.0-Osay/2.0;//ideal combustion O2 amount
  // stoichiometric air
  nair_stoichiometric=xO2*(1.0+79.0/21.0);
  nair=nair_stoichiometric*a_f_ratio;
  double xO2real=xO2*a_f_ratio;  
  double xO2excess;
  if(a_f_ratio>1) xO2excess=xO2*(a_f_ratio-1.0);
  else	          xO2excess=0.0; 
  air=new Gas[airName.length];
  nrt=0;
  for(int i=0;i<airName.length;i++)
  {  air[i]=new Gas(air_name[]);
     nrt+=Nair[i];
  }
  String s[]={"co2","h2o","n2","o2","co","so2"};
  product_name=s;
  nproduct=6;
  Nproduct=new double[nproduct];
  product=new Gas[nproduct];
  for(int i=0;i<nproduct;i++)
  {product[i]=new Gas(s[i]);}
  Nproduct[5]=Ssay;//SO2
  Nproduct[1]=Hsay/2.0;//H2O
  Nproduct[2]=Nsay/2.0+Nair[1];//N2
  Nproduct[3]=xO2excess;//O2
  double xco=0;
  double xco2=0;
  if(a_f_ratio>=1) xco2=Csay;
  else {xco2=Osay+2*xO2real-Csay-Nproduct[1]-2.0*Nproduct[3]-2.0*Nproduct[5];xco=Csay-xco2;}
  Nproduct[0]=xco2;//CO2
  Nproduct[4]=xco;//CO
  int i;
  for(i=0;i<Nproduct.length;i++) npt+=Nproduct[i];
  for(i=0;i<Nfuel.length;i++) 
  {if(fuel_name.equals("o2")||fuel_name.equals("O2"))isO2=true;
  else if(fuel_name.equals("n2")||fuel_name.equals("N2"))isN2=true;
  }
  if(isO2 && isN2) 
  {Nreactant=new double[nfuel];
   reactant_name=new String[nfuel];
   reactant=new Gas[nfuel];
  for(i=0;i<Nfuel.length;i++)
    {if(fuel_name.equals("o2")||fuel_name.equals("O2")) {Nreactant[i]=Nfuel[i]+Nair[0];}
     else if(fuel_name.equals("n2")||fuel_name.equals("N2")) {Nreactant[i]=Nfuel[i]+Nair[1];}
     else {Nreactant[i]=Nfuel[i];}
     reactant_name[i]=fuel_name[i];
     reactant[i]=new Gas(reactant_name[i]);
    }    
  }
  else if(isO2) 
  {Nreactant=new double[nfuel+1];
   reactant_name=new String[nfuel+1];
   reactant=new Gas[nfuel+1];
  for(i=0;i<Nfuel.length;i++)
    {if(fuel_name.equals("o2")||fuel_name.equals("O2")) {Nreactant[i]=Nfuel[i]+Nair[0];}
     else {Nreactant[i]=Nfuel[i];}
     reactant_name[i]=fuel_name[i];
     reactant[i]=new Gas(reactant_name[i]);
    }
  Nreactant[i]=Nair[1];
  reactant_name[i]="n2"; 
  reactant[i]=new Gas(reactant_name[i]); 
  }
  else if(isN2) 
  {Nreactant=new double[nfuel+1];
   reactant_name=new String[nfuel+1];
   reactant=new Gas[nfuel+1];
  for(i=0;i<Nfuel.length;i++)
    {if(fuel_name.equals("n2")||fuel_name.equals("n2")) {Nreactant[i]=Nfuel[i]+Nair[1];}
     else {Nreactant[i]=Nfuel[i];}
     reactant_name[i]=fuel_name[i];
     reactant[i]=new Gas(reactant_name[i]);
    }
  Nreactant[i]=Nair[0];
  reactant_name[i]="o2"; 
  reactant[i]=new Gas(reactant_name[i]); 
  }
  else
  {Nreactant=new double[nfuel];
   reactant_name=new String[nfuel];
   reactant=new Gas[nfuel];
  for(i=0;i<Nfuel.length;i++)
    {Nreactant[i]=Nfuel[i];
     reactant_name[i]=fuel_name[i];
     reactant[i]=new Gas(reactant_name[i]);
     }  
   }
}	
     public void arrange_fuel_atoms() 
     {
     // atomic structure of the reactants
     int i,j;
        for(i=0;i<nfuel;i++)
        { for(j=0;j<fuel[i].natom;j++)
	         { add_fuel_atom(i,j);}
	    }
     }

public static String toString(double left, int w, int d)
// converts a double to a string with given width and decimals.
{
	NumberFormat df=NumberFormat.getInstance(Locale.US);
	df.setMaximumFractionDigits(d);
    df.setMinimumFractionDigits(d);
    df.setGroupingUsed(false);
    String s = df.format(left);
    while (s.length() < w)
      s = " " + s;
    if (s.length() > w)
    {
        s = "";
        for (int i=0; i<w; i++)
          s = s + "-";
    }
    return s;
}

     public int  add_fuel_atom(int i,int j) 
     {
     //determine and order the atomic structure
     int k;
       for(k=0;k<nfatom;k++)
       {
       if(fuel[i].atomList[j].symbol.equals(fatomList[k].symbol))
         {
         fatomList[k]=new Atom(fatomList[k].symbol,fatomList[k].N+fuel[i].atomList[j].N*Nfuel[i]);
         return 1;
         }
       }
     Atom atomL[];
     atomL=new Atom[nfatom+1];
     for(k=0;k<nfatom;k++)
       atomL[k]=new Atom(fatomList[k]);
     atomL[nfatom]=new Atom(fuel[i].atomList[j].symbol,fuel[i].atomList[j].N*Nfuel[i]);
     fatomList=atomL;
     nfatom++;
     return 2;
     }

     public String fuel_toString()
     { String s="";
	 for(int i=0;i<nfatom;i++)
        s+=fatomList[i].toString()+" ";  
     return s;  
	 }
     public String toString()
     {
       // writes chemical reaction in chemistry norm
       String s="[";
       for(int i=0;i<nfuel;i++)
       {
       if(Nfuel[i]==1)
           { s+=fuel[i].toString(); }
       else if(Nfuel[i]!=0)
           { s+=toString(Nfuel[i],8,4)+" "+fuel[i].toString(); }
       if(i!=(nfuel-1) && Nfuel[i]!=0 ) s+=" + ";
       else s+="] + ";
       }
       s+=toString(Nair[0],8,4)+"*[ "+air[0].toString()+" + (79/21) *"+air[1].toString()+"] --> ";
       for(int i=0;i<nproduct;i++)
       {
       if(Nproduct[i]==1)
           { s+=product[i].toString(); }
       else if(Nproduct[i]!=0)
           { s+=toString(Nproduct[i],8,4)+" "+product[i].toString(); }
       if(i!=(nproduct-1) && Nproduct[i]!=0 ) s+=" + ";
       }
       if(s.endsWith(" + "))
         s=s.substring(0,s.length()-3);
     return s;
     }
 
     public String toString(String ch)
     {
     //return the c
     String s="";
     int i,j;
     if(ch.equals("name"))
        s=s+reacName+"\n";
     else if(ch.equals("formula"))
        {
        s=toString();
        }
     else if(ch.equals("composition"))
        {
        for(i=0;i<natom;i++)
        s=s+atomList[i].toString()+" ";
        }
     return s;
     }    
     public double H(double TF,double TA,double TP)
     {
     // reaction enthalpy
      double HH=0;
      double h1,h2,h3;
      for(int i=0;i<nfuel;i++)
      {
       h1=fuel[i].ht(TF)*Nfuel[i];
       HH-=h1;
      }     
      for(int i=0;i<2;i++)
      {
       h2=air[i].ht(TA)*Nair[i];
       HH-=h2;
      }
      for(int i=0;i<nproduct;i++)
      {
       h3=product[i].ht(TP)*Nproduct[i];
       HH+=h3;
      } 
      return HH;
     }
     public double S(double TF,double TA,double TP,double pt)
     {
     //reaction entropy
      double s1,s2,s3;
      double SS=0;
      for(int i=0;i<nfuel;i++)
      {
       if(Nfuel[i]!=0)
          SS-=fuel[i].si(TF,Nfuel[i]/npt*pt)*Nfuel[i];
      }      
      for(int i=0;i<nair;i++)
      {
       if(Nair[i]!=0)
          SS-=air[i].si(TA,Nair[i]/nrt*pt)*Nair[i];
      }
      for(int i=0;i<nproduct;i++)
      {
       if(Nproduct[i]!=0)
          SS+=product[i].si(TP,Nproduct[i]/nrt*pt)*Nproduct[i];
      }
        return SS;
     }
     public double S(double TF,double TA,double TP)
     {
     //reaction entropy
      double pt=1.0;
      return S(TF,TA,TP,pt);
     }
     public double G(double TF,double TA,double TP,double pt)
     { //reaction gibbs free energy
       double g1,g2,g3;
      double GG=0;
      for(int i=0;i<nfuel;i++)
      {
       if(Nfuel[i]!=0)
          GG-=fuel[i].gt(TF,Nfuel[i]/npt*pt)*Nfuel[i];
      }      
      for(int i=0;i<nair;i++)
      {
       if(Nair[i]!=0)
          GG-=air[i].gt(TA,Nair[i]/nrt*pt)*Nair[i];
      }
      for(int i=0;i<nproduct;i++)
      {
       if(Nproduct[i]!=0)
          GG+=product[i].gt(TP,Nproduct[i]/nrt*pt)*Nproduct[i];
      }
        return GG;
     }
     
     public double G(double TF,double TA,double TP)
     {
     // reaction gibbs energy
      double pt=1.0;
        return G(TF,TA,TP,pt);
     }
     
     public double G0(double TF,double TA,double TP)
     { //reaction gibbs free energy
       double g1,g2,g3;
      double GG=0;
      for(int i=0;i<nfuel;i++)
      {
       if(Nfuel[i]!=0)
          GG-=fuel[i].gt(TF)*Nfuel[i];
      }      
      for(int i=0;i<nair;i++)
      {
       if(Nair[i]!=0)
          GG-=air[i].gt(TA)*Nair[i];
      }
      for(int i=0;i<nproduct;i++)
      {
       if(Nproduct[i]!=0)
          GG+=product[i].gt(TP)*Nproduct[i];
      }
        return GG;
     }
     public double K(double TF,double TA,double TP)
     {
     // Equilibrium costant. Usually TF,TA and TP will be in
     // the same temperature
     return Math.exp(G0(TF,TA,TP)/(-8.3145*TP));
     }
     public double funcH(double a_f_ratio,double TF,double TA,double TP)
     {setcombustion(a_f_ratio);
      return H(TF,TA,TP);
     }
     public double adiabatic_a_f_ratio(double TF,double TA,double TP,double xl,double xu)
     { //Adiabatic flame temperature
     // Adyabatic flame temperature using Bisection method
     // bisection method to find roots of zero energy change
     // Q=0
     // defination of variables :
     // xl : lower guess
     // xu : upper guess
     // xr : root estimate
     // es : stopping criterion
     // ea :approximate error
     // maxit : maximum iterations
     // iter  : number of iteration
     double test;
     double xr=0;
     double es,ea;
     double fxl=0,fxr=0;
     int maxit=500,iter=0;
     es=0.00000001;
     ea=1.1*es;
     while((ea>es)&&(iter<maxit))
     {
     xr=(xl+xu)/2.0;
     iter++;
     if((xl+xu)!=0)
        { ea=Math.abs((xu-xl)/(xl+xu))*100;}
     fxl= funcH(xl,TF,TA,TP);
     fxr= funcH(xr,TF,TA,TP);
     test= fxl*fxr;
     if(test==0.0)       ea=0;
     else if(test<0.0)  xu=xr;
     else
       {
       xl=xr;
       }
     } //end of while
        return xr;
     }
     
     public double Taf(double TF,double TA,double xl,double xu)
     { //Adiabatic flame temperature
     // Adyabatic flame temperature using Bisection method
     // bisection method to find roots of zero energy change
     // Q=0
     // defination of variables :
     // xl : lower guess
     // xu : upper guess
     // xr : root estimate
     // es : stopping criterion
     // ea :approximate error
     // maxit : maximum iterations
     // iter  : number of iteration
     double test;
     double xr=0;
     double es,ea;
     double fxl=0,fxr=0;
     int maxit=500,iter=0;
     es=0.000001;
     ea=1.1*es;
     while((ea>es)&&(iter<maxit))
     {
     xr=(xl+xu)/2.0;
     iter++;
     if((xl+xu)!=0)
        { ea=Math.abs((xu-xl)/(xl+xu))*100;}
     fxl= H(TF,TA,xl);
     fxr= H(TF,TA,xr);
     test= fxl*fxr;
     if(test==0.0)       ea=0;
     else if(test<0.0)  xu=xr;
     else
       {
       xl=xr;
       }
     } //end of while
        return xr;
     }
     
     double Vi(double TF,double TA,double TP,int Ni)
     {return -G(TF,TA,TP)/(96487*Ni);}
	 
	 double Vi0(double TF,double TA,double TP,int Ni)
     {return -G0(TF,TA,TP)/(96487*Ni);}     
}
