# Projectile motion with drag. My problem is in the loop as the number of elements are not equal apparently on the left and right. Id appreciate detailed input ive been lost on this for a couple of days now

조회 수: 17(최근 30일)
편집: James Tursa 22 Dec 2020
function Yoted
clear
clc
format long g;
workspace
format compact
fontSize = 20;
g = 9.81; % acceleration due to gravity in Y
y0 = 0; % intial position y
vHit = 50; % initial velocity magnitude
v0Max = vHit/10 + vHit;
v0Min = vHit - vHit/10;
v0 = randi([v0Min v0Max],1) %randomly generated v0
wHit = 100; % Angular Velocity (rad/s)
wMin = wHit - wHit/10;
wMax = wHit + wHit/10;
Angular_w = randi([wMin wMax],1) % randomly generate w
Angle = 30; % angle of strike
aMin = Angle - Angle/10;
aMax = Angle + Angle/10;
Alpha = randi([aMin aMax],1) % Angle Alpha hit
c = 4.5*10^(-5); % Coefficient of drag
r = 21.5*10^(-3); % Radius of gold ball
m = 46*10^(-3); % Mass of golf ball
rho_air = 1.2; % Density of air
p = [0 0 0];
dt=0.1;
t=[0.1:dt:200];
v = v0; % Velocity Magnitude
Velx = v(1)*cosd(Alpha)
Vely = v(1)*sind(Alpha)- g*t(1)
Vel(1,:) = [Velx Vely 0] % Velocity Vector
FdragX(1) = v(1).*Vel(1).*c; % Drag force in X direction
FdragY(1) = v(1).*Vel(2).*c; % Drag force in Y direction
AccDragX(1) = -FdragX(1)./ m; % F=M*A decceleration due to drag in X
AccDragY(1) = -FdragY(1)./ m; % F=M*A decceleration due to drag in Y
AccDrag(1,:) = [AccDragX(1) AccDragY(1) 0]
xImpact_no_Drag = (sind(2*Alpha)*(v0^2))/g
t(5)
for i=2:length(t)
p(i) = p(i-1) + Vel.*t(i-1) + AccDrag.*(t(i-1).^2);
Vel = Vel + AccDrag.*t(i-1)
v= sqrt(sum((Vel.^2),"all"))
Velx = v*cosd(Alpha)
Vely = v*sind(Alpha)- g.*t(i)
Vel = [Velx Vely 0];
FdragX = v.*Velx*c; % Drag force in X direction
FdragY = v.*Vely*c; % Drag force in Y direction
AccDragX = -FdragX./ m; % F=M*A decceleration due to drag in X
AccDragY = -FdragY./ m; % F=M*A decceleration due to drag in Y
AccDrag = [AccDragX AccDragY 0]
if p<0
break
end
end
plot(t,p)

댓글을 달려면 로그인하십시오.

### 채택된 답변

James Tursa 22 Dec 2020
편집: James Tursa 22 Dec 2020
This makes AccDrag a vector:
AccDrag(1,:) = [AccDragX(1) AccDragY(1) 0]
So here you have a vector on the right hand side and a scalar on the left hand side:
p(i) = p(i-1) + Vel.*t(i-1) + AccDrag.*(t(i-1).^2);
Looks like you could fix this with:
p(i,:) = p(i-1,:) + etc.
BUT, the above equation is wrong. The right hand side needs to be current position + delta position. But you have used formulas that include the absolute time, not the delta time. This will generate an incorrect answer. You need to change how you are updating the position and velocity to use delta times, not absolute times.
Also this test:
if p<0
should be this (if you use the fix above):
if p(i,2)<0

댓글을 달려면 로그인하십시오.

### Community Treasure Hunt

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

Start Hunting!

Translated by