import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# Parameters of the car suspension system
m = 1000  # mass of the car (kg)
c = 1500  # damping coefficient (Ns/m)
k = 20000  # spring constant (N/m)

# Impulse force due to the pothole
def impulse_force(t, duration=0.05, magnitude=5000):
    if 0 <= t <= duration:
        return magnitude / duration
    else:
        return 0

# Define the system of ODEs
def car_suspension_system(t, y):
    x, x_dot = y
    force = impulse_force(t)
    x_ddot = (force - c * x_dot - k * x) / m
    return [x_dot, x_ddot]

# Initial conditions
x0 = 0.0  # initial displacement
x_dot0 = 0.0  # initial velocity
initial_conditions = [x0, x_dot0]

# Time span for the simulation
t_span = (0, 2)  # simulate for 2 seconds
t_eval = np.linspace(t_span[0], t_span[1], 1000)

# Solve the ODE
solution = solve_ivp(car_suspension_system, t_span, initial_conditions, t_eval=t_eval)

# Extract the results
t = solution.t
x = solution.y[0]

# Plot the results
plt.figure(figsize=(10, 6))
plt.plot(t, x, label='Displacement (x)')
plt.title('Vibration of the Car Over a Pothole')
plt.xlabel('Time (s)')
plt.ylabel('Displacement (m)')
plt.legend()
plt.grid(True)
plt.show()