import os.path, glob from os.path import join import json import h5py import pandas as pd import numpy as np import matplotlib.pyplot as pl from scipy.optimize import curve_fit from helper import * from style import * discretization_length = 0.15 # µm xSize = 508 ySize = 33 f = 90 I_m = 0.8 # 0.8 I_c I_dc = 0.9 # 0.9 I_c TimeSteps = 135 TimeStepUnit = 0.246326 # ms meta_data_dict = {'discretization_length':discretization_length, 'xSize':xSize, 'ySize':ySize, 'f':f, 'I_m':I_m, 'I':I_dc} path_data_repo = "../Fig_4" ###### Colormaps cmap = cmcrameri.cm.romaO ################################################# # barrier position # get the barrier position to draw it at the correct positon in the later plots df = pd.read_csv(join(path_data_repo, "barrier_position.csv")) def bar_pos_func(t,vc,xm): """ function giving the barrier positon, following the Shapiro protocol, see paper """ return vc * t + xm * np.sin(2 * np.pi * f * t) + 18*1e-6 times = df["time"].to_numpy()*1e-3 pos = df["bar_pos"].to_numpy() * discretization_length * 1e-6 p0 = [0.41e-3, 0.55e-6] p_opt, _ = curve_fit(bar_pos_func, times,pos ) # fig = pl.figure() pl.plot(df["bar_pos"].to_numpy() * discretization_length, df["time"].to_numpy(), label = 'Data') pl.plot(bar_pos_func(times,*p_opt)*1e6 , df["time"].to_numpy(), label = 'fit') #pl.plot(bar_pos_func(times,*p0)*1e6 , df["time"].to_numpy(), label = 'function') #pl.plot(pos*1e6 , times*1e3, label = 'function') pl.xlabel("x [µm]") pl.ylabel("Time [ms]") pl.legend() #pl.show() pl.close() ############################################################################# ############################################################################# # plot density and phase planes in the same plot ### Density ## set_plot_style() rows = 5 cols = 2 fig, ax = pl.subplots(nrows=rows,ncols=cols , figsize = (5,3), sharex=True, sharey=True, layout="constrained") fontsize = 8 times = [0, 14.3, 19.5, 26.8, 31] # times in ms ### import Thomas_fermi radius of the to get the cloud boundaries boundary = pd.read_csv(join(path_data_repo, "y_Thomas_Fermi_radius.csv")) shift = 0 boundary_x = boundary["x"].to_numpy()[int(shift / discretization_length):] - shift boundary_r_TF = boundary["r_TF"].to_numpy()[int(shift / discretization_length):] # max density vmax = 80 for i, t in enumerate(times): # open image path = os.path.join(path_data_repo, fr'Sim_density_time_evolution_t_{t:.2f}ms_slice.h5') # open files h5_file = h5py.File(path, 'r') image = h5_file['image'] # get meta_data discretization_length = h5_file.attrs['discretization_length'] xSize = h5_file.attrs['xSize'] ySize = h5_file.attrs['ySize'] f = h5_file.attrs['f'] I_m = h5_file.attrs['I_m'] I = h5_file.attrs['I'] # transform image in 1/µm^3 data = np.asarray(image) * 1/discretization_length**3 y_max = ySize * discretization_length extent = (0, data.shape[1] * discretization_length, -y_max /2 ,y_max/2) # plot data im = ax[rows-1-i, 1].imshow(data, aspect='auto', origin='lower', extent=extent, cmap=cmcrameri.cm.vik, vmin = -vmax, vmax = vmax,interpolation='none') #Thomas Fermi boundaries boundary ##upper boundary ax[rows - 1 - i, 1].fill_between(boundary_x,boundary_r_TF ,y_max/2, color = RPTU_COLORS['schiefer'], alpha = 0.9, ec = None) ax[rows - 1 - i, 1].plot(boundary_x,boundary_r_TF, color = RPTU_COLORS['schiefer'], alpha = 0.9)# + (ySize * discretization_length) / 2) ##lower boundary ax[rows - 1 - i, 1].fill_between(boundary_x, (-1) * boundary_r_TF, -y_max/2 , color = RPTU_COLORS['schiefer'], alpha = 0.9, ec = None) ax[rows - 1 - i, 1].plot(boundary_x, (-1)*boundary_r_TF,color = RPTU_COLORS['schiefer'], alpha = 0.9)# + (ySize * discretization_length) / 2, color = RPTU_COLORS['schiefer'], alpha = 0.1) # barrier time = times[i] bar_pos = bar_pos_func(time * 1e-3, *p_opt) * 1e6 print(np.round(time,1)) ax[rows - 1 - i, 1].axvline(bar_pos, ls = '--', lw = 1 , color='white', alpha = 1) ax[rows-1-i, 1].set_xlim(0,65) # comment in to get a colorbar # cbar = fig.colorbar(im, ax = ax[:, 1], shrink = 0.7) # cbar.ax.set_ylabel(r'$\Delta n$ [1/µm³]',fontsize =fontsize ) # ticks = [vmax,int(vmax/2),0,- int(vmax/2),-vmax] # cbar.ax.set_yticks(ticks) # cbar.ax.set_yticklabels(ticks, fontsize = fontsize) ###################################################################################### ###################################################################################### # Phase # max phase v_max = 1 for i, t in enumerate(times): # open image path = os.path.join(path_data_repo, fr'Sim_phase_time_evolution_t_{t:.2f}ms_slice.h5') # open files h5_file = h5py.File(path, 'r') image = h5_file['image'] # get meta_data discretization_length = h5_file.attrs['discretization_length'] xSize = h5_file.attrs['xSize'] ySize = h5_file.attrs['ySize'] f = h5_file.attrs['f'] I_m = h5_file.attrs['I_m'] I = h5_file.attrs['I'] # normalize it to pi data = np.asarray(image) / np.pi data = data[:, int(shift / discretization_length):] y_max = ySize * discretization_length extent = (0, data.shape[1] * discretization_length, -y_max /2 ,y_max/2) # plot data im = ax[rows-1-i, 0].imshow(data, aspect='auto', origin='lower', extent=extent, cmap=cmap, vmax=v_max, vmin=-v_max,interpolation='none') # boundary # upper boundary ax[rows - 1 - i, 0].fill_between(boundary_x,boundary_r_TF ,y_max/2, color = RPTU_COLORS['schiefer'], alpha = 0.9, ec = None) ax[rows - 1 - i, 0].plot(boundary_x,boundary_r_TF, color = RPTU_COLORS['schiefer'], alpha = 0.9)# + (ySize * discretization_length) / 2) ##lower boundary ax[rows - 1 - i, 0].fill_between(boundary_x, (-1) * boundary_r_TF, -y_max/2 , color = RPTU_COLORS['schiefer'], alpha = 0.9, ec = None) ax[rows - 1 - i, 0].plot(boundary_x, (-1)*boundary_r_TF,color = RPTU_COLORS['schiefer'], alpha = 0.9) time = times[i] # barrie position bar_pos = bar_pos_func(time*1e-3,*p_opt)*1e6 print(np.round(time,1)) ax[rows - 1 - i, 0].axvline(bar_pos, ls = '--', lw = 1.5 , color= 'white', alpha = 1) ax[i,0].set_ylabel(r"y [µm]", fontsize=fontsize) ticks = [-2,0,2] ax[i, 0].set_yticks(ticks) ax[i, 0].set_yticklabels(ticks, fontsize = fontsize) ax[rows - 1 - i, 0].text(x=76, y=0, s=f"{time:.1f} ms", ha='center', verticalalignment='center', fontsize=8, bbox={'boxstyle': 'round', 'fc': RPTU_COLORS["mango25pct"], 'ec': RPTU_COLORS["mango"], 'ls': '-', 'lw': 1.5}, zorder=5) ax[rows-1-i, 0].set_xlim(0,65) ax[i, 1].set_xlabel(r"x [µm]", fontsize=fontsize) ax[i, 0].set_xlabel(r"x [µm]", fontsize=fontsize) ax[i, 1].set_xticks([0, 20, 40, 60]) ax[i, 1].set_xticklabels([0, 20, 40, 60], fontsize = fontsize) ax[i, 0].set_xticks([0, 20, 40, 60]) ax[i, 0].set_xticklabels([0, 20, 40, 60], fontsize = fontsize) #fig.subplots_adjust(hspace=0.5, right=0.85) #fig.tight_layout() #fig.savefig("phase_density_profile.png") #################################################################################### #################################################################################### # Phase evolution overview # open image path = os.path.join(path_data_repo, fr'Sim_phase_time_evolution_full.h5') # open files h5_file = h5py.File(path, 'r') image = h5_file['image'] # get meta_data discretization_length = h5_file.attrs['discretization_length'] TimeSteps = h5_file.attrs['TimeSteps'] TimeStepUnit = h5_file.attrs['TimeStepUnit'] # normalize it to pi data = np.asarray(image) / np.pi extent = (0, data.shape[1] * discretization_length, 0, TimeSteps * TimeStepUnit) set_plot_style() fig = pl.figure(figsize=(2.5, 3)) fontsize = 8 # plot data im = pl.imshow(data, aspect='auto', origin='lower', extent=extent, cmap=cmap, vmax=v_max, vmin=-v_max, interpolation = 'none') # plot barrier ax = pl.gca() ax2 = ax.twinx() ax.plot(df["bar_pos"].to_numpy() * discretization_length - shift, df["time"].to_numpy(), color='white', lw = 2, ls = '--' ) pl.xlim(0,65) # comment in to get a colorbar # cbar = pl.colorbar(im, shrink = 0.7, location = 'left', orientation = 'vertical', pad = 0.25) # cbar.ax.set_ylabel(r'$\delta \phi / \pi $', fontsize=fontsize) # cbar.ax.set_yticks([1, 0.5, 0, -0.5, -1]) # cbar.ax.set_yticklabels([1, 0.5, 0, -0.5, -1], fontsize = fontsize) ax.set_yticks([0, 10, 20, 30]) ax.set_yticklabels([0, 10, 20, 30], fontsize = fontsize) ax.set_xticks([0, 20, 40, 60]) ax.set_xticklabels([0, 20, 40, 60], fontsize = fontsize) # set seconad y-axis labels ax2.set_yticks([0, 14.3, 19.5, 26.8, 31]) ax2.set_yticklabels(["0ms", "14.3ms", "19.5ms", "26.8ms", "31ms"], fontsize = fontsize-2) ax.set_ylim(0, TimeSteps * TimeStepUnit) pl.ylim(0, TimeSteps * TimeStepUnit) ax2.set_ylim(0, TimeSteps * TimeStepUnit) ax.set_xlabel(r"x [µm]", fontsize = fontsize) ax.set_ylabel(r"Time [ms]", fontsize = fontsize) pl.tight_layout() pl.savefig("phase_time_evolution.png") pl.show()