Projectile Motion

import numpy as np
import matplotlib.pyplot as plt

# Projectile Motion Simulation

# ==========================================
# PARAMETERS
# ==========================================

v0 = 25.0        # Initial velocity (m/s)
angle = 45.0     # Launch angle (degrees)
h0 = 0.0         # Initial height (m)
g = 9.81         # Gravitational acceleration (m/s²)

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

# Velocity components
vx = v0 * np.cos(angle_rad)  # Horizontal component
vy = v0 * np.sin(angle_rad)  # Vertical component

print(f"Initial velocity: {v0} m/s at {angle}°")
print(f"Horizontal velocity: {vx:.2f} m/s")
print(f"Vertical velocity: {vy:.2f} m/s")

# ==========================================
# KINEMATIC EQUATIONS
# ==========================================

# Position as function of time:
# x(t) = vx * t
# y(t) = h0 + vy * t - 0.5 * g * t²

# Velocity as function of time:
# vx(t) = vx (constant)
# vy(t) = vy - g * t

# ==========================================
# CALCULATE KEY VALUES
# ==========================================

# Time of flight (when y = 0)
# Solve: 0 = h0 + vy*t - 0.5*g*t²
# Using quadratic formula: t = (vy ± sqrt(vy² + 2*g*h0)) / g

discriminant = vy**2 + 2*g*h0
if discriminant >= 0:
    t_flight = (vy + np.sqrt(discriminant)) / g
else:
    t_flight = 0
    print("Warning: Invalid parameters (negative flight time)")

# Maximum height (when vy = 0)
# At apex: vy - g*t = 0, so t = vy/g
t_max_height = vy / g
max_height = h0 + vy * t_max_height - 0.5 * g * t_max_height**2

# Range (horizontal distance)
range_distance = vx * t_flight

print(f"\nTime of flight: {t_flight:.2f} s")
print(f"Maximum height: {max_height:.2f} m")
print(f"Range: {range_distance:.2f} m")

# ==========================================
# SIMULATE TRAJECTORY
# ==========================================

# Create time array
dt = 0.01  # Time step (s)
time_array = np.arange(0, t_flight + dt, dt)

# Calculate position at each time step
x_positions = vx * time_array
y_positions = h0 + vy * time_array - 0.5 * g * time_array**2

# Filter out negative heights (after landing)
valid_indices = y_positions >= 0
x_positions = x_positions[valid_indices]
y_positions = y_positions[valid_indices]
time_array = time_array[valid_indices]

print(f"\nSimulated {len(time_array)} trajectory points")

# ==========================================
# ANALYZE SPECIFIC TIME POINTS
# ==========================================

def analyze_time(t):
    """Analyze projectile state at time t"""
    x = vx * t
    y = h0 + vy * t - 0.5 * g * t**2
    
    vx_t = vx  # Constant
    vy_t = vy - g * t
    
    speed = np.sqrt(vx_t**2 + vy_t**2)
    angle_t = np.degrees(np.arctan2(vy_t, vx_t))
    
    print(f"\n--- Time t = {t:.2f} s ---")
    print(f"Position: ({x:.2f}, {y:.2f}) m")
    print(f"Velocity: ({vx_t:.2f}, {vy_t:.2f}) m/s")
    print(f"Speed: {speed:.2f} m/s")
    print(f"Angle: {angle_t:.1f}°")
    
    return x, y, vx_t, vy_t, speed

# Analyze at launch, apex, and landing
print("\n" + "="*40)
print("TRAJECTORY ANALYSIS")
print("="*40)

analyze_time(0)  # Launch
analyze_time(t_max_height)  # Apex
analyze_time(t_flight)  # Landing

# ==========================================
# ENERGY ANALYSIS
# ==========================================

def calculate_energy(t, mass=1.0):
    """Calculate kinetic and potential energy"""
    y = h0 + vy * t - 0.5 * g * t**2
    vy_t = vy - g * t
    
    KE = 0.5 * mass * (vx**2 + vy_t**2)
    PE = mass * g * y
    total_E = KE + PE
    
    return KE, PE, total_E

print("\n" + "="*40)
print("ENERGY CONSERVATION (m = 1 kg)")
print("="*40)

for t in [0, t_max_height, t_flight]:
    KE, PE, E_total = calculate_energy(t)
    print(f"\nt = {t:.2f} s:")
    print(f"  KE = {KE:.2f} J")
    print(f"  PE = {PE:.2f} J")
    print(f"  Total = {E_total:.2f} J")

# ==========================================
# VISUALIZATION
# ==========================================

plt.figure(figsize=(12, 8))

# Subplot 1: Trajectory
plt.subplot(2, 2, 1)
plt.plot(x_positions, y_positions, 'b-', linewidth=2, label='Trajectory')
plt.plot([x_positions[0]], [y_positions[0]], 'go', markersize=10, label='Launch')
plt.plot([x_positions[-1]], [y_positions[-1]], 'ro', markersize=10, label='Landing')
plt.xlabel('Distance (m)')
plt.ylabel('Height (m)')
plt.title('Projectile Trajectory')
plt.grid(True, alpha=0.3)
plt.legend()
plt.axis('equal')

# Subplot 2: Height vs Time
plt.subplot(2, 2, 2)
plt.plot(time_array, y_positions, 'g-', linewidth=2)
plt.axhline(y=max_height, color='r', linestyle='--', alpha=0.5, label=f'Max: {max_height:.1f} m')
plt.xlabel('Time (s)')
plt.ylabel('Height (m)')
plt.title('Height vs Time')
plt.grid(True, alpha=0.3)
plt.legend()

# Subplot 3: Velocity Components
plt.subplot(2, 2, 3)
vx_array = np.full_like(time_array, vx)
vy_array = vy - g * time_array
plt.plot(time_array, vx_array, 'b-', linewidth=2, label='vx (horizontal)')
plt.plot(time_array, vy_array, 'r-', linewidth=2, label='vy (vertical)')
plt.axhline(y=0, color='k', linestyle='-', alpha=0.3)
plt.xlabel('Time (s)')
plt.ylabel('Velocity (m/s)')
plt.title('Velocity Components')
plt.grid(True, alpha=0.3)
plt.legend()

# Subplot 4: Energy
plt.subplot(2, 2, 4)
KE_array = 0.5 * (vx**2 + (vy - g * time_array)**2)
PE_array = g * y_positions
E_total_array = KE_array + PE_array

plt.plot(time_array, KE_array, 'b-', linewidth=2, label='Kinetic Energy')
plt.plot(time_array, PE_array, 'r-', linewidth=2, label='Potential Energy')
plt.plot(time_array, E_total_array, 'k--', linewidth=2, label='Total Energy')
plt.xlabel('Time (s)')
plt.ylabel('Energy (J/kg)')
plt.title('Energy vs Time')
plt.grid(True, alpha=0.3)
plt.legend()

plt.tight_layout()
plt.show()

# ==========================================
# OPTIMAL ANGLE FOR MAX RANGE
# ==========================================

print("\n" + "="*40)
print("OPTIMAL LAUNCH ANGLE")
print("="*40)

angles = np.linspace(0, 90, 91)
ranges = []

for ang in angles:
    ang_rad = np.radians(ang)
    vx_temp = v0 * np.cos(ang_rad)
    vy_temp = v0 * np.sin(ang_rad)
    
    disc = vy_temp**2 + 2*g*h0
    if disc >= 0:
        t_temp = (vy_temp + np.sqrt(disc)) / g
        r_temp = vx_temp * t_temp
        ranges.append(r_temp)
    else:
        ranges.append(0)

max_range_idx = np.argmax(ranges)
optimal_angle = angles[max_range_idx]
max_range = ranges[max_range_idx]

print(f"For v0 = {v0} m/s and h0 = {h0} m:")
print(f"Optimal angle: {optimal_angle:.1f}°")
print(f"Maximum range: {max_range:.2f} m")
print(f"\nNote: For h0 = 0, optimal angle is 45°")