Double rake probe spotlight

The following scripts introduces basic aspects of the double rake probe measurements.

The double rake probe diagnostic

This is the double rake probe:

Its wiki page is located here. The double rake probe has two rows of pins ("rakes"), out of which the first 6 in each row are connected to the data acquisition system. The pin signals are called DRP-R1 to DRP-R6 and DRP-L1 to DRP-L6.

Like any Langmuir probe, the rake probe pins can be electrically insulated (measuring the floating potential $V_f$), biased to a negative voltage -100 V (measuring the ion saturated current $I_{sat}$), their biasing voltage may be swept (measuring $V_f$, $I_{sat}$ and the electron temperature $T_e$ if the $I$-$V$ characteristic is fitted with an exponential function) or they can perform other measurements. During this campaign, all the pins measure the ion saturated current $I_{sat}$. This can be processed to gain information about the plasma turbulent transport.

Import basic libraries

To visualise and process the double rake probe data, first we import basic libraries, Numpy and Matplotlib. The line %matplotlib notebook enables the drawn figures to be interactive.

%matplotlib inline
import numpy as np
import matplotlib.pyplot as pl

#Number of the test discharge; 0 means the last one
shot_no = 35881 #35041 #test discharge
shot = shot_no

Load double rake probe data

The data directory of the double rake probe is http://golem.fjfi.cvut.cz/shots/{shot}/Diagnostics/DoubleRakeProbe/, where {shot} stands for the discharge number. Note that in different sessions, the data may be saved elsewhere and it might be needed to update the URL variable in the following text. The data directory may be located from the individual shot homepage, and tends to be the same within a single session.

In the following, we write a function to download the rake probe data.

from urllib.error import HTTPError # recognise the error stemming from missing data
import pandas as pd # for reading csv files

#Define an exception which will be raised if the data is missing and stop the notebook execution
class StopExecution(Exception):
    def _render_traceback_(self):
        pass

def get_data(shot, identifier):
    URL = "http://golem.fjfi.cvut.cz/shots/{shot}/Diagnostics/DoubleRakeProbe/{identifier}.csv"
    url = URL.format(shot=shot, identifier=identifier)
    try:
        df = pd.read_csv(url, names=['time', identifier], index_col='time')
    except HTTPError:
        print('File not found at %s . Aborting notebook execution.' % url)
        raise StopExecution
    t = np.array(df.index)
    data = np.transpose(np.array(df))[0]
    return t, data/46.7

Notice that in the last line, the data values are divided by 46.7. This is the constant which converts the raw DAS signal to the ion saturated current in A. Its physical meaning is the resistance of the measuring resistor in the "silver box".

Load the discharge beginning and end

Prior to plotting the double rake probe data, we'll write a little function which loads the time of the discharge beginning (plasma formation) and end (recombination). We'll use it to set the X axis limits later. Notice that t1 and t2 are in ms.

from urllib.request import urlopen

def get_file(url, shot, silent=False):
    URL = 'http://golem.fjfi.cvut.cz/shots/%i/%s' % (shot, url)
    if not silent:
        print(URL)
    f = urlopen(URL)
    try:
        return np.loadtxt(f)
    except ValueError: # the data couldn't be converted to a row of numbers
        return np.array([np.nan])

t1 = float(get_file('Diagnostics/BasicDiagnostics/Results/t_plasma_start', shot))
t2 = float(get_file('Diagnostics/BasicDiagnostics/Results/t_plasma_end', shot))
print('Discharge %i: %.1f ms - %.1f ms' % (shot, t1, t2))

http://golem.fjfi.cvut.cz/shots/35881/Diagnostics/BasicDiagnostics/Results/t_plasma_start
http://golem.fjfi.cvut.cz/shots/35881/Diagnostics/BasicDiagnostics/Results/t_plasma_end
Discharge 35881: 3.6 ms - 19.0 ms

Plot all the signals from the double rake probe pins

Next, we load the double rake probe data for the current discharge and plot them.

#Load the signals of all the rake probe pins
DRP_R = []
DRP_L = []
for i in range(1, 7):
    t, data = get_data(shot, 'DRP-R%i' % i)
    DRP_R.append(data)
for i in range(1, 7):
    t, data = get_data(shot, 'DRP-L%i' % i)
    DRP_L.append(data)
DRP_R = np.array(DRP_R)
DRP_L = np.array(DRP_L)

#Plot the rake probe pin signals
fig, axs = pl.subplots(6, 2, num=('All rake probe signals in the current discharge'), figsize=(9,9),
                       sharex=True, sharey=True)
for i in range(0, 6):
    axs[i, 0].plot(t*1000, DRP_L[i]*1000, label='DRP-L%i' % (i+1))
for i in range(0, 6):
    axs[i, 1].plot(t*1000, DRP_R[i]*1000, label='DRP-R%i' % (i+1))
fig.tight_layout()
fig.subplots_adjust(hspace=0, wspace=0.3)
for i in range(6):
    for j in range(2):
        axs[i,j].legend(loc=2)
        axs[i,j].grid(True)
        axs[i,j].set_xlim(t1-3, t2+3)
        axs[i,j].set_ylim(-0.5, 4)
        axs[i,j].axhline(c='k')
        axs[i,j].set_ylabel('$I_{sat}$ [mA]')
pl.savefig("icon-fig.png")

Notice how the signal wanes on the pins with a higher number. This is because they are not as deep in the plasma column and they are surrounded by plasma of lower temperature and density, hence lower $I_{sat}$.

Load the rake probe position

#R_probe = 70 #mm
R_probe = get_file('Diagnostics/DoubleRakeProbe/Parameters/r_first_tip', shot)
print('R_first_pin = %i mm' % R_probe)

http://golem.fjfi.cvut.cz/shots/35881/Diagnostics/DoubleRakeProbe/Parameters/r_first_tip
R_first_pin = 85 mm

Plot histogram of the $I_{sat}$ fluctuations

The basic method of investigating the presence of turbulent structures in the edge plasma is plotting the $I_{sat}$ histogram. A histogram is an approximation of the distribution function - the probability of a certain value occurring in the signal. It is plotted by the function pl.hist().

In the following, we plot the histogram of the $I_{sat}$ fluctuations in the test discharge shot. The fluctuations are sampled during a 1ms window from 8 ms to 9 ms (approx. middle of the discharge).

# Choose data of the sample
data = DRP_R[0] #deepest pin in the R row
sample = data[(0.008 < t) & (t < 0.009)]

#Plot the histogram
fig = pl.figure('Histogram of I_sat fluctuations', figsize=(6,4))
fig.tight_layout()
pl.hist(sample*1000, bins=30)
pl.xlabel('$I_{sat}$ [mA]')
pl.ylabel('value count')

#Plot the mean with a thick black line
pl.axvline(sample.mean(), c='k', lw=2)

<matplotlib.lines.Line2D at 0x7f40f68985b0>

Notice that the histogram might be asymmetrical.

Calculate the distribution skewness

Skewness $S$ is the third statistical moment of a distribution function, the first two being the mean and the standard deviation. If the distribution is symmetric, its skewness is zero. For instance, the skewness of the normal (Gaussian) distribution is zero. The skewness of the $I_{sat}$ fluctuations in this case, however, isn't zero and the distribution is not normal. It has an abundance of large positive events (blobs) while lacking the same number of large negative events. This is evidence of the presence of blobs in the SOL.

In the following, we calculate the skewness using the function scipy.stats.skew.

from scipy.stats import skew

S = skew(sample)
print('The sample skewness is %.2f.' % S)

The sample skewness is 0.07.

Plot the skewness profile from all probes

In the following, we use all the probe pins and plot the profile of the $I_{sat}$ fluctuations skewness. The distance between the pins is 2.5 mm.

#Prepare the variables
S = []
T = np.arange(np.ceil(t1+0.5), t2-0.5) #windows centres, spaced by 1 ms
pins = ['L'+str(i) for i in range(1,7)] + ['R'+str(i) for i in range(1,7)]
r = [0, 2.5, 5, 7.5, 10, 12.5] * 2
r = np.array(r) + R_probe #mm
colour_indices = np.linspace(0, 1, T.size)

#Iterate over the pins and time windows and calculate skewness
for pin in pins:
    t, DRP = get_data(shot, 'DRP-'+pin)
    t *= 1000 #ms
    s = []
    for tau in T:
        sample = DRP[(tau-0.5 < t) & (t < tau+0.5)]
        s.append(skew(sample))
    S.append(s)
S = np.array(S)
    
#Plot the skewness profile
fig = pl.figure('Skewness profile', figsize=(8,6))
for i in range(T.size):
    pl.plot(r[:6], S[:6,i], marker='o', color=pl.cm.copper_r(colour_indices[i]),
            lw=0.5, label='%i ms' % T[i])
    pl.plot(r[6:], S[6:,i], marker='o', color=pl.cm.copper_r(colour_indices[i]), lw=0.5)
pl.xlabel('$r$ [mm]')
pl.ylabel('$S$')
pl.axhline(0, c='k')
pl.legend()
pl.grid(True)