""" This module contains classes for importing, filtering and analyzing raw I-V and IF data obtained from SIS mixer experiments. Two classes (``RawData0`` and ``RawData``) are provided to help manage the data. ``RawData0`` is intended for data that was collected with no LO injection (i.e., unpumped data), and ``RawData`` is intended for data that was collected with LO injection (i.e., pumped data). Note: Experimental data can be passed to these classes either in the form of CSV data files or Numpy arrays. In both bases, the data should have two columns: the first for voltage, and the second for current or IF power, depending on the file type. For CSV files, you can define the delimiter using the keyword argument ``delimiter=','``, the number of rows to skip for the header using ``skip_header=1``, and which columns to import using ``usecols=(0,1)``. Take a look at the data in ``QMix/notebooks/eg-230-data/`` for an example. Also, take a look at ``QMix/notebooks/analyze-experimental-data.ipynb`` for an example of how to use this module. """ import glob import os from copy import deepcopy import matplotlib.pyplot as plt import numpy as np import scipy.constants as sc from scipy import special from scipy.interpolate import UnivariateSpline import qmix from qmix.exp.if_data import dcif_data, if_data from qmix.exp.if_response import if_response from qmix.exp.iv_data import dciv_curve, iv_curve from qmix.exp.parameters import params as PARAMS from qmix.harmonic_balance import harmonic_balance from qmix.mathfn.filters import gauss_conv from qmix.mathfn.misc import slope_span_n from qmix.misc.terminal import cprint from qmix.qtcurrent import qtcurrent from qmix.respfn import RespFnFromIVData # Colors for plotting _pale_blue = '#6ba2f9' _pale_green = '#42b173' _pale_red = '#f96b6b' _red = 'r' _blue = '#1f77b4' _dark_blue = '#1f77b4' # Impedance recovery parameters _good_error = 7e-7 _step = 1e-5 # Parameters for saving figures _plot_params = {'dpi': 500, 'bbox_inches': 'tight'} # Note: All plotting functions are excluded from coverage tests # by using: "# pragma: no cover" # FILE HIERARCHY ------------------------------------------------------------- _file_structure = {'DC IV data': '01_dciv/', 'Pumped IV data': '02_iv_curves/', 'IF data': '03_if_data/', 'Impedance recovery': '04_impedance/', 'IF noise': '05_if_noise/', 'Noise temperature': '06_noise_temp/', 'IF spectrum': '07_spectrum/', 'Overall performance': '08_overall_performance/', 'CSV data': '09_csv_data/'} """Default file hierarchy to use when plotting experimental data.""" # CLASSES FOR RAW DATA ------------------------------------------------------- class RawData0(object): """Class for importing and analyzing experimental DC data (with no LO). Note: Experimental data can be passed to this class either in the form of CSV data files or Numpy arrays. In both cases, the data should have two columns: the first for voltage, and the second for current or IF power, depending on the file type. For CSV files, you can define the delimiter using the keyword argument ``delimiter=','``, the number of rows to skip for the header using ``skip_header=1``, and which columns to import using ``usecols=(0,1)``. Take a look at the data in ``QMix/notebooks/eg-230-data/`` for an example. Also, take a look at ``QMix/notebooks/analyze-experimental-data.ipynb`` for an example of how to use this module. See ``qmix.exp.parameters.params`` for all possible keyword arguments. These parameters control how the data is imported and analyzed. Args: dciv: DC I-V curve. Either a CSV data file or a Numpy array. The data should have two columns: the first for voltage, and the second for current. If you are using CSV files, the properties of the CSV file can be set through additional keyword arguments (see below). dcif: DC IF data. Either a CSV data file or a Numpy array. The data should have two columns: the first for voltage, and the second for IF power. If you are using CSV files, the properties of the CSV file can be set through additional keyword arguments (see below). Keyword arguments: delimiter (str): Delimiter for CSV files. usecols (tuple): List of columns to import (tuple of length 2). skip_header (int): Number of rows to skip, used to skip the header. v_fmt (str): Units for voltage ('uV', 'mV', or 'V'). i_fmt (str): Units for current ('uA', 'mA', or 'A'). vmax (float): Maximum voltage to import in units [V]. npts (int): Number of points to have in I-V interpolation. debug (bool): Plot each step of the I-V processing procedure. voffset (float): Voltage offset, in units [V]. ioffset (float): Current offset, in units [A]. voffset_range (float): Voltage range over which to search for offset, in units [V]. voffset_sigma (float): Standard deviation of Gaussian filter when searching for offset. rseries (float): Series resistance in experimental measurement system, in units [ohms]. i_multiplier (float): Multiply the imported current by this value. v_multiplier (float): Multiply the imported voltage by this value. ifdata_vmax (float): Maximum IF voltage to import. ifdata_npts (int): Number of points for interpolation. filter_data (bool): Filter data vgap_guess (float): Guess of gap voltage. Used to temporarily normalize while filtering. Given in units [V]. igap_guess (float): Guess of gap current. Used to temporarily normalize while filtering. Given in units [A]. filter_theta (float): Angle by which to the rotate data while filtering. Given in radians. filter_nwind (int): Window size for Savitsky-Golay filter. filter_npoly (int): Order of Savitsky-Golay filter. ifdata_sigma (float): Standard deviation of Gaussian used for filtering, in units [V] area (float): Area of the junction in um^2. vgap_threshold (float): The current to measure the gap voltage at. rn_vmin (float): Lower voltage range to determine the normal resistance rn_vmax (float): Upper voltage range to determine the normal resistance vrsg (float): The voltage at which to calculate the subgap resistance. vleak (float): The voltage at which to calculate the subgap leakage current. vshot (list): Voltage range over which to fit shot noise slope, in units [V]. Can be a list of lists to define multiple ranges. comment (str): Comment to describe this instance. verbose (bool): Print to terminal. """ def __init__(self, dciv, dcif=None, **kw): # Import keyword arguments tmp = deepcopy(PARAMS) tmp.update(kw) kw = tmp # Unpack keyword arguments comment = kw['comment'] v_smear = kw['v_smear'] vleak = kw['vleak'] area = kw['area'] verbose = kw['verbose'] self.kwargs = kw self.comment = comment self.vleak = vleak if isinstance(dciv, str): # input type: CSV file self.file_path = dciv elif isinstance(dciv, np.ndarray): # input type: Numpy array self.file_path = 'Numpy array' else: raise ValueError('Input data type not recognized.') # Get DC I-V data self.voltage, self.current, self.dc = dciv_curve(dciv, **kw) # Unpack DC I-V metadata self.vgap = self.dc.vgap self.igap = self.dc.igap self.fgap = self.dc.fgap self.rn = self.dc.rn self.rsg = self.dc.rsg self.q = self.rsg / self.rn self.rna = self.rn * area * 1e-12 self.jc = self.vgap / self.rna self.offset = self.dc.offset self.vint = self.dc.vint self.ileak = np.interp(vleak / self.vgap, self.voltage, self.current) * self.igap # Generate response function from DC I-V curve self.resp = RespFnFromIVData(self.voltage, self.current, check_error=False, verbose=False, v_smear=None) # Generate smeared response function from DC I-V curve self.resp_smear = RespFnFromIVData(self.voltage, self.current, check_error=False, verbose=False, v_smear=v_smear) # Import DC IF data (if it exists) if dcif is not None: # Import self.if_data, dcif = dcif_data(dcif, self.dc, **kw) # Unpack DC IF metadata self.dcif = dcif self.if_noise = dcif.if_noise self.corr = dcif.corr self.shot_slope = dcif.shot_slope self.if_fit = dcif.if_fit else: # pragma: no cover self.dcif = None self.if_data = None self.if_noise = None self.corr = None self.shot_slope = None self.if_fit = None if verbose: print(self) def __str__(self): # pragma: no cover message = "\033[35m\nDC I-V data:\033[0m {0}\n".format(self.comment) message += "\tVgap: \t\t{:6.2f}\tmV\n".format(self.vgap * 1e3) message += "\tfgap: \t\t{:6.2f}\tGHz\n".format(self.fgap / 1e9) message += "\n" message += "\tRn: \t\t{:6.2f}\tohms\n".format(self.rn) message += "\tRsg: \t\t{:6.2f}\tohms\n".format(self.rsg) message += "\tQ: \t\t{:6.2f}\n".format(self.q) message += "\n" message += "\tJc: \t\t{:6.2f}\tkA/cm^2\n".format(self.jc / 1e7) message += "\tIleak: \t\t{:6.2f}\tuA\n".format(self.ileak * 1e6) message += "\n" message += "\tOffset:\t\t{:6.2f}\tmV\n".format(self.offset[0] * 1e3) message += "\t \t\t{:6.2f}\tuA\n".format(self.offset[1] * 1e6) message += "\n" message += "\tVint: \t\t{:6.2f}\tmV\n".format(self.vint * 1e3) if self.if_noise is not None: message += "\tIF noise:\t{:6.2f}\tK\n".format(self.if_noise) return message def __repr__(self): # pragma: no cover msg = "DC I-V curve: Vgap = {:.2f} mV, Rn = {:.2f} ohms" return msg.format(self.vgap / 1e-3, self.rn) def print_info(self): # pragma: no cover """Print information about the DC I-V curve. This method is deprecated. Just use ``print(dciv)`` instead, assuming that ``dciv`` is an instance of this class. """ print(self) def plot_dciv(self, fig_name=None, ax=None, vmax_plot=4., **kw): # pragma: no cover """Plot DC I-V curve. Some additional labels will be added as well, including normal-state resistance, subgap resistance, gap voltage, and gap current. Note: If ``fig_name`` is provided, this method will save the plot to the specified folder and then close the plot. This means that the Matplotlib axis object will not be returned in this case. This is done to prevent too many plots from being open at the time. Args: fig_name (str): figure filename ax: Matplotlib axis vmax_plot (float): max voltage to include in plot (in mV) kw: keyword arguments (not used) """ # Check plotting arguments # assert fig_name is None or ax is None # Unnormalize the data mv = self.vgap * 1e3 # norm -> mV ua = self.vgap / self.rn * 1e6 # norm -> uA v_mv = self.voltage * mv i_ua = self.current * ua v_v = v_mv / 1000 # Other values to plot rn_slope = -self.vint / self.rn * 1e6 + self.voltage * ua i_at_gap = np.interp([1.], self.voltage, self.current) * ua i_leak = np.interp(self.vleak / self.vgap, self.voltage, self.current) # Fit sub-gap resistance mask = (self.vleak - 0.2e-3 < v_v) & (v_v < self.vleak + 0.2e-3) psg = np.polyfit(v_mv[mask], i_ua[mask], 1) # Strings for legend labels lgd_str1 = 'DC I-V' lgd_str2 = r'$R_\mathrm{{n}}$ = %.2f $\Omega$' % self.rn lgd_str3 = r'$V_\mathrm{{gap}}$ = %.2f mV' % (self.vgap * 1e3) lgd_str4 = r'$I_\mathrm{{leak}}$ = %.2f $\mu$A' % (i_leak * ua) lgd_str5 = r'$R_\mathrm{{sg}}$ = %.1f $\Omega$' % self.rsg # Plot DC I-V curve if ax is None: fig, ax = plt.subplots() else: fig = ax.get_figure() ax.plot(v_mv, i_ua, label=lgd_str1) # Label gap ax.plot(self.vgap * 1e3, i_at_gap, marker='o', ls='None', color='r', mfc='None', markeredgewidth=1, label=lgd_str3) # Label leakage current ax.plot(2, i_leak * ua, marker='o', ls='None', color='g', mfc='None', markeredgewidth=1, label=lgd_str4) # Fit line to normal resistance slope ax.plot(v_mv, rn_slope, 'k--', label=lgd_str2) # Fit line to sub-gap resistance ax.plot(v_mv, np.polyval(psg, v_mv), 'k:', label=lgd_str5) # Set max voltage idx = np.abs(vmax_plot - v_mv).argmin() ax.set_xlim([0, v_mv[idx]]) ax.set_ylim([0, i_ua[idx]]) # Set figure properties ax.set_xlabel(r'Bias Voltage (mV)') ax.set_ylabel(r'Current ($\mu$A)') ax.minorticks_on() ax.legend(loc=2) if fig_name is not None: fig.savefig(fig_name, **_plot_params) plt.close(fig) return else: return ax def plot_offset(self, fig_name=None, ax=None, **kw): # pragma: no cover """Plot DC I-V curve at the origin to see if there is an offset. Note: If ``fig_name`` is provided, this method will save the plot to the specified folder and then close the plot. This means that the Matplotlib axis object will not be returned in this case. This is done to prevent too many plots from being open at the time. Args: fig_name (str): figure filename ax: Matplotlib axis kw: keyword arguments (not used) """ # Check arguments # assert fig_name is None or ax is None # Unnormalize the data mv = self.vgap * 1e3 # norm -> mV ua = self.vgap / self.rn * 1e6 # norm -> uA v_mv = self.voltage * mv i_ua = self.current * ua # Only plot around the origin mask = (-0.2 < v_mv) & (v_mv < 0.2) # Plot offset (to make sure it was corrected properly) if ax is None: fig, ax = plt.subplots() else: fig = ax.get_figure() ax.plot(v_mv[mask], i_ua[mask]) ax.set_xlim([-0.2, 0.2]) ax.set_xlabel(r'Bias Voltage (mV)') ax.set_ylabel(r'Current ($\mu$A)') ax.minorticks_on() ax.grid() if fig_name is not None: fig.savefig(fig_name, **_plot_params) plt.close(fig) return else: return ax def plot_rdyn(self, fig_name=None, ax=None, vmax_plot=4., **kw): # pragma: no cover """Plot dynamic resistance of the DC I-V curve. The dynamic resistance is the derivative of the I-V data, inverted to get resistance. Note: If ``fig_name`` is provided, this method will save the plot to the specified folder and then close the plot. This means that the Matplotlib axis object will not be returned in this case. This is done to prevent too many plots from being open at the time. Args: fig_name: figure filename ax: Matplotlib axis vmax_plot: max voltage to include in plot (in mV) kw (dict): keyword arguments (not used) """ # Check arguments # assert not fig_name is not None and ax is not None # De-normalize v_mv = self.voltage * self.vgap * 1e3 i_ma = self.current * self.igap * 1e3 # Calculate dynamic resistance rdyn = slope_span_n(v_mv, i_ma, 11) # Plot dynamic resistance if ax is None: fig, ax = plt.subplots() else: fig = ax.get_figure() ax.semilogy(v_mv, 1 / rdyn) ax.set_xlabel(r'Bias Voltage (mV)') ax.set_ylabel(r'Dynamic Resistance ($\Omega$)') ax.set_xlim([0, vmax_plot]) ax.minorticks_on() if fig_name is not None: fig.savefig(fig_name, **_plot_params) plt.close(fig) return else: return ax def plot_rstat(self, fig_name=None, ax=None, vmax_plot=4., **kw): # pragma: no cover """Plot static resistance of DC I-V data. The static resistance is the DC voltage divided by the DC current. Note: If ``fig_name`` is provided, this method will save the plot to the specified folder and then close the plot. This means that the Matplotlib axis object will not be returned in this case. This is done to prevent too many plots from being open at the time. Args: fig_name: figure filename ax: Matplotlib axis vmax_plot: max voltage to include in plot (in mV) kw (dict): keyword arguments (not used) """ # Check arguments # assert fig_name is not None and ax is not None # De-normalize v_mv = self.voltage * self.vgap * 1e3 i_ma = self.current * self.igap * 1e3 # Only plot up to a given voltage mask = (0 < v_mv) & (v_mv < vmax_plot) v_mv, i_ma = v_mv[mask], i_ma[mask] # Calculate static resistance r_stat = v_mv / i_ma # Plot static resistance if ax is None: fig, ax = plt.subplots() else: fig = ax.get_figure() ax.plot(v_mv, r_stat) ax.set_xlabel(r'Bias Voltage (mV)') ax.set_ylabel(r'Static Resistance ($\Omega$)') ax.set_xlim([0, vmax_plot]) ax.set_ylim(bottom=0) ax.minorticks_on() if fig_name is not None: fig.savefig(fig_name, **_plot_params) plt.close(fig) return else: return ax def plot_if_noise(self, fig_name=None, ax=None, **kw): # pragma: no cover """Plot IF noise. The IF noise is calculated from the slope of the shot noise. Note: If ``fig_name`` is provided, this method will save the plot to the specified folder and then close the plot. This means that the Matplotlib axis object will not be returned in this case. This is done to prevent too many plots from being open at the time. Args: fig_name: figure filename ax: Matplotlib axis (should be tuple with two axes) kw (dict): keyword arguments (not used) """ # Check arguments # assert fig_name is not None and ax is not None if self.if_data is None: print('No DC IF data loaded.\n') return # Line that fits to normal-state resistance rslope = (self.voltage * self.vgap / self.rn - self.vint / self.rn) * 1e6 # Denormalize data mv = self.vgap * 1e3 # norm -> mV ua = self.igap * 1e6 # norm -> uA v_mv = self.voltage * mv i_ua = self.current * ua # vmax = v_mv.max() # imax = i_ua.max() if ax is None: fig, (ax1, ax2) = plt.subplots(2, sharex=True, figsize=(6, 9)) else: ax1, ax2 = ax fig = ax1.get_figure() plt.subplots_adjust(hspace=0., wspace=0.) # Plot DC I-V curve ax1.plot(v_mv, i_ua, label='DC I-V') ax1.plot(v_mv, rslope, 'k--', label=r'$R_\mathrm{{n}}^{{-1}}$ slope') ax1.axvline(self.vint * 1e3, c='k', ls=':', lw=0.5, label=r'$V_\mathrm{{int}}$') ax1.set_ylabel(r'Current ($\mu$A)') ax1.set_ylim(bottom=0) ax1.set_xlim(left=0) ax1.legend(loc=4, frameon=False) # Plot DC IF data v_mv = self.if_data[:, 0] * mv ax2.plot(v_mv, self.if_data[:, 1], _pale_red, label='IF (unpumped)') ax2.plot(self.shot_slope[:, 0] * self.vgap * 1e3, self.shot_slope[:, 1], 'k--', label='Shot noise slope') ax2.plot(self.vint * 1e3, self.if_noise, marker='o', ls='None', color='r', mfc='None', markeredgewidth=1, label='IF Noise: {0:.2f} K'.format(self.if_noise)) ax2.axvline(self.vint * 1e3, c='k', ls=':', lw=0.5) ax2.set_xlabel('Bias Voltage (mV)') ax2.set_ylabel('IF Power (K)') ax2.set_xlim([0, v_mv.max()]) ax2.set_ylim(bottom=0) ax2.legend(loc=4, frameon=False) if fig_name is not None: fig.savefig(fig_name, **_plot_params) plt.close(fig) return else: return ax1, ax2 def plot_all(self, fig_folder, sub_folder=None, **kw): # pragma: no cover """Plot all DC data and save it to a specified directory. This method will call the following methods: ``plot_dciv``, ``plot_offset``, ``plot_rdyn`` and ``plot_if_noise``. These figures will be put in ``fig_folder/sub_folder``. Note that if ``sub_folder`` is left as ``None``, this method will use the default file structure (see ``qmix.exp.exp_data._file_structure``). If you would instead like the figures to go into ``fig_folder``, set this argument as an empty string (""). Args: fig_folder (str): directory where the figures go sub_folder (str): sub-directory where the DC I-V figures go kw: keyword arguments that will be passed to the plotting methods """ # Folder for DC data if sub_folder is None: sub_folder = _file_structure['DC IV data'] folder = os.path.join(fig_folder, sub_folder) else: folder = os.path.join(fig_folder, sub_folder) # Make sure folder exists, create it if not if not os.path.exists(folder): os.makedirs(folder) # Names for figures fig1 = os.path.join(folder, 'dciv.png') fig2 = os.path.join(folder, 'dciv-offset.png') fig3 = os.path.join(folder, 'dciv-rdyn.png') fig4 = os.path.join(folder, 'dcif-shot-noise.png') # Generate plots self.plot_dciv(fig1, **kw) self.plot_offset(fig2, **kw) self.plot_rdyn(fig3, **kw) self.plot_if_noise(fig4, **kw) class RawData(object): """Class for importing and analyzing experimental pumped data (LO present). Note: Experimental data can be passed to this class either in the form of CSV data files or Numpy arrays. In both bases, the data should have two columns: the first for voltage, and the second for current or IF power, depending on the file type. For CSV files, you can define the delimiter using the keyword argument ``delimiter=','``, the number of rows to skip for the header using ``skip_header=1``, and which columns to import using ``usecols=(0,1)``. Take a look at the data in ``QMix/notebooks/eg-230-data/`` for an example. Also, take a look at ``QMix/notebooks/analyze-experimental-data.ipynb`` for an example of how to use this module. See ``qmix.exp.parameters.params`` for all possible keyword arguments. These parameters control how the data is imported and analyzed. Args: ivdata: I-V data. Either a CSV data file or a Numpy array. The data should have two columns: the first for voltage, and the second for current. If you are using CSV files, the properties of the CSV file can be set through additional keyword arguments (see below). dciv (qmix.exp.iv_data.DCIVData): DC I-V metadata Keyword arguments: delimiter (str): Delimiter for CSV files. usecols (tuple): List of columns to import (tuple of length 2). skip_header (int): Number of rows to skip, used to skip the header. v_fmt (str): Units for voltage ('uV', 'mV', or 'V'). i_fmt (str): Units for current ('uA', 'mA', or 'A'). vmax (float): Maximum voltage to import in units [V]. npts (int): Number of points to have in I-V interpolation. debug (bool): Plot each step of the I-V processing procedure. voffset (float): Voltage offset, in units [V]. ioffset (float): Current offset, in units [A]. voffset_range (float): Voltage range over which to search for offset, in units [V]. voffset_sigma (float): Standard deviation of Gaussian filter when searching for offset. rseries (float): Series resistance in experimental measurement system, in units [ohms]. i_multiplier (float): Multiply the imported current by this value. v_multiplier (float): Multiply the imported voltage by this value. ifdata_vmax (float): Maximum IF voltage to import. ifdata_npts (int): Number of points for interpolation. filter_data (bool): Filter data vgap_guess (float): Guess of gap voltage. Used to temporarily normalize while filtering. Given in units [V]. igap_guess (float): Guess of gap current. Used to temporarily normalize while filtering. Given in units [A]. filter_theta (float): Angle by which to the rotate data while filtering. Given in radians. filter_nwind (int): Window size for Savitsky-Golay filter. filter_npoly (int): Order of Savitsky-Golay filter. ifdata_sigma (float): Standard deviation of Gaussian used for filtering, in units [V] analyze_iv (bool): Analyze I-V data? analyze_if (bool): Analyze IF data? area (float): Area of the junction in um^2. freq (float): Frequency of LO signal. cut_low (float): only fit over first photon step, start at Vgap - vph + vph * cut_low cut_high: only fit over first photon step, finish at Vgap - vph * cut_high remb_range (tuple): range of embedding resistances to check, normalized to the normal-state resistance xemb_range (tuple): range of embedding reactances to check, normalized to the normal-state resistance alpha_max (float): Maximum drive level. num_b (int): Summation limits for Bessel functions. t_cold (float): Temperature of cold blackbody load. t_hot (float): Temperature of hot blackbody load. vmax_plot (float): Maximum bias voltage for plots. comment (str): Comment to describe this instance. verbose (bool): Print to terminal. """ def __init__(self, ivdata, dciv, if_hot=None, if_cold=None, **kw): # Import keyword arguments tmp = deepcopy(PARAMS) tmp.update(kw) kw = tmp self.kwargs = kw # Unpack keyword arguments comment = kw['comment'] freq = kw['freq'] analyze = kw['analyze'] analyze_if = kw['analyze_if'] analyze_iv = kw['analyze_iv'] verbose = kw['verbose'] # Analyze data? (deprecated, set individually) if analyze is not None: # pragma: no cover analyze_iv = analyze analyze_if = analyze # Sort out file paths if isinstance(ivdata, str): # input type: CSV file self.iv_file = ivdata self.directory = os.path.dirname(ivdata) self.iv_filename = os.path.basename(ivdata) if if_hot is not None and if_cold is not None: self.filename_hot = os.path.basename(if_hot) self.filename_cold = os.path.basename(if_cold) else: self.filename_hot = None self.filename_cold = None elif isinstance(ivdata, np.ndarray): # input type: Numpy array self.iv_file = 'Numpy array' self.directory = 'Numpy array' self.iv_filename = 'Numpy array' if if_hot is not None and if_cold is not None: self.filename_hot = 'Numpy array' self.filename_cold = 'Numpy array' else: self.filename_hot = None self.filename_cold = None else: raise ValueError("Input data type not recognized.") # Unpack DC I-V metadata self.dciv = dciv self.vgap = dciv.dc.vgap self.igap = dciv.dc.igap self.fgap = dciv.dc.fgap self.rn = dciv.rn self.offset = dciv.offset self.vint = dciv.vint self.dc = dciv.dc self.voltage_dc = dciv.voltage self.current_dc = dciv.current # Get LO frequency if isinstance(ivdata, np.ndarray) and freq is None: str1 = 'If input data is in the form of Numpy arrays, ' str2 = 'you must define the frequency of the LO signal.' raise ValueError(str1 + str2) self.freq, self.freq_str = _get_freq(freq, ivdata) kw['freq'] = self.freq self.vph = self.freq / self.fgap * 1e9 # photon voltage # Print to terminal if verbose: cprint('Importing: {}'.format(comment), 'HEADER') print(" -> Files:") print("\tI-V file: \t{}".format(self.iv_file)) if self.filename_hot is not None: print("\tIF hot file: \t{}".format(self.filename_hot)) if self.filename_cold is not None: print("\tIF cold file:\t{}".format(self.filename_cold)) print(" -> Frequency: {:.1f} GHz".format(self.freq)) # Import pumped I-V curve self.voltage, self.current = iv_curve(ivdata, self.dc, **kw) # Dynamic resistance of I-V curve self.rdyn = slope_span_n(self.current, self.voltage, 21) # Impedance recovery if analyze_iv: self._recover_zemb() else: # pragma: no cover self.zt = None self.vt = None self.fit_good = None self.zw = None self.alpha = None # Import and analyze IF data from hot/cold loads self.good_if_noise_fit = True if if_hot is not None and if_cold is not None and analyze_if: # Import and analyze IF data results, self.idx_best, dcif = if_data(if_hot, if_cold, self.dc, dcif=dciv.dcif, **kw) # Unpack results self.if_hot = results[:, :2] self.if_cold = np.vstack((results[:, 0], results[:, 2])).T self.tn = results[:, 3] self.gain = results[:, 4] # DC IF values self.if_noise = dcif.if_noise self.corr = dcif.corr self.shot_slope = dcif.shot_slope self.good_if_noise_fit = dcif.if_fit # Best values self.tn_best = self.tn[self.idx_best] self.gain_best = self.gain[self.idx_best] self.g_db = 10 * np.log10(self.gain[self.idx_best]) self.v_best = self.if_hot[self.idx_best, 0] # Dynamic resistance at optimal bias voltage i = np.abs(self.voltage - self.v_best).argmin() p = np.polyfit(self.voltage[i:i + 10], self.current[i:i + 10], 1) self.zj_if = self.rn / p[0] else: # pragma: no cover self.filename_hot = None self.filename_cold = None self.if_hot = None self.if_cold = None self.tn = None self.gain = None self.idx_best = None self.if_noise = None self.good_if_noise_fit = None self.shot_slope = None self.tn_best = None self.g_db = None if verbose: print("") def _recover_zemb(self): """Recover the embedding circuit (i.e., the Thevenin eqiv. circuit). The technique used here is the RF voltage match method described by Skalare (1989) and Withington et al. (1995). Note: All currents and voltages are normalized to the gap voltage and to the normal resistance, respectively. Keyword Args: fit_range (list): Fit interval for impedance recovery, normalized to the width of the first photon step. cut_low (float): only fit over first photon step, start at Vgap - vph + vph * cut_low (DEPRECATED) cut_high: only fit over first photon step, finish at Vgap - vph * cut_high (DEPRECATED) remb_range (tuple): range of embedding resistances to check, normalized to the normal-state resistance xemb_range (tuple): range of embedding reactances to check, normalized to the normal-state resistance Returns: thevenin impedance, voltage, and fit (boolean) """ # Fit interval for impedance recovery fit_range = self.kwargs.get('fit_range', None) if fit_range is not None: fit_low = fit_range[0] fit_high = fit_range[1] else: fit_low = self.kwargs.get('cut_low', None) fit_high = 1 - self.kwargs.get('cut_high', None) if fit_low is None or fit_high is None: fit_low = PARAMS['fit_range'][0] fit_high = PARAMS['fit_range'][1] fit_high = 1 - fit_high # Range of impedance values to test remb_range = self.kwargs.get('remb_range', PARAMS['remb_range']) xemb_range = self.kwargs.get('xemb_range', PARAMS['xemb_range']) # Force certain value zemb = self.kwargs.get('zemb', PARAMS['zemb']) cprint(" -> Impedance recovery:") # Unpack vgap = self.vgap # mV rn = self.rn # ohms vph = self.freq * sc.giga / self.fgap resp = self.dciv.resp # Only consider linear region of first photon step # Ratio removed at either end of step v_low = 1 - vph + vph * fit_low v_high = 1 - vph * fit_high mask = (v_low <= self.voltage) & (self.voltage <= v_high) exp_voltage = self.voltage[mask] exp_current = self.current[mask] idx_middle = np.abs(exp_voltage - (1 - vph / 2.)).argmin() # Calculate alpha for every bias voltage alpha = _find_alpha(self.dciv, exp_voltage, exp_current, vph, **self.kwargs) ac_voltage = alpha * vph # Calculate AC junction impedance ac_current = _find_ac_current(resp, exp_voltage, vph, alpha, **self.kwargs) ac_impedance = ac_voltage / ac_current zw = ac_impedance[idx_middle] # Calculate error at every embedding impedance in given range zt_real = np.linspace(remb_range[0], remb_range[1], 101) zt_imag = np.linspace(xemb_range[0], xemb_range[1], 201) err_surf = np.empty((len(zt_real), len(zt_imag)), dtype=float) for i in range(len(zt_real)): for j in range(len(zt_imag)): err_surf[i, j] = _error_function(ac_voltage, ac_impedance, zt_real[i] + 1j * zt_imag[j]) ibest, jbest = np.unravel_index(err_surf.argmin(), err_surf.shape) zt_real_best, zt_imag_best = zt_real[ibest], zt_imag[jbest] if zemb is None: zt_best = zt_real_best + 1j * zt_imag_best else: zt_best = zemb vt_best = _find_source_voltage(ac_voltage, ac_impedance, zt_best) err_best = err_surf[ibest, jbest] # Determine whether or not it was a good fit good_fit = err_best <= _good_error # Print to terminal if good_fit: cprint('\t- good fit', 'OKGREEN') else: cprint('\t- bad fit', 'WARNING') print("\t- embedding circuit:") print("\t\t- voltage: \t{:+6.2f}\t\t* Vgap".format(vt_best)) print("\t\t- impedance: \t{:+12.2f}\t* Rn".format(zt_best)) with np.errstate(divide='ignore', invalid='ignore'): power_avail = np.abs(vt_best * vgap)**2 / 8 / np.real(zt_best * rn) print("\t\t- avail. power:\t{:+7.2f}\t\tnW".format(power_avail / 1e-9)) print("\t- junction:") print("\t\t- drive level:\t{:+6.2f}".format(alpha[idx_middle])) print("\t\t- impedance:\t{:+12.2f}\t* Rn".format(zw)) with np.errstate(divide='ignore', invalid='ignore'): power_delivered = np.abs(ac_voltage[idx_middle] * vgap)**2 / 2 / np.real(zw * rn) print("\t\t- deliv. power:\t{:+7.2f}\t\tnW".format(power_delivered / 1e-9)) # Save values as attributes self.zt = zt_best self.vt = vt_best self.fit_good = good_fit self.zw = zw self.alpha = alpha[idx_middle] self.err_surf = err_surf def plot_iv(self, fig_name=None, ax=None, vmax_plot=4.): # pragma: no cover """Plot pumped I-V curve. Note: If ``fig_name`` is provided, this method will save the plot to the specified folder and then close the plot. This means that the Matplotlib axis object will not be returned in this case. This is done to prevent too many plots from being open at the time. Args: fig_name: figure filename ax: Matplotlib axis vmax_plot: max voltage to include in plot (in mV) """ # Check arguments # assert fig_name is not None and ax is not None # DC I-V curve dciv = self.dciv # De-normalize vmv = dciv.vgap * 1e3 # norm -> mV iua = dciv.igap * 1e6 # norm -> uA imax = np.interp(vmax_plot, dciv.voltage * vmv, dciv.current * iua) if ax is None: fig, ax = plt.subplots() else: fig = ax.get_figure() ax.plot(dciv.voltage * vmv, dciv.current * iua, label="Unpumped") ax.plot(self.voltage * vmv, self.current * iua, 'r', label="Pumped") ax.set_xlabel(r'Bias Voltage (mV)') ax.set_ylabel(r'DC Current (uA)') ax.set_xlim([0, vmax_plot]) ax.set_ylim([0, imax]) msg = 'LO: {:.1f} GHz'.format(self.freq) ax.legend(loc=2, title=msg, frameon=False) ax.minorticks_on() if fig_name is not None: fig.savefig(fig_name, **_plot_params) plt.close(fig) return else: return ax def plot_if(self, fig_name=None, ax=None, vmax_plot=4.): # pragma: no cover """Plot IF power from hot and cold loads. Note: If ``fig_name`` is provided, this method will save the plot to the specified folder and then close the plot. This means that the Matplotlib axis object will not be returned in this case. This is done to prevent too many plots from being open at the time. Args: fig_name: figure filename ax: Matplotlib axis vmax_plot: max voltage to include in plot (in mV) """ # Check arguments # assert fig_name is not None and ax is not None # De-normalize voltage v_mv = self.if_hot[:, 0] * self.vgap * 1e3 # Plot IF data if ax is None: fig, ax = plt.subplots() else: fig = ax.get_figure() ax.plot(v_mv, self.if_hot[:, 1], _pale_red, label='Hot') ax.plot(v_mv, self.if_cold[:, 1], _pale_blue, label='Cold') if self.dciv.if_data is not None: v_tmp = self.dciv.if_data[:, 0] * self.vgap * 1e3 ax.plot(v_tmp, self.dciv.if_data[:, 1], 'k--', label='No LO') ax.set_xlabel('Bias Voltage (mV)') ax.set_ylabel('IF Power (K)') ax.set_ylim(bottom=0) ax.set_xlim([0, vmax_plot]) msg = 'LO: {:.1f} GHz'.format(self.freq) ax.legend(loc=1, title=msg, frameon=False) ax.minorticks_on() if fig_name is not None: fig.savefig(fig_name, **_plot_params) plt.close(fig) return else: return ax def plot_ivif(self, fig_name=None, ax=None, vmax_plot=4.): # pragma: no cover """Plot IV and IF data on same plot. Note: If ``fig_name`` is provided, this method will save the plot to the specified folder and then close the plot. This means that the Matplotlib axis object will not be returned in this case. This is done to prevent too many plots from being open at the time. Args: fig_name: figure filename ax: Matplotlib axis (should be tuple with two axes) vmax_plot: max voltage to include in plot (in mV) """ # Check arguments # assert fig_name is not None and ax is not None # De-normalize mv = self.vgap * 1e3 ua = self.vgap / self.rn * 1e6 imax = np.interp(vmax_plot, self.voltage_dc * mv, self.current_dc * ua) if ax is None: fig, ax1 = plt.subplots() else: ax1, ax2 = ax fig = ax1.get_figure() # Plot I-V data ax1.plot(self.voltage_dc * mv, self.current_dc * ua, '#8c8c8c', label="Unpumped") ax1.plot(self.voltage * mv, self.current * ua, 'k', label="Pumped") ax1.set_xlabel('Bias Voltage (mV)') ax1.set_ylabel(r'DC Current ($\mu$A)') ax1.set_ylim([0, imax]) ax1.legend(loc=2, frameon=True, framealpha=1.) ax1.grid(False) # Plot IF data mv = self.if_hot[:, 0] * self.vgap * 1e3 ax2 = ax1.twinx() ax2.plot(mv, self.if_hot[:, 1], '#f96b6b', label='Hot') ax2.plot(mv, self.if_cold[:, 1], '#6ba2f9', label='Cold') ax2.set_ylabel('IF Power (K)') ax2.legend(loc=1, framealpha=1., frameon=True) ax2.grid(False) ax2.set_ylim(bottom=0) ax2.set_xlim([0, vmax_plot]) ax1.minorticks_on() ax2.minorticks_on() if fig_name is not None: fig.savefig(fig_name, **_plot_params) plt.close(fig) return else: return ax1, ax2 def plot_shapiro(self, fig_name=None, ax=None): # pragma: no cover """Plot Shapiro steps. Note: If ``fig_name`` is provided, this method will save the plot to the specified folder and then close the plot. This means that the Matplotlib axis object will not be returned in this case. This is done to prevent too many plots from being open at the time. Args: fig_name: figure filename ax: Matplotlib axis """ # Check arguments # assert fig_name is not None and ax is not None # De-normalize voltage v_mv = self.if_hot[:, 0] * self.vgap * 1e3 # Calculate Shaprio voltage separation f_hz = float(self.freq) * 1e9 vshapiro = f_hz * sc.h / sc.e / 2 / sc.milli mask = (0. < v_mv) & (v_mv < 3.5 * vshapiro) if ax is None: fig, ax = plt.subplots() else: fig = ax.get_figure() ax.plot(v_mv[mask], self.if_hot[mask, 1], '#f96b6b', label='Hot') ax.plot(v_mv[mask], self.if_cold[mask, 1], '#6ba2f9', label='Cold') ax.axvline(vshapiro, label=r'$\omega_\mathrm{LO}h/2e$', c='k', ls='--') ax.axvline(2 * vshapiro, c='k', ls='--') ax.axvline(3 * vshapiro, c='k', ls='--') ax.set_xlim([0, 3.5 * vshapiro]) ax.set_xlabel('Bias Voltage (mV)') ax.set_ylabel('IF Power (K)') ax.minorticks_on() ax.legend() if fig_name is not None: fig.savefig(fig_name, **_plot_params) plt.close(fig) return else: return ax def plot_if_noise(self, fig_name=None, ax=None): # pragma: no cover """Plot IF noise. Note: If ``fig_name`` is provided, this method will save the plot to the specified folder and then close the plot. This means that the Matplotlib axis object will not be returned in this case. This is done to prevent too many plots from being open at the time. Args: fig_name: figure filename ax: Matplotlib axis (should be tuple with two axes) """ # Check arguments # assert fig_name is not None and ax is not None rslope = (self.voltage_dc * self.vgap / self.rn - self.vint / self.rn) * 1e6 vmax = self.voltage.max() * self.vgap * 1e3 imax = self.current.max() * self.igap * 1e6 if ax is None: fig, (ax1, ax2) = plt.subplots(2, sharex=True, figsize=(6, 9)) else: ax1, ax2 = ax fig = ax1.get_figure() plt.subplots_adjust(hspace=0., wspace=0.) ax1.plot(self.dciv.voltage * self.vgap * 1e3, self.dciv.current * self.igap * 1e6, label='Unpumped') ax1.plot(self.voltage * self.vgap * 1e3, self.current * self.igap * 1e6, 'r', label='Pumped') ax1.plot(self.voltage * self.vgap * 1e3, rslope, 'k--', label=r'$R_\mathrm{{n}}^{{-1}}$ slope') ax1.plot(self.vint * 1e3, 0, 'ro', label=r'$V_\mathrm{{int}}$') ax1.set_ylabel(r'Current ($\mu$A)') ax1.set_ylim([0, imax]) ax1.set_xlim([0, vmax]) ax1.legend() v_mv = self.if_hot[:, 0] * self.vgap * 1e3 ax2.plot(v_mv, self.if_hot[:, 1], _pale_red, label='Hot') ax2.plot(v_mv, self.if_cold[:, 1], _pale_blue, label='Cold') ax2.plot(self.shot_slope[:, 0] * self.vgap * 1e3, self.shot_slope[:, 1], 'k--', label='Shot noise slope') ax2.plot(self.vint * 1e3, self.if_noise, 'ro', label='IF Noise: {:.2f} K'.format(self.if_noise)) ax2.set_xlabel('Bias Voltage (mV)') ax2.set_ylabel('IF Power (K)') ax2.set_xlim([0, v_mv.max()]) ax2.set_ylim([0, np.max(self.shot_slope) * 1.1]) ax2.legend(loc=0) if fig_name is not None: fig.savefig(fig_name, **_plot_params) plt.close(fig) return else: return ax1, ax2 def plot_noise_temp(self, fig_name=None, ax=None, vmax_plot=4.): # pragma: no cover """Plot noise temperature. Note: If ``fig_name`` is provided, this method will save the plot to the specified folder and then close the plot. This means that the Matplotlib axis object will not be returned in this case. This is done to prevent too many plots from being open at the time. Args: fig_name: figure filename ax: Matplotlib axis (should be tuple with two axes) vmax_plot: max voltage to include in plot (in mV) """ # Check arguments # assert fig_name is not None and ax is not None v_mv = self.if_hot[:, 0] * self.vgap * 1e3 hot = gauss_conv(self.if_hot[:, 1], 5) cold = gauss_conv(self.if_cold[:, 1], 5) if ax is None: fig, ax1 = plt.subplots() ax2 = ax1.twinx() else: ax1, ax2 = ax fig = ax1.get_figure() # Plot IF power from hot/cold loads l1 = ax1.plot(v_mv, hot, _pale_red, label='Hot IF') l2 = ax1.plot(v_mv, cold, _pale_blue, label='Cold IF') ax1.set_xlabel('Bias Voltage (mV)') ax1.set_ylabel('IF Power (K)') ax1.set_xlim([1, 3.5]) ax1.set_ylim([0, hot.max() * 1.3]) ax1.grid(False) # Plot noise temperature l3 = ax2.plot(v_mv, self.tn, _pale_green, ls='--', label='Noise Temp.') l4 = ax2.plot(v_mv[self.idx_best], self.tn_best, label=r'$T_\mathrm{{n}}={:.1f}$ K'.format(self.tn_best), marker='o', ls='None', color='k', mfc='None', markeredgewidth=1) ax2.set_ylabel('Noise Temperature (K)', color='g') ax2.set_ylim([0, self.tn_best * 5]) ax2.set_xlim([0., vmax_plot]) for tl in ax2.get_yticklabels(): tl.set_color('g') ax2.grid(False) # Build legend lns = l1 + l2 + l3 + l4 labs = [l.get_label() for l in lns] ax2.legend(lns, labs, loc=2, frameon=True, framealpha=1.) if fig_name is not None: fig.savefig(fig_name, **_plot_params) plt.close(fig) return else: return ax1, ax2 def plot_yfac_noise_temp(self, fig_name=None, ax=None, vmax_plot=4.): # pragma: no cover """Plot Y-factor and noise temperature. Note: If ``fig_name`` is provided, this method will save the plot to the specified folder and then close the plot. This means that the Matplotlib axis object will not be returned in this case. This is done to prevent too many plots from being open at the time. Args: fig_name: figure filename ax: Matplotlib axis vmax_plot: max voltage to include in plot (in mV) """ # Check arguments # assert fig_name is not None and ax is not None # Plot y-factor and noise temperature v_mv = self.if_hot[:, 0] * self.vgap * 1e3 hot = gauss_conv(self.if_hot[:, 1], 5) cold = gauss_conv(self.if_cold[:, 1], 5) yfac = hot / cold if ax is None: fig, ax1 = plt.subplots() ax2 = ax1.twinx() else: ax1, ax2 = ax fig = ax1.get_figure() # Plot y factor ax1.plot(v_mv, yfac, _dark_blue, label='Y-factor') ax1.axhline(293. / 77., c=_dark_blue, ls=':') ax1.set_xlabel('Bias Voltage (mV)') ax1.set_ylabel('Y-factor', color=_dark_blue) for tl in ax1.get_yticklabels(): tl.set_color(_dark_blue) ax1.set_ylim([1., 4.]) # Plot noise temperature ax2.plot(v_mv, self.tn, _red, label='Noise Temp.') ax2.plot(v_mv[self.idx_best], self.tn_best, marker='o', ls='None', color='k', mfc='None', markeredgewidth=1) msg = '{0:.1f} K'.format(self.tn_best) ax2.annotate(msg, xy=(v_mv[self.idx_best], self.tn_best), xytext=(v_mv[self.idx_best] + 0.5, self.tn_best + 50), bbox=dict(boxstyle="round", fc="w", alpha=0.5), arrowprops=dict(color='black', arrowstyle="->", lw=1), ) ax2.set_ylabel('Noise Temperature (K)', color=_red) for tl in ax2.get_yticklabels(): tl.set_color(_red) ax2.set_ylim([0, 300.]) ax2.set_xlim([0., vmax_plot]) if fig_name is not None: fig.savefig(fig_name, **_plot_params) plt.close(fig) return else: return ax1, ax2 def plot_gain_noise_temp(self, fig_name=None, ax=None, vmax_plot=4.): # pragma: no cover """Plot gain and noise temperature. Note: If ``fig_name`` is provided, this method will save the plot to the specified folder and then close the plot. This means that the Matplotlib axis object will not be returned in this case. This is done to prevent too many plots from being open at the time. Args: fig_name: figure filename ax: Matplotlib axis vmax_plot: max voltage to include in plot (in mV) """ # Check arguments # assert fig_name is not None and ax is not None # De-normalize voltage v_mv = self.if_hot[:, 0] * self.vgap * 1e3 if ax is None: fig, ax1 = plt.subplots() ax2 = ax1.twinx() else: ax1, ax2 = ax fig = ax1.get_figure() # Plot gain ax1.plot(v_mv, self.gain, label=r'Gain', color=_dark_blue) ax1.plot(v_mv[self.idx_best], self.gain[self.idx_best], marker='o', ls='None', color='k', mfc='None', markeredgewidth=1) msg = r'$G_\mathrm{{c}}={0:.2f}$'.format(self.gain[self.idx_best]) ax1.annotate(msg, xy=(v_mv[self.idx_best], self.gain[self.idx_best]), xytext=(v_mv[self.idx_best] + 0.75, self.gain[self.idx_best] - 0.1), arrowprops=dict(color='black', arrowstyle="->", lw=0.5), va="center", ha="left", fontsize=16, ) ax1.set_xlabel('Bias Voltage (mV)') ax1.set_ylabel('Gain', color=_dark_blue) for tl in ax1.get_yticklabels(): tl.set_color(_dark_blue) ax1.set_ylim(bottom=0) ax1.minorticks_on() # Plot noise temperature ax2.plot(v_mv, self.tn, _red, label='Noise Temp.') ax2.plot(v_mv[self.idx_best], self.tn_best, marker='o', ls='None', color='k', mfc='None', markeredgewidth=1) msg = r'$T_\mathrm{{N}}={0:.1f}$ K'.format(self.tn_best) ax2.annotate(msg, xy=(v_mv[self.idx_best], self.tn_best), xytext=(v_mv[self.idx_best] + 0.75, self.tn_best + 50), arrowprops=dict(color='black', arrowstyle="->", lw=0.5), va="center", ha="left", fontsize=16) ax2.set_ylabel('Noise Temperature (K)', color=_red) for tl in ax2.get_yticklabels(): tl.set_color(_red) ax2.set_ylim([0, self.tn_best * 5]) ax2.set_xlim([0., vmax_plot]) ax2.minorticks_on() if fig_name is not None: fig.savefig(fig_name, **_plot_params) plt.close(fig) return else: return ax1, ax2 def plot_rdyn(self, fig_name=None, ax=None, vmax_plot=4.): # pragma: no cover """Plot dynamic resistance. Note: If ``fig_name`` is provided, this method will save the plot to the specified folder and then close the plot. This means that the Matplotlib axis object will not be returned in this case. This is done to prevent too many plots from being open at the time. Args: fig_name: figure filename ax: Matplotlib axis vmax_plot: max voltage to include in plot (in mV) """ # Check arguments # assert fig_name is not None and ax is not None # Unnormalize current/voltage v_mv = self.voltage * self.vgap * 1e3 # Determine dynamic resistance (remove 0 values to avoid /0 errors) rdyn = self.rdyn * self.rn # Position of steps steps = np.r_[-1 + self.vph * np.arange(-3, 4, 1), 1 - self.vph * np.arange(3, -4, -1)] v_steps = steps * self.vgap * 1e3 r_steps = np.interp(v_steps, v_mv, rdyn) # Dynamic resistance at 'best' bias point (where TN is best) vb_best = (self.if_hot[:, 0] * self.vgap * 1e3)[self.idx_best] rdyn_bias = np.interp(vb_best, v_mv, rdyn) if ax is None: fig, ax = plt.subplots() else: fig = ax.get_figure() ax.plot(v_mv, rdyn, label=r'$R_\mathrm{dyn}$') ax.plot(vb_best, rdyn_bias, 'r^', label=r'%.1f $\Omega$' % rdyn_bias) ax.plot(v_steps, r_steps, 'k+', label=r'$V_\mathrm{gap} + nV_\mathrm{ph}$') plt.axvline(-1 * self.vgap * 1e3, c='k', ls='--', lw=0.5) plt.axvline(0, c='k', ls='--', lw=0.5) plt.axvline(1 * self.vgap * 1e3, c='k', ls='--', lw=0.5) ax.set_xlabel('Bias Voltage (mV)') ax.set_ylabel(r'Dynamic Resistance ($\Omega$)') ax.set_xlim([0, vmax_plot]) ax.set_ylim(bottom=0) ax.legend(loc=0, title='LO: ' + str(self.freq) + ' GHz') ax.minorticks_on() if fig_name is not None: fig.savefig(fig_name, **_plot_params) plt.close(fig) return else: return ax def plot_gain(self, fig_name=None, ax=None, vmax_plot=4.): # pragma: no cover """Plot gain. Note: If ``fig_name`` is provided, this method will save the plot to the specified folder and then close the plot. This means that the Matplotlib axis object will not be returned in this case. This is done to prevent too many plots from being open at the time. Args: fig_name: figure filename ax: Matplotlib axis vmax_plot: max voltage to include in plot (in mV) """ # Check arguments # assert fig_name is not None and ax is not None # De-normalize voltage v_mv = self.if_hot[:, 0] * self.vgap * 1e3 if ax is None: fig, ax = plt.subplots() else: fig = ax.get_figure() ax.plot(v_mv, self.gain*100, label=r'$G_{{c}}$') ax.set_xlabel('Bias Voltage (mV)') ax.set_ylabel('Gain (%)') ax.set_xlim([0, vmax_plot]) ax.set_ylim([0, self.gain.max() * 105]) ax.minorticks_on() if fig_name is not None: fig.savefig(fig_name, **_plot_params) plt.close(fig) return else: return ax def plot_error_surface(self, fig_name=None, ax=None): # pragma: no cover """Plot error surface (from impedance recovery). Note: If ``fig_name`` is provided, this method will save the plot to the specified folder and then close the plot. This means that the Matplotlib axis object will not be returned in this case. This is done to prevent too many plots from being open at the time. Args: fig_name: figure filename ax: Matplotlib axis """ # Check arguments # assert fig_name is not None and ax is not None # Range of impedance values (normalized) remb_range = self.kwargs.get('remb_range', PARAMS['remb_range']) xemb_range = self.kwargs.get('xemb_range', PARAMS['xemb_range']) # Range of impedance values (de-normalized) zt_real = np.linspace(remb_range[0], remb_range[1], 101) * self.rn zt_imag = np.linspace(xemb_range[0], xemb_range[1], 201) * self.rn zt_real_range = zt_real[-1] - zt_real[0] zt_imag_range = zt_imag[-1] - zt_imag[0] # Recovered impedance zt_best = self.zt * self.rn zt_re_best, zt_im_best = zt_best.real, zt_best.imag # Mesh (for plotting) xx, yy = np.meshgrid(zt_real, zt_imag) zz = np.log10(self.err_surf) if ax is None: fig, ax = plt.subplots() else: fig = ax.get_figure() pc = ax.pcolor(xx, yy, zz.T, cmap='viridis') # Add color bar cbar = plt.colorbar(pc, ax=ax) cbar.ax.set_ylabel(r'$\log_{10}(\varepsilon)$', rotation=90) # Left or right half of plot? zt_re_mid = (zt_real[0] + zt_real[-1]) / 2. if zt_re_best > zt_re_mid: text_posx = zt_re_best - zt_real_range / 10. text_ha = "right" else: text_posx = zt_re_best + zt_real_range / 10. text_ha = "left" # Top or bottom half of plot? zt_im_mid = (zt_imag[0] + zt_imag[-1]) / 2. if zt_im_best > zt_im_mid: text_posy = zt_im_best - zt_imag_range / 10. text_va = "top" else: text_posy = zt_im_best + zt_imag_range / 10. text_va = "bottom" # Annotate best value err_str1 = 'Minimum Error at\n' err_str2 = r'$Z_\mathrm{{T}}$={0:.2f} $\Omega$'.format(zt_best) err_str = err_str1 + err_str2 text_pos = text_posx, text_posy bbox_props = dict(boxstyle="round", fc="w", alpha=0.5) ax.annotate(err_str, xy=(zt_re_best, zt_im_best), xytext=text_pos, bbox=bbox_props, va=text_va, ha=text_ha, fontsize=16, arrowprops=dict(color='black', arrowstyle="->", lw=2)) ax.set_xlabel(r'$R_\mathrm{{T}}$ ($\Omega$)') ax.set_ylabel(r'$X_\mathrm{{T}}$ ($\Omega$)') # Add text box textstr1 = 'Embedding impedance:\n' textstr2 = r'$Z_\mathrm{{T}}=R_\mathrm{{T}}+j\,X_\mathrm{{T}}$' textstr = textstr1 + textstr2 ax.text(0.05, 0.95, textstr, transform=ax.transAxes, verticalalignment='top', bbox=bbox_props) if fig_name is not None: fig.savefig(fig_name, **_plot_params) plt.close(fig) return else: return ax def plot_simulated(self, fig_name=None, ax=None, vmax_plot=4.): # pragma: no cover """Plot simulated I-V curve (from impedance recovery). Note: If ``fig_name`` is provided, this method will save the plot to the specified folder and then close the plot. This means that the Matplotlib axis object will not be returned in this case. This is done to prevent too many plots from being open at the time. Args: fig_name: figure filename ax: Matplotlib axis vmax_plot: max voltage to include in plot (in mV) """ # Check arguments # assert fig_name is not None and ax is not None # Unpack vph = self.freq * sc.giga / self.dciv.fgap resp = self.dciv.resp_smear # Fit interval for impedance recovery fit_range = self.kwargs.get('fit_range', None) if fit_range is not None: fit_low = fit_range[0] fit_high = fit_range[1] else: fit_low = self.kwargs.get('cut_low', None) fit_high = self.kwargs.get('cut_high', None) if fit_low is None or fit_high is None: fit_low = PARAMS['fit_range'][0] fit_high = PARAMS['fit_range'][1] fit_high = 1 - fit_high v_min = (1 - vph + vph * fit_low) * self.vgap * 1e3 v_max = (1 - vph * fit_high) * self.vgap * 1e3 # Build embedding circuit cct = qmix.circuit.EmbeddingCircuit(1, 1) cct.freq[1] = vph cct.zt[1, 1] = self.zt cct.vt[1, 1] = self.vt # Simulate pumped I-V curve vj = harmonic_balance(cct, resp, num_b=30, verbose=False) vph_list = [0, cct.freq[1]] current = qtcurrent(vj, cct, resp, vph_list, num_b=30, verbose=False) # De-normalize mv = self.vgap * 1e3 # norm -> mV ua = self.igap * 1e6 # norm -> uA if ax is None: fig, ax = plt.subplots() else: fig = ax.get_figure() ax.plot(self.dciv.voltage * mv, self.dciv.current * ua, label='Unpumped', c='gray') ax.plot(self.voltage * mv, self.current * ua, label='Pumped') ax.plot(cct.vb * mv, current[0].real * ua, label='Simulated', c='r', ls='--') ax.plot([v_min, v_max], np.interp([v_min, v_max], cct.vb * mv, current[0].real * ua), 'k+', label='Fit Interval') ax.set_xlim([0, vmax_plot]) ax.set_ylim([0, np.interp(vmax_plot, self.dciv.voltage * mv, self.dciv.current * ua)]) ax.set_xlabel(r'Bias Voltage (mV)') ax.set_ylabel(r'DC Current (uA)') msg1 = 'LO: {:.1f} GHz'.format(self.freq) msg2 = r'$V_T$ = {:.2f} mV'.format(self.vt * self.vgap * 1e3) msg3 = r'$Z_T$ = {:.2f} $\Omega$'.format(self.zt * self.rn) msg = msg1 + '\n' + msg2 + '\n' + msg3 ax.legend(title=msg, frameon=False) if fig_name is not None: fig.savefig(fig_name, **_plot_params) plt.close(fig) return else: return ax def plot_all(self, fig_folder, file_struc=None, **kw): # pragma: no cover """Plot all pumped data and save to specified directory. This method will call the following methods: ``plot_iv``, ``plot_ivif``, ``plot_error_surface``, ``plot_simulated``, ``plot_simulated``, ``plot_noise_temp``, and ``plot_gain_noise_temp``. The plots generated by the different methods will be placed in different sub-directories. This is set by the ``file_struc`` argument. This argument is a dictionary with the following keys: ``'Pumped IV data'``, ``'IF data'``, ``'Impedance recovery'``, and ``'Noise temperature'``. Set ``file_struc`` to an empty string ("") if you want all the figures to go into ``fig_folder``. Args: fig_folder (str): folder where the figures go file_struc (dict): dictionary listing all the sub-folders kw: keyword arguments that will be passed to the plotting methods """ # Folders for new plots if file_struc is None: iv_folder = _file_structure['Pumped IV data'] if_folder = _file_structure['IF data'] zr_folder = _file_structure['Impedance recovery'] tn_folder = _file_structure['Noise temperature'] elif isinstance(file_struc, str) and file_struc == '': iv_folder = '' if_folder = '' zr_folder = '' tn_folder = '' elif isinstance(file_struc, dict): iv_folder = file_struc['Pumped IV data'] if_folder = file_struc['IF data'] zr_folder = file_struc['Impedance recovery'] tn_folder = file_struc['Noise temperature'] else: raise ValueError iv_folder = os.path.join(fig_folder, iv_folder) if_folder = os.path.join(fig_folder, if_folder) zr_folder = os.path.join(fig_folder, zr_folder) tn_folder = os.path.join(fig_folder, tn_folder) # Make sure folders exist, create them if not if not os.path.exists(iv_folder): os.makedirs(iv_folder) if not os.path.exists(if_folder): os.makedirs(if_folder) if not os.path.exists(zr_folder): os.makedirs(zr_folder) if not os.path.exists(tn_folder): os.makedirs(tn_folder) # Names for figures f = str(self.freq) fig1 = os.path.join(iv_folder, f + '-iv.png') fig2 = os.path.join(if_folder, f + '-ivif.png') fig3 = os.path.join(zr_folder, f + '-err-surf.png') fig4 = os.path.join(zr_folder, f + '-sim.png') fig5 = os.path.join(tn_folder, f + '-tn.png') fig6 = os.path.join(tn_folder, f + '-tn-gain.png') # Generate plots self.plot_iv(fig1, **kw) self.plot_ivif(fig2, **kw) self.plot_noise_temp(fig5, **kw) self.plot_gain_noise_temp(fig6, **kw) if self.zt is not None: self.plot_error_surface(fig3) self.plot_simulated(fig4, **kw) # ANALYZE IF SPECTRUM DATA --------------------------------------------------- def plot_if_spectrum(data_folder, fig_folder=None, figsize=None): # pragma: no cover """Plot all IF spectra within data_folder. Args: data_folder: data folder fig_folder: figure folder figsize: figure size, in inches """ pstr = "\nImporting and plotting IF data:" cprint(pstr, 'HEADER') if_spectra_files = glob.glob(os.path.join(data_folder, '*comb*.dat')) fig1, ax1 = plt.subplots(figsize=figsize) fig2, ax2 = plt.subplots(figsize=figsize) for if_file in if_spectra_files: filename = os.path.basename(if_file)[:-4] print(" - {}".format(filename)) base = filename.split('_')[0][1:] freq, t_n, p_hot_db, p_cold_db = if_response(if_file) fig2, ax2 = plt.subplots(figsize=figsize) ax2.plot(freq, t_n) ax1.plot(freq, t_n, label="{} GHz".format(base)) ax2.plot(freq, gauss_conv(t_n, sigma=1), label="{} GHz".format(base)) ax2.set_ylabel('Noise Temperature (K)') ax2.set_xlabel('Frequency (GHz)') ax2.set_ylim([0, 400]) ax2.set_xlim([0, 20]) if fig_folder is not None: figname = os.path.join(fig_folder, filename) fig2.savefig(figname + '.png', **_plot_params) ax2.set_ylim([0, 2000]) fig2.savefig(figname + '2.png', **_plot_params) else: fig2.show() ax1.set_ylabel('Noise Temperature (K)') ax1.set_xlabel('Frequency (GHz)') ax1.set_ylim([0, 500]) ax1.set_xlim([0, 20]) ax1.legend() fig1.savefig(os.path.join(fig_folder, 'if_spectra.png'), **_plot_params) ax1.set_ylim([0, 2000]) fig1.savefig(os.path.join(fig_folder, 'if_spectra2.png'), **_plot_params) ax2.set_ylabel('Noise Temperature (K)') ax2.set_xlabel('Frequency (GHz)') ax2.set_ylim([0, 500]) ax2.set_xlim([0, 20]) ax2.legend() fig2.savefig(os.path.join(fig_folder, 'if_spectra_smooth.png'), **_plot_params) ax2.set_ylim([0, 2000]) fig2.savefig(os.path.join(fig_folder, 'if_spectra_smooth2.png'), **_plot_params) print("") # Plot overall results ------------------------------------------------------- def plot_overall_results(dciv, data_list, fig_folder, vmax_plot=4., figsize=None, tn_max=None, f_range=None): # pragma: no cover """Plot all results. This function is somewhat messy, but it will take in a list of RawData class instances, and plot the overall figures of merit (e.g., noise temperature vs LO frequency). Args: dciv: DC I-V data (instance of RawData0) data_list: list of pumped data (instances of RawData) fig_folder: figure destination """ cprint("\nPlotting results.", 'HEADER') # Initialize results directory (if needed) for sub_folder in _file_structure.values(): path = os.path.join(fig_folder, sub_folder) if not os.path.exists(path): os.makedirs(path) plotparam = dict(ls='--', marker='o') csv_folder = os.path.join(fig_folder, '09_csv_data/') fig_folder = os.path.join(fig_folder, '08_overall_performance/') num_data = float(len(data_list)) # Gather data as a function of LO frequency freq, t_n, gain, rdyn = [], [], [], [] f_z, z, v, aemb = [], [], [], [] f_z_all, z_all, v_all, aemb_all = [], [], [], [] if_noise_f, if_noise = [], [] aj, zj = [], [] for data in data_list: freq.append(data.freq) t_n.append(data.tn_best) gain.append(data.g_db) aj.append(data.alpha) zj.append(data.zw) rdyn.append(data.zj_if) if data.fit_good: f_z.append(data.freq) z.append(data.zt * data.rn) v.append(data.vt * data.vgap) aemb.append(data.vt / data.vph) f_z_all.append(data.freq) z_all.append(data.zt * data.rn) v_all.append(data.vt * data.vgap) aemb_all.append(data.vt / data.vph) if data.good_if_noise_fit: if_noise_f.append(data.freq) if_noise.append(data.if_noise) f_z = np.array(f_z) z = np.array(z) f_z_all = np.array(f_z_all) z_all = np.array(z_all) t_n = np.array(t_n) gain = np.array(gain) v = np.array(v) # For normalizing data mv = dciv.vgap * 1e3 ua = dciv.igap * 1e6 imax_plot = np.interp(vmax_plot, dciv.voltage * mv, dciv.current * ua) # Save data in text format ----------------------------------------------- # Save DC I-V curve as csv output_text = np.vstack((dciv.voltage, dciv.current)).T np.savetxt(os.path.join(csv_folder, 'dciv-data.txt'), output_text) # Save impedance as csv with open(os.path.join(csv_folder, 'recovered-emb.txt'), 'w') as fout: for i in range(len(data_list)): vt_tmp = data_list[i].vt * dciv.vgap * 1e3 zt_tmp = data_list[i].zt * dciv.rn zt_fmt = "{:6.2f} + 1j * ({:6.2f})".format(zt_tmp.real, zt_tmp.imag) pstring = '{0}\t{1:5.2f}\t{2}\t{3}\n'.format(data_list[i].freq, vt_tmp, zt_fmt, data_list[i].fit_good) fout.write(pstring) # Write all DC I-V data to a file with open(os.path.join(csv_folder, 'dciv-info.txt'), 'w') as fout: fout.write('Gap voltage \t{:6.2f} [mV]\n'.format(dciv.vgap * 1e3)) fout.write('Normal resistance\t{:6.2f} [ohms]\n'.format(dciv.rn)) fout.write('Gap frequency \t{:6.2f} [GHz]\n'.format(dciv.fgap / 1e9)) # Write all pumped data to a file with open(os.path.join(csv_folder, 'results.txt'), 'w') as fout: headers = ['Frequency (GHz)', 'IV Filename', 'IF Filename (Hot)', 'IF Filename (Cold)', 'Noise Temperature (K)', 'Gain (dB)', 'Drive Level', 'Embedding Impedance (ohms)', 'Embedding Voltage (mV)', 'Embedding Circuit Recovered', 'IF Noise (K)', 'IF Noise Recovered'] fout.write(', '.join(headers) + '\n') for data in data_list: _list = [data.freq, data.iv_file, data.filename_hot, data.filename_cold, "{:.2f}".format(data.tn_best), "{:.2f}".format(data.g_db), "{:.2f}".format(data.alpha), "{:.4f}".format(data.zt * dciv.rn), "{:.4f}".format(data.vt * dciv.vgap * 1e3), data.fit_good, data.if_noise, data.good_if_noise_fit] string = ', '.join([str(item) for item in _list]) fout.write(string + '\n') # Plot all pumped iv curves ---------------------------------------------- fig, ax = plt.subplots(figsize=figsize) ax.plot(dciv.voltage * mv, dciv.current * ua, 'k') for i, data in enumerate(data_list): ax.plot(data.voltage * mv, data.current * ua, color=plt.cm.winter(i / num_data), label=data.freq) if f_range is not None: ax.set_xlim([f_range[0], f_range[1]]) ax.set_xlabel(r'Voltage (mV)') ax.set_ylabel(r'Current (uA)') ax.set_xlim([0, vmax_plot]) ax.set_ylim([0, imax_plot]) ax.legend(title='LO (GHz)', frameon=False) fig.savefig(os.path.join(fig_folder, 'iv_curves.png'), dpi=500) plt.close(fig) # Plot dynamic resistance ------------------------------------------------ fig, ax = plt.subplots(figsize=figsize) ax.plot(freq, rdyn, **plotparam) if f_range is not None: ax.set_xlim([f_range[0], f_range[1]]) ax.set_xlabel('Frequency (GHz)') ax.set_ylabel(r'Dynamic resistance ($\Omega$)') ax.set_ylim(bottom=0) fig.savefig(os.path.join(fig_folder, 'rdyn.png'), dpi=500) plt.close(fig) # Plot noise temperature results ----------------------------------------- fig, ax = plt.subplots(figsize=figsize) ax.plot(freq, t_n, color=_blue, **plotparam) ax.set_xlabel('Frequency (GHz)') ax.set_ylabel('Noise Temperature (K)') if f_range is not None: ax.set_xlim([f_range[0], f_range[1]]) if tn_max is None: ax.set_ylim(bottom=0) else: ax.set_ylim([0, tn_max]) ax.grid() fig.savefig(os.path.join(fig_folder, 'noise_temperature.png'), dpi=500) plt.close(fig) # Plot noise temperature with spline fit --------------------------------- fig, ax = plt.subplots(figsize=figsize) freq_t = np.linspace(np.min(freq), np.max(freq), 1001) sp_1 = UnivariateSpline(freq, t_n) ax.plot(freq, t_n, 'o', color=_blue) ax.plot(freq_t, sp_1(freq_t), '--', color=_blue) ax.set_xlabel('Frequency (GHz)') ax.set_ylabel('Noise Temperature (K)') if f_range is not None: ax.set_xlim([f_range[0], f_range[1]]) if tn_max is None: ax.set_ylim(bottom=0) else: ax.set_ylim([0, tn_max]) ax.grid() fname = os.path.join(fig_folder, 'noise_temperature_spline_fit.png') fig.savefig(fname, dpi=500) plt.close(fig) # Plot noise temperature results on log scale ---------------------------- fig, ax = plt.subplots(figsize=figsize) ax.semilogy(freq, t_n, color=_blue, **plotparam) ax.set_xlabel('Frequency (GHz)') ax.set_ylabel('Noise Temperature (K)') if f_range is not None: ax.set_xlim([f_range[0], f_range[1]]) ax.grid() fig.savefig(os.path.join(fig_folder, 'noise_temperature_logy.png'), dpi=500) plt.close(fig) # Plot noise temperature and gain ---------------------------------------- fig, ax1 = plt.subplots(figsize=figsize) ax1.plot(freq, t_n, c=_pale_red, **plotparam) ax1.set_xlabel('Frequency (GHz)') ax1.set_ylabel('Noise Temperature (K)', color=_pale_red) if f_range is not None: ax1.set_xlim([f_range[0], f_range[1]]) if tn_max is None: ax1.set_ylim(bottom=0) else: ax1.set_ylim([0, tn_max]) for tl in ax1.get_yticklabels(): tl.set_color(_pale_red) ax2 = ax1.twinx() ax2.plot(freq, gain, c=_blue, **plotparam) ax2.set_ylabel('Gain (dB)', color=_blue) for tl in ax2.get_yticklabels(): tl.set_color(_blue) fig.savefig(os.path.join(fig_folder, 'noise_temperature_and_gain.png'), dpi=500) plt.close(fig) # Plot IF noise contribution results ------------------------------------- fig, ax = plt.subplots(figsize=figsize) ax.plot(if_noise_f, if_noise, 'o--', color=_pale_red) if f_range is not None: ax.set_xlim([f_range[0], f_range[1]]) ax.set_xlabel('Frequency (GHz)') ax.set_ylabel(r'IF Noise Contribution (K)') ax.set_ylim(bottom=0) fname = os.path.join(fig_folder, 'if_noise.png') fig.savefig(fname, dpi=500) plt.close(fig) # Plot embedding impedance results --------------------------------------- fig, ax = plt.subplots(figsize=figsize) ax.plot(f_z, z.real, c=_pale_blue, label='Real', **plotparam) ax.plot(f_z, z.imag, c=_pale_red, label='Imaginary', **plotparam) if f_range is not None: ax.set_xlim([f_range[0], f_range[1]]) ax.set_xlabel('Frequency (GHz)') ax.set_ylabel(r'Embedding Impedance ($\Omega$)') ax.legend(frameon=False) ax.minorticks_on() fig.savefig(os.path.join(fig_folder, 'embedding_impedance.png'), dpi=500) plt.close(fig) # Plot embedding impedance results --------------------------------------- fig, ax = plt.subplots(figsize=figsize) ax.plot(f_z_all, z_all.real, c=_pale_blue, label='Real', **plotparam) ax.plot(f_z_all, z_all.imag, c=_pale_red, label='Imaginary', **plotparam) if f_range is not None: ax.set_xlim([f_range[0], f_range[1]]) ax.set_xlabel('Frequency (GHz)') ax.set_ylabel(r'Embedding Impedance ($\Omega$)') ax.legend(frameon=False) ax.minorticks_on() fig.savefig(os.path.join(fig_folder, 'embedding_impedance_all.png'), dpi=500) plt.close(fig) # Plot embedding impedance results --------------------------------------- fig, ax = plt.subplots(figsize=figsize) ax.plot(f_z, v * 1e3, c=_pale_green, ls='--', marker='o') if f_range is not None: ax.set_xlim([f_range[0], f_range[1]]) ax.set_xlabel('Frequency (GHz)') ax.set_ylabel(r'Embedding Voltage (mV)') ax.set_ylim(bottom=0) ax.minorticks_on() fig.savefig(os.path.join(fig_folder, 'embedding_voltage.png'), dpi=500) plt.close(fig) # Plot embedding impedance results --------------------------------------- fig, ax = plt.subplots(figsize=figsize) ax.plot(freq, aj, c=_pale_green, **plotparam) if f_range is not None: ax.set_xlim([f_range[0], f_range[1]]) ax.set_xlabel('Frequency (GHz)') ax.set_ylabel(r'Drive Level, $\alpha$') ax.set_ylim([0, 1.2]) ax.minorticks_on() ax.grid() fig.savefig(os.path.join(fig_folder, 'drive_level.png'), dpi=500) plt.close(fig) # Plot the impedance of the SIS junction --------------------------------- fig, ax1 = plt.subplots(figsize=figsize) zj = np.array(zj) * dciv.rn ax1.plot(freq, zj.real, c=_pale_blue, label=r'Re$\{Z_J\}$', **plotparam) ax1.plot(freq, zj.imag, c=_pale_red, label=r'Im$\{Z_J\}$', **plotparam) if f_range is not None: ax1.set_xlim([f_range[0], f_range[1]]) ax1.set_xlabel('Frequency (GHz)') ax1.set_ylabel(r'Junction Impedance ($\Omega$)') ax1.set_ylim() ax1.legend(loc=0) ax1.minorticks_on() plt.savefig(os.path.join(fig_folder, 'junction_impedance.png'), dpi=500) plt.close(fig) print(" -> Done\n") # IMPEDANCE RECOVERY HELPER FUNCTIONS (PRIVATE) ------------------------------ def _error_function(vwi, zwi, zs): """Calculate error function. Equation 26 from: S. Withington, K. G. Isaak, S. Kovtonyuk, R. Panhuyzen, and T. M. Klapwijk, "Direct detection at submillimetre wavelengths using superconducting tunnel junctions," Infrared Phys. Technol., vol. 36, no. 7, pp. 1059-1075, Dec. 1995. """ err1 = np.sum(np.abs(vwi)**2) err2 = np.sum(np.abs(vwi * zwi / (zs + zwi))) err3 = np.sum(np.abs(zwi / (zs + zwi))**2) return (err1 - err2**2 / err3) / len(vwi) def _find_source_voltage(vwi, zwi, zs): """Calculate source voltage (i.e., embedding voltage). Equation 27 from: S. Withington, K. G. Isaak, S. Kovtonyuk, R. Panhuyzen, and T. M. Klapwijk, "Direct detection at submillimetre wavelengths using superconducting tunnel junctions," Infrared Phys. Technol., vol. 36, no. 7, pp. 1059-1075, Dec. 1995. """ v1 = np.sum(np.abs(vwi * zwi / (zs + zwi))) v2 = np.sum(np.abs(zwi / (zs + zwi))**2) return v1 / v2 def _find_ac_current(resp, vb, vph, alpha, num_b=20, **kw): """Calculate AC tunneling current. This is the large-signal equation from Tucker theory. """ ac_current = np.zeros_like(vb, dtype=complex) for n in range(-num_b, num_b + 1): # Bessel functions j_n = special.jv(n, alpha) j_minus = special.jv(n - 1, alpha) j_plus = special.jv(n + 1, alpha) ac_current += j_n * (j_minus + j_plus) * resp.idc(vb + n * vph) ac_current += 1j * j_n * (j_minus - j_plus) * resp.ikk(vb + n * vph) return ac_current def _find_pumped_iv_curve(resp, vb, vph, alpha, num_b=20, **kw): """Calculate DC tunneling current (from Tucker theory). """ dc_current = np.zeros_like(vb, dtype=float) for n in range(-num_b, num_b + 1): dc_current += special.jv(n, alpha)**2 * resp.idc(vb + n * vph) return dc_current def _find_alpha(dciv, vdc_exp, idc_exp, vph, alpha_max=1.5, num_b=20, **kw): """Find the drive level (alpha) at each bias voltage. """ resp = dciv.resp # Guess initial alpha value using the Bisection Method idc_tmp = _find_pumped_iv_curve(resp, vdc_exp, vph, alpha_max, num_b=num_b, **kw) idciv = resp.idc(vdc_exp) alpha = (idc_exp - idciv) / (idc_tmp - idciv) * alpha_max alpha[alpha < 0] = 0 # Refine alpha using an iterative technique alpha_step = alpha_max / 4. for it in range(15): idc_tmp = _find_pumped_iv_curve(resp, vdc_exp, vph, alpha, num_b=num_b, **kw) idc_err_tmp = idc_tmp - idc_exp alpha[idc_err_tmp > 0] -= alpha_step alpha[idc_err_tmp < 0] += alpha_step alpha[alpha < 0] = 0 alpha_step /= 2. return np.array(alpha) # FILE MANAGEMENT ------------------------------------------------------------ def initialize_dir(fig_folder): # pragma: no cover """Initialize a new directory for storing results. If you use either Args: fig_folder: desired location """ folder_list = [''] folder_list += list(_file_structure.values()) for folder in folder_list: if not os.path.exists(fig_folder + folder): os.makedirs(fig_folder + folder) print(' - Created: ' + folder) print(" ") # FILE HELPER FUNCTIONS ------------------------------------------------------ def _get_freq_from_filename(file_path): """Get frequency from filename. This is used by the ``RawData`` class if a frequency is not provided as an argument. Note: This function assumes that the only numbers in the filename are there to represent the frequency. E.g., ``f230_0_iv.csv`` will be analyzed as 230.0 GHz, but ``f230_0_iv12.csv`` will be analyzed as 230.012 GHz. More importantly, ``no15_f230_iv.csv`` will be analyzed as 152.30 GHz. This function also assumes that the first digit in the file name represents the hundreds (x100), so to represent a frequency below 100 GHz, you should use a leading zero. E.g., 85 GHz should be saved as ``f085_0_iv.csv``, or something along those lines. Args: file_path: file path Returns: float: Frequency, in units [GHz] """ filename = os.path.basename(file_path) freq_nums = [int(s) for s in list(filename) if s.isdigit()] mult = 100. freq = 0. for c in freq_nums: freq += c * mult mult /= 10. return freq def _get_freq(freq, filepath): """Get frequency. If ``freq`` is not ``None``, return ``freq``. Otherwise, if ``freq`` is ``None``, try to get it from the filename (see ``_get_freq_from_filename()`` function). Also, return the frequency as a string, so that it can be used to name output files. But, replace periods with underscores. For example, represent a frequency of 230.0 GHz as ``230_0``. Args: freq: frequency, in units GHz filepath: filename of pumped I-V data Returns: tuple: frequency float and frequency string """ if freq is None: freq = float(_get_freq_from_filename(filepath)) else: freq = float(freq) freq_str = "{0:05.1f}".format(freq) freq_str = freq_str.replace('.', '_') return freq, freq_str