## ENGG1811 Lab 9: Simulation and its applications, and Matlab functions

### Objectives

After completing this lab, students should be able to
• Write Matlab scripts and functions
• Apply the Matlab programming concepts of for-loop and matrix indexing
• Use built-in functions find and max
• Gain a better understanding of simulation and how it can be used for design and decision making

### Assessment

This lab has four parts: Parts A to D. The four parts together form a sequence of tasks. These parts are not independent (like the previous labs) and should be done one after another. You need to show your tutors all the four parts. All scripts and functions should be reasonably documented.

The on-line assessment is about for-loop, and can be done at any time.

You should make a directory called lab09 to store your files for this lab.

### Introduction:

Imagine a mobile phone game where you use a strawberry launcher to feed an alien. Your goal is to feed as many strawberries to the alien in as little time as possible. Technically speaking, the problem to solve is, for a given location of the alien, what launch speed and launch angle will minimise the travel time of the strawberries from the launcher to the alien.

Remark: The real-life application of this problem is missiles. For this lab, let us be innocent and consider strawberries and aliens instead.

The starting point of simulation is mathematical modelling. In ENGG1811, we will always give you the mathematical models. After you have learned enough knowledge in your own discipline, you will be able to come out with your own mathematical model for the work that you do.

The motion of the strawberries in the air is affected by two factors: gravitational pull and air resistance. A common model for air resistance is that the resistive force is proportional to total speed of the object. We assume the launcher, strawberries and alien live in a two-dimension space. The horizontal direction is the x-direction and the vertical direction is the y-direction. The position of the strawberry is therefore defined by its x- and y-coordinates. We use px(t) and py(t) to denote the x- and y-coordinates of the strawberry at time t. In addition, we use vx(t) and vy(t) to denote the velocity of the strawberry in x- and y-directions. The motion of the strawberry obeys the following set of ordinary differential equations:

where m and c are respectively the mass and drag coefficient, and g is the acceleration due to gravity. Since the strawberry has an x- and a y-component in velocity, its total speed is given by √(vx(t)2 + vy(t)2). If you are interested, Equations (3) and (4) are derived from Newton's Second Law. Equation (3) says that the strawberry experiences resistance as it moves in the air in the x-direction. Equation (4) says that the strawberry experiences both gravitational pull and air resistance when it moves in the y-direction.

Given Equations (1)-(4), you can use Euler's forward rule to derive how to calculate the position and velocity at time t + Δ from that at time t. We have derived the equations for you, as follows:

Let us assume that you use the Matlab vector variables vecPositionX,  vecPositionY, vecVelocityX and vecVelocityY to store the coordinates and velocities of the strawberry at different time instances. Equations (5) and (7) mean that you will be updating the elements of the vectors vecPositionX and vecVelocityX using these Matlab statements:

vecPositionX(k) = vecPositionX(k-1) + vecVelocityX(k-1) * dt;
speedTotal = sqrt(vecVelocityX(k-1)^2+vecVelocityY(k-1)^2);
vecVelocityX(k) = vecVelocityX(k-1) - c*speedTotal*vecVelocityX(k-1)/m*dt;

where dt in the Matlab statements is the time increment, which is denoted by Δ in the mathematical equations. Essentially you identify px(t) with the vecPositionX(k-1), px(t+ Δ) with the vecPositionX(k) etc.

You are asked to work out the other Matlab statements for updating the elements of the vectors vecPositionY and vecVelocityY.

### Part A: Simulation

Download the zip file which contains the files simulateProjectile_PartiallyDone.m and testSimulateProjectile.m, and move them to the directory lab09 that you have created earlier. The file simulateProjectile_PartiallyDone.m contains a Matlab function to simulate the motion of the strawberry. First, rename the file as simulateProjectile.m. This is so that if you, for some reason, need to download the file again, the newly downloaded file will not overwrite the work that you have.

The function in simulateProjectile.m has the following form:

[vecPositionX,vecPositionY] = simulateProjectile(vecTime,m,c,v0,theta0d)

where vecTime is a Matlab vector containing time instances, m is mass, c is drag coefficient, v0 is initial launch speed, theta0d is launch angle in degrees. This function returns two outputs, which are the x- and y-coordinates at the corresponding time instances specified by the vector vecTime. The function is missing two lines of code (within the for-loop located at the end of the file), which are the two Matlab statements for updating positionY and velocityY. You are asked to fill in these two statements at the places indicated. You may wish to read the function to see what it does. In particular, note that v0 is the initial total speed, which means the initial speed in the x- and y-directions, are, respectively v0*cosd(theta0d) and v0*sind(theta0d) where cosd() and sind() are Matlab trigonometric functions which take angles expressed in degrees.

After you have completed this, you can test the function by using the Matlab script testSimulateProjectile.m. Take a look at what the script does. It defines the problem and simulation parameters, as follows:
• time starts at 0 and ends at 10 with an increment of 0.01
• m = 0.145
• c = 0.0013
• v0 = 50
We will not change these parameters for this lab.

For this part, you can use theta0d = 35. We will use a range of values for theta0d later on.

If you've added the correct Matlab statements, then a plot of the x- and y-coordinates of the strawberry should look like this. (Remark: In Matlab, you can generate picture files for graphs using the command print. See the help page on how to use it.)

### Part B: Determine the landing distance and landing time

The simulation in Part A allows you to determine the trajectory of the strawberry launched at a given speed and angle. Let us now assume that the alien is located at a height (y-coordinate) of +10, and the alien's can only swallow the strawberries when they are falling down. The trajectory plot in Part A shows that, with a launch angle of 35 degrees, a rough estimate of the landing position of the strawberries will be at the coordinates (96,10). This position is obtained from using visual inspection from the graph. Your task for this Part is to use Matlab to automatically determine the landing position (assuming a landing height of 10) as well as the time to reach the landing position, from the trajectories computed in Part A.We will define the exact meaning of landing using an example in a moment.

We will give you some idea how you can use Matlab to find the landing time and position. The trajectory has 501 numbers per vector so they may be hard to work with. So we will use a small example to illustrate the method. This is a typical problem solving technique. If you can solve the problem for a small example, then you can transfer the solution to a larger instance of the problem. Let us assume that you have three vectors px, py and tv, which have analogous meaning to vecPositionX, vecPositionY and vecTime that we have used above. We made up some values for the elements of these three vectors:

tv =      [   0.6    0.7        1.5     1.6     1.7    1.8     1.9]

px =     [   2.4    2.9       5.7     6.6      7.5    8.4    9.3]

py =      [  9.7    10.1   10.5   10.2   10.1   9.9    9.8]

Note that this example contains the main feature of the problem with the positions in the y-direction dropping from a value above 10 to a value below 10. It means that at some point in time, the strawberry has a height of exactly 10. Instead of determining the exact time of crossing, we use an approximation here. First, we consider only the time instances in the vector tv. Second, we say that landing occurs at time t if the vertical position of the strawberry at time t is closest to the level of 10, but is still at or above the level of 10. For this example, landing occurs at time instant 1.7 because the strawberry is just above the level 10, at a height of 10.1 to be specific. At the next time instance, at time 1.8, the strawberry is below the level of 10, at 9.9 to be specific. This is not exact, but if the time increment is small, the error will not be big. (This is also on the condition that the alien opens its mouth wide!) We say that the landing distance (which measures position in the x-direction) is 7.5 and the landing time is 1.7.

We suggest that you work with this example first to find out what Matlab statements you need to determine the landing distance and landing time from tv, px and py. Once you get this going, you can test whether it works for the vecTime, vecPositionX and vecPositionY vectors.

You can use any method to solve this problem as long as it can be automated. Having done some programming, you may think that you need some loops and if statements to solve this problem. Yes, you can do that, but I want to say that this problem has some special features and you can solve this problem using the Matlab built-in functions find and max. If you do not know how to use these functions, you can click on the help button (white question mark on a blue circle) to look at the help pages of these functions. If you want some more hint, you can click here.

You should try to do this Part using a script, because you will need to write a function to do these steps in Part C. You can keep the script to show your tutor.

### Part C: Write a function to determine the landing distance and landing time

Once you have got Part B working, you should write a function to determine the landing distance and landing time (these are the outputs) from the Matlab vectors vecTime, vecPositionX and vecPositionY and a scalar landingLevel (these are the inputs). Note that landingLevel is 10 but we should try to make the function general so that it can work with any landing level.

The beginning of your function should look like this:

function [landingPosition,landingTime] = ...
findLandingPositionAndTime(vecTime,vecPositionX,vecPositionY,landingLevel)

### Part D: Determining the launch angle to feed the alien in a minimum time

First, let us recap what you have achieved by now.

• For a given launch angle, the function simulateProjectile gives you vecPositionX and vecPositionY (Part A)
• From vecTime, vecPositionX and vecPositionY, the function findLandingPositionAndTime gives you landingPosition and landingTime (Part C)

By combining what you have in Parts A and C, you can determine, for a given launch angle, what the landing position and landing time are. What this means is that, you have all the functions that you need to study how the launch angle impacts on landing position and landing time.

In this part, your task is to write a Matlab script that does the following:

• Create a vector that starts with 20 and ends at 60 at an increment of 1. This vector contains all the launch angles that you want to use. We will call this vector vecAngles but you should come out with a better variable name.
• Create two zero vectors that have the same size as the vector vecAngles. These vectors will be used to store the landing distances and landing times for different launch angles.
• For each of the launch angles, you want to:
• Simulate to obtain the trajectory. You can use the same parameters as Part A.
• Determine the landing distance and landing time for the trajectory assuming a landing height of 10
• Store the calculated landing distance and landing time in vectors that you have created

This script does some calculations repeatedly so you need some type of loops (which type?). Once you get your script working, you can plot a graph of landing distance versus landing time. You should get a graph similar to that below.

Let us assume that the alien is at a horizontal distance of 90 from the launcher. You can see from the graph that there are two launch angles that will allow you to feed the alien and one of them will allow you to feed the alien quicker. Your final task is to determine the launch angle to use by using Matlab statements. You can add these statements to the script and write a Matlab disp statement to display the answer. We do not need a very accurate answer. As long as the launch angle is within 1 degree of the correct value, that is acceptable. You can use any method you like to determine the launch angle. The type of method you used in Part B is also helpful here.

--- Happy feeding! ---