Trajectory of a 3d body

Problem Statement: –
I have a 3d body which i have launched from cruising aircraft travelling at 200m/sec at 10km altitude. The launch speed is 50m/sec and it is launched at 45degrees angle pointing downwards. I want to create the whole trajectory of the body using everything like change in pressure, density, drag, lift etc. I created a code but it is not giving the accurate answer.

the code: –
function aircraft_launch_trajectory_6DOF_advanced()
% Constants
g = 9.81; % Acceleration due to gravity (m/s^2)

% Convert given units to SI
m = 49.964 * 0.0283495; % Mass (kg) (1 ouncemass = 0.0283495 kg)
A = 68.166 * 0.00064516; % Area (m^2) (1 in^2 = 0.00064516 m^2)
Izz = 19.58 * 0.0283495 * 0.0254^2; % Moment of inertia (kg*m^2)

% Body dimensions
length = 1.80 * 0.0254; % Length (m)
width = 3.30 * 0.0254; % Width (m)
height = 8.00 * 0.0254; % Height (m)

% Center of mass location (assuming from the base of the body)
cm_z = 1.708 * 0.0254; % (m)

% Assumed center of pressure location (you may need to adjust this)
cp_z = 4.5 * 0.0254; % (m)

% Distance between center of mass and center of pressure
d_cp_cm = cp_z – cm_z;

% Initial conditions
h0 = 10000; % Initial height (m)
v0_aircraft = 200; % Aircraft velocity (m/s)
v0_launch = 50; % Launch velocity (m/s)
angle = -45; % Launch angle (degrees)
omega0 = 0; % Initial angular velocity (rad/s)

% Calculate initial velocity components
v0x = v0_aircraft – v0_launch * cosd(angle);
v0y = -v0_launch * sind(angle);

% Initial state vector [x, y, vx, vy, theta, omega]
y0 = [0; h0; v0x; v0y; deg2rad(angle); omega0];

% Time settings
tspan = [0 300]; % Time span for integration

% ODE solver options
options = odeset(‘Events’, @ground_impact, ‘RelTol’, 1e-6, ‘AbsTol’, 1e-6);

% Solve ODE using ode45 (4th order Runge-Kutta)
[t, y] = ode45(@(t, y) equations_of_motion_6DOF_advanced(t, y, m, Izz, A, d_cp_cm, length), tspan, y0, options);

% Extract results
x = y(:, 1);
h = y(:, 2);
vx = y(:, 3);
vy = y(:, 4);
theta = y(:, 5);
omega = y(:, 6);

% Calculate velocity magnitude
v = sqrt(vx.^2 + vy.^2);

% Plot trajectory
figure;
plot(x/1000, h/1000);
xlabel(‘Horizontal distance (km)’);
ylabel(‘Altitude (km)’);
title(‘2D Trajectory of Body (Advanced 6DOF Simulation)’);
grid on;

% Plot velocity
figure;
plot(t, v);
xlabel(‘Time (s)’);
ylabel(‘Velocity (m/s)’);
title(‘Velocity vs Time’);
grid on;

% Plot angular position and velocity
figure;
subplot(2,1,1);
plot(t, rad2deg(theta));
ylabel(‘Angular Position (degrees)’);
title(‘Angular Motion’);
grid on;
subplot(2,1,2);
plot(t, rad2deg(omega));
xlabel(‘Time (s)’);
ylabel(‘Angular Velocity (deg/s)’);
grid on;

% Display results
disp([‘Time of flight: ‘, num2str(t(end)), ‘ s’]);
disp([‘Horizontal distance traveled: ‘, num2str(x(end)/1000), ‘ km’]);
disp([‘Impact velocity: ‘, num2str(v(end)), ‘ m/s’]);
disp([‘Final angle: ‘, num2str(rad2deg(theta(end))), ‘ degrees’]);
disp([‘Max altitude: ‘, num2str(max(h)/1000), ‘ km’]);
disp([‘Max velocity: ‘, num2str(max(v)), ‘ m/s’]);

% Output state at specific times for debugging
debug_times = [0, 5, 10, 30, 60, 120];
for debug_time = debug_times
idx = find(t >= debug_time, 1);
if ~isempty(idx)
disp([‘State at t = ‘, num2str(t(idx)), ‘ s:’]);
disp([‘ Position: (‘, num2str(x(idx)/1000), ‘ km, ‘, num2str(h(idx)/1000), ‘ km)’]);
disp([‘ Velocity: ‘, num2str(v(idx)), ‘ m/s’]);
disp([‘ Angle: ‘, num2str(rad2deg(theta(idx))), ‘ degrees’]);
end
end
end

function dydt = equations_of_motion_6DOF_advanced(t, y, m, Izz, A, d_cp_cm, length)
% Unpack state vector
x = y(1);
h = y(2);
vx = y(3);
vy = y(4);
theta = y(5);
omega = y(6);

% Calculate velocity magnitude and angle of attack
v = sqrt(vx^2 + vy^2);
alpha = atan2(vy, vx) – theta;

% Calculate air density and speed of sound using improved atmospheric model
[rho, a] = atmosphere_model(h);

% Calculate Mach number
M = v / a;

% Get aerodynamic coefficients based on angle of attack and Mach number
[Cd, Cl, Cm] = aero_coefficients(alpha, M);

% Calculate aerodynamic forces and moment
q = 0.5 * rho * v^2; % Dynamic pressure
Fd = q * Cd * A; % Drag force
Fl = q * Cl * A; % Lift force
M_aero = q * Cm * A * length; % Aerodynamic moment

% Add moment due to center of pressure offset
M_total = M_aero + Fl * d_cp_cm;

% Apply initial boost for first 7 seconds
boost_duration = 7; % seconds
if t <= boost_duration
boost_force = m * 200 / boost_duration; % Force to achieve 200 m/s in 7 seconds
boost_fx = boost_force * cos(theta);
boost_fy = boost_force * sin(theta);
else
boost_fx = 0;
boost_fy = 0;
end

% Calculate accelerations
ax = (-Fd * cos(theta) + Fl * sin(theta) + boost_fx) / m;
ay = -9.81 + (-Fd * sin(theta) – Fl * cos(theta) + boost_fy) / m;
alpha_dot = M_total / Izz; % Angular acceleration

% Return state derivatives
dydt = [vx; vy; ax; ay; omega; alpha_dot];
end

function [rho, a] = atmosphere_model(h)
% Constants
R = 287.05; % Specific gas constant for air (J/(kg*K))
g0 = 9.80665; % Standard gravity (m/s^2)
gamma = 1.4; % Ratio of specific heats for air

% Sea level conditions
T0 = 288.15; % Temperature at sea level (K)
P0 = 101325; % Pressure at sea level (Pa)

if h < 11000 % Troposphere
T = T0 – 0.0065 * h;
P = P0 * (T / T0)^(g0 / (R * 0.0065));
elseif h < 20000 % Lower Stratosphere
T = 216.65;
P = 22632.1 * exp(-g0 * (h – 11000) / (R * T));
else % Upper Stratosphere
T = 216.65 + 0.001 * (h – 20000);
P = 5474.89 * (216.65 / T)^(34.1632);
end

rho = P / (R * T);
a = sqrt(gamma * R * T); % Speed of sound
end

function [Cd, Cl, Cm] = aero_coefficients(alpha, M)
% This function returns aerodynamic coefficients based on angle of attack and Mach number
% Note: These are simplified models and should be replaced with actual data or more complex models

% Convert alpha to degrees for easier calculations
alpha_deg = rad2deg(alpha);

% Drag coefficient
Cd0 = 0.1 + 0.1 * M; % Base drag increases with Mach number
Cdi = 0.01 * alpha_deg^2; % Induced drag
Cd = Cd0 + Cdi;

% Lift coefficient
Cl_slope = 0.1; % Lift curve slope
Cl = Cl_slope * alpha_deg * (1 – M^2/4)^(-0.5); % Prandtl-Glauert correction for compressibility

% Moment coefficient
Cm0 = -0.05; % Zero-lift moment coefficient
Cm_alpha = -0.01; % Moment curve slope
Cm = Cm0 + Cm_alpha * alpha_deg;

% Limit coefficients to realistic ranges
Cd = max(0, min(Cd, 2));
Cl = max(-2, min(Cl, 2));
Cm = max(-1, min(Cm, 1));
end

function [position, isterminal, direction] = ground_impact(t, y)
position = y(2); % Height
isterminal = 1; % Stop integration when ground is reached
direction = -1; % Negative direction (approaching ground)
endProblem Statement: –
I have a 3d body which i have launched from cruising aircraft travelling at 200m/sec at 10km altitude. The launch speed is 50m/sec and it is launched at 45degrees angle pointing downwards. I want to create the whole trajectory of the body using everything like change in pressure, density, drag, lift etc. I created a code but it is not giving the accurate answer.