Computational Astronomy Assignment 3
====================================
Due date: Friday 17 April.


1) Neither scipy or numpy seems to have a routine to fit a straight line to a set of points which have differing uncertainties. So part (1) of this week's assignment is to write your own in python. I just wrote my own version to check that it works. In order to set you on the right theoretical track, I reproduce here the 'def' line and the docstring from my own code:

def generic_line_fit(x_vals, y_vals, y_uncs=None):
  """This function fits a straight line of the form

	y = intercept + slope * x

to the input data. According to section 6-3, Bevington (1969), the best-fit values of intercept and slope obey the following equation:

	( sum[1/s/s]    sum[x/s/s]  )   ( intercept )   (  sum[y/s/s]  )
	(                           ) * (           ) = (              )
	( sum[x/s/s]   sum[x*x/s/s] )   (  slope    )   ( sum[x*y/s/s] )

where I've written x, y and s respectively for x_vals, y_vals and y_uncs. The function solves this matrix equation by using the solve() function of the numpy.linalg module.

Note that use of numpy is unavoidable, so I have imported within the function definition.
"""


2) You are observing an unpolarized, flat-spectrum, pointlike source with a radio telescope of system temperature (which I take in this case to include all the non-source contributions) 30 K, an effective area of 2800 m^2 and a bandwidth of 32 MHz. You spend half your time observing the source and the other half pointing at source-free sky background. You subtract the average on-source temperature from the average off-source temperature to obtain the nett source temperature, then convert that to janskys.

	a) What will be the uncertainty in your measurement of average off-source temperature if your total observing time is 1 hour? (Ie 1/2 hour on-source, 1/2 hour off-source.)

	b) What will be the uncertainty in the calculated nett source temperature? (Assume the source temperature is negligible compared to the Tsys.)

	c) What then is the approximate minimum source flux (in janskys) observable in that time? (Boltzmann's constant k is 1.38e-23 J K^-1.)


3) Measurements of a partially polarized radio source give the following sets of Stokes parameters (all in mJy):

	Wavelength:	   I		  Q	  U	  V
	  (cm)		(+- 0.1)             (all +- 0.05)
	======================================================
	    6		 29.1		-0.09	 0.57	-0.03
	   10		 42.7		 0.13	 1.04	-0.00
	   21		 74.6		 2.65	-0.22	 0.05
	   30		 97.4		-3.01	-4.46	-0.01
	======================================================

	a) Plot log(I) against log(lambda).
		i) Why are the values of I consistent with a power-law spectrum?
		ii) Use your line-fitting function from part (1) to estimate the spectral index.
		iii) What is the most likely physical process powering the radio emission?

	b) Calculate the polarization fraction of the source at 6 cm.

	c) Is there a significant amount of circular polarization at any wavelength? Give your reasons.

	d) Calculate the polarization angle for each wavelength, and the uncertainties in these angles.

	e) Plot the polarization angle against the square of the wavelength.

	f) Use your line-fitting function from part (1) to deduce the rotation measure (in degrees cm^-2) and the 'true' polarization angle of the source.

	Some hints:
		- The 'math' module has functions you will need such as sqrt() and atan2(), plus constants such as pi.
		- The polarization angle (ie the angle between the direction of vibration of the electric field, measured anticlockwise from horizontal) is just 1/2 the angle made by a Stokes vector in the Q-U plane of the Poincare sphere. Obviously the same will be true of uncertainties.
		- Because the Q-U angle is uncertain modulo 2*pi, the polarization angle is uncertain modulo pi. Hence, on your plot in part (e), you may have to add or subtract n*pi to one or more values in order to get all the points to 'line up' properly.
		- If sigma_x is approximately equal to sigma_y, sigma_angle (in radians) is approx equal to sigma_x / sqrt(x*x+y*y).