Discussion :

[ThomasAU]

Bonjour

Je cherche à mettre en équation la trajectoire apparente du soleil pendant une journée.

Pour cela j'ai ce code :

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
 
clf;
clear all;
putenv('GNUTERM','x11');
 
jour = 100;
heure = linspace(0,24,100);
 
[declinaison, angleHoraire, hauteur, azimut] = Soleil(46, jour, heure);
 
%figure(1);
%plot(heure, hauteur);
%figure(2);
%plot(heure, azimut);
 
figure(3);
hold on;
 
l = cos(hauteur);
plot3(l.*cos(azimut), l.*sin(azimut), sin(hauteur));
xlabel('x');
ylabel('y');
zlabel('z');
 
 
[declinaison, angleHoraire, hauteur, azimut] = Soleil(46, jour, 12);
 
plot3([0 0.5*cos(azimut)], [0 sin(azimut)], [0 0], 'r');

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
 
function [declinaison, angleHoraire, hauteur, azimut] = Soleil(latitude, jour, heure)
 
latitude = latitude * pi/180;
declinaison = declinaison(jour);
angleHoraire = angleHoraire(heure);
 
hauteur = asin( sin(latitude)*sin(declinaison) + cos(latitude)*cos(declinaison)*cos(angleHoraire));
 
azimut = asin( cos(declinaison)*sin(angleHoraire)./cos(hauteur) );

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
 
function dec = declinaison(jour)
 
sinDec = sin(23.45*pi/180)*sin(2*pi/365.25*(jour-81));
dec = 180/pi * asin(sinDec);

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
 
function Ah = angleHoraire(tSV)
 
Ah = pi*(tSV/12 - 1);

J'obtient bien une ellipse incliné, mais elle me parait étrange.
Pouvez-vous y jeter un coup d'œil et m'aider à trouver l'erreur ?

Merci et Bonne Année

[Merel]

Salut,

Qu'est ce qui te parait étrange dans le résultat?
J'obtiens bien une eclipse inclinée qui parait correcte.

[ThomasAU]

Bonjour

En faite, c'est complètement faux.
Il y avait déjà une erreur d'unité dans la fonction "declinaison"

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
 
function dec = declinaison(jour)
 
sinDec = sin(23.45*pi/180)*sin(2*pi/365.25*(jour-81));
dec = 180/pi * asin(sinDec);

La déclinaison doit-être en radian et non pas en degré.

Correction

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
 
function dec = declinaison(jour)
 
sinDec = sin(23.45*pi/180)*sin(2*pi/365.25*(jour-81));
dec = asin(sinDec);

Mais dans sa cas j'ai une ellipse qui "rebondit" sur le plan x=0
C'est pas terrible. Il faudrait que l'ellipse puisse passer de l'autre côté.
Alors es une erreur mathématique ou de programmation ?
Je panche plutôt pour un défaut de programmation. Pour les math, je me suis inspiré de plusieurs sites qui donnent tous les même formules.

A noté que chaque courbe est décalé d'une vingtaine de jour.

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
 
figure(1);clf;
figure(2);clf;
 
clear all;
putenv('GNUTERM','x11');
 
 
figure(1);
hold on;
 
figure(2);
hold on;
 
debutJour = 1;
finJour =81*2;
pas = 20;
 
latitude = 46;
 
for jour = debutJour:pas:finJour
 
		heure = linspace(0,24,1000);
 
		[declinaison, angleHoraire, hauteur, azimut] = Soleil(latitude, jour, heure);
 
		figure(1);
		plot(heure, hauteur, 'b');
		plot(heure, azimut, 'g');
 
		figure(2);
		plot3(cos(hauteur).*cos(azimut), cos(hauteur).*sin(azimut), sin(hauteur));
 
endfor
 
 
figure(1);
grid on;
title('Hauteur & azimut');
legend('hauteur', 'azimut');
print('Hauteur et Azimut du Soleil.png','-dpng');
%print('Hauteur et Azimut.png', '-png');
 
figure(2);
xlabel('x');
ylabel('y');
zlabel('z');
print('Trajectoire Apparente du Soleil.png', '-dpng');
%print('Trajectoire dans le ciel.png', '-png');

[Jerome Briot]

Il faut faire attention avec les fonctions trigonométriques inverses...

ASIN n'est définie que sur l'intervalle [-pi/2 pi/2]

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
th = linspace(0,2*pi,9)
 
s = sin(th)
 
th2 = asin(s)

renvoie

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
th =
 
         0    0.7854    1.5708    2.3562    3.1416    3.9270    4.7124    5.4978    6.2832
 
 
s =
 
         0    0.7071    1.0000    0.7071    0.0000   -0.7071   -1.0000   -0.7071   -0.0000
 
 
th2 =
 
         0    0.7854    1.5708    0.7854    0.0000   -0.7854   -1.5708   -0.7854   -0.0000

LA fonction ATAN2 est beaucoup plus robuste

[ThomasAU]

je comprend pour asin, mais je vois pas en quoi atan2 peut m'aider ...

[Jerome Briot]

je comprend pour asin, mais je vois pas en quoi atan2 peut m'aider ...

Je voulais juste dire que la fonction inverse ATAN2 ne pose pas de problème alors qu'avec les autres fonctions trigonométriques inverses, il faut tenir compte du quadrant dans lequel l'angle se trouve

[ThomasAU]

Merci pour vos réponse.

J'ai résolu mon problème avec le code suivant :

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
 
function [declinaison, angleHoraire, hauteur, azimute] = Soleil(latitude, jour, heure)
 
 
latitude = latitude * pi/180;
declinaison = declinaison(jour);
angleHoraire = angleHoraire(heure);
 
hauteur = asin( sin(latitude).*sin(declinaison) + cos(latitude).*cos(declinaison).*cos(angleHoraire));
 
azimute = asin( cos(declinaison).*sin(angleHoraire)./cos(hauteur) );
 
s = size(azimute);
 
for i = s(2)/2:s(2)
 
	if(azimute(i-1) > azimute(i))
		azimute(i) = pi-azimute(i);
	endif
endfor
 
for i = 1:s(2)/2
 
	if(azimute(i) > azimute(i+1))
		azimute(i) = -pi-azimute(i);
	endif
 
endfor