Landau-Zener Transition.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Import "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using matplotlib backend: Qt5Agg\n"
     ]
    }
   ],
   "source": [
    "%matplotlib \n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from math import pi, sqrt\n",
    "from qutip import *\n",
    "import time\n",
    "\n",
    "from numpy import min, max\n",
    "from scipy import linspace\n",
    "import scipy.signal"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "code_folding": [
     24
    ]
   },
   "outputs": [],
   "source": [
    "class ExpParams(object):\n",
    "    \"\"\" A placehoder for parameters \"\"\"\n",
    "    pass\n",
    "\n",
    "def plot_eps_vs_time(tlist, eps_list, params):\n",
    "    \"\"\"\n",
    "    Plot energy vs time\n",
    "    \"\"\"\n",
    "    \n",
    "    H0 = - params.delta/2.0 * sigmax()\n",
    "    eigen = [(H0 - eps/2.0 * sigmaz()).eigenstates() for eps in eps_list]\n",
    "\n",
    "    upper_band = [eig[0][0] for eig in eigen]\n",
    "    lower_band = [eig[0][1] for eig in eigen]\n",
    "\n",
    "    fig, ax = plt.subplots(figsize=(12,8))\n",
    "    ax.plot(tlist, upper_band, tlist, lower_band)\n",
    "    ax.set_xlabel('Time (us)', fontsize=20)\n",
    "    ax.set_ylabel('eps (MHz)', fontsize=20)\n",
    "    ax.tick_params(axis='x', labelsize=18)\n",
    "    ax.tick_params(axis='y', labelsize=18)\n",
    "    ax.set_title('Landau-Zener transition')\n",
    "    ax.grid()\n",
    "\n",
    "# def qubit_evolve(psi0, eps_list, tlist, H_args, n_th=0):\n",
    "#     \"\"\"\n",
    "#     eps_t : Function that represents eps(t) in -eps(t)/2*sz, \n",
    "#             i.e., defines how energy difference changes\n",
    "#     eps_list : numpy array of eps list with float type\n",
    "#     Return: expectation values\n",
    "#     \"\"\"\n",
    "#     # Hamiltonian\n",
    "#     sz = sigmaz()\n",
    "#     sm = destroy(2)\n",
    "\n",
    "#     H0 = - H_args['delta']/2.0 * sigmax() # - eps0/2.0 * sigmaz()     \n",
    "\n",
    "#     # collapse operators\n",
    "#     c_op_list = []\n",
    "# #     n_th = 0.0 # zero temperature\n",
    "\n",
    "#     # relaxation\n",
    "#     rate1 = H_args['gamma1'] * (1 + n_th)\n",
    "#     # excitation\n",
    "#     rate2 = H_args['gamma1'] * n_th\n",
    "#     # dephasing \n",
    "#     rate3 = H_args['gamma2']\n",
    "    \n",
    "#     c_op_list = [sqrt(rate1) * sm, sqrt(rate2) * sm.dag(), sqrt(rate3) * sz]\n",
    "# #     eps_list = np.array([eps_t(t, H_args) for t in tlist], dtype=float) \n",
    "#     H = [H0, [-1/2.0 * sz, eps_list]]\n",
    "#     output = mesolve(H, psi0, tlist, c_op_list, [sm.dag() * sm], args=H_args)  \n",
    "\n",
    "#     return output.expect[0]\n",
    "\n",
    "def qubit_evolve(tlist, eps_list, params):\n",
    "    \"\"\"\n",
    "    eps_t : Function that represents eps(t) in -eps(t)/2*sz, \n",
    "            i.e., defines how energy difference changes\n",
    "    eps_list : numpy array of eps list with float type\n",
    "    Return: expectation values\n",
    "    \"\"\"\n",
    "    # Hamiltonian\n",
    "    sz = sigmaz()\n",
    "    sm = destroy(2)\n",
    "\n",
    "    H0 = - params.delta / 2.0 * sigmax() # - eps0/2.0 * sigmaz()     \n",
    "\n",
    "    # collapse operators\n",
    "    c_op_list = []\n",
    "#     n_th = 0.0 # zero temperature\n",
    "\n",
    "    # relaxation\n",
    "    rate1 = params.gamma1 * (1 + params.n_th)\n",
    "    # excitation\n",
    "    rate2 = params.gamma1 * params.n_th\n",
    "    # dephasing \n",
    "    rate3 = params.gamma2\n",
    "    \n",
    "    c_op_list = [sqrt(rate1) * sm, sqrt(rate2) * sm.dag(), sqrt(rate3) * sz]\n",
    "#     eps_list = np.array([eps_t(t, H_args) for t in tlist], dtype=float) \n",
    "    H = [H0, [-1/2.0 * sz, eps_list]]\n",
    "    output = mesolve(H, params.psi0, tlist, c_op_list, [sm.dag() * sm])  \n",
    "\n",
    "    return output.expect[0]\n",
    "\n",
    "# def plot_probality_time(eps_t):\n",
    "#     eps_list = [eps_t(t, H_args) for t in tlist]\n",
    "#     fig, (ax,ax1) = plt.subplots(2,1,figsize=(12,10))\n",
    "#     ax.plot(tlist, np.real(p_ex), 'b', tlist, np.real(1-p_ex), 'r')\n",
    "#     ax.set_xlabel('Time (us)', fontsize=20)\n",
    "#     ax.set_ylabel('Occupation probability', fontsize=20)\n",
    "#     ax.tick_params(axis='x', labelsize=18)\n",
    "#     ax.tick_params(axis='y', labelsize=18)\n",
    "#     ax.set_title(f'Pulse out and idle, \\n f1={f1:.3f}, df01/dt={df01df*dfdt:.1f} (MHz/us), \\n df/dt={dfdt*1e3:.3f} (mPhi0/us), delta={delta:.3f} MHz, Gamma1={gamma1:.3f}, Gamma2={gamma2:.3f}', fontsize=16)\n",
    "#     ax.legend((\"R state\", \"L state\"), loc=0, fontsize=20);\n",
    "#     ax.grid()\n",
    "\n",
    "#     ax1.plot(tlist, eps_list)\n",
    "#     ax1.set_xlabel('Time (us)', fontsize=20)\n",
    "#     ax1.set_ylabel('eps (MHz)', fontsize=20)\n",
    "#     ax1.tick_params(axis='x', labelsize=18)\n",
    "#     ax1.tick_params(axis='y', labelsize=18)\n",
    "#     ax1.legend((\"R state\", \"L state\"), loc=0);\n",
    "#     ax1.grid() "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pulse out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "time elapsed = 0.17228388786315918\n"
     ]
    }
   ],
   "source": [
    "# Linear sweep from f1 to -f1\n",
    "\n",
    "# parameters\n",
    "params =ExpParams()\n",
    "params.delta = 10   #  MHz\n",
    "params.df01df= 5000/0.015 # define energy vs flux slope MHz/Phi0\n",
    "params.dfdt = 0.01/10 # in Phi0/us\n",
    "params.f1 = -0.002\n",
    "params.dt = 0.002    # us\n",
    "params.gamma1 = 0  # relaxation rate\n",
    "params.gamma2 = 0  # dephasing  rate\n",
    "params.n_th = 0 \n",
    "params.psi0 = basis(2,0)      # initial state\n",
    "\n",
    "\n",
    "tlist = np.linspace(0,2*np.absolute(params.f1)/params.dfdt\n",
    "                    , int(np.round(2*np.absolute(params.f1)/params.dfdt/params.dt)))\n",
    "eps_t = lambda t, params : (params.dfdt * t + params.f1) * params.df01df\n",
    "eps_list = np.array([eps_t(t, params) for t in tlist], dtype=float)\n",
    "\n",
    "start_time = time.time()\n",
    "p_ex = qubit_evolve(tlist, eps_list, params)\n",
    "print('time elapsed = ' + str(time.time() - start_time))\n",
    "\n",
    "# plot\n",
    "fig, ax = plt.subplots(figsize=(12,8))\n",
    "ax.plot(tlist, np.real(p_ex), 'b', tlist, np.real(1-p_ex), 'r')\n",
    "ax.plot(tlist, 1 - np.exp(-np.pi * params.delta **2 / (2 * params.df01df*params.dfdt)) * np.ones(shape(tlist)), 'k')\n",
    "ax.set_xlabel('Time (us)', fontsize=20)\n",
    "ax.set_ylabel('Occupation probability', fontsize=20)\n",
    "ax.set_title('Landau-Zener transition', fontsize=20)\n",
    "ax.tick_params(axis='x', labelsize=18)\n",
    "ax.tick_params(axis='y', labelsize=18)\n",
    "ax.legend((\"R state\", \"L state\", \"Landau-Zener formula\"), loc=0, fontsize=20);\n",
    "ax.grid()\n",
    "\n",
    "plot_eps_vs_time(tlist, eps_list, params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.round?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fa66a101fa0>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.constants import h, hbar \n",
    "t_list = np.linspace(10,1,101)*1e-9\n",
    "P_LZ = np.exp(-2 * np.pi * (h*200e6)**2 / (hbar * (h*500e6/t_list)))\n",
    "plt.plot(t_list/1e-9, P_LZ, label='s')\n",
    "plt.xlabel('time (ns)')\n",
    "plt.ylabel('P_LZ')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pulse out and idle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "code_folding": [],
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "time elapsed = 0.1523418426513672\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "333.33333333333337"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# pulse and idle\n",
    "def eps_t(t, params):\n",
    "    if t <= params.t1:\n",
    "        f = t * params.dfdt + params.f1\n",
    "    else:\n",
    "        f = params.f2\n",
    "        \n",
    "    eps = params.df01df * f\n",
    "    \n",
    "    return eps\n",
    "\n",
    "# parameters\n",
    "params =ExpParams()\n",
    "params.delta = 20   #  MHz\n",
    "params.df01df= 5000/0.015 # define energy vs flux slope MHz/Phi0\n",
    "params.dfdt = 0.01/10 # in Phi0/us\n",
    "params.f1 = -0.002\n",
    "params.f2 = -0.000\n",
    "params.dt = 0.001    # us\n",
    "params.gamma1 = 0  # relaxation rate\n",
    "params.gamma2 = 0  # dephasing  rate\n",
    "params.n_th = 0 \n",
    "params.psi0 = basis(2,0)      # initial state\n",
    "params.idle = 0.5\n",
    "\n",
    "# compute\n",
    "params.t1 = (params.f2-params.f1)/ params.dfdt\n",
    "t1 = params.t1\n",
    "dt = params.dt\n",
    "tlist1 = np.linspace(0, t1, int(np.round(t1/dt)+1))\n",
    "tlist2 = np.linspace(0, params.idle, int(np.round(params.idle/dt))+1)\n",
    "tlist = np.concatenate((tlist1, tlist2[1:]+t1))\n",
    "eps_list = np.array([eps_t(t, params) for t in tlist], dtype=float)\n",
    "\n",
    "start_time = time.time()\n",
    "p_ex = qubit_evolve(tlist, eps_list, params)\n",
    "print('time elapsed = ' + str(time.time() - start_time))\n",
    "\n",
    "# plot\n",
    "fig, (ax,ax1) = plt.subplots(2,1,figsize=(12,10))\n",
    "ax.plot(tlist, np.real(p_ex), 'b', tlist, np.real(1-p_ex), 'r')\n",
    "ax.set_xlabel('Time (us)', fontsize=20)\n",
    "ax.set_ylabel('Occupation probability', fontsize=20)\n",
    "ax.tick_params(axis='x', labelsize=18)\n",
    "ax.tick_params(axis='y', labelsize=18)\n",
    "ax.set_title(f'Pulse out and idle, \\n f1={params.f1:.3f}, df01/dt={params.df01df*params.dfdt:.1f} (MHz/us), \\n df/dt={params.dfdt*1e3:.3f} (mPhi0/us), delta={params.delta:.3f} MHz, Gamma1={params.gamma1:.3f}, Gamma2={params.gamma2:.3f}', fontsize=16)\n",
    "ax.legend((\"R state\", \"L state\"), loc=0, fontsize=20);\n",
    "ax.grid()\n",
    "\n",
    "ax1.plot(tlist, eps_list)\n",
    "ax1.set_xlabel('Time (us)', fontsize=20)\n",
    "ax1.set_ylabel('eps (MHz)', fontsize=20)\n",
    "ax1.tick_params(axis='x', labelsize=18)\n",
    "ax1.tick_params(axis='y', labelsize=18)\n",
    "ax1.legend((\"R state\", \"L state\"), loc=0);\n",
    "ax1.grid()\n",
    "\n",
    "plot_eps_vs_time(tlist, eps_list, params)\n",
    "params.df01df * params.dfdt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "# f2 vs idle , 2D\n",
    "\n",
    "# pulse and idel\n",
    "def eps_t(t, args):\n",
    "    if t <= args['t1']:\n",
    "        f = t * args['dfdt'] + args['f1']        \n",
    "    else:\n",
    "        f = args['f2']\n",
    "        \n",
    "    eps = args['df01df'] * f\n",
    "    \n",
    "    return eps\n",
    "\n",
    "# parameters\n",
    "delta = 10   #  MHz\n",
    "df01df= 5000/0.02 # define energy vs flux slope MHz/Phi0\n",
    "dfdt = 0.01/10  # in Phi0/us\n",
    "f1 = -0.001\n",
    "dt = 0.005    # us\n",
    "gamma1 = 0    # relaxation rate\n",
    "gamma2 = 0     # dephasing  rate\n",
    "\n",
    "# compute\n",
    "psi0 = basis(2,0)      # initial state\n",
    "f2_list = np.linspace(-0.0002, 0.0002, 11) \n",
    "data = np.ones((f2_list.size, np.int(np.round(idle/dt)+1)))\n",
    "for i, f2 in enumerate(f2_list):\n",
    "    t1 = (f2-f1)/ dfdt\n",
    "    tlist1 = np.linspace(0, t1, np.round(t1/dt)+1)\n",
    "    tlist2 = np.linspace(0, idle, np.round(idle/dt)+1)\n",
    "    tlist = np.concatenate((tlist1, tlist2[1:]+t1))\n",
    "    exp_args ={'dfdt':dfdt, 'df01df':df01df, 't1':t1, 'f1':f1,'f2':f2, 't1':t1}\n",
    "    eps_list = np.array([eps_t(t, exp_args) for t in tlist], dtype=float)\n",
    "    H_args = {'delta':delta, 'gamma1':gamma1, 'gamma2':gamma2}\n",
    "    p_ex = qubit_evolve(psi0, eps_list, tlist, H_args)\n",
    "    data[i,:] = p_ex[-np.int(np.round(idle/dt))-2:-1]\n",
    "\n",
    "# plot\n",
    "idle_mat, f2_mat = np.meshgrid(tlist2, f2_list)\n",
    "fig, ax = plt.subplots(figsize=(12,8))\n",
    "ax.pcolor(idle_mat, f2_mat*1e3, data)\n",
    "ax.set_xlabel('Idle time (us)',fontsize=20)\n",
    "ax.set_ylabel(r'$\\Delta$Flux (m$\\Phi_0$)',fontsize=20)\n",
    "ax.tick_params(axis='x', labelsize=18)\n",
    "ax.tick_params(axis='y', labelsize=18)\n",
    "ax.set_title(f'Pulse out and idle, delta={delta} (MHz), f1={f1:.3f}(Phi0), df01/dt={df01df*dfdt:.1f} (MHz/us), \\n df/dt={dfdt:.3f} (Phi0/us), Gamma1={gamma1:.3f}, Gamma2={gamma2:.3f}', fontsize=16)\n",
    "\n",
    "eps_list = [eps_t(t, exp_args) for t in tlist]\n",
    "fig, ax1 = plt.subplots(figsize=(12,8))\n",
    "ax1.plot(tlist, eps_list)\n",
    "ax1.set_xlabel('Time (us)', fontsize=20)\n",
    "ax1.set_ylabel('eps (MHz)', fontsize=20)\n",
    "ax.tick_params(axis='x', labelsize=18)\n",
    "ax.tick_params(axis='y', labelsize=18)\n",
    "ax1.legend((\"R state\", \"L state\"), loc=0, fontsize=20);\n",
    "ax1.grid()\n",
    "\n",
    "plot_eps_vs_time(delta, eps_list, tlist)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Out and back"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "# def eps_t_outandback_linear(t, p):\n",
    "#     \"\"\" Ideal linear ramp, idle and ramp down\"\"\"\n",
    "#     t1 = (p.f2 - p.f1) / p.dfdt\n",
    "#     if t <= t1:\n",
    "#         f = t * p.dfdt + p.f1\n",
    "#     elif t<= t1 + p.idle:\n",
    "#         f = p.f2\n",
    "#     else:\n",
    "#         f = -(t-(t1 + p.idle)) * p.dfdt + p.f2\n",
    "        \n",
    "#     eps = p.df01df * f\n",
    "    \n",
    "#     return eps\n",
    "\n",
    "# def eps_t_exponential(t, args):\n",
    "#     \"\"\" Not very good. Exponential ramp up and linear ramp-down\"\"\"\n",
    "#     factor=0.5\n",
    "\n",
    "#     if t <= args['t1'] + args['idle']:        \n",
    "#         s = (1-np.exp(-t/factor))*(args['f2']-args['f1']) + args['f1']\n",
    "        \n",
    "#     else:\n",
    "#         t2 = args['t1'] + args['idle']\n",
    "#         s = -args['dfdt']*(t-t2) + (1-np.exp(-t2/factor))*(args['f2']-args['f1']) + args['f1']\n",
    "                        \n",
    "       \n",
    "#     return s *args['df01df'] \n",
    "\n",
    "\n",
    "def eps_list_no_filter(tlist, eps_list, filter_args):\n",
    "    return eps_list\n",
    "\n",
    "def eps_list_RLC(tlist, eps_list, params):\n",
    "    \"\"\"\n",
    "    Create nonideal eps_t.\n",
    "    Args: \n",
    "        tlist: time list\n",
    "        eps_list: ideal eps list\n",
    "    Return: nonideal_eps_list\n",
    "    \"\"\"\n",
    "\n",
    "    t = tlist\n",
    "    u = eps_list - eps_list[0]\n",
    "    L, C, factor = (params.filter_params['L']\n",
    "                    , params.filter_params['C']\n",
    "                    , params.filter_params['factor'])\n",
    "    \n",
    "    RLC = np.array([1, L*factor, C*factor])   # R,L,C\n",
    "    num, den = [0,0,1], [RLC[1]*RLC[2], RLC[0]*RLC[2],1]\n",
    "    # num, den = [0,RLC[0]*RLC[2],0], [RLC[1]*RLC[2], RLC[0]*RLC[2],1]\n",
    "\n",
    "    # z, p, k = scipy.signal.besselap(1)\n",
    "    # z, p, k = scipy.signal.buttap(1)\n",
    "    # num, den = scipy.signal.zpk2tf(z, p, k)\n",
    "\n",
    "    system = scipy.signal.lti(num, den)\n",
    "    tout, y, x = scipy.signal.lsim(system, u, t)\n",
    "\n",
    "    return np.array(y + eps_list[0], dtype=float)\n",
    "\n",
    "# def eps_t_with_extension(t, args):\n",
    "#     if t <= args['t1']:\n",
    "#         f = t * args['dfdt'] + args['f1']        \n",
    "#     elif t<= args['t1'] + args['idle']:\n",
    "#         f = args['f2']\n",
    "#     elif t<=args['t1']*2 + args['idle']:\n",
    "#         f = -(t-(args['t1']+args['idle'])) * args['dfdt'] + args['f2']\n",
    "#     else:\n",
    "#         f = args['f1']\n",
    "        \n",
    "#     eps = args['df01df'] * f\n",
    "    \n",
    "#     return eps\n",
    "\n",
    "# def get_tlist_outandback(p):\n",
    "#     \"\"\" p: ExpParams instance\"\"\"\n",
    "#     t1 = (p.f2 - p.f1) / p.dfdt\n",
    "#     idle = p.idle\n",
    "#     n1, n2, n3 = np.round(t1/p.dt), np.round(p.idle/p.dt), np.round(t1/p.dt)\n",
    "    \n",
    "#     tlist1 = np.linspace(0, t1, n1+1)\n",
    "#     idle_list_dt = np.linspace(0, idle, n2+1)\n",
    "#     tlist2 = np.linspace(0, t1, n3+1)\n",
    "#     if p.extension != 0:\n",
    "#         n4 = np.round(p.extension/p.dt)\n",
    "#         tlist3 = np.linspace(0, p.extension, n4+1)\n",
    "#         tlist = np.concatenate((tlist1, idle_list_dt[1:]+t1, tlist2[1:]+t1+idle, tlist3[1:]+2*t1+idle))\n",
    "#     else:\n",
    "#         tlist = np.concatenate((tlist1, idle_list_dt[1:]+t1, tlist2[1:]+t1+idle))\n",
    "    \n",
    "#     return tlist\n",
    "\n",
    "def get_tlist_eps_list_outandback_linear(params):\n",
    "    \"\"\" Ideal linear ramp, idle and ramp down\"\"\"\n",
    "\n",
    "    t1 = (params.f2 - params.f1) / params.dfdt\n",
    "    idle = params.idle\n",
    "    dt = params.dt\n",
    "    \n",
    "    # get tlist\n",
    "    n1, n2, n3 = np.round(t1/dt).astype(int), np.round(idle/dt).astype(int), np.round(t1/dt).astype(int)\n",
    "    \n",
    "    tlist1 = np.linspace(0, t1, n1+1)\n",
    "    tlist2 = np.linspace(0, idle, n2+1)\n",
    "    tlist3 = np.linspace(0, t1, n3+1)\n",
    "    if params.extension != 0:\n",
    "        n4 = np.round(params.extension/dt)\n",
    "        tlist4 = np.linspace(0, params.extension, n4+1)\n",
    "        tlist = np.concatenate((tlist1, tlist2[1:]+t1, tlist3[1:]+t1+idle, tlist4[1:]+2*t1+idle))\n",
    "    else:\n",
    "        tlist = np.concatenate((tlist1, tlist2[1:]+t1, tlist3[1:]+t1+idle))\n",
    "    \n",
    "    # get eps_list\n",
    "    f1 = tlist1 * params.dfdt + params.f1\n",
    "    f2 = params.f2 * np.ones(tlist2.size-1)\n",
    "    f3 = -tlist3[1:] * params.dfdt + params.f2\n",
    "    if params.extension !=0:\n",
    "        f4 = params.f1 * np.ones(tlist4.size-1)\n",
    "        f_list = np.concatenate((f1, f2, f3, f4))\n",
    "    else:\n",
    "        f_list = np.concatenate((f1, f2, f3))\n",
    "    eps_list = params.df01df * f_list\n",
    "    \n",
    "    return (tlist, eps_list)\n",
    "\n",
    "def out_idle_back(params, eps_filter):\n",
    "    \"\"\" Calculate R state probability in out-idle-back measurement\n",
    "        Args:\n",
    "            params: ExpParams instance\n",
    "            eps_list_filter: a filter function that converts ideal eps_list to \n",
    "            new eps_list\n",
    "        Return:\n",
    "            numpy array with float datatype.\n",
    "    \"\"\"\n",
    "    # compute\n",
    "    p_ex_idle_list = np.ones(params.idle_list.size)\n",
    "    for j, idle in enumerate(params.idle_list):\n",
    "        params.idle = idle\n",
    "        tlist, eps_list = get_tlist_eps_list_outandback_linear(params)        \n",
    "        eps_list_filter = eps_filter(tlist, eps_list, params) # after conversion\n",
    "        p_ex = qubit_evolve(tlist, eps_list_filter, params)\n",
    "        p_ex_idle_list[j] = p_ex[-1]\n",
    "        \n",
    "    return p_ex_idle_list, tlist, eps_list, eps_list_filter\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.rint(1.8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## a fixed idle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "code_folding": [],
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "time elapsed = 0.07844996452331543\n"
     ]
    }
   ],
   "source": [
    "# pulse out and back with with a fixed idle, 1D\n",
    "\n",
    "# parameters\n",
    "params =ExpParams()\n",
    "params.delta = 10   #  MHz\n",
    "params.df01df= 5000/0.015 # define energy vs flux slope MHz/Phi0\n",
    "params.dfdt = 0.01/1 # in Phi0/us\n",
    "params.f1 = -0.002\n",
    "params.f2 = 0.0000\n",
    "params.dt = 0.002    # us\n",
    "params.gamma1 = 0  # relaxation rate\n",
    "params.gamma2 = 0  # dephasing  rate\n",
    "params.n_th = 0 \n",
    "params.extension = 0 # us, extra time after ramp-down.\n",
    "params.idle = 2\n",
    "params.psi0 = basis(2,0)      # initial state\n",
    "\n",
    "start_time = time.time()\n",
    "tlist, eps_list = get_tlist_eps_list_outandback_linear(params)\n",
    "p_ex = qubit_evolve(tlist, eps_list, params)\n",
    "print('time elapsed = ' + str(time.time() - start_time))\n",
    "\n",
    "# plot\n",
    "fig, (ax,ax1) = plt.subplots(2,1,figsize=(12,10))\n",
    "ax.plot(tlist, np.real(p_ex), 'b', tlist, np.real(1-p_ex), 'r')\n",
    "ax.set_xlabel('Time (us)', fontsize=20)\n",
    "ax.set_ylabel('Occupation probability', fontsize=20)\n",
    "ax.tick_params(axis='x', labelsize=18)\n",
    "ax.tick_params(axis='y', labelsize=18)\n",
    "ax.set_title(f'Pulse out and back, idle={params.idle:.3f} (us) \\n'\n",
    "             + f', delta={params.delta} (MHz), f1={params.f1*1e3:.3f}(mPhi0),  f2={params.f2*1e3:.3f}f(mPhi0), df01/dt={params.df01df*params.dfdt:.3f} (MHz/us)\\n'\n",
    "             + f', df/dt={params.dfdt:.3f} (Phi0/us), Gamma1={params.gamma1:.3f}, Gamma2={params.gamma2:.3f}', fontsize=16)\n",
    "\n",
    "ax.legend((\"R state\", \"L state\"), loc=0, fontsize=20);\n",
    "ax.grid()\n",
    "\n",
    "ax1.plot(tlist, eps_list)\n",
    "ax1.set_xlabel('Time (us)', fontsize=20)\n",
    "ax1.set_ylabel('eps (MHz)', fontsize=20)\n",
    "ax1.tick_params(axis='x', labelsize=18)\n",
    "ax1.tick_params(axis='y', labelsize=18)\n",
    "ax1.legend((\"R state\", \"L state\"), loc=0, fontsize=20);\n",
    "ax1.grid()\n",
    "\n",
    "# plot_eps_vs_time(delta, eps_list, tlist)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## vs idle, 1D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "code_folding": [
     4
    ]
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tlist not uniform?\n",
      "tlist not uniform?\n",
      "tlist not uniform?\n",
      "tlist not uniform?\n",
      "tlist not uniform?\n",
      "tlist not uniform?\n",
      "tlist not uniform?\n",
      "tlist not uniform?\n",
      "tlist not uniform?\n",
      "tlist not uniform?\n",
      "tlist not uniform?\n",
      "tlist not uniform?\n",
      "tlist not uniform?\n",
      "tlist not uniform?\n",
      "tlist not uniform?\n",
      "tlist not uniform?\n",
      "tlist not uniform?\n",
      "tlist not uniform?\n",
      "tlist not uniform?\n",
      "tlist not uniform?\n",
      "tlist not uniform?\n",
      "tlist not uniform?\n",
      "tlist not uniform?\n",
      "tlist not uniform?\n",
      "tlist not uniform?\n",
      "tlist not uniform?\n",
      "tlist not uniform?\n",
      "tlist not uniform?\n",
      "tlist not uniform?\n",
      "tlist not uniform?\n",
      "tlist not uniform?\n",
      "tlist not uniform?\n",
      "tlist not uniform?\n",
      "tlist not uniform?\n",
      "tlist not uniform?\n",
      "tlist not uniform?\n",
      "tlist not uniform?\n",
      "tlist not uniform?\n",
      "tlist not uniform?\n",
      "tlist not uniform?\n",
      "tlist not uniform?\n",
      "tlist not uniform?\n",
      "tlist not uniform?\n",
      "tlist not uniform?\n",
      "tlist not uniform?\n",
      "tlist not uniform?\n",
      "tlist not uniform?\n",
      "tlist not uniform?\n",
      "tlist not uniform?\n",
      "tlist not uniform?\n"
     ]
    },
    {
     "ename": "AttributeError",
     "evalue": "'numpy.ndarray' object has no attribute 'norm'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mAttributeError\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-44-5037f79c87ea>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     50\u001b[0m \u001b[0max1\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mgrid\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     51\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 52\u001b[0;31m \u001b[0mplot_eps_vs_time\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtlist\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0meps_list_filter\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mparams\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m<ipython-input-11-06d23781de2b>\u001b[0m in \u001b[0;36mplot_eps_vs_time\u001b[0;34m(tlist, eps_list, params)\u001b[0m\n\u001b[1;32m      9\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     10\u001b[0m     \u001b[0mH0\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mparams\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdelta\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;36m2.0\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0msigmax\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 11\u001b[0;31m     \u001b[0meigen\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mH0\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0meps\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;36m2.0\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0msigmaz\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0meigenstates\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0meps\u001b[0m \u001b[0;32min\u001b[0m \u001b[0meps_list\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     12\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     13\u001b[0m     \u001b[0mupper_band\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0meig\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0meig\u001b[0m \u001b[0;32min\u001b[0m \u001b[0meigen\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-11-06d23781de2b>\u001b[0m in \u001b[0;36m<listcomp>\u001b[0;34m(.0)\u001b[0m\n\u001b[1;32m      9\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     10\u001b[0m     \u001b[0mH0\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mparams\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdelta\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;36m2.0\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0msigmax\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 11\u001b[0;31m     \u001b[0meigen\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mH0\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0meps\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;36m2.0\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0msigmaz\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0meigenstates\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0meps\u001b[0m \u001b[0;32min\u001b[0m \u001b[0meps_list\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     12\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     13\u001b[0m     \u001b[0mupper_band\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0meig\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0meig\u001b[0m \u001b[0;32min\u001b[0m \u001b[0meigen\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/.local/lib/python3.8/site-packages/qutip/qobj.py\u001b[0m in \u001b[0;36meigenstates\u001b[0;34m(self, sparse, sort, eigvals, tol, maxiter)\u001b[0m\n\u001b[1;32m   1697\u001b[0m         ekets = np.array([Qobj(vec, dims=new_dims) for vec in evecs],\n\u001b[1;32m   1698\u001b[0m                          dtype=object)\n\u001b[0;32m-> 1699\u001b[0;31m         \u001b[0mnorms\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mket\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnorm\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mket\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mekets\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1700\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0mevals\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mekets\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0mnorms\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1701\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/.local/lib/python3.8/site-packages/qutip/qobj.py\u001b[0m in \u001b[0;36m<listcomp>\u001b[0;34m(.0)\u001b[0m\n\u001b[1;32m   1697\u001b[0m         ekets = np.array([Qobj(vec, dims=new_dims) for vec in evecs],\n\u001b[1;32m   1698\u001b[0m                          dtype=object)\n\u001b[0;32m-> 1699\u001b[0;31m         \u001b[0mnorms\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mket\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnorm\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mket\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mekets\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1700\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0mevals\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mekets\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0mnorms\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1701\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mAttributeError\u001b[0m: 'numpy.ndarray' object has no attribute 'norm'"
     ]
    }
   ],
   "source": [
    "# out and back\n",
    "# vs idle , 1D\n",
    "\n",
    "# parameters\n",
    "params =ExpParams()\n",
    "params.delta = 100   #  MHz\n",
    "params.df01df= 5000/0.015 # define energy vs flux slope MHz/Phi0\n",
    "params.dfdt = 0.04/0.1 # in Phi0/us\n",
    "params.f1 = -0.001\n",
    "params.f2 = 0.0000\n",
    "params.dt = 0.002    # us\n",
    "params.gamma1 = 50  # relaxation rate\n",
    "params.gamma2 = 50 # dephasing  rate\n",
    "params.n_th = 0 \n",
    "params.extension = 0 # us, extra time after ramp-down.\n",
    "params.idle_list = np.linspace(0, 0.5, 51)\n",
    "params.psi0 = basis(2,0)      # initial state\n",
    "\n",
    "filter_params = {'L':1, 'C':10 , 'factor':0.01} # L, C, factor = 1, 10, 0.005\n",
    "params.filter_params = filter_params\n",
    "\n",
    "# choose eps_t\n",
    "# TODO add 3: eps_list_exponential}\n",
    "eps_filters= {1:eps_list_no_filter, 2:eps_list_RLC}\n",
    "eps_filter = eps_filters[1]\n",
    "\n",
    "# compute\n",
    "p_ex_idle_list, tlist, eps_list, eps_list_filter = out_idle_back(params, eps_filter)\n",
    "\n",
    "# plot\n",
    "fig, (ax,ax1) = plt.subplots(2,1,figsize=(12,10))\n",
    "ax.plot(params.idle_list, p_ex_idle_list, 'b')\n",
    "ax.set_xlabel('Idle time (us)',fontsize=20)\n",
    "ax.set_ylabel('Ocupation Probability',fontsize=20)\n",
    "ax.set_ylim((-0.05,1.05))\n",
    "ax.tick_params(axis='x', labelsize=18)\n",
    "ax.tick_params(axis='y', labelsize=18)\n",
    "ax.legend((\"R state\",), loc=0, fontsize=20);\n",
    "ax.set_title(f'Pulse out and back \\n'\n",
    "             + f', delta={params.delta} (MHz), f1={params.f1*1e3:.3f}(mPhi0),  f2={params.f2*1e3:.3f}f(mPhi0), df01/dt={params.df01df*params.dfdt:.3f} (MHz/us)\\n'\n",
    "             + f', df/dt={params.dfdt:.3f} (Phi0/us), Gamma1={params.gamma1:.3f}, Gamma2={params.gamma2:.3f}', fontsize=16)\n",
    "\n",
    "ax1.plot(tlist, eps_list)\n",
    "if eps_filter != eps_filters[1]:\n",
    "    ax1.plot(tlist, eps_list_filter)\n",
    "ax1.set_xlabel('Time (us)', fontsize=20)\n",
    "ax1.set_ylabel('eps (MHz)', fontsize=20)\n",
    "ax1.tick_params(axis='x', labelsize=18)\n",
    "ax1.tick_params(axis='y', labelsize=18)\n",
    "ax1.grid()\n",
    "\n",
    "plot_eps_vs_time(tlist, eps_list_filter, params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 146,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f36bc4d3250>]"
      ]
     },
     "execution_count": 146,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eps_list_test = [eps_t_outandback_linear(t, params) for t in tlist]\n",
    "plt.plot(tlist, eps_list_test,'bo-')\n",
    "\n",
    "tlist_t, eps_list_t = get_tlist_eps_list_outandback_linear(params)\n",
    "plt.plot(tlist_t, eps_list_t, 'r*-')\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## f2 vs idle, 2D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "code_folding": [],
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "flux=1/41\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/lib/python3/dist-packages/scipy/signal/filter_design.py:1618: BadCoefficients: Badly conditioned filter coefficients (numerator): the results may be meaningless\n",
      "  warnings.warn(\"Badly conditioned filter coefficients (numerator): the \"\n"
     ]
    },
    {
     "ename": "TypeError",
     "evalue": "len() of unsized object",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                              Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-45-15c1d51b81d6>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     38\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     39\u001b[0m     \u001b[0mparams\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mf2\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mf2\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 40\u001b[0;31m     \u001b[0mp_ex_idle_list\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtlist\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0meps_list\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0meps_list_filter\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mout_idle_back\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mparams\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0meps_filter\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     41\u001b[0m     \u001b[0mdata\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mp_ex_idle_list\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     42\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-42-8a81dcae7805>\u001b[0m in \u001b[0;36mout_idle_back\u001b[0;34m(params, eps_filter)\u001b[0m\n\u001b[1;32m    139\u001b[0m         \u001b[0mtlist\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0meps_list\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mget_tlist_eps_list_outandback_linear\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mparams\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    140\u001b[0m         \u001b[0meps_list_filter\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0meps_filter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtlist\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0meps_list\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mparams\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;31m# after conversion\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 141\u001b[0;31m         \u001b[0mp_ex\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mqubit_evolve\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtlist\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0meps_list_filter\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mparams\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    142\u001b[0m         \u001b[0mp_ex_idle_list\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mj\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mp_ex\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    143\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-11-06d23781de2b>\u001b[0m in \u001b[0;36mqubit_evolve\u001b[0;34m(tlist, eps_list, params)\u001b[0m\n\u001b[1;32m     81\u001b[0m \u001b[0;31m#     eps_list = np.array([eps_t(t, H_args) for t in tlist], dtype=float)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     82\u001b[0m     \u001b[0mH\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mH0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;36m2.0\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0msz\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0meps_list\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 83\u001b[0;31m     \u001b[0moutput\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmesolve\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mH\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mparams\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpsi0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtlist\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mc_op_list\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0msm\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdag\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0msm\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     84\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     85\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0moutput\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mexpect\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/.local/lib/python3.8/site-packages/qutip/mesolve.py\u001b[0m in \u001b[0;36mmesolve\u001b[0;34m(H, rho0, tlist, c_ops, e_ops, args, options, progress_bar, _safe_mode)\u001b[0m\n\u001b[1;32m    252\u001b[0m         \u001b[0mss\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mH\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    253\u001b[0m     \u001b[0;32melif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mH\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mlist\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mQobj\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mQobjEvo\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 254\u001b[0;31m         \u001b[0mss\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_mesolve_QobjEvo\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mH\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mc_ops\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtlist\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0moptions\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    255\u001b[0m     \u001b[0;32melif\u001b[0m \u001b[0mcallable\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mH\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    256\u001b[0m         \u001b[0mss\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_mesolve_func_td\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mH\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mc_ops\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mrho0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtlist\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0moptions\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/.local/lib/python3.8/site-packages/qutip/mesolve.py\u001b[0m in \u001b[0;36m_mesolve_QobjEvo\u001b[0;34m(H, c_ops, tlist, args, opt)\u001b[0m\n\u001b[1;32m    283\u001b[0m     \u001b[0mPrepare\u001b[0m \u001b[0mthe\u001b[0m \u001b[0msystem\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mthe\u001b[0m \u001b[0msolver\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mH\u001b[0m \u001b[0mcan\u001b[0m \u001b[0mbe\u001b[0m \u001b[0man\u001b[0m \u001b[0mQobjEvo\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    284\u001b[0m     \"\"\"\n\u001b[0;32m--> 285\u001b[0;31m     \u001b[0mH_td\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mQobjEvo\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mH\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtlist\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mtlist\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    286\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0missuper\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mH_td\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcte\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    287\u001b[0m         \u001b[0mL_td\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mliouvillian\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mH_td\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/.local/lib/python3.8/site-packages/qutip/qobjevo.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, Q_object, args, copy, tlist, state0, e_ops)\u001b[0m\n\u001b[1;32m    484\u001b[0m                 \u001b[0mQ_object\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mQ_object\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    485\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 486\u001b[0;31m         \u001b[0mop_type\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_td_format_check_single\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mQ_object\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtlist\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    487\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mops\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    488\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/.local/lib/python3.8/site-packages/qutip/qobjevo.py\u001b[0m in \u001b[0;36m_td_format_check_single\u001b[0;34m(self, Q_object, tlist)\u001b[0m\n\u001b[1;32m    591\u001b[0m                         \u001b[0;32melif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mop_k\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mndarray\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    592\u001b[0m                             \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtlist\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mndarray\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mor\u001b[0m \u001b[0;32mnot\u001b[0m\u001b[0;31m \u001b[0m\u001b[0;31m\\\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 593\u001b[0;31m                                         \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mop_k\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtlist\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    594\u001b[0m                                 \u001b[0;32mraise\u001b[0m \u001b[0mTypeError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"Time list does not match\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    595\u001b[0m                             \u001b[0mop_type\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m3\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mTypeError\u001b[0m: len() of unsized object"
     ]
    }
   ],
   "source": [
    "# out and back\n",
    "# f2 vs idle , 2D\n",
    "\n",
    "import nmmn.plots\n",
    "from IPython.display import clear_output\n",
    "\n",
    "parula=nmmn.plots.parulacmap() \n",
    "turbo=nmmn.plots.turbocmap() \n",
    "\n",
    "# parameters\n",
    "params =ExpParams()\n",
    "params.delta = 10   #  MHz\n",
    "params.df01df= 5000/0.015 # define energy vs flux slope MHz/Phi0\n",
    "params.dfdt = 0.001/0.1 # in Phi0/us\n",
    "params.f1 = -0.001\n",
    "params.dt = 0.002    # us\n",
    "params.gamma1 = 1  # relaxation rate\n",
    "params.gamma2 = 1  # dephasing  rate\n",
    "params.n_th = 0 \n",
    "params.extension = 0 # us, extra time after ramp-down.\n",
    "params.idle_list = np.linspace(0, 1, 21)\n",
    "params.f2_list = np.linspace(-0.001, 0.003, 41) \n",
    "params.psi0 = basis(2,0)      # initial state\n",
    "\n",
    "filter_params = {'L':1, 'C':10 , 'factor':0.01} # L, C, factor = 1, 10, 0.005\n",
    "params.filter_params = filter_params\n",
    "\n",
    "# choose eps_t\n",
    "# TODO add 3: eps_list_outbandback_exponential}\n",
    "eps_filters= {1:eps_list_no_filter, 2:eps_list_RLC}\n",
    "eps_filter = eps_filters[2]\n",
    "\n",
    "# compute\n",
    "data = np.ones((params.f2_list.size, params.idle_list.size))\n",
    "for i, f2 in enumerate(params.f2_list):\n",
    "    clear_output()\n",
    "    print(f'flux={i+1}/{params.f2_list.size}')\n",
    "   \n",
    "    params.f2 = f2\n",
    "    p_ex_idle_list, tlist, eps_list, eps_list_filter = out_idle_back(params, eps_filter)    \n",
    "    data[i,:] = p_ex_idle_list\n",
    "\n",
    "# plot\n",
    "idle_mat, f2_mat = np.meshgrid(params.idle_list, params.f2_list)\n",
    "fig, ax = plt.subplots(figsize=(12,8))\n",
    "color = ax.pcolormesh(idle_mat, f2_mat*1e3, data, cmap=parula)\n",
    "ax.set_xlabel('Idle time (us)',fontsize=20)\n",
    "ax.set_ylabel(r'$\\Delta$Flux (m$\\Phi_0$)',fontsize=20)\n",
    "ax.tick_params(axis='x', labelsize=18)\n",
    "ax.tick_params(axis='y', labelsize=18)\n",
    "ax.set_title(f'Pulse out and back, delta={delta} (MHz), f1={f1*1e3:.3f}(mPhi0),\\n df01/df={df01df:.3f} (MHz/Phi), df01/dt={df01df*dfdt:.3f} (MHz/us), df/dt={dfdt:.3f} (Phi0/us),\\n Gamma1={gamma1:.3f}, Gamma2={gamma2:.3f}', fontsize=16)\n",
    "\n",
    "cb = fig.colorbar(color)\n",
    "cb.set_label('|R> occupation', fontsize=20)\n",
    "cb.ax.tick_params(labelsize=18)\n",
    "\n",
    "plot_eps_vs_time(tlist, eps_list_filter, params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import min, max\n",
    "from scipy import linspace\n",
    "from scipy.signal import lti, step, impulse\n",
    "\n",
    "# making transfer function\n",
    "# example from Ogata Modern Control Engineering \n",
    "# 4th edition, International Edition page 307\n",
    "\n",
    "# num and den, can be list or numpy array type\n",
    "num = [6.3223, 18, 12.811] \n",
    "den = [1, 6, 11.3223, 18, 12.811]\n",
    "\n",
    "tf = lti(num, den)\n",
    "\n",
    "# get t = time, s = unit-step response\n",
    "t, s = step(tf)\n",
    "\n",
    "# recalculate t and s to get smooth plot\n",
    "t, s = step(tf, T = linspace(min(t), t[-1], 500))\n",
    "\n",
    "# get i = impulse\n",
    "t, i = impulse(tf, T = linspace(min(t), t[-1], 500))\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "plt.plot(t, s, t, i)\n",
    "plt.title('Transient-Response Analysis')\n",
    "plt.xlabel('Time(sec)')\n",
    "plt.ylabel('Amplitude')\n",
    "plt.hlines(1, min(t), max(t), colors='r')\n",
    "plt.hlines(0, min(t), max(t))\n",
    "plt.xlim(xmax=max(t))\n",
    "plt.legend(('Unit-Step Response', 'Unit-Impulse Response'), loc=0)\n",
    "plt.grid()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.signal import butter, lfilter, freqz\n",
    "import scipy\n",
    "def butter_lowpass(cutoff, fs, order=5):\n",
    "    nyq = 0.5 * fs\n",
    "    normal_cutoff = cutoff / nyq\n",
    "    b, a = butter(order, normal_cutoff, btype='low', analog=False)\n",
    "    return b, a\n",
    "\n",
    "def butter_lowpass_filter(data, cutoff, fs, order=5):\n",
    "    b, a = butter_lowpass(cutoff, fs, order=order)\n",
    "    y = lfilter(b, a, data)\n",
    "    return y\n",
    "\n",
    "t = np.linspace(0,1,100)\n",
    "sig = np.tanh(np.linspace(0,2,100))\n",
    "# sig = np.exp(-t**2)\n",
    "# sig = np.tanh(t/0.4)\n",
    "cutoff = 0.01\n",
    "fs = 0.1\n",
    "order =5\n",
    "# proc = butter_lowpass_filter(sig, cutoff, fs, order)\n",
    "# f = scipy.fft(sig)\n",
    "# f[25:275] = 0\n",
    "# ff = scipy.ifft(f)\n",
    "# proc = sig * np.exp(-((t-1.5)/2)**2)\n",
    "plt.plot(t, sig)\n",
    "# plt.plot(t,sig, t,np.absolute(ff))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "L, C = 2, 4\n",
    "factor = 0.01\n",
    "nonideal =  nonideal_eps_t(tlist, eps_list, L, C, factor )\n",
    "\n",
    "# plot\n",
    "fig, (ax1) = plt.subplots(1, 1, figsize=(12,8))\n",
    "\n",
    "ax1.plot(tlist, eps_list, tlist, nonideal)\n",
    "ax1.set_xlabel('t',fontsize=20)\n",
    "ax1.set_ylabel(r'Response',fontsize=20)\n",
    "ax1.tick_params(axis='x', labelsize=18)\n",
    "ax1.tick_params(axis='y', labelsize=18)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## out and back with exponential"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "def eps_t_exponential(idle, params):\n",
    "    \"\"\" Not very good. Exponential ramp up and linear ramp-down\"\"\"\n",
    "    \n",
    "    idle_list = np.arange(0, idle+params.dt, params.dt)\n",
    "    f_list = ((1 - params.exp_A * np.exp(-idle_list/params.exp_tau))-(1-params.exp_A) )*(params.f2 - params.f1)/params.exp_A + params.f1\n",
    "    \n",
    "    f0 = f_list[-1] \n",
    "    tlist_trailing = np.arange(params.dt, params.exp_tau *3 , params.dt)\n",
    "    f_trailing = np.exp(-tlist_trailing/params.exp_tau) * (f0 - params.f1) + params.f1\n",
    "    \n",
    "    f_list = np.concatenate( (f_list, f_trailing))\n",
    "                        \n",
    "    eps_list = f_list  * params.df01df\n",
    "    tlist = np.concatenate((idle_list, tlist_trailing + idle))\n",
    "    \n",
    "    return tlist, eps_list\n",
    "def eps_t_exponential_double(idle, params):\n",
    "    \"\"\" double exponential. fast rise + slow rise for ramp up. fast ramp-down\"\"\"\n",
    "    idle_list = np.arange(0, idle+params.dt, params.dt)\n",
    "    f_list = ((1 - (params.exp_A * np.exp(-idle_list/params.exp_tau)\n",
    "                    + (1-params.exp_A) * np.exp(-idle_list/params.exp_tau1))))*(params.f2 - params.f1) + params.f1\n",
    "    \n",
    "    f0 = f_list[-1] \n",
    "    tlist_trailing = np.arange(params.dt, params.exp_tau *3 , params.dt)\n",
    "    f_trailing = np.exp(-tlist_trailing/params.exp_tau) * (f0 - params.f1) + params.f1\n",
    "    \n",
    "    f_list = np.concatenate( (f_list, f_trailing))\n",
    "                        \n",
    "    eps_list = f_list  * params.df01df\n",
    "    tlist = np.concatenate((idle_list, tlist_trailing + idle))\n",
    "    \n",
    "    return tlist, eps_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "code_folding": [
     3,
     19
    ]
   },
   "outputs": [],
   "source": [
    "# out and back\n",
    "# vs idle , 1D\n",
    "\n",
    "# parameters\n",
    "params =ExpParams()\n",
    "params.delta = 10   #  MHz\n",
    "params.df01df= 5000/0.015 # define energy vs flux slope MHz/Phi0\n",
    "# params.dfdt = 0.001/0.3 # in Phi0/us\n",
    "params.f1 = -0.002\n",
    "params.f2 = 0.001\n",
    "params.dt = 0.002    # us\n",
    "params.gamma1 = 0.3  # relaxation rate\n",
    "params.gamma2 = 10  #dephasing  rate\n",
    "params.n_th = 0 \n",
    "params.extension = 0 # us, extra time after ramp-down.\n",
    "params.idle_list = np.linspace(0, 4, 51)\n",
    "params.psi0 = basis(2,0)      # initial state\n",
    "# params.exp_tau = 0.2\n",
    "# params.exp_A = 1\n",
    "params.exp_tau = 0.2\n",
    "params.exp_tau1 = 3\n",
    "params.exp_A = 0.8  #  0 <=exp_A <=1\n",
    "\n",
    "\n",
    "# compute\n",
    "p_ex_idle_list = np.ones(params.idle_list.size)\n",
    "for j, idle in enumerate(params.idle_list):\n",
    "    tlist, eps_list = eps_t_exponential_double(idle, params)\n",
    "    p_ex = qubit_evolve(tlist, eps_list, params)\n",
    "    p_ex_idle_list[j] = p_ex[-1]\n",
    "\n",
    "# plot\n",
    "fig, (ax,ax1) = plt.subplots(2,1,figsize=(12,10))\n",
    "ax.plot(params.idle_list, p_ex_idle_list, label=f'gamma2={params.gamma2}')\n",
    "ax.set_xlabel('Idle time (us)',fontsize=20)\n",
    "ax.set_ylabel('0R ocupation Probability',fontsize=20)\n",
    "ax.set_ylim((-0.05,1.05))\n",
    "ax.tick_params(axis='x', labelsize=18)\n",
    "ax.tick_params(axis='y', labelsize=18)\n",
    "ax.legend((\"R state\",), loc=0, fontsize=20);\n",
    "ax.set_title(f'Pulse out and back \\n'\n",
    "             + f', delta={params.delta} (MHz), f1={params.f1*1e3:.3f}(mPhi0),  f2={params.f2*1e3:.3f}f(mPhi0), df01/df={params.df01df:.1f} (MHz/us)\\n'\n",
    "             + f', Gamma1={params.gamma1:.3f}, Gamma2={params.gamma2:.3f},exp_tau={params.exp_tau} us', fontsize=16)\n",
    "ax.legend(fontsize=18, loc=1)\n",
    "\n",
    "ax1.plot(tlist, eps_list)\n",
    "ax1.set_xlabel('Time (us)', fontsize=20)\n",
    "ax1.set_ylabel('eps (MHz)', fontsize=20)\n",
    "ax1.tick_params(axis='x', labelsize=18)\n",
    "ax1.tick_params(axis='y', labelsize=18)\n",
    "ax1.grid()\n",
    "\n",
    "        \n",
    "# plot_eps_vs_time(tlist, eps_list, params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 131,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f4e48df8dc0>"
      ]
     },
     "execution_count": 131,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ax.legend(fontsize=14, loc=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "code_folding": [
     0,
     16
    ],
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Out-idle-back pulse\\n tau =0.05 us')"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def eps_t_exponential(idle, params):\n",
    "       \n",
    "    idle_list = np.arange(0, idle+params.dt, params.dt)\n",
    "    f_list = ((1 - params.exp_A * np.exp(-idle_list/params.exp_tau))-(1-params.exp_A) )*(params.f2 - params.f1)/params.exp_A + params.f1\n",
    "    \n",
    "    f0 = f_list[-1] \n",
    "    tlist_trailing = np.arange(params.dt, params.exp_tau *3 , params.dt)\n",
    "    f_trailing = np.exp(-tlist_trailing/params.exp_tau) * (f0 - params.f1) + params.f1\n",
    "    \n",
    "    f_list = np.concatenate( (f_list, f_trailing))\n",
    "                        \n",
    "    eps_list = f_list  * params.df01df\n",
    "    tlist = np.concatenate((idle_list, tlist_trailing + idle))\n",
    "    \n",
    "    return tlist, eps_list\n",
    "\n",
    "def eps_t_exponential_double(idle, params):\n",
    "    \"\"\" double exponential. fast rise + slow rise for ramp up. fast ramp-down\"\"\"\n",
    "    idle_list = np.arange(0, idle+params.dt, params.dt)\n",
    "    f_list = ((1 - (params.exp_A * np.exp(-idle_list/params.exp_tau)\n",
    "                    + (1-params.exp_A) * np.exp(-idle_list/params.exp_tau1))))*(params.f2 - params.f1) + params.f1\n",
    "    \n",
    "    f0 = f_list[-1] \n",
    "    tlist_trailing = np.arange(params.dt, params.exp_tau *3 , params.dt)\n",
    "    f_trailing = np.exp(-tlist_trailing/params.exp_tau) * (f0 - params.f1) + params.f1\n",
    "    \n",
    "    f_list = np.concatenate( (f_list, f_trailing))\n",
    "                        \n",
    "    eps_list = f_list  * params.df01df\n",
    "    tlist = np.concatenate((idle_list, tlist_trailing + idle))\n",
    "    \n",
    "    return tlist, eps_list\n",
    "\n",
    "params.exp_tau = 0.05\n",
    "params.exp_tau1 = 1\n",
    "params.exp_A = 1  #  0 <=exp_A <=1\n",
    "\n",
    "# a, b=eps_t_exponential(2, params)\n",
    "a, b=eps_t_exponential_double(2, params)\n",
    "fig, ax = plt.subplots(1,1, figsize=(9,6))\n",
    "ax.plot(a,b)\n",
    "ax.set_xlabel('Time (us)',fontsize=20)\n",
    "ax.set_ylabel(r'eps (MHz)',fontsize=20)\n",
    "ax.tick_params(axis='x', labelsize=18)\n",
    "ax.tick_params(axis='y', labelsize=18)\n",
    "ax.set_title(f'Out-idle-back pulse\\n tau ={params.exp_tau} us', fontsize=16)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "code_folding": [
     9
    ]
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "flux=81/81\n",
      "0.003\n"
     ]
    }
   ],
   "source": [
    "# out and back\n",
    "# f2 vs idle , 2D\n",
    "\n",
    "import nmmn.plots\n",
    "from IPython.display import clear_output\n",
    "\n",
    "parula=nmmn.plots.parulacmap() \n",
    "turbo=nmmn.plots.turbocmap() \n",
    "\n",
    "# parameters\n",
    "params =ExpParams()\n",
    "params.delta = 10   #  MHz\n",
    "params.df01df= 5000/0.015 # define energy vs flux slope MHz/Phi0\n",
    "# params.dfdt = 0.01/0.1 # in Phi0/us\n",
    "params.f1 = -0.001\n",
    "params.dt = 0.002    # us\n",
    "params.gamma1 = 1  # relaxation rate\n",
    "params.gamma2 = 10  # dephasing  rate\n",
    "params.n_th = 0 \n",
    "params.extension = 0 # us, extra time after ramp-down.\n",
    "params.idle_list = np.linspace(0, 1, 51)\n",
    "params.f2_list = np.linspace(-0.001, 0.003, 81) \n",
    "# params.f2_list = np.linspace(-0.001, 0.001, 41) \n",
    "params.psi0 = basis(2,0)      # initial state\n",
    "# params.exp_tau = 0.4\n",
    "# params.exp_A = 1\n",
    "params.exp_tau = 0.05\n",
    "params.exp_tau1 = 1\n",
    "params.exp_A = 1  #  0 <=exp_A <=1\n",
    "\n",
    "\n",
    "# compute\n",
    "data = np.ones((params.f2_list.size, params.idle_list.size))\n",
    "for i, f2 in enumerate(params.f2_list):\n",
    "    clear_output()\n",
    "    print(f'flux={i+1}/{params.f2_list.size}')\n",
    "      \n",
    "    params.f2 = f2\n",
    "    print(f2)\n",
    "    p_ex_idle_list = np.ones(params.idle_list.size)\n",
    "    for j, idle in enumerate(params.idle_list):\n",
    "        tlist, eps_list = eps_t_exponential_double(idle, params)\n",
    "        p_ex = qubit_evolve(tlist, eps_list, params)\n",
    "        p_ex_idle_list[j] = p_ex[-1]\n",
    "    \n",
    "#     p_ex_idle_list, tlist, eps_list, eps_list_filter = out_idle_back(params, eps_filter)    \n",
    "    data[i,:] = p_ex_idle_list\n",
    "\n",
    "# plot\n",
    "idle_mat, f2_mat = np.meshgrid(params.idle_list, params.f2_list)\n",
    "fig, ax = plt.subplots(figsize=(12,8))\n",
    "color = ax.pcolormesh(idle_mat, f2_mat*1e3, data, cmap=parula)\n",
    "ax.set_xlabel('Idle time (us)',fontsize=20)\n",
    "ax.set_ylabel(r'$\\Delta$Flux (m$\\Phi_0$)',fontsize=20)\n",
    "ax.tick_params(axis='x', labelsize=18)\n",
    "ax.tick_params(axis='y', labelsize=18)\n",
    "ax.set_title(f'Pulse out and back, delta={params.delta} (MHz), f1={params.f1*1e3:.3f}(mPhi0)\\n, exp_tau={params.exp_tau} (us),exp_tau1={params.exp_tau1} (us), exp_A={params.exp_A}\\n df01/df={params.df01df:.3f} (MHz/Phi),\\n Gamma1={params.gamma1:.3f}, Gamma2={params.gamma2:.3f}', fontsize=16)\n",
    "\n",
    "cb = fig.colorbar(color)\n",
    "cb.set_label('|R> occupation', fontsize=20)\n",
    "cb.ax.tick_params(labelsize=18)\n",
    "\n",
    "plot_eps_vs_time(tlist, eps_list, params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1, 'Pulse out and back, delta=10 (MHz), f1=-1.000(mPhi0), exp_tau=0.4 (us)\\n df01/df=333333.333 (MHz/Phi), (Phi0/us),\\n Gamma1=1.000, Gamma2=10.000')"
      ]
     },
     "execution_count": 82,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ax.set_title(f'Pulse out and back, delta={params.delta} (MHz), f1={params.f1*1e3:.3f}(mPhi0), exp_tau={params.exp_tau} (us)\\n df01/df={params.df01df:.3f} (MHz/Phi), (Phi0/us),\\n Gamma1={params.gamma1:.3f}, Gamma2={params.gamma2:.3f}', fontsize=16)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:matplotlib.legend:No handles with labels found to put in legend.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fe2ae5f53a0>"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t = np.linspace(0,1, 101)\n",
    "tau1 = 0.4\n",
    "tau2 = 0.01\n",
    "y = -(0.2*np.exp(-t/tau1) + 1.8*np.exp(-t/tau2))/2\n",
    "plt.plot(t, y)\n",
    "plt.grid()\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# old"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.unravel_index(1,(2,2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "options = Options()\n",
    "print(options)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Rabi oscillations\n",
    "H0 = 5000/2 * sigmaz()\n",
    "H1 = sigmax()\n",
    "\n",
    "psi0 = basis(2, 0)\n",
    "tlist = np.linspace(0,0.01, 300)\n",
    "H = [H0, [H1,'w*cos(5000*t)']]\n",
    "args = {'w':100}\n",
    "# H = [H0]\n",
    "out = mesolve(H, psi0,tlist, [], [sigmax(), sigmay(),sigmaz()], args=args)\n",
    "exp_x, exp_y, exp_z = out.expect[0], out.expect[1], out.expect[2]\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(tlist, exp_x, tlist, exp_y, tlist, exp_z)\n",
    "ax.legend(('x', 'y', 'z'))\n",
    "\n",
    "# plot in Bloch sphere\n",
    "bloch = Bloch()\n",
    "bloch.add_points([exp_x, exp_y, exp_z], meth='l')\n",
    "bloch.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# into rotating frame\n",
    "H_R = lambda t: (-1j*H0*t).expm()\n",
    "\n",
    "# Rabi oscillations\n",
    "H0 = 5000/2 * sigmaz()\n",
    "H1 = sigmax()\n",
    "\n",
    "psi0 = basis(2, 0)\n",
    "tlist = np.linspace(0,0.01, 300)\n",
    "H = [H0, [H1,'w*cos(5000*t)']]\n",
    "args = {'w':100}\n",
    "# H = [H0]\n",
    "out = mesolve(H, psi0,tlist, [], [], args=args)\n",
    "\n",
    "psi= [H_R(t)*out.states*psi0 for t in tlist]\n",
    "\n",
    "exp_x, exp_y, exp_z = expect(sigmax(), psi), expect(sigmay(), psi),expect(sigmaz(), psi)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(tlist, exp_x, tlist, exp_y, tlist, exp_z)\n",
    "ax.legend(('x', 'y', 'z'))\n",
    "\n",
    "# plot in Bloch sphere\n",
    "bloch = Bloch()\n",
    "bloch.add_points([exp_x, exp_y, exp_z], meth='l')\n",
    "bloch.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hi Ganesh\n",
    "\n",
    "If you calculate the operator for unitary evolution you will get a 2x2 qobj propagator that you can multiply your initial state (ket) to get the final state:\n",
    "\n",
    "T = ...\n",
    "U = propagator(H, T, [], H_args={..})\n",
    "psi0 = basis(2.0)\n",
    "psi_T = U * psi0\n",
    "\n",
    "If you include collapse operators you will get a 4x4 qobj propagator that describes how the density matrix in vector form is transform. So you need to convert your initial state into a density matrix, like this: rho0 = ket2dm(psi0), and then convert it to vector form. We don't have any pretty functions for this yet so its a bit of a trick, which goes something like:\n",
    "\n",
    "U = propagator(H, T, c_ops, H_args={...})\n",
    "rho0 = ket2dm(psi0)\n",
    "rho_vec = mat2vec(rho0.full())\n",
    "rho_T = Qobj(vec2mat(U.data * rho_vec))\n",
    "\n",
    "You could also do the same with with Qobj instead of the underlaying sparse matrix data by doing something like:\n",
    "\n",
    "def rho_vector_form(rho):\n",
    "    return Qobj(mat2vec(rho.full()), dims=[[rho.dims[0]] * 2, [1]])\n",
    "def rho_matrix_form(rho_vec):\n",
    "    return Qobj(vec2mat(rho_vec.full()))\n",
    "\n",
    "rho_T_vec = U * rho_vector_form(rho0)\n",
    "rho_T = rho_matrix_form(rho_T_vec)\n",
    "\n",
    "\n",
    "Hope this helps,\n",
    "Rob"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "406.667px",
    "left": "454px",
    "top": "240px",
    "width": "231.188px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}