//====================================================//
// This program is developed as a part of solar       //
// energy aplications in Polser Limited, İzmir-Turkey //
// Reference : Solar Engineering of Thermal Processes //
// John A. Duffie & William A. Beckman                //
// Third edition John Willey and Sons                 //
// ISBN: 13 978-0-471-69867-8                         //
// program developer : Dr. M. Turhan Coban            //
//====================================================//

import java.util.*;
import javax.swing.*;
import static java.lang.Math.*;
 
public class solar extends Mathd
{	
// Solar constant  page 5
public static final double Gsc=1367.0;//W/m^2;

//  extraterresterial radiation incident on the plane 
//  normal to the radiation on the nth day of the year
//  page 9
public static double Gon(int n)
{double g=solar.Gsc*(1.0+0.033*cosd(360.0*n/365.0));
return g;
}

// page 9:  1.4.1a
public static double Gon(double n)
{double g=solar.Gsc*(1.0+0.033*cosd(360.0*n/365.0));
return g;
}
// page 9:  1.4.1a
public static double Gon(date calendar)
{   double n=n(calendar);
	double g=solar.Gsc*(1.0+0.033*cosd(360.0*n/365.0));
    return g;
}

//  Go for a horizontal surface at anytime between sinrise and sunset W/m^2
public static double G0tetaz(date calendar,double fi,double gamma)
{ double n=n(calendar);
  double G0=Gon(n);
  double tetaz=tetaz(calendar,fi,gamma);
  double Gotetaz=G0*cosd(tetaz);
  return Gotetaz;
}
//  Ho for a horizontal surface at anytime between sinrise and sunset
//  J/m^2 daily total radiation
public static double H0tetaz(date calendar,double fi,double gamma)
{ double n=n(calendar);
  double delta=delta(calendar);
  double ws=wss(calendar,fi);
  double H0=24.0*3600*Gsc/PI*(1.0+0.033*cosd(360.0*n/365.0))*
  (cosd(fi)*cosd(delta)*sind(ws)+PI*ws/180.0*sind(fi)*sind(delta));
  //System.out.println("n="+n+"delta="+delta+"ws="+ws+"H0="+H0);
  return H0;
}

// Io for a horizontal surface between two given times
//  W/m^2 daily total enegy input
public static double I0tetaz(date calendar1,date calendar2,double fi,double gamma)
{ double n=n(calendar1);
  double delta=delta(calendar1);
  double ws1=w(calendar1);
  double ws2=w(calendar2);
  System.out.println("n="+n+"ws1="+ws1+"ws2="+ws2);
  double I0=12.0*3600*Gsc/PI*(1.0+0.033*cosd(360.0*n/365.0))*
  (cosd(fi)*cosd(delta)*(sind(ws2)-sind(ws1))+PI*(ws2-ws1)/180.0*sind(fi)*sind(delta));
  System.out.println("I0="+I0);
  return I0;
}
// page 9 : eqn 1.4.2  
public static double B(double n)
{   // derece
	double g=360.0/365.0*(n-1);
    return g;
}
//air mass page 10 :equation 1.5.1 and Kasten and Young equation(1989) combined
public static double m(double tetaz,double h)
{  // tetaz zenith angle degree
   // h altitude
   // m air mass
   double g;
   if(tetaz<=70) 	
	 g=1.0/cosd(tetaz);
   else
     g=Math.exp(-0.0001184*h)/(cosd(tetaz)+0.5057*Math.pow((96.080-tetaz),-1.634));	
   return g;
}

//The day of the year
public static double n(date calendar)
{int nn[]={0,31,59,90,120,151,181,212,243,273,304,334};
int day=calendar.gc.get(Calendar.DAY_OF_MONTH);	
int month=calendar.gc.get(Calendar.MONTH);
int hour=calendar.gc.get(Calendar.HOUR_OF_DAY);
int minute=calendar.gc.get(Calendar.MINUTE);
double n=day+nn[month]+hour/24.0+minute/(24*60);
return n;
}

public static String direction_beam_radiation()
{String s="";
 
 char fi='\u03A6';
 char smaller='\u2264';
 char derece='\u00B0';
 char delta='\u03B4';
 char beta='\u03B2';
 char gamma='\u03B3';
 char omega='\u03C9';
 char theta='\u03B8';
 char alpha='\u03B1';
 s+=" "+fi+" Latitude, the angular location north or south of equator, north positive -90"+derece+smaller+fi+smaller+"90"+derece+"\n\n";
 s+=" "+delta+" Declination, the angular position of the sun at solar noon (i.e. when the sun is on the local meridian)\nwith respect to plane of the equator -23.45"+derece+smaller+delta+smaller+"23.45"+derece+"\n\n";
 s+=" "+beta+" Slope, the angle between the plane of the surface in question and the horizontal 0"+derece+smaller+beta+smaller+"180"+derece+"\n("+beta+">90"+derece+" means that the surface has a downward component.)\n\n";
 s+=" "+gamma+" Surface azimuth angle, the deviation of the projection on a horizontal plane of the normal to the \nsurface from the local meridian, with zero due south, east negative and west positive -180"+derece+smaller+gamma+smaller+"180"+derece+"\n("+beta+">90"+derece+" means that the surface has a downward component.)\n\n";
 s+=" "+omega+" hour angle, the angular displacement of the sun east or west of the local meridian due to rotation \nof the earth on its axis at 15"+derece+"an hour; morning negative, afternoon positive\n\n";
 s+=" "+theta+" angle of incidence, the angle betwen the beam radiation on the surface and the normal to that surface\n\n";
 s+=" "+theta+"z Zenith angle, the angle between the vertical and the line of the sun, that is the angle of incidence\n of beam radiation on a horizontal surface\n\n";
 s+=" "+alpha+"s, Solar altitude angle, the angle between the horizontal and the line to the sun, that is compliment\n of the zenith angle\n\n";
 s+=" "+gamma+"s, Solar azimuth angle, the angular displacement from south of the projection of beam radiation on \nthe horizontal plane. Displacements east of south are negative and west of south are positive\n\n";
 return s;
}

public static void print_direction_beam_radiation()
{   ImageIcon picture=new ImageIcon("sun.jpg");
	JOptionPane.showMessageDialog(null,direction_beam_radiation(),"Direction of Beam Radiation",JOptionPane.INFORMATION_MESSAGE,picture);}

//solar time general
public static double solar_time(date calendar, double Lloc,double Lstandard)
{int nn[]={0,31,59,90,120,151,181,212,243,273,304,334};
int day=calendar.gc.get(Calendar.DAY_OF_MONTH);	
int month=calendar.gc.get(Calendar.MONTH);
int hour=calendar.gc.get(Calendar.HOUR_OF_DAY);
int minute=calendar.gc.get(Calendar.MINUTE);
double n=day+nn[month]+hour/24.0+minute/(24*60);
double B=B(n);
double E=229.2*(0.000075+0.001868*cosd(B)-0.032077*sind(B)-0.014615*cosd(2.0*B)-0.04089*sind(2.0*B));
double a=4.0*(Lstandard-Lloc)+E;
return a;
}

//solar time in turkey
public static double solar_time(date calendar, double Lloc)
{ 
int nn[]={0,31,59,90,120,151,181,212,243,273,304,334};
double Lstandard=30.0;
int day=calendar.gc.get(Calendar.DAY_OF_MONTH);	
int month=calendar.gc.get(Calendar.MONTH);
int hour=calendar.gc.get(Calendar.HOUR_OF_DAY);
int minute=calendar.gc.get(Calendar.MINUTE);
double n=day+nn[month]+hour/24.0+minute/(24*60);
double B=B(n);
double E=229.2*(0.000075+0.001868*cosd(B)-0.032077*sind(B)-0.014615*cosd(2.0*B)-0.04089*sind(2.0*B));
double a=4.0*(Lstandard-Lloc)+E;
 return a;
}

// solar time in turkey now
public static double solar_time(double Lloc)
{  
int nn[]={0,31,59,90,120,151,181,212,243,273,304,334};
date calendar=new date();
double Lstandard=30.0;
int day=calendar.gc.get(Calendar.DAY_OF_MONTH);	
int month=calendar.gc.get(Calendar.MONTH);
int hour=calendar.gc.get(Calendar.HOUR_OF_DAY);
int minute=calendar.gc.get(Calendar.MINUTE);
double n=day+nn[month]+hour/24.0+minute/(24*60);
double B=B(n);
double E=229.2*(0.000075+0.001868*cosd(B)-0.032077*sind(B)-0.014615*cosd(2.0*B)-0.04089*sind(2.0*B));
double a=4.0*(Lstandard-Lloc)+E;
return a;
}

//solar time in any location now!
public static double solar_time(double Lloc,double Lstandard)
{  
int nn[]={0,31,59,90,120,151,181,212,243,273,304,334};
date calendar=new date();
int day=calendar.gc.get(Calendar.DAY_OF_MONTH);	
int month=calendar.gc.get(Calendar.MONTH);
int hour=calendar.gc.get(Calendar.HOUR_OF_DAY);
int minute=calendar.gc.get(Calendar.MINUTE);
double n=day+nn[month]+hour/24.0+minute/(24*60);
double B=B(n);
double E=229.2*(0.000075+0.001868*cosd(B)-0.032077*sind(B)-0.014615*cosd(2.0*B)-0.04089*sind(2.0*B));
double a=4.0*(Lstandard-Lloc)+E;
 return a;
}

//declination angle
public static double delta(date calendar)
{  
int nn[]={0,31,59,90,120,151,181,212,243,273,304,334};
int day=calendar.gc.get(Calendar.DAY_OF_MONTH);	
int month=calendar.gc.get(Calendar.MONTH);
int hour=calendar.gc.get(Calendar.HOUR_OF_DAY);
int minute=calendar.gc.get(Calendar.MINUTE);
double n=day+nn[month]+hour/24.0+minute/(24*60);
double B=B(n);
//double d=(0.00691-0.399912*cosd(B)+0.070257*sind(B)-0.006758*cosd(2.0*B)+0.000907*sind(2.0*B))-0.002679*cosd(3.0*B)+0.00148*sind(3.0*B);
double d1=23.45*sind(360*(284+n)/365.0); 
return d1;
}
//hour angle
// noon 0
// morning -
// afternoon positive
// 15 degree per hour
public static double w(date calendar)
{
int day=calendar.gc.get(Calendar.DAY_OF_MONTH);	
int month=calendar.gc.get(Calendar.MONTH);
int hour=calendar.gc.get(Calendar.HOUR_OF_DAY);
int minute=calendar.gc.get(Calendar.MINUTE);
double h=hour+(double)(minute/60.0)-12.0;
h*=15.0;
System.out.println("day="+day+"month="+month+"hour="+hour+"minute="+minute+"h="+h);
return h;
}
//angle of incidence
public static double teta(date calendar,double beta,double fi,double gamma)
{
// fi latitude 90: north pole, 0 equator, -90: south pole
// beta slope : angle between horizontal surface and the plane of the surface
// gamma : surface azimuth angle: angle from south
//         east minus; west : plus
// w     : hours angle 
double delta=delta(calendar);
double w=w(calendar);
System.out.println("delta="+delta+"w="+w);
double costeta=sind(delta)*sind(fi)*cosd(beta)-sind(delta)*cosd(fi)*sind(beta)*cosd(gamma)+
cosd(delta)*cosd(fi)*cosd(beta)*cosd(w)+cosd(delta)*sind(fi)*sind(beta)*cosd(gamma)*cosd(w)+
cosd(delta)*sind(beta)*sind(gamma)*sind(w); 	
return acosd(costeta);
}
// zenith angle
public static double tetaz(date calendar,double fi,double gamma)
{
// fi latitude 90: north pole, 0 equator, -90: south pole
// beta slope : 0 angle between horizontal surface and the plane of the surface
// gamma : surface azimuth angle: angle from south
//         east minus; west : plus
// w     : hours angle 
return teta(calendar,0,fi,gamma);
}

public static int sign(double x)
{return (int)(x/abs(x)+0.0001);}

//solar azimuth angle
public static double gammas(date calendar,double fi,double gamma)
{
double w=w(calendar);
double tetaz=tetaz(calendar,fi,gamma);
double delta=delta(calendar);
double cosx=(cosd(tetaz)*sind(fi)-sind(delta))/(sind(tetaz)*cosd(fi));
double x=acosd(cosx);
double gammas=sign(w)*abs(x);
return gammas;
}

//sunset hour angle
public static double wss(date calendar,double fi)
{double delta=delta(calendar);
 double ws=acosd(-sind(fi)*sind(delta)/(cosd(fi)*cosd(delta))); 
 return ws;
}

//sunset hour
public static double sunset_double(date calendar,double fi)
{double delta=delta(calendar); 
 double ws=wss(calendar,fi)/15.0+12.0;
 return ws;
}

// sunset hour
public static date sunset_date(date calendar,double fi)
{double delta=delta(calendar); 
 double ws=wss(calendar,fi)/15.0+12.0;
 int hour=(int)ws;
 int minute=(int)((ws-hour)*60.0);
 int day=calendar.gc.get(Calendar.DAY_OF_MONTH);
 int year=calendar.gc.get(Calendar.YEAR)+1900;	
 int month=calendar.gc.get(Calendar.MONTH)+1;
 date d=new date(day,month,year,hour,minute);
 return d;
}

//sunrise hour angle
public static double wsr(date calendar,double fi)
{ return -wss(calendar,fi);}

//sunrise hour
public static double sunrise(date calendar,double fi)
{double delta=delta(calendar); 
 double ws=wsr(calendar,fi)/15.0+12.0;
 return ws;
}

public static date sunrise_date(date calendar,double fi)
{double delta=delta(calendar); 
 double ws=wsr(calendar,fi)/15.0+12.0;
 int hour=(int)ws;
 int minute=(int)((ws-hour)*60.0);
 int day=calendar.gc.get(Calendar.DAY_OF_MONTH);
 int year=calendar.gc.get(Calendar.YEAR)+1900;	
 int month=calendar.gc.get(Calendar.MONTH)+1;
 date d=new date(day,month,year,hour,minute);
 return d;
}

//Number of daylight hours
public static double N_hour(date calendar,double fi)
{double delta=delta(calendar);
 double nh=2.0/15.0*acosd(-tand(fi)*tand(delta));
 return nh;
}

//Solar altitude angle
// 90-tetaz
public static double alfastetaz(date calendar,double fi,double gamma)
{
 return 90-tetaz(calendar,fi,gamma);
}

// Ratio of beam radiation on tilted surface to that on horizontal surface
public double Rb(date calendar,double beta,double fi,double gamma)
{return teta(calendar,beta,fi,gamma)/tetaz(calendar,fi,gamma);}

// TOb : athmospheric scattering factor 
// A  altitude  km
public static double TOb(String climatetype,double A,date calendar,double fi,double gamma)
{
double r0,r1,rk;
String ce[]={"Tropical","Midlatitude summer","Subarctic summer","Midlatitude winter"};
String ct[]={"Tropik","ortaenlem yaz","kutup yaz","ortaenlem kış"};
if(climatetype.equals(ct[0])||climatetype.equals(ce[0]))
{r0=0.95;r1=0.98;rk=1.02;}
else if(climatetype.equals(ct[1])||climatetype.equals(ce[1]))
{r0=0.97;r1=0.99;rk=1.02;}	
else if(climatetype.equals(ct[2])||climatetype.equals(ce[2]))
{r0=0.99;r1=0.99;rk=1.01;}	
else if(climatetype.equals(ct[3])||climatetype.equals(ce[3]))
{r0=1.03;r1=1.01;rk=1.0;}
else
{System.out.println("Error wrong climate type");
r0=1.0;r1=1.0;rk=1.0;
}
double a0s=0.4237-0.00821*(6-A)*(6-A);
double a1s=0.5055+0.00595*(6.5-A)*(6.5-A);
double ks=0.2711+0.01858*(2.5-A)*(2.5-A);
double a0=a0s*r0;
double a1=a1s*r1;
double k=ks*rk;
double tetaz=tetaz(calendar,fi,gamma);
double to=a0+a1*exp(-k/Mathd.cosd(tetaz));
return to;
}

// Ratio of diffuse sky radiation to total radiation
// Erbs at all correlation
public static double Id_over_I(double kt)
{double y=0;
 if(kt<=0.22) {y=1.0-0.09*kt;}
 else if(kt>0.22 && kt<=0.8) {y=0.9511-0.1604*kt+4.388*kt*kt-16.638*kt*kt*kt+12.336*kt*kt*kt*kt; }
 else if(kt>0.8 & kt<=1.0) {y=0.165;}
 else y=0;
 return y;
}

// Beam and diffuse components of daily radiation
// Collares-Pereira and Rabl, Solar Energy,22 155(1979a) "The Average Distribution of Solar 
// Radiation-Corrolations Between Diffuse and emispherical and Between Daily and Hourly 
// Insolation Values ws sunset hour angle
public static double Hd_over_H(double kt)
{double y=0;
 if(kt<=0.17) {y=0.99;}
 else if(kt>0.17 && kt<=0.75) {y=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.8) {y=-0.54*kt+0.632;}
 else if(kt>0.8 & kt<=1.0) {y=0.2;}
 else y=0;
 return y;
}

// beam and diffuse components of daily radiation
//page 78
public static double Hd_over_H(double kt,double ws)
{	double y=0;
 if(ws<=81.4)
 {
 if(kt<=0.715) {y=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 && kt<=1.0) {y=0.143; }
 else y=0;
 }
 else
 {
 if(kt<=0.722) {y=1.0+0.2832*kt-2.5557*kt*kt+0.8448*kt*kt*kt;}
 else if(kt>0.722 && kt<=1.0) {y=0.175;}
 else y=0;
 }
 return y;
}
 
public static void main(String arg[])
{   int year=2014;
    String month="April";
    int day=15;
    int hour1=13;
    int minute=25;
    int hour2=16;
    date gd1=new date(day,month,year,hour1,minute);
    double G0=Gon(gd1);
    double fi=32.75;//latitude angle
    double gamma=0.0;
    double ws=wss(gd1,fi);
    System.out.println("ws="+ws); 
    double H0=H0tetaz(gd1,fi,gamma);
    double H=2.3e7;
    double KT=H/H0;
    System.out.println("KT="+KT+" H0 ="+H0);
    double HDOH=Hd_over_H(KT,ws);
    System.out.println(" Hd/H ="+HDOH);
    System.out.println("sunset hour : "+sunset_date(gd1,fi).toStringT());
    System.out.println("sunrise hour : "+sunrise_date(gd1,fi).toStringT());
    print_direction_beam_radiation();
}

}