//======================================================
// 
// Class roSUN
// Mehmet Aydın SOYDAŞ , öğrenci
// Danışman: Dr. Turhan Coban
// Ege University, School of Engineering, 
// Department of Mechanical Engineering
// Bornova İZMİR - TURKEY
// Program yazılırken kullanılan kaynak: Duffie, J., A., Beckman, W., A.,
//										 "Solar Engineering of Thermal Processes",
//										 John Wiley and Sons, New York, 1991.
//

import javax.swing.JOptionPane;


public class radoSUN
{
	int n,d;
	int m;
	int day_mounth[]={0,31,59,90,120,151,181,212,243,273,304,334};
	double enlem,boylam;
	double B,o;
	double E;  // E zaman düzeltme faktörü
	double st;  // standart time - yerel saat zamanı (dakika)
	double solarT; // solar time - güneş zamanı (dakika)
	double Lst; // standart saat boylamı( Batı boylamları (+), doğu boylamları (-))
	double Lloc; // yerel boylam( Batı boylamları (+), doğu boylamları (-))
	
	
  // *** Açılar *** 
    double w; // hour - saat açısı	
	double fi; // Latitude - enlem açısı; kuzey yarım küre için pozitif, güney yarım küre için negatif (-90'< fi <90')
	double delta; // Declination - Deklinasyon açısı; -günlük hesaplanır
	double tetaZ,costetaZ; // Zeith açısı - Güneş ışınlarının yatay bir düzleme gelişini ölçer.
	double gama=0; // Yüzey azimut açısı
	double beta=0; // Slope - Eğimli bir yüzeyin yatayla yapığı açı
	double teta,costeta;

	double YSZgd,w1;
	double YSZgb,w2;
	double gu;
 
  // * katsayılar *
    double a00,a0;
	double a11,a1;
	double kk,k;
	
	double r0;
	double r1;
	double rk;
	
	double A; // deniz seviyesinden (km) olarak yükseklik
	int iklim; // climate type
	String s_iklim;
	double TOb; // yeryüzündeki doğrudan ışınımın, atmosfer dışındaki ışınıma oranı
	double TOd; // yatay bir yüzeye gelen yaygın ışınım, atmosfer dışındaki ışınıma oranı 
	
	double Gsc=1367; //1370; // (W/m^2) güneş sabiti
	double Gan; // yatay bir yüzeye gelen ışınımın hesaplanması
	double Gcnb; // güneş ışınlarına dik bir yüzey için açık doğrudan ışınım şiddeti
	double Gcb; // yatay bir yüzeye gelen doğrudan ışınım şiddeti
	double Gd;
	double Gtoplam; //yatay yüzeye gelen toplam ışınım miktarı
	double Getoplam; //eğik yüzeye gelen toplam ışınım miktarı
	double Rb; // Eğik yüzeye gelen doğrudan ışınım yatay yüzeye gelen doğrudan ışınıma oranı
	double a,b; // Rb nin hasas modeli için katsayılar
	double q;
	double ad,am; // günlük ortalama ve aylık ortalama ışınımlar
	
	
	public radoSUN(  double enlem, double boylam, double ssboylamı, double yukseklik, int climate)
	{
		
		
		fi=enlem;
		Lloc=boylam;
		Lst=ssboylamı; // standart saat boylamı
		A=yukseklik; // deniz seviyesinden (km) olarak yükseklik
		iklim=climate;
		
		
		if( iklim==1 )  // Tropical
		{r0=0.95; r1=0.98; rk=1.02; s_iklim="Tropical"; } 
			else if( iklim==2 )  //  Midlatitude summer
			 {r0=0.97; r1=0.99; rk=1.02; s_iklim="Midlatitude summer"; }
				else if( iklim==3 )  //  Subarctic summer
				{r0=0.99; r1=0.99; rk=1.01; s_iklim="Subarctic summer"; }
				 else if( iklim==4 ) //  Midlatitude winter
				  {r0=1.03; r1=1.01; rk=1.00; s_iklim="Midlatitude winter"; }
		
		
	}

	public int day_sayac(int day, int mounth)
	{
		d=day;
		m=mounth;
		n=d+day_mounth[m-1];
		
		return n;
	}
	
	public int[] day_sayac(int n)
	{
		int xDay=0;
        int xMounth=0;
			
			for(int i=0; i<=n; i++)
			{
				
				for(int j=0; j<=11; j++)
				{
				if( i==day_mounth[j] )
					{xMounth=(j+1);
					 } 
				 }
					 
				 
				xDay=i-day_mounth[(xMounth-1)];
			}
			
			d=xDay+1;
			m=xMounth;
			
			int k[] = new int[2];
				k[0] = xDay+1;
				k[1] = xMounth;
		
		return k;
	}
	
	public double solar_time(int yearsday, double standartTime)
	{
		st=standartTime; // standart time - yerel saat zamanı (dakika)
		n=yearsday;
		B=((n-1)*(360/365))*(Math.PI/180);
		E=229.2*(0.000075+0.001868*Math.cos(B)-0.032077*Math.sin(B)-0.014615*Math.cos(2*B)-0.04089*Math.sin(2*B));  // dakika cinsinden
		solarT=st+4*(Lst-Lloc)+E;
		
		return solarT;
	}
	
	public double radiation(int nn,double standartTime)
	{
		
		solarT=solar_time(nn,standartTime);
		
		w=15*((solarT/60)-12);  // hour - saat açısı
		
		// deklinasyon açısı hesabı
		delta=23.45*Math.sin(2*Math.PI*(284+n)/365);
		
		// yatay bir yüzey için zenith açısı
		costetaZ=Math.cos(delta*Math.PI/180)*Math.cos(fi*Math.PI/180)*Math.cos(w*Math.PI/180)
						+Math.sin(delta*Math.PI/180)*Math.sin(fi*Math.PI/180);
		tetaZ=Math.acos(costetaZ)*180/Math.PI;				
						
		// güneş doğuş ve batış saatlerinin hesaplanması
		
		double ws=(Math.acos(-Math.tan(delta*(Math.PI/180))*Math.tan(fi*(Math.PI/180))))*(180/Math.PI);
		
		YSZgd=12-((4*(Lst-Lloc)+E)/60)-ws/15;
		YSZgb=12-((4*(Lst-Lloc)+E)/60)+ws/15;	
		
		gu=2*ws/15;
		
		// yatay bir yüzeye gelen ışınımın hesaplanması
		
		Gan=Gsc*(1+0.033*Math.cos(2*Math.PI*nn/365));
		
		a00=0.4237-0.00821*(6-A)*(6-A);
		a11=0.5055+0.00595*(6.5-A)*(6.5-A);
		kk=0.2711+0.01858*(2.5-A)*(2.5-A);	
		
		a0=r0*a00;
		a1=r1*a11;
		k=rk*kk;
		
		
		TOb=a0+a1*Math.exp(-(k/costetaZ));
		
		if ( costetaZ == 0 )
		{ 
			TOb=0;
		}
		
		// güneş ışınlarına dik bir yüzey için açık doğrudan ışınım şiddeti
		Gcnb=Gan*TOb;
		
		// yatay bir yüzeye gelen doğrudan ışınım şiddeti
		Gcb=Gan*TOb*costetaZ;
		
		// yatay bir yüzeye gelen yaygın ışınım şiddeti
		TOd=(0.271-0.294*TOb)*costetaZ;
		
		Gd=Gan*TOd;
		
		//yatay yüzeye gelen toplam ışınım miktarı
		
		Gtoplam=(Gcb+Gd+Math.abs(Gcb+Gd))/2;
		
		//double op=1/0;
		if ( n == 337 ||  n == 333 ||  n == 226 ||  n == 150 ||  n == 70 )
		{ 
			Gtoplam=0;
		}
		
		if ( costetaZ == 0 )
		{ 
			Gtoplam=0;
		}
		
		return Gtoplam;
	}
	
	public double radiation(int day, int mounth , double standartTime)
	{
		d=day;
		m=mounth;
		n=day_sayac(d,m);
		Gtoplam=radiation(n,standartTime);
		
		return Gtoplam;
	}
		
	
	public double angles_surface(int n,double standartTime, double yatay_egim, double azimuth, double ro)
	{
		
		beta=yatay_egim;
		gama=azimuth;
		q=ro;
		
		
		Gtoplam=radiation(n,standartTime);
		
		///*  -> Rb yi bulmak için birincil model <-
		
		costeta=Math.sin(delta*Math.PI/180)*Math.sin(fi*Math.PI/180)*Math.cos(beta*Math.PI/180)-
				Math.sin(delta*Math.PI/180)*Math.cos(fi*Math.PI/180)*Math.sin(beta*Math.PI/180)*Math.cos(gama*Math.PI/180)+
				Math.cos(delta*Math.PI/180)*Math.cos(fi*Math.PI/180)*Math.cos(beta*Math.PI/180)*Math.cos(w*Math.PI/180)+
				Math.cos(delta*Math.PI/180)*Math.sin(fi*Math.PI/180)*Math.sin(beta*Math.PI/180)*Math.cos(gama*Math.PI/180)*Math.cos(w*Math.PI/180)+
				Math.cos(delta*Math.PI/180)*Math.sin(beta*Math.PI/180)*Math.sin(gama*Math.PI/180)*Math.sin(w*Math.PI/180);
		
		teta=Math.acos(costeta)*180/Math.PI;
		
		// Eğik yüzeye gelen doğrudan ışınım yatay yüzeye gelen doğrudan ışınıma oranı
		Rb=(costeta/costetaZ); 
		
		
		
		// radiation on sloped surface: the isotropic diffuse model - by liu and jordan
		
	    Getoplam=Gcb*Rb+Gd*((1+Math.cos(beta*Math.PI/180))/2)+Gtoplam*q*((1-Math.cos(beta*Math.PI/180))/2);
		if( Gtoplam<=0  ) 
		{Getoplam=0;} 
		if( Getoplam<=0  ) 
		{Getoplam=0;} 
		
	
	
		return Getoplam;
		
	}
	
	public double angles_surface(int day, int mounth,double standartTime, double yatay_egim, double azimuth, double ro)
	{
		d=day;
		m=mounth;
		n=day_sayac(d,m);
		Getoplam=angles_surface(n,standartTime, yatay_egim, azimuth,ro);
		
		return Getoplam;
	}
	
	
	
	
	// Concentrating Solar Power - Aktive Panel modeli  =>  Yoğunlaştırılmış Güneş Enerjisi - aktif panel(güneşi takip ediyor)
	
	public double CSP_APmodel(int n, double standartTime, double ro)
	{
		
		solarT=solar_time(n,standartTime);
		w=15*((solarT/60)-12);  // hour - saat açısı
		
		gama=90.0;  //yüzey azimut açısı - Sistem Kuzey-güney ekseninde sabit tutuluyor.
		
		double Gctoplam=angles_surface(n,standartTime,w,gama,ro);
	

		return Gctoplam;
		
	}
	
	// Concentrating Solar Power - Aktive Panel modeli  =>  Yoğunlaştırılmış Güneş Enerjisi - aktif panel(güneşi takip ediyor)
	public double CSP_APmodel(int day, int mounth, double standartTime, double ro)
	{
		
		d=day;
		m=mounth;
		n=day_sayac(d,m);
		
		
		double Gctoplam=CSP_APmodel(n,standartTime, ro);
	

		return Gctoplam;
		
	}
										
	
	
	
	public double average_day_CSP_APmodel(int day, int mounth, double ro)  //CSP_APmodel için günlük ortalama ışınım miktarı(W/m^2 )
	{
		
		double x0=0;
		double toplam=0;
		int oran=2;	
        double	orann=oran;
		
		for(int i=0; i<=(24*oran); i++)
			{
				x0=i*(60/oran);
				
				double t=CSP_APmodel(day,mounth,x0, ro);
				
      		toplam+=t;
 			}
 		double average=toplam/(24*orann);
		
		return average;
	}
	
	public double average_day_CSP_APmodel(int n, double ro)  //CSP_APmodel için günlük ortalama ışınım miktarı(W/m^2 )
	{
		
		double x0=0;
		double toplam=0;
		int oran=2;	
        double	orann=oran;
		
		for(int i=0; i<=(24*oran); i++)
			{
				x0=i*(60/oran);
				
				double t=CSP_APmodel(n,x0, ro);
				
				
      		toplam+=t;
 			}
 		double average=toplam/(24*orann);
		
		return average;
	}

			// Depo edilmesi modeli_ Concentrating Solar Power - Aktive Panel modeli    
			// Yoğunlaştırılmış Güneş Enerjisi - aktif panel(güneşi takip ediyor)
	
			public double CSP_APmodel_Storage(int n, double standartTime, double ro)
			{
				
				//...solarT=solar_time(n,standartTime);
				//...w=15*((solarT/60)-12);  // hour - saat açısı
				//...gama=90.0;  //yüzey azimut açısı - Sistem Kuzey-güney ekseninde sabit tutuluyor.
				double storageTime=standartTime-9.0*60;
				
				if( storageTime<0 )
					{storageTime=24*60-storageTime; } 
				
				double Gc_Stoplam=0.35*CSP_APmodel(n,standartTime, ro)+0.65*average_day_CSP_APmodel(n,ro);
				                
	
				return Gc_Stoplam;
				
			}
			
			// Depo edilmesi modeli_ Concentrating Solar Power - Aktive Panel modeli 
			//Yoğunlaştırılmış Güneş Enerjisi - aktif panel(güneşi takip ediyor)
			public double CSP_APmodel_Storage(int day, int mounth, double standartTime, double ro)
			{
				
				d=day;
				m=mounth;
				n=day_sayac(d,m);
				
				
				double Gc_Stoplam=CSP_APmodel_Storage(n,standartTime, ro); 
			
		
				return Gc_Stoplam;
				
									}
	
	
	public double[] maximum_CSP_APmodel( double ro)  //CSP_APmodel için 365 gün içerisinde maksimum ışınım miktarı(W/m^2 )
	{
		
		double time;
		double t=0;
		int nn=0;
		double kontrol=0;
		
		
		
		for(int n=1; n<=365; n++)
			{
				
			for(int i=0; i<=24; i++)
				{
					time=i*60.0;

					if( kontrol<CSP_APmodel(n,time,ro)  ) 
					{
						kontrol=CSP_APmodel(n,time,ro);
						t=time;
						nn=n;
					} 
				
 				}
			}
			
			
		int xDay=0;
        int xMounth=0;
		
			xDay=day_sayac(nn)[0];
			xMounth=day_sayac(nn)[1];
			
 		double m[] = new double[5];
		m[0] = nn;
		m[1] = t;
		m[2] = kontrol;
		m[3] = xMounth;
		m[4] = xDay;
		
		return m;
	
	}
	
	public double average_year_CSP_APmodel(double ro)  //CSP_APmodel için yıllık ortalama ışınım miktarı(W/m^2 )
	{
		
		double x0=0;
		double toplam=0;
		double t;
		
		int oran=1;	
        double	orann=oran;
		
        for(int n=1; n<=365; n++)
			{
				for(int i=0; i<(24*oran); i++)
					{
						x0=i*(60/orann);
				
						t=CSP_APmodel(n, x0, ro);
			

      					toplam+=t;
 					}
			}
			
		double averagem=toplam/((365-5)*24*oran);  // hesaplamada yılın sadece beş günü için NaN hatası alınmakta
													// o yüzden ortalama hesabında bunlar çıkartılmakdadır.
		return averagem;
	}

	public double average_year_angles_surface(double yatay_egim, double azimuth, double ro)
	  //CSP_APmodel için yıllık ortalama ışınım miktarı(W/m^2 )
	{
		
		double x0=0;
		double toplam=0;
		double t;
		
		int oran=1;	
        double	orann=oran;
		
        for(int n=1; n<=365; n++)
			{
				for(int i=0; i<(24*oran); i++)
					{
						x0=i*(60/orann);
				
						t=angles_surface(n,x0, yatay_egim, azimuth, ro);
			

      					toplam+=t;
 					}
			}
			
		double averagem=toplam/((365-5)*24*oran);  // hesaplamada yılın sadece beş günü için NaN hatası alınmakta
													// o yüzden ortalama hesabında bunlar çıkartılmakdadır.
		return averagem;
	}

	public String cikti()
{

String s="Girdiler(Verilenler) : \n";
s+=" Enlem : "+fi+"' \n";
s+=" Boylam : "+Lloc+"' \n";
s+=" İklim Türü : "+s_iklim+" \n";
s+=" Yerel Saat : "+st/60+" saat \n";
s+=" Standart Saat Boylamı: "+Lst+"' \n";
s+=" Deniz seviyesinden yükseklik : "+A+" km \n";
s+=" .............................. \n";
s+=" Ay : "+m+". \n";
s+=" Gün : "+d+".\n";
s+=" Yüzey Eğim Açısı : "+beta+"' \n";
s+=" Yüzey Azimut Açısı : "+gama+"' \n\n\n";


s+="Sonuçlar : \n";
s+=" Güneşin doğuş saati = "+YSZgd+" saat (yerel saat zamanı ile) \n";
s+=" Güneşin batış saati = "+YSZgb+" saat (yerel saat zamanı ile) \n";
s+=" Gün uzunluğu = "+gu+" saat\n";	
s+=" Güneş ışınlarının Yatay yüzeye Geliş açısı = "+tetaZ +"' (derece)\n";

//s+=" Yatay bir yüzeye gelen doğrudan ışınım şiddeti = "+Gcb +" W/m^2 \n";
//s+=" Yatay bir yüzeye gelen yaygın ışınım şiddeti = "+Gd +" W/m^2 \n";
s+=" Yatay yüzeye gelen toplam ışınım miktarı = "+Gtoplam +" W/m^2 \n";
s+=" Güneş ışınlarının Eğik yüzeye geliş açısı miktarı = "+teta +"' (derece) \n";
s+=" Eğimli yüzeye gelen toplam ışınım miktarı = "+Getoplam +" W/m^2 \n";

return s;
}	
	
	public String cikti2()
{

String s="Sonuçlar\n";
s+=" Güneşin doğuş saati = "+YSZgd+" saat (yerel saat zamanı ile) \n";
s+=" Güneşin batış saati = "+YSZgb+" saat (yerel saat zamanı ile) \n";
s+=" Gün uzunluğu = "+gu+" saat\n";	
s+=" Yatay bir yüzey için zenith açısı = "+tetaZ +" (derece)\n";
s+=" Yatay bir yüzeye gelen doğrudan ışınım şiddeti = "+ Gcb+" W/m^2 \n";
s+=" Yatay bir yüzeye gelen yaygın ışınım şiddeti = "+Gd +" W/m^2 \n";
s+=" Yatay yüzeye gelen toplam ışınım miktarı = "+Gtoplam +" W/m^2 \n";

return s;
}
	
public String ciktiEY()
{

String s="Sonuçlar\n";
s+=" Güneşin doğuş saati = "+YSZgd+" saat (yerel saat zamanı ile) \n";
s+=" Güneşin batış saati = "+YSZgb+" saat (yerel saat zamanı ile) \n";
s+=" Gün uzunluğu = "+gu+" saat\n";	
s+=" Yatay yüzey için zenith açısı = "+tetaZ +" (derece)\n";
s+=" Yatay yüzeye gelen doğrudan ışınım şiddeti = "+Gcb+" W/m^2 \n";
s+=" Yatay yüzeye gelen yaygın ışınım şiddeti = "+Gd +" W/m^2 \n";
s+=" Yatay yüzeye gelen toplam ışınım miktarı = "+Gtoplam +" W/m^2 \n";
s+=" Eğimli yüzeye gelen toplam ışınım miktarı = "+Getoplam +" W/m^2 \n";

return s;
}
}