""" File to evaluate the Shapiro curves we have measured to get the step width depending on the modulation amplitude in Fig 2 B. The same measurements are used as in Fig 2 A """ # 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 import numpy as np import matplotlib.pyplot as pl from scipy.signal import find_peaks from scipy.special import jv from helper import * from style import * meas_list = [ { "path": "Exp_Shapiro_steps_f_91_Hz_I_m_0.0.csv", "mod_amp": 0., # in µm "color": "grey", "mod_freq": 0, 'v_max': 0.8 }, { "path":"Exp_Shapiro_steps_f_91_Hz_I_m_0.3.csv", "mod_amp": 0.2, # in µm "mod_freq": 90, "p0": [90, 50, 100, 100, 0.05, 0.25, 0.55, 5], "color": "orange", "steps": 3, 'v_max': 0.8 }, { "path": "Exp_Shapiro_steps_f_91_Hz_I_m_0.4.csv", "mod_amp": 0.3, "mod_freq": 90, "p0": [90, 50, 100, 100, 0.05, 0.25, 0.55, 5], "color": "steelblue", "steps": 3, 'v_max': 0.8 }, { "path":"Exp_Shapiro_steps_f_91_Hz_I_m_0.7.csv", "mod_amp": 0.5, "mod_freq": 90, "p0": [90, 40, 45, 40, 0.3, 0.2, 0.7, 5], "color": "darkgreen", "steps": 3, 'v_max': 0.8 }, { "path":"Exp_Shapiro_steps_f_91_Hz_I_m_0.8.csv", "mod_amp": 0.6, "mod_freq": 90, "p0": [90, 50, 100, 100, 0.05, 0.25, 0.55, 5], "color": "purple", "steps": 3, 'v_max': 0.8 }, { "path":"Exp_Shapiro_steps_f_91_Hz_I_m_1.0.csv", "mod_amp": 0.7, "mod_freq": 90, "p0": [65, 100, 100, 100, 0.108, 0.25, 0.55, 5], "color":"firebrick", "steps":2, 'v_max': 0.8 }, { "path": "Exp_Shapiro_steps_f_91_Hz_I_m_1.2.csv", "mod_amp": 0.9, "mod_freq": 90, "p0": [90, 50, 100, 100, 0.05, 0.25, 0.55, 5], "color": "royalblue", "steps": 3, 'v_max': 0.9 }, { "path":"Exp_Shapiro_steps_f_91_Hz_I_m_1.4.csv", "mod_amp": 1, "mod_freq": 90, "p0": [90, 50, 100, 100, 0.05, 0.25, 0.55, 5], "color": "navy", "steps": 3, 'v_max': 0.8 }, { "path": "Exp_Shapiro_steps_f_91_Hz_I_m_1.8.csv", "mod_amp": 1.3, "mod_freq": 90, "p0": [90, 50, 100, 100, 0.05, 0.25, 0.55, 5], "color": "steelblue", "steps": 3, 'v_max': 1 }, { "path": "Exp_Shapiro_steps_f_91_Hz_I_m_1.9.csv", "mod_amp": 1.4, "mod_freq": 90, "p0": [90, 50, 100, 100, 0.05, 0.25, 0.55, 5], "color": "darkgreen", "steps": 3, 'v_max': 1 }, { "path":"Exp_Shapiro_steps_f_91_Hz_I_m_1.5.csv", "mod_amp": 1.1, "mod_freq": 90, "p0": [90, 50, 100, 100, 0.05, 0.25, 0.55, 5], "color": "orange", "steps": 3, 'v_max': 0.8 }, { "path":"Exp_Shapiro_steps_f_91_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.8 }, { "path": "Exp_Shapiro_steps_f_91_Hz_I_m_1.0.csv", "mod_amp": 0.7, "mod_freq": 90, "p0": [90, 50, 100, 100, 0.05, 0.25, 0.55, 5], "color": "orange", "steps": 3, 'v_max': 1 }, ] # data path path_data_repo = "../Fig_2/shapiro_curves" # parameters set for evaluation, see paper data_supplementary material epsilon = 45 delta_f = 90 #fig1, (ax1, ax2) = pl.subplots(2, 1,figsize=(8,6)) # init result array result = [] for meas in meas_list: # get data vel, I, mu, mu_err, I_m, f, mod_amp, dz, dz_err = get_meas_data(meas["path"], path_data_repo=path_data_repo, mod_amp_return = True) style = get_style(meas["color"]) v_max = 1 #meas['v_max'] #ax1.errorbar(vel[vel <= v_max],mu[vel <= v_max],mu_err[vel <= v_max], **style, label=f'$A_m = {meas["mod_amp"]:.1f}\mathrm{{µm}}$, f= {meas["mod_freq"]} Hz') # plot each shapiro curve in a single plot fig3, (ax3, ax4) = pl.subplots(2, 1, figsize=(4, 6)) ax3.errorbar(vel[vel <= v_max],mu[vel <= v_max],mu_err[vel <= v_max], **style, label=f'$A_m = {meas["mod_amp"]:.1f}\mathrm{{µm}}$, f= {meas["mod_freq"]} Hz') ax3.set_xlabel("Velocity [mm/s]") ax3.set_ylabel("$\Delta \mu$ ") ax3.legend() # fig3.tight_layout() # fig3.savefig(join(f"single_trajectory_{mod_amp}.png")) # calculate the derivative # calculate the derivative x = vel # Estimate the smoothing and error on the number of data points. sigma = 5 # Smooth the data and calculate the numerical derivative. smooth = gaussian_filter1d(dz, sigma=sigma) derivative = np.gradient(smooth) / np.gradient(x) # Search for maxima in the derivative, i.e., flanks in the actual data. peaks, properties = find_peaks(derivative, distance=5, height=10e-7) # plot the derivative ax4.plot(x[vel <= v_max], derivative[vel <= v_max], ls=' ', marker='d', color=meas["color"], label=f'$A_m = {meas["mod_amp"]:.1f}\mathrm{{µm}}$, f= {meas["mod_freq"]} Hz') ax4.plot(x[peaks], derivative[peaks], marker='X', ls='none', color='red') # Error estimation v_err = np.mean(np.diff(vel)) * sigma ############ step 0 ############### step_number = 0 mask = (mu <= (step_number * delta_f + epsilon)) & (mu >= (step_number * delta_f - epsilon)) ax3.axhline(step_number * delta_f, ls='--', color='grey') # first plateau try: # search maximum of the derivative to get the position of first step v_0 = vel[mask][np.argmax(derivative[mask])] ax4.plot(vel[mask][np.argmax(derivative[mask])], derivative[np.argmax(derivative[mask])], marker='X', ls='none', color='blue') ax3.axvline(x[mask][np.argmax(derivative[mask])], ls='--', color='red') ax4.axvline(x[mask][np.argmax(derivative[mask])], ls='--', color='red') except IndexError: v_0 = 0 ############ step 1 ############### step_number = 1 mask = (mu <= (step_number * delta_f + epsilon)) & (mu >= (step_number * delta_f - epsilon)) # print(mask) ax3.axhline(step_number * delta_f, ls='--', color='grey') ax3.axhline(step_number * delta_f + epsilon, ls='--', color='green') ax3.axhline(step_number * delta_f - epsilon, ls='--', color='green') try: # find the two maxima of the derivative within the epsilon region, # giving the beginning and the end of the shapiro plateau length = int(len(x[mask]) / 2) left_side_derivative = derivative[mask][:length] left_side_x = x[mask][:length] right_side_derivative = derivative[mask][length + 1:] right_side_x = x[mask][length+1:] v_1 = np.abs(right_side_x[np.argmax(right_side_derivative)] - left_side_x[np.argmax(left_side_derivative)]) ax4.plot(right_side_x[np.argmax(right_side_derivative)], right_side_derivative[np.argmax(right_side_derivative)], marker='X', ls='none', color='blue') ax4.plot(left_side_x[np.argmax(left_side_derivative)], left_side_derivative[np.argmax(left_side_derivative)], marker='X', ls='none', color='blue') ax3.axvline(left_side_x[np.argmax(left_side_derivative)], ls='--', color='green') ax4.axvline(left_side_x[np.argmax(left_side_derivative)], ls='--', color='green') ax3.axvline(right_side_x[np.argmax(right_side_derivative)], ls='--', color='green') ax4.axvline(right_side_x[np.argmax(right_side_derivative)], ls='--', color='green') ax4.plot(vel[mask], derivative[mask], color='red') except IndexError: # v_1 = np.nan v_1 = 0 ############ step 2 ############### step_number = 2 mask = (mu <= (step_number * delta_f + epsilon)) & (mu >= (step_number * delta_f - epsilon)) # print(mask) ax3.axhline(step_number * delta_f, ls='--', color='grey') ax3.axhline(step_number * delta_f + epsilon, ls='--', color='blue') ax3.axhline(step_number * delta_f - epsilon, ls='--', color='blue') # widths = widths[1:] - widths[:-1] try: # search also the position of the second step even though i tis not used in the paper length = int(len(x[mask]) / 2) left_side_derivative = derivative[mask][:length] left_side_x = x[mask][:length] right_side_derivative = derivative[mask][length + 1:] right_side_x = x[mask][length+1:] v_2 = np.abs(right_side_x[np.argmax(right_side_derivative)] - left_side_x[np.argmax(left_side_derivative)]) ax4.plot(right_side_x[np.argmax(right_side_derivative)] , right_side_derivative[np.argmax(right_side_derivative)], marker='X', ls='none', color='blue') ax4.plot(left_side_x[np.argmax(left_side_derivative)] , left_side_derivative[np.argmax(left_side_derivative)], marker='X', ls='none', color='blue') ax3.axvline(left_side_x[np.argmax(left_side_derivative)], ls='--', color='blue') ax4.axvline(left_side_x[np.argmax(left_side_derivative)], ls='--', color='blue') ax3.axvline(right_side_x[np.argmax(right_side_derivative)], ls='--', color='blue') ax4.axvline(right_side_x[np.argmax(right_side_derivative)], ls='--', color='blue') ax4.plot(vel[mask], derivative[mask], color='red') except IndexError: # v_1 = np.nan v_2 = 0 if meas["mod_amp"] == 0: # set the values for I_M = 0, where I0 is basically the critical velocity, set the step width of step 1 and 2 # to zero since there is no step. I0 = v_0 v_1 = 0 v_2 = 0 # put the values in a list result.append([mod_amp,v_0,v_1,v_2,v_err]) ax4.set_xlabel("Velocity [mm/s]") ax4.set_ylabel("Derivative $dz/dv$") ax4.legend() fig3.tight_layout() #fig3.savefig(join( f'single_trajectory_{mod_amp}.png')) # pl.close() res = np.asarray(result) mod_amp = res[:,0] # normalize the step widths by I0 /(I_c) v_0 = res[:,1] / I0 v_1 = res[:,2] / I0 v_2 = res[:,3] / I0 # sort the modulation amplitudes sortArg = np.argsort(mod_amp) mod_amp = mod_amp[sortArg] v_0 = v_0[sortArg] v_1 = v_1[sortArg] v_2 = v_2[sortArg] v_err = res[:, 4] / I0 I_err = v_err ################################################# ################################################# # Fit the results n_max = 3 shape = (n_max, len(meas_list)) 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() fig, (ax1, ax2) = pl.subplots(2, 1)#,figsize=(5.5,4.5)) I_c = current(1,v_c=1) # get the critical current from the func keyword parameter R = 0.9e-3 # resistance prefac = I_c * R / f #2.27 see paper def bes_0(x): return np.abs(jv(0, prefac * x)) def bes_1(x): return np.abs(jv(1, prefac * x)) def bes_2(x): return np.abs(jv(2, prefac * x)) x_values = np.tile(mod_amp, n_max) y_values = np.concatenate((v_0,v_1,v_2)) # 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) # save results in a Dataframe data = {'mod_amp':mod_amp, 'I_0':v_0,'I_0_err':I_err, 'I_1':v_1,'I_1_err':v_err} DATA = pd.DataFrame(data) ################################################ ############################################### # plot the results # modulation amplitudes in current I(v_m) ; v_m = omega_m * x_m I_m = current(mod_amp * 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 # plot the fit vaues ax1.plot(I_fit, fit_values[0, :],color = 'darkgreen', ls='-', marker='', label=f'Fit $n = {0}$') ax2.plot(I_fit, fit_values[1, :],color = 'navy', ls='-', marker='', label=f'Fit $n = {1}$') # plot theoretical considerations ax1.plot(I_fit, bes_0(I_fit),color = 'orange', ls='-', marker='', label=r'$\mathrm{J}_{0}(\alpha_{\mathrm{th}} \cdot I_\mathrm{m}/I_c)$') ax2.plot(I_fit, bes_1(I_fit),color = 'orange', ls='-', marker='', label=r'$\mathrm{J}_{1}(\alpha_{\mathrm{th}} \cdot I_\mathrm{m}/I_c)$') # plot the data ax1.errorbar(I_m,v_0,label=r'$I_0$', yerr=I_err,**get_style('darkgreen',ls = '--')) ax2.errorbar(I_m,v_1,label=r'$I_1$', yerr=I_err,**get_style('navy',ls = '--')) ax2.set_xlabel("Modulation Amplitude $I_\mathrm{m}/I_c$") ax1.set_ylabel(r"$I_0/I_c$") ax2.set_ylabel(r"$I_1/I_c$") ax1.legend() ax2.legend(loc='lower right') fig.tight_layout() pl.show()