Hey @hawa hassan, please refer to below code where i have plotted the above reqirements, take it as reference and modify accordingly:
% ProjectileMotionWithAirResistance.m
% This script simulates projectile motion with air resistance using ode45.
clear; close all; clc;
%% 1. Define Problem Parameters and Initial Conditions
% Your ID and angle calculation
ID = [1, 8, 3, 9, 8, 1, 0];
theta_degrees = sum(ID); % Sum of ID is 30, so launch angle is 30 degrees
% Initial conditions for the projectile
initial_x = 0; % Initial horizontal position
initial_y = mean(ID); % Initial vertical position (mean(ID) = 30/7 = 4.2857)
% Assume an initial speed. Your original code implies an initial speed of 1.
% For a more visible trajectory, a higher speed is usually needed.
initial_speed = 10; % meters/second (adjust as needed for your problem)
% Calculate initial velocity components based on initial_speed and theta
initial_vx = initial_speed * cosd(theta_degrees);
initial_vy = initial_speed * sind(theta_degrees);
% Air resistance coefficient (c)
% This value determines the strength of air resistance.
% If c=0, there is no air resistance. A small positive value will show its effect.
c_air_resistance = 0.01; % Example coefficient (adjust based on your problem's specifics)
% Gravity (g)
% Your original dydt had '-1' for gravity. This implies g=1 acting downwards.
% For real-world physics, use 9.81. We'll stick to 1 to match your setup.
gravity_val = 1; % Value for gravity (e.g., 9.81 for Earth, or 1 as in your original code)
% Initial state vector y0 = [x; vx; y; vy]
y0 = [initial_x; initial_vx; initial_y; initial_vy];
% Time span for the simulation
% We'll simulate for a sufficiently long time and then trim the data
% once the projectile hits the ground.
tspan = [0 20]; % Start at t=0, simulate up to 20 seconds (adjust if needed)
%% 2. Solve the Ordinary Differential Equation (ODE)
% Call ode45. We use an anonymous function @(t,y_vec) to pass additional
% parameters (c_air_resistance, gravity_val) to our 'projectile' function.
[t, y] = ode45(@(t, y_vec) projectile(t, y_vec, c_air_resistance, gravity_val), tspan, y0);
%% 3. Post-processing: Trim data to when projectile hits the ground
% Find the first index where vertical position (y) is less than or equal to 0
% (assuming y(3) is the vertical position)
ground_hit_idx = find(y(:,3) <= 0, 1, 'first');
if isempty(ground_hit_idx)
warning('Projectile did not hit the ground within the simulated time span. Consider increasing tspan.');
t_plot = t;
y_plot = y;
else
% Trim the time and state vectors to the point of impact
t_plot = t(1:ground_hit_idx);
y_plot = y(1:ground_hit_idx,:);
end
%% 4. Plotting the Results
% Plot 1: Trajectory (y vs x)
figure;
plot(y_plot(:,1), y_plot(:,3), 'b-', 'LineWidth', 2);
hold on;
plot(y_plot(1,1), y_plot(1,3), 'go', 'MarkerSize', 8, 'MarkerFaceColor', 'g', 'DisplayName', 'Start Point'); % Start point
plot(y_plot(end,1), y_plot(end,3), 'rx', 'MarkerSize', 8, 'LineWidth', 2, 'DisplayName', 'End Point'); % End point
grid on;
xlabel('Horizontal Distance (x)');
ylabel('Vertical Height (y)');
title('Projectile Motion Trajectory with Air Resistance');
legend('Trajectory', 'Start Point', 'End Point', 'Location', 'best');
axis tight; % Adjust axes to fit the data
hold off;
% Plot 2: State Variables vs. Time
figure;
plot(t_plot, y_plot(:,1), 'DisplayName', 'x (Horizontal Position)');
hold on;
plot(t_plot, y_plot(:,2), 'DisplayName', 'vx (Horizontal Velocity)');
plot(t_plot, y_plot(:,3), 'DisplayName', 'y (Vertical Position)');
plot(t_plot, y_plot(:,4), 'DisplayName', 'vy (Vertical Velocity)');
grid on;
xlabel('Time (s)');
ylabel('Value');
title('Projectile State Variables vs. Time');
legend('Location', 'best');
hold off;
% Display key results
disp(['Initial Angle (theta): ', num2str(theta_degrees), ' degrees']);
disp(['Initial Speed (V0): ', num2str(initial_speed), ' m/s']);
disp(['Air Resistance Coefficient (c): ', num2str(c_air_resistance)]);
disp(['Initial Height: ', num2str(initial_y)]);
disp(['Total Flight Time: ', num2str(t_plot(end)), ' s']);
disp(['Total Horizontal Distance (Range): ', num2str(y_plot(end,1)), ' units']);
%% --- Projectile ODE Function ---
% Save this as 'projectile.m' in the same directory as your main script,
% or include it directly below the main script in the same .m file.
function dydt = projectile(t, y, c_coeff, g_val)
% PROJECTILE ODE function for projectile motion with air resistance.
% dydt = projectile(t, y, c_coeff, g_val) calculates the derivatives
% of the state vector y at time t.
%
% State vector y:
% y(1) = x (horizontal position)
% y(2) = vx (horizontal velocity)
% y(3) = y (vertical position)
% y(4) = vy (vertical velocity)
%
% Inputs:
% t : Current time (scalar)
% y : Current state vector (column vector)
% c_coeff : Air resistance coefficient (scalar)
% g_val : Value of gravitational acceleration (scalar, positive)
% Parameters
g = g_val; % Gravitational acceleration (e.g., 9.81 m/s^2)
c = c_coeff; % Air resistance coefficient
% Calculate the magnitude of the velocity vector
v = sqrt(y(2)^2 + y(4)^2); % v = sqrt(vx^2 + vy^2)
% Equations of motion (dydt)
% 1. dx/dt = vx
dxdt = y(2);
% 2. dvx/dt = ax (acceleration in x-direction due to air resistance)
% Air resistance force is proportional to v^2 and acts opposite to velocity.
% F_drag = -c * v^2 * (velocity_vector / v) = -c * v * velocity_vector
% ax = F_drag_x / m = -c * v * vx (assuming mass m=1 or c incorporates mass)
dvxdt = -c * y(2) * v;
% 3. dy/dt = vy
dydt_val = y(4);
% 4. dvy/dt = ay (acceleration in y-direction due to gravity and air resistance)
% ay = -g - c * v * vy (assuming mass m=1 or c incorporates mass)
dvydt = -g - (c * y(4) * v);
% Combine derivatives into a column vector
dydt = [dxdt; dvxdt; dydt_val; dvydt];
end
