# Lens Modeling: How to
## 0. Handling FITS Files and Generating Needed Data
### FITS Format
`.fits` is the most used digital file format in astronomy. These files are opened with several different software, but we'll use **[SAOImageDS9](https://sites.google.com/cfa.harvard.edu/saoimageds9)**. The main functionalities we'll need are the viewing of data and the creation of regions. See: [DS9Guide](https://www.jb.man.ac.uk/~gbendo/Sci/Pict/DS9guide.pdf)
A `.fits` file consists of 2 main parts, the `header` and the `data`.
A typical `header` looks like this
```
IMPLE = T / conforms to FITS standard
BITPIX = -32 / array data type
NAXIS = 2 / number of array dimensions
NAXIS1 = 2362
NAXIS2 = 2423
EXTEND = T
COMMENT FITS (Flexible Image Transport System) format is defined in 'Astronomy
COMMENT and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H
FILENAME= '../reduced_data/WG0214-2105_F160W_drz_sci.fits' / name of file
FILETYPE= 'SCI ' / type of data found in data file
TELESCOP= 'HST' / telescope used to acquire data
INSTRUME= 'WFC3 ' / identifier for instrument used to acquire data
EQUINOX = 2000.0 / equinox of celestial coord. system
/ DATA DESCRIPTION KEYWORDS
ROOTNAME= '../reduced_data/WG0214-2105_F160W' / rootname of the observation set
IMAGETYP= 'EXT ' / type of exposure identifier
PRIMESI = 'WFC3 ' / instrument designated as prime
/ TARGET INFORMATION
TARGNAME= 'WG0214-2105 ' / proposer's target name
RA_TARG = 3.356812500000E+01 / right ascension of the target (deg) (J2000)
DEC_TARG= -2.109304722222E+01 / declination of the target (deg) (J2000)
/ PROPOSAL INFORMATION
PROPOSID= 15652 / PEP proposal identifier
LINENUM = 'Z3.009 ' / proposal logsheet line number
PR_INV_L= 'Treu ' / last name of principal investigator
PR_INV_F= 'Tommaso ' / first name of principal investigator
PR_INV_M= 'L. ' / middle name / initial of principal investigator
[...]
[...]
[...]
```
It includes the metadata corresponding to the observation done. Among others you can find useful information about exposure time, pixel size, WCS, etc.
On the other side, the `data` is multidimensional array, where each entry corresponds to a pixel In essence, the array
```
data = np.array([[1, 0, 0], [0, 1, 0], [0,0,1]])
```
represents the following image:
<center>
</center>
### Cutout
Observation data includes a lot of information which is not needed, so in order to model we would like to have a cutout of the system we're trying to model.
In order to do this we use the `Cutout2D` function from `astropy`
```
class Cutout2D(
data: data,
position: tuple[int, int],
size: int,
wcs: Any | None = None,
mode: str = "trim",
fill_value: float = np.nan,
copy: bool = False
)
```
where `data` is the data from the initial `.fits` file. Note: It is good practice to update the header (particularly the WCS) accordingly. An example of the script to do a cutout can be seen here:
```
import astropy
from astropy.io import fits
import matplotlib.pyplot as plt
from astropy.visualization import astropy_mpl_style
from astropy.wcs import WCS
from astropy.nddata import Cutout2D
import tkinter, tkinter.filedialog
import matplotlib.pyplot as plt
tkinter.Tk().withdraw()
# Open the SCI file
print("Select the SCI File to be cropped")
file_path = tkinter.filedialog.askopenfilename()
hdulist = fits.open(file_path)
#Initialize variables for header, data and WCS
header = hdulist['PRIMARY'].header
data = hdulist['PRIMARY'].data
wcs = WCS(header)
# Read the coordinates for the center of the cutout
c0x =int(input('Cutout center, x-axis, pixels: \n'))
c0y =int(input('Cutout center, y-axis, pixels: \n'))
center = (c0x, c0y)
# Read size of cutout
size=int(input('Size of the (square) cutout, pixels: \n'))
#Use Cutout2D to generate a cutout
cutout = Cutout2D(data, center, size, wcs=wcs)
#Keep the header but modified for the cutout - this is good practice
cheader = cutout.wcs.to_header()
primaryhdu = fits.PrimaryHDU(cutout.data, cheader)
hdulist = fits.HDUList([primaryhdu])
#plot the file
plt.figure()
plt.style.use(astropy_mpl_style)
plt.imshow(cutout.data, cmap='gray',origin= 'lower')
plt.colorbar()
plt.show()
# Save the file
hdulist.writeto('Cutout.fits', overwrite=True)
```
### Background Reduction
Similarly, we'll use `astropy` to reduce the background of the image. We will subtract a floor that leaves a very small (close to zero, they can be negative!) value in every pixel. To obtain the floor value we make a region in ds9 on an empty patch. In the region statistics we look at the computed median and use this as our floor.
```
import astropy
import numpy as np
from astropy.io import fits
import matplotlib.pyplot as plt
from astropy.visualization import astropy_mpl_style
from astropy.wcs import WCS
from astropy.nddata import Cutout2D
import tkinter, tkinter.filedialog
import matplotlib.pyplot as plt
tkinter.Tk().withdraw()
# Open the SCI file
print("Select the SCI File to reduce background from")
file_path = tkinter.filedialog.askopenfilename()
hdulist = fits.open(file_path)
#Initialize variables for header, data and WCS
header = hdulist['PRIMARY'].header
data = hdulist['PRIMARY'].data
wcs = WCS(header)
# Reduce the value in each pixel
bg= float(input("Input the background value: \n"))
newdata = data - bg
#Keep the header
primaryhdu = fits.PrimaryHDU(newdata, header)
#Create fits file
hdulist = fits.HDUList([primaryhdu])
# Save the file
hdulist.writeto('Noisereduced.fits', overwrite=True)
```
### PSF
The Point Spread Function (PSF) is the geometric reaction an optical system has to a point source. In Hubble it is characterized by the diffraction spikes on star. In JWST, newer, it is characterized by the diamont shape that stars have.
It affects the observation of the system, as spikes and *airy* rings are visible around the quasars. This airy rings are not to be confused with the arc generated by the gravitational lenses. The latter will be discused further down this guide.
To generate the PSF please refer to ``starred``, a python package to generate psf from stars in the sci file by Michalewicz et al.
Files needed for the reconstruction of the PSF are:
- Cutout of the stars in the SCI File. (Make them all the same size).
- Patching the cutout of the stars, using a random normal distribution with a mean and RMS obtained from the BGR Sci file. Make sure the BGR is used.
- Error Map of each star. A basic $\sigma$ map of the star cutout. Use the weight files for this.
In the case of WG0214 in IR band (160), four stars with a decent amount of spikes/rings were used.
After feeding these images, error maps, and extra model information, the final PSF reconstruction looks like this:
In which the Airy rings around the point source are quite visible. The spikes are also visible if the contrast is adjusted.
**Important:** A PSF generated for one set of observations is NOT compatible with another set of obs. The PSF is wavelenght and construction dependent. Since Hubble is in LEO, the effects of the earth atmosphere, solar induced contraction and expansion, this varies. It also varies in other optical systems, such as the ELT, VLA, JWST, etc. It is always better to obtain the psf from the observed data rather than other sources.
[GitHub Repo](https://gitlab.com/cosmograil/starred)
### Masks
We initially need to do 3 masks. We first make the region files for this.
#### 1-Pixel mask
This is just a mask with all but one single corner pixel set to 0.
#### Arc Mask
This one covers both, the arc and the light from the source galaxy - leaving the lens galaxy uncovered. Make sure the region file is saved in the image coordinate system (pixels), DS9 format, and double check there is no other figures in the file.
<center>
 
</center>
#### Object Mask
This one covers all the light objects not belonging to the system to be modelled. Save each region individually. Make sure each region file is saved with the image coordinate system (pixels), DS9 format and double check there is no other figures in the file.
<center>
 
</center>
Now, we use the script to generate the fits files for the masks.
### Error Map
For this we also need a region, which is more or less "empty", so that we can get a sample of the noise in the whole image.
<center>
 
</center>
Then we use this region (Once again, make sure the region file is saved in the image coordinate system, DS9 format, and double check there is no other figures in the file) and the errormap script to generate an errormap.
<center>
</center>
## 1. Mass Modelling
### The configfile
We begin with a `configfile` this text file (extensionless) contains the data of the model. Our initial version looks like this:
```
chi2type 1
minimiser siman
seed 1
siman_iter 1
siman_nT 1000
siman_dS 10
siman_Sf 1.1
siman_k 2
siman_Ti 10
siman_Tf 1.1
siman_Tmin 1
lenses_vary 2
spemd
4.80 #x-coord flat:4.7,4.9 step:0.01
4.80 #y-coord flat:4.7,4.9 step:0.01
0.90 #b/a flat:0.8,1 step:0.01
0.0 #theta flat:0,6.283: step:0.01
0.828 #theta_e flat:0.5,1.2 scale: step:0.01
0.0000001E-4 #r_core exact:
0.500522 #gam flat:0.3,0.7 step:0.01
shear
0.0000 flat:0,0.3 step:0.01 #magnitude
0.0000 noprior: step:0.01 #theta
sources 1
Dds/Ds 1.000000 exact:
source weighted
srcx 4.82 exact:
srcy 4.90 exact:
4
4.3200000 3.9200000 error:0.04
5.6000000 4.7200000 error:0.04
5.1200000 5.6000000 error:0.04
4.24000000 5.3600000 error:0.04
```
From top to bottom:
This first part
```
chi2type 1
minimiser siman
seed 1
siman_iter 1
siman_nT 1000
siman_dS 10
siman_Sf 1.1
siman_k 2
siman_Ti 10
siman_Tf 1.1
siman_Tmin 1
```
contains information on how the modelling is gonna take place. Since we're initially modelling first the `chi2type` we'll use is 1, which corresponds to the source plane positions. We'll be optimizing that first. We'll be initially the using simulated annealing (siman) agorithm with the parameters, Temperature, Step Size, etc., as defined above.
Now, this part
```
lenses_vary 2
spemd
4.80 #x-coord flat:4.7,4.9 step:0.01
4.80 #y-coord flat:4.7,4.9 step:0.01
0.90 #b/a flat:0.8,1 step:0.01
0.0 #theta flat:0,6.283: step:0.01
0.828 #theta_e flat:0.5,1.2 scale: step:0.01
0.0000001E-4 #r_core exact:
0.500522 #gam flat:0.3,0.7 step:0.01
shear
0.0000 flat:0,0.3 step:0.01 #magnitude
0.0000 noprior: step:0.01 #theta
```
defines the initial values for the lens galaxy. The mass model assumed in the above case is a single power elliptical mass distribution (SPEMD). `x/y-coord` are the coordinates (in arcseconds) of the center of the lens, `b/a` is the eccentricity of the galaxy, `theta` is the angle of the galaxy measured in radians from the +x-axis, `theta_e` is the einstein radius in `arcsecongd`, `r_core` is the radius of the galaxy core. Usually a very small number., lastly `gamm` is a factor in the mass distribution (see literature).
Shear is the distortion caused due to the large scale structures of the universe (0.1% to 1%) and causes some minor effects, which are factored into the model by defining it.
We define the priors to be the interval inbetween which the actual value will be in, not everything has to have a prior, but it is easier with one. `flat:` means the probability of every value is the same, but other options like `gaussian:` are possible (see glee wiki).
Lastly, the final part of the initial config file
```
sources 1
Dds/Ds 1.000000 exact:
source weighted
srcx 4.82 exact:
srcy 4.90 exact:
4
4.3200000 3.9200000 error:0.04
5.6000000 4.7200000 error:0.04
5.1200000 5.6000000 error:0.04
4.24000000 5.3600000 error:0.04
```
This correspond to the sources. `srcx` and `srcy` are the weighted averages if the position of the images, followed by the number of images (this case 4) and their position.
### Source Plane Position Optimization (chi2type=1)
We run `glee -m -i configfile` and this returns a `configfile.001` file. (For organization reasons, I'll refer to it as `configfile01`). This new updated file contains the first optimized version of the file. On this new file, you can run `glee -c configfile01` to see the current `chi2` value. This should be samll (less than TBC)
### Image Plane Position Optimization (chi2type=2)
We change the `chi2type` to 2 in `configfile01` and rerun an annealing (`glee -m -i configfile1`). This returns `configfile02`. Again, check that the chi2 is below the value.
### Visualization
We can see the results of this optimization running the `gleeview.py` tool. Before running this, we need to generate the `.fits` cube with the required data (caustics, critical curves, contours, ...), we do this with the `glee -s configfile02`. The result should look like this. The sources (green diamonds) should be bundled up together. The caustics are the red lines, the critical curve is the blue line, and the images are the circle. The center of the lens is the blue square.
## 2. Lens Plane Modelling
### Lens Light
To begin modeling the light of the lens, we have to modify the mass-optimized `configfile`. We first set the `chi2type` to 16. Then we set all the mass parameters to `exact:` (leave the `scale:` in `theta_e`). Now we add the sources and we do so by adding the `esource` (extended sources) section. We also add a new algorithm, the MCMC Chains. The file now looks like:
```
chi2type 16
minimiser siman
seed 1
siman_iter 1
siman_nT 1000
siman_dS 10
siman_Sf 1.1
siman_k 2
siman_Ti 10
siman_Tf 1.1
siman_Tmin 1
mcmc_n 100000
mcmc_dS 0.04375
mcmc_dSini 0
mcmc_k 2
lenses_vary 2
spemd
4.794076 #x-coord exact:
4.785908 #y-coord exact:
0.809780 #b/a exact:
0.105977 #theta exact:
0.881267 #theta_e exact: scale:
1.0000e-11 #r_core exact:
0.540466 #gam exact:
shear
0.083024 #magnitude exact:
0.744317 #theta exact:
sources 1
Dds/Ds 1.000000 exact:
source weighted
srcx 4.335643 exact:
srcy 4.377210 exact:
4
4.320000 3.920000 error:0.04
5.600000 4.720000 error:0.04
5.120000 5.600000 error:0.04
4.240000 5.360000 error:0.04
esources 1
Dds/Ds 1.000000 exact:
esource_ngy 20
esource_ngx 120
esource_dx 0.080000
esource_data /afs/mpa/home/allansch/Lens/WG0214/WG0214_F160W_sci_CO_NR.fits
esource_err /afs/mpa/home/allansch/Lens/WG0214/ErrorMap.fits
esource_arcmask /afs/mpa/home/allansch/Lens/WG0214/Mask_Arc1.fits
esource_lensmask /afs/mpa/home/allansch/Lens/WG0214/Mask_Lens1.fits
esource_mod_light LensOnly
esource_psf /afs/mpa/home/allansch/Lens/WG0214/PSF.fits
esource_sub_agn_psf /afs/mpa/home/allansch/Lens/WG0214/PSF_agn.fits
esource_sub_agn_psf_factor 3
esource_sub_esr_psf /afs/mpa/home/allansch/Lens/WG0214/PSF_esr.fits
esource_sub_esr_psf_factor 3
esource_regopt SpecRegPrecSigFigOnce
esource_reglampre 1
esource_reglamnup 1000
esource_regtype curv
esource_reglam 1
esource_reglamlo 1e-05
esource_reglamhi 100000
esource_light 2
sersic
4.800000 #x-coord flat:4.72,4.88 label:sersic_x step:0.01
4.800000 #y-coord flat:4.72,4.88 label:sersic_y step:0.01
0.800000 #q flat:0.6,1 step:0.01
3.000000 #PA noprior: step:0.01
0.300000 #amp noprior: min:0 step:0.01
0.900000 #r_eff flat:0,1.5 step:0.01
0.700000 #n_sersic flat:0.5,9 step:0.01
sersic
4.810615 #x-coord link:sersic_x a:0,1,1
4.791608 #y-coord link:sersic_y a:0,1,1
0.800000 #q flat:0.6,1 step:0.01
3.000000 #PA noprior: step:0.01
0.300000 #amp noprior: min:0 step:0.01
0.900000 #r_eff flat:0,1.5 step:0.01
0.700000 #n_sersic flat:0.5,9 step:0.01
esource_end
```
Specific details on what everything means in the GLEE wiki. Here we’re modelling the lens galaxy using 2 sersic profiles. Notice how the centroids of these profiles are linked to each other. The triplet following `a: a1,a2,a3` means the value will be transformed according to: $y = a_1+a_2\cdot x^{a_3}$.
Once everything is set up, you run one siman to jump into parameter space. This is followed by several MCMC chains (done with `glee -M -i configfile0x`). Pay special care that the acceptance has a value of 0.25±0.05. You can adjust this by changing the value of `mcmc_dS`. After each chain (generally 100k-300k iterations), save the covariance matrix using the command `glee -h configfile0x.mcmc` and add it to the next `configfile` by adding the lines
```
sampling_f gaussian
sampling_cov /afs/mpa/home/allansch/Lens/WG0214/WG0214_0x.cov
```
after the `mcmc_k` line. Then repeat.
To visualize the returned file, `configfile0x.mcmc` we use `gleeviewsample.py`. here we can see the sampling and it's `log_evidence` value, which tells us how good the fit is and also gives us a corner plot of how the parameters vary with each other.
<center>
</center>
after several chains:
Furthermore, by running `glee -f 2 configfile0x` we can export `fits` cube that contains (1) the original image, (2) the model, (3) residuals, (4) normalized residuals.
### Lens Light + Quasar Light
Now we set all the parameters regarding the light of the lens to `exact:`, we add the number of images (in this case 4) to the number of profiles in `esource_light`. Also, we now remove the `arcmask` and replace it with a mask of the exact same dimensions, but in which the data is just zeros except for one pixel in any corner, set to 1. This because glee requires a mask, so we put in this minimally invasive mask to make up for it.
The file now looks like
```
chi2type 16
minimiser siman
seed 1
siman_iter 1
siman_nT 1000
siman_dS 10
siman_Sf 1.1
siman_k 2
siman_Ti 10
siman_Tf 1.1
siman_Tmin 1
mcmc_n 100000
mcmc_dS 0.3
mcmc_dSini 0
mcmc_k 2
lenses_vary 2
spemd
4.794076 #x-coord exact:
4.785908 #y-coord exact:
0.809780 #b/a exact:
0.105977 #theta exact:
0.881267 #theta_e exact: scale:
1.0000e-11 #r_core exact:
0.540466 #gam exact:
shear
0.083024 #magnitude exact:
0.744317 #theta exact:
sources 1
Dds/Ds 1.000000 exact:
source weighted
srcx 4.335643 exact:
srcy 4.377210 exact:
4
4.320000 3.920000 error:0.04
5.600000 4.720000 error:0.04
5.120000 5.600000 error:0.04
4.240000 5.360000 error:0.04
esources 1
Dds/Ds 1.000000 exact:
esource_ngy 20
esource_ngx 120
esource_dx 0.080000
esource_data /afs/mpa/home/allansch/Lens/WG0214/WG0214_F160W_sci_CO_NR2.fits
esource_err /afs/mpa/home/allansch/Lens/WG0214/ErrorMap.fits
esource_arcmask /afs/mpa/home/allansch/Lens/WG0214/Mask_1pix.fits
esource_lensmask /afs/mpa/home/allansch/Lens/WG0214/Mask_Lens1.fits
esource_mod_light LensOnly
esource_psf /afs/mpa/home/allansch/Lens/WG0214/PSF.fits
esource_sub_agn_psf /afs/mpa/home/allansch/Lens/WG0214/PSF_agn.fits
esource_sub_agn_psf_factor 3
esource_sub_esr_psf /afs/mpa/home/allansch/Lens/WG0214/PSF_esr.fits
esource_sub_esr_psf_factor 3
esource_regopt SpecRegPrecSigFigOnce
esource_reglampre 1
esource_reglamnup 1000
esource_regtype curv
esource_reglam 9.99288
esource_reglamlo 1e-05
esource_reglamhi 100000
esource_light 6
...
psf
5.205020 #x-coord noprior: step:0.01
5.495106 #y-coord noprior: step:0.01
170.869832 #amp noprior: min:0 step:2.61399
psf
5.580431 #x-coord noprior: step:0.01
4.644163 #y-coord noprior: step:0.01
225.230848 #amp noprior: min:0 step:2.88811
psf
4.386284 #x-coord noprior: step:0.01
3.922756 #y-coord noprior: step:0.01
351.470253 #amp noprior: min:0 step:2.63407
psf
4.193880 #x-coord noprior: step:0.01
5.309741 #y-coord noprior: step:0.01
149.845956 #amp noprior: min:0 step:1.54551
esource_end
```
Google Drive
Gist
Clipboard
Sign in with Wallet
