from math import radians, cos, sin, asin, sqrt
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams.update({'font.size': 14})
plt.rc('text', usetex=True)
plt.rc('text.latex', preamble=r'\usepackage{amsmath}')
plt.rcParams["font.family"] = "Times New Roman"

###############################################################################
# Distance London to Sydney                                                   #
###############################################################################

# https://www.latlong.net/
coord_london = (51.470020,-0.454295) # Heathrow Airport, London, UK
coord_sydney = (-33.947346,151.179428) # Sydney Airport, NSW, Australia

def haversine(coord1, coord2):
    """source: https://stackoverflow.com/a/4913653"""
    lon1, lat1, lon2, lat2 = map(radians, [coord1[1], coord1[0], coord2[1], coord2[0]])
    dlon = lon2 - lon1
    dlat = lat2 - lat1
    a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2
    c = 2 * asin(sqrt(a))
    r = 6371 # Earth's radius (km)
    return c * r

distance = haversine(coord_london, coord_sydney) # km
print(f"Distance from London to Sydney: {distance:.2f} km")


###############################################################################
# Travel time                                                                 #
###############################################################################

a = 1.23 # maximum comfortable m/s^2 (https://doi.org/10.1016/j.apergo.2022.103881, p.8)
distance *= 1000 # convert distance to meters
t = 0 # s
v = 0 # m/s
s = 0 # m
data = [] # (velocity, distance) for plotting

# Acceleration
while s < distance/2:
    t += 1
    v += a
    s += v
    data.append((v*3.6, s/1000))

print(f"Maximum velocity: {v*3.6:.2f} km/h at {t/3600:.2f} h ({t:.0f} s)")

# Deceleration
while s < distance:
    t += 1
    v -= a
    s += v
    data.append((v*3.6, s/1000))

print(f"Travel time from London to Sydney: {t/3600:.2f} hours ({t:.0f} s)")


###############################################################################
# Plotting                                                                    #
###############################################################################

fig, axs = plt.subplots(2, 2, figsize=(12, 6))
fig.suptitle("Travel from London to Sydney")

# Plot velocity vs time
axs[0,0].plot(np.arange(0, t, 1), list(zip(*data))[0], 'b-')
axs[0,0].set_xlabel('Time [s]')
axs[0,0].set_ylabel('Velocity [km/h]')
axs[0,0].set_ylim(top=17500)
axs[0,0].set_yticks(np.arange(0, 17501, 2500))
axs[0,0].set_xlim(right=8000)
axs[0,0].set_xticks(np.arange(0, 8001, 1000))
axs[0,0].grid(True, axis='both', which='both', linestyle='--', linewidth=0.5)

# Plot distance vs velocity
axs[0,1].plot(list(zip(*data))[1], list(zip(*data))[0], 'r-')
axs[0,1].set_xlabel('Distance [km]')
axs[0,1].set_ylabel('Velocity [km/h]')
axs[0,1].set_ylim(top=17500)
axs[0,1].set_yticks(np.arange(0, 17501, 2500))
axs[0,1].set_xlim(right=17500)
axs[0,1].set_xticks(np.arange(0, 17501, 2500))
axs[0,1].grid(True, axis='both', which='both', linestyle='--', linewidth=0.5)

# Plot acceleration vs time
axs[1,0].plot(np.arange(0, round(t/2), 1), np.full(round(t/2), a), 'g-')
axs[1,0].plot(np.arange(round(t/2), t, 1), np.full(round(t/2)+1, -a), 'g-')
axs[1,0].set_xlabel('Time [s]')
axs[1,0].set_ylabel('Acceleration [m/s$^2$]')
axs[1,0].set_ylim((-1.5,1.5))
axs[1,0].set_yticks(np.arange(-1.5, 1.51, 0.5))
axs[1,0].set_xlim(right=8000)
axs[1,0].set_xticks(np.arange(0, 8001, 1000))
axs[1,0].grid(True, axis='both', which='both', linestyle='--', linewidth=0.5)

# Plot distance vs time
axs[1,1].plot(np.arange(0, t, 1), list(zip(*data))[1], 'm-')
axs[1,1].set_xlabel('Time [s]')
axs[1,1].set_ylabel('Distance [km]')
axs[1,1].set_ylim(top=17500)
axs[1,1].set_yticks(np.arange(0, 17501, 2500))
axs[1,1].set_xlim(right=8000)
axs[1,1].set_xticks(np.arange(0, 8001, 1000))
axs[1,1].grid(True, axis='both', which='both', linestyle='--', linewidth=0.5)

plt.tight_layout()
plt.subplots_adjust(top=0.9)
plt.savefig("distance_london_sydney.png", dpi=300)
plt.show()