Fix SOLPSTotalRadiatedPower and make solps_total_radiated_power() more universal

[vsnever]

When calculating the spectral emissivity in SOLPSTotalRadiatedPower.emission_function(), the total emissivity is divided by 4π, wavelength range and the number of spectral bins. It seems to me that dividing by the number of spectral bins is unnecessary here and breaks the normalisation.

I propose a change from:

spectrum.samples_mv[:] = self.total_rad.evaluate(point.x, point.y, offset_z) / (4 * pi * wvl_range * spectrum.bins)

to:

spectrum.samples_mv[:] = self.total_rad.evaluate(point.x, point.y, offset_z) / (4 * pi * wvl_range)

Note, however, that even after this correction, this emitter will not work properly if dispersive rendering is enabled.

I also propose to rename SOLPSTotalRadiatedPower to SOLPSRadiatedPower, solps_total_radiated_power() to solps_radiated_power() and convert the latter to a more universal function. The new function will look like this:

def solps_radiated_power(world, solps_mesh, solps_radiated_power, step=0.01):

    if not isinstance(solps_mesh, SOLPSMesh):
        raise TypeError('Argument solps_mesh must be a SOLPSMesh instance.')
    if not isinstance(solps_radiated_power, SOLPSFunction2D):
        raise TypeError('Argument solps_radiated_power must be a SOLPSFunction2D instance.')

    radiated_power_f3d = AxisymmetricMapper(solps_radiated_power)
    outer_radius = solps_mesh.mesh_extent['maxr']
    inner_radius = solps_mesh.mesh_extent['minr'] - 0.001
    plasma_height = solps_mesh.mesh_extent['maxz'] - solps_mesh.mesh_extent['minz']
    lower_z = solps_mesh.mesh_extent['minz']

    main_plasma_cylinder = Cylinder(outer_radius, plasma_height, parent=world,
                                    material=SOLPSRadiatedPower(radiated_power_f3d, vertical_offset=lower_z, step=step),
                                    transform=translate(0, 0, lower_z))

    return main_plasma_cylinder

After this change, the solps_radiated_power() can be used for example to calculate the wavelength-integrated H-alpha signals (issue #48 and PR #49).

@Mateasek, @jacklovell, what do you think?

[Mateasek]

Hi @vsnever,

good catch with the normalisation!

I also like the proposed changes to the class and methods. Adds more variability, so I think you can go ahead, if @jacklovell is not against. We can always keep the solps_total_radiated_power as a utility function which would simply call the proposed solps_radiated_power.

[jacklovell]

I think the idea behind dividing by the number of spectral bins is to distribute the power evenly among each bin. The division by spectrum.nbins is necessary to conserve power: if you didn't then summing the spectrum over all bins (i.e. integrating over wavelengths) would give you nbins times the total power. To me it looks like the existing function is correct, though I would be happy to see a demonstration which proves otherwise.

I'm not against providing another utility function to create an emitter from the SOLPS data. I would prefer some more instructive names though: make_solps_emitter for the function, and radiation_function instead of the solps_radiated_power argument. Otherwise you have a function and an argument to the function with the same name, which will be confusing. I always through the solps_total_radiated_power function was oddly named anyway: it implies computation of a property of the simulation, not the creation of an emitting primitive.

To be honest, I'm struggling to see what benefit SOLPSTotalRadiatedPower gives over the standard Cherab core idiom of using a RadiationFunction material with a VolumeTransform modifier, and just passing AxisymmetricMapper(whatever-solps-2d-field-you-want) to RadiationFunction. See https://cherab.github.io/documentation/demonstrations/radiation_loads/radiation_function.html#radiation-function. Was there a reason we didn't use this? Maybe historical?

[vsnever]

Hi @jacklovell, Thanks for the tip with VolumeTransform modifier. I didn't know about it. With this modifier, the RadiationFunction is enough indeed.

I think the idea behind dividing by the number of spectral bins is to distribute the power evenly among each bin. The division by spectrum.nbins is necessary to conserve power: if you didn't then summing the spectrum over all bins (i.e. integrating over wavelengths) would give you nbins times the total power. To me it looks like the existing function is correct, though I would be happy to see a demonstration which proves otherwise.

The power is divided by the wavelength range of the ray, so the integral over the wavelength gives the power again, not the nbins times the power.

[vsnever]

Looks like we don't need SOLPSTotalRadiatedPower anymore. Then the make_solps_emitter() function will be like this:

def make_solps_emitter(solps_mesh, radiation_function, parent=None, step=0.01):

    if not isinstance(solps_mesh, SOLPSMesh):
        raise TypeError('Argument solps_mesh must be a SOLPSMesh instance.')
    if not isinstance(radiation_function, SOLPSFunction2D):
        raise TypeError('Argument radiation_function must be a SOLPSFunction2D instance.')

    radiation_function_3d = AxisymmetricMapper(radiation_function)
    outer_radius = solps_mesh.mesh_extent['maxr']
    inner_radius = solps_mesh.mesh_extent['minr']
    plasma_height = solps_mesh.mesh_extent['maxz'] - solps_mesh.mesh_extent['minz']
    lower_z = solps_mesh.mesh_extent['minz']
    emitter = RadiationFunction(radiation_function_3d, step=step)
    material = VolumeTransform(emitter, transform=translate(0, 0, -lower_z))

    plasma_volume = Subtract(Cylinder(outer_radius, height), Cylinder(inner_radius, height),
                             material=material, parent=parent, transform=translate(0, 0, lower_z))

    return plasma_volume

[jacklovell]

The power is divided by the wavelength range of the ray, so the integral over the wavelength gives the power again, not the nbins times the power.

Ah yes, I see now.

[vsnever]

Fixed in #52.