H.E.S.S. I (CT3)

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

from iactrace import Camera, Telescope, show_image, show_sensor_chain, show_telescope

Load telescope and camera

telescope = Telescope.from_yaml(
    "../../configs/HESS/CT3.yaml",
    n_samples=128,
    key=jax.random.key(42),
)
camera = Camera.from_yaml("../../configs/HESS/HESS1U.yaml")

E0804 10:37:17.800954 3111357 numa_hwloc.cc:121] Call to hwloc_set_cpubind() failed: Invalid argument [22]
E0804 10:37:17.801153 3111358 numa_hwloc.cc:121] Call to hwloc_set_cpubind() failed: Invalid argument [22]
E0804 10:37:17.801196 3111359 numa_hwloc.cc:121] Call to hwloc_set_cpubind() failed: Invalid argument [22]

telescope = telescope.apply_roughness(0, 12)

3D scene

We can visualize the telescope and the detector chain via trimesh scenes:

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

scene = show_sensor_chain(camera)
scene.show(viewer="jupyter")

Star field image

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

fov_deg = 6
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(ncols=1, figsize=(8, 8))
ax1 = show_image(image, camera.sensor_groups[0], ax=ax)

CPU times: user 6.08 s, sys: 569 ms, total: 6.65 s
Wall time: 6.5 s

Bokeh from a finite-distance point source

N_points = 1
key1, key2 = jax.random.split(jax.random.key(45))

x = jax.random.uniform(key1, N_points, minval=-1, maxval=1)
y = jax.random.uniform(key2, N_points, minval=-1, maxval=1)
z = jnp.ones(N_points) * 250

points = jnp.array([x, y, z]).T
f_points = jnp.ones(len(points))

rays = telescope.render(points, f_points, source_type="point")
image = camera.image(rays)

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

Spot diagram from explicit ray geometry

from iactrace.analysis import FlatFocalPlane

n_rays = 1000
key1, key2 = jax.random.split(jax.random.key(123))

r = 6.0 * jnp.sqrt(jax.random.uniform(key1, (n_rays,)))
theta = jax.random.uniform(key2, (n_rays,)) * 2 * jnp.pi
ray_origins = jnp.stack(
    [r * jnp.cos(theta), r * jnp.sin(theta), jnp.ones(n_rays) * 20.0],
    axis=1,
)

tilt_angle = 0 * jnp.pi / 180
ray_directions = jnp.broadcast_to(
    jnp.array([jnp.sin(tilt_angle), 0.0, -jnp.cos(tilt_angle)]),
    (n_rays, 3),
)
ray_values = jnp.ones(n_rays)

rays, _ = telescope.trace(ray_origins, ray_directions, ray_values)
image = camera.image(rays)

focal_plane = FlatFocalPlane()
hits = focal_plane.intersect(rays)

fig, ax = plt.subplots(ncols=2, figsize=(14, 6))
show_image(image, camera.sensor_groups[0], ax=ax[0])
h2 = ax[1].hist2d(
    hits.xy_local[hits.alive, 0],
    hits.xy_local[hits.alive, 1],
    bins=30,
    norm="log",
)