How to properly create for loop with if statement

7 次查看（过去 30 天）

Light 2023-4-2

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1939384-how-to-properly-create-for-loop-with-if-statement

编辑： Torsten 2023-4-22

for loop.png

Hi, how do you create a for or while loop with an if statement which adds 360 to an existing variable l, if l is negative?

采纳的回答

Star Strider 2023-4-2

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1939384-how-to-properly-create-for-loop-with-if-statement#answer_1206604

在 MATLAB Online 中打开

It is necessary to use the index on ‘L’ on both sides of the equation.

L = -180:20:180
L = 1×19
  -180  -160  -140  -120  -100   -80   -60   -40   -20     0    20    40    60    80   100   120   140   160   180
for l = 1:length(L)
    if L(l) < 0
        L(l) = L(l) + 360;
    end
end
L
L = 1×19
   180   200   220   240   260   280   300   320   340     0    20    40    60    80   100   120   140   160   180

A one-line version (avoiding the loop) would be —

L = -180:20:180
L = 1×19
  -180  -160  -140  -120  -100   -80   -60   -40   -20     0    20    40    60    80   100   120   140   160   180
L(L<0) = L(L<0) + 360
L = 1×19
   180   200   220   240   260   280   300   320   340     0    20    40    60    80   100   120   140   160   180

The one-line version is more efficient.

91 个评论
显示 89更早的评论隐藏 89更早的评论

Light 2023-4-2

编辑：Walter Roberson 2023-4-12

在 MATLAB Online 中打开

IMG-20230331-WA0004.jpg

Okay, I will paste it:

% l= longitude, b= latitude
l1 = -96.81
b1 = 32.78
% translating polar coordinates of the first city (Dallas) to cartesian coordinates
x1 = cosd (-96.81)*cosd (32.78)
y1 = sind (-96.81)*cosd (32.78)
z1 = sind (32.78)
l2 = 14.42
b2 = 50.07
% translating polar coordinates of the second city (Prague) to cartesian coordinates
x2 = cosd (14.42)*cosd (50.07)
y2 = sind (14.42)*cosd (50.07)
z2 = sind (50.07)
% calculating angular distance between the two cities (Dallas and Prague)
psi = acosd ((x1*x2)+(y1*y2)+(z1*z2))
x3 = ((x2 - x1)*cosd(psi))/ sind(psi)
y3 = ((y2 - y1)*cosd(psi))/ sind(psi)
z3 = ((z2 - z1)*cosd(psi))/ sind(psi)
% calculating the distance between the two cities (Dallas and Prague)
r = 6371 Distance = 2*pi*r*(psi/360)
% Plotting the great circle path between the two cities
phi = linspace(0, psi, 100)
    X = x1*cos(phi) + x3*sin(phi);
    Y = y1*cos(phi) + y3*sin(phi);
    Z = z1*cos(phi) + z3*sin(phi);
    % converting cartesian coordinates x, y, z to polar coordinates l, b % (longitude and latitude)
    b = asind(Z)
    l = atan2d (Y,X)
    l(l < 0) = l(l < 0) + 360
    load('topo.mat', 'topo', 'topomap1');
    % load earth topography
    whos topo topomap1;
    % set coordinates of points
    Longitude = 0:15:360; Latitude = -72:6:72;
    % set background parameters
    contour(0:359, -89:90, topo, [0,0], 'w');
    axis equal box on
    set(gca, 'Color', [0 0 1], 'XLim', [0 360], 'YLim', [-90 90], 'XTick',... [0: 180: 360], 'YTick', [-90: 90: 90]);
        %set(gca, 'Color', [0 0 1], 'XLim', [0 360], 'YLim', [-90 90]);
    grid; hold on
    plot(l, b, 'y.', 'linewidth', 1);
    hold off

Star Strider 2023-4-3

编辑：Star Strider 2023-4-3

在 MATLAB Online 中打开

Running it —

% l= longitude, b= latitude

l1 = -96.81

l1 = -96.8100

b1 = 32.78

b1 = 32.7800

% translating polar coordinates of the first city (Dallas) to cartesian coordinates

x1 = cosd (-96.81)*cosd (32.78)

x1 = -0.0997

y1 = sind (-96.81)*cosd (32.78)

y1 = -0.8348

z1 = sind (32.78)

z1 = 0.5414

l2 = 14.42

l2 = 14.4200

b2 = 50.07

b2 = 50.0700

% translating polar coordinates of the second city (Prague) to cartesian coordinates

x2 = cosd (14.42)*cosd (50.07)

x2 = 0.6216

y2 = sind (14.42)*cosd (50.07)

y2 = 0.1598

z2 = sind (50.07)

z2 = 0.7668

% calculating angular distance between the two cities (Dallas and Prague)

psi = acosd ((x1*x2)+(y1*y2)+(z1*z2))

psi = 77.3049

x3 = ((x2 - x1)*cosd(psi))/ sind(psi)

x3 = 0.1625

y3 = ((y2 - y1)*cosd(psi))/ sind(psi)

y3 = 0.2241

z3 = ((z2 - z1)*cosd(psi))/ sind(psi)

z3 = 0.0508

% calculating the distance between the two cities (Dallas and Prague)

r = 6371;

Distance = 2*pi*r*(psi/360)

Distance = 8.5959e+03

% Plotting the great circle path between the two cities

phi = linspace(0, psi, 100)

phi = 1×100

0 0.7809 1.5617 2.3426 3.1234 3.9043 4.6851 5.4660 6.2469 7.0277 7.8086 8.5894 9.3703 10.1512 10.9320 11.7129 12.4937 13.2746 14.0554 14.8363 15.6172 16.3980 17.1789 17.9597 18.7406 19.5214 20.3023 21.0832 21.8640 22.6449

X = x1*cosd(phi) + x3*sind(phi);

Y = y1*cosd(phi) + y3*sind(phi);

Z = z1*cosd(phi) + z3*sind(phi);

% converting cartesian coordinates x, y, z to polar coordinates l, b % (longitude and latitude)

b = asind(Z)

b = 1×100

32.7800 32.8237 32.8606 32.8907 32.9138 32.9301 32.9395 32.9420 32.9376 32.9263 32.9081 32.8831 32.8511 32.8123 32.7667 32.7142 32.6550 32.5890 32.5162 32.4367 32.3505 32.2576 32.1582 32.0522 31.9396 31.8206 31.6951 31.5632 31.4250 31.2804

l = atan2d (Y,X)

l = 1×100

-96.8100 -96.6843 -96.5576 -96.4298 -96.3008 -96.1706 -96.0392 -95.9063 -95.7721 -95.6363 -95.4990 -95.3600 -95.2193 -95.0768 -94.9323 -94.7859 -94.6374 -94.4867 -94.3337 -94.1783 -94.0204 -93.8599 -93.6966 -93.5304 -93.3612 -93.1888 -93.0131 -92.8339 -92.6511 -92.4644

l(l < 0) = l(l < 0) + 360

l = 1×100

263.1900 263.3157 263.4424 263.5702 263.6992 263.8294 263.9608 264.0937 264.2279 264.3637 264.5010 264.6400 264.7807 264.9232 265.0677 265.2141 265.3626 265.5133 265.6663 265.8217 265.9796 266.1401 266.3034 266.4696 266.6388 266.8112 266.9869 267.1661 267.3489 267.5356

load('topo.mat', 'topo', 'topomap1');

% load earth topography

whos topo topomap1;

Name Size Bytes Class Attributes topo 180x360 518400 double topomap1 64x3 1536 double

% set coordinates of points

Longitude = 0:15:360; Latitude = -72:6:72;

figure

% set background parameters

contour(0:359, -89:90, topo, [0,0], 'w');

axis equal

box on

set(gca, 'Color', [0 0 1], 'XLim', [0 360], 'YLim', [-90 90], 'XTick',...

[0: 180: 360], 'YTick', [-90: 90: 90]);

%set(gca, 'Color', [0 0 1], 'XLim', [0 360], 'YLim', [-90 90]);

grid; hold on

plot(l, b, 'y.', 'linewidth', 1);

hold off

I had to insert some linefeeds to clarify the code lines, however it runs without error. I defer to you to determine if it produces the correct result.

Star Strider 2023-4-3

在 MATLAB Online 中打开

As always, my pleasure!

Here is my analysis, and it appears to provide the correct flight path (that should provide a nice view of tha aurora borealis about midway through the flight for those on the left side of the airplane) —

% % % % % GREAT CIRCLE NAVIGATION

% % % great_circle.m

% % % 03-Mar-2023

% % % Star Strider

b1 = 32.78; % StartLat

l1 = -96.81; % StartLon

b2 = 50.07; % EndLat

l2 = 14.42; % EndLon

r = 6371; % Earth Radius (km)

[x1,y1,z1] = sp2ct(b1,l1);

[x2,y2,z2] = sp2ct(b2,l2);

psi = acosd(x1.*x2+y1.*y2+z1.*z2) % Cities Angular Distance (°)

psi = 77.3049

Distance = 2*pi*r*(psi/360);

fprintf('Distance = %.2f km\n',Distance)

Distance = 8595.92 km

phiv = linspace(0,psi,250).'; % Angular Flight Path

[xv,yv,zv] = angdist(x1,y1,z1,x2,y2,z2,phiv,psi);

[lv,bv] = ct2sp(xv,yv,zv);

lv(lv<0) = lv(lv<0) + 360;

figure

load('topo.mat', 'topo', 'topomap1');

% load earth topography

whos topo topomap1;

Name Size Bytes Class Attributes topo 180x360 518400 double topomap1 64x3 1536 double

% set coordinates of points

Longitude = 0:15:360; Latitude = -72:6:72;

% set background parameters

contour(0:359, -89:90, topo, [0,0], 'w');

axis equal

box on

set(gca, 'Color', [0 0 1], 'XLim', [0 360], 'YLim', [-90 90], 'XTick',...

[0: 180: 360], 'YTick', [-90: 90: 90]);

%set(gca, 'Color', [0 0 1], 'XLim', [0 360], 'YLim', [-90 90]);

grid; hold on

plot(lv, bv, 'y.', 'linewidth', 1);

hold off

function [lon,lat] = ct2sp(x,y,z)

lon = atan2d(y,x);

lat = asind(z);

end

function [x,y,z] = angdist(x1,y1,z1,x2,y2,z2,phi,psi)

% x90 = (x2-x1.*cosd(psi))./sind(psi);

% y90 = (y2-y1.*cosd(psi))./sind(psi);

% z90 = (z2-z1.*cosd(psi))./sind(psi);

% x = x1.*cosd(phi).*x90.*sind(phi);

% y = y1.*cosd(phi).*y90.*sind(phi);

% z = z1.*cosd(phi).*z90.*sind(phi);

x = (x1.*sind(psi-phi)+x2.*sind(phi))./sind(psi);

y = (y1.*sind(psi-phi)+y2.*sind(phi))./sind(psi);

z = (z1.*sind(psi-phi)+z2.*sind(phi))./sind(psi);

end

function [x,y,z] = sp2ct(lat,lon)

x = cosd(lon).*cosd(lat);

y = sind(lon).*cosd(lat);

z = sind(lat);

end

I created functions for the repetitive calculations, and also to simplify the code design. (This is simply a personal preference, and I do not intend it to be criticism of your code.) There are also MATLAB functions cart2sph and sph2cart that will do the conversions. I coded the versions of them here that correspond to the provided great circle documentation. I will let you sort this to determine where the errors may be in your code.

For my part, I now have a great circle .m file if I want to use it in the future!

Star Strider 2023-4-4

‘Okay,its mainly following the guide given by the calculations guide.’

Pretty much completely following it.

‘Could you just explain the the 'angdist' function and the 'function' subplot at the end?’

I thought I already explained ‘angdist’ as being

applied to each ‘(x1,x2)’, ‘(y1,y2)’, and ‘(z1,z2)’ pair in turn. It is the same relation applied to each pair. It returns the Cartesian coordinates ‘xv’, ‘yv’, and ‘zv’ that are vectors because φ is the vector ‘phiv’. Those are then transformed into the longitude coordinates ‘lv’, and corresponding latitude coordinates ‘bv’, in the ‘ct2sp’ call immediately after that, both of them being vectors as well, because their arguments are vectors. The ‘lv’ vector is further processed by adding 360 to its negative elements. Those are then plotted on the map.

There is no function subplot. The code after the figure call I copied from your existing code because it loads the map and sets its parameters. (I cannot improve on that, and besides it is consistent with the desired result with respect to displaying the flight path.) The only difference between my code and yours is:

plot(lv, bv, 'y.', 'linewidth', 1);

The three functions at the end of my script are the functions I use in the code.

I just now noticed that you did not post the MATLAB version/release that you are using in either of your threads. My code should run without error in all recent releases. If you are having problems with it because of version/release compatibility issues, I might be able to help with that. The online Run feature here and my offline computers all run R2023a.

Light 2023-4-12

编辑：Walter Roberson 2023-4-21

在 MATLAB Online 中打开

PRGA2001 Drag Coefficient.pdf

Okay no problem, thanks very much. I know I didnt declare the differential equations properly, but this is all I have right now:

%Differential Equations
% m (d2x / dt^2) = -k|v| (dx/dt)
% m (d2y / dt^2) = -mg - k|v| (dy/dt) 
% |v| = sqrt ( (dx/ dt)^2 + (dy/ dt)^2 )
= v0*cos,  = v0*sin
% k = air resistance coefficient
% 𝜃 (theta) = initial angle of elevation of the ball when it was kicked
% x and y are respectively the horizontal and vertical spatial coordinates
% t is the time measured from the moment the ball is kicked
% |v| = mag_v = is the magnitude of the velocity (speed, changes with time)
% theta = ?
m = 0.45 % in kg
g = 9.81 % gravitational acceleration in m/s^2
v0 = 108 % in km/h
function dBalldt = Football_F(t, Ball)
m = 0.45 % in kg
g = 9.81 % gravitational acceleration in m/s^2
v0 = 108 % in km/h
% theta = ?
ICs when t = 0
x = 0, y = 0
mag_v = sqrt((v0*cos*theta).^2 + (v0*sin*theta).^2)
% dxdt = @ (t,x)[-k*(mag_v*v0*cos*theta); (-m*g)*-k*(mag_v*v0*sin*theta)];
% m_dxdt = (m*dxdt);
dx/dt = Vx
dy/dt = Vy
dVx= -k*(mag_v*v0*cos*theta)
dVy = (-m*g)*-k*(mag_v*v0*sin*theta)
m_dVx = m*(-k*(mag_v*v0*cos*theta))
m_dVy = m*(-m*g)*-k*(mag_v*v0*sin*theta)
[t,x] = ode45 (m_dxdt, [0,20], [0,0,0]);
return 

Star Strider 2023-4-13

在 MATLAB Online 中打开

Today turned out to be an unscheduled ‘errand day’ so I was away for several hours.

If I am correct, the objective here is to develop a set of differnetial equations and then estimate parameters using the supplied data. If that is the situation, it would be relatively straightforward to code that. See Parameter Estimation for a System of Differential Equations for an example.

Extracting ‘VAR DATA’ to a table and saving it to a .mat file is straightforward.

VAR_DATA = array2table([0 108

0.5 107

1.1 106

1.7 105

2.3 104

2.9 104

3.5 103

4.0 102

4.6 101

5.2 100

5.7 99

6.3 99

6.8 98

7.4 97

7.9 97

8.4 96

9.0 95

9.5 94

10.0 94

10.6 93

11.1 92

11.6 92

12.1 91

12.6 91

13.1 90

13.6 89

14.1 89

14.6 88

15.1 88

15.6 87

16.1 87

16.5 86

17.0 86

17.5 85

18.0 85

18.4 84

18.9 83

19.4 82], 'VariableNames',{'Distance (m)','Speed (kph)'})

VAR_DATA = 38×2 table

Distance (m) Speed (kph) ____________ ___________ 0 108 0.5 107 1.1 106 1.7 105 2.3 104 2.9 104 3.5 103 4 102 4.6 101 5.2 100 5.7 99 6.3 99 6.8 98 7.4 97 7.9 97 8.4 96

save('VAR_DATA.mat','VAR_DATA')

L = size(VAR_DATA,1);

t = linspace(0, L-1, L);

figure

stem3(t, VAR_DATA{:,1}, VAR_DATA{:,2}, 'filled')

grid on

xlabel('t')

ylabel('Distance')

zlabel('Velocity')

axis('equal')

figure

plot(VAR_DATA{:,1}, VAR_DATA{:,2}, '.-')

grid on

xlabel('Distance')

ylabel('Velocity')

axis('equal')

I will address this tomorrow. In the interim, using this information and my prototype code, experiment with coding it. Ideally, it would be preferable if there was also a time vector. in addition to distance and velocity. It it possible to get that, or infer it from the data? For example, are the VAR data sampled at a constant sampling interval or frequency? (I don’t know, and I’ve never encountered this previously, so I’ve never explored it. I just enjoy watching the matches.)

However, it might be possible to use the Symbolic Math Toolbox to solve the differential equations in a closed form , with that solution being the objective function for a nonlinear parameter estimation problem. Try that as well.

Light 2023-4-19

编辑：Walter Roberson 2023-4-21

在 MATLAB Online 中打开

Hi good day Star, I have been going over the code and I have been having some trouble with it. I'm getting errors which im not sure how to resolve, any idea how the code could be fixed? Im really trying my best

Differential Equations

%Differential Equations
% m (d2x / dt^2) = -k|v| (dx/dt)
 % m (d2y / dt^2) = -mg - k|v| (dy/dt) 
% |v| = sqrt ( (dx/ dt)^2 + (dy/ dt)^2 )
% |v| = mag_v = sqrt((v0*cos(theta).^2 + (v0*sin(theta).^2)))
% k = air resistance coefficient
% 𝜃 (theta) = initial angle of elevation of the ball when it was kicked
% x and y are respectively the horizontal and vertical spatial coordinates
% t is the time measured from the moment the ball is kicked
% |v| = mag_v = is the magnitude of the velocity (speed, changes with time)
global g
%function dBalldt = Football_F(t, Ball)
%Football_F (1,[1 2 3 4])
m = 0.45 % in kg
g = 9.81 % gravitational acceleration in m/s^2
v0 = 108 % in km/h
mag_v = sqrt((v0*cos(theta).^2 + (v0*sin(theta).^2)))
t_range = 0: 0.1: 2;
theta = 30*pi/180;
% v = 10
% k = ?
Vx = v0*cos(theta)
Vy = v0*sin(theta)
ICs = [0 0 v0*cos(theta) v0*sin(theta)]
% when t = 0
% x = 0, y = 0
% dxdt = @ (t,x)[-k*(mag_v*v0*cos*theta); (-m*g)*-k*(mag_v*v0*sin*theta)];
% m_dxdt = (m*dxdt);
dxdt = Vx
dydt = Vy
%dVx= -k*(mag_v*v0*cos*theta)
%dVy = (-m*g)*-k*(mag_v*v0*sin*theta)
m_dVx = m*(-k*(mag_v*v0*cos*theta))
m_dVy = m*(-m*g)*-k*(mag_v*v0*sin*theta)
[t,xyVxVy] = ode45 (Football_F, t_range, ICs);
%extra arrays
x = xyVxVy (:, 1);
y = xyVxVy (:, 2);
Vx = xyVxVy (:, 3);
Vy = xyVxVy (:, 4);
plot (x,y)
grid
[x y]
y_x5_theoretical = interp1 (x, y, 5);

Star Strider 2023-4-19

在 MATLAB Online 中打开

VAR_DATA.mat

I created and uploaded ‘VAR_DATA.mat’ to make this a bit easier.

I am still somewhat lost on this, and have not given it much thought.

syms k t x(t) y(t) Y T

m = sym(450); % g

g = sym(9.81); % m/s^2

v0 = sym(108); % km/h

v = sqrt(diff(x)^2+diff(y)^2);

Eq1 = m*diff(x,2) == -k*v*diff(x);

Eq2 = m*diff(y,2) == -m*g - k*v*diff(y);

[VF,Subs] = odeToVectorField(Eq1,Eq2)

VF =

Subs =

deqfcn = matlabFunction(VF,'Vars',{T,Y,k})

deqfcn = function_handle with value:

@(T,Y,k)[Y(2);k.*sqrt(Y(2).^2+Y(4).^2).*Y(2).*(-1.0./4.5e+2)-9.81e+2./1.0e+2;Y(4);k.*sqrt(Y(2).^2+Y(4).^2).*Y(4).*(-1.0./4.5e+2)]

load('VAR_DATA.mat')

VAR_DATA

VAR_DATA = 38×2 table

Distance (m) Speed (kph) ____________ ___________ 0 108 0.5 107 1.1 106 1.7 105 2.3 104 2.9 104 3.5 103 4 102 4.6 101 5.2 100 5.7 99 6.3 99 6.8 98 7.4 97 7.9 97 8.4 96

VAR_MTX = table2array(VAR_DATA);

t = linspace(0, 3, size(VAR_MTX,1)).';

B0 = [1; pi/2; VAR_MTX(1,2)];

B = lsqcurvefit(@objfcn, B0, t, VAR_MTX)

Local minimum possible. lsqcurvefit stopped because the final change in the sum of squares relative to its initial value is less than the value of the function tolerance.

B = 3×1

50.5159 1.4735 103.6831

xy = objfcn(B,t)

xy = 38×2

10.0711 103.1928 10.0711 103.1606 10.0711 103.0644 10.0711 102.9057 10.0711 102.6865 10.0711 102.4099 10.0711 102.0791 10.0711 101.6982 10.0711 101.2711 10.0711 100.8022

figure

plot3(VAR_MTX(:,1), VAR_MTX(:,2), t, '.')

hold on

plot3(xy(:,1), xy(:,2), t, '-', 'LineWidth',1.5)

grid

xlabel('Distance')

ylabel('Velocity')

zlabel('t')

legend('Data','Fit', 'Location','best')

function xy = objfcn(b,t)

k = b(1);

theta = b(2);

v0 = b(3);

deqfcn = @(T,Y,k)[Y(2);k.*sqrt(Y(2).^2+Y(4).^2).*Y(2).*(-1.0./4.5e+2)-9.81e+2./1.0e+2;Y(4);k.*sqrt(Y(2).^2+Y(4).^2).*Y(4).*(-1.0./4.5e+2)];

ic = [v0*sin(theta) 0 v0*cos(theta) 0];

[t,y] = ode45(@(t,y)deqfcn(t,y,k), t, ic);

xy = y(:,[3 1]);

end

This is how I usually solve these sorts of problems, however I genuinely do not have a good feeling about what the variables are here. We do not have a time vector (although the exact time vector would be helpful), and I am not certain how to model this. I used the Symbolic Math Toolbox to derive the code for the differential equaations simply because it avoids the algebraic errors that I usually make (and that can be frustrating to detect and correct). Experiment with this code to estimate the parameters and correctly determine what the data are that are being modeled.

Note that you can select the variables that ‘objfcn’ returns by defining what columns the ‘xy’ variable selects from the four columns of ‘y’ returned by the ode45 call. These values will be used by lsqcurvefit to fit the data. The ‘names’ of the variables are in the ‘Subs’ output of the odeToVectorField call.

The earlier problem was straightforward aand relatively easy to solve. This one is not.

Light 2023-4-20

编辑：Light 2023-4-20

在 MATLAB Online 中打开

PRGA2001 Drag Coefficient (2).pdf

Do you know how steps 1-6 in the "Steps'' section of the problem is to be carried out, using trapz or quad for number 3?

This is what I have thus far:

% k = air resistance coefficient
% 𝜃 (theta) = initial angle of elevation of the ball when it was kicked
% x and y are respectively the horizontal and vertical spatial coordinates
% t is the time measured from the moment the ball is kicked
% |v| = mag_v = is the magnitude of the velocity (speed, changes with time)
global g
%function dBalldt = Football_F(t, Ball)
%Football_F (1,[1 2 3 4])
m = 0.45 % in kg
g = 9.81 % gravitational acceleration in m/s^2
v0 = 108 % in km/h
mag_v = sqrt((v0*cos(theta).^2 + (v0*sin(theta).^2)))
t = 0: 0.1: 2;
theta = 30*pi/180;  %initial guess for theta
k = 3               %initial guess for air resistance in N
% v = 10
Vx = v0*cos(theta)
Vy = v0*sin(theta)
ICs = [0 0 v0*cos(theta) v0*sin(theta)]
% when t = 0
% x = 0, y = 0
% dxdt = @ (t,x)[-k*(mag_v*v0*cos*theta); (-m*g)*-k*(mag_v*v0*sin*theta)];
% m_dxdt = (m*dxdt);
dxdt = Vx
dydt = Vy
%dVx= -k*(mag_v*v0*cos*theta)
%dVy = (-m*g)*-k*(mag_v*v0*sin*theta)
m_dVx = m*(-k*(mag_v*v0*cos*theta))
m_dVy = m*(-m*g)*-k*(mag_v*v0*sin*theta)
[t,xyVxVy] = ode45 (Football_F, t, ICs);
%extra arrays
x = xyVxVy (:, 1);
y = xyVxVy (:, 2);
Vx = xyVxVy (:, 3);
Vy = xyVxVy (:, 4);
plot (x,y)
grid
[x y]
y_x5_theoretical = interp1 (x, y, 5)