Removing the background

Written by Matin Sharkawi

Why do we need to remove the background?

Unfortunately, London is plagued by a constant orange haze of light pollution. For the purposes of our project, light pollution harms our data as it adds a seemingly random value to all areas of the image, making it impossible to read the actual brightness of celestial bodies.

In order to remove the noise pollution, there are a few possible methods one could take. We tried two methods with varying degrees of success.

When taking the photos of each celestial body, we included an extra offset shot of the nearby sky. This shot was to provide us with a reference from which we could determine values to remove from our main photos, in order to remove the noise.

Method 1 – Minima removal

The simplest method of using the offset image is to compare each image with the offset image and generating an array of minima between the two images, as shown below.

Code in order to generate a "mud" and remove it from the original image.
Figure 1 – Code in order to generate a “mud” and remove it from the original image.

As described, the code generates an image showing the minima between the main images and the offset image, then subtracts this ‘mud’ from the original image to remove the noise pollution. The function backround_removal also needs error handling for underflow errors, as when values go below zero, they wrap back around to the maximum 32 bit integer. This can be avoided by implementing an if statement that checks if any values drop below zero, and returns them to zero.

However, as is now visible, artifacts in the form of dark holes appear on the image with the mud removed that line up with the stars visible in the offset image. In an ideal scenario, you would attempt to find a patch of sky with as few stars as possible in order to avoid this.

Image 1 - The "Mud". Note how the holes line up with the brightest stars in Figure 2 below.
Figure 2 – The “Mud”. Note how the holes line up with the brightest stars in Figure 2 below.
The offset image, which lines up with the holes in the mud
Figure 3 – The offset image, which lines up with the holes in the mud

Method 2 – Gradient Removal

n order to get rid of the artifacts, we instead found the gradient between the offset image and the 10 images, and from this we generated a noise map that represented the background in each image. Adding another error handler to avoid underflow errors, we then subtract the noise map from the original image, resulting in a clean, artifact free image from which we can extract data.

The alternate method, generating a 2D background image that can then be subtracted from the original
Figure 4 -The alternate method, generating a 2D background image that can then be subtracted from the original.
The noise map generated by Fig. 3
Figure 5 – The noisemap generated by figure 3. As you can see, there are no visible stars left to cause artifacts, just a map of the skies brightness levels.
One of the final images with no background.
Figure 6 – One of the final images with no background.

Star clusters and their HR diagrams

Written by Elvi Haaramaki and Sof Hachem

Star clusters

In this project, we are looking at objects called star clusters, which are groups of stars that are gravitationally bound together. There are two types – globular clusters and open clusters.

Globular clusters are large collections of thousands of millions of red stars, most of which are population II stars formed a few million years after the Big Bang. The stars are densely packed together, and the shape of the cluster is roughly spherical1

Figure 1: Globular cluster M3.

Open clusters are made up of a few hundred stars which are much less densely packed together than globular clusters, and are without a distinct shape. These stars are younger blue stars, tens of millions of years old, and are weakly gravitationally bound2. For this project we took pictures of M35, M36, and M37, which are all open clusters. 

Open clusters are ideal for the purposes of the project, since the stars are further apart from each other, making them easier to resolve individually. This is necessary to be able to find the apparent magnitudes of the stars. Additionally, it is assumed that all the stars in the cluster are the same distance from us, allowing us to calculate the absolute magnitudes.

Figures 2 and 3: Open clusters M36 and M37.

HR diagrams

HR diagrams are a useful tool in astronomy, as plotting stars within a cluster based on their magnitude and colour-temperature allows the study of the stars’ relative stellar evolutionary stages and provides an idea about the age and the relative distances of given clusters3, which was also the objective during the project. 

For the data here, green and blue filters4 were used to attain the difference in the magnitudes corresponding to a colour index in the observational HR diagram, being plotted against the apparent (green) magnitude.  

Figures 4 and 5: HR diagrams for M36 and M37.

The upper left region in the diagrams corresponds to the bluer and brighter stars, while the redder and dimmer stars are located at the bottom right. The stars in a HR diagram for a relatively young open cluster, which mostly contains stars burning hydrogen into helium, are expected to lie along the main sequence branch – this is seen in the HR diagram for M36.  In older clusters, the main sequence cut-off becomes more visible as the more massive stars burn hydrogen at faster rates and are consequently spending less time in the main sequence, leading older clusters to have more prominent red-giant branches. In globular clusters, more stellar evolutionary stages are present due to the clusters age, which can be seen in the HR diagram for M3.   

Figure 6: HR diagram of M3, showing more distinct branching and stars in variety of evolutionary stages6.

Identifying the turnoff point from the diagram and comparing the mass and brightness of the most massive remaining star in the main sequence, allows us to further estimate the age of the star and – as stars within the cluster are estimated to have formed around the same time – the age of the cluster5







6: (Figure 6)

Removing the Hot Pixels

Written by Stylianos Spyrou

In this post, the concept of the hot pixels will be analyzed, the importance of removing the hot pixels from the FITS files will be emphasized and parts of the coding and the reasoning behind the steps followed will be explained.

What are the hot pixels and how are they obtained?

Hot pixels are single sharp pixels located at random locations of common images taken by digital cameras. They are points that do not react linearly to incident light captured by the lens. An important feature of the hot pixels is that they appear in the same location regardless the frame. This means they do not move, and they remain at a fixed position.

Hot pixels are caused due to electrical charges which appear into the sensor well of the camera’s lens. Since they are very sharp, and they are individual extremely bright pixels, they cannot be observed while looking at the image taken. However, they are very easy to determine while zooming at the images closely during processing, especially if the background of the image is very dark, such as the images taken from a telescope. They can also be visible when the sensor’s temperature increases or at very high ISOs[1]. The weather conditions during the photoshoot also have impact on the existence of the hot pixels since, if the temperature of the surroundings is high, it is ideal for their formation of hot pixels on the lens during the shoot. Finally, they tend to appear far more often in long exposure images. The reasoning behind the formation of the hot pixels is that while capturing less light from the scenery in the given moment, the patterns obtained by the camera sensor are comparatively stronger in that specific moment.  Lenses and camera sensors get hotter and hotter as they use long exposures.

Why & how to remove hot pixels?

It is very important to remove the hot pixels from an image. The main reason is that it affects the image significantly during close viewing or printing. The only way they can be removed is by following a processing method after the image is taken. This can either happen through photo editing or coding. In this particular article, the second method will be analyzed and explained in detail.

Figure 1- defining the function that will find the hot pixels of the FITS files. In this part of the code, median filter is being used, however, it can be replaced by a gaussian filter instead.

On Figure 1, the function that determines the hot pixels of each FITS file must be defined. In the first part of the code, the edges must be ignored, and they will be found separately later on. In order to understand the purpose behind it, we need to imagine the filter being a square with the center passing from every point of the image obtained. As the detector approaches the edges of the image, the corners of the square filter will exceed the limits and therefore the values they will obtain will be set by Python to zero, which affects the rest of the data (median/Gaussian filter).

Figure 2- Finding the hot pixels at the edges.

As indicated in the Figure 2 above, the next step of the code is to find the pixels on the edges of the image, but not the corners. Repeat the process for the left and right edges, top and bottom.

The next step is to find the hot pixels of the image at its four corners. The code used is very similar to the one used to determine those at the edges and can be indicated in the Figure 3 below:

Figure 3- Finding the hot pixels at the four corners of the image.

Using the previous Figures, the hot pixels can be obtained. Since the code would be required to be repeated for all images in order for them to get stacked, a part that repeats the following method can be introduced, which keeps the code simple and neat. However, the easiest method is to simply repeat the code for each FITS file of each cluster. A suggestion can be found below on Figure 4:

Figure 4 – Removing the hot pixels for or images before stacking.


Fits File Handling

Written by Sarthak Maji


When the data of the star clusters are collected using a CCD camera, the data is saved as a FITS file. To handle this FITS file in python, first, we need to understand what FITS files are. FITS (Flexible Image Transport System) is a file format designed to store, transmit, and manipulate the data saved on the file. The data stored in a fits file is multi-dimensional, like a 2D array. The pixel data from the camera is stored as a 2D or 3D array depending on the type of image. For example, if the data is monochrome, it is stored as a 2D array; if it is RGB, it is stored as a 3D array. The image metadata is stored as headers in ASCII format, which is readable by humans.

Handling FITS file in Python

To handle the FITS file in python, we first need to install a package called AstroPy. AstroPy is installed using the command line function on our compiler, which is python. We also install other valuable packages like photutils which detect stars and their positions on the image.

Figure 1: Installing AstroPy Package using command line function.

In the command line, we use the pip command to install the latest version of AstroPy and other useful packages. Once the needed packages are installed, we must import all necessary packages onto the jupyter notebook. For example, to open the fits file, we need a specific function of astropy to handle fits; the fits function is imported from, as shown in the code below.

import numpy as np
import matplotlib.pyplot as plt
from import fits
from matplotlib.colors import LogNorm as log
from astropy.stats import SigmaClip
from photutils.background import Background2D, MedianBackground

Once the packages are imported, we import the fits file onto our jupyter notebook and save it to a variable, as shown in the code below. Every fits file is saved into different variables.

Img1 = ('')
Img2 = ('')
Img3 = ('')
Img4 = ('')
Img5 = ('')
Img6 = ('')
Img7 = ('')
Img8 = ('')
Img9 = ('')
Img10 =('')
Imgo = ('')

To ensure the fits files are imported correctly, we must also ensure that the fits files are saved in the same directory as the jupyter notebook. After importing the fits file onto a variable, as shown in the code above, we now need to open them to manipulate the data, which is done using the function, as shown in the code below.

hdul1 =
hdul2 =
hdul3 =
hdul4 =
hdul5 =
hdul6 =
hdul7 =
hdul8 =
hdul9 =
hdul10 =
hdulo =

Once the fits file is opened onto a variable, we now need to extract the data and store it onto another variable so that the data can be manipulated. The code to extract data is shown in the code below.

Data1 = hdul1[0].data.astype('int32')
Data2 = hdul2[0].data.astype('int32')
Data3 = hdul3[0].data.astype('int32')
Data4 = hdul4[0].data.astype('int32')
Data5 = hdul5[0].data.astype('int32')
Data6 = hdul6[0].data.astype('int32')
Data7 = hdul7[0].data.astype('int32')
Data8 = hdul8[0].data.astype('int32')
Data9 = hdul9[0].data.astype('int32')
Data10 = hdul10[0].data.astype('int32')
Datao = hdulo[0].data.astype('int32')

We extract the data as a 32-bit file (as shown in the code above) so that while manipulating the data, we avoid problems like stack overflow while processing the images. The extracted data can remove hot pixels, background, and stacking. This extracted data is a 2D array where every data point is the pixel value. We can also determine the shape of the array, which gives us an idea of the no. of pixels in the camera sensor.

Figure 2: Extracted data as an array and the shape of the array.
Figure 3: Raw data of M35B.

Figure 3 is the image representation of the extracted raw data of M35B in the log scale. This raw data contains hot pixels and a background, which is removed before stacking all the raw images collected.

Telescope Group 2023: Measuring the Age of a Stellar Cluster

Written by Sara Fekri


Hello everyone! We are the Telescope Group 2023: Sof, Stel, Sara, Sarthak, Elvi, Matin, and Malcom. As part of our 3rd year research project, we are looking at different star clusters in London’s dark skies to measure the age of the stars within the clusters and produce a Hertzsprung-Russell (HR) diagram for them. 

Here, we will be presenting how the FITS files, an astronomy file format, are handled in Python and the steps taken to get the HR diagram of the star clusters. The objects we have captured are the open clusters Messier 35, Messier 36, and Messier 37 on the night of 08/February/2023. We also observed the Orion Nebulae M42, and globular cluster M3.

The steps taken for our project each have a detailed blog post and are as follows:

HR Diagram

A Hertzsprung-Russell (HR) diagram is plot of stars showing the relationship between the stars’ absolute magnitude versus their effective temperatures. The absolute magnitude corresponds to the stars’ luminosities with negative values being brighter. The effective temperature can be shown using the stellar classification, where the stars from O to M correspond to decreasing temperate on the x-axis.1

An important feature in the HR diagram is the main sequence (MS), the region where the primary hydrogen-burning takes place within a star’s lifetime. Looking into detail where the MS ends, the “MS turn-off”, gives astronomers information about the age of the star cluster.

The general layout of the HR diagram can be seen below, where stars of greater luminosity are located at the top of the diagram, and stars with higher surface temperature are towards the left side of the diagram. 

Hertzberg-Russel diagram. (Image taken from here)

Our Cluster Data – M35

The data captured using the telescope at King’s gives us a FITS file, the most common digital file format in astronomy, and were captured with the blue and green filters attached. From that, we were able to open the files and get on with extracting the data. For this post, data from the M35 blue filter will be used. The raw image for M35_blue from Python can be seen:

Raw image of one of the data files for M35 Blue filter.

However, this raw image does not only contain the stars; it includes hot pixels captured by the filters and a background. The histogram below shows the numbers of counts vs the pixels of the image above:

Histogram of the data set above with the number of counts on the x-axis and the number pixels on the y-axis.

The brightest stars would be located on the histogram where the sharp peaks occur at the right of the curve. The highest peak of the histogram on the left side corresponds to the background of the image.

Next, the hot pixels needed to be removed from the image (detailed post here). Below is a picture showing the stars detected on the M35_Blue image after removing hot pixels using a code in Python:

Stars from M35_Blue are detected with a red circle around them.

The following procedure involved finding the background for it to be subtracted from each data set. The detailed methodology of background removal can be accessed here. The same steps were taken for each of the 10 data sets captured for M35_blue and M35_green. Now, the 10 images can be stacked to form:

Stacked image of M35 after removal of background and hot pixels in the blue filter.

The step afterwards was to identify the stars in the blue and green stacked images by comparing the two and see which stars overlap, i.e. correspond to the same star in the different filters.

After pairing the stars, the pixel counts of each star in both the blue and green images are found to compute the star’s apparent magnitudes (B and G). Finally, G vs. B-G is plotted, and this would form the observational HR-diagram of M35.

Observational HR diagram of open cluster M35.

The plot above displays HR diagram of M35, where the Main Sequence branch is clearly shown. The MS turn-off point is where the stars leave the MS and form “Giants” which are seen on the top right of the diagram. Since most stars are on the MS branch, and there are very few giants, this concludes that the open cluster M35 is fairly young.


This project aimed to produce an HR diagram of M35 starting from extracting the data in the FITS files to comparing the stars in the stacked images of the blue and green filter. The same methodology was applied to M36 and M37. This post provided a summary of the steps taken to reach the final result where details can be found in the appropriate separate posts.

Hope you enjoyed reading!

2020 Vision: Spectral Classification

By Ceinwen Cheng

The most typical way to classify stars is based on their spectral profiles. Electromagnetic radiation is passed through filters displaying lines on a spectrum. The intensity of each spectral line provides information on the abundance of the element and the temperature of the star’s photosphere.

The system for classification is named the Morgan-Keenan (MK) system, using letters O, B, A, F, G, K and M, where O-type is the hottest and M-type is the coolest. Each letter class has a further numerical subclass from 0 to 9 where 0 is the hottest and 9 is the coolest. In addition to this, a luminosity class (Yerkes Spectral Classification) can be added using roman numerals, based on the width of the absorption lines in the spectrum which is affected but the density of the star’s surface.

The aim of  our very last imaging session was to image spectra, including stars from all the classes in the MK-system:

  • 0-Class: Alnitak
  • B-Class: Regulus and Alnilam
  • A-Class: Alhena and Sirius
  • F-Class: Procyon
  • G-Class: Capella
  • K-Class: Alphard
  • M-Class: Betelgeuse

For each star, two sets of spectra were recorded at a 60 s exposure set on 2×2 binning. As the exposure time is relatively long, it was important to manually keep the star between the cross-hairs of the eyepiece while imaging, as there is a tendency for the star to stray from the center of the telescope as the earth rotates.

There is a trend of worsening spectra as you move down the spectral types. As luminosity is proportional to the fourth power of temperature, cooler stars in the lower end of the spectral classes such as K and M-classes are dimmer than an O-class. Displayed below are the spectra we collected for each star, excluding Alphard and Betelgeuse as no identifiable spectral lines were seen for both stars.

spectra alnilam B

Alnilam B-class: For this class we expect medium amounts of Hydrogen Balmer lines, and some neutral Helium lines.  We observe spectral lines of    H-γ, H-δ,  He I, and C III in our experimental spectra analysis.

spectra alnitak O

Alnitak O-class:  In the hottest class of stars, we expect to see ionized helium features. In this graph, we see prominent absorption lines of H-β, Hγ, Hδas well as He I, He II and He III, which are in correlation with our expectations.

spectra capella G

Capella G-class: We expect heavier elements such as calcium and for the Hydrogen Balmer lines to be less prominent. Compared to stars in A or F-classes, our experimental Capella spectra have less defined H-ε, Hγ, Hδ absorption lines, but a more obvious Ca II absorption line.

spectra sirius A

Sirius A-class: This class has the strongest features of the Hydrogen Balmer series. We observe H-β, Hγ, Hδ, Hε strong spectral lines, and slight H-ζ and H-η absorption lines.  Additionally, characteristically of A-class stars, Sirius displays spectral lines of heavier elements such as Mg II and Ca II.

spectra alhena

Alhena A-class: Alhena’s spectra are very similar to Sirius’s, being in the same class.

spectra regulus B

Regulus B-class: Displaying prominent H-β, Hγ, and Hδ Balmer lines, it is in the lower temperature end of B-class, almost an A-class. It has strong Balmer features but still a characteristic He I absorption line that places it in the B-class.

spectra procyon F

Procyon F-class: We expect to see ionized metals and weaker hydrogen lines than A-classes. Our data shows weaker H-β, Hγ, and Hδ Balmer lines than Alhena or Sirius, but more prominent absorption lines for ionized elements such as Ca II .

Stars can be assumed to be black bodies, and its absorption spectra will overall take the shape of its black body radiation. After extracting the absorption spectra, we attempted to normalize it by dividing the spectra data by its polynomial interpolation to the fourth degree. The polynomial interpolation is our estimate of the black body curve that would be displayed by the star.

Characteristically, an ideal normalize spectrum will not deviate from 1.0 on the y-axis, other than peaks where the absorption lines are. We can directly compare that the peaks in the normalized graphs are where the troughs in the initial spectra data are.

normalized alhena

Normalized spectra, Alhena A-class

normalized alnilam B

Normalized spectra, Alnilam B-class

normalized alnitak O

Normalized spectra, Alnitak O-class

normalized capella G

Normalized spectra, Capella G-class

normalized procyon F

Normalized spectra, Procyon F-class

normalized regulus B

Normalized spectra, Regulus B-class

normalized sirius A

Normalized spectra, Sirius A-class

Code: Extracting spectra and normalizing data

From here, we will go through how we extracted the spectra from the .fit files using Python, as well as the normalization of the spectra through a polynomial interpolation. This is for 1 star, specifically Alhena, where the code is repeated for the rest of the stars as well, where the only deviation is the addition of absorption lines manually.

Step 1:

#importing relevant modules
import numpy as np
import os
import math
from PIL import Image
from PIL import ImageOps
import astropy
%config InlineBackend.rc = {}
import matplotlib
import matplotlib.pyplot as plt
%matplotlib inline
from import fits
from astropy.nddata import Cutout2D
import os
from numpy import asarray
import matplotlib.image as mpimg
from specutils.spectra import Spectrum1D, SpectralRegion
from specutils.fitting import fit_generic_continuum
from numpy import *
import numpy.polynomial.polynomial as poly
from scipy import interpolate

Step 2:

#opening both files to check stats, and positioning a cut out around both spectra
hdulista1 ="C:\Users\Ceinwen\Desktop\ProjectPics\Spectra\")
dataa1 = ((hdulista1[0].data)/256)

positiona1 = (838.5, 660)
sizea1 = (150,1677)
cutouta1 = Cutout2D(np.flipud(dataa1), positiona1, sizea1)
plt.imshow(np.flipud(dataa1), origin='lower', cmap='gray')

hdulista2 ="C:\Users\dhruv\Desktop\Project Pics\Spectra\")
dataa2 = ((hdulista2[0].data)/256)

positiona2 = (838.5, 690)
sizea2 = (150,1677)
cutouta2 = Cutout2D(np.flipud(dataa2), positiona2, sizea2)
plt.imshow(np.flipud(dataa2), origin='lower', cmap='gray')

We open both .fit files to check the stats of the data, and cut out the part of the image with the spectra for both data sets, as you can see, the spectra are very faint.

Step 3:

#plotting the cut out

plt.imshow(, origin='lower', cmap='gray')

plt.imshow(, origin='lower', cmap='gray')

Plotting the cut-out, taking a closer look, the spectra lines are more visible now.

Step 4:

#converting the file's data into an array and plotting it

xa1 =
a1 = np.trapz(xa1,axis=0)

xa2 =
a2 = np.trapz(xa2,axis=0)

Step 5:

#finding the mean of both data and plotting it, alongside spectral lines.

a = np.mean([a1, a2], axis=0)

array with lines

Step 6:

#doing a line of best fit through the spectral lines's positions.

m1, b1 = np.polyfit(x1, y1, 1)

straight line

Step 7:

#mapping the pixel number to wavelength
xo1 = np.arange(1,1678)
func1 = lambda t: (m1*t)+b1
xn1 = np.array([func1(xi) for xi in xo1])
plt.xlabel ('Wavelength ($\AA$)')

in wavelenght

Step 8:

#plotting the spectra in term of wavelength and including named spectral lines.
plt.xlabel ('Wavelength ($\AA$)')
plt.axvline(x=4481.325,color='magenta',label="Mg II",ls=":")
plt.axvline(x=3933.66,color='lime',label="Ca II",ls=":")

in wavelenth plus lines

Step 9:

#doing a polynomial fit to the 4th degree, plotting it on the same graph as the spectra
coefs1 = poly.polyfit(xn1, a, 4)
ffit1 = poly.Polynomial(coefs1)
plt.plot(xn1, ffit1(xn1))
plt.xlabel ('Wavelength ($\AA$)')

polynomial interpolation

Step 10:

#diving the spectra by the polynomial fit to normalise, and plotting normalised graph

plt.xlabel ('Wavelength ($\AA$)')

Finally, dividing the spectra by the polynomial interpolation to normalize the spectra.



Brief Note:

The end of the 2020 third year project came quicker than most years before us, unexpected situations such as strikes and social distancing due to coronavirus have unfortunately prevented us from collecting any more data on the variations in the apparent magnitude of Betelgeuse, detailed in the previous blog post. We have had to cancel our poster and presentation as well.

We are all saddened by the fact we can’t clamber about the dark rooftop that overlooks the Thames and half of London anymore, but I believe this project has inspired many in the group to go further in the field of Astronomy. We all want to thank our absolutely brilliant supervisor, Prof. Malcolm Fairbairn, for his guidance and leadership. He has been looking out for us in more than this project, being an amazing mentor.

2020 Vision: Tracking of Betelgeuse – Recent dimming and brightening

Written by Ceinwen Cheng

Betelgeuse, what used to be the 10th brightest star in the night sky, lies on the left shoulder of Orion’s constellation. This reddish star attribute’s its color to its classification of being a red supergiant around 15 times the mass of the Sun.

However, ever since October 2019, Betelgeuse’s apparent magnitude has had a drastic decrease that can be observed with the naked eye, dimming from an apparent magnitude of 0.6 in October 2019 to 1.8 in early February. The apparent magnitude scale is reverse logarithmic, i.e the lower the magnitude, the brighter the star. A star of apparent magnitude 1.0 is 2.512 times brighter than a star of apparent magnitude 2.0.


A graph published by Nature news article in February 26th c 2020 displaying a graphical representation of Betelgeuse’s apparent magnitude

Additionally, as recently as February 22nd, Betelgeuse has started brightening again, an increase of almost 10% from the dimmest point. There are many proposed explanations, one being a change in extinction and therefore not a change in the actual luminosity of the star, which is supported by the fact that the infrared observations have of Betelgeuse has barely changed since the dimming phenomenon. In astronomy, extinction is the absorption and scattering of light particles by cosmic dust between the observer and the light emitter.

We acknowledge several problems measuring the apparent magnitude of Betelgeuse in Central London, one of the most light polluted cities in Europe. If we measure Betelgeuse directly, our data will be skewed by a range of obscure interference such as clouds, atmospheric seeing, thermal turbulence, air and light pollution. These factors are hard to mitigate and could alter or change the trend of our data. Therefore, we approach this project by also tracking a star close to Betelgeuse that apparent magnitude isn’t changing much, Betelgeuse’s sister star Bellatrix. We are interested in the ratio of luminosity between Bellatrix and Betelgeuse.

Imaging both Bellatrix and Betelgeuse in 5 different sessions in 2020, on the 28th Jan, 3rd Feb, 11th Feb, 17th Feb, 2nd  March 2020,  we collected the images to be processed into a trend. Using Python 3, Dhruv Gandotra, our python expert, extracted the data from our images to produce the following trend in the graph below, showing the trend of the ratio of magnitudes between Betelgeuse and Bellatrix. Additionally, he also found a source online [twitter handle @Betelbot] that tracks the apparent magnitude of Betelgeuse over the same time frame as ours, allowing us to peer review our data.

Data points betelbot and ours (5 sessions)

@Betelbot ‘s measure of Betelgeuse apparent magnitude (in red), compared to our experimentally calculated difference in apparent magnitude between Betelgeuse and Bellatrix inverted.  Where Day 0 is the 28th Jan.


trendline betel bot and ours

@Betelbot ‘s data and ours combined, with trend line calculated for both sets of data.


Brief Description of our code:

Importing our .fit files (a file format used in astronomy) into python, there are first and foremost several important points we need to take into account. All the images have noise, some images have hot pixels, and the position of the stars differ for each image. To solve these problems at the same time, Dhruv came up with an idea, to select a picture form a session, use an array to find the position of the maximum pixel brightness, then define a box around the star of fixed width, and only take data from that box. This way, any hot pixels can be cut out, the amount of noise will be relatively the same for the pictures in each session, and as the position of the stars do not move around too much, the star will still be in the box of data.

Note: We will be attaching parts of the code alongside its explanation for future references

After importing all our relevant modules into our code, we must first preview a single Betelgeuse image from the session, looking at the file’s information and then the data available in the file. From this we can extract the minimum, maximum, standard deviation, or mean of the pixel brightness.


#importing relevant modules
import numpy as np
import os
import math
from PIL import Image
import astropy
%config InlineBackend.rc = {}
import matplotlib
import matplotlib.pyplot as plt
%matplotlib inline
from import fits

#opening single image and it's data
from astropy.nddata import Cutout2Dhdulist1a ="C:\Users\dhruv\Desktop\Project Pics\Betelgeuse\Session 1\")
data1a = ((hdulist1a[0].data)/256)

print('Min:', np.min(data1a))
print('Max:', np.max(data1a))
print('Mean:', np.mean(data1a))
print('Stdev:', np.std(data1a))

The next part of the code involves finding the position of the maximum pixel brightness i.e. where the star will be located, through an array. From this, we create a 50×50 cut out around the star, and graphed the cut out. Now that we have a preview on an image, we can repeat this easily for the 10-15 other images taken per session, add all their data and find the mean average of their data.


#locating position of star and creating a cut out
position1 = (223, 248)
size1 = (50,50)
cutout1a = Cutout2D(np.flipud(data1a), position1, size1)
plt.imshow(np.flipud(data1a), origin='lower', cmap='gray') #using np.flipud as the image is inverted
cutout1a.plot_on_original(color='white') #displacing the position of the cut out in the image

#plotting cut out
plt.imshow(, origin='lower', cmap='gray')

#importing the rest of the images taken in the session
hdulist2a ="C:\Users\dhruv\Desktop\Project Pics\Betelgeuse\Session 1\")
hdulist3a ="C:\Users\dhruv\Desktop\Project Pics\Betelgeuse\Session 1\")
hdulist4a ="C:\Users\dhruv\Desktop\Project Pics\Betelgeuse\Session 1\")
hdulist5a ="C:\Users\dhruv\Desktop\Project Pics\Betelgeuse\Session 1\")
hdulist6a ="C:\Users\dhruv\Desktop\Project Pics\Betelgeuse\Session 1\")
hdulist7a ="C:\Users\dhruv\Desktop\Project Pics\Betelgeuse\Session 1\")
hdulist8a ="C:\Users\dhruv\Desktop\Project Pics\Betelgeuse\Session 1\")
hdulist9a ="C:\Users\dhruv\Desktop\Project Pics\Betelgeuse\Session 1\")
hdulist10a ="C:\Users\dhruv\Desktop\Project Pics\Betelgeuse\Session 1\")
hdulist11a ="C:\Users\dhruv\Desktop\Project Pics\Betelgeuse\Session 1\")
hdulist12a ="C:\Users\dhruv\Desktop\Project Pics\Betelgeuse\Session 1\")
hdulist13a ="C:\Users\dhruv\Desktop\Project Pics\Betelgeuse\Session 1\")
hdulist14a ="C:\Users\dhruv\Desktop\Project Pics\Betelgeuse\Session 1\")

#extracting the data of each image
data2a = (hdulist2a[0].data)/256
data3a = (hdulist3a[0].data)/256
data4a = (hdulist4a[0].data)/256
data5a = (hdulist5a[0].data)/256
data6a = (hdulist6a[0].data)/256
data7a = (hdulist7a[0].data)/256
data8a = (hdulist8a[0].data)/256
data9a = (hdulist9a[0].data)/256
data10a = (hdulist10a[0].data)/256
data11a = (hdulist11a[0].data)/256
data12a = (hdulist12a[0].data)/256
data13a = (hdulist13a[0].data)/256
data14a = (hdulist14a[0].data)/256

#stacking the images and plotting
sum1 = data1a+data2a+data3a+data4a+data5a+data6a+data7a+data8a+data9a+data10a+data11a+data12a+data13a+data14a
average1 = sum1/14
ax1a = plt.plot(average1)

average data plot session1 betel

All 14 of the images of Betelgeuse taken in session 1 stacked, where the x axis is the number of pixels???, and the y-axis is the relative mean average brightness of?? the stacked images.

stacked cut out session1 betel

The graph above shows the cut out relative to the actual image taken of Betelgeuse in Session 1.

stacked box session1 betel

This graph displays the actual cut out, where we shall be focusing on taking our data from.

An interesting thing we can find in our stacked image is the amount of background noise there is, which can be seen if we reduce the y-axis stacked data to a limit of 1.10 to 1.50. Background noise is light registered by the CCD from a space in the sky where there are no light sources. Background noise is often caused by light diffusion or scattering in our atmosphere, and telescopes in space such as the Hubble Space Telescope do not encounter such a problem.


#changing the y-axis of the graph to show noise
ax1a = plt.plot(average1)

noise session1 betel

A zoomed in graph of the earlier stacked images to a y-axis limit of 1.10 to 1.50, displaying the noise present in the images.

stacked cut out data session1 betel

This is the sum of the stacked cut out data after converting the image into an array. One thing to point out is that there are to be two peaks, i.e. the position of the star moved. This does not change anything as we are finding the summed maximum pixel brightness.


#repeating finding the position of the star and creating cut out for the stacked image
position2 = (224, 253)
size2 = (50,50)
cutoutf1 = Cutout2D(np.flipud(average1), position2, size2)
plt.imshow(np.flipud(average1), origin='lower', cmap='gray')
plt.imshow(, origin='lower', cmap='gray')
xf1 =
ax1b = plt.plot(xf1)
mn1 = np.sum(xf1)

print('Min:', np.min(xf1))
print('Max:', np.max(xf1))
print('Mean:', np.mean(xf1))
print('Stdev:', np.std(xf1))

Finally, we convert the stacked images into an array, and sum all the data, as seen in the graph above. The curve we produced is seen to be a nice Gaussian distribution. As we are looking for the ratio of apparent brightness between Bellatrix and Betelgeuse, we do not need the apparent magnitude. The whole process is repeated for Bellatrix, and then the sum of the respective data is placed into a ratio and normalised 2.512 (due to the reverse logarithmic scale of apparent magnitudes), to give the experimental difference in apparent magnitudes.


#where mn1 is the sum of brightness for Betelgeuse and mn2 for Belletrix 
#finding the difference in apparent magnitude
a = math.log(mn1/mn2,2.512)
print a
((a-1.22)/1.22)*100  #percentage error

Graphs of Interest

stacked data session 1 betel

Discarded stacked data graph of failed session, with a hot pixel seen as a green line on the left.

noise session 1 betel

The same graph with the y-axis reduced to 1.0-3.5 mean photons received in 0.1 s seconds. Still observing the hot pixel as well as an increase in noise past the x-axis 535 pixels mark.

When using the telescope, the CCD (Charged-Coupled device) needs to be cooled, up to around -31ºC. However, the graphs above shows the result of taking images when the CCD is not sufficiently cooled, in this image, the camera was cooled to around -7ºC. Hot pixels occur when the camera sensor heats up, an electrical charge is leaked, resulting in a pixel that is much brighter than expected. In our data, there is an obvious hot pixel, as there is a pixel that is very bright but not located at where the star is, additionally, there is no Gaussian spread in it neighbouring pixels. The hot pixel should go away after the CCD is sufficiently cooled.

Rapidly cooling the CCD may cause ice crystals to form in the detector, often dry ceramic pills placed in the detector to reduce the amount of moisture. On the far right of our data graphs, we see one of the stacked photos have a large amount of noise compared to the other photos. A few plausible explanations for such a phenomenon: ice crystals forming due to the CCD rapidly cooling, or a naughty undergrad using their phone near the telescope, increasing the background luminosity.


The span of our allocated project time is minuscule compare to many astronomical events. The dimming of Betelgeuse occurred over 4-5 months and should still be monitored in the future due to it sudden brightening. By the time this blog is published, we have only taken 5 session of images of Betelgeuse and Bellatrix, and to produce a trend or a timeline worth looking at, we need many more sessions and images. Clouds remain our worst enemy. Springtime in London brings rain and cloudy days, and we cannot image on a cloudy night.

However, even though Malcolm dislikes the company of optimists, there are still a few more weeks before the end of the allocated project time, where we will continue imaging Betelgeuse and Bellatrix in hopes of providing you with a more complete trend of Betelgeuse’s dimming and brightening.

2020 Vision: Imaging of Messier-42, Orion Nebula

Orion’s constellation remains one of the most identifiable set of stars in the night sky, and right below the belt, resides Orion’s nebula, also known as Messier-42. We took several images of Messier-42 at 20 second exposures, using 5 different filters. Each filter corresponds to a specific wavelength on the spectrum of light, where filter 1 to 5 are blue, hydrogen-alpha, oxygen-3, Sulphur, and green, respectively.

We stacked the photos using the image editor Gimp, where we colourised the photos from each filter, and then stacked them again so they formed the final image. The photos corresponding to filter 1 and 5 were not included in the final stack.

aadama.png Aadam’s edit


Alex’s Edit


Conor’s edit

Stacked (RGB)

Dhruv’s edit


Vlad’s edit


Ceinwen’s edit

Our photos can hardly compare to the startling images taken by the Hubble Space Telescope, the many drawbacks of imaging in central London. There are 3 main factors impacting the resolution of our images:

  • Air pollution: the absorption and back-scattering of cosmic radiation, reducing the intensity of the received light
  • Light pollution: artificial light from nearby buildings causes the night sky to take on a glow (due to scatter from molecules in the air), worsening contrast in images.
  • Atmospheric Turbulence: differences in air pressure and temperature changes the refractive index of air, induces scintillation effects such as changes in visual position and fluctuations in brightness due to refraction. Makes the stars seem like they are “dancing” or “twinkling”

These atmospheric distortions are somewhat mitigated by prolonging the exposure and stacking of images to remove background.

Quick introduction to the project group of 2020! We are a team of third year (mostly) physicists consisting of Ceinwen, Dhruv, Aadam, Conor (formerly biochemistry!?), Alex and Vlad.


Alex, Conor, Dhruv and Malcolm spelling out LOVE with their hands and internally their hearts


Aadam enjoying the cold frigid January weather at midnight and definitely not thinking about his warm bed


Malcolm fixing our mistakes


Wild nights on the roof

Written by Ceinwen Cheng


Telescopes Project 2019: Simulating an Exoplanet Transit

This post was written by Tim who is taking part in the project this year:-

One of the tasks we hoped to achieve in our project was to observe the transit of an exoplanet. By taking images of a star over an extended period, we should see a small dip in the luminosity as an exoplanet passes in front of it.

The dip in flux of photons tends to be very small (in the order of 0.2-0.5%), so before embarking on this task, we needed to ask the question: What transits are we likely to be able to detect? Tim went about simulating a transit to see how reliable our data is likely to be.

Note: At the end of each section, the Python code that was used is included for future reference.

Step 1: Finding the error in flux from our images

The first step was to find out how much the flux varies in our images of stars that are not undergoing a transit. From working on our HR diagram, we had 9 photos (of each filter) of Messier 47. A star was chosen at random and the standard deviation of photon counts for that star across all 9 B images was calculated, giving us a standard deviation of 0.094516 (as a fraction of total photon count) for a 20 second exposure.


### 0. Locating objects in FIT file ###

import numpy as np, pandas as pd
import imageio
import matplotlib.pyplot as plt

from scipy import ndimage
from scipy import misc
from import fits
from scipy.spatial import distance

from collections import namedtuple
import os

def total_intensity_image(input_voxels, threshold):
    voxels_copy = np.copy(input_voxels)
    threshold_indices = voxels_copy < threshold
    voxels_copy[threshold_indices] = 0

def threshold_image(input_voxels, threshold):
    voxels_copy = np.copy(input_voxels)
    threshold_indices = voxels_copy < threshold     voxels_copy[threshold_indices] = 0     plt.imshow(voxels_copy)     plt.title('Image with threshold at '+str(threshold))      def get_optimal_threshold(input_voxels, start=2, stop=9, num=300):     """     The optimal threshold is defined as the value for which the      second derivative of the total intensity is maximum     """     array = np.linspace(start,stop,num)     intensities = [total_intensity_image(input_voxels, i) for i in array]     array_deriv1 = array[2:]     intensities_deriv1 = [(y2-y0)/(x2-x0) for x2, x0, y2, y0 in zip(array_deriv1, array, intensities[2:], intensities)]     array_deriv2 = array_deriv1[2:]     intensities_deriv2 = [(y2-y0)/(x2-x0) for x2, x0, y2, y0 in zip(array_deriv2, array_deriv1, intensities_deriv1[2:], intensities_deriv1)]     return(array_deriv2[np.where(intensities_deriv2 == max(intensities_deriv2))]) GFiltering = namedtuple('GFiltering', ['filteredImage','segrObjects', 'nbrObjects' ,'comX', 'comY']) #defining named tuple containing (1) labeled voxels,\ #(2) integer number of different objetcs, (3) array of x-ccord of centers of mass for all objects, (4) same for y def gaussian_filtering(input_voxels, blur_radius):     # Gaussian smoothing     img_filtered = ndimage.gaussian_filter(input_voxels, blur_radius)     threshold = get_optimal_threshold(input_voxels)     # Using NDImage     labeled, nr_objects = ndimage.label(img_filtered > threshold) 

    # Compute Centers of Mass
    centers_of_mass = np.array(ndimage.center_of_mass(input_voxels, labeled, range(1, nr_objects+1)))
    array_com_x = [i[1] for i in centers_of_mass]
    array_com_y = [i[0] for i in centers_of_mass]
    return(GFiltering(img_filtered, labeled, nr_objects, array_com_x, array_com_y))

### 1. Retrieving coordinates and photons counts from a folder of fits files ###

# Take a single fits file, return data from with coordinates, sizes and photon counts of each star
def get_stars_from_fits(file):
    voxels =[0].data
    data = gaussian_filtering(np.log10(voxels), 7) # This just find object locations, we won't actually use log10
    stars = []
    for star_label in range (1, data.nbrObjects + 1):
        objectCoord = np.where(data.segrObjects == star_label)
        zipped_list = np.array(list(zip(objectCoord[1], objectCoord[0], voxels[objectCoord])))
        coords = np.average(zipped_list[:,:2], axis=0, weights=zipped_list[:,2])
        size = len(objectCoord[0])
        photon_count = voxels[objectCoord].sum()
        stars.append({'coords': coords, 'size': size, 'photon_count': photon_count})
    return pd.DataFrame(stars)

# Repeat above function for an entire folder (returns list of DataFrame, each item being data from one image)
def get_stars_from_folder(folder):
    filenames = os.listdir(folder)
    data = []
    for filename in filenames:
        if '.fit' in filename:
    return data

### 2. Aligning images based on largest star so we can later match stars between files ###

# Get data from largest star in an image
def get_largest_star(image):
    return image[image['size'] == image['size'].max()].iloc[0]

# Get coordinates of largest star in image
def get_largest_star_coords(image):
    return get_largest_star(image).coords

# Take the full data set and return the image with most stars for calibrating against
def choose_calib_image(images):
    max_stars = max([image.shape[0] for image in images])
    return [image for image in images if image.shape[0]==max_stars][0]

# Take an image and align so largest star matches location with calibration image
def align_against_calib_image(image, calib_image):
    largest_star_coords = get_largest_star_coords(image)
    calib_image_coords = get_largest_star_coords(calib_image)
    image.coords = image.coords + calib_image_coords - largest_star_coords
# Take set of images and align them all against calibration image
def align_all_against_calib_image(images, calib_image):
    for image in images: align_against_calib_image(image, calib_image)
# Take the full data set, get all images, choose calibration images automatically and align all against it
def align_all_auto(images):
    align_all_against_calib_image(images, choose_calib_image(images))
    return images

### 3. Matching stars between multiple images to group into one DataFrame ###

# Get distance to all stars in file from specified coords. Outputs list of tuples (star object, distance)
def get_star_distances(coords, image):
    star_distances = [] # list of tuples: (star, distance)
    for i, star in image.iterrows():
        star_distances.append((star, distance.euclidean(coords, star.coords)))
    return star_distances

# Get closest star to particular coords in an image. Returns tuple (star object, distance)
def get_closest_star(coords, image):
    star_distances = get_star_distances(coords, image)
    smallest_distance = min([star_distance[1] for star_distance in star_distances])
    return [star_distance for star_distance in star_distances if star_distance[1]==smallest_distance][0]

# Take coordinates (x,y), and go through images selecting nearest stars. Ignore if too far away.
def find_matching_stars(coords, images, max_distance=10):
    closest_stars = []
    for image in images:
        (star, distance) = get_closest_star(coords, image)
        if distance < max_distance:
    return {'coords': coords, 'closest_stars': pd.DataFrame(closest_stars)}

# Returns list of tuples (x,y) of star coordinates from an image
def list_all_coords(image):
    return [star.coords for i, star in image.iterrows()]

# Repeat find_matching_stars for all stars in a file, output pandas DataFrame
def match_all_stars(image, images, max_distance=10):
    full_star_list = []
    for coords in list_all_coords(image):
        matching_stars = find_matching_stars(coords, images, max_distance)['closest_stars']
        if len(matching_stars) == len(images): # Only include stars found in all files
            full_star_list.append({'coords': coords, 'matched_stars': matching_stars})
    return pd.DataFrame(full_star_list)

# Same as match_all_stars but takes entire data set and chooses calibration files by itself
def match_all_stars_auto(images, max_distance=10):
    images = align_all_auto(images)
    return match_all_stars(choose_calib_image(images), images, max_distance)

### 4. Finally retrieving the data we actually want ###

def aggregate_star_data(images):
    data = match_all_stars_auto(images)
    stars = []
    for i, star in data.iterrows():
        coords = tuple(star.coords)
        photon_counts = star.matched_stars.photon_count.values
        sizes = star.matched_stars['size'].values # Had to use [''] notation because of my stupid naming
            'coords': coords,
            'avg_photon_counts': np.mean(photon_counts),
            'avg_size': np.mean(sizes),
            'standard_deviation_photon_counts': np.std(photon_counts),
            'coeff_of_var_photon_counts': np.std(photon_counts) / np.mean(photon_counts)
    return pd.DataFrame(stars)

folder = 'FIT_Files/' # Make sure to include '/' after folder name
data = get_stars_from_folder(folder)
aggregate_star_data = aggregate_star_data(data)
aggregate_star_data.to_csv('CSV/combined_star_data.csv') # Save the data for use later

Step 2: Finding a functional form for light curves

Next, to help generate functions for simulating a transit, data was taken from a transit of HIP 41378, which was observed by the Kepler space telescope (, dataset K2SFF211311380-C18). This shows how flux changes over a few days. At around 3439 days, there is a clear dip in the data where a transit took place:


Zooming in on this dip, we start to see the shape of the “light curve” we are looking for:


After a bit of trial and error and help from Malcolm, a curve was fitted over this which could be used for simulations.


The formula for the curved part finally came out as

Screenshot from 2019-03-19 16-42-55


  • Φrel = Relative flux (1.0 being flux when no transit)
  • A = Amplitude (maximum dip in relative flux)
  • T = Length of time of transit (days)
  • t = time (days)
  • t0 = Time of minimum flux (i.e. time halfway through transit)

With this equation, the parameters can now be varied as we like.



### 0. Imports ###

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import os

### 1. Get some sample data from an actual exoplanet transit ###

# Read the data and plot the relative flux over time
filepath = 'Kepler_Data/test-light-curve-data.csv'
sample_data = pd.read_csv(filepath)
sample_data = sample_data.rename(index=str, columns={"BJD - 2454833": "time", "Corrected Flux": "relative_flux"})
sample_data.plot(x='time', y='relative_flux')

# zooming in to the transit (secton with big dip) and normalising so flux starts at 1.0
start_time = 3437.55
end_time = 3438.61
transit_data_length = end_time - start_time#

sample_data = sample_data.drop(sample_data[(sample_data['time'] < 3437.55) | (sample_data['time'] > 3438.61)].index)
sample_data['relative_flux'] = sample_data['relative_flux'] / sample_data['relative_flux'][0]
sample_data.plot.scatter(x='time', y='relative_flux', ylim=[0.994, 1.001])

# Adjusting to start at time=0
sample_data['time_adjusted'] = sample_data.time - start_time
sample_data.plot.scatter(x='time_adjusted', y='relative_flux', ylim=[0.994, 1.001])

### 2. Creating a functional form for this data ###

## 2.1 Fitting a curve to the data##

# Just a quick function to help in the formula below
# Return x^exponent if x>=0, (-x)^exponent if x<0 def real_power(x, exponent):     if x >= 0: return x**exponent
    else: return -(-x)**exponent

relative_flux_range = relative_flux_range = max(sample_data.relative_flux) - min(sample_data.relative_flux)
time = np.arange(0.0, transit_data_length, 0.01)
time_centre_of_transit = transit_data_length / 2

# After fiddling around the below formula seems to fit for during the transit (but not outside)
relative_flux_curve = [1.0 - relative_flux_range*(real_power(np.cos((4.2*(x-time_centre_of_transit))), 0.2)) for x in time]

transit_curve = pd.DataFrame({'time': time, 'relative_flux': relative_flux_curve})
# Set all points outside of the transit to 1.0
transit_curve.loc[transit_curve['relative_flux'] > 1, 'relative_flux'] = 1.0

ax = transit_curve.plot(x='time', y='relative_flux')
sample_data.plot.scatter(ax=ax, x='time_adjusted', y='relative_flux', color='orange')

## 2.2 Switching to units of seconds and photons per second ##

expected_photons_per_second = 80000 # This is a random choice for now but not far off true data we'll use
seconds_in_day = 86400

# Convert time from days to seconds
transit_curve['time_seconds'] = list(map(int, round(transit_curve.time * seconds_in_day)))

# We only generated data for each 0.01 days, let's fill the curve for every second
def interpolate_transit_curve(curve):
    # Start with list of empty data points
    time_range = max(curve.time_seconds)
    data = pd.DataFrame({'time_seconds': list(range(0, time_range + 1))})
    relative_flux = [None] * (time_range + 1)

    # Fill each row that we have data for
    for index, row in curve.iterrows(): # Populate data we have
        relative_flux[] = row.relative_flux 
    if not relative_flux[0]: relative_flux[0] = 1.0
    if not relative_flux[-1]: relative_flux[-1] = 1.0
    data['relative_flux'] = relative_flux
    # Fill in the gaps
    data = data.interpolate()
    return data

adjusted_transit_curve = interpolate_transit_curve(transit_curve)

# Changing y axis to photons per second
adjusted_transit_curve['photons_per_second'] = adjusted_transit_curve.relative_flux * expected_photons_per_second

adjusted_transit_curve.plot(x='time_seconds', y='photons_per_second')

### 3. Putting this all into python functions ###

# Change units to seconds and photons per second
def change_transit_curve_units(transit_curve, expected_photons_per_second):
    transit_curve['time_seconds'] = list(map(int, round(transit_curve.time * seconds_in_day))) # Convert time from days to seconds
    adjusted_transit_curve = interpolate_transit_curve(transit_curve) # Fill in missing data
    adjusted_transit_curve['photons_per_second'] = adjusted_transit_curve.relative_flux * expected_photons_per_second # Changing y axis to photons per second
    return adjusted_transit_curve

def curve_formula(x, relative_flux_range, transit_length_days):
    time_centre_of_transit = transit_length_days / 2
    cos_factor = (1.06 / transit_length_days) * 4.2 # This is just variable used in the equation for the curve
    # After fiddling around the below formula seems to fit for during the transit (but not outside)
    return 1.0 - relative_flux_range*(real_power(np.cos(cos_factor*(x-time_centre_of_transit)), 0.2))

# transit_length (in seconds) actually slightly longer than actual transit - includes some flat data either side
# expected_photons_per_second is for when no transit taking place
# relative_flux_range is the amount we want the curve to dip

def get_transit_curve(transit_length, expected_photons_per_second, relative_flux_range):
    second_in_day = 86400
    transit_length_days = transit_length / seconds_in_day
    time = np.arange(0.0, transit_length_days, 0.001) # List of times to fill with data points for curve
    relative_flux_curve = [curve_formula(x, relative_flux_range, transit_length_days) for x in time]
    transit_curve = pd.DataFrame({'time': time, 'relative_flux': relative_flux_curve})
    # Set all points outside of the transit to 1.0
    transit_curve.loc[transit_curve['relative_flux'] > 1, 'relative_flux'] = 1.0
    return change_transit_curve_units(transit_curve, expected_photons_per_second)

# Just a test. This is closer to what we'll be observing
get_transit_curve(10000, 80000, 0.0025).plot(x='time_seconds', y='photons_per_second')

Step 3: Simulating the observation of a transit

The first simulation was done using the parameters for the transit above: T = 1.06 days, A = 0.0053475. A data point was generated for every 30 seconds of data, simulating the photon count for a 20 second exposure based on the given curve. The data point was then shifted randomly up or down based on a normal distribution around that point with the standard deviation taken from part 1 above. In theory, we should have seen the data points dip and rise (in similar fashion to the curve they were generated from), but given the very large standard deviation relative to the change in flux, they came out looking quite random:


Step 4: Testing the data

The question must now be asked: Does this data show anything, or is it just a random mess? To test this, we take a selection of curves of varying amplitudes (from zero change in flux to twice that of the expected change), and test which curve our data fits best against. Ideally, we end up with it fitting closest to the one in the middle (the amplitude we started with for our simulation).


For each curve, a chi-square test was done, comparing the simulated data to the curve. The chi-square statistic is defined as

Screenshot from 2019-03-19 16-55-26

where observedk are our simulated data, and expectedk are the expected data points taken from a curve. σ is the standard deviation.

In general, the lower the chi-square stat, the closer our data fits the model. So, what we hope to see as we vary our amplitude is that the chi square stat decreases (as amplitude increases from zero), hits a minimum at our “good” amplitude (about 0.05, which the data was generated from), and then begins to increase again as the amplitude grows larger. It turns out this is exactly what we saw! The orange line represents the amplitude we were hoping for:


It appears that the best fit for the random-looking simulated data is indeed the curve we were hoping for. This is good news for our chances of observation, but unfortunately this simulation was a bit unrealistic for our purposes: The upcoming transits that Mohyudding has found for us that we are hoping to observe are less than three hours long (the one above is closer to 25), and the dip in flux is around 0.25%, rather than 0.5%. Rerunning the simulation and chi-square testing with a 170 minute transit and amplitude of 0.0025, we get a slightly different result:


This isn’t too bad – it indicates that we would still detect the dip in flux, but perhaps would not be able to accurately measure the amount of dip (in this case it would have been over-estimated).

Unfortunately, there was still one final kick in the teeth awaiting us: The stars we are likely to observe are of around magnitude 11 (very dim). It seems that the dimmer the star, the worse our error in measurement. When retesting with the standard deviation of a magnitude 9 star (less dim, but still quite dim), we get this result:


Oh dear! This result implied that no amplitude (i.e. no transit at all) would be a better fit than any dip in flux, let alone the exact dip we were hoping for. Unfortunately, it seems that our errors are just too large to detect a transit accurately.

In conclusion, it seems that detected an exoplanet transit is going to be quite a challenge. We haven’t completely given up hope though: If we are lucky enough to have a clear sky on the right night, we might try to observe a transit using much longer exposure photographs. This will hopefully allow us to reduce the error in our readings.


### 0. Imports and constants ###

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import os

seconds_in_day = 86400

# Data taken from star 14 in file A-luminosity-standard-deviations.ipynb
test_expected_photons_per_second = round(15158415.7/20)
test_standard_deviation = 1432719

# Data from the imported Kepler data in file B-exoplanet-transit-curve-function.ipynb
test_transit_time = 91584
test_relative_flux_range = 0.0053475457296916495

### 1. Functions to generate curves for simulating transits ###

# (mostly just copied from previous section) #

def interpolate_transit_curve(curve):
    # Start with list of empty data points
    time_range = max(curve.time_seconds)
    data = pd.DataFrame({'time_seconds': list(range(0, time_range + 1))})
    relative_flux = [None] * (time_range + 1)

    # Fill each row that we have data for
    for index, row in curve.iterrows(): # Populate data we have
        relative_flux[] = row.relative_flux 
    if not relative_flux[0]: relative_flux[0] = 1.0
    if not relative_flux[-1]: relative_flux[-1] = 1.0
    data['relative_flux'] = relative_flux
    # Fill in the gaps
    data = data.interpolate()
    return data

# Change units to seconds and photons per second
def change_transit_curve_units(transit_curve, expected_photons_per_second):
    transit_curve['time_seconds'] = list(map(int, round(transit_curve.time * seconds_in_day))) # Convert time from days to seconds
    adjusted_transit_curve = interpolate_transit_curve(transit_curve) # Fill in missing data
    adjusted_transit_curve['photons_per_second'] = adjusted_transit_curve.relative_flux * expected_photons_per_second # Changing y axis to photons per second
    return adjusted_transit_curve

# Used for the formula version of transit curve. Returns x^n if x>=0, -(-x)^n if x<0 def real_power(x, exponent):     if x >= 0: return x**exponent
    else: return -(-x)**exponent
def curve_formula(x, relative_flux_range, transit_length_days):
    time_centre_of_transit = transit_length_days / 2
    cos_factor = (1.06 / transit_length_days) * 4.2 # This is just variable used in the equation for the curve
    # After fiddling around the below formula seems to fit for during the transit (but not outside)
    return 1.0 - relative_flux_range*(real_power(np.cos(cos_factor*(x-time_centre_of_transit)), 0.2))

# transit_length (in seconds) actually slightly longer than actual transit - includes some flat data either side
# expected_photons_per_second is for when no transit taking place
# relative_flux_range is the amount we want the curve to dip

def get_transit_curve(transit_length, expected_photons_per_second, relative_flux_range):
    transit_length_days = transit_length / seconds_in_day
    time = np.arange(0.0, transit_length_days, 0.001) # List of times to fill with data points for curve
    relative_flux_curve = [curve_formula(x, relative_flux_range, transit_length_days) for x in time]
    transit_curve = pd.DataFrame({'time': time, 'relative_flux': relative_flux_curve})
    # Set all points outside of the transit to 1.0
    transit_curve.loc[transit_curve['relative_flux'] > 1, 'relative_flux'] = 1.0
    return change_transit_curve_units(transit_curve, expected_photons_per_second)

### 2. Simulating observing a transit ###

def get_simulated_data(transit_curve, standard_deviation, observation_frequency=30, exposure_time=20):
    observation_times = list(range(0, transit_curve.shape[0]-exposure_time, observation_frequency))
    expected_photons, simulated_photons = [], []
    coeff_of_variation = standard_deviation / (max(transit_curve.photons_per_second)*exposure_time)

    for observation_time in observation_times:
        end_observation_time = observation_time + exposure_time
        curve_rows = transit_curve.loc[(transit_curve['time_seconds'] >= observation_time) & (transit_curve['time_seconds'] < end_observation_time)]         expected_photons_observed = int(np.sum(curve_rows.photons_per_second))         simulated_photons_observed = round(np.random.normal(expected_photons_observed, expected_photons_observed*coeff_of_variation))         expected_photons.append(expected_photons_observed)         simulated_photons.append(simulated_photons_observed)              return pd.DataFrame({         'observation_time': observation_times,         'expected_photons': expected_photons,         'simulated_photons': simulated_photons     }) test_curve = get_transit_curve(test_transit_time, test_expected_photons_per_second, test_relative_flux_range) test_simulated_data = get_simulated_data(test_curve, test_standard_deviation) test_simulated_data.plot.scatter(x='observation_time', y='simulated_photons', s=1) ### 3. Chi-square testing ### ## 3.1 The functions for chi-square stat and generating curves of varying amplitude ## # Get chi-square stat def get_chi_square(data, transit_curve, standard_deviation, observation_frequency=30, exposure_time=20):     observation_times = data.observation_time          expected_photon_counts_from_curve = []       for observation_time in observation_times:         end_observation_time = observation_time + exposure_time         curve_rows = transit_curve.loc[(transit_curve['time_seconds'] >= observation_time) & (transit_curve['time_seconds'] < end_observation_time)]
    simulated_photon_counts = data.simulated_photons
    variance = standard_deviation**2
    chi_square = np.sum(np.square(simulated_photon_counts - expected_photon_counts_from_curve)) / variance
    return chi_square

# Generate list of transit curves, all paramters kept the same except amplitude (demonstrated below)
def get_curves_to_compare_against(transit_length, expected_photons_per_second, expected_relative_flux_range, number_of_curves=40):
    curves = []
    for flux_range in np.arange(0, expected_relative_flux_range*2, (expected_relative_flux_range)*2/number_of_curves):
        curves.append(get_transit_curve(transit_length, expected_photons_per_second, flux_range))
    return curves

# Get list chi squares for simulated data tested against varying curves
def get_multiple_chi_squares(data, transit_length, expected_photons_per_second, relative_flux_range, standard_deviation, number_of_curves=40):
    curves = get_curves_to_compare_against(transit_length, expected_photons_per_second, relative_flux_range)
    chi_square = []
    amplitude = []
    for curve in curves:
        chi_square.append(get_chi_square(data, curve, standard_deviation))
        amplitude.append(max(curve.relative_flux) - min(curve.relative_flux))
    return pd.DataFrame({'chi_square': chi_square, 'amplitude': amplitude})

## 3.2 Testing these out with the paramaters used for the actual transit in file B ##

# A demonstration showing all the curves I am comparing the simulated data against
curves = get_curves_to_compare_against(test_transit_time, test_expected_photons_per_second, test_relative_flux_range)
for curve in curves: plt.plot(curve.time_seconds, curve.relative_flux)
test_chi_square_df = get_multiple_chi_squares(test_simulated_data, test_transit_time, test_expected_photons_per_second, test_relative_flux_range, test_standard_deviation)

test_chi_square_df.plot.scatter(x='amplitude', y='chi_square')
plt.xlim([0, 0.011])
plt.axvline(x=test_relative_flux_range, color='orange')

## 3.2 Testing with more realistic parameters ##

# Similar to my previous simulated data but with much fewer points due to less time to observe
transit_time = 10000 # Rounded up to account for the flat parts of the curve at either end
relative_flux_range = 0.0025
transit_curve = get_transit_curve(transit_time, test_expected_photons_per_second, relative_flux_range)
simulated_data = get_simulated_data(transit_curve, test_standard_deviation)
simulated_data.plot.scatter(x='observation_time', y='simulated_photons', s=1)

chi_square_df = get_multiple_chi_squares(simulated_data, transit_time, test_expected_photons_per_second, relative_flux_range, test_standard_deviation)
chi_square_df.plot.scatter(x='amplitude', y='chi_square')
plt.xlim([0, 0.007])
plt.axvline(x=relative_flux_range, color='orange')

## 3.3 Testing with worse standard deviation ##

standard_dev = 3*test_expected_photons_per_second*20
x_simulated_data = get_simulated_data(transit_curve, standard_dev)
x_simulated_data.plot.scatter(x='observation_time', y='simulated_photons', s=1)

chi_square_df = get_multiple_chi_squares(x_simulated_data, transit_time, test_expected_photons_per_second, relative_flux_range, standard_dev)
chi_square_df.plot.scatter(x='amplitude', y='chi_square')
plt.xlim([0, 0.007])
plt.axvline(x=relative_flux_range, color='orange')

Telescope Project 2019

Hello! We are the telescope team of 2019! We are Aman, Paul, Victor, Ki, Tim, Mohyuddin and of course, Malcolm!

The purpose of this blog will be to keep a log of what we’ve been up to on a regular basis, and hopefully form a nice basis for the dreaded report.

Over the past few weeks we have been working on creating a Hertzsprung-Russell diagram of the Messier 47 star cluster. A HR diagram is a graph that plots stars according to their luminosity and colour. Using this diagram we are able to identify a star’s type depending on its position on the diagram.

Even more, by looking at the turn off point (the point where main sequence stars become red giants) we are able to figure out the cluster’s age.

In order to do this we first needed to get a number of images to work with.  In the observatory we managed to take 20 pictures of Messier 47. 10 of these were taken with a blue filter and the other 10 with a green filter.

The reason for taking these two sets of images are to allow us to calculate a B-V colour index, fulfilling one of the axis requirements for the HR diagram.

The next stage was for us to take the images and merge them into one. But first we had to get rid of ‘hot pixels’. These are bright points in our images caused by defects in the cameras system. Below is an example image – if you look closely you can very bright and very small dots. These are hot pixels, not stars!

An image we have taken of Messier 47. If you look closely enough you can see the tiny, but bright, hot pixels.

We first attempted processing the images by hand.  We removed the hot pixels by painting them out in GIMP, an image manipulation program. We then layered all the corresponding B and G images on top of each other and merged them using the ‘addition’ blending channel. Unfortunately, since the telescope was not perfectly still when it was taking the photos, you could see a significant directional blur. To remedy this we had to align all the images manually by eye. Below is what resulted.

This is the result of aligning and merging all 10 blue filter images of M47 in GIMP

Same as above but with the green filter instead.

While it was sufficient for our purposes, we decided that this approach was far too labour intensive. If we managed to go up again a take even more pictures, we’d have to repeat the whole tedious process. Instead it would be better to handle all this using the magic of python!

By working with the FITS files through python we hoped to automate the whole process.

We first converted the files into arrays and then took the logarithms of the intensity values. This gave the images a better contrast allowing us to see far more stars.

The raw FITS file

The same file but with a logarithm applied

However after taking the logarithm we were left with images with significant noise. The next step was to set a threshold intensity. Any value below this threshold would be changed to 0, and as a result we should be left with only stars and some hot pixels. We determined the optimum threshold value by taking the maximum of the second derivative of the total intensity.

Intensity values within a FITS file of M47

A graph of the 2nd derivative

For this image we found the optimum threshold intensity value to be 3.38.

The determined threshold applied to the log10 FITS file.

Since what we now had was just an array of points with different intensities, the next step was to cluster them into their own star arrays. First we used a Gaussian blur to smooth out the objects. This also had the additional benefit of blurring away the hot pixels.


Here is a Boolean test demonstrating the thresholding applied to the image. Notice all the noise and jagged shape of the stars.

But with the Gaussian blur applied we have a cleaner image with much smoother objects.

Here we have an even larger blurring radius applied.

The next step was to cluster pixels into their own groups of stars. This was accomplished with the ndimage.label function.

Here are the segregated objects. Each star has its own unique colour, indicating that it’s has it’s own distinct label.

Now that we had separated the intensities into their own stars, we could then work out the centre of masses for each one.

Centres of masses applied to the clusters

Centre of masses isolated

With the pre-processing complete, we could now work on automating the alignment of all the images with respect to each other.