# Etendue¶

This example code calculates the etendue of a pinhole and compares it with a simpler analytical calculation. The distance between the pinhole and the observer is varied, showing that the ray-tracer exhibits the correct limits. The method of calculation is described in Carr, M., Meakins, A., et al. “Description of complex viewing geometries of fusion tomography diagnostics by ray-tracing.” Review of Scientific Instruments 89.8 (2018): 083506.


import numpy as np
import matplotlib.pyplot as plt

from raysect.core import Point3D, Vector3D, rotate_basis, translate, Ray as CoreRay
from raysect.core.math.sampler import DiskSampler3D, RectangleSampler3D, TargettedHemisphereSampler
from raysect.optical import World
from raysect.primitive import Box, Cylinder, Subtract
from raysect.optical.material import AbsorbingSurface, NullMaterial

R_2_PI = 1 / (2 * np.pi)

world = World()

# Setup pinhole
target_plane = Box(Point3D(-10, -10, -0.000001), Point3D(10, 10, 0.000001))
hole = Cylinder(0.001, 0.001, transform=translate(0, 0, -0.0005))
pinhole = Subtract(target_plane, hole, parent=world, material=AbsorbingSurface())

target = Cylinder(0.0012, 0.001, transform=translate(0, 0, -0.0011), parent=world, material=NullMaterial())

def analytic_etendue(area_det, area_slit, distance, alpha, gamma):

return area_det * area_slit * np.cos(alpha/360 * (2*np.pi)) * np.cos(gamma/360 * (2*np.pi)) / distance**2

# generate the transform to the detector position and orientation
detector_transform = translate(0, 0, distance) * rotate_basis(Vector3D(0, 0, -1), Vector3D(0, -1, 0))

# generate bounding sphere and convert to local coordinate system
sphere = target.bounding_sphere()

# instance targetted pixel sampler
targetted_sampler = TargettedHemisphereSampler(spheres)

solid_angle = 2 * np.pi
etendue_sampled = solid_angle * detector_area

etendues = []
for i in range(batches):

# sample pixel origins
origins = point_sampler(samples=ray_count)

passed = 0.0
for origin in origins:

# obtain targetted vector sample
direction, pdf = targetted_sampler(origin, pdf=True)
path_weight = R_2_PI * direction.z/pdf

origin = origin.transform(detector_transform)
direction = direction.transform(detector_transform)

while True:

# Find the next intersection point of the ray with the world
intersection = world.hit(CoreRay(origin, direction))

if intersection is None:
passed += 1 * path_weight
break

elif isinstance(intersection.primitive.material, NullMaterial):
hit_point = intersection.hit_point.transform(intersection.primitive_to_world)
origin = hit_point + direction * 1E-9
continue

else:
break

if passed == 0:
raise ValueError("Something is wrong with the scene-graph, calculated etendue should not zero.")

etendue_fraction = passed / ray_count
etendues.append(etendue_sampled * etendue_fraction)

etendue = np.mean(etendues)
etendue_error = np.std(etendues)

return etendue, etendue_error

area = 0.001**2 * np.pi

detector_etendue = 0.001**2 * np.pi * np.pi  # etendue = A * omega * 1/2

distance_samples = [10**i for i in np.arange(-4, -1.1, 0.10)]

analytic_values = []
raytraced_values = []
raytraced_errors = []

for d in distance_samples:
analytic_values.append(analytic_etendue(area, area, d, 0, 0))
value, error = raytraced_etendue(d)
raytraced_values.append(value)
raytraced_errors.append(error)

analytic_values = np.array(analytic_values)
raytraced_values = np.array(raytraced_values)
raytraced_errors = np.array(raytraced_errors)

plt.ion()

plt.figure()
ax = plt.gca()
ax.set_xscale("log", nonposx='clip')
ax.set_yscale("log", nonposy='clip')
plt.axhline(y=detector_etendue, linestyle='--', color='k', label='detector etendue')
plt.plot(distance_samples, analytic_values, label='analytic etendue')
plt.errorbar(distance_samples, raytraced_values, raytraced_errors, label='ray-traced etendue')
plt.xlabel('Distance between slit and detector (m)')
plt.ylabel('Etendue (m^2 str)')
plt.title('Ray-traced VS approximate analytic etendue')
plt.legend()

# plt.figure()
# ax = plt.gca()
# ax.set_xscale("log", nonposx='clip')
# plt.errorbar(distance_samples, np.abs(raytraced_values-analytic_values)/raytraced_values, raytraced_errors/raytraced_values)
# plt.xlim(0.001, 0.1)
# plt.ylim(0, 0.5)
# plt.xlabel('Distance between slit and detector (m)')
# plt.ylabel('Fractional error')
# plt.show()