{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\nComplex Resonator Model\n=======================\n\nThis notebook shows how to fit the parameters of a complex resonator,\nusing `lmfit.Model` and defining a custom `Model` class.\n\nFollowing Khalil et al. (https://arxiv.org/abs/1108.3117), we can model the\nforward transmission of a microwave resonator with total quality factor\n$Q$, coupling quality factor $Q_e$, and resonant frequency\n$f_0$ using:\n\n\\begin{align}S_{21}(f) = 1 - \\frac{Q Q_e^{-1}}{1+2jQ(f-f_0)/f_0}\\end{align}\n\n$S_{21}$ is thus a complex function of a real frequency.\n\nBy allowing $Q_e$ to be complex, this model can take into account\nmismatches in the input and output transmission impedances.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import matplotlib.pyplot as plt\nimport numpy as np\n\nimport lmfit" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since ``scipy.optimize`` and ``lmfit`` require real parameters, we represent\n$Q_e$ as ``Q_e_real + 1j*Q_e_imag``.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def linear_resonator(f, f_0, Q, Q_e_real, Q_e_imag):\n Q_e = Q_e_real + 1j*Q_e_imag\n return (1 - (Q * Q_e**-1 / (1 + 2j * Q * (f - f_0) / f_0)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The standard practice of defining an ``lmfit`` model is as follows:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "class ResonatorModel(lmfit.model.Model):\n __doc__ = \"resonator model\" + lmfit.models.COMMON_DOC\n\n def __init__(self, *args, **kwargs):\n # pass in the defining equation so the user doesn't have to later.\n super().__init__(linear_resonator, *args, **kwargs)\n\n self.set_param_hint('Q', min=0) # Enforce Q is positive\n\n def guess(self, data, f=None, **kwargs):\n verbose = kwargs.pop('verbose', None)\n if f is None:\n return\n argmin_s21 = np.abs(data).argmin()\n fmin = f.min()\n fmax = f.max()\n f_0_guess = f[argmin_s21] # guess that the resonance is the lowest point\n Q_min = 0.1 * (f_0_guess/(fmax-fmin)) # assume the user isn't trying to fit just a small part of a resonance curve.\n delta_f = np.diff(f) # assume f is sorted\n min_delta_f = delta_f[delta_f > 0].min()\n Q_max = f_0_guess/min_delta_f # assume data actually samples the resonance reasonably\n Q_guess = np.sqrt(Q_min*Q_max) # geometric mean, why not?\n Q_e_real_guess = Q_guess/(1-np.abs(data[argmin_s21]))\n if verbose:\n print(\"fmin=\", fmin, \"fmax=\", fmax, \"f_0_guess=\", f_0_guess)\n print(\"Qmin=\", Q_min, \"Q_max=\", Q_max, \"Q_guess=\", Q_guess, \"Q_e_real_guess=\", Q_e_real_guess)\n params = self.make_params(Q=Q_guess, Q_e_real=Q_e_real_guess, Q_e_imag=0, f_0=f_0_guess)\n params['%sQ' % self.prefix].set(min=Q_min, max=Q_max)\n params['%sf_0' % self.prefix].set(min=fmin, max=fmax)\n return lmfit.models.update_param_vals(params, self.prefix, **kwargs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's use the model to generate some fake data:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "resonator = ResonatorModel()\ntrue_params = resonator.make_params(f_0=100, Q=10000, Q_e_real=9000, Q_e_imag=-9000)\n\nf = np.linspace(99.95, 100.05, 100)\ntrue_s21 = resonator.eval(params=true_params, f=f)\nnoise_scale = 0.02\nnp.random.seed(123)\nmeasured_s21 = true_s21 + noise_scale*(np.random.randn(100) + 1j*np.random.randn(100))\n\nplt.figure()\nplt.plot(f, 20*np.log10(np.abs(measured_s21)))\nplt.ylabel('|S21| (dB)')\nplt.xlabel('MHz')\nplt.title('simulated measurement')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Try out the guess method we added:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "guess = resonator.guess(measured_s21, f=f, verbose=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And now fit the data using the guess as a starting point:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "result = resonator.fit(measured_s21, params=guess, f=f, verbose=True)\n\nprint(result.fit_report() + '\\n')\nresult.params.pretty_print()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we'll make some plots of the data and fit. Define a convenience function\nfor plotting complex quantities:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def plot_ri(data, *args, **kwargs):\n plt.plot(data.real, data.imag, *args, **kwargs)\n\n\nfit_s21 = resonator.eval(params=result.params, f=f)\nguess_s21 = resonator.eval(params=guess, f=f)\n\nplt.figure()\nplot_ri(measured_s21, '.')\nplot_ri(fit_s21, 'r.-', label='best fit')\nplot_ri(guess_s21, 'k--', label='inital fit')\nplt.legend(loc='best')\nplt.xlabel('Re(S21)')\nplt.ylabel('Im(S21)')\n\nplt.figure()\nplt.plot(f, 20*np.log10(np.abs(measured_s21)), '.')\nplt.plot(f, 20*np.log10(np.abs(fit_s21)), 'r.-', label='best fit')\nplt.plot(f, 20*np.log10(np.abs(guess_s21)), 'k--', label='initial fit')\nplt.legend(loc='best')\nplt.ylabel('|S21| (dB)')\nplt.xlabel('MHz')" ] } ], "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.7.6" } }, "nbformat": 4, "nbformat_minor": 0 }