%matplotlib inline
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
import math
shot_no = 40485
The data is read from ISF files, using the isfreader
library. It is installed bellow automatically if not already installed.
!pip3 install isfreader
Requirement already satisfied: isfreader in /opt/anaconda3/lib/python3.8/site-packages (1.0.0) Requirement already satisfied: numpy>=1.9.1 in /opt/anaconda3/lib/python3.8/site-packages (from isfreader) (1.19.2)
import isfreader
ds = np.DataSource('/tmp') # temporary storage for downloaded files
data_URL = 'http://golem.fjfi.cvut.cz/shots/{shot_no}/Diagnostics/Interferometry/LukLob/DAS_raw_data_dir/ch{ch_id}.isf'
scalars_URL = 'http://golem.fjfi.cvut.cz/shots/{shot_no}/Diagnostics/BasicDiagnostics/DetectPlasma/Results/{name}'
def get_scalar(shot_no, name):
return float(ds.open(scalars_URL.format(shot_no=shot_no, name=name)).read())
t_plasma_start = get_scalar(shot_no, 't_plasma_start')
t_plasma_end = get_scalar(shot_no, 't_plasma_end')
is_plasma = get_scalar(shot_no, 'is_plasma')
def load_channel(shot_no, chanel_id):
fname = ds.open(data_URL.format(shot_no=shot_no, ch_id=chanel_id)).name
data = isfreader.read_file(fname)
data[:, 0] = data[:, 0] * 1e3
return data
mixer = load_channel(shot_no, 1)
ref_saw = load_channel(shot_no, 3)
x, y = mixer.shape
f_s = x / (mixer[-1, 0] - mixer[0, 0]) # ms -> kHz
f_s
250000.05000001
The mixer signals 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.
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
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()]
f_base # [kHz]
540.0001080000216
f_base_ref_saw
540.0001080000216
fig, ax = plt.subplots()
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.legend()
<matplotlib.legend.Legend at 0x7f4ec00e6880>
The difference between the phase of the reference and mixer signal is proportional to the line-integrated electron density. Therefore, it is necessary to obtain the "phase" of both signals.
The instantaneous phase and amplitude of the base signal can be inferred only for the base band, i.e. by removing higher and lower frequencies (i.e. extracting the base sine wave from the reference saw-tooth signal)
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')
def freq_filter(signal, sos_filter):
signal[:, 1] = sigproc.sosfiltfilt(sos_filter, signal[:, 1])
return signal
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
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])]
def find_peaks(data):
peaks_indexes, _ = sigproc.find_peaks(data[:, 1])
return np.vstack((data[peaks_indexes, 0], data[peaks_indexes, 1])).T
def initial_phase_shift(peaks, peaks_ref, mixer_filtered):
phase_mean = peaks[0, 0] - peaks_ref[0, 0]
mixer_filtered[:, 0] -= phase_mean
peaks[:, 0] -= phase_mean
return peaks, mixer_filtered
def differences(signal):
data = signal.copy()
x, y = data.shape
for i in range(1, int(x) - 1):
data[i, 0] = signal[i, 0] - signal[i - 1, 0]
data[i, 1] = signal[i + 1, 0] - signal[i, 0]
return np.hstack((signal, data))
def cut_edges(diff_peaks_mixer, diff_peaks_ref, peaks, peaks_ref):
diff_peaks_mixer = diff_peaks_mixer[(diff_peaks_mixer[0, 0] < diff_peaks_mixer[:, 0]) & (diff_peaks_mixer[:, 0] < diff_peaks_mixer[-1, 0])]
diff_peaks_ref = diff_peaks_ref[(diff_peaks_ref[0, 0] < diff_peaks_ref[:, 0]) & (diff_peaks_ref[:, 0] < diff_peaks_ref[-1, 0])]
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 diff_peaks_mixer, diff_peaks_ref, peaks, peaks_ref
def remove_basic_diff(diff_peaks_mixer, diff_peaks_ref):
basic_diff = np.mean(diff_peaks_ref[:, 2])
diff_peaks_mixer[:, 2] -= basic_diff
diff_peaks_ref[:, 2] -= basic_diff
return diff_peaks_mixer, diff_peaks_ref
def sort_differencies(diff_peaks_mixer):
data = diff_peaks_mixer.copy()
x, y = data.shape
for i in range(0, x):
if data[i, 2] < 0:
data[i, 2] *= 0.9
data[:, 2] = abs(data[:, 2])
return data[data[:, 2].argsort()[::-1]]
def smooth(y, box_pts):
box = np.ones(box_pts)/box_pts
y_smooth = np.convolve(y, box, mode='same')
return y_smooth
def without_correction(mixer_filtered, ref_saw_filtered):
mixer_data = mixer_filtered.copy()
ref_saw_data = ref_saw_filtered.copy()
peaks = find_peaks(mixer_data)
peaks_ref = find_peaks(ref_saw_data)
peaks, mixer_data = initial_phase_shift(peaks, peaks_ref, mixer_data)
diff_peaks_mixer = differences(peaks)
diff_peaks_ref = differences(peaks_ref)
diff_peaks_mixer, diff_peaks_ref, peaks, peaks_ref = cut_edges(diff_peaks_mixer, diff_peaks_ref, peaks, peaks_ref)
diff_peaks_mixer, diff_peaks_ref = remove_basic_diff(diff_peaks_mixer, diff_peaks_ref)
return mixer_data, peaks, peaks_ref, diff_peaks_mixer, diff_peaks_ref
def basic_cycle(mixer_filtered, ref_saw_filtered):
peaks = find_peaks(mixer_filtered)
peaks_ref = find_peaks(ref_saw_filtered)
peaks, mixer_filtered = initial_phase_shift(peaks, peaks_ref, mixer_filtered)
diff_peaks_mixer = differences(peaks)
diff_peaks_ref = differences(peaks_ref)
diff_peaks_mixer, diff_peaks_ref, peaks, peaks_ref = cut_edges(diff_peaks_mixer, diff_peaks_ref, peaks, peaks_ref)
diff_peaks_mixer, diff_peaks_ref = remove_basic_diff(diff_peaks_mixer, diff_peaks_ref)
sorted_diff_mixer = sort_differencies(diff_peaks_mixer)
return peaks, peaks_ref, mixer_filtered, diff_peaks_mixer, diff_peaks_ref, sorted_diff_mixer
def find_nearest(array, value):
array = array.copy()
x, y = array.shape
array_data = np.zeros((x, 1))
for i in range(0, x):
array_data[i, 0] = array[i, 0]
idx = (np.abs(array_data[:, 0] - value)).argmin()
return array_data[idx]
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] # TODO: zmenit na peaks_ref?
dphase[i, 1] = peaks[i, 0] - peaks_ref[i, 0]
return dphase
def calc_lost_phase(peaks, peaks_ref, diff_peaks_mixer, diff_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 2*math.pi*f_base*np.average(data[:, 1]), dphase
def repair(mixer_filtered, ref_saw_filtered, lost_phase, dphase):
peaks, peaks_ref, mixer_filtered, diff_peaks_mixer, diff_peaks_ref, sorted_diff_mixer = basic_cycle(mixer_filtered, ref_saw_filtered)
# make repairing cycle
number_of_bad_peaks = round(lost_phase/(2*math.pi))
remainder = lost_phase - 2*math.pi*number_of_bad_peaks
x_peaks_mixer, y_peaks_mixer = peaks.shape
x_peaks_ref, y_peaks_ref = peaks_ref.shape
deriv = dphase.copy()
deriv[:, 1] = np.gradient(dphase[:, 1])
deriv[:, 1] = smooth((deriv[:, 1]), 10)
deriv_sort = deriv.copy()
# deriv_sort[:, 1] = abs(deriv_sort[:, 1])
deriv_sort = deriv_sort[deriv_sort[:, 1].argsort()[::-1]]
bad_peaks_indexes = np.empty((0, 1))
k = 0
l = 0
while k < number_of_bad_peaks:
index = np.argwhere(diff_peaks_mixer[:, 0] == deriv_sort[l, 0])
if (l > 0) and (abs((index[0, 0] - find_nearest(bad_peaks_indexes, index[0, 0]))) < 10):
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
if number_of_bad_peaks == 0:
for i in range(x_peaks_mixer):
if remainder > math.pi:
index = np.argwhere(diff_peaks_mixer[:, 0] == deriv_sort[i, 0])
peaks_ref = np.delete(peaks_ref, index, 0)
bad_peaks_indexes = np.vstack((bad_peaks_indexes, index[0, 0]))
return peaks, peaks_ref, bad_peaks_indexes, diff_peaks_mixer, deriv
def calc_dphase(peaks, peaks_ref, repaired_mixer_peaks, repaired_ref_peaks):
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))
dphase[:, 1] *= 2*math.pi*f_base
for i in range(0, int(len(dphase))):
dphase[i, 0] = peaks_ref[i, 0]
dphase[i, 1] = peaks[i, 0] - peaks_ref[i, 0]
x_peaks, y_peaks = repaired_mixer_peaks.shape
x_ref_peaks, y_ref_peaks = repaired_ref_peaks.shape
dphase1 = np.ones((min(x_peaks, x_ref_peaks), 2))
for i in range(0, int(len(dphase1))):
dphase1[i, 0] = repaired_ref_peaks[i, 0]
dphase1[i, 1] = repaired_mixer_peaks[i, 0] - repaired_ref_peaks[i, 0]
dphase1[:, 1] *= 2*math.pi*f_base
return dphase, dphase1
mixer_filtered_no_corr, peaks_no_corr, peaks_ref_no_corr, diff_peaks_mixer_no_corr, diff_peaks_ref_no_corr = without_correction(mixer_filtered, ref_saw_filtered)
lost_phase, dphase_zero = calc_lost_phase(peaks_no_corr, peaks_ref_no_corr, diff_peaks_mixer_no_corr, diff_peaks_ref_no_corr)
repaired_mixer_peaks, repaired_ref_peaks, bad_peaks_indexes, diff_peaks_mixer, deriv = repair(mixer_filtered, ref_saw_filtered, lost_phase, dphase_zero)
dphase, dphase1 = calc_dphase(peaks_no_corr, peaks_ref_no_corr, repaired_mixer_peaks, repaired_ref_peaks)
dphase1[:, 1] = smooth(dphase1[:, 1], 100)
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.
omega_carrier = 2 * np.pi * 71e9 # 71 GHz microwave generator
a = 0.085 # limiter radius
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 = dphase1.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))]
fig, ax = plt.subplots(figsize=(6, 4.5), dpi=300)
ax.plot(ne_lav[:, 0], ne_lav[:, 1]/1e18)
for t in (t_plasma_start, t_plasma_end):
plt.axvline(t, color='k', linestyle='--')
ax.set(xlabel='time [ms]', ylim=(0, None), ylabel='$\\bar n_e$ [10$^{18}$ m$^{-3}$]')
plt.savefig('icon-fig.png')
plt.show()
np.savetxt('ne_lav.csv', ne_lav, delimiter=",")
ne_lav_plasma = ne_lav[((t_plasma_start) < ne_lav[:, 0]) & (ne_lav[:, 0] < (t_plasma_end))]
for m in ('mean', 'max'):
v = getattr(ne_lav_plasma[:, 1], m)()
with open(f'ne_lav_{m}', 'w') as f: f.write(f'{v:.3e}')