public class brayton4
{
Gas g;
double PP1[],vv1[];
double PP2[],vv2[];
double TT1[],ss1[];
double TT2[],ss2[];
double Wturbine,Wcompressor1,Wcompressor2,Wcompressor,Wnet,Q1,Qregenerator,Qintercooler,Qnet,eta;
double T[],P[],h[],s[];
double m; // air mass flow rate;
combustion cc;
double a_f; //air to fuel ratio
double mair;
double mratio;
double mproduct; //proct mass flow rate kg/s  
double mfuel;
double ds5p; //entropy difference of two air and product gases
double dh5p; //enthalpy difference of two air and product gases
public brayton4(String fueli[],double Nfueli[],double mi,double T1,double P1,double Tfuel,double P2,double P4,double T3,double T6,double dT,double etacompressor1,double etacompressor2,double etaturbine1,double etaturbine2)
{  
   m=mi;
   g=new Gas("air");	
   g.base("mass");
   T=new double[11];
   P=new double[11];
   h=new double[11];
   s=new double[11];
   TT1=new double[25];ss1=new double[25];
   TT2=new double[2];ss2=new double[2];
// 1 compressor 1 inlet
P[1]=P1;T[1]=T1;
P[2]=P2; 
P[3]=P2;T[3]=T3;
P[4]=P4;
P[5]=P4;
P[6]=P4;T[6]=T6;
P[7]=P1;
P[8]=P1;
// 1 inlet
double a[]=g.property(T[1],P[1]);
h[1]=a[3];
s[1]=a[5];
TT1[0]=T[1];ss1[0]=s[1];
double pr1=a[17];
// 2i = 11 isentropic outlet of compressor 1
double pr2=P2/P1*pr1;
double T2i=g.Ti(pr2);
a=g.property(T2i,P2);
double h2i=a[3];
double s2i=s[1];//a[5];
double Cp2i=a[9];
T[9]=T2i;
s[9]=s2i;
TT1[1]=T2i;ss1[1]=s2i;
// 2 outlet of compressor 1
h[2]=(h2i-h[1])/etacompressor1+h[1];
T[2]=(h[2]-h2i)/Cp2i+T2i;
a=g.property(T[2],P[2]);
s[2]=a[5];
// 3 compressor 2 inlet
a=g.property(T[3],P[3]);
h[3]=a[3];
s[3]=a[5];
double pr3=a[17];
// 4i = 12 isentropic process compressor 2 outlet
double pr4=P[4]/P[3]*pr3;
double T4i=g.Ti(pr4);
a=g.property(T4i,P[4]);
double h4i=a[3];
double s4i=a[5];
double Cp4i=a[9];
T[10]=T4i;
s[10]=s4i;
//4 outlet of compressor 2
h[4]=(h4i-h[3])/etacompressor2+h[3];
T[4]=(h[4]-h4i)/Cp4i+T4i;
a=g.property(T[4],P[4]);
s[4]=a[5];
double Cp4=a[9];
//6 combustion chamber outlet
   a_f=1.0;
   cc=new combustion("combustion chamber",fueli,Nfueli,a_f);
   a_f=cc.adiabatic_a_f_ratio(T1,Tfuel,T6,1.0,10.0);
   mair=cc.nair/28.97;//mass flow rate of air from reaction
   mratio=m/mair;//air mass flow rate correction ratio compare to input air flow rate
   mproduct=(cc.productmix.N/cc.productmix.M)*mratio; //proct mass flow rate kg/s  
   mfuel=cc.fuelmix.N/cc.fuelmix.M*mratio;//fuel mass flow rate kg/s
   cc.fuelmix.base("mass");
   cc.productmix.base("mass");//productmix calculations will be mass bases kJ/kg..
   a=cc.productmix.property(T[6],P[6]);
   h[6]=a[3];
   s[6]=a[5];
   double Cp6=a[9];
   double pr6=a[17];
    a=g.property(T[6],P[6]);
   //fuel inlet enthalpy
   a=cc.fuelmix.property(Tfuel,P[6]);  
   double hfuel=a[3];
   //7i isentropic turbine outlet
   double pr7=P[7]/P[6]*pr6;
   double T7i=cc.productmix.Ti(pr7);
   a=cc.productmix.property(T7i,P[7]);
   double h7i=a[3];
   double s7i=a[5];
   double Cp7i=a[9];
   T[0]=T7i;
   s[0]=s7i;
   //7 turbine outlet
   h[7]=h[6]-etaturbine1*(h[6]-h7i);
   T[7]=(h[7]-h7i)/Cp7i+T7i;
   a=cc.productmix.property(T[7],P[7]);
   s[7]=a[5];
   double Cp7=a[9]/cc.productmix.M;
   //8 heat exchanger outlet
   T[8]=T[4]+dT;
   a=cc.productmix.property(T[8],P[8]);
   h[8]=a[3];
   s[8]=a[5];
   double Cp8=a[8]/cc.productmix.M;
   //5 heat exchanger exit
   Qregenerator=mproduct*(h[7]-h[8]);
   T[5]=T[8]+Qregenerator/mproduct/Cp8;
   a=g.property(T[5],P[5]);
   h[5]=a[3];
   s[5]=a[5];
   //5' heat exchanger exit calculated as product gases
   //normalisation of two gas data
   a=cc.productmix.property(T[5],P[5]);
   double h5p=a[3];
   double s5p=a[5];
   ds5p=s5p-s[5];
   //System.out.println("ds5p="+ds5p);
   s[6]-=ds5p;
   s[7]-=ds5p;
   s[8]-=ds5p;
   s[0]-=ds5p;
Wturbine=mproduct*(h[6]-h[7]);
Wcompressor1=m*(h[2]-h[1]);
Wcompressor2=m*(h[4]-h[3]);
Wcompressor=Wcompressor1+Wcompressor2;
Q1=mproduct*h[6]-m*h[5]-mfuel*hfuel;
double Qregenerator1=m*(h[5]-h[4]);
Qintercooler=m*(h[2]-h[3]);
Qnet=Q1;
Wnet=Wturbine-Wcompressor;
eta=Wnet/Qnet;
}
public double[][] line(double T1,double T2,double Pi,int n)
{
double TT1[]=new double[n];
double ss1[]=new double[n];
double aa[];
for(int i=0;i<n;i++)
{TT1[i]=T1+i*(T2-T1)/(n-1);
 aa=g.property(TT1[i],Pi);
 ss1[i]=aa[5];
}
double a[][]={TT1,ss1};
return a;
}
public double[][] line1(double T1,double T2,double Pi,int n)
{
double TT1[]=new double[n];
double ss1[]=new double[n];
double aa[];
for(int i=0;i<n;i++)
{TT1[i]=T1+i*(T2-T1)/(n-1);
 aa=cc.productmix.property(TT1[i],Pi);
 ss1[i]=aa[5]-ds5p;
}
double a[][]={TT1,ss1};
return a;
}
public void plot()
{double a[][]=new double[2][61];
 double c1[][]=new double[2][3];
 double c3[][]=new double[2][3];
 double c6[][]=new double[2][3];
 c1[0][0]=T[1];
 c1[1][0]=s[1];
 c1[0][1]=T[9];
 c1[1][1]=s[9];
 c1[0][2]=T[2];
 c1[1][2]=s[2];
 
 c3[0][0]=T[3];
 c3[1][0]=s[3];
 c3[0][1]=T[10];
 c3[1][1]=s[10];
 c3[0][2]=T[4];
 c3[1][2]=s[4];
 
 c6[0][0]=T[6];
 c6[1][0]=s[6];
 c6[0][1]=T[0];
 c6[1][1]=s[0];
 c6[0][2]=T[7];
 c6[1][2]=s[7];
 
 int i=0;
 a[0][0]=T[1];
 a[1][0]=s[1];
 double b[][]=line(T[2],T[3],P[2],15);
 for(i=0;i<15;i++)
 {a[0][i+1]=b[0][i];
  a[1][i+1]=b[1][i];
 }
 b=line(T[4],T[5],P[4],15);
 for(i=0;i<15;i++)
 {a[0][i+16]=b[0][i];
  a[1][i+16]=b[1][i];
 }
 b=line1(T[5],T[6],P[4],15);
 for(i=0;i<15;i++)
 {a[0][i+31]=b[0][i];
  a[1][i+31]=b[1][i];
 }
 b=line1(T[7],T[8],P[7],15);
 for(i=0;i<15;i++)
 {a[0][i+46]=b[0][i];
  a[1][i+46]=b[1][i];
 }

 Plot p=new Plot(a[1],a[0]);
 p.addData(c1[1],c1[0],1);
 p.addData(c3[1],c3[0],1);
 p.addData(c6[1],c6[0],1);
 p.addData(s,T,22);
 p.setPlabel("Brayton (Gas turbine) cycle with combustion process"); 
 p.setXlabel("s entropy kJ/(kgK)") ;
 p.setYlabel("T degree C");
 p.plot();
}

public String toString()
{
String ss=" Brayton cycle \n";
ss+="Combustion Reaction : \n";
ss+=cc.toString()+"\n";
ss+="mair ="+m+"kg/s mproduct="+mproduct+" kg/s mfuel\n";
ss+="Turbine power output W turbine = "+Wturbine+" kW \n";
ss+="Low pressure compressor power input W turbine = "+Wcompressor1+" kW \n";
ss+="High pressure compressor power input W turbine = "+Wcompressor2+" kW \n";
ss+="Compressor power input W turbine = "+Wcompressor+" kW \n";
ss+="Net power output  Wnet = "+Wnet+" kW\n";
ss+="Combustion chamber heat input Q combustion = "+Qnet+" kW\n";
ss+="Cycle efficiency eta  = "+eta+"\n";
ss+="Equivalent temperature carnot cycle efficiency eta  = "+(1.0-T[1]/T[6])+"\n";
ss+="Regenerator heat transfer   = "+Qregenerator+" kW\n";
ss+="Termodynamic states : \n";
for(int i=1;i<9;i++)
{ss+="P["+i+"] = "+P[i]+"bar T["+i+"] = "+(T[i]-273.15)+" degree C h["+i+"] = "+h[i]+" kJ/kg s["+i+"] = "+s[i]+" kJ/kgK\n";}
return ss;
}
}