HW 3

 

1)      Write a function called findel.pro that will return the index of an array that is closest to a given value. For example:

                IDL> array = [1.3, 2.4, 3.5, 4.2]

     IDL> value = 3.1           

       IDL> result = findel(value, array)           

       IDL> print,result

                  2

       IDL> print,array[result]

                  3.50000

 

HINT: You may have to experiment at the command line first in order to figure this one out.

 

2)      Write a function called integrate.pro that will determine the integral of an array F defined at x-values X. Use the Trapezoidal Rule (ignore the composite term) for your numerical integration method and avoid all loops! Here’s an example:

 

IDL> x = my_fillarr(0.25, -5, 5)

IDL> f = 3*x^2 + 2*x

IDL> result = integrate(x, f)

IDL> print, result

               250.312

 

                  Where the true answer is 250.

 

HINT: Use shift() and be careful about what you do at the edges since shifting wraps array elements around to the other side.

 

3)      Write a wrapper procedure for the histogram() function and plot procedure called plot_hist.pro. Your procedure should do the following:

1.       If no keywords are set, take a distribution of points and create a histogram using the default settings of histogram()

2.       Plot the histogram with the appropriate x-axis values.

3.       Return the x-axis values through a keyword called xaxis

4.       Have keyword options that allow the user to specify the minimum and maximum values to be included in the histogram, as well as the binsize.

5.       Use the _extra keyword to pass plot options, such as the plot color, axis labels and title to the plot command.

 

Here’s an example:

 

IDL> dist = randomn(seed, 1000)*4

IDL> plot_hist, dist, binsize=1.5, min=-20, max=20, xtitle='X Values', ytitle='Number', title='My Histogram', charsize=2

 

And the resulting plot should look like this.

 

HINT: Use the psym=10 keyword to plot to get the “blocky” plot style.

 

BONUS:

 

In Lecture 3 I showed a method of finding the peaks in a vector of points using numerical derivatives (subtracting a shifted version of a vector from itself). This is something that you’ll find yourself doing over and over again as an observational astronomer, so I couldn’t pass up the opportunity to assign it. So here it is!

 

Write a function called peaks.pro that finds the indices of all of the peaks in a spectrum (vector of points). For example, for a spectrum consisting of two Gaussians:

 

IDL> x = my_fillarr(0.25, -10, 10)

IDL> spec = gaus(x, [2., -5, 1.2])

IDL> spec += gaus(x, [3., 4.5, 1.8])

IDL> plot, x, spec

IDL> pk = peaks(spec)

IDL> print,pk

          20          58

IDL> print,x[pk]

     -5.00000      4.50000

 

HINT: If you plan on using this for your research in the future, you’ll need to be wary of what happens at the extreme ends of the array.