3: Writing and including custom Effects

In this tutorial, we will load the model of MICADO (including Armazones, ELT, MORFEO) and then turn off all effect that modify the spatial extent of the stars. The purpose here is to see in detail what happens to the distribution of the stars flux on a sub-pixel level when we add a plug-in astrometric Effect to the optical system.

For real simulation, we will obviously leave all normal MICADO effects turned on, while still adding the plug-in Effect. Hopefully this tutorial will serve as a refernce for those who want to see how to create Plug-ins and how to manipulate the effects in the MICADO optical train model.

Create and optical model for MICADO and the ELT

from tempfile import TemporaryDirectory

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm

import scopesim as sim
from scopesim_templates.stellar import stars, star_grid

Scopesim works by using so-called instrument packages, which have to be downloaded separately. For normal use, you would set the package directory (a local folder path, local_package_folder in this example), download the required packages once, and then remove the download command.

local_package_folder = "./inst_pkgs"

However, to be able to run this example on the Readthedocs page, we need to include a temporary directory.

Do not copy and run this code locally, it is only needed to set things up for Readthedocs!

from tempfile import TemporaryDirectory
local_package_folder = TemporaryDirectory().name

Download the required instrument packages for an observation with MICADO at the ELT.

Again, you would only need to do this once, not every time you run the rest of the script, assuming you set a (permanent) instrument package folder.

sim.download_packages(["Armazones", "ELT", "MICADO", "MORFEO"])

astar.scopesim.server.database - Gathering information from server ...
astar.scopesim.server.database - Connection successful, starting download ...

Downloading Armazones: 100%|██████████| 74.3k/74.3k [00:00<00:00, 639kB/s]
Extracting  Armazones: 100%|██████████| 10/10 [00:00<00:00, 3171.74it/s]
Downloading ELT      : 100%|██████████| 53.4k/53.4k [00:00<00:00, 25.8MB/s]
Extracting  ELT      : 100%|██████████| 16/16 [00:00<00:00, 5812.81it/s]
Downloading MICADO   : 100%|██████████| 14.4M/14.4M [00:27<00:00, 549kB/s]
Extracting  MICADO   : 100%|██████████| 121/121 [00:00<00:00, 843.72it/s]
Downloading MORFEO   : 100%|██████████| 1.38M/1.38M [00:02<00:00, 562kB/s]
Extracting  MORFEO   : 100%|██████████| 12/12 [00:00<00:00, 953.34it/s]

[PosixPath('/home/docs/checkouts/readthedocs.org/user_builds/scopesim/checkouts/latest/docs/source/examples/inst_pkgs/Armazones.zip'),
 PosixPath('/home/docs/checkouts/readthedocs.org/user_builds/scopesim/checkouts/latest/docs/source/examples/inst_pkgs/ELT.zip'),
 PosixPath('/home/docs/checkouts/readthedocs.org/user_builds/scopesim/checkouts/latest/docs/source/examples/inst_pkgs/MICADO.zip'),
 PosixPath('/home/docs/checkouts/readthedocs.org/user_builds/scopesim/checkouts/latest/docs/source/examples/inst_pkgs/MORFEO.zip')]

We can see which Effects are already included by calling micado.effects:

cmd = sim.UserCommands(use_instrument="MICADO", set_modes=["SCAO", "IMG_1.5mas"])
micado = sim.OpticalTrain(cmd)

micado.effects

Table length=23

element	name	class	included
str13	str23	str31	bool
armazones	skycalc_atmosphere	SkycalcTERCurve	True
ELT	telescope_reflection	SurfaceList	True
MICADO	micado_static_surfaces	SurfaceList	True
MICADO	micado_ncpas_psf	NonCommonPathAberration	True
MICADO	filter_wheel_1 : [open]	FilterWheel	True
MICADO	filter_wheel_2 : [Ks]	FilterWheel	True
MICADO	pupil_wheel : [open]	FilterWheel	True
MICADO_DET	full_detector_array	DetectorList	False
MICADO_DET	detector_window	DetectorWindow	True
MICADO_DET	qe_curve	QuantumEfficiencyCurve	True
...	...	...	...
MICADO_DET	detector_linearity	LinearityCurve	True
MICADO_DET	border_reference_pixels	ReferencePixelBorder	True
MICADO_DET	readout_noise	PoorMansHxRGReadoutNoise	True
MICADO_DET	source_fits_keywords	SourceDescriptionFitsKeywords	False
MICADO_DET	extra_fits_keywords	ExtraFitsKeywords	True
default_ro	relay_psf	FieldConstantPSF	True
default_ro	relay_surface_list	SurfaceList	True
default_ro	extra_fits_keywords_ro	ExtraFitsKeywords	True
MICADO_IMG_HR	zoom_mirror_list	SurfaceList	True
MICADO_IMG_HR	micado_adc_3D_shift	AtmosphericDispersionCorrection	False

Now we turn off all Effects that cause spatial aberrations:

for effect_name in ["full_detector_array", "micado_adc_3D_shift",
                    "micado_ncpas_psf", "relay_psf"]:
    micado[effect_name].include = False
    print(micado[effect_name])

DetectorList: "full_detector_array"
AtmosphericDispersionCorrection: "micado_adc_3D_shift"
NonCommonPathAberration: "micado_ncpas_psf"
FieldConstantPSF: "relay_psf"

The normal detector window is set to 1024 pixels square. Let’s reduce the size of the detector readout window:

micado["detector_window"].data["x_cen"] = 0          # [mm] distance from optical axis on the focal plane
micado["detector_window"].data["y_cen"] = 0
micado["detector_window"].data["x_size"] = 64          # [pixel] width of detector
micado["detector_window"].data["y_size"] = 64

By default ScopeSim works on the whole pixel level for saving computation time. However it is capable of integrating sub pixel shift. For this we need to turn on the sub-pixel mode:

micado.cmds["!SIM.sub_pixel.flag"] = True

We can test what’s happening by making a grid of stars and observing them:

src = star_grid(n=9, mmin=20, mmax=20.0001, separation=0.0015 * 15)
src.fields[0]["x"] -= 0.00075
src.fields[0]["y"] -= 0.00075

micado.observe(src, update=True)

plt.figure(figsize=(8,8))
plt.imshow(micado.image_planes[0].data, origin="lower")

 FOVs:   0%|          | 0/1 [00:00<?, ?it/s]
 FOV effects: 0it [00:00, ?it/s]
 FOVs: 100%|██████████| 1/1 [00:00<00:00, 13.62it/s]
 Image Plane effects: 100%|██████████| 1/1 [00:00<00:00, 5249.44it/s]

<matplotlib.image.AxesImage at 0x7f2f6012c790>

micado["detector_window"].data

Table length=1

id	x_cen	y_cen	x_size	y_size	angle	gain	pixel_size
int64	str6	str6	str10	str11	int64	int64	float64
0	0	0	64	64	0	1	0.015

Writing a custom Effect object

The following code snippet creates a new Effect class.

import numpy as np
from astropy.table import Table

from scopesim.effects import Effect
from scopesim.base_classes import SourceBase


class PointSourceJitter(Effect):
    def __init__(self, **kwargs):
        super(PointSourceJitter, self).__init__(**kwargs)   # initialise the underlying Effect class object
        self.meta["z_order"] = [500]                        # z_order number for knowing when and how to apply the Effect
        self.meta["max_jitter"] = 0.001                     # [arcsec] - a parameter needed by the effect
        self.meta.update(kwargs)                            # add any extra parameters passed when initialising

    def apply_to(self, obj):                                # the function that does the work
        if isinstance(obj, SourceBase):
            for field in obj.fields:
                if isinstance(field, Table):
                    dx, dy = 2 * (np.random.random(size=(2, len(field))) - 0.5)
                    field["x"] += dx * self.meta["max_jitter"]
                    field["y"] += dy * self.meta["max_jitter"]

        return obj

Lets break it down a bit (THIS IS JUST A STEP-BY-STEP EXPLANATION OF THE CODE ABOVE, NOT SOMETHING NEW!):

class PointSourceJitter(Effect):
    ...

Here we are subclassing the Effect object from ScopeSim. This has the basic functionality for reading in ASCII and FITS files, and for communicating with the OpticsManager class in ScopeSim.

The initialisation function looks like this:

def __init__(self, **kwargs):
    super(PointSourceJitter, self).__init__(**kwargs)   # initialise the underlying Effect class object
    self.meta["z_order"] = [500]

Here we make sure to activate the underlying Effect object. The z_order keyword in the meta dictionary is used by ScopeSim to determine when and where this Effect should be applied during a simulations run. The exact z-order numbers are described in [insert link here].

The main function of any Effect is the apply_to method:

def apply_to(self, obj):
    if isinstance(obj, SourceBase):
        ...

    return obj

It should be noted that what is passed in via (obj) must be returned in the same format. The contents of the obj can change, but the obj object must be returned.

All the code which enacts the results of the physical effect are contained in this method. For example, if we are writing a redshifting Effect, we could write the code to shift the wavelength array of a Source object by z+1 here.

There are 4 main classes that are cycled through during an observation run: * SourceBase: contains the original 2+1D distribtion of light, * FieldOfViewBase: contains a (quasi-)monochromatic cutout from the Source object, * ImagePlaneBase: contains the expectation flux image on the detector plane * DetectorBase: contains the electronic readout image

An Effect object can be applied to any number of objects based on one or more of these base classes. Just remember to segregate the base-class-specific code with if statements.

One further method should be mentioned: def fov_grid(). This method is used by FOVManager to estimate how many FieldOfView objects to generate in order to best simulation the observation. If your Effect object might alter this estimate, then you should include this method in your class. See the code base for further details.

Note: The fov_grid method will be depreciated in a future release of ScopeSim. It will most likely be replaced by a FOVSetupBase class that will be cycled through the apply_to function. However this is not yet 100% certain, so please bear with us.

Including a custom Effect

First we need to initialise an instance of the Effect object:

jitter_effect = PointSourceJitter(max_jitter=0.001, name="random_jitter")

Then we can add it to the optical model:

micado.optics_manager.add_effect(jitter_effect)

micado.effects

Table length=24

element	name	class	included
str13	str23	str31	bool
armazones	skycalc_atmosphere	SkycalcTERCurve	True
armazones	random_jitter	PointSourceJitter	True
ELT	telescope_reflection	SurfaceList	True
MICADO	micado_static_surfaces	SurfaceList	True
MICADO	micado_ncpas_psf	NonCommonPathAberration	False
MICADO	filter_wheel_1 : [open]	FilterWheel	True
MICADO	filter_wheel_2 : [Ks]	FilterWheel	True
MICADO	pupil_wheel : [open]	FilterWheel	True
MICADO_DET	full_detector_array	DetectorList	False
MICADO_DET	detector_window	DetectorWindow	True
...	...	...	...
MICADO_DET	detector_linearity	LinearityCurve	True
MICADO_DET	border_reference_pixels	ReferencePixelBorder	True
MICADO_DET	readout_noise	PoorMansHxRGReadoutNoise	True
MICADO_DET	source_fits_keywords	SourceDescriptionFitsKeywords	False
MICADO_DET	extra_fits_keywords	ExtraFitsKeywords	True
default_ro	relay_psf	FieldConstantPSF	False
default_ro	relay_surface_list	SurfaceList	True
default_ro	extra_fits_keywords_ro	ExtraFitsKeywords	True
MICADO_IMG_HR	zoom_mirror_list	SurfaceList	True
MICADO_IMG_HR	micado_adc_3D_shift	AtmosphericDispersionCorrection	False

When we want to observe, we need to include the update=True flag so that the optical model is updated to include the instance of our new Effect:

micado.observe(src, update=True)

plt.figure(figsize=(8,8))
plt.imshow(micado.image_planes[0].data, origin="lower")

 FOVs:   0%|          | 0/1 [00:00<?, ?it/s]
 FOV effects: 0it [00:00, ?it/s]
 FOVs: 100%|██████████| 1/1 [00:00<00:00, 13.62it/s]
 Image Plane effects: 100%|██████████| 1/1 [00:00<00:00, 4655.17it/s]

<matplotlib.image.AxesImage at 0x7f2f5f470f10>

We can update the parameters of the object on-the-fly by accessing the meta dictionary:

micado["random_jitter"].meta["max_jitter"] = 0.005

micado.observe(src, update=True)

plt.figure(figsize=(8,8))
plt.imshow(micado.image_planes[0].data, origin="lower")

 FOVs:   0%|          | 0/1 [00:00<?, ?it/s]
 FOV effects: 0it [00:00, ?it/s]
 FOVs: 100%|██████████| 1/1 [00:00<00:00, 13.66it/s]
 Image Plane effects: 100%|██████████| 1/1 [00:00<00:00, 5236.33it/s]

<matplotlib.image.AxesImage at 0x7f2f5f35c650>

Here we can see that there is a certain amount of sub-pixel jitter being introduced into each observation. However this bare-bones approach is not very realistic. We should therefore turn the PSF back on to get a more realistic observation:

micado["relay_psf"].include = True

micado.observe(src, update=True)
hdus = micado.readout()

plt.figure(figsize=(8,8))
plt.imshow(hdus[0][1].data, origin="lower", norm=LogNorm())

 FOVs:   0%|          | 0/1 [00:00<?, ?it/s]
 FOV effects:   0%|          | 0/1 [00:00<?, ?it/s]

astar.scopesim.effects.psf_utils - WARNING: PSF center off

 FOVs: 100%|██████████| 1/1 [00:00<00:00, 10.29it/s]
 Image Plane effects: 100%|██████████| 1/1 [00:00<00:00, 5210.32it/s]

astar.scopesim.detector.detector_array - Extracting from 1 detectors...

<matplotlib.image.AxesImage at 0x7f2f5f37d6d0>