---
jupyter:
  jupytext:
    formats: ipynb,md
    text_representation:
      extension: .md
      format_name: markdown
      format_version: '1.2'
      jupytext_version: 1.6.0
  kernelspec:
    display_name: Python 3
    language: python
    name: python3
---

<!-- #region hide_input=false -->
# Load libraries + setup functions
<!-- #endregion -->

```python hide_input=true
import warnings
warnings.filterwarnings('ignore')

import sys
import os
from pyrpl import Pyrpl
import numpy as np
import time
import matplotlib.pyplot as plts

import IPython
import ipywidgets as widgets

from bokeh.plotting import figure, show
from bokeh.io import output_notebook, push_notebook
from bokeh.layouts import column, row
from bokeh.models import ColumnDataSource, Toggle, Range1d

from scipy.signal import welch

from datetime import datetime

#Arduino Serial connection
import serial

# Red Pitaya Stuff

HOSTNAME = "rp-f06897.local"
p = Pyrpl(hostname=HOSTNAME, config='test',gui=False)
r = p.rp
s = r.scope

# Some good default settings
s.input1 = 'in1'
s.input2 = 'in2'
s.decimation = 1024 # or s.duration =0.01
s.average = True
s.trigger_source = 'immediately'

# For the gui updates during the live view loops
ipython = IPython.get_ipython()

# For Bokeh
output_notebook()

# The full scale ranges
I_full_scale = 520 #uA
V_full_scale = 800 #uV

# An autoscaling routine
def auto_scale(r, d):
    d_max = np.max(d)
    d_min = np.min(d)
    d_range = d_max - d_min 
    if d_range < 0.01:
        d_range = 0.01
    r.start = d_min - d_range*0.1
    r.end = d_max + d_range*0.1

def get_IV_t(decimation = 1024):
    s.decimation = decimation
    s._start_acquisition()
    time.sleep(s.duration)
    I,V = s._get_curve()
    I *= 20 # we have RP in HV mode
    V *= 20 # we have RP in HV mode
    I /= 0.91 # scaling correction for my RP
    V /= 0.91 # scaling correction for my RP
    I /= 1e4 # Conversion factor Mr SQUID 10000 A/V
    V /= 1e4 # Conversion factor Mr SQuID 10000 V/V
    t = s.times
    return I*1e6,V*1e6,t
```

<!-- #region hide_input=false -->
# Your information

Please update the practicum round
<!-- #endregion -->

```python hide_input=false
academic_year = "2020_2021"
practicum_round = 4
```

Run this code to create the folder where your data will be saved

```python
folder = academic_year + "_round_%02d" % practicum_round
if not os.path.exists(folder):
    os.makedirs(folder)
```

<!-- #region hide_input=false -->
# IV vs Time Live Display
<!-- #endregion -->

```python hide_input=true
I = []
V = []
t = []

def update_display():
    global I,V,t
    I,V,t = get_IV_t()
    I_avg = np.average(I)
    I_ptp = np.max(I) - np.min(I)
    p1.title.text = "Current          Average = %.2f uA     Peak-to-Peak %.2f uA" % (I_avg, I_ptp)
    V_avg = np.average(V)
    V_ptp = np.max(V) - np.min(V)
    p2.title.text = "Voltage          Average = %.2f uV     Peak-to-Peak %.2f uV" % (V_avg, V_ptp)
    R_ptp.value = "%.2f" % (V_ptp / I_ptp)
    R_avg.value = "%.2f" % (V_avg / I_avg)
    source1.data = dict(x=t, y=I)
    source2.data = dict(x=t, y=V)
    if scale_mode.value == "Autoscale":# and False:
        auto_scale(p1.y_range,I)
        auto_scale(p2.y_range,V)
    elif scale_mode.value == "Full range":
        p1.y_range.start = -I_full_scale
        p1.y_range.end = I_full_scale
        p2.y_range.start = -V_full_scale
        p2.y_range.end = V_full_scale
        
range1 = Range1d()
range2 = Range1d()
source1 = ColumnDataSource()
source2 = ColumnDataSource()

p1 = figure(title="Current", plot_height=300, plot_width=900,toolbar_location=None, y_range = range1)
#p1 = figure(title="Current", plot_height=300, plot_width=900,toolbar_location=None)
p1.yaxis.axis_label = 'Curret (uA)'

p2 = figure(title="Voltage", plot_height=300, plot_width=900,toolbar_location=None, y_range = range2)
#p2 = figure(title="Voltage", plot_height=300, plot_width=900,toolbar_location=None)

p2.xaxis.axis_label = 'Time (s)'
p2.yaxis.axis_label = 'Voltage (uV)'

p1.line('x', 'y', source=source1)
p2.line('x', 'y', source=source2)

style = {'description_width': 'initial'}

stop_button = widgets.ToggleButton(description='Stop')
pause_button = widgets.ToggleButton(description='Pause')

def save_data(w):
    fmt = f"mrsquid_{mrsquid_mode.value}_%Y-%m-%d-%H_%M_%S.dat"
    outname = folder + "/" + datetime.now().strftime(fmt)
    np.savetxt(outname, np.array([t,I,V]).T)
    filename.value = outname
    
save_button = widgets.Button(description='Save data')
save_button.on_click(save_data)
filename = widgets.HTML(description="Last filename: ",style=style)
filename.value = "(none)"

scale_mode = widgets.RadioButtons(options=['Autoscale', 'Full range', 'Freeze range'])
mrsquid_mode = widgets.RadioButtons(options=['V-I', 'V-Phi'])

R_ptp = widgets.Text(description="V/I ptp (ohms)",style=style)
R_avg = widgets.Text(description="V/I avg (ohms)",style=style)
R_ptp.layout.width = '150px'
R_avg.layout.width = '150px'

update_display()
target = show(column(p1,p2), notebook_handle=True)
rows = []
rows.append(widgets.HBox([stop_button,pause_button,scale_mode, mrsquid_mode]))
rows.append(widgets.HBox([save_button, filename]))
rows.append(widgets.HBox([R_ptp, R_avg]))
display(widgets.VBox(rows))

while True:
    ipython.kernel.do_one_iteration()
    if not pause_button.value:
        update_display()
    ipython.kernel.do_one_iteration()
    if stop_button.value:
        break
    push_notebook(handle=target)

print("Live view done")
```

# IV Live Display

```python hide_input=true
def update_display():
    global I,V,t
    ipython.kernel.do_one_iteration()
    I,V,t = get_IV_t()
    ipython.kernel.do_one_iteration()
    source1.data = dict(x=I, y=V)
    
    # Update widgets
    I_avg = np.average(I)
    I_ptp = np.max(I) - np.min(I)
    V_avg = np.average(V)
    V_ptp = np.max(V) - np.min(V)
    widgets = (I_avg_w, V_avg_w,  R_avg_w, I_ptp_w, V_ptp_w, R_ptp_w)
    vals = (I_avg, V_avg, V_avg / I_avg, I_ptp, V_ptp, V_ptp / I_ptp)
    for w,v in zip(widgets,vals):
        w.value = "%.2f" % v
    
    # Update scaling
    if scale_mode.value == "Autoscale":
        auto_scale(p1.x_range,I)
        auto_scale(p1.y_range,V)
    elif scale_mode.value == "Full range":
        p1.x_range.start = -x_full_scale
        p1.x_range.end = x_full_scale
        p1.y_range.start = -y_full_scale
        p1.y_range.end = y_full_scale
        
    x = np.array([np.min(I), np.max(I)])
    R = line_slope.value
    y = np.array([V_avg - I_ptp*R/2, V_avg+I_ptp*R/2])
    source2.data = dict(x=x, y=y)
    
    if plot_line_button.value:
        slope_line.visible = True
    else:
        slope_line.visible = False    
        
        
range1 = Range1d()
range2 = Range1d()
source1 = ColumnDataSource()
source2 = ColumnDataSource()

p1 = figure(title="IV", plot_height=500, plot_width=600,toolbar_location=None, x_range = range1, y_range = range2)
p1.yaxis.axis_label = 'Measured Voltage (uV)'
p1.xaxis.axis_label = 'Applied Current (uA)'
p1.line('x', 'y', source=source1)
slope_line = p1.line('x','y',source=source2, color='red')
slope_line.visible = False

x_full_scale = I_full_scale #uA
y_full_scale = V_full_scale # uV

stop_button = widgets.ToggleButton(description='Stop')
pause_button = widgets.ToggleButton(description='Pause')

scale_mode = widgets.RadioButtons(options=['Autoscale', 'Full range', 'Freeze range'])
mrsquid_mode = widgets.RadioButtons(options=['V-I', 'V-Phi'])

style = {'description_width': 'initial'}

I_avg_w = widgets.Text(description="I_avg (uA)",style=style)
V_avg_w = widgets.Text(description="V_avg (uV)",style=style)
I_ptp_w = widgets.Text(description="I_ptp (uA)",style=style)
V_ptp_w = widgets.Text(description="V_ptp (uV)",style=style)
R_ptp_w = widgets.Text(description="V/I ptp (ohms)",style=style)
R_avg_w  = widgets.Text(description="V/I avg (ohms)",style=style)

plot_line_button = widgets.ToggleButton(description="Plot line")
line_slope = widgets.FloatText(description="R_line (ohms)")
line_slope.value = 50
fit_slope_button = widgets.Button(description="Fit")

def fit_slope(w):
    coef = np.polyfit(I,V,1)
    line_slope.value = coef[0]
fit_slope_button.on_click(fit_slope)

for w in (I_avg_w, V_avg_w, I_ptp_w, V_ptp_w, R_ptp_w, R_avg_w):
    w.layout.width = '150px'
    
def save_data(w):
    fmt = f"mrsquid_{mrsquid_mode.value}_%Y-%m-%d-%H_%M_%S.dat"
    outname = folder + "/" + datetime.now().strftime(fmt)
    np.savetxt(outname, np.array([t,I,V]).T)
    filename.value = outname
    
save_button = widgets.Button(description='Save data')
save_button.on_click(save_data)
filename = widgets.HTML(description="Last filename: ",style=style)
filename.value = "(none)"

update_display()
target = show(p1, notebook_handle=True)
rows = []

rows.append(widgets.HBox([I_avg_w, V_avg_w, R_avg_w]))
rows.append(widgets.HBox([I_ptp_w, V_ptp_w, R_ptp_w]))
rows.append(widgets.HBox([plot_line_button, line_slope, fit_slope_button]))
rows.append(widgets.HBox([save_button, filename]))
rows.append(widgets.HBox([stop_button,pause_button,scale_mode,mrsquid_mode]))
controls = widgets.VBox(rows)
display(controls)

while True:
    ipython.kernel.do_one_iteration()
    if not pause_button.value:
        update_display()
    if stop_button.value:
        break
    push_notebook(handle=target)

print("Live view stopped")
```

# Temperature trace recorder

```python hide_input=true
# Global variables 
R = []
T = []
tm = []
Iavg = []
Iptp = []
Vavg = []
Vptp = []
t0 = time.time()
timer_t0 = time.time()

# The functions

def update_iv():
    I,V,_ = get_IV_t()
    source1.data = dict(x=I, y=V)

def get_temperature(serial_port, baudrate, timeout):
    # serial_port='COM7', baudrate=9600, timeout=1
    
    #variables (could be adjusted to inputs?)
    alpha = 4.1139e-3; R0 = 1000; R1=1081; Vin=3.272        #should be according to literature: alpha = 3.85e-3 & R0 = 1000
    
    #open Serial port to Arduino + flush what's still on it
    ser = serial.Serial(serial_port,baudrate,timeout=0.5)
    ser.flush()
    
    #sent command to Arduino so we can read it out:
    time.sleep(0.5)                       #needs to sleep, otherwise we get a: 'can not convert float... at the Vout-line'
    ser.reset_input_buffer()
    ser.write(b'g')
    
    #Serial read + Temperature calculations:
    Vout = float(str(ser.readline().decode("utf-8")))*(Vin/1023.0)          #calculates the voltage from the voltage divider (over the PT1000)
    Rpt = (Vout*R1)/(Vin-Vout)                                              #calculates the resistance of the PT1000
    Temp_C = ((Rpt/R0)-1)/alpha                                             #with the above calculated resistance and 'known' variable we calculate the temperature in Celsius
    Temp_K = Temp_C + 273.15                                                #Temperature in Kelvin
    
    return Temp_K


def take_measurement():
    global timer_t0 
    timer_t0 = time.time()
    T.append(get_temperature('COM7',9600,1)) 
    I,V,_ = get_IV_t()
    R.append(np.polyfit(I,V,1)[0])
    Iavg.append(np.average(I))
    Iptp.append(np.max(I)-np.min(I))
    Vavg.append(np.average(V))
    Vptp.append(np.max(V)-np.min(V))    
    tm.append(time.time()-t0)
    source2.data = dict(x=tm,y=R)
    source3.data = dict(x=tm,y=T)
    source4.data = dict(x=T,y=R)

# The plots

p1 = figure(title="IV", plot_height=300, plot_width=300,toolbar_location=None)
p1.yaxis.axis_label = 'Measured Voltage (uV)'
p1.xaxis.axis_label = 'Applied Current (uA)'
p2 = figure(title="R vs time", plot_height=300, plot_width=300,toolbar_location=None)
p2.yaxis.axis_label = 'Resistance (Ohm)'
p2.xaxis.axis_label = 'Time (s)'
p3 = figure(title="T vs time", plot_height=300, plot_width=300,toolbar_location=None)
p3.yaxis.axis_label = 'Temperature (K)'
p3.xaxis.axis_label = 'Time (s)'
p4 = figure(title="R vs T", plot_height=300, plot_width=900,toolbar_location=None)
p4.yaxis.axis_label = 'Resistance (Ohm)'
p4.xaxis.axis_label = 'Temperature (K)'

init_data = dict(x=[0,1],y=[0,1])

source1 = ColumnDataSource()
source1.data = init_data
p1.line('x', 'y', source=source1)
source2 = ColumnDataSource()
source2.data = init_data
p2.line('x', 'y', source=source2)
p2.circle('x', 'y', source=source2, fill_color="white", size=8)
source3 = ColumnDataSource()
source3.data = init_data
p3.line('x', 'y', source=source3)
p3.circle('x', 'y', source=source3, fill_color="white", size=8)
source4 = ColumnDataSource()
source4.data = init_data
p4.circle('x', 'y', source=source4, fill_color="white", size=8)

take_measurement()

target = show(column(row(p1,p2,p3),p4), notebook_handle=True)

# The widgets

save_button = widgets.Button(description='Save data')
def save_data(w):
    fmt = f"mrsquid_R_vs_T_%Y-%m-%d-%H_%M_%S.dat"
    outname = folder + "/" + datetime.now().strftime(fmt)
    np.savetxt(outname, np.array([R,T,tm,Iavg,Iptp,Vavg,Vptp]).T)
    filename.value = outname
save_button.on_click(save_data)

filename = widgets.HTML(description="Last filename: ",style=style)
filename.value = "(none)"

reset_button = widgets.Button(description='Restart data')
def reset(w):
    global R,T,t,Iavg,Iptp,Vavg,Vptp
    R = []
    T = []
    tm = []
    Iavg = []
    Iptp = []
    Vavg = []
    Vptp = []
    t0 = time.time()
    take_measurement()
reset_button.on_click(reset)

meas_button = widgets.Button(description='Measure now')
def meas(w):
    take_measurement()
meas_button.on_click(meas)

stop_button = widgets.ToggleButton(description='Stop')
pause_button = widgets.ToggleButton(description='Pause')
interval = widgets.FloatText(description="Measurement interval (s)",style={'description_width': 'initial'})
interval.value = 5

rows = []
rows.append(widgets.HBox([interval,meas_button]))
rows.append(widgets.HBox([pause_button, save_button, reset_button, stop_button]))
rows.append(widgets.HBox([filename]))
controls = widgets.VBox(rows)
display(controls)

while True:
    ipython.kernel.do_one_iteration()
    update_iv()
    if (time.time()-timer_t0) > interval.value and not pause_button.value:
        take_measurement()
    if stop_button.value:
        break
    push_notebook(handle=target)
    
print("Live view stopped")
```

# Power Spectral Density Live Display

```python hide_input=true
# The traces we download
I = []
V = []
t = []

# The calculated spectral densities
I_psd = []
V_psd = []
f = []

# Each PSD trace is 8193 points long
n_psd = 8193
nmax  =  100
I_psd_array = np.zeros([nmax,n_psd])
V_psd_array = np.zeros([nmax,n_psd])

# For keeping track of how many we've filled and where the next goes
ind_next = 0
n_filled = 0

def update_display():
    global I,V,t,f,I_psd,V_psd,last_trace_time,ind_next,n_filled
    I,V,t = get_IV_t(decimation=decimation.value)
    I_avg = np.average(I)
    I_ptp = np.max(I) - np.min(I)
    p1.title.text = "Current          Average = %.2f uA     Peak-to-Peak %.2f uA" % (I_avg, I_ptp)
    V_avg = np.average(V)
    V_ptp = np.max(V) - np.min(V)
    p2.title.text = "Voltage          Average = %.2f uV     Peak-to-Peak %.2f uV" % (V_avg, V_ptp)
    f, I_psd_trace = welch(I,1/t[1],nperseg=len(I))
    f, V_psd_trace = welch(V,1/t[1],nperseg=len(V))
    
    if enable_average.value:
        # ind_next: the running index
        # n_filled: the number of traces taken since the last reset
        
        I_psd_array[ind_next,:] = I_psd_trace
        V_psd_array[ind_next,:] = V_psd_trace
        
        if n_filled  == 0:
            I_psd = I_psd_trace
            V_psd = V_psd_trace
        elif n_filled < nmax:
            I_psd = I_psd * n_filled / (n_filled+1) + I_psd_trace
            V_psd = V_psd * n_filled / (n_filled+1) + V_psd_trace
        else:
            # this should work because of numpy array index wrapping
            I_psd += I_psd_trace / nmax - I_psd_array[ind_next-nmax]/nmax
            V_psd += V_psd_trace / nmax - V_psd_array[ind_next-nmax]/nmax  
        ind_next += 1
        if ind_next == 100:
            ind_next = 0
        if n_filled < 100:
            n_filled += 1
            averaged_num.value = "%d" % n_filled
            
    else:
        I_psd = I_psd_trace
        V_psd = V_psd_trace
        n_filled = 0
        ind_next = 0
        if averaged_num.value != "1":
            averaged_num.value = "1"
    
    source1.data = dict(x=f, y=I_psd)
    source2.data = dict(x=f, y=V_psd)
    if manual_freq.value:
        p1.x_range.start = freq_start.value
        p1.x_range.end = freq_stop.value
        p2.x_range.start = freq_start.value
        p2.x_range.end = freq_stop.value
    else:
        p1.x_range.start = f[0]
        p1.x_range.end = f[-1]
        p2.x_range.start = f[0]
        p2.x_range.end = f[-1]
        
range1 = Range1d()
range2 = Range1d()
source1 = ColumnDataSource()
source2 = ColumnDataSource()

p1 = figure(title="Current", plot_height=300, y_axis_type="log",
            plot_width=900, toolbar_location=None, x_range = range1)
p1.yaxis.axis_label = 'Curret PSD dB(A^2/Hz)'

p2 = figure(title="Voltage", plot_height=300, y_axis_type="log",
            plot_width=900, toolbar_location=None, x_range = range2)
p2.xaxis.axis_label = 'Frequency (Hz)'
p2.yaxis.axis_label = 'Voltage PSD dB(V^2/Hz)'

p1.line('x', 'y', source=source1)
p2.line('x', 'y', source=source2)

stop_button = widgets.ToggleButton(description='Stop')
pause_button = widgets.ToggleButton(description='Pause')

def reset_average(w):
    global n_filled, ind_next
    n_filled = 0
    ind_next = 0

# For the decimation GUI
options = []
for d,t in zip(s.decimation_options[9:], s.duration_options[9:]):
    if t < 1e-3:
        t*= 1e6
        u = "us"
    elif t<1:
        t*= 1e3
        u = "ms"
    else: 
        u = "s"
    options.append((str(d) + ", " "{0:.3g}".format(t) + " " + u,d))
options

def update_freq(w):
    s.decimation = decimation.value
    trace_info.value = "Frequency resolution: %.2f Hz &nbsp; &nbsp; &nbsp; &nbsp;" % (1/s.duration)
    trace_info.value += "Max frequency: %.2f Hz" % (1/s.duration*n_psd)

decimation = widgets.Dropdown(description="Decimation:", 
                              options=options)
decimation.observe(reset_average)
decimation.observe(update_freq)
trace_info = widgets.HTML()

enable_average = widgets.Checkbox(description="Enable averaging (max 100 traces)", style=style)
averaged_num = widgets.HTML(value="1", description="Number averaged:", style=style)
reset_button = widgets.Button(description="Reset averaging")
reset_button.on_click(reset_average)

manual_freq = widgets.Checkbox(description="Manual frequency range")
freq_start = widgets.FloatText(value=0,description="Start (Hz)")
freq_stop = widgets.FloatText(value=1e5,description="Stop (Hz)")

def save_data(w):
    fmt = f"mrsquid_PSD_%Y-%m-%d-%H_%M_%S.dat"
    outname = folder + "/" + datetime.now().strftime(fmt)
    np.savetxt(outname, np.array([f,I_psd,V_psd]).T)
    filename.value = outname
    
save_button = widgets.Button(description='Save data')
save_button.on_click(save_data)
filename = widgets.HTML(description="Last filename: ",style=style)
filename.value = "(none)"

update_display()
update_freq(0)
target = show(column(p1,p2), notebook_handle=True)
rows = []
rows.append(widgets.HBox([decimation,trace_info]))
rows.append(widgets.HBox([enable_average, averaged_num, reset_button]))
rows.append(widgets.HBox([manual_freq,freq_start,freq_stop]))
rows.append(widgets.HBox([save_button, filename]))
rows.append(widgets.HBox([stop_button,pause_button]))
display(widgets.VBox(rows))

while True:
    ipython.kernel.do_one_iteration()
    if not pause_button.value:
        update_display()
    ipython.kernel.do_one_iteration()
    if stop_button.value:
        break
    push_notebook(handle=target)
    time.sleep(0.1)

print("Live view done")
```