Lecture 2 - Logic, Loops, and Arrays

This iPython notebook covers some of the most important aspects of the Python language that is used daily by real Astronomers and Physicists. Topics will include:

1. The logic of Python, including while loops and if/else statements
2. Function definitions and how to make your own
3. A review of numpy arrays and a discussion of their usefulness in solving real problems
4. Reading in data from text and numpy file formats, along with creating your own outputs to be used later

A. Logic, If/Else, and Loops

You can make conditional (logical) statements in Python, which return either "True" or "False", also known as "Booleans."

A basic logic statement is something that we've used already: x < y. Here is this one again, and a few more.

If you let a and b be conditional statements (like the above statements, e.g. a = x < y), then you can combine the two together using logical operators, which can be thought of as functions for conditional statements.

There are three logical operators that are handy to know:

• And operator: a and b
• Or operator: a or b
• Not operator: not(a)

Now, these might not seem especially useful at first, but they're the bread and butter of programming. Even more importantly, they are used when we are doing if/else statements or loops, which we will now cover.

An if/else statement (or simply an if statement) are segments of code that have a conditional statement built into it, such that the code within that segment doesn't activate unless the conditional statement is true.

Here's an example. Play around with the variables x and y to see what happens.

The idea here is that Python checks to see if the statement (in this case "x < y") is True. If it is, then it will do what is below the if statement. The else statement tells Python what to do if the condition is False.

Note that Python requires you to indent these segments of code, and WILL NOT like it if you don't. Some languages don't require it, but Python is very particular when it comes to this point. (The parentheses around the conditional statement, however, are optional.)

You also do not need an "else" segment, which effectively means that if the condition isn't True, then that segment of code doesn't do anything, and Python will just continue on past the if statement.

Here is an example of such a case. Play around with it to see what happens when you change the values of x and y.

While-loops are similar to if statements, in the sense that they also have a conditional statement that is built into it and it executes when the conditional is True. However, the only difference is, it will KEEP executing that segment of code until the conditional statement becomes False.

This might seem a bit strange, but you can get the hang of it!

For example, let's say we want Python to count from 1 to 10.

Note here that we tell Python to print the number x (x starts at 1) and then redefining x as itself +1 (so, x=1 gets redefined to x = x+1 = 1+1 = 2). Python then executes the loop again, but now x has been incremented by 1. We continue this process from x = 1 to x = 10, printing out x every time. Thus, with a fairly compact bit of code, you get 10 lines of output.

It is sometimes handy to define what is known as a DUMMY VARIABLE, whose only job is to count the number of times the loop has been executed. Let's call this dummy variable i.

But...what if the conditional statement is always true?

B. Defining Your Own Functions

So far, we have really focused on using built-in functions (such as from numpy and matplotlib), but what about defining our own? This is easy to do, and can be a way to not only clean up your code, but also allows you to apply the same set of operations to multiple variables without having to explicitly write it out every time.

For example, let's say we want to define a function that takes the square root of a number. It's probably a good idea to check if the number is positive first, otherwise we'll end up with an imaginary answer.

So the outline for a function is

def <function name> (<input variable>):
    <some code here>
    return <output variable>

In general, many common mathematical functions like sqrt, log, exp, sin, cos can be found in the math module. So we don't have to write our own - phew!

When defining your own functions, you can also use multiple input variables. For example, if we want to find the greatest common divisor (gcd) of two numbers, we could apply something called Euclid's algorithm.

We define gcd(a,0) = a. Then we note that the gcd(a,b) = gcd(b,r), where r is the remainderwhen a is divided by b. So we can repeat this process until we end up with zero remainder. Then, we return whatever number is left in a as the greatest common divisor.

A command you might not have encountered yet is %. The expression x % y returns the remainder when x is divided by y.

Challenge 1 - Fibonacci Sequence and the Golden Ratio

The Fibonacci sequence is defined by $f_{n+1} =f_n + f_{n-1}$ and the initial values $f_0 = 0$ and $f_1 = 1$.

The first few elements of the sequence are: $0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377 ...$

Using what you just learned about functions, define a function fib, which calculates the $n$-th element in the Fibonacci sequence. It should have one input variable, $n$ and its return value should be the $n$-th element in the Fibonacci sequence.

The ratio of successive elements in the Fibonacci sequence converges to $$\phi = (1 + \sqrt{5})/2 ≈ 1.61803\dots$$ which is the famous golden ratio. Your task is to approximate $\phi$.

Define a function phi_approx that calculates the approximate value of $\phi$ obtained by the ratio of the $n$-th and $(n−1)$-st elements, $$f_n /f_{n-1} \approx \phi$$

phi_approx should have one variable, $n$. Its return value should be the $n$-th order approximation of $\phi$.

C. Numpy Arrays - Review of Basics and Some More Advanced Topics

Recall in the first lecture that we introduced a python module known as numpy and type of variable known as a numpy array. For review, we will call numpy to be imported into this notebook so we can use its contents.

Here we are calling in the contents of numpy and giving it the shorthand name 'np' for convenience.

To create an array variable (let's call it 'x'), we simply assign 'x' to be equal to the output of the np.array() function, using a list as an input. You can then verify its contents by using the print() function.

As we learned in Lecture 1, numpy arrays are convenient because they allow us to do math across the whole array and not just individual numbers.

For example, let's say we want to make a new variable 'y' such that y = $x^{2}$, then this is done simply as

The documentation of possible functions that can be applied to integers and floats (i.e. single numbers), as well as numpy arrays, can be found here: https://docs.scipy.org/doc/numpy/reference/routines.math.html

As discussed previously, there are numerous ways to create arrays beyond np.numpy(). These include:

• np.arange()
• np.linspace()

These create arrays of numbers within a range with a specific step-size between each consecutive number in the array.

It is sometimes convenient to have Python create other arrays for you, depending on the problem that you are going to solve. For example, sometimes it is handy to create an array of all zeros, which can then be replaced later with data. This can be done by using np.zeros().

Going back to the Fibonacci example, let's say we want to store the first 10 elements of the Fibonacci sequence in an array for easy access in the future. To ready such an array, you simply do the following.

Now how do we assign a new value to an element of the array? We use the following "square bracket" notation:

array_name[index_number] = value

In this, the array (with the name "array_name" or whatever it is you have named it) will have "value" replace whatever is in the position corresponding to "index_number."

Arrays are numbered starting from 0, such that

• First position = 0
• Second position = 1
• Third position = 2
• etc.

It is a bit confusing, but after a bit of time, this becomes quite natural. Let's practice with the Fibonacci example.

First, let's store the first Fibonacci number in our array. We use the brackets to store that value in the first position (0 index number) in the data array we made above.

Now you try it. Store the second Fibonacci number in the second position of your array and use a print statement to verify that you have done so.

Python array indexing is fairly straightforward once you get the hang of it.

Let's say you wanted the last element of the array, but you don't quite recall the size of the array. One of the easiest ways to access that element is to use negative indexing.

Negative indexing is the same as normal indexing, but backward, in the sense that you start with the last element of the array and count forward. More explicitly, for any array:

• array[-1] = last element of array
• array[-2] = second to last element of the array
• array[-3] = third to last element of the array
• etc

Now then, let's create an array using np.arange() with 10 elements, and see if you can access the last element and the second to last element using negative indexing. Print out these values.

Now, sometimes its useful to access more than one element of an array. Let's say that we have an array with 100 elements in the range [0,10] (including endpoints). If you recall, this can be done via the np.linspace() function.

Now then, in order to get a range of elements rather than simply a single one, we use the notation:

x[i_start,i_end+1]

For example, let's say you want the 1st, 2nd, and 3rd element, then you'd have to do

x[0:3]

In this notation, ":" represents you want everything between 0 and 3, and including 0. Let's test this.

If you want everything passed a certain point of the array (including that point), then you would just eliminate the right number, for example

x[90:]

would give you everything after (and including) the 90 index element. Similarly, if you want everything before a certain index

x[:90]

would give you everything before the 90 index element.

So, let's say that you would want everything up to and including the tenth element of the array $x$. How would you do that?

(Remember, the tenth element has an index of 9)

Finally, simply using the ":" gives you all the elements in the array.

Challenge 2 - Projectile Motion

In this challenge problem, you will be building what is known as a NUMERICAL INTEGRATOR in order to predict the projectiles trajectory through a gravitational field (i.e. what happens when you throw a ball through the air)

Let's say that you have a projectile (let's say a ball) in a world with 2 spatial dimensions (dimensions x and y). This world has a constant acceleration due to gravity (call it simply g) that points in the -y direction and has a surface at y = 0.

Can we calculate the motion of the projectile in the x-y plane after the projectile is given some initial velocity vector v? In particular, can we predict where the ball will land? With loops, yes we can!

Let's first define all of the relevant variables so far. Let g = -9.8 (units of m/s, so an Earth-like world), the initial velocity vector being an array v = [3.,3.], and an initial position vector (call it r) in the x-y plane of r = [0.,1.]. For ease, let's use numpy arrays for the vectors.

Now then, remember that,

$a = \frac{dv}{dt}$ and thus, $dv = a\ dt$

So, the change of an objects velocity ($dv$) is equal to the acceleration ($a = g$ in this case) multiplied by the change in time ($dt$)

Likewise:

$v_{x} = \frac{dx}{dt}$ and $v_{y} = \frac{dy}{dt}$, or

$v_{x}\ dt = dx$ and $v_{y}\ dt = dy$

Now, in this case, since there is only downward acceleration, the change of $v_{x}$ is 0 until the projective hits the ground.

Now, we're going to define two functions, one that will calculate the velocity vector components and the other, the position vector components, and returning a new vector with the new components. I'll give you the first one.

Now that we've defined intV (short for "integrate v"), let's use it real quick, just to test it out. Let dt = 0.1 (meaning, your taking a step forward in time by 0.1 seconds).

As you can see, $V_{x}$ hasn't changed, but $V_{y}$ has decreased, representing the projectile slowing down as it's going upward.

I'll let you define the function now for the position vector. Call it intR, and it should be a function of (r,v,dt), and remember that now both $r_{x}$ and $r_{y}$ are changing. Remember to return an array.

Now we have the functions that calculate the changes in the position and velocity vectors. We're almost there!

Now, we will need a while-loop in order to step the projectile through its trajectory. What would the condition be?

Well, we know that the projectile stops when it hits the ground. So, one way we can do this is to have the condition being (r[1] >= 0), since the ground is defined at y = 0.

So, having your intV and intR functions, along with a while-loop and a dt = 0.1 (known as the "step-size"), can you use Python to predict where the projectile will land?

Now, note that we've defined the while-loop such that it doesn't stop exactly at 0. Firstly, this was strategic, since the initial y = 0, and thus the while-loop wouldn't initialize to begin with (you can try to change it). One way you can overcome this issue is to decrease dt, meaning that you're letting less time pass between each step. Ideally, you'd want dt to be infinitely small, but we don't have that convenience in reality. Re-run the cells, but with dt = 0.01 and we will get much closer to the correct answer.

So, we know where it's going to land...can we plot the full trajectory? Yes, but this is a bit complicated, and requires one last function: np.append(). https://docs.scipy.org/doc/numpy/reference/generated/numpy.append.html

The idea is to use np.append() to make an array that keeps track of where the ball has been.

Let's define two empty arrays that will store our information (x and y). This is an odd idea, defining an array variable without any elements, so instead think of it as a basket without anything inside of it yet, and we will np.append() to fill it.

Now, all you have to do, is each time the while-loop executes, you use np.append() for the x and y arrays, adding the new values to the end of them.

How do you do that? Well, looking at the np.append() documentation, for x, you do x = np.append(x,[r[0]])

The same syntax is used for the y array.

After that, you simply use plt.plot(x,y,'o') to plot the trajectory of the ball (the last 'o' is used to change the plotting from a line to points).

Good luck! Also, don't forget to reset your v and r arrays (otherwise, this will not work)

If everything turns out alright, you should get the characteristic parabola.

Also, if you're going to experiment with changing the intial position and velocity, remember to re-run the cell where we define the x and y arrays in order to clear the plot.

Now you've learned how to do numerical integration. This technique is used all throughout Physics and Astronomy, and while there are more advanced ways to do it in order to increase accuracy, the heart of the idea is the same.

Here is a figure made by Corbin Taylor (head of the Python team) that used numerical integration to track the position of a ray of light as it falls into a spinning black hole.

D. Loading And Saving Data Arrays

So, we have learned a lot about data arrays and how we can manipulate them, either through mathematics or indexing. However, up until this point, all we've done is use arrays that we ourselves created. But what happens if we have data from elsewhere? Can Python use that?

The answer is of course yes, and while there are ways to import data that are a bit complicated at times, we're going to teach you some of the most basic, and most useful, ways.

For this section, we will be using plotting to visualize the data as you're working with it, and as such, we will be loading in the package "matplotlib.pyplot" which you used in Lecture 1

Now then, let's say we are doing a timing experiment, where we look at the brightness of an object as a function of time. This is actually a very common type of measurement that you may do in research, such as looking for dips in the brightness of stars as a way to detect planets.

This data is stored in a text file named timeseries_data.txt in the directory lecture2_data. Let's load it in.

Now we have the data loaded into Python as a numpy array, and one handy thing you can do is to use Python to find the dimensions of the array. This is done by using ".shape" as so.

In this format, we know that this is a 2x1000 array (two rows, 1000 columns). Another way you can think about this is that you have two 1000-element arrays contained within another array, where each of those arrays are elements (think of it as an array of arrays).

The first row is the time stamps when each measurement was taken, while the second row is that of the value of the measurement itself.

For ease of handling this data, one can in principle take each of these rows and create new arrays out of them. Let's do just that.

Here, you have 2 dimensions with the array timeseriesData, and as such much specify the row first and then the column. So,

• array_name[n,:] is the n-th row, and all columns within that row.
• array_name[:,n] is the n-th column, and all rows within that particular column.

Now then, let's see what the data looks like using the plot() function that you learned last time. Do you remember how to do it? Why don't you try! Plot t as your x-axis and signal as your y-axis. Don't forget to show your plot.

Looking at our data, you see clear spikes that jump well above most of the signal. (I've added this to the data to represent outliers that may sometimes appear when you're messing with raw data, and those must be dealt with). In astronomy, you sometimes have relativistic charged particles, not from your source, that hit the detector known as cosmic rays, and we often have to remove these.

There are some very complex codes that handle cosmic rays, but for our purposes (keeping it easy), we're going to just set a hard cut off of, let's say 15.

In order to do this, we can use conditional indexing in place of normal indices. This involves taking a conditional statement (more on those later) and testing whether it evaluates to True on each element in the array.

This gives an array of Booleans, which we can use as logical indices to select only the entries for which the logical statement is True.

In this case, the conditional statement that we have used is signal < cutOff.

Here, conditional indexing keeps the data that we have deemed "good" by this criteria. We can also do the same for the corresponding time stamps, since t and signal have the same length.

Now let's plot it. You try.

Now that you have your data all cleaned up, it would be nice if we could save it for later and not have to go through the process of cleaning it up every time. Fear not! Python has you covered.

There are two formats that we are going to cover, one that is Python-specific, and the other a simple text format.

First, we must package our two cleaned up arrays into one again. This can be done simply with the np.array() function.

Then, we can use either the np.save() function or the np.savetxt function, the first saving the array into a '.npy' file and the other, into a '.txt' file. The syntax is pretty much the same for each.

Now that your data files are saved, you can load them up again, using np.loadtxt() and np.load() for .txt and .npy files respectively. We used np.loadtxt() above, and np.load works the same way. So, let's load in the .npy file and see if our data was saved correctly.

Now, let's see if you can do the same thing, but with the .txt file that we saved.

So, to summarize, not only can you manipulate arrays, but now you can save them and load them. In a way, those are some of the most important skills in scientific computing. Almost everything you'll be doing requires you know this, and now that you've mastered it, you're well on your way to being an expert in computational physics and astronomy!