import java.io.*;

public class flatplate
{
	/////////////////////////
	//                     //
	//  OĞUZ EMRAH TURGUT  //
	//                     // 
	/////////////////////////
 
	 public static double [] thomas(double f[],double e[],double g[],double r[])
  {
  // 3 band matrix system of equation solving algorithm
  int n=f.length;
  double x[]=new double[n];
  for(int k=1;k<n;k++)
     {e[k]=e[k]/f[k-1];
      f[k]=f[k]-e[k]*g[k-1];
     }
  for(int k=1;k<n;k++)
     {r[k]=r[k]-e[k]*r[k-1];
     }
  x[n-1]=r[n-1]/f[n-1];   
  for(int k=(n-2);k>=0;k--)
     {x[k]=(r[k]-g[k]*x[k+1])/f[k];}
  return x;      
  }
  
  
  public static double [][] cubic_spline(double xi[],double yi[],double c0,double cn)
  {
  int n=xi.length;
  double h[]=new double[n];
  double w[]=new double[n];
  double f[]=new double[n];
  double e[]=new double[n];
  double g[]=new double[n];
  double d[]=new double[n];
  double x[]=new double[n];
  double S[][]=new double[4][n];
  int k;	  
  for(k=0;k<(n-1);k++)
     {h[k]=xi[k+1]-xi[k];
      w[k]=(yi[k+1]-yi[k])/h[k];      
	 }
  d[0]=c0;
  d[n-1]=cn;
  for(k=1;k<(n-1);k++)
     {d[k]=6.0*(w[k]-w[k-1]);}
  f[0]=1.0;
  f[n-1]=1.0;
  g[0]=0.0;
  g[n-1]=0.0;
  e[0]=0.0;
  e[n-1]=0.0;
  for(k=1;k<(n-1);k++)
  {f[k]=2.0*(h[k]+h[k-1]);e[k]=h[k-1];g[k]=h[k];}
  S[2]=thomas(f,e,g,d);
  S[3]=xi;
  for(k=0;k<(n-1);k++)
    {S[0][k]=(6.*yi[k+1]-h[k]*h[k]*S[2][k+1])/(6.0*h[k]);
     S[1][k]=(6.*yi[k]-h[k]*h[k]*S[2][k])/(6.0*h[k]);
    } 
  return S;    	 	     	 
  }
  
  public static double funcSpline(double S[][],double x)
  {
  int n=S[0].length;
  double xx1=0;
  double xx2=0;
  double y=0;
  double hk=0;
  for(int k=0;k<(n-1);k++)
     {if(S[3][k]<=x && x<=S[3][k+1]) 
      {hk=(S[3][k+1]-S[3][k]);
       xx1=(x-S[3][k]);
       xx2=(S[3][k+1]-x);
       y=S[0][k]*xx1+S[1][k]*xx2+(xx1*xx1*xx1*S[2][k+1]+xx2*xx2*xx2*S[2][k])/(6.0*hk);
       break;
       }
     }
       if(y==0 && S[3][n-2]<=x )
       {    
	   int k=n-2;
	   hk=(S[3][k+1]-S[3][k]);
       xx1=(x-S[3][k]);
       xx2=(S[3][k+1]-x);
       y=S[0][k]*xx1+S[1][k]*xx2+(xx1*xx1*xx1*S[2][k+1]+xx2*xx2*xx2*S[2][k])/(6.0*hk);
       }
  return y;
  }     
	
	public static double[] properties(double t)//in C
	{
	// properties of air	
	double[] T={-150,-100,-50,0,20,40,60,80,100,120,140,160,180,200,250,300,350,400};//C
	double[] ro={2.793,1.980,1.534,1.293,1.205,1.127,1.067,1.000,0.946,0.898,0.854,0.815,0.779,0.746,0.675,0.616,0.566,0.524};//kg/m3
	double[] cp={1026,1009,1005,1005,1005,1005,1009,1009,1009,1013,1013,1017,1022,1026,1034,1047,1055,1068}; //J/kgK	
	double[] k={0.0116,0.0160,0.0204,0.0243,0.0257,0.0271,0.0285,0.0299,0.0314,0.0328,0.0343,0.0358,0.0372,0.0368,0.0421,0.0454,0.0485,0.0515};//W/mK
	double[] kvis={3.08E-06,5.95E-06,9.55E-06,13.30E-06,15.11E-06,16.97E-06,18.90E-06,20.94E-06,23.06E-06,25.23E-06,27.55E-06,29.85E-06,32.29E-06,34.63E-06,41.17E-06,47.85E-06,55.05E-06,62.53E-06};//m2/s
	double[] Pr={0.76,0.74,0.725,0.715,0.713,0.711,0.709,0.708,0.703,0.700,0.695,0.69,0.69,0.69,0.685,0.680,0.680,0.680,0.680};
	double[] expansion={8.21E-03,5.82E-03,4.51E-03,3.67E-03,3.43E-03,3.20E-03,3.00E-03,2.83E-03,2.68E-03,2.55E-03,2.43E-03,2.32E-03,2.21E-03,2.11E-03,1.91E-03,1.75E-03,1.61E-03,1.49E-03};
	
	double[][] roint=cubic_spline(T,ro,0,0);
	double newro=funcSpline(roint,t);
	double[][] cpint=cubic_spline(T,cp,0,0);
	double newcp=funcSpline(cpint,t);
	double[][] kint=cubic_spline(T,k,0,0);
	double newk=funcSpline(kint,t);	
	double[][] kvisint=cubic_spline(T,kvis,0,0);
	double newkvis=funcSpline(kvisint,t);	
	double[][] Print=cubic_spline(T,Pr,0,0);
	double newPr=funcSpline(Print,t);	
	double[][] expansionint=cubic_spline(T,expansion,0,0);
	double newexpansion=funcSpline(expansionint,t);	
	
	double[] proper=new double[6]; 	
	proper[0]=newro;//density
	proper[1]=newcp;//specific heat
	proper[2]=newk;// heat cond coeff
	proper[3]=newkvis;// kinematic vis
	proper[4]=newPr;// Prandtl
	proper[5]=newexpansion;//Expansion coeff;
	return proper;
	}
	
	
	public static double Ut()
	{
	double N=2;// cam örtü sayısı	
    double beta=45;// kolektör eğimi 	
    double epsg=0.88;//emmitance of glass
    double epsp=0.1;//emmitance of plate	
    double Tamb=10+273;;//ambient temp	(K)
    double Tpm=100+273;//mean plate temp (K) 
    double hw=10; //wind heat transfer coefficient(W/m2C)	
	double sigma=5.670373E-08;// stefan--boltzmann coeff
    	
	double C=0.0;
	
	if((0<beta)&&(70>beta))
	{C=520*(1-0.000051*beta*beta);}
	else if((70<beta)&&(90>beta))
	{beta=70;
	C=520*(1-0.000051*beta*beta);
	}	
		
	double f=(1+0.089*hw-0.1166*hw*epsp)*(1+0.07866*N);
	double e=0.430*(1-(100/Tpm));
	double AA=1/((N/((C/Tpm)*Math.pow(((Tpm-Tamb)/(N+f)),e)))+(1/hw));
	double BB=sigma*(Tpm+Tamb)*(Tpm*Tpm+Tamb*Tamb);
	double CC=(1/(epsp+0.00591*N*hw))+(((2*N)+f-1+(0.133*epsp))/(epsg))-N;	
	double U=AA+(BB/CC);	
	return U;	
	}
	
	public static double toRadian(double beta)
	{return beta*Math.PI/180;}
	
	public static double[] iterations(double tc1,double tc2)
	{
	 double space=0.025;//m	
	 double	epsp=0.1;// plate emmitance
	 double epsg=0.88;//glass emmitance
	 double Tamb=10+273;//ambient temp K
	 double Tsky=10+273;//sky temp K
	 double hw=10;//W/m2C	
	 double Tmp=100+273;//mean plate temperature K
	 double sigma=5.670373E-08;// stefan--boltzmann coeff	
	 double g=9.801;
	 double Tc1=tc1+273;
	 double Tc2=tc2+273;
	 double beta=45;//tilt angle
	 double length=2;
	 double width=1;
	 double Ac=2*(length+width);//section perimeter
	 double insulationthickness=0.5E-3;//m
	 double A=length*width;
	 
	 double kinsulation=0.045;//W/m2C
	 double bckinsthckness=0.050;//m
	 
	 double collectorthck=2.3E-03;
	 
	 double tol=0.001;
	 double iter=1;
	 double Ut=0.0;
	 for(int i=0;i<3000000;i++)
	 {
	 double Tort1=((Tc1-273)+(Tmp-273))/2;
	 double Tort2=((Tc1-273)+(Tc2-273))/2;
	 
	 
	 double[] prop1=properties(Tort1);
	 double k1=prop1[2]; 
	 double ro1=prop1[0];
	 double kvis1=prop1[3];
	 double Pr1=prop1[4];
	 double expan1=prop1[5];
	 
	 double[] prop2=properties(Tort2);
	 double k2=prop2[2]; 
	 double ro2=prop2[0];
	 double kvis2=prop2[3];
	 double Pr2=prop2[4];
	 double expan2=prop2[5];
	 
	 double hrp_c1=(sigma*(Tmp*Tmp+Tc1*Tc1)*(Tmp+Tc1))/((1/epsp)+(1/epsg)-1);	
	 double hrc1_c2=(sigma*(Tc1+Tc2)*(Tc1*Tc1+Tc2*Tc2))/((1/epsg)+(1/epsg)-1);
	 double hrc2_ca=epsg*sigma*(Tc2*Tc2+Tsky*Tsky)*(Tc2+Tsky)*(Tc2-Tsky)/(Tc2-Tamb);
	 
	 double Ra1=(g*expan1*(Tmp-Tc1)*space*space*space)/(kvis1*kvis1/Pr1);	
	 double beta_rd=toRadian(beta);
	 // Nusselt number calculation(natural convection)
	 double value=1.0;
	 double value2=1.0;
	 double decider1=(1-(1708/(Ra1*Math.cos(beta_rd))));
	 double decider2=(Math.pow((Ra1*Math.cos(beta_rd)/5830),1/3)-1);
	 if (decider1>0)
	 {value=decider1;}
	 else
	 {value=0.0;}	
	 if(decider2>0)
	 {value2=decider2;}	
	 else
	 {value2=0.0;}	
	 double Nu=1+((1.44*(1-((1708*Math.pow(Math.sin(1.8*beta_rd),1.6))/(Ra1*Math.cos(beta_rd)))))*value)+value2;
	 double hncp_c1=Nu*k1/space;
	 
	 double Ra2=(g*expan2*(Tc1-Tc2)*space*space*space)/(kvis2*kvis2/Pr2);	
	 double value3=1.0;
	 double value4=1.0;
	 double decider3=(1-(1708/(Ra2*Math.cos(beta_rd))));
	 double decider4=(Math.pow((Ra2*Math.cos(beta_rd)/5830),1/3)-1);
	 if (decider3>0)
	 {value3=decider3;}
	 else
	 {value=0.0;}	
	 if(decider4>0)
	 {value4=decider4;}	
	 else
	 {value4=0.0;}	
	 double Nu2=1+((1.44*(1-((1708*Math.pow(Math.sin(1.8*beta_rd),1.6))/(Ra2*Math.cos(beta_rd)))))*value3)+value4;
	 double hnc1_c2=Nu2*k2/space;
	 
	 double R1=1/(hw+hrc2_ca);
	 double R2=1/(hnc1_c2+hrc1_c2);
	 double R3=1/(hrp_c1+hncp_c1);
	  Ut=1/(R1+R2+R3);
	 double Tdummy1=Tmp-((Ut*(Tmp-Tamb))/(hncp_c1+hrp_c1));
	 double Tdummy2=Tdummy1-(Ut*(Tmp-Tamb)/(hnc1_c2+hrc1_c2));
	 iter+=1;
	 if((Math.abs(Tdummy1-Tc1)<tol)&&(Math.abs(Tdummy2-Tc2)<tol))
	 {break;}
	 Tc1=Tdummy1;
	 Tc2=Tdummy2;
	 }	
	 
	 double Ub=kinsulation/bckinsthckness;
	 double Ue=(kinsulation/insulationthickness)*Ac*collectorthck/A;
	 double Utop=Ut+Ub+Ue;
     
	 double[] returnval=new double[7];
     returnval[0]=Ut;
     returnval[1]=Ub;
     returnval[2]=Ue;
     returnval[3]=Utop;
     returnval[4]=Tc1;
     returnval[5]=Tc2;
     returnval[6]=iter;
     
     return returnval;	
  }
	
	public static void main(String args[])
	{
	double firstguessTc1=60;
	double firstguessTc2=50;
	double[] values=iterations(firstguessTc1,firstguessTc2);
	double Ucalc1=values[0];
	double Ubottom=values[1];
	double Uedge=values[2];
	double Usum=values[3];
	double Tcam1=values[4]-273;
	double Tcam2=values[5]-273;
	double iterasyon=values[6];
	
	
	double U=Ut();
	System.out.println("****************************************************" );
	System.out.println("* 1.cam örtünün ilk tahmin sıcaklığı = "+firstguessTc1+" C");
	System.out.println("* 2.cam örtünün ilk tahmin sıcaklığı = "+firstguessTc2+" C");
	System.out.println("* İterasyon sonucu bulunan Ut değeri = "+Ucalc1+" W/m2K" );
	System.out.println("* Ubottom değeri = "+Ubottom+" W/m2K" );
	System.out.println("* Uedge değeri = "+Uedge+" W/m2K" );
	System.out.println("* Toplam U değeri = "+Usum+" W/m2K" );
	System.out.println("* 1.cam örtünün sıcaklığı = "+Tcam1+" C");
	System.out.println("* 2.cam örtünün sıcaklığı = "+Tcam2+" C " );	
	System.out.println("* İterasyon sayısı = "+iterasyon);		
	System.out.println("* Hesaplanan Ut = "+U);
	System.out.println("****************************************************" );
	}
}