public class Example6_5
{
	double k;
	double ro;
	double cp;
	double L;
	double alfa;
	public Example6_5(double ik,double iro,double icp,double iL)
	{
	L=iL;	
	cp=icp;	
	ro=iro;	
	k=ik;	
	alfa=k/(ro*cp);
	}
	
	double betam(double m)
	{
    return m*Math.PI/L;		
	}  

	double ft1(double t)
	{
	return 100.0*Math.sin(3.1415*t/180.0);	
	}
	
	double ft2(double t)
	{
	return 100.0*Math.sin(3.1415*t/180.0);	
	}
	
	double Fx(double x)
	{
	return 100*(3.0*x*x+5*x);	
	}

	double Gx(double x,double to)
	{
	return 3.0*x*x+4.0*to;	
	} 
	
	double func1int(double m,double x)
	{
	return Math.sin(betam(m)*x)*Fx(x);	
	}
	
	public  double newton_cotes1(double a,double b, int n,double m)
    {       //Newton-cotes integral 5 points
        //Boole rule
        // ff integrating function
        // a,b integral limits
        // n : number of sub division h=(b-a)/n;
	    double h=(b-a)/n;
        double h1=h/4.0;
        double sum=0;
        int p,i;
        double x[]=new double[5];
        double f[]=new double[5];
        for(p=0;p<n;p++)
        {  {for(i=0;i<5;i++) {x[i]=a+i*h1+p*h;f[i]=func1int(m,x[i]);}}
           sum+=h*(7.0*f[0]+32.0*f[1]+12.0*f[2]+32.0*f[3]+7.0*f[4])/90.0; 
        }
        return sum;
   }

    double func1(double t,double xx)
	{ 
		int n=20;
		double a=0.0;
		double b=L;
		double sum=0.0;
		for(double m=1;m<20.0;m=m+1.0)
		{
		sum+=Math.exp(-alfa*betam(m)*betam(m)*t)*Math.sin(betam(m)*xx)*newton_cotes1(a,b,n,m);	
		}
		return 2.0*sum/L;
	}
	
	///////////////////////////////////////////////////////////////////////////////////////////
	
	
	double intfunc2(double x,double to,double m)
	{
	return Math.sin(betam(m)*x)*Gx(x,to);	
	}
	
	public  double newton_cotes2ic(double m,double to)
    {       //Newton-cotes integral 5 points
        //Boole rule
        // ff integrating function
        // a,b integral limits
        // n : number of sub division h=(b-a)/n;
        double a=0.0;
        double b=L;
        int n=20;
	    double h=(b-a)/n;
        double h1=h/4.0;
        double sum=0;
        int p,i;
        double x[]=new double[5];
        double f[]=new double[5];
        for(p=0;p<n;p++)
        {  {for(i=0;i<5;i++) {x[i]=a+i*h1+p*h;f[i]=intfunc2(x[i],to,m);}}
           sum+=h*(7.0*f[0]+32.0*f[1]+12.0*f[2]+32.0*f[3]+7.0*f[4])/90.0; 
        }
        return sum;
   }
	
    double intfunc2dıs(double to,double m)
	{
	return Math.exp(alfa*betam(m)*betam(m)*to)*newton_cotes2ic(m,to);	
	}
	
	public  double newton_cotes2dıs(double m,double t)
    {       //Newton-cotes integral 5 points
        //Boole rule
        // ff integrating function
        // a,b integral limits
        // n : number of sub division h=(b-a)/n;
        double a=0.0;
        double b=t;
        int n=20;
	    double h=(b-a)/n;
        double h1=h/4.0;
        double sum=0;
        int p,i;
        double x[]=new double[5];
        double f[]=new double[5];
        for(p=0;p<n;p++)
        {  {for(i=0;i<5;i++) {x[i]=a+i*h1+p*h;f[i]=intfunc2dıs(x[i],m);}}
           sum+=h*(7.0*f[0]+32.0*f[1]+12.0*f[2]+32.0*f[3]+7.0*f[4])/90.0; 
        }
        return sum;
   }
	
   double func2(double t,double xx)
   {
   double sum=0.0;
   for(double m=1;m<20;m=m+1.0)
   {
   sum+=Math.exp(-alfa*betam(m)*betam(m)*t)*Math.sin(betam(m)*xx)*newton_cotes2dıs(m,t);	   
   }	   
   return ((alfa*2.0)/(k*L))*sum;	   
   }
   
   /////////////////////////////////////////////////////////////////////////////////////////
   
   
    double intfunc3(double to,double m)
	{
	return Math.exp(alfa*betam(m)*betam(m)*to)*ft1(to);	
	}
	
	public  double newton_cotes3(double m,double t)
    {       //Newton-cotes integral 5 points
        //Boole rule
        // ff integrating function
        // a,b integral limits
        // n : number of sub division h=(b-a)/n;
        double a=0.0;
        double b=t;
        int n=20;
	    double h=(b-a)/n;
        double h1=h/4.0;
        double sum=0;
        int p,i;
        double x[]=new double[5];
        double f[]=new double[5];
        for(p=0;p<n;p++)
        {  {for(i=0;i<5;i++) {x[i]=a+i*h1+p*h;f[i]=intfunc3(x[i],m);}}
           sum+=h*(7.0*f[0]+32.0*f[1]+12.0*f[2]+32.0*f[3]+7.0*f[4])/90.0; 
        }
        return sum;
   }
	
   double func3(double t,double xx)
   {
	 double sum=0.0;
	 for(double m=1;m<21.0;m=m+1.0)
	 {
	 sum+=Math.exp(-alfa*betam(m)*betam(m)*t)*betam(m)*Math.sin(betam(m)*xx)*newton_cotes3(m,t);	 
	 }  
	 return (alfa*2.0/L)*sum; 
   }
	
  ///////////////////////////////////////////////////////////////////////////////////////////////////
    double intfunc4(double to,double m)
	{
	return Math.exp(alfa*betam(m)*betam(m)*to)*ft2(to);	
	}
  
	public  double newton_cotes4(double m,double t)
    {       //Newton-cotes integral 5 points
        //Boole rule
        // ff integrating function
        // a,b integral limits
        // n : number of sub division h=(b-a)/n;
        double a=0.0;
        double b=t;
        int n=20;
	    double h=(b-a)/n;
        double h1=h/4.0;
        double sum=0;
        int p,i;
        double x[]=new double[5];
        double f[]=new double[5];
        for(p=0;p<n;p++)
        {  {for(i=0;i<5;i++) {x[i]=a+i*h1+p*h;f[i]=intfunc4(x[i],m);}}
           sum+=h*(7.0*f[0]+32.0*f[1]+12.0*f[2]+32.0*f[3]+7.0*f[4])/90.0; 
        }
        return sum;
   }
   
   double func4(double t,double xx)
   {
	 double sum=0.0;
	 for(double m=1;m<21.0;m=m+1.0)
	 {
	 sum+=Math.pow(-1.0,m)*Math.exp(-alfa*betam(m)*betam(m)*t)*betam(m)*Math.sin(betam(m)*xx)*newton_cotes3(m,t);	 
	 }  
	 return (alfa*2.0/L)*sum; 
   }
   
   double T(double t,double xx)
   {
	return func1(t,xx)+func2(t,xx)+func3(t,xx)-func4(t,xx);   
   }
	
	public static void main(String args[])
	{
	double k=10;
	double ro=100;
	double cp=2000;
	double L=1.0;	
	
		
		Example6_5 e65=new Example6_5(k,ro,cp,L);
		
		double Tdep[][]=new double[30][11];
	
		int i=0;
	  	for(double t=1.0;t<31.0;t=t+1.0)
	  	{
		  int j=0;
		  for(double xx=0.0;xx<=1.0;xx=xx+0.1)
		  {
		  Tdep[i][j]=e65.T(t,xx);
		  j++;
          }
          i++;
	    }
        System.out.println(Matrix.toString(Tdep));
	    timeplot tp=new timeplot(Tdep,0.0,L,1.0,1.0);
	    tp.start();
        
	} 
}