Calculate RTM using multiple spectral lines

In this example the radiance from multiple spectral lines is calculated simultaniously, which significantly reduces the computational time compared to the case when the radiance for each spectral line is calculated one by one.

Note

The adaptive sampling can be used effectively only if the emission profiles are not very different for different spectral lines. Therefore it is not recomended to mix the emissions from impurities and hydrogen isotopes in a single Raysect simulation.

import numpy as np
from raysect.optical import World, translate, rotate
from raysect.optical.material import AbsorbingSurface
from raysect.optical.observer import PinholeCamera, SpectralRadiancePipeline2D
from raysect.primitive import Cylinder, Subtract
from cherab.iter.machine import load_pfc_mesh
from raysect.core.workflow import MulticoreEngine
from cherab.iter.divimp import read_divimp_reg
from cherab.iter.observer import FovCamera, FoV, SpectralAdaptiveSampler2D
from cherab.iter.emitter import RegularRPhiZIntegrator, DiscreteSpectralRPhiZEmitter


# === Configuration ===========================================================
case = "2226"
view = "DIM_EP1"
dphi = 0.1
dtheta = 0.1
step = 0.0025
cutoff = 0.1
cutoff_frac = 0.01
min_samples = 1000
nr = 460
nz = 945
rmin, rmax = 3.95, 8.55
zmin, zmax = -4.65, 4.8
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()

# === wall creation ===========================================================
central_column = Cylinder(rmin, zmax - zmin, material=AbsorbingSurface(), parent=world, transform=translate(0, 0, zmin))
load_pfc_mesh(world)

# === plasma creation =========================================================
askukush_filename = "../cherab/iter/divimp/data/{}/lines_{}.txt".format(case, case)
with open(askukush_filename, "r") as f:
    keys_all = f.readlines()[2].split()[4:]
profcols = {key: i + 4 for i, key in enumerate(keys_all)}
divimpprof = read_divimp_reg(askukush_filename, r, z, profcols=profcols, skiprows=3)
keys = [key for key in keys_all if divimpprof[key].max() > 0]  # removing all zero-only profiles
# Note that all deuterium (or all non-deuterium) spectral lines also can be removed if needed
# keys = [key for key in keys if key[0] != 'D']
# keys = [key for key in keys if key[0] == 'D']
lines_wavelength = np.array([0.1 * float(key[2:]) for key in keys])  # spectral lines wavelengths in nm
min_wavelength = lines_wavelength.min()
max_wavelength = lines_wavelength.max() + 0.1
delta_wavelength = np.diff(np.sort(lines_wavelength)).min()
bins = int((max_wavelength - min_wavelength) / delta_wavelength) + 1
emissivity = np.zeros((nr, nz, len(keys)))
for i, key in enumerate(keys):
    # SpectralRadiancePipeline2D does not integrate the radiance over the wavelength, so the conversion to spectral radiance is not required
    emissivity[:, :, i] = divimpprof[key]
emitt_material = DiscreteSpectralRPhiZEmitter(
    emissivity,
    lines_wavelength,
    min_wavelength,
    max_wavelength,
    bins,
    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 = (int(round(phi / dphi)), int(round(theta / dtheta)))
pipeline = SpectralRadiancePipeline2D()
sampler = SpectralAdaptiveSampler2D(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 = min_wavelength
camera.max_wavelength = max_wavelength
camera.spectral_bins = bins
camera.spectral_rays = 1
camera.pixel_samples = min_samples
camera.ray_extinction_prob = 0
camera.render_engine = MulticoreEngine(8)

# === raytracing ==============================================================
imsize = camera.pixels[0] * camera.pixels[1]
npix = imsize
while npix > cutoff_frac * imsize:
    camera.observe()
    errors = pipeline.frame.errors()[:, :, emitt_material.indx]
    mean = pipeline.frame.mean[:, :, emitt_material.indx]
    ind = np.where((mean > 1.0e-4 * mean.max((0, 1))[None, None, :]).sum(2))
    rel_err = (errors[ind] / mean[ind]).max(1)
    npix = np.sum(rel_err > cutoff)
    print("Max. error: %.2e. Pixels to add samples for: %d" % (rel_err.max(), npix))
    for i, key in enumerate(keys):
        fname = "%s_%s_case_%s" % (key, view, case)
        np.savetxt(fname + ".txt", mean[:, :, i])
        np.savetxt(fname + "_errors.txt", errors[:, :, i])