Cz flux arc fix

pycqed/measurement/waveform_control/pulse_library.py

@@ -3,7 +3,6 @@
 Library containing pulse shapes.
 '''
 from pycqed.measurement.waveform_control.pulse import Pulse, apply_modulation
-from pycqed.measurement.waveform_control_CC.waveform import martinis_flux_pulse
 from pycqed.utilities.general import int_to_bin
 
 class MW_IQmod_pulse(Pulse):
@@ -341,78 +340,7 @@ def chan_wf(self, chan, tvals):
 
 
 class MartinisFluxPulse(Pulse):
-
-    def __init__(self, channel=None, channels=None, name='Martinis flux pulse',
-                 **kw):
-        Pulse.__init__(self, name)
-        if channel is None and channels is None:
-            raise ValueError('Must specify either channel or channels')
-        elif channels is None:
-            self.channel = channel  # this is just for convenience, internally
-            # this is the part the sequencer element wants to communicate with
-            self.channels.append(channel)
-        else:
-            self.channels = channels
-
-        self.amplitude = kw.pop('amplitude', 1)
-
-        self.length = kw.pop('length', None)
-        self.lambda_coeffs = kw.pop('lambda_coeffs', None)
-
-        self.theta_f = kw.pop('theta_f', None)
-        self.g2 = kw.pop('g2', None)
-        self.E_c = kw.pop('E_c', None)
-        self.f_bus = kw.pop('f_bus', None)
-        self.f_01_max = kw.pop('f_01_max', None)
-        self.dac_flux_coefficient = kw.pop('dac_flux_coefficient', None)
-
-        self.kernel_path = kw.get('kernel_path', None)
-        if self.kernel_path is not None:
-            kernelvec = np.loadtxt(self.kernel_path)
-            self.kernel = np.zeros((len(kernelvec), len(kernelvec)))
-            for i in range(len(kernelvec)):
-                for j in range(len(kernelvec)):
-                    self.kernel[i, j] = kernelvec[i-j]
-            del(kernelvec)
-        else:
-            ValueError('Must specify kernel path')
-
-    def __call__(self, **kw):
-        self.amplitude = kw.pop('amplitude', self.amplitude)
-
-        self.pulse_buffer = kw.pop(
-            'pulse_buffer', self.pulse_buffer)
-
-        self.length = kw.pop(
-            'length', self.length)
-        self.lambda_coeffs = kw.pop('lambda_coeffs', self.lambda_coeffs)
-        self.theta_f = kw.pop('theta_f', self.theta_f)
-        self.g2 = kw.pop('g2', self.g2)
-        self.E_c = kw.pop('E_c', self.E_c)
-        self.f_bus = kw.pop('f_bus', self.f_bus)
-        self.f_01_max = kw.pop('f_01_max', self.f_01_max)
-        self.dac_flux_coefficient = kw.pop(
-            'dac_flux_coefficient', self.dac_flux_coefficient)
-
-        self.channel = kw.pop('channel', self.channel)
-        self.channels = kw.pop('channels', self.channels)
-        self.channels.append(self.channel)
-        return self
-
-    def chan_wf(self, chan, tvals):
-        t = tvals - tvals[0]
-        martinis_pulse = martinis_flux_pulse(
-            length=self.length,
-            lambda_coeffs=self.lambda_coeffs,
-            theta_f=self.theta_f, g2=self.g2,
-            E_c=self.E_c, f_bus=self.f_bus,
-            f_01_max=self.f_01_max,
-            dac_flux_coefficient=self.dac_flux_coefficient,
-            return_unit='V')
-
-        # amplitude is not used because of parameterization but
-        # we use the sign of the amplitude to set the flux compensation
-        return martinis_pulse*np.sign(self.amplitude)
+    pass
 
 
 class QWG_Codeword(Pulse):

pycqed/measurement/waveform_control_CC/waveform.py

@@ -6,6 +6,14 @@
     Prerequisites:
     Usage:
     Bugs:
+
+Contains most basic waveforms, basic means having a few parameters and a
+straightforward translation to AWG amplitude, i.e., no knowledge of qubit
+parameters.
+
+Examples of waveforms that are too advanced are flux pulses that require
+knowledge of the flux sensitivity and interaction strengths and qubit
+frequencies. See e.g., "waveform_control_CC/waveforms_flux.py".
 '''
 
 import logging
@@ -292,8 +300,9 @@ def mod_gauss_VSM(amp, sigma_length, f_modulation, axis='x', phase=0,
                                  Q_phase_delay=Q_phase_delay)
     return G_I_mod, G_Q_mod, D_I_mod, D_Q_mod
 
+
 def mod_square(amp, length, f_modulation,  phase=0,
-              motzoi=0, sampling_rate=1e9):
+               motzoi=0, sampling_rate=1e9):
     '''
     Simple modulated gauss pulse. All inputs are in s and Hz.
     '''
@@ -320,140 +329,3 @@ def mod_square_VSM(amp_G, amp_D, length, f_modulation,
     D_I_mod, D_Q_mod = mod_pulse(D_I, D_Q, f_modulation,
                                  sampling_rate=sampling_rate)
     return G_I_mod, G_Q_mod, D_I_mod, D_Q_mod
-
-
-#####################################################
-# Flux pulses
-#####################################################
-
-
-def martinis_flux_pulse(length: float, lambda_2: float, lambda_3: float,
-                        theta_f: float,
-                        f_01_max: float, J2: float,
-                        V_offset: float=0, V_per_phi0: float=1,
-                        E_c: float=250e6, f_bus: float =None,
-                        f_interaction: float =None,
-                        asymmetry: float =0, sampling_rate: float =1e9,
-                        return_unit: str='V'):
-    """
-    Returns the pulse specified by Martinis and Geller
-    Phys. Rev. A 90 022307 (2014).
-
-    \theta = \theta _0 + \sum_{n=1}^\infty  (\lambda_n*(1-\cos(n*2*pi*t/t_p))/2
-
-    note that the lambda coefficients are rescaled to ensure that the center
-    of the pulse has a value corresponding to theta_f.
-
-    length          (float) lenght of the waveform (s)
-    lambda_2
-    lambda_3
-
-    theta_f         (float) final angle of the interaction in degrees.
-                    Determines the Voltage for the center of the waveform.
-
-    f_01_max        (float) qubit sweet spot frequency (Hz).
-    J2              (float) coupling between 11-02 (Hz),
-                    approx sqrt(2) J1 (the 10-01 coupling).
-    E_c             (float) Charging energy of the transmon (Hz).
-    f_bus           (float) frequency of the bus (Hz).
-    f_interaction   (float) interaction frequency (Hz).
-    asymmetry       (float) qubit asymmetry
-
-    sampling_rate   (float) sampling rate of the AWG (Hz)
-    return_unit     (enum: ['V', 'eps', 'f01', 'theta']) whether to return the
-                    pulse expressed in units of theta: the reference frame of
-                    the interaction, units of epsilon: detuning to the bus
-                    eps=f12-f_bus
-    """
-    # Define number of samples and time points
-
-
-    # Pulse is generated at a denser grid to allow for good interpolation
-    # N.B. Not clear why interpolation is needed at all... -MAR July 2018
-    fine_sampling_factor = 1  # 10
-    nr_samples = int(np.round((length)*sampling_rate * fine_sampling_factor))
-    rounded_length = nr_samples/(fine_sampling_factor * sampling_rate)
-    tau_step = 1/(fine_sampling_factor * sampling_rate)  # denser points
-    # tau is a virtual time/proper time
-    taus = np.arange(0, rounded_length-tau_step/2, tau_step)
-    # -tau_step/2 is to make sure final pt is excluded
-
-    # Derived parameters
-    if f_interaction is None:
-        f_interaction = f_bus + E_c
-    theta_i = np.arctan(2*J2 / (f_01_max - f_interaction))
-    # Converting angle to radians as that is used under the hood
-    theta_f = 2*np.pi*theta_f/360
-    if theta_f < theta_i:
-        raise ValueError(
-            'theta_f ({:.2f} deg) < theta_i ({:.2f} deg):'.format(
-                theta_f/(2*np.pi)*360, theta_i/(2*np.pi)*360)
-            + 'final coupling weaker than initial coupling')
-
-    # lambda_1 is scaled such that the final ("center") angle is theta_f
-    lambda_1 = (theta_f - theta_i) / (2 + 2 * lambda_3)
-
-    # Calculate the wave
-    theta_wave = np.ones(nr_samples) * theta_i
-    theta_wave += lambda_1 * (1 - np.cos(2 * np.pi * taus / rounded_length))
-    theta_wave += (lambda_1 * lambda_2 *
-                   (1 - np.cos(4 * np.pi * taus / rounded_length)))
-    theta_wave += (lambda_1 * lambda_3 *
-                   (1 - np.cos(6 * np.pi * taus / rounded_length)))
-
-    # Clip wave to [theta_i, pi] to avoid poles in the wave expressed in freq
-    theta_wave_clipped = np.clip(theta_wave, theta_i, np.pi-.01)
-    if not np.array_equal(theta_wave, theta_wave_clipped):
-        logging.warning(
-            'Martinis flux wave form has been clipped to [{}, 180 deg]'
-            .format(theta_i))
-
-    # Transform from proper time to real time
-    t = np.array([np.trapz(np.sin(theta_wave)[:i+1], dx=1/(10*sampling_rate))
-                  for i in range(len(theta_wave))])
-
-    # Interpolate pulse at physical sampling distance
-    t_samples = np.arange(0, length, 1/sampling_rate)
-    # Scaling factor for time-axis to get correct pulse length again
-    scale = t[-1]/t_samples[-1]
-    interp_wave = scipy.interpolate.interp1d(
-        t/scale, theta_wave_clipped, bounds_error=False,
-        fill_value='extrapolate')(t_samples)
-
-    # Return in the specified units
-    if return_unit == 'theta':
-        # Theta is returned in radians here
-        return np.nan_to_num(interp_wave)
-
-    # Convert to detuning from f_interaction
-    delta_f_wave = 2 * J2 / np.tan(interp_wave)
-    if return_unit == 'eps':
-        return np.nan_to_num(delta_f_wave)
-
-    # Convert to parametrization of f_01
-    f_01_wave = delta_f_wave + f_interaction
-    if return_unit == 'f01':
-        return np.nan_to_num(f_01_wave)
-
-    # Convert to voltage
-    voltage_wave = Qubit_freq_to_dac(
-        frequency=f_01_wave,
-        f_max=f_01_max,
-        E_c=E_c,
-        dac_sweet_spot=V_offset,
-        V_per_phi0=V_per_phi0,
-        asymmetry=asymmetry,
-        branch='positive')
-
-    # why sometimes the last sample is nan is not known,
-    # but we will surely figure it out someday.
-    # (Brian and Adriaan, 14.11.2017)
-    # This may be caused by the fill_value of the interp_wave (~30 lines up)
-    # that was set to 0 instead of extrapolate. This caused
-    # the np.tan(interp_wave) to divide by zero. (MAR 10-05-2018)
-    voltage_wave = np.nan_to_num(voltage_wave)
-
-    return voltage_wave
-############################################################################
-#
-############################################################################

pycqed/measurement/waveform_control_CC/waveforms_flux.py

@@ -0,0 +1,129 @@
+'''
+    File:               waveform_flux.py
+    Author:             Filip Malinowski and Adriaan Rol
+    Purpose:            generate flux waveforms that require knowledge
+        of the qubit Hamiltonian.
+    Prerequisites:
+    Usage:
+    Bugs:
+
+For more basic waveforms see e.g., waveforms.py
+'''
+import logging
+import scipy.interpolate
+import numpy as np
+
+
+def martinis_flux_pulse(length: float,
+                        theta_i: float, theta_f: float,
+                        lambda_2: float, lambda_3: float=0, lambda_4: float=0,
+                        sampling_rate: float =2.4e9):
+    """
+    Returns the pulse specified by Martinis and Geller as θ(t) specified in
+        Phys. Rev. A 90 022307 (2014).
+    Note that θ still needs to be transformed into detuning from the
+    interaction and into AWG amplitude V(t).
+
+        θ(τ) = θ_i + Σ_{n=1}^N  λ_n*(1-cos(n*2*pi*τ/τ_p))
+
+    Args:
+        length      :   lenght of the waveform (s)
+        lambda_2    :   lambda coeffecients
+        lambda_3    :
+        lambda_3    :
+        theta_i     :   initial angle of interaction (rad).
+        theta_f     :   final angle of the interaction (rad).
+        sampling_rate : sampling rate of AWG in (Hz)
+
+    This waveform is generated in several steps
+        1. Generate a time grid, may include fine sampling.
+        2. Generate θ(τ) using eqs 15 and 16
+        3. Transform from proper time "τ" to real time "t" using interpolation
+
+    """
+    if theta_f < theta_i:
+        raise ValueError(
+            'theta_f ({:.2f} deg) < theta_i ({:.2f} deg):'.format(
+                np.rad2deg(theta_f), np.rad2deg(theta_i))
+            + 'final coupling weaker than initial coupling')
+
+    # 1. Generate a time grid, may include fine sampling.
+
+    # Pulse is generated at a denser grid to allow for good interpolation
+    # N.B. Not clear why interpolation is needed at all... -MAR July 2018
+    fine_sampling_factor = 2  # 10
+    nr_samples = int(np.round((length)*sampling_rate * fine_sampling_factor))
+    rounded_length = nr_samples/(fine_sampling_factor * sampling_rate)
+    tau_step = 1/(fine_sampling_factor * sampling_rate)  # denser points
+    # tau is a virtual time/proper time
+    taus = np.arange(0, rounded_length-tau_step/2, tau_step)
+    # -tau_step/2 is to make sure final pt is excluded
+
+    # lambda_1 is scaled such that the final ("center") angle is theta_f
+    # Determine lambda_1 using the constraint set by eq 16 from Martinis 2014
+    lambda_1 = (theta_f - theta_i) / (2) - lambda_3
+
+    # 2. Generate θ(τ) using eqs 15 and 16
+    theta_wave = np.ones(nr_samples) * theta_i
+    theta_wave += lambda_1 * (1 - np.cos(2 * np.pi * taus / rounded_length))
+    theta_wave += lambda_2 * (1 - np.cos(4 * np.pi * taus / rounded_length))
+    theta_wave += lambda_3 * (1 - np.cos(6 * np.pi * taus / rounded_length))
+    theta_wave += lambda_4 * (1 - np.cos(8 * np.pi * taus / rounded_length))
+
+    # Clip wave to [theta_i, pi] to avoid poles in the wave expressed in freq
+    theta_wave_clipped = np.clip(theta_wave, theta_i, np.pi-.01)
+    if not np.array_equal(theta_wave, theta_wave_clipped):
+        logging.warning(
+            'Martinis flux wave form has been clipped to [{}, 180 deg]'
+            .format(np.rad2deg(theta_i)))
+
+    # 3. Transform from proper time τ to real time t using interpolation
+    t = np.array([np.trapz(np.sin(theta_wave_clipped)[:i+1],
+                           dx=1/(fine_sampling_factor*sampling_rate))
+                  for i in range(len(theta_wave_clipped))])
+
+    # Interpolate pulse at physical sampling distance
+    t_samples = np.arange(0, length, 1/sampling_rate)
+    # Scaling factor for time-axis to get correct pulse length again
+    scale = t[-1]/t_samples[-1]
+    interp_wave = scipy.interpolate.interp1d(
+        t/scale, theta_wave_clipped, bounds_error=False,
+        fill_value='extrapolate')(t_samples)
+
+    # Theta is returned in radians here
+    return np.nan_to_num(interp_wave)
+
+
+def eps_to_theta(eps: float, g: float):
+    """
+    Converts ε into θ as defined in Phys. Rev. A 90 022307 (2014)
+
+        θ = arctan(Hx/Hz)
+        θ = arctan(2*g/ε)
+
+    args:
+        ε: detuning
+        g: coupling strength
+    returns:
+        θ: interaction angle (radian)
+    """
+    # Ignore divide by zero as it still gives a meaningful angle
+    with np.errstate(divide='ignore'):
+        theta = np.arctan(np.divide(2*g, eps))
+    return theta
+
+
+def theta_to_eps(theta: float, g: float):
+    """
+    Converts θ into ε as defined in Phys. Rev. A 90 022307 (2014)
+        Hz = Hx/tan(θ)
+        ε = 2g/tan(θ)
+
+    args:
+        θ: interaction angle (radian)
+        g: coupling strength
+    returns:
+        ε: detuning
+    """
+    eps = 2 * g / np.tan(theta)
+    return eps