Et le petit script qui va avec.
#!/usr/bin/env ruby
# Constantes universelles
$G = 6.6742 * (10**(-11)) # N.m2.kg-2
# Constantes terrestres
$earth_r = 6_371_000 # Rayon en metres
$earth_µ = 5510 # Masse volumique en Kg/m3
$earth_m = 5.972*(10**24) # Masse en Kg
# Voxel pour intégration
$ddim = 20000 # Taille du voxel en m
$dx = $ddim
$dy = $ddim
$dz = $ddim
$hdx = $dx/2 # Pour centrage
$hdy = $dy/2 # Idem
$hdz = $dz/2 # Idem
$dvol = $dx * $dy * $dz
$mcube = $dvol * $earth_µ # Masse d'un micro cube
$mg = $mcube * $G
def voxel_attraction(xv, yv, zv, x, y, z)
distx = (x-xv)
disty = (y-yv)
distz = (z-zv)
dsq = distx**2 + disty**2 + distz**2
dst = Math:qrt(dsq)
vforce = $mg/dsq # m1.m2.G/r² (on prend une masse de m2=1 pour avoir la valeur du champ)
vforceuni = vforce/dst # Pour décomposer sur x,y,z
[vforceuni * distx, vforceuni * disty, vforceuni * distz] # Projection de la force
end
def voxel_loop(a,b,c,x,y,z)
xf = yf = zf = 0
# On approxime le début de la boucle au voxel le plus proche
xs = (a/$ddim)*$ddim
ys = (b/$ddim)*$ddim
zs = (c/$ddim)*$ddim
vokc=0
# ON intègre sur x sur [-a;a], sur y sur [-b;b], sur z de [-c à c]
xi = -xs # Arrondi de la borne de départ au voxel le plux proche
while xi < xs
# puts xi
yi = -ys # Idem
while yi < ys
zi = -zs # Idem
while zi < zs
xv = xi+$hdx # Centre du voxel
yv = yi+$hdy # Centre du voxel
zv = zi+$hdz # Centre du voxel
if yield(xv,yv,zv)
f = voxel_attraction(xv,yv,zv,x,y,z)
xf += f[0]
yf += f[1]
zf += f[2]
vokc +=1
end
zi += $dz
end
yi += $dy
end
xi += $dx
end
{ x: xf, y: yf, z: zf, vox: vokc }
end
def pres(res)
puts "Pesanteur: (%.03f, %.03f, %0.03f). Voxels utilisés: #{res[:vox]}" % [res[:x],res[:y],res[:z]]
end
# Champ de gravité en x, y, z
def flat_earth_gravity(a, b, c, x, y, z)
voxel_loop(a,b,c,x,y,z) {|xv,yv,zv|
true # Toujours garder le voxel
}
end
def round_earth_gravity(r, x, y, z)
r2 = r * r
voxel_loop(r,r,r,x,y,z) {|xv, yv, zv|
xv*xv + yv*yv + zv*zv < r2 # Garder le voxel si dans la sphere
}
end
def cylindrical_earth_gravity(r, c, x, y, z)
r2 = r * r
voxel_loop(r,r,c,x,y,z) {|xv, yv, zv|
xv*xv + yv*yv < r2
}
end
# r : Rayon du cercle de la calotte (pas de la sphère !!)
# h : hauteur de la calotte
def spherical_cap_gravity(r, c, x, y ,z)
h = c*2
# Calcul du rayon de la sphere
rsphere = (h*h + r*r)/(2*h)
r2 = rsphere * rsphere
# On place l'origine de la bounding box d'iteration à h/2, il faut faire un peu gaffe
# zv se balade entre -h/2 et h/2
voxel_loop(r,r,c,x,y,z) {|xv, yv, zv|
zs = -rsphere + zv + c # -rsphere + h à la surface de la terre, -rsphere au pole
xv*xv + yv*yv + zs*zs < r2 # Garder le voxel si dans la sphere
}
end
def bench_spherical_earth
puts "Terre ronde"
puts "==========="
pres round_earth_gravity($earth_r,0,0,$earth_r)
puts ""
end
def bench_flat_earth
puts "Terre plate, carrée"
puts "==================="
l = 22_584_672 #
ll = 22_584_672 #
h = 2_123_654
# On prend les demi dimensions
a = l/2
b = ll/2
c = h/2
puts "1) Milieu, en surface"
pres flat_earth_gravity(a,b,c,0,0,c) # Centre, surface
puts "2) Extrémité Est, en surface"
pres flat_earth_gravity(a,b,c,a,0,c) # Extrémité Est, surface
puts ""
end
def bench_cylindrical_earth
puts "Terre plate, circulaire"
puts "======================="
r = 12_742_036
h = 2_123_654
c = h/2
puts "1) Milieu, en surface"
pres cylindrical_earth_gravity(r,c,0,0,c) # Centre, surface
puts "2) Extrémité Est, en surface"
pres cylindrical_earth_gravity(r,c,r,0,c) # Extrémité Est, surface
puts ""
end
def bench_spherical_cap_earth
puts "Terre calotte sphérique"
puts "======================="
r = 12_742_036
h = 4_105_286
c = h/2
puts "1) Milieu, en surface"
pres spherical_cap_gravity(r,c,0,0,c) # Centre, surface
puts "2) Extrémité Est, en surface"
pres spherical_cap_gravity(r,c,r,0,c) # Extrémité Est, surface
puts ""
end
bench_spherical_earth
#bench_flat_earth
#bench_cylindrical_earth
bench_spherical_cap_earth