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"])
httpx - HTTP Request: GET https://scopesim.univie.ac.at/InstPkgSvr/example_data/ "HTTP/1.1 200 OK"
httpx - HTTP Request: GET https://scopesim.univie.ac.at/InstPkgSvr/instruments/ "HTTP/1.1 200 OK"
httpx - HTTP Request: GET https://scopesim.univie.ac.at/InstPkgSvr/locations/ "HTTP/1.1 200 OK"
httpx - HTTP Request: GET https://scopesim.univie.ac.at/InstPkgSvr/psfs/ "HTTP/1.1 200 OK"
httpx - HTTP Request: GET https://scopesim.univie.ac.at/InstPkgSvr/telescopes/ "HTTP/1.1 200 OK"
httpx - HTTP Request: GET https://scopesim.univie.ac.at/InstPkgSvr/instruments/ "HTTP/1.1 200 OK"
httpx - HTTP Request: GET https://scopesim.univie.ac.at/InstPkgSvr/locations/ "HTTP/1.1 200 OK"
httpx - HTTP Request: GET https://scopesim.univie.ac.at/InstPkgSvr/telescopes/ "HTTP/1.1 200 OK"
httpx - HTTP Request: GET https://scopesim.univie.ac.at/InstPkgSvr/ "HTTP/1.1 200 OK"
httpx - HTTP Request: GET https://scopesim.univie.ac.at/InstPkgSvr/example_data/ "HTTP/1.1 200 OK"
httpx - HTTP Request: GET https://scopesim.univie.ac.at/InstPkgSvr/instruments/ "HTTP/1.1 200 OK"
httpx - HTTP Request: GET https://scopesim.univie.ac.at/InstPkgSvr/locations/ "HTTP/1.1 200 OK"
httpx - HTTP Request: GET https://scopesim.univie.ac.at/InstPkgSvr/psfs/ "HTTP/1.1 200 OK"
httpx - HTTP Request: GET https://scopesim.univie.ac.at/InstPkgSvr/telescopes/ "HTTP/1.1 200 OK"
astar.scopesim.server.database - Connection successful, starting download ...
httpx - HTTP Request: GET https://scopesim.univie.ac.at/InstPkgSvr/locations/Armazones.2023-07-11.zip "HTTP/1.1 200 OK"
httpx - HTTP Request: GET https://scopesim.univie.ac.at/InstPkgSvr/telescopes/ELT.2024-02-29.zip "HTTP/1.1 200 OK"
httpx - HTTP Request: GET https://scopesim.univie.ac.at/InstPkgSvr/instruments/MICADO.2024-02-29.zip "HTTP/1.1 200 OK"
---------------------------------------------------------------------------
KeyboardInterrupt Traceback (most recent call last)
Cell In[4], line 1
----> 1 sim.download_packages(["Armazones", "ELT", "MICADO", "MORFEO"])
File ~/checkouts/readthedocs.org/user_builds/scopesim/checkouts/stable/scopesim/server/database.py:414, in download_packages(pkg_names, release, save_dir)
412 for pkg_name in pkg_names:
413 try:
--> 414 pkg_path = _download_single_package(
415 client, pkg_name, release, all_versions, folders_dict,
416 save_dir, padlen)
417 except PkgNotFoundError as error:
418 # Whole stack trace not useful for enduser.
419 # Could log it to file though...
420 # logger.exception(error)
421 logger.error(error)
File ~/checkouts/readthedocs.org/user_builds/scopesim/checkouts/stable/scopesim/server/database.py:333, in _download_single_package(client, pkg_name, release, all_versions, folders_dict, save_dir, padlen)
330 pkg_url = f"{folders_dict[pkg_name]}/{zip_name}"
332 save_path = save_dir / f"{pkg_name}.zip"
--> 333 handle_download(client, pkg_url, save_path, pkg_name, padlen)
334 handle_unzipping(save_path, save_dir, pkg_name, padlen)
336 return save_path.absolute()
File ~/checkouts/readthedocs.org/user_builds/scopesim/checkouts/stable/scopesim/server/download_utils.py:130, in handle_download(client, sub_url, save_path, name, padlen, chunk_size, disable_bar, params)
124 total = int(response.headers.get("Content-Length", 0))
126 with (save_path.open("wb") as file_outer,
127 tqdm.wrapattr(
128 file_outer, "write", miniters=1, total=total,
129 **tqdm_kwargs, disable=disable_bar) as file_inner):
--> 130 for chunk in response.iter_bytes(chunk_size=chunk_size):
131 file_inner.write(chunk)
133 except httpx.HTTPStatusError as err:
File ~/checkouts/readthedocs.org/user_builds/scopesim/envs/stable/lib/python3.11/site-packages/httpx/_models.py:829, in Response.iter_bytes(self, chunk_size)
827 chunker = ByteChunker(chunk_size=chunk_size)
828 with request_context(request=self._request):
--> 829 for raw_bytes in self.iter_raw():
830 decoded = decoder.decode(raw_bytes)
831 for chunk in chunker.decode(decoded):
File ~/checkouts/readthedocs.org/user_builds/scopesim/envs/stable/lib/python3.11/site-packages/httpx/_models.py:887, in Response.iter_raw(self, chunk_size)
884 chunker = ByteChunker(chunk_size=chunk_size)
886 with request_context(request=self._request):
--> 887 for raw_stream_bytes in self.stream:
888 self._num_bytes_downloaded += len(raw_stream_bytes)
889 for chunk in chunker.decode(raw_stream_bytes):
File ~/checkouts/readthedocs.org/user_builds/scopesim/envs/stable/lib/python3.11/site-packages/httpx/_client.py:124, in BoundSyncStream.__iter__(self)
123 def __iter__(self) -> typing.Iterator[bytes]:
--> 124 for chunk in self._stream:
125 yield chunk
File ~/checkouts/readthedocs.org/user_builds/scopesim/envs/stable/lib/python3.11/site-packages/httpx/_transports/default.py:111, in ResponseStream.__iter__(self)
109 def __iter__(self) -> typing.Iterator[bytes]:
110 with map_httpcore_exceptions():
--> 111 for part in self._httpcore_stream:
112 yield part
File ~/checkouts/readthedocs.org/user_builds/scopesim/envs/stable/lib/python3.11/site-packages/httpcore/_sync/connection_pool.py:367, in PoolByteStream.__iter__(self)
365 except BaseException as exc:
366 self.close()
--> 367 raise exc from None
File ~/checkouts/readthedocs.org/user_builds/scopesim/envs/stable/lib/python3.11/site-packages/httpcore/_sync/connection_pool.py:363, in PoolByteStream.__iter__(self)
361 def __iter__(self) -> Iterator[bytes]:
362 try:
--> 363 for part in self._stream:
364 yield part
365 except BaseException as exc:
File ~/checkouts/readthedocs.org/user_builds/scopesim/envs/stable/lib/python3.11/site-packages/httpcore/_sync/http11.py:349, in HTTP11ConnectionByteStream.__iter__(self)
347 with ShieldCancellation():
348 self.close()
--> 349 raise exc
File ~/checkouts/readthedocs.org/user_builds/scopesim/envs/stable/lib/python3.11/site-packages/httpcore/_sync/http11.py:341, in HTTP11ConnectionByteStream.__iter__(self)
339 try:
340 with Trace("receive_response_body", logger, self._request, kwargs):
--> 341 for chunk in self._connection._receive_response_body(**kwargs):
342 yield chunk
343 except BaseException as exc:
344 # If we get an exception while streaming the response,
345 # we want to close the response (and possibly the connection)
346 # before raising that exception.
File ~/checkouts/readthedocs.org/user_builds/scopesim/envs/stable/lib/python3.11/site-packages/httpcore/_sync/http11.py:210, in HTTP11Connection._receive_response_body(self, request)
207 timeout = timeouts.get("read", None)
209 while True:
--> 210 event = self._receive_event(timeout=timeout)
211 if isinstance(event, h11.Data):
212 yield bytes(event.data)
File ~/checkouts/readthedocs.org/user_builds/scopesim/envs/stable/lib/python3.11/site-packages/httpcore/_sync/http11.py:224, in HTTP11Connection._receive_event(self, timeout)
221 event = self._h11_state.next_event()
223 if event is h11.NEED_DATA:
--> 224 data = self._network_stream.read(
225 self.READ_NUM_BYTES, timeout=timeout
226 )
228 # If we feed this case through h11 we'll raise an exception like:
229 #
230 # httpcore.RemoteProtocolError: can't handle event type
(...)
234 # perspective. Instead we handle this case distinctly and treat
235 # it as a ConnectError.
236 if data == b"" and self._h11_state.their_state == h11.SEND_RESPONSE:
File ~/checkouts/readthedocs.org/user_builds/scopesim/envs/stable/lib/python3.11/site-packages/httpcore/_backends/sync.py:126, in SyncStream.read(self, max_bytes, timeout)
124 with map_exceptions(exc_map):
125 self._sock.settimeout(timeout)
--> 126 return self._sock.recv(max_bytes)
File ~/.asdf/installs/python/3.11.9/lib/python3.11/ssl.py:1295, in SSLSocket.recv(self, buflen, flags)
1291 if flags != 0:
1292 raise ValueError(
1293 "non-zero flags not allowed in calls to recv() on %s" %
1294 self.__class__)
-> 1295 return self.read(buflen)
1296 else:
1297 return super().recv(buflen, flags)
File ~/.asdf/installs/python/3.11.9/lib/python3.11/ssl.py:1168, in SSLSocket.read(self, len, buffer)
1166 return self._sslobj.read(len, buffer)
1167 else:
-> 1168 return self._sslobj.read(len)
1169 except SSLError as x:
1170 if x.args[0] == SSL_ERROR_EOF and self.suppress_ragged_eofs:
KeyboardInterrupt:
cmd = sim.UserCommands(use_instrument="MICADO", set_modes=["SCAO", "IMG_1.5mas"])
micado = sim.OpticalTrain(cmd)
micado.effects
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])
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
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")
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")
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")
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")
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())
