public class solar_eng
{
	
	
	
    
	public static double hour_angle(double saatfarkı,boolean b1) 
	{   // false saat 12 den önce 
	    // true  saat 12 den sonra  
		//saat farkı in hour ex 1.5 ,2.0....
		double f=0;
		if(b1==true){f=saatfarkı*15;}
		else if(b1==false){f=-saatfarkı*15;}
		return f;
		
	}
	
	
	public static  double toRadian(double C)
	{ 
		  return (Math.PI/180)*C;
	
	}
	
	
	public static double gunsayısı(String s1,double aygunu)
	{
		double n=0;
		double i=0;
		if(s1=="Ocak"){n=aygunu;}
		else if(s1=="Şubat"){n=aygunu+31;}
		else if(s1=="Mart"){n=aygunu+59;}
		else if(s1=="Nisan"){n=aygunu+90;}
		else if(s1=="Mayıs"){n=aygunu+120;}
		else if(s1=="Haziran"){n=aygunu+151;}
		else if(s1=="Temmuz"){n=aygunu+181;}
		else if(s1=="Ağustos"){n=aygunu+212;}
		else if(s1=="Eylül"){n=aygunu+243;}
		else if(s1=="Ekim"){n=aygunu+273;}
		else if(s1=="Kasım"){n=aygunu+304;}
		else if(s1=="Aralık"){n=aygunu+334;}
		return n;
	} 

	   public static double decl_angle(String s1,double aygunu)
    { 
	   double n1=gunsayısı(s1,aygunu);
	   double dummy=360*(284+n1)/365;
	   double rad=toRadian(dummy);
	   return 23.45*Math.sin(rad);    
	}

       public static double horiz_zenith_angle(double enlem,String s1,double aygunu,boolean b1,double saatfarkı)
    {
	    double decl_degre=decl_angle(s1,aygunu);
	    double decl_rad=toRadian(decl_degre);
	    double rad_enlem=toRadian(enlem);
	    double hour_angle_degree=hour_angle(saatfarkı,b1);
	    double hour_angle_rad=toRadian(hour_angle_degree);
	    double dum1=(Math.cos(rad_enlem)*Math.cos(decl_rad)*Math.cos(hour_angle_rad))+(Math.sin(decl_rad)*Math.cos(rad_enlem));
	    double zenith_rad=Math.acos(dum1);
	    return zenith_rad;
	    //return (180/Math.PI)*zenith_rad;
	}
	
	  public static double beam_radiation_transmittance(String s2,double rakım,double enlem,String s1,double aygunu,boolean b1,double saatfarkı) 
	{   // clear sky beam radiation transmittance
		//rakım in km 
		//2.8.1a-b-c-d
		//
		double ro=0;
		double r1=0;
		double rk=0;
		
		if(s2=="Tropikal"){ro=0.95;r1=0.98;rk=1.02;}
		else if(s2=="Orta_enlem_yaz"){ro=0.97;r1=0.99;rk=1.02;}
		else if(s2=="Supartik_yaz"){ro=0.99;r1=0.99;rk=1.01;}
		else if(s2=="Orta_enlem_kış"){ro=1.03;r1=1.01;rk=1.00;}
		
		double a00=0.4237-0.00821*(6-rakım)*(6-rakım);
		double a10=0.5055+0.00595*(6.5-rakım)*(6.5-rakım);
		double a20=0.2711+0.01858*(2.5-rakım)*(2.5-rakım);
		
		double a0=a00*ro;
		double a1=a10*r1;
		double k=a20*rk;
		
		double cosazi=Math.cos(horiz_zenith_angle( enlem, s1, aygunu, b1, saatfarkı));
		return a0+a1*Math.exp(-k/cosazi);
	}
	
	    public static double diffuse_radiation_transmittance(String s2,double rakım,double enlem,String s1,double aygunu,boolean b1,double saatfarkı) 	
	   {   //clear sky diffuse radiation
		   //(2.8.5)
	       double td=0.271-0.294*beam_radiation_transmittance(s2,rakım,enlem,s1,aygunu,b1,saatfarkı);    
		   return td;
	   }
		
	    public static double beam_radiation_horizontal(String s1,double aygunu,double enlem,boolean b1,double saatfarkı) 
	{
		
	    double Gsc=	1367;//W/m2
	    double n=gunsayısı(s1,aygunu);
	    double dum2=360*n/365;
	    double raddum=toRadian(dum2); 
	    double zenithrd=horiz_zenith_angle( enlem, s1, aygunu, b1, saatfarkı);
	    double Go=Gsc*(1+0.033*Math.cos(raddum))*Math.cos(zenithrd);
		return Go;//W/m2
	}
		
	    public static double G_on(String s1,double aygunu)
    {  
	    //1.4.1
	    // extrateressial radiation  on the nth day of the year
	    double Gsc=1367;//W/m2	    
	    double n=gunsayısı(s1,aygunu);	    
	    double Gon =Gsc*(1+(0.033*Math.cos(n*360/365))) ;   
		return Gon;    
    }
		
	public static double daily_solar_radiation(double enlem,String s1,double aygunu,double incline_angle,double surf_azimuth)
	{
		//daily solar radiation
		// H0
		//1.10.3
		double Gsc=1367;//W/m2
		double enlem_rd=toRadian(enlem);
		double decl_angle_rad=toRadian(decl_angle(s1,aygunu));
		double n=gunsayısı(s1,aygunu);
		double[] sunsetrize=sun_set_rise(s1,aygunu,enlem,incline_angle,surf_azimuth);
		double ws=sunsetrize[1];//güneşin batış açısı
		double ws_rd=toRadian(ws);
		double H0=(24*3600*Gsc/Math.PI)*(1+0.033*Math.cos(360*n/365))*(Math.cos(enlem_rd)*Math.cos(decl_angle_rad)*Math.sin(ws_rd)+(Math.PI*ws/180)*Math.sin(enlem_rd)*Math.sin(decl_angle_rad));
		return H0;//in J/m2
	}    
	
	    public static double hour_period_solar_radiation(double enlem,String s1,double aygunu,boolean b1,boolean b2,double saatfarkı1,double saatfarkı2)//
	{   //eq 1.10.4 
	    //I0
	    double Gsc=1367;//W/m2
	    double enlem_rd=toRadian(enlem);	    
		double decl_angle_rad=toRadian(decl_angle(s1,aygunu));
		double n=gunsayısı(s1,aygunu);
		double w1=hour_angle(saatfarkı1,b1);
		double w2=hour_angle(saatfarkı2,b2); 
		double w1_rd=toRadian(w1);
		double w2_rd=toRadian(w2); // w2 is larger than w1  
		double I0=(12*3600*Gsc/Math.PI)*(1+0.033*Math.cos(360*n/365))*((Math.cos(enlem_rd)*Math.cos(decl_angle_rad)*(Math.sin(w2_rd)-Math.sin(w1_rd)))+(Math.PI*(w2-w1)/180)*(Math.sin(enlem_rd))*Math.sin(decl_angle_rad)); 
        return I0;
		
	}
	
	public static double daylighthours(double enlem,String s1,double aygunu)
	{
	// number of daylight hours	
	double da=decl_angle( s1, aygunu);
	double enlem_rd=toRadian(enlem);
	double da_rd=toRadian(da);	
	return 2*(Math.acos(-Math.tan(enlem_rd)*Math.tan(da_rd)))/15;			
	}
	
	
	
	
	
	
	
	     // saatlik ışınım direct ve yayılı
	    // Erbs ve Orgill_Hollands  
	
	
		public static double Orgill_and_Hollands(double kt)
	{
	   double Id_I=0.0;
		 
	     if(kt<0.35){Id_I=1.0-0.249*kt;}
		 else if((kt>=0.35)&&(kt<=0.75)){Id_I=1.557-1.84*kt;}
		 else if((kt>0.75)){Id_I=0.177;}		
	   	
		 return Id_I;	
	}
	
	    
		
	    public static double Erbs(double kt)
	{
	    double IdI=0.0	;
	    if(kt<=0.22){IdI=1.0-0.09*kt;}
	    else if((kt>=0.22)&&(kt<=0.80)){IdI=0.9511-0.1604*kt+4.388*kt*kt-16.638*kt*kt*kt+12.336*kt*kt*kt*kt;}    
	  	else if(kt>0.80){IdI=0.165;}
	  	return IdI;    
		    
	}
	   // günlük ortalama ışınım
	
	    public static double Erbs(double Kt,double ws)
	{
	    double HdH=0.0;
		
		if(ws<=81.4)
		{if(Kt<0.715){HdH=1.0-0.2727*Kt+2.4495*Kt*Kt-11.9514*Kt*Kt*Kt+9.3879*Kt*Kt*Kt*Kt;}
		 else if(Kt>=0.715){HdH=0.143;}
		}     
	    else if(ws>81.4)
	    {if(Kt<=0.722){HdH=1.0+0.2832*Kt-2.5557*Kt*Kt+0.8448*Kt*Kt*Kt;}
	    else if(Kt>0.722){HdH=0.175;}
	    }
	    
	    return HdH;
			    
     }
     
       public static double Collares(double Kt)
     {
	     double Hdh=0.0;
	     
	     if(Kt<=0.17){Hdh=0.99;}
	     else if ((Kt>0.17)&&(Kt<0.75)){Hdh=1.188-2.272*Kt+9.473*Kt*Kt-21.865*Kt*Kt*Kt+14.648*Kt*Kt*Kt*Kt;}
	     else if ((Kt>=0.75)&&(Kt<0.80)){Hdh=-0.54*Kt+0.632;}
	     else if (Kt>=0.80){Hdh=0.2;}
	     
	     return Hdh;    
	 }
     
	 
	 // Güneşin doğuş ve batış dereceleri
	 
	 public static double[] sun_set_rise(String s1,double aygunu,double enlem,double incline_angle,double surf_azimuth)      	
	{ 
	 
	 double dec_angle_radian=toRadian(decl_angle( s1,aygunu));
	 double enlem_radian=toRadian(enlem);
	 double incline_angle_radian=toRadian(incline_angle);
	 double surf_azimuth_radian=toRadian(surf_azimuth);
	 double A=Math.sin(dec_angle_radian)*Math.sin(enlem_radian)*Math.cos(incline_angle_radian);
	 double B=-Math.sin(dec_angle_radian)*Math.cos(enlem_radian)*Math.sin(incline_angle_radian)*Math.cos(surf_azimuth_radian);
	 double C=Math.cos(dec_angle_radian)*Math.cos(enlem_radian)*Math.cos(incline_angle_radian);
	 double D=Math.cos(dec_angle_radian)*Math.sin(enlem_radian)*Math.sin(incline_angle_radian)*Math.cos(surf_azimuth_radian);
	 double E=Math.cos(dec_angle_radian)*Math.sin(incline_angle_radian)*Math.sin(surf_azimuth_radian)	;	
	 double K=A-B;
	 double L=C+D;
	 double M=L*L+E*E;
	 double N=2*K*L;
	 double P=K*K-E*E;
	 double ws1=-Math.tan(enlem_radian)*Math.tan(dec_angle_radian);
	 double ws=Math.acos(ws1);
	 double ws_degree=ws*(180/Math.PI);
	 double dummy1=((-N-Math.sqrt(N*N-4*M*P))/(2*M));
	 double wb1=Math.acos(dummy1);
	 double dummy2=((-N-Math.sqrt(N*N-4*M*P))/(2*M));
	 double wb2=Math.acos(dummy2);
	 double wb1_degree=wb1*(180/Math.PI);
	 double wb2_degree=wb2*(180/Math.PI);
	  // wbb[0]=güneş doğuş açısı
	  // wbb[1]=güneş batış açısı     	
	 double[] wbb=new double[2];
	 if (surf_azimuth==0)
	 {wbb[0]=Math.max(-wb1_degree,ws_degree);
	  wbb[1]=Math.min(wb1_degree,ws_degree);}
	 else if((surf_azimuth<0)&&(surf_azimuth>-180))
	 {wbb[0]=Math.max(-Math.max(wb1_degree,wb2_degree),ws_degree);
	  wbb[1]=Math.min(Math.min(wb1_degree,wb2_degree),ws_degree);
	 }
	 else if((surf_azimuth>0)&&(surf_azimuth<180))
	 {wbb[0]=Math.max(-Math.min(wb1_degree,wb2_degree),-ws_degree);
	  wbb[1]=Math.min(Math.max(wb1_degree,wb2_degree),ws_degree);	 
	 }
	 else if((surf_azimuth==180)||(surf_azimuth==0))
	 {wbb[0]=Math.max(-wb1_degree,-ws_degree);
	  wbb[1]=Math.min(wb1,ws);
	 }
	
     return wbb;	 
 }
	   
           public static void main(String args[])
       {
	       double enlem=38; // in degree
	       String s1="Nisan"; 
	       
	       double aygunu=14;
	       double rakım=27*(1E-03); // in km 
	       String s2="Orta_enlem_yaz";
	       double incline_angle=0.0001;//in degree
	       double surf_azimuth=0;//in degree
	       boolean b11=false;
	       boolean b12=false;
	       double saatfarkı2=3.0;
	       double saatfarkı1=4.0;
	       double[] sunsetrise=sun_set_rise(s1, aygunu,enlem,incline_angle,surf_azimuth);      	
	       double ws=sunsetrise[1];// güneşin batış açısı 
	       double Ho=daily_solar_radiation(enlem,s1,aygunu,incline_angle,surf_azimuth); // in MJ/m2
	       double I0=hour_period_solar_radiation(enlem,s1,aygunu,b11,b12,saatfarkı1,saatfarkı2);//
	       
	       double zen_rd1=horiz_zenith_angle(enlem,s1,aygunu,false,0.5);
	       double rb1=beam_radiation_transmittance(s2,rakım,enlem,s1,aygunu, false,0.5) ;
	       double Gon=G_on(s1,aygunu)*3600/((1E06));
	       double Gcnb1=Gon*rb1;
	    
	       double Gcb1=Gon*rb1*Math.cos(zen_rd1);
	       double rd1=diffuse_radiation_transmittance(s2,rakım,enlem,s1,aygunu,false,0.5) ;
	       double Gcd1=Gon*Math.cos(zen_rd1)*rd1;
	       double Gctop1=Gcb1+Gcd1;
	       //*********************************************************************
	       double zen_rd2=horiz_zenith_angle(enlem,s1,aygunu,false,1.5);
	       double rb2=beam_radiation_transmittance(s2,rakım,enlem,s1,aygunu, false,1.5) ;
	       double Gon2=G_on(s1,aygunu)*3600/((1E06));
	       double Gcnb2=Gon*rb2;
	    
	       double Gcb2=Gon*rb2*Math.cos(zen_rd2);
	       double rd2=diffuse_radiation_transmittance(s2,rakım,enlem,s1,aygunu,false,1.5) ;
	       double Gcd2=Gon*Math.cos(zen_rd2)*rd2;
	       double Gctop2=Gcb2+Gcd2;
	       //************************************************************************
	       double zen_rd3=horiz_zenith_angle(enlem,s1,aygunu,false,2.5);
	       double rb3=beam_radiation_transmittance(s2,rakım,enlem,s1,aygunu, false,2.5) ;
	       double Gon3=G_on(s1,aygunu)*3600/((1E06));
	       double Gcnb3=Gon*rb3;
	    
	       double Gcb3=Gon*rb3*Math.cos(zen_rd3);
	       double rd3=diffuse_radiation_transmittance(s2,rakım,enlem,s1,aygunu,false,2.5) ;
	       double Gcd3=Gon*Math.cos(zen_rd3)*rd3;
	       double Gctop3=Gcb3+Gcd3;
	       //*************************************************************************
	       double zen_rd4=horiz_zenith_angle(enlem,s1,aygunu,false,3.5);
	       double rb4=beam_radiation_transmittance(s2,rakım,enlem,s1,aygunu, false,3.5) ;
	       double Gon4=G_on(s1,aygunu)*3600/((1E06));
	       double Gcnb4=Gon*rb4;
	    
	       double Gcb4=Gon*rb4*Math.cos(zen_rd4);
	       double rd4=diffuse_radiation_transmittance(s2,rakım,enlem,s1,aygunu,false,3.5) ;
	       double Gcd4=Gon*Math.cos(zen_rd4)*rd4;
	       double Gctop4=Gcb4+Gcd4;
	       //***************************************************************************
	       double zen_rd5=horiz_zenith_angle(enlem,s1,aygunu,false,4.5);
	       double rb5=beam_radiation_transmittance(s2,rakım,enlem,s1,aygunu, false,4.5) ;
	       double Gon5=G_on(s1,aygunu)*3600/((1E06));
	       double Gcnb5=Gon*rb5;
	    
	       double Gcb5=Gon*rb5*Math.cos(zen_rd5);
	       double rd5=diffuse_radiation_transmittance(s2,rakım,enlem,s1,aygunu,false,4.5) ;
	       double Gcd5=Gon*Math.cos(zen_rd5)*rd5;
	       double Gctop5=Gcb5+Gcd5;
	       //****************************************************************************
	       double zen_rd6=horiz_zenith_angle(enlem,s1,aygunu,false,5.5);
	       double rb6=beam_radiation_transmittance(s2,rakım,enlem,s1,aygunu, false,5.5) ;
	       double Gon6=G_on(s1,aygunu)*3600/((1E06));
	       double Gcnb6=Gon*rb6;
	    
	       double Gcb6=Gon*rb6*Math.cos(zen_rd6);
	       double rd6=diffuse_radiation_transmittance(s2,rakım,enlem,s1,aygunu,false,5.5) ;
	       double Gcd6=Gon*Math.cos(zen_rd6)*rd6;
	       double Gctop6=Gcb6+Gcd6;
	       //****************************************************************************
	       double zen_rd7=horiz_zenith_angle(enlem,s1,aygunu,false,6.5);
	       double rb7=beam_radiation_transmittance(s2,rakım,enlem,s1,aygunu, false,6.5) ;
	       double Gon7=G_on(s1,aygunu)*3600/((1E06));
	       double Gcnb7=Gon*rb7;
	    
	       double Gcb7=Gon*rb7*Math.cos(zen_rd7);
	       double rd7=diffuse_radiation_transmittance(s2,rakım,enlem,s1,aygunu,false,6.5) ;
	       double Gcd7=Gon*Math.cos(zen_rd7)*rd7;
	       double Gctop7=Gcb7+Gcd7;
	       
	       
	       
	       
	       double totalrad= 2*(Gctop1+Gctop2+Gctop3+Gctop4+Gctop5+Gctop6+Gctop7);
	       //double total_radiation=27.92*(1E6);//J/m2
	       //double Kt=total_radiation/Ho;
	       //double Hd_Collares=Collares(Kt);
	       //double Hd_Erbs=Erbs(Kt,ws);
	       //double H_diffuse_herbs=Hd_Erbs*total_radiation;
	       //double H_diffuse_Collares=Hd_Collares*total_radiation;
	       // System.out.println("beam rad trans") 
	        System.out.println("*********************************************************************************************************************************");
	        System.out.println("* -- Hours-------Beam rad trans----Beam radiation-----Hrzntl pl rad------Diff rad trans-----Diffuse radiation-----Total radiation");
	        System.out.println("* 11-12//12-13 "+rb1+"  "+Gcnb1+"  "+Gcb1+"  "+rd1+"   "+Gcd1+"   "+Gctop1);
	        System.out.println("* 10-11//13-14 "+rb2+"  "+Gcnb2+"  "+Gcb2+"  "+rd2+"   "+Gcd2+"   "+Gctop2);
	        System.out.println("* 9-10//14-15  "+rb3+"  "+Gcnb3+"  "+Gcb3+"  "+rd3+"   "+Gcd3+"   "+Gctop3);
	        System.out.println("* 8-9//15-16   "+rb4+"  "+Gcnb4+"  "+Gcb4+"  "+rd4+"   "+Gcd4+"   "+Gctop4);
	        System.out.println("* 7-8//16-17   "+rb5+"  "+Gcnb5+"  "+Gcb5+"  "+rd5+"   "+Gcd5+"   "+Gctop5);
	        System.out.println("* 6-7//17-18   "+rb6+"  "+Gcnb6+"  "+Gcb6+"  "+rd6+"   "+Gcd6+"   "+Gctop6); 
	        System.out.println("* 5-6//18-19   "+rb7+"  "+Gcnb7+"  "+Gcb7+"  "+rd7+"   "+Gcd7+"   "+Gctop7);
	        System.out.println("*************************************************************************************************************************************");
	        System.out.println("Daily total radiation=   "+totalrad +"   MJ/m2");
	        
	        //System.out.println(H_diffuse_Collares+" "+H_diffuse_herbs);
	       
	       
	   }






}


