Calculate RTM using kukushkin profile

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


# === Configuration ===========================================================
case = "2226"
view = "DIM_EP1"
element = "D_4860"
step = 0.0025
cutoff = 0.1
cutoff_frac = 0.01
min_samples = 100
ratio = 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
fraction = 0.01
r, dr = np.linspace(rmin, rmax, nr + 1, retstep=True)
z, dz = np.linspace(zmin, zmax, nz + 1, retstep=True)

world = World()

# === Wall  ===================================================================
central_column = Cylinder(rmin, zmax - zmin, material=AbsorbingSurface(), parent=world, transform=translate(0, 0, zmin))
mesh = import_iter_mesh(world)

# === Plasma  =================================================================
divimpprof = read_divimp_mesh(2226, r, z, [element], unit="photon/m3")
emissivity = divimpprof[element]

emitt_material = RegularRPhiZEmitter(emissivity, rmin, dr, dz, integrator=RegularRPhiZIntegrator(step))
emitter = Subtract(
    Cylinder(rmax - eps_rz, zmax - zmin - eps_rz), Cylinder(rmin + eps_rz, zmax - zmin - eps_rz), parent=world
)
emitter.material = emitt_material
emitter.transform = translate(0, 0, zmin)

# === Camera ==================================================================
fov = FoV(view)
pupils = fov.view["PUPILS"]
rotate_angles = fov.euler()
phi, theta = fov.view["FOV"]
pixels = fov.view["PIXELS"]

pipeline = RadiancePipeline2D()
sampler = MonoAdaptiveSampler2D(pipeline, fraction=fraction, cutoff=cutoff, min_samples=min_samples, ratio=ratio)
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 = 10
camera.spectral_rays = 1
camera.pixel_samples = min_samples
camera.ray_extinction_prob = 0

# === Ray Tracing =============================================================
plt.ion()
camera.observe()
plt.ioff()
plt.show()

np.savetxt("{}_{}_case_{}.txt".format(element, view, case), pipeline.frame.mean)
np.savetxt("{}_{}_case_{}_errors.txt".format(element, view, case), pipeline.frame.errors())