Редактировать: Извините, я не думал о написании тестов.Я сделаю это и посмотрю, не смогу ли я выяснить, что я сделал неправильно.Спасибо человеку, который предложил мне написать тесты!
Я пытаюсь написать компьютерную симуляцию на Python, которая симулирует электрическую силу и то, как атомы взаимодействуют с ней.Для тех, кто не знает, по существу, вещи с противоположными (положительными и отрицательными) зарядами притягиваются, и подобно зарядам отталкиваются, и величина силы падает как 1 / (квадрат расстояния).Я пытаюсь поместить отрицательно заряженную частицу (ион кислорода) и положительно заряженную частицу (ион водорода) в систему координат и ожидаю, что они будут притягиваться и сближаться, но вместо этого они отталкиваются!Естественно, я думал, что это просто опечатка, поэтому я добавил отрицательный знак к моей модели электрической силы, но они все еще отражают!Я понятия не имею, что происходит, и надеялся, что некоторые из вас, ребята, могут иметь представление о том, что здесь происходит не так.Ниже приведен код.Я включил все, чтобы вы могли просто запустить его из собственного терминала, чтобы увидеть для себя, что происходит.
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np
import math
Å = 10 ** (-10)
u = 1.660539040 * (10 ** (-27))
k = 8.987551787368 * (10 ** 9)
e = 1.6021766208 * (10 ** (-19))
hydrogen1 = dict(
name='Hydrogen 1',
charge=1 * e,
mass=1.00793 * u,
position=[0.5 * Å, 0.5 * Å, 0.5 * Å],
velocity=[0, 0, 0],
acceleration=[0, 0, 0],
force=[0, 0, 0]
)
oxygen1 = dict(
name='Oxygen 1',
charge=-2 * e,
mass=15.9994 * u,
position=[0, 0, 0],
velocity=[0, 0, 0],
acceleration=[0, 0, 0],
force=[0, 0, 0]
)
atoms = [hydrogen1, oxygen1]
def magnitude(vector):
magnitude = 0
for coordinate in vector:
magnitude += (coordinate ** 2)
return math.sqrt(magnitude)
def scale_vector(vector, scalefactor):
scaled_vector = vector
i = 0
while i < len(vector):
vector[i] *= scalefactor
i += 1
return scaled_vector
def sum_vectors(vectors):
resultant_vector = [0, 0, 0]
for vector in vectors:
i = 0
while i < len(vector):
resultant_vector[i] += vector[i]
i += 1
return resultant_vector
def distance_vector(point1, point2):
if type(point1) is list and type(point2) is list:
pos1 = point1
pos2 = point2
elif type(point1) is dict and type(point2) is dict:
pos1 = point1['position']
pos2 = point2['position']
vector = []
i = 0
while i < len(pos1):
vector.append(pos2[i] - pos1[i])
i += 1
return vector
def distance(point1, point2):
return magnitude(distance_vector(point1, point2))
def direction_vector(point1, point2):
vector = distance_vector(point1, point2)
length = magnitude(vector)
return scale_vector(vector, 1 / length)
def eletric_force(obj1, obj2):
length = k * obj1['charge'] * \
obj2['charge'] / ((distance(obj1, obj2)) ** 2)
force_vector = scale_vector(direction_vector(obj1, obj2), length)
return force_vector
def force_to_acceleration(force, mass):
scalefactor = 1 / (mass)
return scale_vector(force, scalefactor)
time = 10
t = 0
period = 1 / 1000
while t < time:
i = 0
while i < len(atoms):
atom = atoms[i]
position = atom['position']
velocity = atom['velocity']
acceleration = atom['acceleration']
# Moving the atom
atom['position'] = sum_vectors(
[position, scale_vector(velocity, period)])
# Accelerating the atom using its current acceleration vector
atom['velocity'] = sum_vectors([
velocity, scale_vector(acceleration, period)])
# Calculating the net force on the atom
force = [0, 0, 0]
j = 0
while j < len(atoms):
if j != i:
force = sum_vectors([force, eletric_force(atoms[i], atoms[j])])
j += 1
# Updating the force and acceleration on the atom
atoms[i]['force'] = [force[0], force[1], force[2]]
atom['acceleration'] = force_to_acceleration(
[force[0], force[1], force[2]], atom['mass'])
i += 1
t += period
np.random.seed(19680801)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
for atom in atoms:
name = atom['name']
position = atom['position']
X = position[0]
Y = position[1]
Z = position[2]
print(
f'Position of {name}: [{X}, {Y}, {Z}]')
color = 'green'
if 'Oxygen' in atom['name']:
color = 'red'
ax.scatter(X, Y, Z, color=color)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
plt.show()
Ниже приведены тесты:
from functions import *
from constants import *
from atoms import hydrogen1, hydrogen2, oxygen1
import math
def test_magnitude():
assert magnitude({1, 3, -5}) == math.sqrt(35)
def test_sum_vectors():
assert sum_vectors([[1, 2, 3], [0, -4, 8]]) == [1, -2, 11]
def test_scale_vector():
assert scale_vector([1, 4, -3], -2) == [-2, -8, 6]
def test_distance_vector():
assert distance_vector([1, 4, 3], [0, 1, 1]) == [-1, -3, -2]
assert distance_vector(hydrogen2, hydrogen1) == [Å, Å, Å]
def test_distance():
assert distance([1, 2, 3], [3, 2, 1]) == math.sqrt(8)
assert distance(hydrogen1, oxygen1) == math.sqrt(0.75) * Å
assert distance(hydrogen1, hydrogen2) == Å * math.sqrt(3)
def test_direction_vector():
assert direction_vector([1, 1, 1], [7, 5, -3]) == [6 /
math.sqrt(68), 4 / math.sqrt(68), -4 / math.sqrt(68)]
m = 1 / math.sqrt(3)
for component in direction_vector(hydrogen2, hydrogen1):
assert abs(component - m) < 10 ** (-12)
def test_electric_force():
m = 4.439972744 * 10 ** (-9)
for component in electric_force(hydrogen1, hydrogen2):
assert abs(component - m) < 10 ** (-12)
def test_force_to_acceleration():
assert force_to_acceleration(
[4, 3, -1], 5.43) == [4 / 5.43, 3 / 5.43, -1 / 5.43]