I’m trying to simulate the flight of a glider I 3D printed and tested.
My model is the following :
The glider is launched with an initial velocity and angle of attack (θ in the code), the forces at play here are its weight, drag and lift, no thrust at all.
I’ve put in equation the model and here is what I found :

dX/dt = -B(X^2+Y^2)(Cl*sin(θ)+Cd*cos(θ) dY/dt = - g + B(X^2+Y^2)(Cl*cos(θ)-Cd*sin(θ) where X = Vx and Y = Vy, B = 0.5*µ*S/m and Cl, Cd the lift and drag coefficients

Now for the code :

First I import Cl/Cd data from a CSV file and interpolate them.
I define my dX/dt, dY/dt functions
I define my Euler function to solve

import numpy as np
import matplotlib.pyplot as plt
import scipy.interpolate as interpolate

# Data
rho = 1.2
m = 0.06
S = 0.05
g = 9.81

B = 0.5*rho*S*(1/m)

# Data on angle of attack, lift coefficient, and drag coefficient
angles = []  
CL_raw = []  
CD_raw = []  

file = open("C:/Cours Prépa/TIPE/xf-naca0015-il-50000.csv")
L = file.readlines()
for i in range(11, 100):
    angles.append(float(L[i].split(",")[0]))
    CL_raw.append(float(L[i].split(",")[1]))
    CD_raw.append(float(L[i].split(",")[2]))

# Interpolation of the data
f_CD = interpolate.interp1d(angles, CD_raw)
f_CL = interpolate.interp1d(angles, CL_raw)

# Creating arrays with higher resolution values
angles_interp = np.linspace(min(angles), max(angles), 1000)
CD_interp = f_CD(angles_interp)  # Using interpolation function for CD
CL_interp = f_CL(angles_interp)  # Using interpolation function for CL


# Functions to calculate derivatives dX/dt and dY/dt based on the inclination
def dX_dt(X, Y, θ):
    CD = f_CD(θ)
    CL = f_CL(θ)
    return -B*(X**2 + Y**2)*(CL*np.sin(np.radians(θ)) + CD*np.cos(np.radians(θ)))

def dY_dt(X, Y, θ):
    CD = f_CD(θ)
    CL = f_CL(θ)
    return -g + B*(X**2 + Y**2)*(CL*np.cos(np.radians(θ)) - CD*np.sin(np.radians(θ)))


# Main function for Euler's method
def euler(dt):
    # Initial conditions
    X = [15]  # Initial horizontal velocity in m/s
    Y = [0]   # Initial vertical velocity in m/s
    x = [0]   # Initial horizontal position in m
    y = [10]  # Initial vertical position in m
    θ = [np.rad2deg(np.arctan(Y[0]/X[0]))] # Initial angle of attack
    v = [np.sqrt(X[0]**2 + Y[0]**2)]
    elapsed_time = 0

    i = 0

    # Euler's method
    while y[i] >= 0:
        X.append(X[i] + dX_dt(X[i], Y[i], θ[i])*dt)
        Y.append(Y[i] + dY_dt(X[i], Y[i], θ[i])*dt)
        x.append(x[i] + X[i]*dt)
        y.append(y[i] + Y[i]*dt)
        v.append(np.sqrt(X[i]**2 + Y[i]**2))
        elapsed_time += dt
        θ.append(np.rad2deg(np.arctan(Y[i]/X[i])))
        if θ[i] > angle_limit or θ[i] < -angle_limit:
            print("Stall")
            return x, y, v, elapsed_time, θ
        i += 1

    return x, y, v, elapsed_time, θ

# Simulation parameters
dt = 0.001   # Time step in seconds
angle_limit = 10  # Stall angle limit in degrees


x, y, v, elapsed_time, θ = euler(dt)
print(elapsed_time)
time_array = np.arange(0, elapsed_time, dt)

# Plotting the graph
fig = plt.figure(figsize=(15, 10))

plt.subplot(3, 2, 1)
plt.plot(angles, CD_raw, 'o', label='Raw data')
plt.plot(angles_interp, CD_interp, '-', label='Interpolation')
plt.xlabel('Angle of Attack (degrees)')
plt.ylabel('Drag Coefficient (CD)')
plt.title('Variation of CD with Angle of Attack')
plt.legend()
plt.grid(True)

plt.subplot(3, 2, 2)
plt.plot(angles, CL_raw, 'o', label='Raw data')
plt.plot(angles_interp, CL_interp, '-', label='Interpolation')
plt.xlabel('Angle of Attack (degrees)')
plt.ylabel('Lift Coefficient (CL)')
plt.title('Variation of CL with Angle of Attack')
plt.legend()
plt.grid(True)

plt.subplot(3, 2, 3)
plt.plot(x, y)
plt.ylim([0, 11])
plt.xlabel('x (m)')
plt.ylabel('y (m)')
plt.title('Aircraft Trajectory')
plt.grid(True)

# Plotting the graph of the aircraft's inclination relative to the ground
plt.subplot(3, 2, 4)
plt.plot(time_array, θ)
plt.xlabel('Time')
plt.ylabel('Inclination (degrees)')
plt.title('Aircraft Inclination Relative to Ground')
plt.grid(True)

plt.subplot(3, 2, 5)
plt.plot(time_array, v)
plt.xlabel('Time')
plt.ylabel('Speed magnitude v')
plt.title('Speed Magnitude Over Time')
plt.grid(True)

plt.tight_layout()
plt.show()

When I run the code, the angle of attack suddenly drops no matter what I try.
I tried a lot of initial conditions and time steps etc… but no matter what I do the angle of attack always diverges rapidly, which is not what happens in real testing

Let me know if you have any idea
Thanks

When I run the code, the angle of attack suddenly drops no matter what I try.
I tried a lot of initial conditions and time steps etc… but no matter what I do the angle of attack always diverges rapidly, which is not what happens in real testing

Let me know if you have any idea
Thanks