# file to evaluate a measurement import os.path, glob from os.path import join import json import pandas as pd from scipy.ndimage import gaussian_filter1d from scipy.optimize import curve_fit from scipy.signal import find_peaks from scipy.special import jv import scipy.constants as c import numpy as np import matplotlib.pyplot as pl from helper import * from style import * # path to data path_data_repo = "../Fig_2" ### Import the data ### ###### step height ###### # Experiment tab_step_height = pd.read_csv(os.path.join(path_data_repo,"Exp_step_height_data.csv")) ## single imbalance curves for plotting meas_list = [ { "path": "Exp_Shapiro_steps_f_135f_Hz_I_m_1.4.csv", "mod_amp": 0.7, "mod_freq": 135, "p0": [135, 100, 100, 100, 0.108, 0.25, 0.55, 5], "color": "darkgreen", "steps": 3, 'v_max': 1, 'v_plot': 0.4, 'plot': True }, { "path": "Exp_Shapiro_steps_f_90f_Hz_I_m_1.6.csv", "mod_amp": 1.2, "mod_freq": 90, "p0": [90, 50, 100, 100, 0.05, 0.25, 0.55, 5], "color": "orange", "steps": 3, 'v_max': 0.65, 'v_plot': 0.35, 'plot': True }, { "path": "Exp_Shapiro_steps_f_60f_Hz_I_m_1.1.csv", "mod_amp": 1.2, "mod_freq": 60, "color": "navy", "p0": [60, 50, 50, 50, 0.15, 0.2, 0.2, -10], "steps": 3, 'v_max': 0.65, 'v_plot': 0.32, 'plot': True }, ] # Simulations tab_step_height_simulation = pd.read_csv(os.path.join(path_data_repo,"Sim_step_height_data.csv")) ## Step Width ##### # Experiment tab_step_width = pd.read_csv(os.path.join(path_data_repo,"Exp_step_width_data.csv")) # Simulations tab_step_width_sim = pd.read_csv(os.path.join(path_data_repo,"Sim_step_width_data.csv")) ##### Plot the data ### fig, ax = pl.subplot_mosaic( [["b", "a"], ["c", "d"]], layout="constrained", figsize=(7, 7), dpi=300, gridspec_kw={'height_ratios': [0.8, 1]}) set_plot_style() # Step height experiment # get data mod_freq = tab_step_height["mod_freq"] step_height = tab_step_height["step_height"] step_height_err = tab_step_height["step_height_error"] # plot ax["a"].errorbar(mod_freq, step_height, yerr=step_height_err, **get_style(RPTU_COLORS["nacht"]), label='Exp.', markersize=6, markeredgewidth=1.5,zorder = 3) # µ = f curve x = np.arange(0, 200, 0.01) ax["a"].plot(x, x, color=RPTU_COLORS['mango'], label=fr'$\mu / h = f_\mathrm{{m}}$', lw=2.5) # Simulations # get data freq = tab_step_height_simulation["freq"].to_numpy() d_mu = tab_step_height_simulation["d_mu"].to_numpy() d_mu_err = tab_step_height_simulation["d_mu_err"].to_numpy() # plot data ax["a"].errorbar(freq, d_mu, yerr=d_mu_err, **get_style(RPTU_COLORS["himbeere"], marker='v'), label='Sim.', markersize=5.5, markeredgewidth=1.5) ax["a"].legend(fontsize=10) ax["a"].set_xlabel("$f_\mathrm{m}$ [Hz]") ax["a"].set_ylabel(r"$\Delta \mu \, / \, h$ [Hz]") # Plot the Shapiro example curves in A n = 3 colors = [RPTU_COLORS['himbeere'],RPTU_COLORS['nacht'],RPTU_COLORS['mango']] for i, meas in enumerate(meas_list): # get meas data vel, I, mu, mu_err, I_m, f = get_meas_data(meas['path'],path_data_repo=path_data_repo) I_c = current(1, v_c=1)# get the critical current from the func keyword parameter # fit the data : try: if meas["steps"] == 2: sigmoid = sigmoid_2 elif meas["steps"] == 1: sigmoid = sigmoid_1 else: sigmoid = sigmoid_3 v_max = meas['v_max'] popt, pcov = curve_fit(sigmoid, vel[vel <= v_max], mu[vel <= v_max], p0=meas["p0"], bounds=(-10, np.inf)) v = np.arange(0, vel[-1], 0.001) p_err = np.sqrt(np.diag(pcov)) I = current(vel) / I_c I_fit = current(v) / I_c color = colors[i] style = get_style((cl.to_rgba(color, 0.9))) ax["b"].errorbar(I[vel <= v_max], mu[vel <= v_max], mu_err[vel <= v_max], **style, label=fr'$f_\mathrm{{m}}= {meas["mod_freq"]}$ Hz',) ax["b"].plot(I_fit, sigmoid(v, *popt), color=color,lw=2) except RuntimeError: pass ax["b"].legend(fontsize=10) ax["b"].set_xlim(-0.02, 1.5) ax["b"].set_ylim(-10, 270) ax["b"].set_xlabel(r"Current $I/I_\mathrm{c}$ ") ax["b"].set_ylabel(r"$\Delta \mu \, / \, h$ [Hz]") ###### Step Width ##### # Experiment mod_amp = tab_step_width['mod_amp'].to_numpy() # in µm I_0 = tab_step_width['I_0'].to_numpy() # in I_0/I_c I_1 = tab_step_width['I_1'].to_numpy() # in I_1/I_c I_err = tab_step_width['I_err'].to_numpy() # in I_err/I_c f = 90 # modulation amplitudes in current I(v_m) ; v_m = omega_m * x_m in mm/s I_m = current(mod_amp * 2 * np.pi * f * 1e-3) / I_c color_n0 = RPTU_COLORS["nacht"] color_n1 = RPTU_COLORS["nacht"] color_Jth = RPTU_COLORS["mango"] # plot the data ax["c"].errorbar(I_m, I_0, label=r'$I_0$', yerr=I_err, **get_style(color_n0, ls=''), markersize=6, markeredgewidth=1.5, zorder = 3) ax["d"].errorbar(I_m, I_1, label=r'$I_1$', yerr=I_err, **get_style(color_n1, ls=''), markersize=6, markeredgewidth=1.5, zorder = 3) ### Fit the data ##### n_max = 2 shape = (n_max, len(mod_amp)) def bessel(x, alpha, shape=shape): # Sum of Bessel functions for the different plateau numbers. x_values = np.reshape(x, shape) result = np.zeros(shape) for i, A in enumerate([1, ] * n_max): result[i, :] = A * np.abs(jv(i, alpha * x_values[i])) ** 1 return result.flatten() x_values = np.tile(mod_amp, n_max) y_values = np.concatenate((I_0, I_1)) # The data needs to be transposed to match the x-values. p_opt, p_cov = curve_fit(bessel, x_values, y_values, p0=[0.1]) p_err = np.sqrt(np.diag(p_cov)) # Get the fitted values, which are already flattened and reshape them. mod_amp_fit = np.arange(0, mod_amp[-1], 0.01) x_values_fit = np.tile(mod_amp_fit, n_max) shape_fit = (n_max, len(mod_amp_fit)) fit_values = bessel(x_values_fit, *p_opt, shape=shape_fit).reshape(shape_fit) # modulation amplitudes in current for the fit v_m = omega_m * x_m I_fit = current(mod_amp_fit * 2 * np.pi * f * 1e-3) / I_c # plot the fit results ax["c"].plot(I_fit, fit_values[0, :], color=color_n0, ls='-', marker='', lw=2, label=f'Fit $n = {0}$') ax["d"].plot(I_fit, fit_values[1, :], color=color_n1, ls='-', marker='', lw=2, label=f'Fit $n = {1}$') ##### Theory values ##### I_c = current(1, v_c=1) # get the critical current from the func keyword parameter R = 0.907e-3 # resistance see paper prefac = I_c * R / f #2.27 print(prefac) # plot the theory estimate def bes_0(x): return np.abs(jv(0, prefac * x ** 1)) def bes_1(x): return np.abs(jv(1, prefac * x ** 1)) ax["c"].plot(I_fit, bes_0(I_fit), color=color_Jth, ls='--', marker='', lw=2, label=r'$\mathrm{J}_{0}(\alpha_{\mathrm{th}} \cdot I_\mathrm{m}/I_c)$') ax["d"].plot(I_fit, bes_1(I_fit), color=color_Jth, ls='--', marker='', lw=2, label=r'$\mathrm{J}_{1}(\alpha_{\mathrm{th}} \cdot I_\mathrm{m}/I_c)$') #### Simulations #### # open simulation data mod_amp_sim = tab_step_width_sim["mod_amp"].to_numpy() I0_sim = tab_step_width_sim["I0"].to_numpy() I0_sim_err = tab_step_width_sim["I0_err"].to_numpy() I1_sim = tab_step_width_sim["I1"].to_numpy() - I0_sim I1_sim_err = tab_step_width_sim["I1_err"].to_numpy() mod_amp_fit = np.arange(0, mod_amp[-1], 0.01) # modulation amplitudes in current I(v_m) ; v_m = omega_m * x_m I_m_sim = current(mod_amp_sim * 2 * np.pi * f * 1e-3) / I_c # modulation amplitudes in current for the fit v_m = omega_m * x_m I_fit = current(mod_amp_fit * 2 * np.pi * f * 1e-3) / I_c ax["c"].errorbar(I_m_sim, I0_sim, yerr=I0_sim_err, **get_style(RPTU_COLORS["himbeere"], marker='v', ls=''), label="Sim.", markersize=6, markeredgewidth=1.5) ax["d"].errorbar(I_m_sim, I1_sim, yerr=I1_sim_err, **get_style(RPTU_COLORS["himbeere"], marker='v', ls=''), label="Sim.", markersize=6, markeredgewidth=1.5) ax["c"].legend(fontsize=10) ax["d"].legend(fontsize=10) ax["c"].set_xlabel("Modulation Amplitude $I_\mathrm{m}/I_c$") ax["d"].set_xlabel("Modulation Amplitude $I_\mathrm{m}/I_c$") ax["c"].set_ylabel(r"$I_0 \, / \, I_c$") ax["d"].set_ylabel(r"$I_1 \, / \, I_c$") fig.tight_layout() #pl.savefig(("Fig2.png"), dpi=300) pl.show()