plasma properties not affecting emission from cherab.core Beam()

[Cam42069]

Dear Cherab Team,

I am new to using cherab and have been trying to setup a CXRS emission sim for the HSX stellarator but I'm running into issues where varying plasma parameters has no effect on the observed emission. I've been through the docs, examples, and repo to no avail yet.

Specifically, changing the plasma's geometry is what is bothering me most. To my understanding you can define a raysect primitive like Box(), add it to the plasma object's geometry field and the profiles you defined in species will fit to your geometry along some axis. However, no matter the plasma size or shape the observed emission stays the same. It also remains the same when I change the plasma species' properties like dens, temp, etc.

The only way I can change the beam's emission is by changing the beam parameters, which makes me think the observation is purely of the beam with no regard to the plasma.

I have been basing my work off of the Active Spectroscopy CXS example and I have attached a sample code that shows this problem with a Box() and Cylinder() primitive.

Is this just a misunderstanding on my part or is this a bug in the code?

cxrs_beam_issue.py

Thank you,
Cam

[Mateasek]

Hi @Cam42069

Thanks for giving Cherab a try. Also thanks for the example. I hope I understand you problem correcty:

1. You expect to see some change in the image properties when you modify your plasma parameters, but you don't.
2. You expected to see radiation from plasma, but you don't.
3. You would expect to see difference in the active CX component when you change plasma properties

If yes, then here are the main reasons:

No plasma emission

From the example I see that you didn't add any radiation model to plasma. If you want to see plasma radiation, you need to specify what radiation plasma should be working with, in your example it is probably the 'passive lines' in CXRS spectra:

exc_cv = ExcitationLine(cVI_8_7, atomic_data=adas)
rec_cvi  = RecombinationLine(cVI_8_7, atomic_data=adas)
plasma.models = [exc_cv, rec_cvi]

Missing geometry transform in plasma

You correctly shifted the plasma object within your scene, but you also need to specify Plasma.geometry_transform. As you pointed out correctly, plasma "lives" in a geometric primitive (e.g. cylinder, box, etc.). This object you shift within you scene using plasma.transform. But there is a catch. The geometric object internally uses its own coordinate base, so you also need to properly shift it by correctly assigning Plasma.geometry_transform. I will give you examples from the tokamak universe since it is simple and you will have to figure out the right application for your stellarator. You would usually use Cylinder primitive for a tokamak plasma geometry. Since most tokamaks assign Z=0 at midplane, you would shift your plasma downwards to cover the bottom part of your plasma, e.g.:

z_bottom = -2.5
plasma_translate = translate(0, 0, z_bottom)

This shifts the primitive downwards within your scene. However, this causes a point with Z=0 in the "world" basis to have Z=2.5 in the "plasma" basis. If you specify your plasma profiles with midplane at Z=0, you will see a vertical shift caused by the translation. To counteract this, you need to apply the plasma.geometry_transform which internally applies affine transform within the plasma primitive to fix this problem. For stellarators there might be a bit more than a vertical shift and I will have to leave it up to you to come up with the right expression.

Debugging Plots

Moving to the third problem. I think you should look at the spectra to see difference in the active emission when changing plasma properties. The 2D colormap of total intensity is probably hiding the changes. Simulate more spectral bins with your camera and plot a spectrum for some pixel and you will see the spectrum reacts as you would assume.

Modified Example

I took the liberty of rewriting your example and applied the fixes I suggested. You can see changes in plasma radiation and spectra. Please have in mind the plasma.geometry_transform I applied is probably not correct, you need to figure out how to do it correctly.

### imports ###
from itertools import accumulate
import numpy as np
import matplotlib.pyplot as plt
from scipy.constants import electron_mass, atomic_mass

from raysect.core import translate, rotate_basis, Point3D, Vector3D
from raysect.primitive import Box, Cylinder
from raysect.optical import World, Ray
from raysect.optical.observer import PinholeCamera, SpectralRadiancePipeline2D

from cherab.core import Plasma, Beam, Maxwellian, Species
from cherab.core.math import sample3d, ScalarToVectorFunction3D
from cherab.core.atomic import hydrogen, deuterium, carbon, Line
from cherab.core.model import (
    SingleRayAttenuator,
    BeamCXLine,
    ExcitationLine,
    RecombinationLine,
)
from cherab.tools.plasmas.slab import NeutralFunction, IonFunction
from cherab.openadas import OpenADAS


def plot_pipeline(
    pipeline: SpectralRadiancePipeline2D,
    central_pixel_x: int = 40,
    central_pixel_y: int = 55,
) -> None:
    _, ax = plt.subplots(figsize=(10, 10))  # define the figure
    ax.imshow(pipeline.frame.mean[:, :].sum(axis=-1))
    ax.axvline(central_pixel_x, 0, 100, color="red")
    ax.axhline(central_pixel_y, 0, 100, color="red")
    ax.set_title("Total Radiance Image")
    ax.figure.tight_layout()
    ax.figure.show()

    _, ax = plt.subplots(figsize=(10, 10))  # define the figure
    ax.plot(
        pipeline.wavelengths, SpecRad2D.frame.mean[central_pixel_y, central_pixel_x, :]
    )
    ax.set_title("Central Pixel Spectral Radiance")
    ax.figure.tight_layout()
    ax.figure.show()


###############  Main Cherab Inputs  ###############

cylinder = True

# Plasma volume geometry
ppos = np.array([0.1301, 1.4442, 0.2365])
plasma_translate = translate(*ppos)
plasma_rotate_basis = rotate_basis(
    Vector3D(-0.48249505265067577, 0.0, -0.8758986951512268), Vector3D(0, 0, 1)
)
plasma_transform = plasma_translate * plasma_rotate_basis
radius = 1
if cylinder:
    plasma_geometry = Cylinder(radius, height=5, transform=plasma_transform)
else:
    width = 1.0
    length = 1.0
    height = 3.0
    plasma_geometry = Box(
        Point3D(0, -width / 2, -height / 2), Point3D(length, width / 2, height / 2)
    )


# Pinhole camera translation
cpos = np.array([-39.367, 106.097, 38.354]) / 100
cam_translate = translate(*cpos)
cam_rotate_basis = rotate_basis(
    Vector3D(0.74046228, 0.58483411, -0.33118675), Vector3D(0, 0, 1)
)
cam_transform = cam_translate * cam_rotate_basis


# NB properties
beam_energy = 30_000  # eV
beam_power = 50_000  # J?
beam_sigma = 0.05
beam_divergence_x = 0.5
beam_divergence_y = 0.5
beam_length = 4

### Species params ###
# Neutral Species
peak_h0_density = 4.4e13
h0_NeutralFunction_sigma = 0.1
h0_temperature = 0.5  # [eV]
#  - vel profile -Assume-> follows T_h1

# Ion Species:
peak_h1_density = 5e19  # [m^-3]
peak_h1_temperature = 2500  # [eV]
peak_h1_vel = 120e3  # [m/s]?

# Electron Species:
#  - n Profile -Assume-> follows n_h1
#  - Ti Profile -Assume-> follows T_h1
#  - vel profile -Assume-> follows T_h1

# Impurity Species:
#  - n Profile -Assume-> scaled by n_h1
#  - Ti Profile -Assume-> follows T_h1
#  - vel profile -Assume-> follows T_h1
impurities = [
    (carbon, 6, 0.005)
]  # impurity, ionization, concentration (to main ion species)


## Globally Assigned ##
peak_density = peak_h1_density  # peak density of bulk plasma
peak_temperature = peak_h1_temperature  # peak temperature of bulk plasma
peak_velocity = peak_h1_vel
#
pedestal_top = 1  # Either the distance the plasma is from the walls or a scalar value that scales the distance the plasma is from the walls/bounding box
impurities = [
    (carbon, 6, 0.005),
    (carbon, 5, 0.001),
]  # impurity, ionization, concentration (to main ion species)


### Setup Scene ###
world = World()
adas = OpenADAS(permit_extrapolation=True, missing_rates_return_null=True)
#
plasma = Plasma(parent=world)
plasma.atomic_data = adas
plasma.geometry = plasma_geometry  ### plasma.geometry = Box(Point3D(0, -width / 2, -height / 2), Point3D(length, width / 2, height / 2))
plasma.geometry_transform = plasma_translate  # The plasma transform has to be applied internally within the plasma volume
#
species = []


### Generate Plasma Profiles per species, add to plasma
# Plasma Velocity Profile
vy_profile = IonFunction(peak_velocity, 0, pedestal_top=pedestal_top)
velocity_profile = ScalarToVectorFunction3D(0, vy_profile, 0)

# define neutral species distribution
h0_density = NeutralFunction(
    peak_h0_density, sigma=h0_NeutralFunction_sigma, pedestal_top=pedestal_top
)
h0_distribution = Maxwellian(
    h0_density, h0_temperature, velocity_profile, hydrogen.atomic_weight * atomic_mass
)
species.append(Species(hydrogen, 0, h0_distribution))

# define hydrogen ion species distribution
h1_density = IonFunction(peak_h1_density, 0, pedestal_top=pedestal_top)
h1_temperature = IonFunction(peak_h1_temperature, 0, pedestal_top=pedestal_top)
h1_distribution = Maxwellian(
    h1_density, h1_temperature, velocity_profile, hydrogen.atomic_weight * atomic_mass
)
species.append(Species(hydrogen, 1, h1_distribution))

# add impurities
if impurities:
    for impurity, ionisation, concentration in impurities:
        imp_density = IonFunction(
            peak_density * concentration, 0, pedestal_top=pedestal_top
        )
        imp_temperature = IonFunction(peak_temperature, 0, pedestal_top=pedestal_top)
        imp_distribution = Maxwellian(
            imp_density,
            imp_temperature,
            velocity_profile,
            impurity.atomic_weight * atomic_mass,
        )
        species.append(Species(impurity, ionisation, imp_distribution))

# define the electron distribution
e_density = IonFunction(peak_density, 0, pedestal_top=pedestal_top)
e_temperature = IonFunction(peak_temperature, 0, pedestal_top=pedestal_top)
e_distribution = Maxwellian(e_density, e_temperature, velocity_profile, electron_mass)

# add plasma species and other properties
plasma.b_field = Vector3D(0, 0, 0)
plasma.electron_distribution = e_distribution
plasma.composition = species


## Define Beam() & add params ###
cVI_8_7 = Line(carbon, 5, (8, 7))

### Add passive plasma radiation
exc_cv = ExcitationLine(cVI_8_7)
rec_cvi = RecombinationLine(cVI_8_7)
plasma.models = [exc_cv, rec_cvi]


### Observer ###
SpecRad2D = SpectralRadiancePipeline2D(
    accumulate=False
)  # don't accumulate observations
camera = PinholeCamera(
    (100, 100), parent=world, transform=cam_transform, pipelines=[SpecRad2D]
)
camera.spectral_rays = 1
camera.spectral_bins = 1
camera.pixel_samples = 10
camera.spectral_bins = 100
camera.max_wavelength = 533
camera.min_wavelength = 527

### Observe plasma emission only ###
camera.observe()
plt.ion()
plot_pipeline(camera.pipelines[0])
plt.ioff()


## Add beam to the scene ##
integration_step = 0.0025
beam_transform = plasma_transform
beam_components = []

world2 = World()
# Initialize beam for all 3 energy components
for i in range(3):
    beam_i = Beam(parent=world, transform=beam_transform)
    beam_i.plasma = plasma
    beam_i.atomic_data = adas
    beam_i.energy = beam_energy / (i + 1)
    beam_i.power = beam_power
    beam_i.element = deuterium
    beam_i.sigma = beam_sigma
    beam_i.divergence_x = beam_divergence_x
    beam_i.divergence_y = beam_divergence_y
    beam_i.length = beam_length
    beam_i.attenuator = SingleRayAttenuator(clamp_to_zero=True)
    beam_i.models = [BeamCXLine(cVI_8_7)]
    beam_i.integrator.step = integration_step
    beam_i.integrator.min_samples = 10
    #
    beam_components.append(beam_i)


### Observe plasma and beam emission ###
camera.observe()
plt.ion()
plot_pipeline(camera.pipelines[0])
plt.ioff()