Программа моделирования гравитации
Переделанная из GAS.PY Язык PYTHON библиотека Visual
from visual import * from visual.graph import * from random import random # A model of an ideal gas with hard-sphere collisions # Program uses numpy arrays for high speed computations # Bruce Sherwood win=800 Natoms =2400 # change this to have more or fewer atoms # Typical values L = 0.8# container is a cube L on a side gray = (0.7,0.7,0.7) # color of edges of container Matom = 4E-3/6E23 # helium mass Ratom = 0.01 # wildly exaggerated size of helium atom #k = 1.4E-23 # Boltzmann constant k = 5.0E-20 # Boltzmann constant T = 0.01 # around room temperature kd=1.0 #dt = 0.0 dt = 1E-5 scene = display(title="Gas", width=win, height=win, x=0, y=0, center=(L/2.,L/2.,L/2.)) deltav = 100. # binning for v histogram vdist = gdisplay(x=0, y=win, ymax = Natoms*deltav/1000., width=win, height=0.6*win, xtitle='v', ytitle='dN') theory = gcurve(color=color.cyan) dv = 10. for v in arange(0.,3001.+dv,dv): # theoretical prediction theory.plot(pos=(v, (deltav/dv)*Natoms*4.*pi*((Matom/(2.*pi*k*T))**1.5) *exp((-0.5*Matom*v**2)/(k*T))*v**2*dv)) observation = ghistogram(bins=arange(0.,3000.,deltav), accumulate=1, average=1, color=color.red) xaxis = curve(pos=[(0,0,0), (L,0,0)], color=gray) yaxis = curve(pos=[(0,0,0), (0,L,0)], color=gray) zaxis = curve(pos=[(0,0,0), (0,0,L)], color=gray) xaxis2 = curve(pos=[(L,L,L), (0,L,L), (0,0,L), (L,0,L)], color=gray) yaxis2 = curve(pos=[(L,L,L), (L,0,L), (L,0,0), (L,L,0)], color=gray) zaxis2 = curve(pos=[(L,L,L), (L,L,0), (0,L,0), (0,L,L)], color=gray) Atoms = [] colors = [color.red, color.green, color.blue, color.yellow, color.cyan, color.magenta] poslist = [] xyzclist = [] xlist = [] ylist = [] zlist = [] plist = [] mlist = [] rlist = [] rrlist = [] cclist = [] d1=12 xd1=0.34 yd1=0.34 zd1=0.34 d2=6 xd2=0.48 yd2=0.48 zd2=0.48 dd1=Ratom*d1*kd dd2=Ratom*d2*kd na=0 nna=0 while na < Natoms-2: rrlist.append(Ratom) cclist.append(colors[0]) Lmin = 1.1*Ratom Lmax = L-Lmin lmm= Lmin+(Lmax-Lmin) xr=lmm*random() yr=lmm*random() zr=lmm*random() if sqrt((xd1-xr)*(xd1-xr)+(yd1-yr)*(yd1-yr)+(zd1-zr)*(zd1-zr))<(dd1/2) or sqrt((xd2-xr)*(xd2-xr)+(yd2-yr)*(yd2-yr)+(zd2-zr)*(zd2-zr))<(dd2/2): nna=nna+1 else: na=na+1 xyzclist.append((xr,yr,zr,colors[0],Ratom)) print na, nna,Natoms xyzclist.append((xd1,yd1,zd1,colors[1],Ratom*d1)) xyzclist.append((xd2,yd2,zd2,colors[2],Ratom*d2)) for i in range(Natoms): x,y,z,vv,r =xyzclist[i] Atoms = Atoms+[sphere(pos=(x,y,z), radius=r, color=vv)] mass = Matom*r**3/Ratom**3 pavg = sqrt(2.*mass*1.5*k*T) # average kinetic energy p**2/(2mass) = (3/2)kT if i Исправил код потом почистил пространство на 10% перед большими шарами это скриншот
это график
Тенденция однако