HW 5

NO LOOPS ALLOWED!

1) Write a function called makecirc.pro that creates a 2-dimensional circle-mask consisting of 0’s outside of a given radius and 1’s inside. The inputs are a 2-D array (perhaps a sub-image of a CCD frame containing a star), and an array containing both the center of the circle (in pixel units) and the radius. For example, for a circle centered on pixel (3,3) and a radius of 2:

IDL> data = round(randomu(seed, 7, 7)*10+1)

IDL> mask = makecirc(data, [3,3,2])

IDL> print,data

5 5 5 9 6 9 10

6 2 10 7 5 11 9

9 8 6 11 1 9 7

9 3 11 10 9 3 6

3 9 3 9 10 3 1

5 7 9 7 7 8 10

2 10 3 3 5 2 9

IDL> print,mask

0 0 0 0 0 0 0

0 0 0 1 0 0 0

0 0 1 1 1 0 0

0 1 1 1 1 1 0

0 0 1 1 1 0 0

0 0 0 1 0 0 0

0 0 0 0 0 0 0

IDL> print,mask*data

0 0 0 0 0 0 0

0 0 0 7 0 0 0

0 0 6 11 1 0 0

0 3 11 10 9 3 0

0 0 3 9 10 0 0

0 0 0 7 0 0 0

0 0 0 0 0 0 0

HINT: You’ll need one array of x-values and another for the y-values. This is your chance to use the rebin() and reform() tricks I showed you in class.

2) Write a function called gauss2d.pro that creates a 2-dimensional Gaussian profile

using the equation:

G(x, y) = A₀ + A₁ exp(-U/2)

where

U = (x’/a)² + (y’/b)²

and

x’= (x-h)cos T – (y-k)sin T

y’= (x-h)sin T + (y-k)cos T

The inputs are a 1-D arrays containing the x- and y-axis, and an array A containing the parameters defining the Gaussian, where

A[0] = A0 = constant offset

A[1] = A1 = amplitude

A[2] = a = s width of Gaussian in the X direction

A[3] = b = s width of Gaussian in the Y direction

A[4] = h = X centroid

A[5] = k = Y centroid

A[6] = T = Theta, the rotation of the ellipse from the X axis in radians.

Make the rotation, T, optional. If the user sends in a 6-element array, set the rotation to zero.

For example:

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

IDL> y = x

IDL> a = [0., 1.0, 1.5, 0.75, 0, 0, -!pi/4]

IDL> g = gauss2d(x, y, a)

IDL> shade_surf, g, x, y, charsize=3

The output should look like this.

3) Write a function called sum_gauss.pro that computes N Gaussians and totals them together to form a final 1-D spectrum. The N Gaussians are specified by a 4-by-N (column, row) array of parameters. The columns contain the amplitudes, centroids, 1-sigma widths. The final column contains the y-offsets for the N Gaussians, but the y-offsets are optional and is set to zero if the input array contains only 3 columns. Pass the Npix-by-N array of Gaussians out through an optional keyword called gaussarray. For N=4:

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

IDL> b = transpose([[1., 0.1, 0.3, 0.8], [-0.5, -1.1, 0.4, 1.0], [1.0, 0.7, 0.3, 1.5]])

IDL> print,b

1.00000 -0.500000 1.00000

0.100000 -1.100000 0.700000

0.300000 0.400000 0.300000

0.800000 1.00000 1.500000

IDL> g = sum_gauss(x, b, gaussarray=garr)

IDL> title = 'Funky Spectral Line'

IDL> yt = ‘Intensity’

IDL> xt = 'Velocity [km/s]'

IDL> plot, x, g, xtitle=xt, ytitle=yt, title=title, chars=1.5

IDL> for i = 0, 3 do oplot, x, garr[*, i], lines=2

The plot looks like this.