//======================================================
// Numerical Analysis package in java
// example to show differential equation solution
// and differentiation (derivative) functions
// Dr. Turhan Coban
// =====================================================

import java.io.*;

class fm1 extends f_xi
{
//multivariable function
 
 double func(double x[],int x_ref)
 {
   double p=1.0;     //atm
   double R=82.05;   //cm^3 atm/K/gmol)
   double T=1500.0;  //K
   double xfi=1/24.8;//original fuel
   double n=p/(R*T); //mol çarpanı
   double nfi=xfi*n;
   double nCO2i=0.0;
   double nCO2f=3*nfi;   
   double nO2=(3.5*x[1])+(1.5*nfi-0.5*x[2]);
   double nCO=3.0*(nfi-x[1])-x[2];
   //System.out.println("nfi="+nfi+"ncO2f="+nCO2f+"nco="+(nCO/nCO2f));
   double nH2O=4.0*(nfi-x[1]);
   double rf1=-1.0e12*Math.exp(-30000.0/(1.987*T))*Math.pow(x[1],0.1)*Math.pow(nO2,1.65);
   double rf2=Math.pow(10.0,14.6)*Math.exp(-40000.0/(1.987*T))*nCO*Math.pow(nO2,0.25)*Math.pow(nH2O,0.5)-
   5.0e8*Math.exp(-40000/(1.987*T))*x[2];
   //System.out.println("t="+x[0]+"x="+x[1]+"no2="+nO2+"x[2]="+x[2]+"rf="+rf1+"rf2="+rf2);
   if(x_ref==0)  return rf1;
   if(x_ref==1)  return rf2; 
   else return 0;
 }
}

public class diferansiyel3_2
{

public static double[][] RK6(f_xi fp,double x0,double xn,double f0[],int N)
{
//6th order Runge Kutta Method
//fp  : given set of derivative functions dfj/dxi(fj,x)
// xo : initial value of the independent variable
// xn : final value of the independent variable
// f0 : initial value of the dependent variable
// N  : number of dependent variable to be calculated
// fi : dependent variable
double h=(xn-x0)/N;
int M=f0.length;
double fi[][];
fi=new double[M][N+1];
double xi[]=new double[M+1];
double k[]=new double[6];
int i,j;
double x;
for(j=0;j<M;j++)
  {
  fi[j][0]=f0[j];
  xi[j+1]=f0[j];
  }
for(x=x0,i=0;i<N;x+=h,i++)
  {
  for(j=1;j<=M;j++)
    {
    xi[0]=x;
    xi[j]=fi[j-1][i];
    k[0]=h*fp.func(xi,j-1);
    xi[0]=x+h/2.0;
    xi[j]=fi[j-1][i]+k[0]/2;
    k[1]=h*fp.func(xi,j-1);
    xi[0]=x+h/2.0;
    xi[j]=fi[j-1][i]+k[0]/4.0+k[1]/4.0;
    k[2]=h*fp.func(xi,j-1);
    xi[0]=x+h;
    xi[j]=fi[j-1][i]-k[1]+2.0*k[2];
    k[3]=h*fp.func(xi,j-1);
    xi[0]=x+2.0/3.0*h;
    xi[j]=fi[j-1][i]+7.0/27.0*k[0]+10.0/27.0*k[1]+1.0/27.0*k[3];
    k[4]=h*fp.func(xi,j-1);
    xi[0]=x+1.0/5.0*h;
    xi[j]=fi[j-1][i]+28.0/625.0*k[0]-1.0/5.0*k[1]+546.0/625.0*k[2]+54.0/625.0*k[3]-378/625.0*k[4];
    k[5]=h*fp.func(xi,j-1);
    fi[j-1][i+1]=fi[j-1][i]+k[0]/24.0+5.0*k[3]/48.0+27.0*k[4]/56.0+125.0*k[5]/336.0;
    xi[j]=fi[j-1][i];
    }
  }
double a[][]=new double[M+1][N+1];
for(x=x0,i=0;i<=N;x+=h,i++)
{
  a[0][i]=x;
  for(j=1;j<=M;j++)
  {a[j][i]=fi[j-1][i];}
}

return a;
}

public static void main(String args[]) throws IOException
  {

  fm1 b3=new fm1();
   double y[];
   y=new double[2];
   double p=1.0;     //atm
   double R=82.05;   //cm^3 atm/K/gmol)
   double T=1500.0;  //K
   double xfi=1/24.8;
   double n=p/(R*T);
   double nfi=xfi*n;
   double nCO2i=0.0;
   double nCO2f=3.0*nfi;
   System.out.println("n="+n+"nfi="+nfi+"nco2i="+nCO2i+"nCO2f="+nCO2f);
   y[0]=nfi;
   y[1]=nCO2i;
   double a[][]=RK6(b3,0.0,3.02e-3,y,2000);
   double x1[]=a[0];
   double y1[]=a[1];
   double y2[]=a[2];
   double y3[]=new double[a[0].length];
   double y4[]=new double[a[0].length];
   double y5[]=new double[a[0].length];   
   for(int i=0;i<a[0].length;i++) 
   {x1[i]=a[0][i];
    y1[i]=a[1][i];        // fuel
    y2[i]=a[2][i];        // CO2 
    y3[i]=((3.5*a[1][i])+(1.5*nfi-0.5*a[2][i])); //O2
    y4[i]=(3.0*(nfi-a[1][i])-a[2][i]); //CO
    y5[i]=4.0*(nfi-a[1][i]); //H2O
    System.out.println(x1[i]+" "+y1[i]+" "+y2[i]+" "+y3[i]+" "+y4[i]+" "+y5[i]);
   }
   //normalize the results
   for(int i=0;i<a[0].length;i++) 
   {
    y1[i]=y1[i]/nfi; //fuel
    y2[i]=y2[i]/nCO2f;        //CO2 
    y3[i]=y3[i]/nfi;//O2
    y4[i]=y4[i]/nCO2f;//CO
    y5[i]=y5[i]/nfi;//H2O
    System.out.println(x1[i]+" "+y1[i]+" "+y2[i]+" "+y3[i]+" "+y4[i]+" "+y5[i]);
   }
   
   Plot pp=new Plot(x1,y1);
   pp.addData(x1,y2);
   //pp.addData(x1,y3); 
   pp.addData(x1,y4);
   //pp.addData(x1,y5);    
   pp.plot();
 }
}
