Sunday, September 9, 2018

Some science


import numpy as np
import matplotlib.pyplot as plt

# --- Physical Constants ---
g = 9.81          # Gravity (m/s^2)
dt = 0.01         # Time step (seconds)
v0 = 50.0         # Initial velocity (m/s)
angle = 45.0      # Launch angle (degrees)
k = 0.01          # Air resistance coefficient (drag)

# Convert angle to radians
theta = np.radians(angle)

# Initial velocity components
vx = v0 * np.cos(theta)
vy = v0 * np.sin(theta)

# Lists to store the trajectory
x_points = [0.0]
y_points = [0.0]

# --- Numerical Integration (Euler's Method) ---
while y_points[-1] >= 0:
    # Current velocity magnitude
    v = np.sqrt(vx**2 + vy**2)
    
    # Calculate accelerations (F = ma, so a = F/m)
    # Drag force acts opposite to velocity: -k * v * vx
    ax = -k * v * vx
    ay = -g - (k * v * vy)
    
    # Update velocities
    vx += ax * dt
    vy += ay * dt
    
    # Update positions
    new_x = x_points[-1] + vx * dt
    new_y = y_points[-1] + vy * dt
    
    x_points.append(new_x)
    y_points.append(new_y)

# --- Visualization ---
plt.figure(figsize=(10, 5))
plt.plot(x_points, y_points, label=f"Trajectory (Drag k={k})", color='blue', lw=2)
plt.title("Projectile Motion with Air Resistance")
plt.xlabel("Distance (m)")
plt.ylabel("Height (m)")
plt.axhline(0, color='black', lw=1) # Ground line
plt.grid(True, linestyle='--', alpha=0.6)
plt.legend()
plt.show()

print(f"Total Distance Covered: {x_points[-1]:.2f} meters")

No comments: