Merge 6e9f0b6 into f0b3b85

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
@@ -329,12 +337,10 @@ def mod_square_VSM(amp_G, amp_D, length, f_modulation,
 
 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'):
+                        f_q0: float,
+                        f_q1: float,
+                        anharmonicity_q0: float,
+                        sampling_rate: float =1e9):
     """
     Returns the pulse specified by Martinis and Geller
     Phys. Rev. A 90 022307 (2014).
@@ -344,116 +350,100 @@ def martinis_flux_pulse(length: float, lambda_2: float, lambda_3: float,
     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)
+    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
+    theta_f             (float) final angle of the interaction in degrees.
+                        Determines the Voltage for the center of the waveform.
+
+    f_q0                (float) frequency of the fluxed qubit (Hz).
+    f_q1                (float) frequency of the lower, not fluxed qubit (Hz)
+    anharmonicity_q0    (float) anharmonicity of the fluxed qubit (Hz)
+                        (negative for conventional transmon)
+    J2                  (float) coupling between 11-02 (Hz),
+                        approx sqrt(2) J1 (the 10-01 coupling).
+
+    sampling_rate       (float) sampling rate of the AWG (Hz)
+
     """
-    # 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
-############################################################################
-#
-############################################################################
+    raise NotImplementedError("This class is deprecated")
+#     # Define number of samples and time points
+#     logging.warning("Deprecated, use waveforms_flux.martinis_flux_pulse")
+
+#     # 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
+#     f_initial = f_q0
+#     f_interaction = f_q1 - anharmonicity_q0
+#     detuning_initial = f_q0 + anharmonicity_q0 - f_q1
+#     theta_i = np.arctan(2*J2 / detuning_initial)
+
+#     # 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_clipped)[:i+1], dx=1/(10*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)
+
+#     # 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 between f_20 and f_11
+#     # It is equal to detuning between f_11 and interaction point
+#     delta_f_wave = 2 * J2 / np.tan(interp_wave)
+
+#     # Convert to parametrization of f_01
+#     f_01_wave = delta_f_wave + f_interaction
+
+#     return np.nan_to_num(f_01_wave)
+#     # 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)
+# ############################################################################
+# #
+# ############################################################################

pycqed/measurement/waveform_control_CC/waveforms_flux.py

@@ -0,0 +1,128 @@
+'''
+    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*t/t_p))/2
+
+    Args:
+        length      :   lenght of the waveform (s)
+        lambda_2    :   lambda coeffecients
+        lambda_3    :
+        lambda_3    :
+        theta_f     :   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))
+
+    # 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/(10*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)
+        ε: detuning
+    returns:
+        g: coupling strength
+    """
+    eps = 2 * g / np.tan(theta)
+    return eps