ASTR 310: Computing in Astronomy Homework #7
Due online Thursday, April 7, 2022, by 11:00 AM
For this assignment, your solutions should consist of a PDF file listing your answers and program code and showing your interaction with your program. Please clearly label your answers. For programming problems, you can either write a standalone program and run it from a terminal window, or you can use a Jupyter notebook. If you are using a notebook, please clean it up so that it runs straight through from the beginning and does not include extraneous material. Use Jupyter’s export function to create a PDF file and submit that.
1. Making a map of a volcano on Mars (35 pts)
Copyright By PowCoder代写 加微信 powcoder
Download the image “Olympus_Mons_MOLA_64ppd.jpg” from the exercise page on Moodle. This is an image containing elevation data for the region around the Olympus Mons volcano on Mars, derived from the MOLA instrument on the Surveyor spacecraft.
a) Read the image into Python. The result will be a 3D array containing (R, G, B) values for each pixel. Sum over the axis with length 3 to get total values, then scale the result to the range (-5664, 21185). The result will be a 2D array of elevation values in meters. (5 pts)
b) Determine the grid spacing. The image is centered on latitude 11.875 degrees N, longitude 235.5 degrees E. The radius of Mars is 3389.5 km. The image samples the MOLA dataset at 64 pixels per degree. (5 pts)
c) Make a color map of the data with appropriate image extents and produce a vertical color bar. Label the color bar “Elevation (m)”. (10 pts)
d) Produce an array of contour levels spaced by 1000 m, starting with -5000 m, and ending with 21000m. Using these levels, overlay a contour plot on the image, making the contours black with a line width of 0.5. You will need to use the argument origin=’upper’ to make the contours match up with the image (which has a different default for this keyword). Add appropriate x and y labels and title, then show the plot. (15 pts)
2. Working with FITS X- (65 pts)
Download the file “rp900025n00_bas.fits” from the course website. This is an X-ray photon event file from an 8 ksec (8,000 second) ROSAT PSPC observation of the wind-blown bubble NGC6888. (ROSAT [Röntgensatellit] was a German-led X-ray telescope that flew for about 20 years beginning in 1990. The Position Sensitive Proportional Counters [PSPC] were among its principal instruments). In X- rays most astronomical sources are faint enough that we receive one photon at a time, so rather than recording an image with intensities, detectors usually produce files giving information about the arrival of each photon.
In this assignment you will read in, analyze, and plot results from this file. The file contains an empty primary HDU and five binary table extension HDUs. Each row of the second binary table (STDEVT)
corresponds to a single X-ray photon. The table columns that we care about (see the header for this HDU for the full list) are X, Y, and PI(X, Y) is the position of the photon (see the header for the values of X and Y axis reference pixel [TCRPX1/2], sky coordinate value of reference pixel [TCRVL1/2], and pixel size needed to derive sky positions from table values [TCDLT1/2]). PI is the gain-corrected pulse height value for the photon; the photon energy in keV is approximately equal to the PI value times 0.01. Note: we are going to skip some of the standard X-ray data analysis steps, such as limiting counts to time intervals considered to produce reliable observations, accounting for the detector response matrix, and correcting for dead time. So, the results may look a bit ratty and shouldn’t be published!
a) Write code to read the file and select the data by photon energy so that you have a NumPy record array containing (x,y) positions for photons with energies between 0.4 and 2.0 keV (pulse height values 40-200). You will need to convert the (X,Y) positions from the file into sky coordinates in degrees using the WCS header information discussed above:
x sky coord = TCRVL1 + TCDLT1*(X – TCRPX1)
y sky coord = TCRVL2 + TCDLT2*(Y – TCRPX2)
Use Astropy commands to get the necessary header information from the file. (10 pts)
b) Write code to construct an image array from this selection; that is, produce a 2D mesh covering the region of sky from which photons were received, and set the array values equal to the number of photons that fell within each mesh cell. Use 32 mesh cells. Divide the photon counts by the area of the cell in square arcminutes, the effective area of the detector (assume it to be 150 cm^2), and the “live time” of the observation (use the LIVETIME keyword in the header to get this). At the end of this operation, you will have an image array containing values in units of photons per arcmin^2 per cm^2 per second. (10 pts)
c) Use the image array to produce a color image plot. The plot should have appropriate axis ranges (remember units!) and labels and a title, and a color bar legend should be included. Use whatever color palette you like. If you have done this successfully, you should see photons coming only from a circular region with dark radial “spokes” (the ribs on the detector window) as shown below: (15 pts)
https://heasarc.gsfc.nasa.gov/docs/rosat/pspc.html
d) Modify your code to designate two special regions “A” and “B” on the map. The “A” region should be a circle of radius 4′ centered on J2000 coordinates (20h12m6.9s, +38d21m18s) (convert these to decimal degrees to compare to photon positions). The “B” region should be an annulus centered on the same point with an inner radius of 4′ and an outer radius of 12′. Have your code draw the outline of each region on the image and label the outlines appropriately (“A”, “B”). (15 pts)
e) Return to the original file data and select photons by position so that you have two NumPy record arrays, one containing photon energies for photons arriving in region “A” and the other containing energies for photons in region “B.” Apart from using the positions to select photons from the two re- gions we will now be interested only in their energies.
Write code to create a second plot showing spectra for the two regions (two curves on one plot). The spectra should be shown as “step”-type histograms of photon counts per cm^2 per second per keV versus photon energy in keV. Both x and y axes should be logarithmic, with ranges chosen to allow the curves to fill most of the plot. Include a legend, appropriate axis labels (remember units!), and an appropriate title. In choosing the bin size for each curve (they do not have to be the same), ensure that at least 10 photons contribute to each bin. You will need to experiment with the bin size to find the appropriate value.
When constructing photon energy histograms, keep in mind that the NumPy or Matplotlib histogram routines will return data that are effectively in units of photons per bin. So, to get the spectra in the desired units you need to divide the histograms by effective area in cm^2, live time in seconds, and bin size in keV. (15 pts)
程序代写 CS代考 加微信: powcoder QQ: 1823890830 Email: powcoder@163.com