#################################################
n =172 # day number
Rdc=(360*(pi/180))*((n-1)/365) # day thing

# equation of time
ET=2.2918*(0.00075+0.1868*cos(Rdc)-3.2077*sin(Rdc)-1.4615*cos(2*Rdc)-4.089*sin(2*Rdc))
print(paste0("Equation of time ", ET))

#################################################
LON=-0.45 # local longitude
LST=12 # local standard time
TZ=0 # time zone
LSM=15*TZ # standard meridian longitude

# apparent solar time
AST=LST+(ET/60)+((LON-LSM)/15) 
print(paste0("Apparent solar time ", AST))

#################################################
# solar declination
DEC=23.45*sin(360*((n+284)/365)*(pi/180))
print(paste0("solar declination ", DEC))

# hour angle
H=15*(AST-12) 
print(paste0("hour angle ", H))

L=51.48 # local latitude

# solar altitude angle
B=asin((sin(L*(pi/180))*sin(DEC*(pi/180))+cos(L*(pi/180))*cos(DEC*(pi/180))*cos(H*(pi/180))))
print(paste0("solar altitude angle ", B*(180/pi)))

# solar azimuth
AZ=asin(sin(H*(pi/180))*(cos(DEC*(pi/180))/cos(B*(pi/180))))
print(paste0("solar azimuth angle ", AZ*(180/pi)))
#################################################
# relative air mass
m=1/(sin(B)+0.50572*(6.07995+B)^(-1.6364))
print(paste0("relative air mass ", m))

# optical depths

days=c(0,21,52,80,111,141,172,202,233,264,294,325,355,365)

taubX=c(0.344,0.352,0.361,0.399,0.382,0.376,0.370,0.396,0.383,0.368,0.352,0.337,0.336,0.344)
taub=approx(days,taubX,n = 366)$y[n]

taudX=c(2.508,2.481,2.437,2.272,2.255,2.271,2.286,2.237,2.320,2.935,2.476,2.561,2.535,2.508)
taud=approx(days,taudX,n = 366)$y[n]

print(paste0("taub ", taub))
print(paste0("taud ", taud))

#air mass exponents
b=1.219-((0.043*taub)-(0.151*taud)-(0.204*taub*taud))
d=0.202-((0.852*taub)-(0.007*taud)-(0.357*taub*taud))

Esc = 1367 # solar constant
# extraterrestrial solar irradiance incident on a surface normal
Eo = Esc*(1+0.033*cos((360*((n-3)/365))*(pi/180)))
print(paste0("Extraterrestrial solar irradiance incident on a surface normal ", Eo))

# beam normal irradiance
Eb=Eo*exp(-taub*m^b)
print(paste0("beam normal irradiance ", Eb))

# diffuse horizontal irradiance
Ed=Eo*exp(-taud*m^d)
print(paste0("diffuse horizontal irradiance ", Ed))

# surface azimuth
surfAZ=52*(pi/180)

# surface-solar azimuth angle
ssAZ=surfAZ-AZ

# tilt angle (the angle between the surface and the horizontal plane)
TILT=90*(pi/180)

# angle of incidence
INC=acos(cos(B)*cos(ssAZ)*sin(TILT)+sin(B)*cos(TILT))
print(paste0("angle of incidence ", INC*(180/pi)))

# beam component
if (cos(INC) > 0){
    Etb=Eb*cos(INC) 
} else {
   Etb=0
}

print(paste0("beam component ", Etb))

# ratio of clear-sky diffuse irradiance on a vertical surface to clear-sky diffuse irradiance on the horizontal
Y=max(0.45,0.55+0.437*cos(INC*(pi/180))+0.313*cos(INC*(pi/180))^2)

# diffuse component
if (TILT <= 90*(pi/180)){
    Etd=Ed*(Y*sin(TILT)+cos(TILT)) 
} else {
    Etd=Ed*Y*sin(TILT) 
}
print(paste0("diffuse component ", Etd*(180/pi)))

# ground reflectance
rho=0.2
Etr=(Eb*sin(B)+Ed)*rho*((1-cos(TILT))/2)
print(paste0("ground reflectance ", Etr*(180/pi)))

# total clear-sky irradiance
Et=Etb+Etd+Etr
print(paste0("total clear-sky irradiance ", Et*(180/pi)))
#print(Eb)
#print(Ed)