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

public class NC_vertical_channel
 { double a,b,tickness,H,dH,L,dL;
   double t[][];
   double tininf[],toutinf;
   int n,m;
   double Q[];      // heat transfer in each layer
   double Qair[];   // heat transfer in each layer inside the channel
   double Qbase[];
   double Qt;       // Total heat transfer
   double alfa,h,ro,nu,Cp,alpha,Pr,Nu,Ra,Gr,Re; // Nusselt, Prandtl,Grashoff, Rayleigh and Reynolds numbers
   double U,deltat; // velocity and thermal boundary layer thickness
   double Aref;     //referance channel area a*b/2;
   double Aboundary;//thermal boundary layer flow area
   // a : channel length (x)
   // b : channel with   (y)
   // tickness : fin thickness
   // H :channel height  (z)
   // dH: finite difference dz
   // L : a+b/2; channel will be considered as a single fin of length L with isolated end point
   // dL finite difference  
  public NC_vertical_channel(double ai,double bi,double ti, double Hi,double dHi,double dLi)
  {a=ai;b=bi;tickness=ti;H=Hi;L=a+b/2.0;dH=dHi;dL=dLi;n=(int)(H/dH+1.0);m=(int)(L/dL+1.0);
  System.out.println("n="+n+"m="+m);
  t=new double[n+1][m];
  tininf=new double[n+1];
  Q=new double[n+1];
  Qair=new double[n+1];
  Qbase=new double[n+1];
  }
  public static double kmetal(double t)
  {double T=t+273.15;
   double a=71.62587939698594-0.03447361809045638*T-7.927135678388277E-6*T*T;
  return a;
  }
  public double k(double t)
  {double T=t+273.15;
  Gas g=new Gas("air");
   double ka=g.k(T);
  return ka;
  }	 	   	 
  public double Rayleigh(double t,double tinf,double x) 
  { // Rayleigh number (natural convection)
    // t degree c fin surface temperature
    // tinf degree C outside air temperature
    // L height of plate m
    // h convective heat transfer coefficient W/(m^2K)
    double P=1.01325; //bar
    Gas g=new Gas("air");
   double tf=(t+tinf)/2.0;
   double Tf=tf+273.15;
   double T=t+273.15; 
   double beta=1/Tf;
   double gyc=9.806;
   ro=g.M/g.v(Tf,P);
   nu=g.vis(Tf)/ro;
   alfa=g.alfa(Tf,P);
   Pr=g.Prandtl(Tf);
   Cp=g.Cp(Tf)/g.M;
   Gr=gyc*beta*(t-tinf)*x*x*x/(nu*nu);
   Re=Math.sqrt(Gr);
   U=Re*nu/(ro*x);
   Ra=Gr*Pr;
   deltat=5.5*x/Math.pow(Gr,0.25);
   //System.out.println("Gr="+Gr+"Ra="+Ra+"t="+t+"tinf="+tinf+"tf="+tf+"vis="+g.vis(Tf)+"ro="+ro+"nu="+nu);
   //System.out.println("gy="+gyc+"nu="+nu+"beta="+beta+"alfa="+alfa);
   return Ra;
  }
  
  public double h(double t,double tinf,double x ) 
  {  // Natural convection heat transfer coefficient 
     // vertical plate
  Gas g=new Gas("air");
  double tf=(t+tinf)/2.0;
  double Tf=tf+273.15;
  double Pr=g.Prandtl(Tf);
  double Ra=Rayleigh(t,tinf,x);
  double k=g.k(Tf);
  double Nu;
  double x1,x2,x3;
  x2=Math.pow((0.492/Pr),(9.0/16.0));
  if(Ra<1e9)
  {
  x1=Math.pow(Ra,0.25);
  Nu=0.68+0.67*x1/Math.pow((1+x2),(4/9));
  } 
  else  
  {
  x1=Math.pow(Ra,(1.0/6.0));
  x2=Math.pow((0.492/Pr),(9/16));
  x3=0.825+0.387*x1/Math.pow((1.0+x2),(8.0/27.0));
  Nu=x3*x3;
  }  
  double hh=Nu*k/x;
  return hh;
  }
  
  public void tempprofile(double tb,double t2inf,int i)
  {t[i][0]=tb;   
   double t1f,t2f,h1=0,h2=0,h3=0;
   double dA=dH*dL;
   double dA1=b/2.0*dH;
   double dA2=0;
   double dAc=tickness*dH;
   double dQ=0;
   double dQair=0;
   double U0=0;
   double ro0=0;
   double Cp0=0;
   double dt=0;
   int j;
   Q[i]=0;
   double x=dH*(i+1);
   for(j=0;j<m-1;j++)
   { t1f=(t[i][0]+tininf[i])/2;
     t2f=(t[i][0]+t2inf)/2;
     h1=h(t[i][0], tininf[i],x ); 
 
     dA2+=dH*dL*deltat;
     U0+=U*dH*dL*deltat;
     ro0+=ro*dH*dL*deltat;
     Cp0+=Cp*dH*dL*deltat;
     h2=h(t[i][0], t2inf,x ); 
     
     if(t[i][0]>tininf[i]) {dQ=(h1*dA*(t[i][0]-tininf[i])+h2*dA*(t[i][0]-t2inf));dQair+=h1*dA*(t[i][0]-tininf[i]);}
     else                  {dQ=0;dQair+=0.0;}
     dt=dQ*dL/kmetal(t[i][0])/dAc;
     t[i][j+1]=t[i][j]-dt;
     Q[i]+=dQ;
     Qair[i]+=dQair;
     //System.out.println("i="+i+"dt="+dt+"Qi="+Q[i]+"Qair="+Qair+"dQ="+dQ+"tininf="+tininf[i]+"Pr="+Pr+"Gr="+Gr+"x="+x+"h1="+h1+"h2="+h2+"Ra="+Ra);
   } 
   h3=h(tb,tininf[i],x);
   dA2+=b/2.0*deltat;
   U0+=U*b/2.0*deltat;
   ro0+=ro*b/2.0*deltat;
   Cp0+=Cp*b/2.0*deltat;   
   U0/=dA2;  // average velocity
   ro0/=dA2; // average dendity
   Cp0/=dA2; // average specific heat
   Qbase[i]=h3*dA1*(tb-tininf[i]);
   Qair[i]+=Qbase[i];
   //if(dA2<dA1) {tininf[i+1]=tininf[i];System.out.println("dA2="+dA2+"dA1="+dA1);}
   //else   
   double dt1=Qair[i]/(ro0*U0*dA2*Cp0*1.0e3);
   double div=(ro0*U0*dA2*Cp0*1.0e3);
   if(dA2<dA1) tininf[i+1] = tininf[i];
   else{
   double t1=((dt1+tininf[i])*dA2+tininf[i]*(dA1-dA2))/dA1;//System.out.println("second region fully developed flow dA2="+dA2+"dA1="+dA1);
   tininf[i+1] = tininf[i]+t1;
   System.out.println("i="+i+"Qair="+Qair[i]+"div="+div+"U0="+U0+"dA2="+dA2+"dt="+dt1+"t1="+t1+"tininf[i]="+tininf[i]+"tininf[i+1]="+tininf[i+1]+"h1="+h1+"h2="+h2+"h3="+h3+"U0="+U0+"ro0="+ro0+"Cp0="+Cp0+"Qair="+Qair[i]); 
       }
} 
   
  public void  Q(double tb1,double tb2,double tinf)
  {   double tb;
      Qt=0;
      tininf[0]=tinf;
	  for(int i=0;i<n;i++)
	  {tb=tb1-(tb1-tb2)/(n-1)*i;
	   tempprofile(tb,tinf,i);
	   Qt+=Q[i]+Qbase[i];
	  }
	  System.out.println("Qt="+Qt);
  } 

   
  public static void main(String args[]) throws IOException
  {
  NC_vertical_channel nc=new NC_vertical_channel(2e-2,1.0e-2,0.4e-3,1.0,0.05,2.5e-3);
  double t1=90;
  double t2=70;
  double tinf=20.0;
  nc.Q(t1,t2,tinf);
  Text.print(nc.t);
  Text.print(nc.tininf);

   double a[][]=Text.readDouble("cc.txt");
      f2a f=new f2a(a,6);
	    surfacePlot p=new surfacePlot(f,0.0,1.0,0.0,0.05,70.0,90.0,"X","Y","Z",50);
        p.plot();
   
  }
 }