import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.animation as animation

def simulate_particle_in_magnetic_field(initial_position, initial_velocity, magnetic_field, time_step, num_steps):
    positions = [initial_position]
    velocities = [initial_velocity]

    for _ in range(num_steps):
        velocity = velocities[-1]
        position = positions[-1]

        # Calculate the acceleration using the Lorentz force equation
        acceleration = np.cross(velocity, magnetic_field)

        # Update the velocity and position using the Euler method
        new_velocity = velocity + acceleration * time_step
        new_position = position + velocity * time_step

        velocities.append(new_velocity)
        positions.append(new_position)

    return np.array(positions), np.array(velocities)

# Constants
initial_position = np.array([0, 0, 0])  # Initial position (x, y, z)
initial_velocity = np.array([1, 0, 0])  # Initial velocity (vx, vy, vz)
magnetic_field = np.array([0, 0, 1])  # Magnetic field (Bx, By, Bz)
time_step = 0.01  # Time step size
num_steps = 100  # Number of simulation steps

# Simulate the particle
positions, velocities = simulate_particle_in_magnetic_field(initial_position, initial_velocity, magnetic_field, time_step, num_steps)

# Create the figure and axes
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

# Set up the initial plot
particle, = ax.plot([], [], [], 'ro')

# Set the plot limits
ax.set_xlim3d(-2, 2)
ax.set_ylim3d(-2, 2)
ax.set_zlim3d(-2, 2)

# Animation update function
def update(frame):
    x, y, z = positions[frame]

    # Update the particle's position
    particle.set_data([x], [y])
    particle.set_3d_properties([z])

    return particle,

# Create the animation
ani = animation.FuncAnimation(fig, update, frames=len(positions), interval=100, blit=True)

# Display the animation
plt.show()

aW1wb3J0IG51bXB5IGFzIG5wCmltcG9ydCBtYXRwbG90bGliLnB5cGxvdCBhcyBwbHQKZnJvbSBtcGxfdG9vbGtpdHMubXBsb3QzZCBpbXBvcnQgQXhlczNECmltcG9ydCBtYXRwbG90bGliLmFuaW1hdGlvbiBhcyBhbmltYXRpb24KCmRlZiBzaW11bGF0ZV9wYXJ0aWNsZV9pbl9tYWduZXRpY19maWVsZChpbml0aWFsX3Bvc2l0aW9uLCBpbml0aWFsX3ZlbG9jaXR5LCBtYWduZXRpY19maWVsZCwgdGltZV9zdGVwLCBudW1fc3RlcHMpOgogICAgcG9zaXRpb25zID0gW2luaXRpYWxfcG9zaXRpb25dCiAgICB2ZWxvY2l0aWVzID0gW2luaXRpYWxfdmVsb2NpdHldCgogICAgZm9yIF8gaW4gcmFuZ2UobnVtX3N0ZXBzKToKICAgICAgICB2ZWxvY2l0eSA9IHZlbG9jaXRpZXNbLTFdCiAgICAgICAgcG9zaXRpb24gPSBwb3NpdGlvbnNbLTFdCgogICAgICAgICMgQ2FsY3VsYXRlIHRoZSBhY2NlbGVyYXRpb24gdXNpbmcgdGhlIExvcmVudHogZm9yY2UgZXF1YXRpb24KICAgICAgICBhY2NlbGVyYXRpb24gPSBucC5jcm9zcyh2ZWxvY2l0eSwgbWFnbmV0aWNfZmllbGQpCgogICAgICAgICMgVXBkYXRlIHRoZSB2ZWxvY2l0eSBhbmQgcG9zaXRpb24gdXNpbmcgdGhlIEV1bGVyIG1ldGhvZAogICAgICAgIG5ld192ZWxvY2l0eSA9IHZlbG9jaXR5ICsgYWNjZWxlcmF0aW9uICogdGltZV9zdGVwCiAgICAgICAgbmV3X3Bvc2l0aW9uID0gcG9zaXRpb24gKyB2ZWxvY2l0eSAqIHRpbWVfc3RlcAoKICAgICAgICB2ZWxvY2l0aWVzLmFwcGVuZChuZXdfdmVsb2NpdHkpCiAgICAgICAgcG9zaXRpb25zLmFwcGVuZChuZXdfcG9zaXRpb24pCgogICAgcmV0dXJuIG5wLmFycmF5KHBvc2l0aW9ucyksIG5wLmFycmF5KHZlbG9jaXRpZXMpCgojIENvbnN0YW50cwppbml0aWFsX3Bvc2l0aW9uID0gbnAuYXJyYXkoWzAsIDAsIDBdKSAgIyBJbml0aWFsIHBvc2l0aW9uICh4LCB5LCB6KQppbml0aWFsX3ZlbG9jaXR5ID0gbnAuYXJyYXkoWzEsIDAsIDBdKSAgIyBJbml0aWFsIHZlbG9jaXR5ICh2eCwgdnksIHZ6KQptYWduZXRpY19maWVsZCA9IG5wLmFycmF5KFswLCAwLCAxXSkgICMgTWFnbmV0aWMgZmllbGQgKEJ4LCBCeSwgQnopCnRpbWVfc3RlcCA9IDAuMDEgICMgVGltZSBzdGVwIHNpemUKbnVtX3N0ZXBzID0gMTAwICAjIE51bWJlciBvZiBzaW11bGF0aW9uIHN0ZXBzCgojIFNpbXVsYXRlIHRoZSBwYXJ0aWNsZQpwb3NpdGlvbnMsIHZlbG9jaXRpZXMgPSBzaW11bGF0ZV9wYXJ0aWNsZV9pbl9tYWduZXRpY19maWVsZChpbml0aWFsX3Bvc2l0aW9uLCBpbml0aWFsX3ZlbG9jaXR5LCBtYWduZXRpY19maWVsZCwgdGltZV9zdGVwLCBudW1fc3RlcHMpCgojIENyZWF0ZSB0aGUgZmlndXJlIGFuZCBheGVzCmZpZyA9IHBsdC5maWd1cmUoKQpheCA9IGZpZy5hZGRfc3VicGxvdCgxMTEsIHByb2plY3Rpb249JzNkJykKCiMgU2V0IHVwIHRoZSBpbml0aWFsIHBsb3QKcGFydGljbGUsID0gYXgucGxvdChbXSwgW10sIFtdLCAncm8nKQoKIyBTZXQgdGhlIHBsb3QgbGltaXRzCmF4LnNldF94bGltM2QoLTIsIDIpCmF4LnNldF95bGltM2QoLTIsIDIpCmF4LnNldF96bGltM2QoLTIsIDIpCgojIEFuaW1hdGlvbiB1cGRhdGUgZnVuY3Rpb24KZGVmIHVwZGF0ZShmcmFtZSk6CiAgICB4LCB5LCB6ID0gcG9zaXRpb25zW2ZyYW1lXQoKICAgICMgVXBkYXRlIHRoZSBwYXJ0aWNsZSdzIHBvc2l0aW9uCiAgICBwYXJ0aWNsZS5zZXRfZGF0YShbeF0sIFt5XSkKICAgIHBhcnRpY2xlLnNldF8zZF9wcm9wZXJ0aWVzKFt6XSkKCiAgICByZXR1cm4gcGFydGljbGUsCgojIENyZWF0ZSB0aGUgYW5pbWF0aW9uCmFuaSA9IGFuaW1hdGlvbi5GdW5jQW5pbWF0aW9uKGZpZywgdXBkYXRlLCBmcmFtZXM9bGVuKHBvc2l0aW9ucyksIGludGVydmFsPTEwMCwgYmxpdD1UcnVlKQoKIyBEaXNwbGF5IHRoZSBhbmltYXRpb24KcGx0LnNob3coKQoK