Usage of RTM

import numpy as np
from raysect.optical import World, translate, rotate
from raysect.optical.material import AbsorbingSurface
from raysect.optical.observer import PinholeCamera, RadiancePipeline2D, MonoAdaptiveSampler2D
from raysect.primitive import Cylinder, Subtract
from cherab.iter.machine import import_iter_mesh
from cherab.iter.divimp import read_divimp_mesh
from cherab.iter.observer import FovCamera, FoV
from cherab.iter.emitter import RegularRPhiZIntegrator, RegularRPhiZEmitter
from matplotlib import pyplot as plt
from cherab.iter.tools import RTM, plot_fov_2D


# === Configuration ===========================================================
case = "2226"
view = "DIM_EP1"
dphi = 0.1
dtheta = 0.1
step = 0.0025
cutoff = 0.1
min_samples = 1000
nr = 450
nz = 950
rmin, rmax = 3.95, 8.55
zmin, zmax = -4.65, 4.8
wlmin, wlmax = 655.6, 656.6
eps_rz = 1.0e-4
r, dr = np.linspace(rmin, rmax, nr + 1, retstep=True)
z, dz = np.linspace(zmin, zmax, nz + 1, retstep=True)


world = World()

# === Load ITER PFC ===========================================================
wall = import_iter_mesh(world)

# === Plasma ==================================================================
askukush_filename = f"../cherab/iter/divimp/lines/lines_{case}.txt"
divimpprof = read_divimp_mesh(askukush_filename, r, z, ["D_6561"])
emissivity = divimpprof["D_6561"] / (wlmax - wlmin)  # spectral emissivity

# === Camera ==================================================================
fov = FoV(view)
pupils = fov.view["PUPILS"]
rotate_angles = fov.euler()
phi, theta = fov.view["FOV"]
pixels = (int(round(phi / dphi)), int(round(theta / dtheta)))

pipeline = RadiancePipeline2D()
sampler = MonoAdaptiveSampler2D(pipeline, fraction=1.0, cutoff=cutoff, min_samples=min_samples, ratio=10000.0)
camera = FovCamera(pixels, fov=(phi, theta), parent=world, pipelines=[pipeline], frame_sampler=sampler)
camera.transform = translate(*pupils) * rotate(*rotate_angles)
camera.min_wavelength = wlmin
camera.max_wavelength = wlmax
camera.spectral_bins = 1
camera.spectral_rays = 1
camera.pixel_samples = min_samples
camera.ray_extinction_prob = 0

# === Create RTM ==============================================================
rtm = RTM()
# rtm.load("test.npy")
rtm.create(camera, r, z, wall=wall, pixel_samples=10)
# rtm.create(camera, r, z, wall=wall, pixel_samples=1, mask=emissivity)
rtm.save("test.npy")
radiation = rtm.construct(emissivity)


# === visualization ===========================================================
phi_max = phi / 2.0
theta_max = theta / 2.0
phi_horizontal = np.linspace(-phi_max, phi_max, pixels[0], endpoint=True)
theta_vertical = np.linspace(theta_max, -theta_max, pixels[1], endpoint=True)

plot_fov_2D(phi_horizontal, theta_vertical, radiation)