prad_peaking returns NaNs prior to shot termination

[yumouwei]

Problem

The prad_peaking parameter method currently implemented in disruption-py returns NaNs prior to the termination of a shot. This issue occurs in both disruptive and non-disruptive shots, although it's more pronounced in disruptive ones.

1150805012: flat-top disruption

1150805013: non-disruptive shot

Analysis of the problem

1. The prad_peaking method first fetches signals from all available chords in the AXA and AXJ bolometers. At each time point, the function calculate the "distance of the chord" (assuming it's the chord's sight line) with respect to the plasma's magnetic axis. If the chord's is within 20% of the plasma's minor radius then the chord's signal is considered as part of the core radiation. The peaking factor is then calculated using mean(core radiation)/mean(all radiation)

disruption-py/disruption_py/shots/parameter_methods/cmod/basic_parameter_methods.py

Lines 1251 to 1271 in 87df64e

	for i in range(len(params.shot_props.times)):
	core_radiation = np.array([])
	all_radiation = np.array([])
	if got_axa:
	axa_dist = np.sqrt((r_axa - r0[i]) ** 2 + (z0[i] - z_axa) ** 2)
	axa_core_index = axa_dist < 0.2 * a_minor[i]
	core_radiation = np.append(
	core_radiation, axa_interp[axa_core_index, i]
	)
	all_radiation = np.append(all_radiation, axa_interp[:, i])
	if got_axj:
	axj_dist = np.sqrt((r_axj - r0[i]) ** 2 + (z0[i] - z_axj) ** 2)
	axj_core_index = axj_dist < 0.2 * a_minor[i]
	core_radiation = np.append(
	core_radiation, axj_interp[axj_core_index, i]
	)
	all_radiation = np.append(all_radiation, axj_interp[:, i])
	try:
	prad_peaking[i] = np.nanmean(core_radiation) / np.nanmean(all_radiation)
	except:
	prad_peaking[i] = np.nan

2. In those NaN time slices, the core radiation array does not contain any value. This indicates the method could not find any chord that satisfies the above criteria. This could potentially be caused by the plasma's motion and shrinkage prior to a disruption.

Potential solutions

0. Determine if this issue is worth investigating -- if data at these time slices are going to be useful for disruption studies.
1. Understand the geometry of the bolometer diagnostic(s)
2. Understand the design choices in the method:
   • Why chose 20% of a_minor as the cutoff?
   • Could we possibly fit the bolometer profile into a smooth function (e.g. a gaussian) and then determine the peaking factor?
   • We should also check the bolometer's radiation profile and compare it with the plasma's magnetic axis & minor radius

Related issues & PRs

• a_minor tree source in calculating Prad_peaking #92
• interp1 Fails in _get_prad_peaking #192
• Fix EFIT source tree in prad peaking computation #200

[hwietfeldt]

Closing this issue after the July 9th meeting. Eventually, we'll need to consider whether the NaNs in prad_peaking representing when the plasma core is out of view of any chords should instead be output as 0s. We should also fill in Alessandro Pau on how Prad Peaking is calculated on C-Mod.

Details of the AXUV arrays (including geometry in Fig. 1 shown below) can be found in Matt Reinke's 2008 paper in Nuclear Fusion: https://iopscience.iop.org/article/10.1088/0029-5515/48/12/125004

Here are notes from the meeting:

• C-Mod peaking factors are different to other machines due to different diagnostics. For example, DIII-D has vertical fans.
• AXUV arrays are photodiode arrays. The UV and visible radiated power might not be peaked in the center and is not Gaussian, so profile fitting doesn't make sense
• To get absolute values we would really need to perform an Abel transform but we're just interested in differences in Prad Peaking
• If we're interested in edge radiation, perhaps we could calculate a mean(edge)/mean(all) statistic
• We could consider making the "core radius" factor (currently 0.2) as a parameter w/ 0.2 as the default value.
• Note this is a bug in the Matlab scripts, which open the analysis tree instead of the efit18 tree.