H_deg = [ 6 5 4 3 2 1 0 -1 -2 -3 -4 -5 ].*15; %the month is November
L = 31.77; %in jerusalem
r = 0.18;
Ic_monthly = 0;

for  n = 305:334
zeta = 23.45.*(sind((360/365).*(n-81))); % solar declination angle
beta = asind((cosd(L)).*(cosd(zeta)).*(cosd(H_deg)) + (sind(L)).*(sind(zeta))); % altitude angle
azimuh = asind(((cosd(zeta)).*(sind(H_deg)))./(cosd(beta))); %sun's azimuth angle
azi = 180 - azimuh;
check = tand(zeta)./tand(L);
if(cosd(H_deg) <= check )
azimuh = azi;
end 
m = 1./sind(beta); % air mass ratio
a = 1160 + 75.*(sind((360/365).*(n-275))); % apparent extraterrestrial flux
k = 0.174 + 0.035.*(sind((360/365).*(n-100))); % optical depth

Ib = a.*exp(-k*m);
q = cosd(beta).*cosd(azimuh - 0 ).*sind(L) + sind(beta).*cosd(L); % cos(theta)
Ibc = Ib.*q % direct bean radiation
C = 0.095 + 0.04.*(sind((360/365).*(n-100))); % sky diffuse factor
Idc = Ib.*C.*((1 + cosd(L))/2) % diffuse radiation
Irc = Ib.*r.*(C + sind(beta)).*((1 - cosd(L))/2) % reflected radiation
Ic = Ibc + Idc + Irc % total radiation
Ic_monthly = Ic_monthly + Ic;

plot(((H_deg)/15)+12,Ic)
xlabel('Hours/(1-24)');
ylabel('IC(solar Irradiance in Watt/m^2)');
title('Solar Irradiance all November ');
hold on

end 
PSH = sum(Ic_monthly);
PSH = PSH/1000/30

IC_monthly_avg = sum(Ic_monthly/30)
IC_daily_avg = (IC_monthly_avg)/13 




