Matlab Tutorial

We can save this plot as a png file by print  -dpng

or as an encapsulated postscript file by print  -depsc2

These commands created the files “figure1.png” and “figure1.eps”, which we could rename (e.g. “mv figure1.png plot_x-squared.png” on the linux command line) if we wanted to keep them.

We will now read in some data from the file “partdiv.txt”. Looks what's inside this file by typing “less partudiv.txt” on the LINUX command line. Use the up and down arrows and page-up & page-down to navigate the file. When you're done press “q”.

Load the data into matlab with:

            [congress, year, total, dems, reps] = textread('partydiv.txt', '%s %f %f %f %f',78,'headerlines',3)

Be careful here, there is a space before each % sign.

Let's make some plots:

           plot(year,seats)

           xlabel('Year')

           ylabel('Total House Seats')

Now try

          plot(year,dems./total,'-b') 

          xlabel('Year')

          ylabel('Percentage')

We can overplot by turning the "hold on"

          hold on

          plot(year,reps./total,'-r')

          hold off

We can add a legend with

          legend('democrats', 'republicans')

You can move the legend around using the GUI by first clicking the arrow, then selecting the legend box and dragging it.

You can get a sense for how correlated (or anti-correlated) a dataset is by looking at the “correlation coefficient”:

corrcoef(dems./total,reps./total)

gives you a 2x2 matrix of coefficients relating each column to itself (along the diagonal) and to the other column(s) (off-diagonals).

TASK

Now that you have completed your tutorial, let’s do some science!

Inside your home directory, you should find a file called “RecentIndices.txt”. Use Linux to see what is inside the file and what every column of the data is. The first column is the year and the second is the month of the observation. The third column is the observed sunspot number, the 8th is the observed radio flux, and the 10th is the observed AP index (a geomagnetic index).

To complete this task, you must

1)generate an image (png) files that contain the sunspot number as a function of time, the radio flux as a function of time, and the observed AP index as a function of time. Plot the lines in different colors and make a legend

2)generate 2 postscript files, one that contains the scatter plot of the sunspot number on the x-axis and radio flux on the y axis and one that contains the sunspot number vs AP index. Calculate the correlation coefficient between the two values and include that information in the title of the plot.

Be sure to label all axes of your plots.

When you are done with this task, email the results (3 files) to jacob.heerikhuisen@uah.edu

Loops are useful for manipulating a list (array) of data one element at a time. Some of the above tricks of using matrix multiplication are more efficient, but loops are useful for complicated data. Consider the program:

n = 10;

data = rand(n,1); heads = 0;

tails = 0; for k = 1:n

dat = data(k); if (dat > 0.5)

heads = heads + 1;

else

tails = tails + 1;

end

end

disp(['fraction of heads = ', num2str(heads/n)]) disp(['fraction of tails = ', num2str(tails/n)])

TASK: Write a program to generate a histogram

A histogram is a way of sorting and describing date. For instance, if there are 10 people in the room and you sort them into ages (3 20 year olds, 4 21 year olds, and 3 22 year olds) you have made a histogram. The input to the program would be the ages of people, i.e. ages= [20, 21, 22, 21, 20, 20, 21, 22, 22, 21], the output would be the array that describes the ages (bin= [20,21,22]) and the numbers of people in each age bin (hist =[3,4,3]). In your program, you should have at least 1 input (like ages) and at least 2 outputs (like bin and hist). The other thing you have to consider is the size of the bins. For instance, what if your bin was in 10 year increments, bin = [10, 20, 30]? All the ages would be in the 20 bin (h = [0,10,0]) and your histogram would be pretty boring. Alternatively, the bin size could be too small, you could have increments of days or even minutes. In that case, the bin array would be huge and there would only be one person in each one bin so h would be mostly 0 with 10 1s. So binsize could also be an input into your program or you should figure out a way to make the program to calculate the correct binsize.

Once you have written you histogram program, try it out on some real data. Make a histogram of the sunspot number data we looked at in the first tutorial. Make a histograms of the number of months as a function of sunspot number.

Write the plot into a eps file and send it to jacob.heerikhuisen@uah.edu along with your histogram program.

MATLAB tutorial 3: writing codes to solve a problem

An Orbit Around the Sun

We will try to compute the orbit of an object (assumed massless) around the Sun. Open a new text file and start by clearing all variables:

clear;

Next let's define the length of an astronomical unit (AU) and the number of seconds in a day:

We are going to integrate the equations of motion in thexy-plane: dxdt =vx ; dydt =v y ; dvdtx =ax ; dvdty =ay

The acceleration can be computed from the gravitational force:

a −GM ⃗r

⃗= r2 ̂

where the vector r is a unit vector, and the mass of the object cancels out. A unit vector is a vector of length 1, obtained by taking the position vector ⃗r and dividing it by its length |r|. Using this information, and the fact that ⃗r =( x , y) gives:

The four differential equations require four initial conditions in order to have a unique solution. Let's choose the following:

xx = 1 * au; yy = 0; vx = 0; vy = 29800;

As we integrate over the orbit, we need to keep track of the position of the object, so that we can plot out the orbit once it's been calculated. It's most efficient to “declare” the arrays where we will store the x & y positions up-front:

npoints = 1000;

x_vect = zeros(1,npoints); y_vect = zeros(1,npoints);

Finally, we need to choose a time-step for the integration procedure. dt = 1*day

We are going to use the simplest integration method, known as “Euler's method”. The basic idea is that at each point in time you take a small step in the direction of the local velocity, and then update the velocity according to the local force. The rest of the program looks something like:

x_vect(1) = xx/au; y_vect(1) = yy/au; for kk = 2:npoints