from IPython.display import Image, display
%%bash
open petersonjr_events_r.mpg
This requires modeling the entire time-dependent system, typically in the following order:
In some sense, you "back out" all the effects that impact your data, to yield your final measurements of the sky. These data are then used to test models (frequentist approach), or as ground truth to understand which models are consistent with the data (Bayesian approach).
from IPython.display import Image, display
Image(filename='psf.png')
Astronomical data are typically acquired using Charge Coupled Devices (CCDs), mounted on the back-end of a telescope.
CCDs are silicon-based detectors whose main purpose is to convert incoming photons into electrons via the photoelectric effect. This conversion has a wavelength-dependent sensitivity that must be calibrated in the laboratory.
Electrons are stored in discrete regions within the silicon called "pixels", which are defined by electrodes embedded in the silicon, and prevented from drifting by electric fields.
An "image" is essentially a matrix of pixels, typically in units of 1024 or 2048-ish pixels ("2k by 4k").
A "camera" is a mosaic of CCD devices, frequently edge-buttable to provide a high fill factor of silicon in the focal plane.
Image(filename='Keplerspacecraft-FocalPlane-cutout.png')
To get the electrons out of the silicon and onto a computer, the device must be "read out". During this process, electrons in a row of pixels are shifted down 1 row, with the bottom row being read-out by electronics.
The conversion of electrons to "analog-to-digital units" (ADUs) happens with a gain factor that is characteristic of the readout amplifier, which converts charge into voltage. This process also happens with a certain amount of random per-pixel readnoise, typically a couple of counts.
Because the arrival of the photons can be considered a Poisson process, the uncertainty in the number of photons N received is $\sqrt{N}$, and the uncertainty in the number of ADUs that are stored in your image is $\sqrt{N \times gain}$. Accordingly, the baseline uncertainty on the "counts" in any measured object is $\sqrt{N \times gain + readnoise^2}$.
Image(filename='ccd_fill_factor_big.png')
However, there is also pattern noise in the readout of the image that comes from the process of shifting the electrons around. All of this means that a zero-second exposures does not have zero counts.
The pattern noise is not apparent in any single image, but can be seen in a stack of images.
The way to model this is to take a set of zero-second exposures, with the shutter closed, and coadd them together to make a master "bias" frame. This bias frame is then subtracted from all subsequent exposures.
Image(filename='bias.jpg', width=512)
There are other, second order, effects that also belong here and include:
* Dark current (a time-dependent bias)
* Bad pixels in the detector with degraded or no response to photons
* Fringing (interference of photons with each other in the silicon)
* Non-linearity of the device (the absence or presence of electrons in the pixel make the storage of other electrons more/less likely)
* Crosstalk between CCDs
All of these are time-dependent effects and must be modeled frequently:
* Bias and dark current should be a slow drift with time and temperature of your CCDs
* Dark current will be exposure-time dependent
* The number of bad pixels increases with time, due to e.g. cosmic ray hits
* Non-linearity will most severely affect bright stars
* Fringing is dependent on the instantaneous emission-line properties of the atmosphere
* The amount of crosstalk between CCDs will have a slow drift with time, but will be different starfield-to-starfield
Having modeled the electronic parameters of our collection device, we must model the "above the silicon" portion of our experiment.
This includes all telescope support structures, mirrors, and lenses on the camera. Because of reflections and dust, certain portions of the system will be more or less sensitive to photons, meaning the response of your system is not spatially "flat".
The standard way to model this is to map out the light path that photons take from the aperture of the telescope.
Typically, the telescope points at an illuminated dome screen or a blank part of the sky and opens the camera shutter. Ideally the resulting image should have the same number of counts per pixel, to within the Poisson noise. However...
Image(filename='flat.png', width=512)
A sequence of these flat field images are acquired, bias corrected, and then coadded into a master flat field frame which is then divided from all science frames, making their response to photons spatially flat.
The atmosphere + telescope + camera fully determines the point-spread-function (Psf) of an image. The Psf is the transfer function of a point source through the system, and defines the fundamental resolution limit of a given image.
Image(filename='Zeta_bootis_short_exposure.png', width=128)
Its shape is driven by high altitude (Kolmogorov) turbulence in the atmosphere that provides $1/r^{5/3}$ broad wings, with a central core with a characteristic scale that is approximately $1/r^{2}$.
Stars may effectively be considered point sources (their angular diameters are $10^{-6}$ seconds of arc) and are considered unresolved. Galaxies are extended objects and frequently have sizes larger than the Psf; if so, they are considered resolved. Accordingly, stars provide a very good empirical measurement of the Psf in any given image.
The size of the Psf (full width at half maximum of the Psf profile) is typically referred to as the image "seeing", with typical seeing values being around one second of arc (1 arcsecond = 1/3600$^{th}$ of a degree).
The median seeing for LSST will be 0.7 arcseconds, with images in bluer passbands having slightly worse seeing than images in redder passbands.
The 0.1 arcsecond seeing for HST is not set by any atmospheric turbulence but by the diffraction limit of the instrument itself.
Image(filename='seeing.jpg')
Resolution and depth
You can build the point spread function of an image by using postage-stamps of all the stars. However, we don't know a priori which objects are stars and which are faint but resolved galaxies.
In addition, each object is pixellated by the detector, leading to difficulties in understanding the shape of a given object (i.e. answering the question "is it resolved or unresolved"). If the detector undersamples the Psf, than this process is difficult because the shapes of all objects are fundamentally impacted by the coarseness of the pixel sampling.
i1 = Image(filename='oversampled1.jpg')
i2 = Image(filename='undersampled.jpg')
i3 = Image(filename='undersample2.png')
display(i1,i2,i3)
If the detector critically samples the Psf, this process is far easier. The notion of critically sampling the Psf comes from the Nyquist theorem, which essentially defines how frequently you have to sample a continuous signal to say you have "recovered" it with perfect fidelity.
More intuitively, you want the Fourier decomposition of your signal to be zero outside of some frequency range, and you want your pixels to sample this signal at higher frequency than this band limit.
If this is the case, you can fully describe your signal with a finite number of terms in a Fourier decomposition, and there is no information lost in the discrete sampling process.
from IPython.display import SVG
SVG(filename='Bandlimited.svg')
As a rule of thumb, if the FWHM of the Psf is 2-3 pixels (or larger) than you have critically sampled your Psf.
For LSST the plate scale of the camera is 0.2"/pixel.
For a given device, your best-seeing nights give you your sharpest (and typically deepest, for a given exposure time) data. However, it is also the most pixellated of your images and thus can be more difficult to analyze.
The process of building a Psf is non-trivial, but will be the subject of the next talk.
Performing optimal measurement on an image actually requires you to answer the question: "What am I trying to measure?". It is this profile that defines the optimal filter for detecting objects your data.
The optimal filter for extracting a signal from noisy data is to convolve the data with the signal you are looking for. This generates a maximum likelihood estimate of there being an object of that nature in your data. I.e. if
$I(x,y) = A \times PSF(x,y) + \epsilon(x,y)$
Thus the maximum likelihood estimate for the detection of a star is
$I(x,y) \ast PSF(x,y)$
NOTE: the optimal detection of galaxies requires a different filter than the optimal detection of stars.
Measurements typically are of two classes: aperture shapes and counts, and model (or "weighted") shapes and fluxes. Aperture fluxes are calculated trivially as, whereas model fluxes using a weighting function $w(x,y)$ for the sum.
For the measurements of stars, Psf-photometry provides the optimally weighted flux measurement.
Calibration of your data is necessary to compare to other surveys that have done the same work.
This involves matching the measured centroids of your objects to an absolute astrometric reference frame ([x,y] -> [ra,decl]), and the measured brightnesses of your objects to physical units (counts -> Janskys).
Absolute astrometric calibration is typically performed by assuming that your image is a tangent plane projection on the spherical celestial coordinate system. This image will have a central pointing, a rotation, and typically other high order terms that are set by optical properties of the telescope. The telescope's "guess" at where is it pointing is typically used as a starting point, and a catalog of astrometric standard stars is used to fit the mapping of image centroids to celestial coordinates. A "world coordinate system" (WCS) transform represents this mapping of pixel to sky coordinates. Blind solvers like astrometry.net may also be used if you don't know where your telescope is pointing.
Second order effects that make this process non-trivial include
The creation of an absolute astrometric reference system is a far more difficult process; thank the USNO for their hard work on the USNO-B and UCAC catalogs. The ideal astrometric reference frame, ICRS, uses a network of ~200 distant quasars to define a non-moving coordinate system.
The calibration of photometry requires converting your brightness measurement (in counts) to physical units.
For single-object observations, this involves observations of a set of "standard stars" throughout the night. By comparing the measured counts of these standard stars in your image to the expected counts, you can set a photometric zeropoint for each image.
You then may fit a model to this zeropoint throughout the night to interpolate to your science (non-standard-star) observations.
The discussion of how these standard stars became standard is again a story. Terrestrial sources (copper and platinum heated to several thousand degrees such that they glow like blackbodys) were used to calibrate the star Vega as a transfer standard.
However, Vega was found to be a non-standard star, both in the fact that its spectrum exhibits excess in the infrared, and the fact that it is intrinsically variable.
Standard star networks, based on the hard work of e.g. Arlo Landolt and Peter Stetson, have been established around the sky, but fundamentally depend on the spectrum of Vega.
Image(filename='VegaSpectrum.png', width=600)
Recent surveys have defined a more sane magnitude system, called the AB system based not upon any single object but on an absolute standard in Janskys.
If you have your data in cgs units (erg s$^{−1}$ cm$^{−2}$ Hz$^{−1}$):
$m_{Vega}$ = -2.5 log$_{10}$ (f$_\nu$/f$_{Vega}$)
$m_{AB}$ = -2.5 log$_{10}$ (f_$\nu$) - 48.6
As a concrete example, the SDSS Survey has an internal (relative) self-calibration that is consistent to 1%, based on the repeat measurement of stars.
Image(filename='dr9_spectro_coverage.png', width=800)
This internal magnitude system is itself calibrated to within 2% of the absolute AB system.