//======================================================
// 
// 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

import javax.swing.JOptionPane;


public class roSUN
{
	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 azimuth;
	
	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=0; // climate type
	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; // (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 roSUN(  double standartTime,double enlem, double boylam, double ssboylamı, double yukseklik, int iklim)
	{
		
		st=standartTime; // standart time - yerel saat zamanı (dakika)
		fi=enlem;
		Lloc=boylam;
		Lst=ssboylamı; // standart saat boylamı
		A=yukseklik; // deniz seviyesinden (km) olarak yükseklik
		
		if( iklim==1 )  // Tropical
		{r0=0.95; r1=0.98; rk=1.02;} 
			else if( iklim==2 )  //  Midlatitude summer
			 {r0=0.97; r1=0.99; rk=1.02; }
				else if( iklim==3 )  //  Subarctic summer
				{r0=0.99; r1=0.99; rk=1.01; }
				 else if( iklim==4 ) //  Midlatitude winter
				  {r0=1.03; r1=1.01; rk=1.00; }
		
		
	}

	public int day_sayac(int day, int mounth)
	{
		d=day;
		m=mounth;
		n=d+day_mounth[m-1];
		
		return n;
	}
	
	public double solar_time(int yearsday)
	{
		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 n)
	{
		
		solarT=solar_time(n);
		
		w=15*((solarT/60)-12);
		
		// 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((n/365)*2*Math.PI));
		
		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));
		
		// 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;
		
		return Gtoplam;
	}
	
	public double radiation(int day, int mounth )
	{
		d=day;
		m=mounth;
		n=day_sayac(d,m);
		Gtoplam=radiation(n);
		
		return Gtoplam;
	}
		
	
	public double angles_surface(int n, double yatay_egim, double azimuth, double ro)
	{
		
		beta=yatay_egim;
		gama=azimuth;
		q=ro;
		
		
		Gtoplam=radiation(n);
		
		///*  -> 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;} 
		
		Gtoplam=Getoplam;	
	
		return Gtoplam;
		
	}
	
	public double angles_surface(int day, int mounth, double yatay_egim, double azimuth, double ro)
	{
		d=day;
		m=mounth;
		n=day_sayac(d,m);
		Gtoplam=angles_surface(n, yatay_egim, azimuth,ro);
		
		return Gtoplam;
	}
	
	
	
	
	// 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 ro)
	{
		
		solarT=solar_time(n);
		w=15*((solarT/60)-12);
		
		azimuth=90.0;  //yüzey azimut açısı - Sistem Kuzey-güney ekseninde sabit tutuluyor.
		
		Gtoplam=angles_surface(d,m,w,azimuth,ro);
	

		return Gtoplam;
		
	}
	
	public double CSP_APmodel(int day, int mounth, double ro)
	{
		
		d=day;
		m=mounth;
		n=day_sayac(d,m);
		
		
		Gtoplam=CSP_APmodel(n, ro);
	

		return Gtoplam;
		
	}
	
	
	
	public double average_day()  //Yatay yüzeye gelen günlük ortalama ışınım miktarı(W/m^2 )
	{
		
		//double x0;
		
		for(int i=0; i<=24; i++)
			{
				//x0=i*60;
				
				
				//ad+=CSP_APmodel();
				
 			}		
 			ad=ad/24;
		
		return ad;
	}

	public double maximum()  //Yatay yüzeye gelen günlük ortalama ışınım miktarı(W/m^2 )
	{
		
		//double x0;
		
		for(int i=0; i<=24; i++)
			{
				//x0=i*60;

				ad+=Gtoplam;
				
 			}		
 			ad=ad/24;
		
		return ad;
	}


	
	public String cikti()
{

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";
s+=" Yatay yüzeye gelen günlük ortalama ışınım miktarı = "+average_day() +" 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";
s+=" Yatay yüzeye gelen günlük ortalama ışınım miktarı = "+average_day() +" W/m^2 \n";


return s;
}
}