3 views (last 30 days)

Show older comments

Hi all, I'm working on a balistic computer of sorts that will take in the initial point, the impact point, and the launch angle and will solve for the required launch velocity. Right now I'm trying to get the trejectory to match expectations without considering drag, but next steps include implementing drag and recursively finding the angle with the minimum launch velocity (which is why theta is listed as an additional parameter.

Lastly, before I get to the code, I'm fairly convinced that the problem is coming from the boundary conditions because if you look at the displayed start and end points, the start point is different than the intended start point and the angle is very different than the one that should be enforced by the fifth boundary condition. The solinit "guess" seems to have very little impact on the solution, and I'm confident with my odesys equations.

Thank you in advance for your help!

d = 100;

theta = pi/4;

solinit = bvpinit(linspace(0,d,20),@guess, theta);

sol = bvp4c(@odesys, @resfun, solinit);

y=deval(sol,linspace(1,d,100));

plot(y(3,:),y(1,:));

fprintf('v_0 (from height) = %f \ntheta = %f \n(%f,%f) -> (%f,%f)\n', ...

sqrt(max(y(1,:))* 32.2 * 2), atand(y(2,1)/y(4,1)), y(3,1), y(1,1),y(3,end),y(1,end));

function dy = odesys(x, y, theta)

% launch parameters

g = 32.2; % gravity, ft/s^2

Cd = 0; % drag coefficient

rho = 0.00237; % density of air, slug/ft^3

c = 10; % pumpkin circumfrence (normal to motion), in

c = c / 12; % convert in to ft

r = c / 2 / pi; % radius, ft

A = pi * r ^ 2; % cross sectional area, ft^2

W = 4; % weight, lbf

m = W / g; % convert lbf to slug

% constants for use in defining dy

k = rho * Cd * A / 2 / m; % constant for drag, 1/sqrt(s*ft)

v = sqrt(y(2).^2+y(4).^2);

%defining dy

dy = y;

dy(1) = y(2);

dy(2) = -g + k * v * y(2);

dy(3) = y(4);

dy(4) = k * v * y(4);

end

function res = resfun(Ya, Yb, theta)

d = 100; % target distance, ft

% 50, 100, 200, and 300

h = 4; % launch height, ft

res = zeros(5,1);

% point a: (0,h)

res(1) = Ya(1)-h;

res(2) = Ya(3);

% point b: (d,0)

res(3) = Yb(1);

res(4) = Yb(3)-d;

%initial launch angle

res(5) = theta - asin(Ya(2)/sqrt(Ya(2).^2+Ya(4).^2));

end

function y = guess(x)

% calculate minimum launch velocity with 0 drag

g = 32.2;

sine = sind(45);

d = 100;

h = 4;

v = d / sine / sqrt( 2 * ( d + h) / g);

% guess follows 0 drag 45º launch

y(1) = -g*x^2/v^2/sine^2+x+h;

y(2) = -g*x/v/sine+v*sine;

y(3) = x;

y(4) = v*sine;

end

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!