//======================================================
// Thermodynamic Package in java
// Class boiler
// Dr. Turhan Coban
// Ege University, School of Engineering, 
// Department of Mechanical Engineering
// Bornova İZMİR - TURKEY
// boiler combustion efficiency

import java.io.*;
import java.net.*;
import java.util.*;
import javax.swing.*;

class boiler_gibbs
{

Gmix fuel;
Gmix air;
Gmix exhaust;
Gmix exhaustanalysis;
double Tg_fuel; // fuel inlet temperature
double Tg_air;  // air inlet temperature
double Tg_exhaust; //Exhaust gases exit temperature
double excess_air_coefficient; //coefficient of excess air
double NT; //
double Q; //boiler heat transfer
double Qmax; //Higher heating value or lower heating value
double eff; //efficiency

//Constructor to set up combustion from fuel and air information only
public boiler(String g_fuelname[],double g_fuelmol[],double Ty,double Th,double Te,double hk,boolean higher_heating_value)
{ set_boiler(g_fuelname,g_fuelmol,Ty,Th,Te,hk,higher_heating_value);}
//Constructor to set up combustion from fuel and air and exhaust analysis (mol percentage) information 
public boiler(String g_fuelname[],double g_fuelmol[],String g_exhaustname[],double g_exhaustmol[],double Ty,double Th,double Te,double hfk,boolean higher_heating_value)
{ set_boiler(g_fuelname,g_fuelmol,g_exhaustname,g_exhaustmol,Ty,Th,Te,hfk,higher_heating_value);}

public void set_boiler(String g_fuelname[],double g_fuelmol[],double Ty,double Th,double Te,double hk,boolean higher_heating_value)
{ // When fuel and air is given analyse the efficiency
  // Ty : fuel inlet temperature degree K
  // Th : air  inlet temperature degree K
  // Te : exhaust temperature degree K
  // higher_heating_value : true  : fuel efficiency will be given as a function of higher heating value
  // higher_heating_value : false : fuel efficiency will be given as a function of lower heating value
  // g_fuelname : names of the component gases in fuel
  // g_fuelmol  : moles(amount) of the component gases in fuel
int n=g_fuelname.length;
Gas g_fuel[];
Gas g_air[];
Gas g_exhaust[];
Gas g_exhaustanalysis[];
excess_air_coefficient=hk;
g_fuel=new Gas[n];
NT=0;
for(int i=0;i<n;i++)
  {NT+=g_fuelmol[i];}
for(int i=0;i<n;i++)
  {g_fuelmol[i]=g_fuelmol[i]/NT;g_fuel[i]=new Gas(g_fuelname[i],g_fuelmol[i]);}      
double Csay=0;
double Nsay=0;
double Osay=0;
double Hsay=0;
double Ssay=0;
double xO2=0;
double xN2=0;
for(int i=0;i<n;i++)
  {int m=g_fuel[i].atomList.length;
  for(int j=0;j<m;j++) 
    {if(g_fuel[i].atomList[j].symbol.equals("C")) Csay+=g_fuel[i].N*g_fuel[i].atomList[j].N;
     else if(g_fuel[i].atomList[j].symbol.equals("H")) Hsay+=g_fuel[i].N*g_fuel[i].atomList[j].N;
     else if(g_fuel[i].atomList[j].symbol.equals("O")) Osay+=g_fuel[i].N*g_fuel[i].atomList[j].N;
     else if(g_fuel[i].atomList[j].symbol.equals("S")) Ssay+=g_fuel[i].N*g_fuel[i].atomList[j].N;
     else if(g_fuel[i].atomList[j].symbol.equals("N")) Nsay+=g_fuel[i].N*g_fuel[i].atomList[j].N;
     else {System.out.println("listede olmıyan bir gaz mevcut : "+g_fuel[i].atomList[j].symbol);}
    }
  }
  xO2=Csay+Ssay+Hsay/4.0-Osay/2.0;//ideal yanma için gerekli O miktarı
  //ideal gaz çıkış atomları
  double hg_fuel=0;
  double hg_air=0;
  double hg_exhaust=0;
  for(int i=0;i<n;i++)
  {hg_fuel+=g_fuel[i].HT(298.0);}
  g_air=new Gas[2];
  g_air[0]=new Gas("o2",xO2);
  g_air[1]=new Gas("n2",(xO2*79.0/21.0));
  hg_air=g_air[0].HT(298.0)+g_air[1].HT(298.0);
  g_exhaust=new Gas[5];
  g_exhaust[0]=new Gas("co2",Csay);
  if(higher_heating_value) g_exhaust[1]=new Gas("h2ol",(Hsay/2.0));
  else g_exhaust[1]=new Gas("h2o",(Hsay/2.0));
  g_exhaust[2]=new Gas("so2",Ssay);
  g_exhaust[3]=new Gas("n2",(xO2*79.0/21.0+Nsay/2.0));
  hg_exhaust=g_exhaust[0].HT(298.0)+g_exhaust[1].HT(298.0)+g_exhaust[2].HT(298.0)+g_exhaust[3].HT(298);
  Qmax=hg_exhaust-hg_air-hg_fuel;
  double xO2gercek=xO2*excess_air_coefficient;
  hg_fuel=0;
  for(int i=0;i<n;i++)
  {hg_fuel+=g_fuel[i].HT(Ty);}
  g_air[0]=new Gas("o2",xO2gercek);
  g_air[1]=new Gas("n2",(xO2gercek*79.0/21.0));
  hg_air=g_air[0].HT(Th)+g_air[1].HT(Th);
  g_exhaust[1]=new Gas("h2o",(Hsay/2.0));
  g_exhaust[3]=new Gas("n2",(xO2gercek*79.0/21.0+Nsay/2.0));
  g_exhaust[4]=new Gas("o2",(xO2gercek-xO2));
  hg_exhaust=g_exhaust[0].HT(Te)+g_exhaust[1].HT(Te)+g_exhaust[2].HT(Te)+g_exhaust[3].HT(Te)+g_exhaust[4].HT(Te);
  Q=hg_exhaust-hg_air-hg_fuel;
  eff=Q/Qmax;
  fuel=new Gmix(g_fuel);
  air=new Gmix(g_air);
  exhaust=new Gmix(g_exhaust);
  fuel=new Gmix(g_fuel);
  air=new Gmix(g_air);
  exhaust=new Gmix(g_exhaust);
  }

public void set_boiler(String g_fuelname[],double g_fuelmol[],String g_exhaustname[],double g_exhaustmol[],double Ty,double Th,double Te,double hfk,boolean higher_heating_value)
{ // When fuel and air and exhaust analysis is given analyse the efficiency
  // Ty : fuel inlet temperature degree K
  // Th : air  inlet temperature degree K
  // Te : exhaust temperature degree K
  // higher_heating_value : true  : fuel efficiency will be given as a function of higher heating value
  // higher_heating_value : false : fuel efficiency will be given as a function of lower heating value
int n=g_fuelname.length;
Gas g_fuel[];
Gas g_air[];
Gas g_exhaust[];
Gas g_exhaustanalysis[];
g_fuel=new Gas[n];
excess_air_coefficient=hfk;
NT=0;
for(int i=0;i<n;i++)
  {NT+=g_fuelmol[i];}
for(int i=0;i<n;i++)
  {g_fuelmol[i]=g_fuelmol[i]/NT;g_fuel[i]=new Gas(g_fuelname[i],g_fuelmol[i]);}      
double Csay=0;
double Nsay=0;
double Osay=0;
double Hsay=0;
double Ssay=0;
double xO2=0;
double xN2=0;
for(int i=0;i<n;i++)
  {int m=g_fuel[i].atomList.length;
  for(int j=0;j<m;j++) 
    {if(g_fuel[i].atomList[j].symbol.equals("C")) Csay+=g_fuel[i].N*g_fuel[i].atomList[j].N;
     else if(g_fuel[i].atomList[j].symbol.equals("H")) Hsay+=g_fuel[i].N*g_fuel[i].atomList[j].N;
     else if(g_fuel[i].atomList[j].symbol.equals("O")) Osay+=g_fuel[i].N*g_fuel[i].atomList[j].N;
     else if(g_fuel[i].atomList[j].symbol.equals("S")) Ssay+=g_fuel[i].N*g_fuel[i].atomList[j].N;
     else if(g_fuel[i].atomList[j].symbol.equals("N")) Nsay+=g_fuel[i].N*g_fuel[i].atomList[j].N;
     else {System.out.println("listede olmıyan bir gaz mevcut : "+g_fuel[i].atomList[j].symbol);}
    }
  }
g_exhaustanalysis=new Gas[g_exhaustmol.length];
double NTg_exhaust=0;
double n_N2=0;  
for(int i=0;i<n;i++)
  {NTg_exhaust+=g_exhaustmol[i];}

for(int i=0;i<n;i++)
  {g_exhaustanalysis[i]=new Gas(g_exhaustname[i],g_exhaustmol[i]);} 
double xco2=0; //co2 mol percentage
double xco=0;  //co mol percentage
for(int i=0;i<g_exhaustanalysis.length;i++)
{if(g_exhaustanalysis[i].gasName.equals("co2")) xco2=g_exhaustanalysis[i].N;
 if(g_exhaustanalysis[i].gasName.equals("co"))  xco=g_exhaustanalysis[i].N;
} 
double CX=0;
double NX=0;
double OX=0;
double HX=0;
double SX=0;
int nn=g_exhaustanalysis.length;

for(int i=0;i<n;i++)
  {int m=g_exhaustanalysis[i].atomList.length;
  for(int j=0;j<m;j++) 
    {if(g_exhaustanalysis[i].atomList[j].symbol.equals("C")) CX+=g_exhaustanalysis[i].N*g_exhaustanalysis[i].atomList[j].N;
     else if(g_exhaustanalysis[i].atomList[j].symbol.equals("H")) HX+=g_exhaustanalysis[i].N*g_exhaustanalysis[i].atomList[j].N;
     else if(g_exhaustanalysis[i].atomList[j].symbol.equals("O")) OX+=g_exhaustanalysis[i].N*g_exhaustanalysis[i].atomList[j].N;
     else if(g_exhaustanalysis[i].atomList[j].symbol.equals("S")) SX+=g_exhaustanalysis[i].N*g_exhaustanalysis[i].atomList[j].N;
     else if(g_exhaustanalysis[i].atomList[j].symbol.equals("N")) NX+=g_exhaustanalysis[i].N*g_exhaustanalysis[i].atomList[j].N;
     else {System.out.println("listede olmıyan bir gaz mevcut : "+g_exhaustanalysis[i].atomList[j].symbol);}
    }
  }
double AC=CX/Csay; //g_exhaust karışımı katsayısı
n_N2=AC*(1.0-NTg_exhaust);
      
  xO2=Csay+Ssay+Hsay/4.0-Osay/2.0;//ideal yanma için gerekli O miktarı
  //ideal gaz çıkış atomları
  double hg_fuel=0;
  double hg_air=0;
  double hg_exhaust=0;
  for(int i=0;i<n;i++)
  {hg_fuel+=g_fuel[i].HT(298.0);}
  g_air=new Gas[2];
  g_air[0]=new Gas("o2",xO2);
  g_air[1]=new Gas("n2",(xO2*79.0/21.0));
  hg_air=g_air[0].HT(298.0)+g_air[1].HT(298.0);
  g_exhaust=new Gas[6];
  g_exhaust[0]=new Gas("co2",Csay*xco2/(xco2+xco));
  g_exhaust[5]=new Gas("co",Csay*xco/(xco2+xco));
  if(higher_heating_value) g_exhaust[1]=new Gas("h2ol",(Hsay/2.0));
  else g_exhaust[1]=new Gas("h2o",(Hsay/2.0));
  g_exhaust[2]=new Gas("so2",Ssay);
  g_exhaust[3]=new Gas("n2",(xO2*79.0/21.0+Nsay/2.0));
  //System.out.println("hh2o="+g_exhaust[1].HT(298.0));
  hg_exhaust=g_exhaust[0].HT(298.0)+g_exhaust[5].HT(298.0)+g_exhaust[1].HT(298.0)+g_exhaust[2].HT(298.0)+g_exhaust[3].HT(298);
  //System.out.println("co2="+Csay+"H2o="+(Hsay/2.0)+"So2"+Ssay+"n2="+(xO2*79.0/21.0+Nsay/2.0)+"hg_exhaust="+hg_exhaust);
  Qmax=hg_exhaust-hg_air-hg_fuel;
  double xO2gercek=xO2*excess_air_coefficient;
  hg_fuel=0;
  for(int i=0;i<n;i++)
  {hg_fuel+=g_fuel[i].HT(Ty);}
  g_air[0]=new Gas("o2",xO2gercek);
  g_air[1]=new Gas("n2",(xO2gercek*79.0/21.0));
  hg_air=g_air[0].HT(Th)+g_air[1].HT(Th);
  g_exhaust[1]=new Gas("h2o",(Hsay/2.0));
  g_exhaust[3]=new Gas("n2",(xO2gercek*79.0/21.0+Nsay/2.0));
  g_exhaust[4]=new Gas("o2",(xO2gercek-xO2));
  hg_exhaust=g_exhaust[0].HT(Te)+g_exhaust[5].HT(Te)+g_exhaust[1].HT(Te)+g_exhaust[2].HT(Te)+g_exhaust[3].HT(Te)+g_exhaust[4].HT(Te);
  Q=hg_exhaust-hg_air-hg_fuel;
  eff=Q/Qmax;
  fuel=new Gmix(g_fuel);
  air=new Gmix(g_air);
  exhaust=new Gmix(g_exhaust);
  exhaustanalysis=new Gmix(g_exhaustanalysis); 
  }
  
public double Qi() {return Q;}
public double Qmax() {return Qmax;}

public double V_NM3_h(double QKW)
{
double Ni=QKW/(-Qi()); //kmol/s
return Ni*8314.5*273.15/1.01325e5*3600.0; //NM^3/saat
}

public double V_M3_h(double QKW,double T,double P)
{//T derece K
 //P basınç bar
 //verilen çevre sıcaklığı ve basınçtaki gaz debisi
double Ni=QKW/(-Qi()); //kmol/s
return Ni*8314.5*T/(P*1e5)*3600.0; //NM^3/saat
}

public double Ni(double QKW)
{
double Ni=QKW/(-Qi()); //kmol/s
return Ni*3600.0; //NM^3/saat
}

public String exhaustanalysis()
{
	String s="Exhaust analysis (dry basis): \n";
for(int i=0;i<exhaustanalysis.ngas;i++)
{s+=exhaustanalysis.gasList[i].gasName+" "+exhaustanalysis.gasList[i].N+"\n";}
return s;
}

public String exhaust()
{
	String s="Exhaust : \n";
for(int i=0;i<exhaust.ngas;i++)
{s+=exhaust.gasList[i].gasName+" "+exhaust.gasList[i].N+"\n";}
return s;
}

public String fuel()
{
String s="Fuel : \n";
for(int i=0;i<fuel.ngas;i++)
{s+=fuel.gasList[i].gasName+" "+fuel.gasList[i].N+"\n";}
return s;
}

public String air()
{
String s="Air : \n ";
for(int i=0;i<air.ngas;i++)
{s+=air.gasList[i].gasName+" "+air.gasList[i].N+"\n";}
return s;
}

public double efficiency() {return eff;}

// root solving function

double fp(String g_fuelname[],double g_fuelmol[],double Ty,double Th,double Te,double hk,boolean higher_heating_value,double Qout)
{   //boiler heat is minus when heat is taken off
    //it is plus if we add heat to the boiler(it is not usual)
	//Qout should be define the same way (sign should ne minus for heat taken out of furnace)
    set_boiler(g_fuelname,g_fuelmol,Ty,Th,Te,hk,higher_heating_value);
	return Qi()-Qout;
}

public  double Ti(String g_fuelname[],double g_fuelmol[],double Ty,double Th,double hk,boolean higher_heating_value,double Qout)
{
double a=Th;
double b=3000.0;
// brent metodu ile kök değerini bulur
double test;
double p=0;
double es,ea;
double f1,f2,f3,fp;
int maxit=500,iter=0;
double tol=1.0e-5;
es=0.001;
double x1=a;f1=fp(g_fuelname,g_fuelmol,Ty,Th,x1,hk,higher_heating_value,Qout);
double x2=b;f2=fp(g_fuelname,g_fuelmol,Ty,Th,x2,hk,higher_heating_value,Qout);
double x3=(x1+x2)/2.0;f3=fp(g_fuelname,g_fuelmol,Ty,Th,x3,hk,higher_heating_value,Qout);
if(f1==0) return x1;
else if(f2==0) return x2;
else if(f3==0) return x3; 
//if(f1*f2>0) System.out.println("verilen bölgede kök yok");
p=-(f2*f3*x1*(f2-f3)+f3*f1*x2*(f3-f1)+f1*f2*x3*(f1-f2))/((f1-f2)*(f2-f3)*(f3-f1));
fp=fp(g_fuelname,g_fuelmol,Ty,Th,p,hk,higher_heating_value,Qout);
//System.out.println("p="+p+"fp="+fp);  
ea=Math.abs(x3-p);
while((ea>es)&&(iter<maxit))
{	if(Math.abs(f3)<tol) return x3; 
    if((p<x1) && (p>x2)) 
    {p=(x1+x2)/2.0;    
     if(p>x3) {x1=x3;f1=f3;x3=p;f3=fp;}
     else if(p<x3) {x2=x3;f2=f3;x3=p;f3=fp;}
    }
    else
    {        
    if(p>x3) {x1=x3;f1=f3;x3=p;f3=fp;}
    else if(p<x3) {x2=x3;f2=f3;x3=p;f3=fp;}
    p=-(f2*f3*x1*(f2-f3)+f3*f1*x2*(f3-f1)+f1*f2*x3*(f1-f2))/((f1-f2)*(f2-f3)*(f3-f1));
    fp=fp(g_fuelname,g_fuelmol,Ty,Th,p,hk,higher_heating_value,Qout); 
    ea=Math.abs(x3-p);
    }
    iter++;
} 
if(iter>=maxit) JOptionPane.showMessageDialog(null,"Uyarı maximum iterasyon sayısı aşıldı \n"+
   " çözüm geçerli olmıyabilir","MAKSİMUM ITERASYON SAYISI UYARISI",JOptionPane.WARNING_MESSAGE); 
return p;
}


public  double hki(String g_fuelname[],double g_fuelmol[],double Ty,double Th,double Te,boolean higher_heating_value,double Qout)
{
double a=1.0;
double b=50.0;
// brent metodu ile kök değerini bulur
double test;
double p=0;
double es,ea;
double f1,f2,f3,fp;
int maxit=500,iter=0;
double tol=1.0e-5;
es=0.001;
double x1=a;f1=fp(g_fuelname,g_fuelmol,Ty,Th,Te,x1,higher_heating_value,Qout);
double x2=b;f2=fp(g_fuelname,g_fuelmol,Ty,Th,Te,x2,higher_heating_value,Qout);
double x3=(x1+x2)/2.0;f3=fp(g_fuelname,g_fuelmol,Ty,Th,Te,x3,higher_heating_value,Qout);
if(f1==0) return x1;
else if(f2==0) return x2;
else if(f3==0) return x3; 
//if(f1*f2>0) System.out.println("verilen bölgede kök yok");
p=-(f2*f3*x1*(f2-f3)+f3*f1*x2*(f3-f1)+f1*f2*x3*(f1-f2))/((f1-f2)*(f2-f3)*(f3-f1));
fp=fp(g_fuelname,g_fuelmol,Ty,Th,Te,p,higher_heating_value,Qout);
//System.out.println("p="+p+"fp="+fp);  
ea=Math.abs(x3-p);
while((ea>es)&&(iter<maxit))
{	if(Math.abs(f3)<tol) return x3; 
    if((p<x1) && (p>x2)) 
    {p=(x1+x2)/2.0;    
     if(p>x3) {x1=x3;f1=f3;x3=p;f3=fp;}
     else if(p<x3) {x2=x3;f2=f3;x3=p;f3=fp;}
    }
    else
    {        
    if(p>x3) {x1=x3;f1=f3;x3=p;f3=fp;}
    else if(p<x3) {x2=x3;f2=f3;x3=p;f3=fp;}
    p=-(f2*f3*x1*(f2-f3)+f3*f1*x2*(f3-f1)+f1*f2*x3*(f1-f2))/((f1-f2)*(f2-f3)*(f3-f1));
    fp=fp(g_fuelname,g_fuelmol,Ty,Th,Te,p,higher_heating_value,Qout); 
    ea=Math.abs(x3-p);
    }
    iter++;
} 
if(iter>=maxit) JOptionPane.showMessageDialog(null,"Uyarı maximum iterasyon sayısı aşıldı \n"+
   " çözüm geçerli olmıyabilir","MAKSİMUM ITERASYON SAYISI UYARISI",JOptionPane.WARNING_MESSAGE); 
return p;
}


}
