human
I am working on the file "Sessions/Session11/Day2/MeasuringCentroidsAndProperMotion.ipynb".
The first task is:
Useful tools
First I'm going to give you one class and one function that will be needed for the first couple of problems.  The class called MoffatPSF describes a PSF with the following circular profile:
$$ I(r) \propto \left[1 + (r/r_0)^2\right]^{-\beta}.$$
Many ground-based PSFs look more like this Moffat function than they do like a Gaussian, so let's use this for our astrometry experiments.  We'll stick with $\beta=2.5$ When you create a MoffatPSF, you give it the  half_light_radius $r_{1/2}$ which encloses half of the light, rather than giving it $r_0.$  You also give the coordinates $(x_0,y_0)$ of the center of the star.  All of these arguments we'll assume to be in units of pixels.
Once you have made a MoffatPSF, you can do a few things with it:
* draw a picture of the PSF as a 2d image (such that it sums to unity)
* realizeXY will return you the (x,y) pixel locations of arrival for one observation of the star - including Poisson noise.
* realizeImage will return you a 2d array that is the pixelized image of the arriving photons.
If you're curious you can look inside the class to see how it's done - I've used a few non-obvious tricks to make this fast since we want to create a lot of star images. There is a common trick of random-number generation when we are picking photon locations, where we transform a uniform deviate (which computers are good at) using a function that creates some desired distribution instead.  But you don't need to understand this to do the problems.
Then there is a function addBackground which will add background noise of a chosen level (denoted as $n$ in the lecture notes) to any image.
Can you write Python code for it?
assistant
#Here is the Python code you requested.

# Make a Moffat class to draw
class MoffatPSF:
    def __init__(self,half_light_radius=4., x0=15.75,y0=15.25):
        # Create a Moffat PSF model that is centered at [y0,x0]
        # with the specified half-light radius (given in pixels)
        self.beta = 2.5
        # Calculate r0 from half-light radius
        tmp = np.power(0.5,1./(1-self.beta)) - 1
        self.r0 = half_light_radius / np.sqrt(tmp)
        self.x0 = x0
        self.y0 = y0
        # Factor that makes the integral of PSF be unity:
        self.norm = (self.beta-1) / (np.pi * self.r0**2)
    def draw(self,n_pix=32):
        # Return an array of shape (n_pix,n_pix) drawing the PSF
        # (using numpy convention
        # of zero-indexing and fastest index, x, last)
        
        # First make arrays holding the x and y coordinates
        xy=np.indices( (n_pix,n_pix),dtype=float)
        x = xy[1].copy()
        y = xy[0].copy()
        rsq = np.square(x-self.x0)+np.square(y-self.y0)
        
        # Here's the basic formla for Moffat:
        result = self.norm * np.power(1+rsq/(self.r0**2), -self.beta)
        # Tune up the normalization to unity:
        result /= np.sum(result)
        return result
    def realizeXY(self,n_photon=1e3):
        # Produce an array of the arrival positions of
        # photons, assuming that they are Poisson sample
        # of a mean of n_photon arrivals.
        # The returned array has shape (2,n_detected)
        # where mean of n_detected is n_photon
        
        # First pick total number of photons:
        n_detected = np.random.poisson(n_photon)
        
        # Then distribute them in radius using the trick
        # of inverse cumulative distribution of uniform distribution.
        r_uniform = np.random.random(n_detected)
        rsq = np.power(r_uniform, 1./(1-self.beta))-1.
        r = self.r0 * np.sqrt(rsq)
        
        # Draw position angles at random and produce output
        theta = np.random.random(n_detected) * (2*np.pi)
        return...