В настоящее время я работаю над простой симуляцией n-тела, которая, по сути, является школьным проектом, но больше просто для развлечения. Я решил использовать субатоми c частиц вместо астрономических тел. Я основывал свой код на другом коде, который нашел здесь http://www.cyber-omelette.com/2016/11/python-n-body-orbital-simulation.html. Итак, вот моя проблема. Когда я добавляю функцию «Force_électrique_externe», я не получаю никаких результатов. подозреваю, что результаты не записываются в "self.Fe". Вот код:
#importer les modules
import numpy as np
import scipy as sci
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import scipy.constants as scic
import random
#variables qui sont utiles à controler
#définir les contantes
k = 1 / (4 * np.pi * scic.epsilon_0) # Constante de Coulomb N*m^2/C^2
#Définir les classes
class point:
def __init__(self, x,y,z):
self.x = x
self.y = y
self.z = z
#Particules
class part:
def __init__(self, position, masse, vitesse, charge, nom = ""):
self.position = position
self.masse = masse
self.vitesse = vitesse
self.charge = charge
self.nom = nom
self.Fe = point(0, 0, 0)
# assigner des variables
#Est- ce que je peux utiliser des symboles par exemple q pour la charge?
#Champ électrique généré par l'intéraction entre les particules
def Force_électrique_particules(parts, part_index):
part_A = parts[part_index]
part_A.Fe.x = 0
part_A.Fe.y = 0
part_A.Fe.z = 0
for index, part_B in enumerate(parts):
if index != part_index:
rcarré = (part_A.position.x - part_B.position.x)**2 + (part_A.position.y - part_B.position.y)**2 + (part_A.position.z - part_B.position.z)**2
r = np.sqrt(rcarré)
E_nv = k * part_A.charge * part_B.charge / r**3
part_A.Fe.x += E_nv * (part_B.position.x - part_A.position.x)
part_A.Fe.y += E_nv * (part_B.position.y - part_A.position.y)
part_A.Fe.z += E_nv * (part_B.position.z - part_A.position.z)
#def Force_électrique_externe(part_A):
#E_ext = point(0, 0, 0)
#part_A.Fe.x += part_A.charge * E_ext.x
#part_A.Fe.y += part_A.charge * E_ext.y
#part_A.Fe.z += part_A.charge * E_ext.z
#def Force_électrique(parts, part_index):
#Fe = point(0, 0, 0)
#Fe_part = Force_électrique_externe(parts, part_index)
#Fe_ext = Force_électrique_externe(parts, part_index)
#part_A.Fe.x += Fe_ext.x + Fe_part.x
#part_A.Fe.y += Fe_ext.y + Fe_part.y
#part_A.Fe.z += Fe_ext.z + Fe_part.z
#return Fe
#accélération
def accélération_part_A(parts, part_index):
accélération = point(0, 0, 0)
for part_index, part_A in enumerate(parts):
#appeler la fonction fe qui appele les deux autres fonctions pour les forces électriques
Force_électrique_particules(parts, part_index)
#Force_électrique_externe(part_A)
accélération.x += part_A.Fe.x / part_A.masse
accélération.y += part_A.Fe.y / part_A.masse
accélération.z += part_A.Fe.z / part_A.masse
return accélération
#vitesse
def vitesse(parts, time_step = 0.001):
for part_index, part_A in enumerate(parts):
accélération = accélération_part_A(parts, part_index)
part_A.vitesse.x += accélération.x * time_step
part_A.vitesse.y += accélération.y * time_step
part_A.vitesse.z += accélération.z * time_step
#position
def var_position(parts, time_step = 0.001):
for part_A in parts:
part_A.position.x += part_A.vitesse.x * time_step
part_A.position.y += part_A.vitesse.y * time_step
part_A.position.z += part_A.vitesse.z * time_step
def champ_elec_step(parts, time_step = 0.001):
vitesse(parts, time_step = time_step)
var_position(parts, time_step = time_step)
def plot_output(parts, outfile = None):
fig = plt.figure()
colours = ['r', 'b', 'g', 'y', 'm', 'c']
ax = fig.add_subplot(1,1,1, projection = '3d')
max_range = 0
for part_et in parts:
max_dim = max(max(part_et["x"]), max(part_et["y"]), max(part_et["z"]))
if max_dim > max_range:
max_range = max_dim
ax.plot(part_et["x"], part_et["y"], part_et["z"], c = random.choice(colours), label = part_et["nom"])
ax.set_xlim([-max_range, max_range])
ax.set_ylim([-max_range, max_range])
ax.set_zlim([-max_range, max_range])
if outfile:
plt.savefig(outfile)
plt.legend(loc = 'upper left')
else:
plt.show()
plt.legend(loc = 'upper left')
def faire_simulation(parts, noms = None, time_step = 0.00001, number_of_steps = 1000, report_freq = 10):
#créé l'endroit pour garder les résultats sauvegarger
hist_position_part = []
for part_et in parts:
hist_position_part.append({"x":[], "y":[], "z":[], "nom":part_et.nom})
for i in range(1, number_of_steps):
champ_elec_step(parts, time_step = 0.001)
if i % report_freq == 0:
for index, position_part in enumerate(hist_position_part):
position_part["x"].append(parts[index].position.x)
position_part["y"].append(parts[index].position.y)
position_part["z"].append(parts[index].position.z)
return hist_position_part
#particules
él1 = {"position": point(1, 2, 1), "masse" : scic.m_e, "vitesse" : point(0, 0, 0), "charge" : -scic.e}
él2 = {"position" : point(2, 2, 3), "masse" : scic.m_e, "vitesse" : point(0,0,0), "charge": -scic.e}
prot1 = {"position" : point(1.5, 4, 1), "masse" : scic.m_p, "vitesse" : point(50,0,0), "charge" : scic.e}
prot2 = {"position" : point(1, 4, 3), "masse" : scic.m_p, "vitesse" : point(0, 0, 0), "charge" : scic.e}
if __name__ == "__main__":
parts = [
part(position = él1["position"], masse = él1["masse"], vitesse = él1["vitesse"], charge = él1["charge"], nom = "él1"),
part(position = prot1["position"], masse = prot1["masse"], vitesse = prot1["vitesse"], charge = prot1["charge"], nom = "prot1"),
part(position = prot2["position"], masse = prot2["masse"], vitesse = prot2["vitesse"], charge = prot2["charge"], nom = "prot2"),
part(position = él2["position"], masse = él2["masse"], vitesse = él2["vitesse"], charge = él2["charge"], nom = "él2"),
]
motions = faire_simulation(parts, time_step = 0.001, number_of_steps = 80000, report_freq = 100)
plot_output(motions, outfile = 'déplacement.png')