Electron density evaluation from the interferometric diagnostics¶
The code calculates a phase shift between a reference sawtooth signal and mixer diode signal. The phase shift is proportional to electron density. A big part of the code are repairing mechanisms used, when the measured signal is damaged, which occurs very often.
(author: L. Lobko)
Import libraries¶
In [1]:
%matplotlib inline
import math
import os
from pathlib import Path
import numpy as np
import matplotlib.pyplot as plt
from scipy import constants
from scipy import signal as sigproc
from scipy.fft import next_fast_len
# NumPy 2 removed binary np.fromstring(); older isfreader still relies on it.
_original_fromstring = np.fromstring
def _fromstring_compat(data, dtype=float, count=-1, sep='', *, like=None):
if sep == '' and isinstance(data, (bytes, bytearray, memoryview)):
return np.frombuffer(data, dtype=dtype, count=count)
return _original_fromstring(data, dtype=dtype, count=count, sep=sep, like=like)
np.fromstring = _fromstring_compat
The data is read from ISF files, using the isfreader library. (Can be installed by "!pip3 install isfreader" command.)
In [2]:
import isfreader
In [3]:
shot_no = 53625
In [4]:
def update_db_interferometry(ne_lav_max, ne_lav_mean):
def pgnum(v):
if v is None:
return 'NULL'
fv = float(v)
if math.isnan(fv) or math.isinf(fv):
return 'NULL'
return repr(fv)
tag = 'true'
if ne_lav_max is None or ne_lav_mean is None:
tag = 'false'
else:
fm = float(ne_lav_max)
fn = float(ne_lav_mean)
if math.isnan(fm) or math.isnan(fn) or math.isinf(fm) or math.isinf(fn):
tag = 'false'
sn = int(shot_no)
mx = pgnum(ne_lav_max)
mn = pgnum(ne_lav_mean)
sql = (
f"INSERT INTO diagnostics.interferometry (shot_no, ne_lav_max, ne_lav_mean, tag_valid) "
f"SELECT {sn}, {mx}, {mn}, {tag} "
f"WHERE NOT EXISTS (SELECT 1 FROM diagnostics.interferometry WHERE shot_no = {sn}); "
f"UPDATE diagnostics.interferometry SET ne_lav_max = {mx}, ne_lav_mean = {mn}, "
f"tag_valid = {tag} WHERE shot_no = {sn}"
)
os.system(
'export PGPASSWORD=`cat /tG/production/psql_password`;psql -c "'
+ sql
+ '" -q -U golem golem_database'
)
Data loading¶
In [5]:
data_URL = "http://golem.fjfi.cvut.cz/shots/{shot_no}/Diagnostics/Interferometry/DAS_raw_data_dir/ch{ch_id}.isf"
scalars_URL = "http://golem.fjfi.cvut.cz/shots/{shot_no}/Diagnostics/PlasmaDetection/Results/{name}"
#HOTFIX
data_URL = "/tG/shots/{shot_no}/Diagnostics/Interferometry/DAS_raw_data_dir/ch{ch_id}.isf"
scalars_URL = "/tG/shots/{shot_no}/Diagnostics/PlasmaDetection/Results/{name}"
In [6]:
def shot_path(url_template, **kwargs):
return Path(url_template.format(**kwargs))
def get_scalar(shot_no, name):
return float(shot_path(scalars_URL, shot_no=shot_no, name=name).read_text().strip())
In [7]:
t_plasma_start = get_scalar(shot_no, 't_plasma_start')
t_plasma_end = get_scalar(shot_no, 't_plasma_end')
In [8]:
def load_channel(shot_no, chanel_id):
local_path = shot_path(data_URL, shot_no=shot_no, ch_id=chanel_id)
data = isfreader.read_file(str(local_path))
data[:, 0] = data[:, 0] * 1e3
return data
In [9]:
mixer = load_channel(shot_no, 1)
ref_saw = load_channel(shot_no, 3)
# phase_det = load_channel(shot_no, 4)
In [10]:
x, y = mixer.shape
f_s = x / (mixer[-1, 0] - mixer[0, 0]) # ms -> kHz
print('Sampling frequency is {} MHz.'.format(round(f_s / 1000)))
Sampling frequency is 50 MHz.
Spectral analysis of the signal
The mixer signal is a base sine wave (the envelope of the mixing) at a frequency close to 500 kHz. The reference saw-tooth frequency sweeping wave has the same base frequency, but with a trail of harmonics forming the sharp saw-tooth shape.
In [11]:
def calculate_spectrum(signal, target_dfreq=10):
nperseg = int(f_s / target_dfreq)
f, psd = sigproc.welch(signal[:, 1], fs=f_s, nperseg=nperseg, nfft=next_fast_len(nperseg))
return f, psd
In [12]:
ref_saw_f, ref_saw_psd = calculate_spectrum(ref_saw)
mixer_f, mixer_psd = calculate_spectrum(mixer)
f_base = mixer_f[mixer_psd.argmax()]
f_base_ref_saw = ref_saw_f[ref_saw_psd.argmax()]
print('The base frequency of the mixer is {} kHz.'.format(f_base))
print('The base frequency of the ref_saw is {} kHz.'.format(f_base_ref_saw))
The base frequency of the mixer is 540.0001080000216 kHz. The base frequency of the ref_saw is 540.0001080000216 kHz.
In [13]:
fig, ax = plt.subplots(dpi=150)
ax.set(xlabel='frequency [kHz]', ylabel='power spectral density [V$^2$]')
ax.plot(ref_saw_f, ref_saw_psd, label='ref_saw')
ax.plot(mixer_f, mixer_psd, label='mixer')
plt.axvline(f_base, label=f'base frequency f={f_base:.0f} kHz', color='C3')
ax.loglog()
plt.grid()
plt.legend();
Extract the baseband from the signal
The instantaneous phase and amplitude of the base signal can be inferred only for the baseband, i.e. by removing higher and lower frequencies (i.e. extracting the base sine wave from the reference saw-tooth signal).
In [14]:
base_band_hwidth = 50 # kHz
base_band_filter = sigproc.iirfilter(8, [f_base - base_band_hwidth, f_base + base_band_hwidth], fs=f_s, btype='bandpass', output='sos')
In [15]:
def freq_filter(signal, sos_filter):
signal[:, 1] = sigproc.sosfiltfilt(sos_filter, signal[:, 1])
return signal
In [16]:
mixer_filtered = freq_filter(mixer, base_band_filter)
ref_saw_filtered = freq_filter(ref_saw, base_band_filter)
cut 0.1 ms from the beginning and from the end for better signal processing
In [17]:
mixer_filtered = mixer_filtered[(mixer_filtered[:, 0] < (mixer_filtered[-1, 0] - 0.1)) & ((mixer_filtered[0, 0] + 0.1) < mixer_filtered[:, 0])]
ref_saw_filtered = ref_saw_filtered[(ref_saw_filtered[:, 0] < (ref_saw_filtered[-1, 0] - 0.1)) & ((ref_saw_filtered[0, 0] + 0.1) < ref_saw_filtered[:, 0])]
Define signal processing functions
func find_peaks - finds peaks with optimal output array
In [18]:
def find_peaks(data):
peaks_indexes, _ = sigproc.find_peaks(data[:, 1])
return np.vstack((data[peaks_indexes, 0], data[peaks_indexes, 1])).T
func initial_phase_shift - mixer and ref_saw signals are in general a bit phase shifted from each other -> it calculates "initial" phase shift and removes it
In [19]:
def initial_phase_shift(peaks, peaks_ref):
phase_mean = peaks[0, 0] - peaks_ref[0, 0]
peaks_ref[:, 0] += phase_mean
return peaks_ref
func cut_edges - cut first and last data point, which is distorted from the spectral filter
In [20]:
def cut_edges(peaks, peaks_ref):
peaks = peaks[(peaks[0, 0] < peaks[:, 0]) & (peaks[:, 0] < peaks[-1, 0])]
peaks_ref = peaks_ref[(peaks_ref[0, 0] < peaks_ref[:, 0]) & (peaks_ref[:, 0] < peaks_ref[-1, 0])]
return peaks, peaks_ref
func smooth - classic func for signal smoothing
In [21]:
def smooth(y, box_pts):
box = np.ones(box_pts) / box_pts
y_smooth = np.convolve(y, box, mode='same')
return y_smooth
func without_correction - a sum of basic operations
In [22]:
def without_correction(mixer_filtered, ref_saw_filtered):
peaks = find_peaks(mixer_filtered)
peaks_ref = find_peaks(ref_saw_filtered)
peaks_ref = initial_phase_shift(peaks, peaks_ref)
peaks, peaks_ref = cut_edges(peaks, peaks_ref)
return peaks, peaks_ref
func find_nearest - finds the nearest peak of the given one
In [23]:
def find_nearest(array, value):
array = np.asarray(array).copy()
idx = np.abs(array[:, 0] - value).argmin()
return array[idx, 0]
func calc_dphase_unchanged - calculates dphase from unrepaired data
In [24]:
def calc_dphase_unchanged(peaks, peaks_ref):
x_peaks, y_peaks = peaks.shape
x_ref_peaks, y_ref_peaks = peaks_ref.shape
dphase = np.ones((min(x_peaks, x_ref_peaks), 2))
for i in range(0, int(len(dphase))):
dphase[i, 0] = peaks[i, 0]
dphase[i, 1] = peaks[i, 0] - peaks_ref[i, 0]
dphase[:, 1] *= 2 * math.pi * f_base
return dphase
func calc_lost_phase - calculates lost phase in the signal - to define how much the signal was damaged
In [25]:
def calc_lost_phase(peaks, peaks_ref):
dphase = calc_dphase_unchanged(peaks, peaks_ref)
time_interval = 0.1 # ms
indexes = np.argwhere(dphase[:, 0] > (dphase[-1, 0] - time_interval))
data = dphase[indexes[:, 0]]
return np.average(data[:, 1]), dphase
func optimizing_cycle - defines the most probable parts of the data, where the signal was damaged and deletes the corresponding waveforms from the reference signal, because these waveforms did not travelled properly through the plasma
In [26]:
def optimizing_cycle(number_of_bad_peaks, deriv_sort, distance, peaks, peaks_ref):
bad_peaks_indexes = np.empty((0, 1))
k = 0 # help variable
l = 0 # help variable
while k < number_of_bad_peaks:
index = np.argwhere(peaks[:, 0] == deriv_sort[l, 0])
if len(bad_peaks_indexes) != 0 and (abs((index[0, 0] - find_nearest(bad_peaks_indexes, index[0, 0]))) < distance): # or (t_plasma_end < peaks[index, 0]) or (peaks[index, 0] < t_plasma_start):
l += 1
else:
bad_peaks_indexes = np.vstack((bad_peaks_indexes, index[0, 0]))
peaks_ref = np.delete(peaks_ref, index, 0)
k += 1
l += 1
return bad_peaks_indexes, peaks_ref
func repair - creates the most plausible repair of the given damaged interferometer data (probably caused by the scatter of the probing wave from the plasma, plasma instability, ...)
In [27]:
def repair(dphase, peaks, peaks_ref):
global dphase_final, bad_peaks_indexes, best_smooth_factor, best_distance
smooth_factors = np.zeros((0, 1))
distances = np.zeros((0, 1))
varieties = np.zeros((0, 1))
for smooth_factor in range(1, 20):
deriv = dphase.copy()
deriv[:, 1] = np.gradient(dphase[:, 1])
deriv[:, 1] = smooth((deriv[:, 1]), smooth_factor)
deriv_sort = deriv.copy()
deriv_sort = deriv_sort[deriv_sort[:, 1].argsort()[::-1]]
for distance in range(1, 20):
bad_peaks_indexes, repaired_ref_peaks = optimizing_cycle(number_of_bad_peaks, deriv_sort, distance, peaks, peaks_ref)
dphase_final = calc_dphase(peaks, repaired_ref_peaks)
dphase_final[:, 1] = smooth(dphase_final[:, 1], 100)
if np.min(dphase_final[:, 1]) > -1:
smooth_factors = np.vstack((smooth_factors, smooth_factor))
distances = np.vstack((distances, distance))
varieties = np.vstack((varieties, calc_curve_length(dphase_final)))
all = np.hstack((smooth_factors, distances, varieties))
if len(varieties) == 0:
print("Cannot repair interferometric data.")
else:
varieties_min = varieties.argmin()
best_smooth_factor = int(all[int(varieties_min), 0])
best_distance = int(all[int(varieties_min), 1])
deriv = dphase.copy()
deriv[:, 1] = np.gradient(dphase[:, 1])
deriv[:, 1] = smooth((deriv[:, 1]), best_smooth_factor)
deriv_sort = deriv.copy()
deriv_sort = deriv_sort[deriv_sort[:, 1].argsort()[::-1]]
bad_peaks_indexes, repaired_ref_peaks = optimizing_cycle(number_of_bad_peaks, deriv_sort, best_distance, peaks, peaks_ref)
dphase_final = calc_dphase(peaks, repaired_ref_peaks)
return dphase_final, bad_peaks_indexes, best_smooth_factor, best_distance
func calc_dphase - calculates dphase
In [28]:
def calc_dphase(mixer_peaks, repaired_ref_peaks):
x_peaks, y_peaks = mixer_peaks.shape
x_ref_peaks, y_ref_peaks = repaired_ref_peaks.shape
dphase_final = np.ones((min(x_peaks, x_ref_peaks), 2))
dphase_final[:, 0] = repaired_ref_peaks[0:len(dphase_final), 0]
dphase_final[:, 1] = (mixer_peaks[0:len(dphase_final), 0] - repaired_ref_peaks[0:len(dphase_final), 0]) * 2 * math.pi * f_base
# for i in range(0, int(len(dphase_final))):
# dphase_final[i, 0] = repaired_ref_peaks[i, 0]
# dphase_final[i, 1] = mixer_peaks[i, 0] - repaired_ref_peaks[i, 0]
# dphase_final[:, 1] *= 2 * math.pi * f_base
return dphase_final
func calc_curve_length - calculates the length of the dphase curve as a key parameter to decide, whether the repair is good enough
In [29]:
def calc_curve_length(dphase_final):
x, y = dphase_final.shape
length = 0
part_length = np.sqrt(((dphase_final[1:-1, 0] - dphase_final[0:-2, 0]) ** 2) + (dphase_final[1:-1, 1] - dphase_final[0:-2, 1]) ** 2)
# for i in range(1, x):
# length += np.sqrt(((dphase_final[i, 0] - dphase_final[i - 1, 0]) ** 2) + (dphase_final[i, 1] - dphase_final[i - 1, 1]) ** 2)
# # distance = np.sqrt((xy_inside[:, 0] - xy_inside[int(((n_i) * (n_i) - 1) / 2), 0]) ** 2 + (xy_inside[:, 1] - xy_inside[int(((n_i) * (n_i) - 1) / 2), 1]) ** 2)
return np.sum(part_length)
Final signal processing cycle¶
In [30]:
# prepare data without any correction
peaks_no_corr, peaks_ref_no_corr = without_correction(mixer_filtered, ref_saw_filtered)
# calculate lost phase and dphase_zero - dphase without any correction
lost_phase, dphase_zero = calc_lost_phase(peaks_no_corr, peaks_ref_no_corr)
repaired_discharge = False
unrepairable_discharge = False
number_of_bad_peaks = round(lost_phase / (2 * math.pi))
# make data repair, if it is necessary
if 1 <= number_of_bad_peaks:
try:
dphase_final, bad_peaks_indexes, best_smooth_factor, best_distance = repair(
dphase_zero, peaks_no_corr, peaks_ref_no_corr
)
dphase_final[:, 1] = smooth(dphase_final[:, 1], 100)
repaired_discharge = True
except Exception as exc:
unrepairable_discharge = True
print(f"Cannot repair data: {exc}")
Estimation of the electron density
The ordinary wave (O-mode) with a carrier frequency $\omega$ traveling through a collisionless plasma with the plasma frequency $\omega_{p} = \sqrt{\frac{n_e e^2}{\epsilon_0 m_e}}$ has a refractive index $$N_O=\sqrt{1-\frac{\omega_p^2}{\omega^2}}$$ Under the assumption that the carrier wave frequency is much larger than the plasma frequency $\omega>>\omega_p$ this formula can be expanded into a Taylor series as $$N_O\approx 1-\frac{\omega_p^2}{2\omega^2}$$ A wave traveling through a medium with a refractive index $N(l)$ accumulates a total phase shift $\varphi = \frac{\omega}{c} \int N(l) \mathrm{d}l$. Therefore, in comparison to a wave traveling in vacuum (or clear air) with $N\approx 1$, the wave traveling through the plasma over a chord with length $L$ accumulates a relative phase shift of $$\Delta \varphi = \frac{e^2}{2\omega c\epsilon_0 m_e}\int\limits_L n_e(l) \mathrm{d}l$$ Therefore, it is possible to estimate the line-averaged density $\bar n_e = \frac{1}{L} \int\limits_L n_e(l) \mathrm{d}l$ from the detected phase shift between the reference and mixer signals.
In [31]:
omega_carrier = 2 * np.pi * 71e9 # 71 GHz microwave generator
a = 0.085 # limiter radius [m]
L = 2 * a # central plasma chord estimate
prop_const = constants.elementary_charge ** 2 / (2 * omega_carrier * constants.speed_of_light * constants.epsilon_0 * constants.m_e)
ne_lav - final line-averaged electron density
In [32]:
if repaired_discharge:
ne_lav = dphase_final.copy()
ne_lav[:, 1] = ne_lav[:, 1] * (1 / (prop_const * L))
# ne_lav = ne_lav[(ne_lav[:, 0] >= 0) & (ne_lav[:, 0] <= (t_plasma_end + 5))]
In [33]:
dphase_zero_raw = dphase_zero.copy()
dphase_zero[:, 1] = smooth(dphase_zero[:, 1], 100)
ne_lav_zero - final line-averaged electron density without any correction
In [34]:
ne_lav_zero = dphase_zero.copy()
ne_lav_zero[:, 1] = ne_lav_zero[:, 1] * (1 / (prop_const * L))
# ne_lav_zero = ne_lav_zero[(ne_lav_zero[:, 0] >= 0) & (ne_lav_zero[:, 0] <= (t_plasma_end + 5))]
ne_lav_zero_raw - final line-averaged electron density without any correction and without any smoothing - used to see more clearly, where the data was damaged and repaired in the last figure
In [35]:
ne_lav_zero_raw = dphase_zero_raw.copy()
ne_lav_zero_raw[:, 1] = ne_lav_zero_raw[:, 1] * (1 / (prop_const * L))
# ne_lav_zero_raw = ne_lav_zero_raw[(ne_lav_zero_raw[:, 0] >= 0) & (ne_lav_zero_raw[:, 0] <= (t_plasma_end + 5))]
Comparison of the final line-averaged electron density between the damaged and repaired signal¶
In [36]:
if repaired_discharge:
fig, ax = plt.subplots(dpi=200)
ax.plot(ne_lav_zero[:, 0], ne_lav_zero[:, 1] / 1e18, label='damaged signal', linestyle='dotted', color='red')
ax.plot(ne_lav[:, 0], ne_lav[:, 1] / 1e18, label='repaired signal', color='deepskyblue')
if (np.max(ne_lav_zero[:, 1] / 1e18) < 0.1):
ax.set_ylim([-0.1, 5])
else:
for t in (t_plasma_start, t_plasma_end):
plt.axvline(t, color='k', linestyle='--')
ax.set(xlabel='time [ms]', ylabel='$\\bar n_e$ [10$^{18}$ m$^{-3}$]')
plt.title('Line-averaged electron density (repaired waveforms: {})'.format(number_of_bad_peaks))
plt.legend()
plt.grid()
plt.show()
Final figure of the temporal evolution of the line-averaged electron density¶
In [37]:
fig, ax = plt.subplots(dpi=200)
if repaired_discharge:
ax.plot(ne_lav[:, 0], ne_lav[:, 1] / 1e18, label='$\\bar n_e$', color='deepskyblue')
else:
ax.plot(ne_lav_zero[:, 0], ne_lav_zero[:, 1] / 1e18, label='$\\bar n_e$', color='deepskyblue')
if (np.max(ne_lav_zero[:, 1] / 1e18) < 0.1):
ax.set_ylim([-0.1, 5])
else:
for t in (t_plasma_start, t_plasma_end):
plt.axvline(t, color='k', linestyle='--')
ax.set(xlabel='time [ms]', ylabel='$\\bar n_e$ [10$^{18}$ m$^{-3}$]')
if unrepairable_discharge:
plt.title('Line-averaged electron density (!!!UNREPAIRED waveforms: {}!!!)'.format(number_of_bad_peaks))
else:
plt.title('Line-averaged electron density (repaired waveforms: {})'.format(number_of_bad_peaks))
plt.legend()
plt.grid()
plt.savefig('icon-fig.png')
plt.show()
FOR TESTING (orange lines are repaired locations)¶
In [38]:
if repaired_discharge:
fig, ax = plt.subplots(dpi=200)
ax.plot(ne_lav_zero_raw[:, 0], ne_lav_zero_raw[:, 1] / 1e18, label='damaged signal', color='red')
ax.plot(ne_lav[:, 0], ne_lav[:, 1] / 1e18, label='$\\bar n_e$', color='deepskyblue')
x, y = bad_peaks_indexes.shape
for i in range(0, x):
plt.axvline(peaks_no_corr[int(bad_peaks_indexes[i, 0]), 0], color='orange', linestyle='--')
ax.set(xlabel='time [ms]', ylabel='$\\bar n_e$ [10$^{18}$ m$^{-3}$]')
plt.title('best smooth factor: {}, best distance: {}, repaired waveforms: {}'.format(best_smooth_factor, best_distance, number_of_bad_peaks))
plt.legend()
plt.grid()
plt.show()
Save data¶
file ne_lav.csv - the final line-averaged electron density data file ne_lav_max.txt - max value of the line-averaged electron density file ne_lav_mean.txt - mean value of the line-averaged electron density
In [39]:
if repaired_discharge:
np.savetxt('ne_lav.csv', ne_lav, delimiter=",") # repaired data
np.savetxt('ne_lav_unrepaired.csv', ne_lav_zero, delimiter=",")
ne_lav_plasma = ne_lav[(0.1 < ne_lav[:, 1] / 1e18)]
ne_lav_mean_v = getattr(ne_lav_plasma[:, 1], 'mean')()
ne_lav_max_v = getattr(ne_lav_plasma[:, 1], 'max')()
with open('ne_lav_mean', 'w') as f:
f.write(f'{ne_lav_mean_v:.3e}')
with open('ne_lav_max', 'w') as f:
f.write(f'{ne_lav_max_v:.3e}')
update_db_interferometry(ne_lav_max_v, ne_lav_mean_v)
else:
np.savetxt('ne_lav.csv', ne_lav_zero, delimiter=",")
if unrepairable_discharge:
ne_lav_mean_v = float('NaN')
ne_lav_max_v = float('NaN')
with open('ne_lav_mean', 'w') as f:
f.write(f'{ne_lav_mean_v:.3e}')
with open('ne_lav_max', 'w') as f:
f.write(f'{ne_lav_max_v:.3e}')
update_db_interferometry(ne_lav_max_v, ne_lav_mean_v)
else:
if np.max(ne_lav_zero[:, 1]/1e18) < 0.1:
ne_lav_mean_v = 0.0
ne_lav_max_v = 0.0
with open('ne_lav_mean', 'w') as f:
f.write(f'{ne_lav_mean_v:.3e}')
with open('ne_lav_max', 'w') as f:
f.write(f'{ne_lav_max_v:.3e}')
update_db_interferometry(ne_lav_max_v, ne_lav_mean_v)
else:
ne_lav_plasma = ne_lav_zero[(0.1 < ne_lav_zero[:, 1] / 1e18)]
ne_lav_mean_v = getattr(ne_lav_plasma[:, 1], 'mean')()
ne_lav_max_v = getattr(ne_lav_plasma[:, 1], 'max')()
with open('ne_lav_mean', 'w') as f:
f.write(f'{ne_lav_mean_v:.3e}')
with open('ne_lav_max', 'w') as f:
f.write(f'{ne_lav_max_v:.3e}')
update_db_interferometry(ne_lav_max_v, ne_lav_mean_v)