Parabolic Telescope

import jax
import jax.numpy as jnp
import matplotlib.pyplot as plt

from iactrace import Camera, Telescope, show_image, show_telescope

Load telescope and camera

telescope = Telescope.from_yaml(
    "../../configs/BASIC/Parabolic.yaml",
    n_samples=8192,
    key=jax.random.key(42),
)
camera = Camera.from_yaml("../../configs/BASIC/default_camera.yaml")

E0804 11:07:52.401499 3121878 numa_hwloc.cc:121] Call to hwloc_set_cpubind() failed: Invalid argument [22]
E0804 11:07:52.401643 3121879 numa_hwloc.cc:121] Call to hwloc_set_cpubind() failed: Invalid argument [22]
E0804 11:07:52.401701 3121880 numa_hwloc.cc:121] Call to hwloc_set_cpubind() failed: Invalid argument [22]

Apply 1 arcsecond of surface roughness:

telescope = telescope.apply_roughness(0, 1)

3D scene

scene = show_telescope(telescope)
scene.show(viewer="jupyter")

Star field image

n_stars = 100
key = jax.random.key(42)
key1, key2 = jax.random.split(key)

fov_deg = 5
fov_rad = fov_deg * jnp.pi / 180

x = jax.random.uniform(key1, (n_stars,), minval=-fov_rad / 2, maxval=fov_rad / 2)
y = jax.random.uniform(key2, (n_stars,), minval=-fov_rad / 2, maxval=fov_rad / 2)
z = -jnp.ones(n_stars)

stars = jnp.stack([x, y, z], axis=1)
stars = stars / jnp.linalg.norm(stars, axis=1, keepdims=True)

f_stars = jax.random.uniform(jax.random.key(4242), shape=(len(stars),))

rays = telescope.render(stars, f_stars, source_type="parallel")
image = camera.image(rays)

fig, ax = plt.subplots()
ax1 = show_image(image, camera.sensor_groups[0], ax=ax)

Spot diagram for on- and off-axis stars

from iactrace.analysis import FlatFocalPlane

stars = jnp.array([[0, 0, -1], [0.004, 0.000, -1]])
stars = stars / jnp.linalg.norm(stars, axis=1, keepdims=True)
f_stars = jnp.array([1.0, 1.0])

rays = telescope.render(stars, f_stars, source_type="parallel")

focal_plane = FlatFocalPlane()
hits = focal_plane.intersect(rays)
points = hits.xy_local
weight = hits.values * hits.hit_mask

fig, ax = plt.subplots(ncols=2, figsize=(8, 4))
chunks = jnp.array_split(points, 2)
wchunks = jnp.array_split(weight, 2)
for i, chunk in enumerate(chunks):
    ax[i].hist2d(chunk[:, 0], chunk[:, 1], weights=wchunks[i], bins=50)
    ax[i].axis("off")