Bonjour,
S'il vous plaît, j'ai besoin d'aide. Avec le programme ci-dessous la figure de H qui s'affiche comporte des valeurs négatives. C'est l'énergie, théoriquement on ne peut pas avoir de valeur négative. Merci d'avance;
Code : Sélectionner tout - Visualiser dans une fenêtre à part
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')