#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed May  7 10:25:47 2025

@author: christian
"""

import numpy as np
import matplotlib.pyplot as plt

# Initiate vectors
t = np.linspace(0, 30, 5)  # simulate for 20 seconds
v = 42.78 * (1 - np.exp(-9.81 * t / 42.78))

# Forward difference method
dvdt_fd = (v[1:] - v[:-1]) / (t[1:] - t[:-1])
dvdt_bd = (v[1:] - v[:-1]) / (t[1:] - t[:-1])
dvdt_cd = (v[2:] - v[:-2]) / (t[2:] - t[:-2])

# Plotting
plt.figure(figsize=(8, 5))
plt.plot(t, v, label='Velocity (m/s)', color='blue')
plt.plot(t[:-1], dvdt_fd, label=r'forward difference')
plt.plot(t[1:], dvdt_bd, label=r'backward difference')
plt.plot(t[1:-1], dvdt_cd, label=r'central difference')
plt.axhline(y=42.78, color='red', linestyle='--', label=f'Terminal velocity ≈ {42.78:.2f} m/s')
plt.title('Falling Object Approaching Terminal Velocity')
plt.xlabel('Time (s)')
plt.ylabel('Velocity (m/s)')
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()