code pour animation
gron...@gmail.com
je trouve le concept de la quasiturbine assez génial
et je voulais essayer de me faire ma propre animation
pour illustrer son fonctionnement.
J'ai finalement réussi à écrire le code ci-dessous,
utilisable sous linux avec bc, plotutils et imagemagick.
pour installer ces paquets sous ubuntu :
sudo apt-get install bc plotutils imagemagick
A toute fin utile :
#!/usr/bin/sh
export BC_LINE_LENGTH=5000
bc -ql <<EOF |sh
pi=3.14159265359
a=0.22 /* fonctionne avec a<0.32 , le cas limite valant le
coup d'oeil ! */
define phi(x) {
return c(pi/4*(1+a*c(2*x)))
}
print "cat << FIN_STATOR > /tmp/stator\n"
for(k=0;k<400;k++) {
theta = pi/2*k/100
cos = c(theta)
sin = s(theta)
print phi(theta)*cos, " ", phi(theta)*sin, "\n"
}
print "FIN_STATOR\n"
for(k=0;k<100;k++) {
theta = pi/2*k/100
cos = c(theta)
sin = s(theta)
print "echo "
print phi(theta)*cos, " ", phi(theta)*sin, " "
print -phi(theta+pi/2)*sin, " ", phi(theta+pi/2)*cos,
" "
print -phi(theta+pi)*cos, " ", -phi(theta+pi)*sin, " "
print phi(theta+3*pi/2)*sin, " ", -phi(theta+3*pi/
2)*cos, " "
print phi(theta)*cos, " ", phi(theta)*sin, " "
print "|graph -T pnm -x -1 1 -y -1 1 --blankout 0.0 -
"
print "--reposition 0 0 1 -x -1 1 -y -1 1 /tmp/stator "
print "> img-", 100+k, ".pnm"
print "\n"
}
quit
EOF
rm /tmp/stator
animate img-*.pnm
rm img-*.pnm
gron...@gmail.com
Je me suis donc démené pour écrire un modèle 3D sous Blender.
Par contre pas d'animation 3D pour l'instant.
Quelques copies d'écran :
http://picasaweb.google.fr/grondilu/Quasiturbine
Le script python correspondant (un peu long, mais tout y est) :
import Blender
from Blender import Scene, Mathutils
from Blender.Mathutils import *
from math import *
scn=Scene.getCurrent()
for ob in scn.objects:
if ob.type == 'Mesh':
scn.unlink(ob)
Blender.Redraw()
# PARAMETRES
a=.3 # "exentricit\'e" du stator
b=.2 # rayon des cylindres de contact
N=200 # r\'esolution angulaire
global_angle=0 # dephasage du moteur
# FONCTIONS
def phi(arg): # angle au sommet du losange
return pi/4*(1+a*cos(2*arg))
def phip(arg): # deriv\'ee premi\`ere de phi
return -a*pi/2*sin(2*arg)
def r(arg): # \'equation polaire des centres des cylindres
return (1-b)*cos(phi(arg))
def rp(arg): # deriv\'ee premi\`ere de r
return -(1-b)*phip(arg)*sin(phi(arg))
# STATOR
angles=[k*2*pi/N for k in range(N)]
statorInterieur2D=[]
statorExterieur2D=[]
for theta in angles:
n=sqrt(r(theta)**2+rp(theta)**2)
statorInterieur2D+=[[cos(theta)*r(theta)*(1+b/n)+b*rp(theta)/
n*sin(theta), \
sin(theta)*r(theta)*(1+b/n)-b*rp(theta)/n*cos(theta)]]
statorExterieur2D+=[[cos(theta), sin(theta)]]
verts=[]
for k in range(N):
verts+=[ statorInterieur2D[k]+[0],\
statorInterieur2D[k]+[1],\
statorExterieur2D[k]+[1],\
statorExterieur2D[k]+[0]]
faces=[]
for k in range(0,4*N,4):
faces+=[ [k, (k-4)%(4*N), (k-3)%(4*N), k+1], \
[k, (k-4)%(4*N), (k-1)%(4*N), k+3], \
[k+3, (k-1)%(4*N), (k-2)%(4*N), k+2], \
[k+1, k+2, (k-2)%(4*N), (k-3)%(4*N)] ]
me=Blender.Mesh.New('stator')
me.verts.extend(verts)
me.faces.extend(faces)
stator=Blender.Scene.GetCurrent().objects.new(me, 'stator')
# ROTOR
rotor=[]
alpha=pi/3
for n in range(4):
repere=[]
verts=[]
faces=[]
i=0
za, zb = .025 + n*.25, .225 + n*.25
theta=-pi+alpha/2
for z in 0, 1, za, zb:
verts.append([b/2*cos(theta), b/2*sin(theta),z])
i+=1
theta=-pi+alpha/2
while(theta<pi+alpha/2):
if theta<=pi-alpha/2:
for z in 0, 1, za, zb:
verts.append([b/2*cos(theta), b/2*sin(theta),z])
i+=1
faces.append([i-8,i-6,i-2,i-4])
faces.append([i-5,i-1,i-3,i-7])
faces.append([3,i-1,i-5])
faces.append([2,i-2,i-6])
faces.append([1,i-3,i-7])
faces.append([0,i-4,i-8])
elif theta<pi:
if len(repere)==0:
repere.append(i-4)
theta_backup=theta-2*pi/N
for z in za, zb:
verts.append([b/2*cos(theta), b/2*sin(theta),z])
i+=1
faces.append([3,i-1,i-3])
faces.append([2,i-2,i-4])
faces.append([i-1,i-2,i-4,i-3])
else:
if len(repere)==1:
theta=pi
repere.append(i)
for z in za, zb:
verts.append([b/2*cos(theta), b/2*sin(theta),z])
i+=1
faces.append([3,i-1,i-3])
faces.append([2,i-2,i-4])
theta+=2*pi/N
theta=theta_backup
for z in 0, 1:
verts.append([b*cos(theta),b*sin(theta),z])
i+=1
faces.append([repere[0],repere[0]+1,i-1,i-2])
# repère 2
repere.append(i-2)
for z in 0, 1:
verts.append([b*cos(theta)-b*(1-sin(theta))*sin(theta)/
cos(theta),b,z])
i+=1
faces.append([i-1,i-2,i-4,i-3])
for z in 0, 1:
verts.append([b-1,b,z])
i+=1
faces.append([i-1,i-2,i-4,i-3])
theta=pi/2
while theta < 2*pi-alpha/2-pi/2*(1+a):
for z in 0,1:
verts.append([b-1+b*cos(theta),b*sin(theta),z])
i+=1
faces.append([i-1,i-2,i-4,i-3])
theta+=2*pi/N
# repère 3
repere.append(i-2)
theta-=2*pi/N
for z in 0, 1:
verts.append([b-1+b/2*cos(theta),b/2*sin(theta),z])
i+=1
faces.append([i-1,i-2,i-4,i-3])
while theta > 0:
for z in 0, 1:
verts.append([b-1+b/2*cos(theta), b/2*sin(theta),z])
i+=1
faces.append([i-1,i-2,i-4,i-3])
theta-=2*pi/N
za, zb = .025 + (n+1)%4*.25, .225 + (n+1)%4*.25
theta=0
for z in 0, 1, za, zb:
verts.append([b-1+b/2*cos(theta), b/2*sin(theta),z])
i+=1
theta-=2*pi/N
faces.append([i-4,i-3,i-5,i-6])
# repère 4
repere.append(i-4)
thetamax = max(alpha/2-pi/2*(1-a),-pi*a)
while theta > thetamax:
for z in 0, 1, za, zb:
verts.append([b-1+b/2*cos(theta), b/2*sin(theta),z])
i+=1
theta-=2*pi/N
faces.append([i-4,i-2,i-6,i-8])
faces.append([i-1,i-3,i-7,i-5])
theta+=2*pi/N
for z in 0, 1, za, zb:
verts.append([b-1+b*cos(theta), b*sin(theta),z])
i+=1
faces.append([i-4,i-2,i-6,i-8])
faces.append([i-1,i-3,i-7,i-5])
# repère 5
repere.append(i-4)
for z in 0,1,za, zb:
verts.append([3*b/2-1+b/2*sin(alpha/2)/tan(pi/4*(1-a)),-b/
2*sin(alpha/2),z])
i+=1
for j in range(4):
faces.append([i-8+j,i-12+j,i-4+j])
k=i-12+j
while k>repere[4]+j:
faces.append([k,k-4,i-4+j])
k-=4
faces.append([i-4,i-2,i-6,i-8])
faces.append([i-1,i-3,i-7,i-5])
faces.append([repere[4]+2,repere[4]+3,i-1,i-2])
za, zb = .025 + n*.25, .225 + n*.25
for z in za, zb:
verts.append([-b/2*(1+sin(alpha/2)/tan(pi/4*(1-a))),-b/2*sin(alpha/
2),z])
i+=1
faces.append([i-1,i-2,repere[1],repere[1]+1])
j=repere[1]+2
while (verts[j][1]<verts[j-2][1]):
faces.append([j,j-2,i-2])
faces.append([j+1,j-1,i-1])
j+=2
faces.append([j-2,2,i-2])
faces.append([j-1,3,i-1])
if n<3:
faces.append([i-2,2,0,i-6])
faces.append([i-6,i-4,i-1,i-2])
faces.append([i-4,i-3,3,i-1])
faces.append([3,i-3,i-5,1])
else:
faces.append([i-1,3,1,i-5])
faces.append([i-3,i-2,i-1,i-5])
faces.append([i-3,i-4,2,i-2])
faces.append([i-4,i-6,0,2])
for k in 0, 1:
faces.append([k,repere[0]+k,repere[2]+k, repere[-1]+4+k])
faces.append([repere[4]+k,repere[2]+2+k,repere[2]+k,
repere[-1]+4+k])
k=0
while repere[3]-k>=repere[2]+6:
faces.append([repere[3]-k+1,repere[3]-k-1,repere[3]+k+7,repere[3]+k
+5])
faces.append([repere[3]-k,repere[3]-k-2,repere[3]+k+6,repere[3]+k
+4])
k+=2
k1=repere[3]+k
k2=repere[4]
while k1+4<k2-2:
faces.append([k1+4, k1+2, repere[2]+6])
faces.append([k1+5, k1+3, repere[2]+7])
k1+=2
faces.append([k2, k2-2, repere[2]+2])
faces.append([k2+1, k2-1, repere[2]+3])
k2-=2
for k in 0, 1:
if k1+4<>k2-2:
faces.append([k1+4+k,k2-2+k,repere[2]+6+k,repere[2]+2+k])
else:
faces.append([k1+4+k,k1+2+k,repere[2]+6+k,repere[2]+2+k])
faces.append([k2+k,k2-2+k,repere[2]+6+k,repere[2]+2+k])
theta=n*pi/2+global_angle
R = RotationMatrix(-360*(theta-phi(theta))/(2*pi),4,'z')
T = TranslationMatrix(Vector([r(theta)*cos(theta),
r(theta)*sin(theta), 0]))
Verts=[]
for v in verts:
Verts.append(Vector(R*Vector(v)*T))
me=Blender.Mesh.New('rotor_' +str(n))
me.verts.extend(Verts)
me.faces.extend(faces)
rotor.append(Blender.Scene.GetCurrent().objects.new(me, 'rotor_'+
str(n)))
za, zb = .025 + n*.25, .225 + n*.25
i=0
verts=[]
faces=[]
theta=-pi+pi/4*(1-a)
for z in za, zb:
verts.append([b/2*cos(theta), b/2*sin(theta),z])
i+=1
theta=-pi+pi/4*(1-a)
while(theta<pi-pi/4*(1-a)):
for z in za, zb:
verts.append([b/2*cos(theta), b/2*sin(theta),z])
i+=1
faces.append([i-3,i-1,i-2,i-4])
faces.append([1,i-1,i-3])
faces.append([0,i-2,i-4])
theta+=2*pi/N
theta=pi-pi/4*(1-a)
for z in za, zb:
verts.append([b/2*cos(theta), b/2*sin(theta),z])
i+=1
faces.append([i-3,i-1,i-2,i-4])
faces.append([1,i-1,i-3])
faces.append([0,i-2,i-4])
for z in za, zb:
verts.append([b-2*(1-b)*cos(pi/4*(1+a)),b/2*sin(theta),z])
verts.append([b-2*(1-b)*cos(pi/4*(1+a)),-b/2*sin(theta),z])
i+=2
faces.append([i-1,i-2,i-4,i-3])
faces.append([i-1,i-3,0,1])
faces.append([i-1,i-2,i-5,1])
faces.append([i-5,i-6,i-4,i-2])
faces.append([i-4,i-3,0,i-6])
theta=n*pi/2+global_angle
R = RotationMatrix(-360*theta/(2*pi),4,'z')
T = TranslationMatrix(Vector([r(theta)*cos(theta),
r(theta)*sin(theta), 0]))
Verts=[]
for v in verts:
Verts.append(Vector(R*Vector(v)*T))
me=Blender.Mesh.New('axe_' +str(n))
me.verts.extend(Verts)
me.faces.extend(faces)
rotor.append(Blender.Scene.GetCurrent().objects.new(me, 'axe_'+
str(n)))
i=0
verts=[]
faces=[]
L=((1-b)*cos(pi/4*(1+a))-b)
l=b/2*sin(pi/4*(1-a))
e=.0125
for x in -L, L:
for y in -l-e, -l, l, l+e:
for z in za-e, za, zb, zb+e:
verts.append([x,y,z])
i+=1
faces.append([0,16,19,3])
faces.append([3,19,31,15])
faces.append([0,12,28,16])
faces.append([12,28,31,15])
faces.append([5,21,22,6])
faces.append([9,25,26,10])
faces.append([5,9,25,21])
faces.append([22,26,10,6])
for k in 0, 16:
faces.append([0+k,3+k,6+k,5+k])
faces.append([0+k,12+k,9+k,5+k])
faces.append([12+k,15+k,10+k,9+k])
faces.append([15+k,3+k,6+k,10+k])
Verts=[]
for v in verts:
Verts.append(Vector(R*Vector(v)))
me=Blender.Mesh.New('centre_' +str(n))
me.verts.extend(Verts)
me.faces.extend(faces)
rotor.append(Blender.Scene.GetCurrent().objects.new(me, 'centre_'+
str(n)))
Blender.Redraw()
Gaetan
ensuite le transférer sous Lightwave !
Penses-tu le terminer?
Merci
Gaetan