From 1cbab37d91493b59892f2ece7440f8065b6f8ddb Mon Sep 17 00:00:00 2001 From: Paolo Franchini Date: Tue, 7 Nov 2023 11:50:59 +0000 Subject: [PATCH] Config file and plotting scripts used for the paper. --- toys/config.py | 7 +- toys/lockin_toy_paper.py | 367 ++++++++++++++++++++++++++++++++++++ toys/plot_for_paper_shot.py | 8 +- 3 files changed, 377 insertions(+), 5 deletions(-) create mode 100755 toys/lockin_toy_paper.py diff --git a/toys/config.py b/toys/config.py index 6a1ec27..3fa7554 100644 --- a/toys/config.py +++ b/toys/config.py @@ -18,7 +18,7 @@ ## Cryostat: ################################################ pressure = 0 # [bar] pressure -ttc = 0.13 # T/Tc +ttc = 0.12 # T/Tc energy = 1000; # [eV] deposited energy @@ -51,4 +51,9 @@ M = 10.0e-9 # [H] v0 = 1e-3 # [m/sec] + +## Simulation limits: ####################################### + +diameter_max = 5*1000e-9 + ############################################################## diff --git a/toys/lockin_toy_paper.py b/toys/lockin_toy_paper.py new file mode 100755 index 0000000..96ac5b8 --- /dev/null +++ b/toys/lockin_toy_paper.py @@ -0,0 +1,367 @@ +''' +Toy MC simulation for the bolometer response +using the error from the LOCKIN circuit + +Input: + + - Energy + - Helium-3 Pressure: pressure + - Wire diameter: diameter + - Base temperature: temperature + + - Reponse time; t_w + - Decay constant: t_b + - Voltage height: v_h + - Voltage noise: v_rms + +Output: + - Preliminary plots: voltage error + - Error vs Pressure, Temperature, Diameter + +P. Franchini, 10/2022 + +''' + +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.mlab as mlab +from scipy.optimize import curve_fit +from scipy.stats import norm + +plt.rcParams.update({'font.size': 14}) +plt.rcParams.update({'lines.linewidth': 3}) +plt.rcParams.update({'xtick.direction': 'in'}) +plt.rcParams.update({'ytick.direction': 'in'}) + +# import Tsepelin code +exec(open("../mod_helium3.py").read()) + +# Configuration file with all the parameters +exec(open("config.py").read()) + +## Parameters ################################################ + +#volume = 1e-6 # [m^3] Helium-3 cell +#density = 6.05e3; # [kg/m^3] Niobium-Titanium (NbTi) + +#============================================================= + +#pressure = 0 # [bar] pressure +#ttc=0.1 # T/Tc +#diameter = 400e-9; # [m] vibrating wire diameter +#energy = 10; # [eV] deposited energy + +#============================================================= + +#t_b = 5.00 # [s] decay constant +#amp=100 # gain of the voltage cold amplifier +#v_h = amp*np.pi/2*1e-7 # [V] Base voltage height for a v=1mm/s +#v_rms = 7.9*1e-9 # [V] Error on voltage measurement for a lock-in amplifier + +#============================================================= + +N = 100 # number of toys +verbose=True # verbosity for plotting + +unused = 0.0 + +# =========================================================== + +## More routines: ########################################### + +def Width_from_Temperature(Temperature,PressureBar,Diameter): + + gap = energy_gap_in_low_temp_limit(PressureBar) + width=np.power(Fermi_momentum(PressureBar),2)*Fermi_velocity(PressureBar)*density_of_states(PressureBar)/(2*density*np.pi*Diameter)*np.exp(-gap/(Boltzmann_const*Temperature)) + return width + +def Temperature_from_Width(Width,PressureBar,Diameter): + + gap = energy_gap_in_low_temp_limit(PressureBar) + temperature=-gap/(Boltzmann_const*np.log( Width*2*density*np.pi*Diameter/(np.power(Fermi_momentum(PressureBar),2)*Fermi_velocity(PressureBar)*density_of_states(PressureBar)))) + return temperature + +def DeltaWidth_from_Energy(E,PressureBar,BaseTemperature,Diameter): + # Find delta width from the input energy deposition for a certain base temperature + + # find fit line for the Width variation vs Deposited energy for the base temperature + W0=Width_from_Temperature(BaseTemperature,PressureBar,Diameter) + + DQ = np.array([]) # delta energy [eV] + DW = np.array([]) # delta width [Hz] + + #for dw in np.arange(0,2.5,0.001): # Delta(Deltaf) + for dw in np.arange(0,2.5,0.01): # Delta(Deltaf) FASTER + T2= Temperature_from_Width(W0+dw,PressureBar,Diameter) + T1= Temperature_from_Width(W0,PressureBar,Diameter) + DQ = np.append(DQ,(heat_capacity_Cv_B_phase_intergral_from_0(T2, PressureBar) - heat_capacity_Cv_B_phase_intergral_from_0(T1, PressureBar)) * volume * 6.242e+18) # [eV] + DW = np.append(DW,dw) + + # Fit line to extract the slope alpha: DQ = alpha * DW + global alpha + alpha, _ = np.polyfit(DW, DQ, 1) + + # Input delta width from input energy + deltawidth = E/alpha + + return deltawidth, alpha + + +# Define signal fit function Dw vs time +def df(_t,_fb,_d): # time, base width, delta (delta width) + _t_w = 1/(np.pi*_fb) + _t1=0 + return _fb + np.heaviside(_t-_t1,1)*(_d*np.power(t_b/_t_w,_t_w/(t_b-_t_w))*(t_b/(t_b-_t_w))*(np.exp(-(_t-_t1)/t_b) - np.exp(-(_t-_t1)/_t_w))) + +# Define signal fit function Dw vs time +def winkelmann(_t,_fb,_d): # time, base width, delta (delta width) + _t_w = 1/(np.pi*_fb) + _t1=0 + return _fb + np.heaviside(_t-_t1,1)*(_d*np.power(t_b/_t_w,_t_w/(t_b-_t_w))*(t_b/(t_b-_t_w))*(np.exp(-(_t-_t1)/t_b))) + +########################################################### + +def Toy(_energy,_pressure,_temperature,_diameter): + + print() + # Input delta(width) from input energy + delta, _ = DeltaWidth_from_Energy(_energy,_pressure,_temperature,_diameter) + + # Base width from the input base temperature + f_base = Width_from_Temperature(_temperature,_pressure,_diameter) + + # Response time + t_w = 1/(np.pi*f_base) + + print("Pressure: ",_pressure, " bar") + print("Temperature: ",_temperature*1e6, " uk") + print("Diameter: ",_diameter*1e9, " nm") + print("T/Tc: ",_temperature/temperature_critical_superfluid(_pressure)) + + print("Base width: ",f_base*1000, " mHz") + print("Width variation: ",delta*1000, " mHz") + print("t_w: ",t_w, "s") + print("Bandwidth: ",np.pi*f_base/2, " Hz") + + t = np.linspace(-5, 50, 1000) # time + + base_toy = np.array([]) # base width distribution + delta_toy = np.array([]) # delta width distribution + + # Repeat the fit N times + for i in range(N): + + # Delta f vs time + #print(np.exp(-(t-5)/t_b)) + #print(np.exp(-(t-5)/t_w)) + t1=0 + deltaf = f_base + np.heaviside(t-t1,1)*(delta*np.power(t_b/t_w,t_w/(t_b-t_w))*(t_b/(t_b - t_w))*(np.exp(-(t-t1)/t_b) - np.exp(-(t-t1)/t_w))) + + # Add noise based on voltage error. The voltage error comes from the LOCKIN circuit. + for j in range(len(deltaf)): + deltaf[j] = np.random.normal(deltaf[j],np.power( deltaf[j],2 ) / (v_h*f_base)*v_rms, 1) + + # Fit the noise'd distribution + popt, pcov = curve_fit(df,t,deltaf) + # errors on base width and width increase + base_fit, delta_fit = popt + + delta_toy = np.append(delta_toy,delta_fit) + base_toy = np.append(base_toy,base_fit) + + + if (verbose): + # Plot deltaf(t) + #plt.title(str(_pressure)+' bar - '+str(diameter*1e9)+' nm - '+str(_temperature*1e6)+" $\mu$K") + plt.figure(figsize=(8, 6)) + plt.plot(t, deltaf*1000,linestyle='',marker='.', color="black") + plt.plot(t, df(t,*popt)*1000) + plt.plot(t, winkelmann(t,*popt)*1000, color='red', linestyle='--') + plt.axhline(y=f_base*1000, color='green', linestyle='--') + plt.xlim([-5,30]) + plt.xlabel('time [s]') + plt.ylabel('$\Delta f$ [mHz]') + plt.savefig('deltaf_toy-example400.0.pdf') + plt.show() + + # Plot toy energy distribution + if verbose and _energy==10000: + plt.hist(delta_toy*alpha/1000,100) + plt.xlabel('Energy [KeV]') + plt.ylabel('Entries') + plt.savefig('energy_distribution'+str(d*1e9)+'.pdf') + plt.show() + + # Plot toy base distribution + if verbose and _energy==10000: + plt.hist(base_toy,100) + plt.xlabel('Base width [Hz]') + plt.ylabel('Entries') + plt.savefig('base_width_distribution'+str(d*1e9)+'.pdf') + plt.show() + + + ## Gaussian fit for base and delta toy distributions + (delta_mu, delta_sigma) = norm.fit(delta_toy) + (base_mu, base_sigma) = norm.fit(base_toy) + + print("Fitted base: ",base_mu," ",base_sigma) + print("Fitted delta: ",delta_mu," ",delta_sigma) + + # Calculates alpha_prime for the error propagation + epsilon=1e-9 + _, alpha1 = DeltaWidth_from_Energy(_energy, _pressure, _temperature - epsilon/2, _diameter) + _, alpha2 = DeltaWidth_from_Energy(_energy, _pressure, _temperature + epsilon/2, _diameter) + gap = energy_gap_in_low_temp_limit(_pressure) + alpha_prime = (alpha2-alpha1)/epsilon * Boltzmann_const/gap * np.power(_temperature,2) * 1/Width_from_Temperature(_temperature,_pressure,_diameter) + print("alpha_prime",alpha_prime) + + print(alpha_prime*delta_mu*base_sigma) + print(alpha*delta_sigma) + + _error= np.sqrt( np.power(alpha_prime*delta_mu*base_sigma,2) + np.power(alpha*delta_sigma,2) ) + #energy_error= np.sqrt( np.power(alpha_prime*delta_mu*base_sigma,2) ) + + return _error + + +def Run_Toy_Pressure(start, end, step, _pressure,_ttc,_diameter, _f): + + global error + global value + + for v in np.arange(start, end, step): + _temperature=_ttc*temperature_critical_superfluid(v) + sigma = Toy(energy,v,_temperature,_diameter) + print("Error:", sigma,"eV, ", sigma/energy*100,"%") + + error = np.append(error,sigma/energy*100) + value = np.append(value,v) + + print(v,*(sigma/energy*100),file=_f) + + +def Run_Toy_Diameter(start, end, step, _pressure,_ttc,_diameter, _f): + + global error + global value + + _temperature=_ttc*temperature_critical_superfluid(_pressure) + for v in np.arange(start, end, step): + sigma = Toy(energy,_pressure,_temperature,v) + print("Error:", sigma,"eV, ", sigma/energy*100,"%") + + error = np.append(error,sigma/energy*100) + value = np.append(value,v) + + print(v,*(sigma/energy*100),file=_f) + + +def Run_Toy_Temperature(start, end, step, _pressure,_temperature,_diameter, _f): + + global error + global value + + for v in np.arange(start, end, step): + sigma = Toy(energy,_pressure,v*temperature_critical_superfluid(_pressure),_diameter) + + print("Error:", sigma,"eV, ", sigma/energy*100,"%") + + error = np.append(error,sigma/energy*100) + value = np.append(value,v) + + print(v,*(sigma/energy*100),file=_f) + + +########################################################### + + +if __name__ == "__main__": + + Toy(1000,0,0.11*temperature_critical_superfluid(0),400e-9) + + exit() + + # Output files + f1 = open("output/lockin_toy-error-pressure.txt", "w") + f2 = open("output/lockin_toy-error-diameter.txt", "w") + f3 = open("output/lockin_toy-error-temperature.txt", "w") + print("# pressure[bar]","error[%]",file=f1) + print("# diameter[m]","error[%]",file=f2) + print("# T/Tc","error[%]",file=f3) + + print("Gap:", energy_gap_in_low_temp_limit(pressure)/temperature_critical_superfluid(pressure)/Boltzmann_const) + + # Parameters used + print("\nEnergy : ",energy, " eV") + #print("Pressure: ",pressure, "bar") + #print("T/Tc: ",temperature/temperature_critical_superfluid(pressure)) + + global error + global value + error = np.array([]) + value = np.array([]) + + # Starts the toy simulation for a range of PRESSURES, fixed T/Tc and diameter + print("\nPRESSURE...") + #diameter = 400e-9; # [m] vibrating wire diameter + #ttc=0.1 # T/Tc + + Run_Toy_Pressure(0, 30, 1, unused, ttc, diameter, f1) + f1.close() + + # Plot results + plt.title(str(diameter*1e9)+' nm - T/Tc='+str(ttc)) + plt.plot(value,error, linestyle='', marker='o', color="black") + plt.xlabel('Pressure [bar]') + plt.ylabel('Error [%]') + plt.yscale('log') + plt.savefig('lockin-error-pressure.pdf') + plt.savefig('lockin-error-pressure.png') + plt.show() + + + + # Starts the toy simulation for a range of DIAMETERS, fixed T/Tc and pressure + print("\nDIAMETERS...") + #pressure = 0 # [bar] pressure + #ttc=0.1 # T/Tc + + error = np.array([]) + value = np.array([]) + #Run_Toy_Diameter(50e-9, 1000e-9, 100e-9, pressure, ttc, unused, f2) # error at 4.5um - 0.17 + Run_Toy_Diameter(150e-9, diameter_max, 100e-9, pressure, ttc, unused, f2) + f2.close() + + # Plot results + plt.title(str(pressure)+' bar - T/Tc='+str(ttc)) + plt.plot(value*1e9,error, linestyle='', marker='o', color="black") + plt.xlabel('diameter [nm]') + plt.ylabel('Error [%]') + plt.yscale('log') + plt.savefig('lockin-error-diameter.pdf') + plt.savefig('lockin-error-diameter.png') + plt.show() + + + + # Starts the toy simulation for a range of T/Tc + print("\nT/Tc...") + #pressure = 0 # [bar] pressure + #diameter = 400e-9; # [m] vibrating wire diameter + + error = np.array([]) + value = np.array([]) + Run_Toy_Temperature(0.1, 0.18, 0.01, pressure, unused, diameter, f3) + f3.close() + + # Plot results + plt.title(str(pressure)+' bar - '+str(diameter*1e9)+' nm') + plt.plot(value,error, linestyle='', marker='o', color="black") + plt.xlabel('T/T$_c$') + plt.ylabel('Error [%]') + plt.yscale('log') + plt.savefig('lockin-error-temperature.pdf') + plt.savefig('lockin-error-temperature.png') + plt.show() + diff --git a/toys/plot_for_paper_shot.py b/toys/plot_for_paper_shot.py index ddb633b..891862d 100755 --- a/toys/plot_for_paper_shot.py +++ b/toys/plot_for_paper_shot.py @@ -53,7 +53,7 @@ def heat_capacity(T,A): ax1.text(-0.2, 1, 'a', ha='left', va='top', transform=ax1.transAxes, fontsize=16, weight='bold') # label outside the plot #ax1.text(0.05, 0.90, 'a', transform=ax1.transAxes, fontsize=16, verticalalignment='top', horizontalalignment='left', weight='bold') # label inside the plot -ax1.plot(x1,y1,label='Lock-in amplifier',color='dodgerblue',linestyle='-') +ax1.plot(x1,y1,label='Conventional readout',color='dodgerblue',linestyle='-') ax1.plot(x2,y2,label='SQUID readout',color='red',linestyle='--') ax1.plot(x3,y3,label='QP shot noise',color='green',linestyle=':') ax1.set_xlabel('Energy [eV]') @@ -84,7 +84,7 @@ def heat_capacity(T,A): ax2.text(-0.2, 1, 'b', ha='left', va='top', transform=ax2.transAxes, fontsize=16, weight='bold') #ax2.text(0.05, 0.90, 'b', transform=ax2.transAxes, fontsize=16, verticalalignment='top', horizontalalignment='left', weight='bold') -ax2.plot(x1,y1,label='Lock-in amplifier',color='dodgerblue',linestyle='-') +ax2.plot(x1,y1,label='Conventional readout',color='dodgerblue',linestyle='-') ax2.plot(x2,y2,label='SQUID readout',color='red',linestyle='--') ax2.plot(x3,y3,label='QP shot noise',color='green',linestyle=':') ax2.plot(x1, heat_capacity(np.array(x1), A),linestyle=':', marker='',color="black") @@ -134,7 +134,7 @@ def heat_capacity(T,A): ax1.text(-0.2, 1, 'a', ha='left', va='top', transform=ax1.transAxes, fontsize=16, weight='bold') #ax1.text(0.05, 0.90, 'a', transform=ax1.transAxes, fontsize=16, verticalalignment='top', horizontalalignment='left', weight='bold') -ax1.plot(x1*1e9,y1,label='Lock-in amplifier',color='dodgerblue',linestyle='-') +ax1.plot(x1*1e9,y1,label='Conventional readout',color='dodgerblue',linestyle='-') ax1.plot(x2*1e9,y2,label='SQUID readout',color='red',linestyle='--') ax1.plot(x3*1e9,y3,label='QP shot noise',color='green',linestyle=':') ax1.set_xlabel('Diameter [nm]') @@ -161,7 +161,7 @@ def heat_capacity(T,A): ax2.text(-0.2, 1, 'b', ha='left', va='top', transform=ax2.transAxes, fontsize=16, weight='bold') #ax2.text(0.05, 0.90, 'b', transform=ax2.transAxes, fontsize=16, verticalalignment='top', horizontalalignment='left', weight='bold') -ax2.plot(x1,y1,label='Lock-in amplifier',color='dodgerblue',linestyle='-') +ax2.plot(x1,y1,label='Conventional readout',color='dodgerblue',linestyle='-') ax2.plot(x2,y2,label='SQUID readout',color='red',linestyle='--') ax2.plot(x3,y3,label='QP shot noise',color='green',linestyle=':') ax2.set_xlabel('Pressure [bar]')