vendredi 3 octobre 2008

Simulation solaire avec povray


Image and video hosting by TinyPic
Une simple maison (maison.pov). (trajectoire du Soleil à l'équinoxe) (à vérifier)

//
// maison.pov
// Fichier povray : description de la scène (éclairage mobile, maison, caméra)
//
#declare Azimut=0; // 0=faîte Est-Ouest
#declare Latitude=50;
global_settings
{
assumed_gamma 1.5
noise_generator 2
}
light_source
{
<0, 0, -5000>, rgb <1, 1, 1>
rotate <0,-clock,0>
rotate <0,Azimut,90-Latitude>
}
//----------------
union
{
prism
{ // maison
0, 12, // z
6, // # pts
<0,0>, <6,0>, <6,3>, <3,6>, <0,3>, <0,0>
rotate <-90, 0, 0>
pigment { color rgb <1, 0, 0> }
}
prism
{ // lucarne
0, 12,
6,
<0,0>, <6,0>, <6,3>, <3,6>, <0,3>, <0,0>
scale <.3, .3, .3>
rotate <-90, 90, 0>
translate <6, 3, -6>
pigment { color rgb <0, 1, 0> }
}
}
camera {
perspective
location <12, 6, 3>
look_at <0, 0, -6>
}


La configuration

//
// Configuration povray : maison.ini
//
Antialias=On
Antialias_Threshold=0.3
Antialias_Depth=3

Input_File_Name=maison.pov

Initial_Frame=1
Final_Frame=32
Initial_Clock=180
Final_Clock=0

Cyclic_Animation=on
Pause_when_Done=off

Un petit script pour générer une .gif animée

#
# maison.sh
# génération des images
povray +ua maison.ini
# assemblage de la .gif animée avec convert (ImageMagick)
convert -delay 20 -loop 0 -dispose Previous maison*.png maison.gif

Avec la version 3.6 de povray, on n'a plus besoin du fichier maison.ini et on peut faire :

#
# maison.sh
# génération des images
povray +ua -KFI0 -KFF30 -KI180.0 -KF0.0 maison.pov
# assemblage de la .gif animée avec convert (ImageMagick)
convert -delay 20 -loop 0 -dispose Previous maison*.png maison.gif

Voir www.Povray.org.
----
voir aussi http://users.skynet.be/chricat/gnomon/gnomon.html
----
Un peu de code 'Python' pour calculer un vecteur solaire dans la base (Ouest, Haut, Sud) en fonction du jour de l'année (1..365) comme argument :

#!/usr/bin/python
#
# sun direction (solar vector) for latitude=50.6
# base = (West, Above, South) as Povray does it
# argv[1] = day of the year 1..365
#
import sys
from Numeric import *
from math import *

def rad(deg): # convert degrees to radians
return deg * math.pi / 180

def h2rad(hour): # convert hours to radians
return (hour -12) * math.pi / 12

def sunh(day): # sun height (day of the year) (in radian)
day -= 81 # origin = march 21
return rad(23.4392) * math.sin(rad(365 * day / 365.25))

def sundir(day): # sun vector (West, Above, South) (~Povray)
return array([[0],[math.sin(sunh(day))], [math.cos(sunh(day))], [1]])

def Rx(theta): # rotation around East-West axis
return array(
[[1, 0, 0, 0],
[0, math.cos(theta), -math.sin(theta), 0],
[0, math.sin(theta), cos(theta), 0],
[0, 0, 0, 1]])

def Ry(theta): # rotation around the vertical axis
return array(
[[math.cos(theta), 0, math.sin(theta), 0],
[0, 1, 0, 0],
[-math.sin(theta), 0, cos(theta), 0],
[0, 0, 0, 1]])

day = int(sys.argv[1]); latitude = 50.6
print "day=",day," --- latitude = ",latitude
for hour in range(6, 18):
print hour, "h, ",
print zip(*matrixmultiply(Rx(-rad(90-latitude)), matrixmultiply(Ry(h2rad(hour)), sundir(day))))

Qui donne, pour le 81-ième jour (21 mars) (heures solaires) :

day= 81 --- latitude = 50.6
6 h, [(-1.0, 3.8864750971621759e-17, 4.7314722195836354e-17, 1.0)]
7 h, [(-0.96592582628906831, 0.16428034532444383, 0.19999816561124292, 1.0)]
8 h, [(-0.8660254037844386, 0.31736525660113385, 0.38636678674867564, 1.0)]
9 h, [(-0.70710678118654746, 0.44882225011134086, 0.54640514987049038, 1.0)]
10 h, [(-0.49999999999999994, 0.54969274899029774, 0.66920690500583568, 1.0)]
11 h, [(-0.25881904510252074, 0.61310259543578471, 0.74640331548173333, 1.0)]
12 h, [(0.0, 0.63473051320226759, 0.77273357349735106, 1.0)]
13 h, [(0.25881904510252074, 0.61310259543578471, 0.74640331548173333, 1.0)]
14 h, [(0.49999999999999994, 0.54969274899029774, 0.66920690500583568, 1.0)]
15 h, [(0.70710678118654746, 0.44882225011134086, 0.54640514987049038, 1.0)]
16 h, [(0.8660254037844386, 0.31736525660113385, 0.38636678674867564, 1.0)]
17 h, [(0.96592582628906831, 0.16428034532444383, 0.19999816561124292, 1.0)]