#QUADRATIC FRICTION - EULER METHOD (Security parable) 
import numpy as np
import matplotlib.pyplot as plt
from math import cos, sin, pi

# Parámetros
m=1  # masa kg
b=0.1  # coeficiente de arrastre
g=9.81  # aceleración debido a la gravedad
vo=100  # velocidad m/s

# Tiempo
t0=0
tfinal=10
dt=0.001
t=np.arange(t0,tfinal,dt)

for ꞵ in range(5,95,5):
    θ=ꞵ*pi/180  # Convertir ángulo a radianes

    # Condiciones iniciales
    x0=0
    y0=0
    vx0=vo*cos(θ)
    vy0=vo*sin(θ)

    # Arrays para almacenar las soluciones
    x=np.zeros_like(t)
    y=np.zeros_like(t)
    vx=np.zeros_like(t)
    vy=np.zeros_like(t)

    # Condiciones iniciales
    x[0]=x0
    y[0]=y0
    vx[0]=vx0
    vy[0]=vy0

    # Resolver las ecuaciones usando el método de Euler
    for i in range(1,len(t)):
        # Calcular aceleraciones
        ax = -b*vx[i-1]*np.sqrt(vx[i-1]**2 + vy[i-1]**2)/m
        ay = -g-b*vy[i-1]*np.sqrt(vx[i-1]**2 + vy[i-1]**2)/m
        
        # Actualizar velocidades
        vx[i]=vx[i-1]+ax*dt
        vy[i]=vy[i-1]+ay*dt
        
        # Actualizar posiciones
        x[i]=x[i-1]+vx[i]*dt
        y[i]=y[i-1]+vy[i]*dt

        if y[i]<0:  # Si el proyectil toca el suelo, detener la simulación
            break

    # Graficar la trayectoria
    plt.plot(x[:i],y[:i],color='blue')

# Configuración de la gráfica
plt.xlabel('$x(m)$',fontsize=12)
plt.ylabel('$y(m)$',fontsize=12)
plt.grid(True)

# Ajustar el tamaño de la letra de los ticks
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
plt.show()