1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56
| clear all
close all
%-----------------Calcul de la Puissance ou énergie------------
n=0:10:365;% le nombre de jour de l'année
phi=45; % la latitude du lieu à Grenoble
gama=0;%angle d'azimut
gamaes=0;%angle exterieur d'azimut
Isc=1367;%Constante Solaire
rho=0.2;%albedo
t0=0;
t=1:100:3600;
i=real(0);
HTIF(37,8)=0;
for i=0:1:4
ii=real(i+1);
beta=15*i;
Y=23.45*sin((pi/180)*360*(284+n)/365);%la declinaison horaire
%omgas=acos(-tan(Y*pi/180)*tan(phi*pi/180))*180/pi;% l'angle horaire du coucher du soleil
omgas=acos(-tan(Y*pi/180)*tan(phi*pi/180));% l'angle horaire du coucher du soleil
A=(86400)*Isc/pi.*(1+0.033*cos(2*pi*n/365.25));%
C=cos(pi*phi/180).*cos(pi*Y/180).*sin(pi*omgas/180);%
%B=(pi/180)*omgas.*sin(pi*phi/180).*sin(pi*Y/180);%
B=omgas.*sin(pi*phi/180).*sin(pi*Y/180);%
H0=(A.*(B+C));% La puissance du rayonnement solaire hors atmosphère
TS=12;% temps solaire vrai
omga=15*(TS-12);%L'angle horaire omega
a=0.409+0.501*sin(pi*omgas/180-pi/3);%la valeur de la cste a de carlos
b=0.6609-0.4767*sin(pi*omgas/180-pi/3);%la valeur de la cste b de carlos
Arg1=(pi/24)*(a+b.*cos(pi*omga/180));
rt=Arg1.*(cos(pi*omga/180)-cos(pi*omgas/180))./(sin(pi*omgas/180)-omgas.*cos(pi*omgas/180));%rapport entre les valeurs horaires et journalières de l'irradiation globale
rd=(pi/24)*(cos(pi*omga/180)-cos(pi*omgas/180))./(sin(pi*omgas/180)-omgas.*cos(pi*omgas/180));%rapport entre les valeurs horaires et journalières de l'irradiation diffuse
H=H0*0.5;%Rayonnement sur une surface horizontale
Hd=H.*(0.0775+0.347*(omgas-pi/2)-(0.505+0.261*(omgas-pi/2))*cos(pi*2*0.5/180-1.8*pi/180));
Hb=H-Hd;%rayonnement global
costetaz=-cos(pi*phi/180).*cos(pi*Y/180).*cos(pi*omga/180)+sin(pi*phi/180).*sin(pi*Y/180);
tetaz=acos(pi*costetaz/180);
costeta=costetaz+sin(pi*tetaz/180)*sin(pi*beta/180)*cos(pi*(gamaes-gama)/180);
Rb=(costeta)/(costetaz);
Hbi=Hb*Rb;
Hdi=Hd*(1+cos(pi*beta/180))/2;
Hr=rho*H*(1-cos(pi*beta/180))/2;
Hti=Hbi+Hdi+Hr;
HTIF(:,ii)=Hti;
end;
Ha=HTIF(:,2);%Energie retenue
Ia=diff(Ha);%Puissance ou l'éclairement
figure
SUBPLOT(3,1,1), plot(n,Y);
title('Declinaison')
SUBPLOT(3,1,2), plot(n,omgas);
title('Omegas')
SUBPLOT(3,1,3), plot(n,H);
title('Rayonnement') |
Partager