.. note:: :class: sphx-glr-download-link-note Click :ref:`here ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_examples_example_complex_resonator_model.py: Complex Resonator Model ======================= This notebook shows how to fit the parameters of a complex resonator, using `lmfit.Model` and defining a custom `Model` class. Following Khalil et al. (https://arxiv.org/abs/1108.3117), we can model the forward transmission of a microwave resonator with total quality factor :math:`Q`, coupling quality factor :math:`Q_e`, and resonant frequency :math:`f_0` using: .. math:: S_{21}(f) = 1 - \frac{Q Q_e^{-1}}{1+2jQ(f-f_0)/f_0} :math:`S_{21}` is thus a complex function of a real frequency. By allowing :math:`Q_e` to be complex, this model can take into account mismatches in the input and output transmission impedances. .. code-block:: default import matplotlib.pyplot as plt import numpy as np import lmfit Since ``scipy.optimize`` and ``lmfit`` require real parameters, we represent :math:`Q_e` as ``Q_e_real + 1j*Q_e_imag``. .. code-block:: default def linear_resonator(f, f_0, Q, Q_e_real, Q_e_imag): Q_e = Q_e_real + 1j*Q_e_imag return (1 - (Q * Q_e**-1 / (1 + 2j * Q * (f - f_0) / f_0))) The standard practice of defining an ``lmfit`` model is as follows: .. code-block:: default class ResonatorModel(lmfit.model.Model): __doc__ = "resonator model" + lmfit.models.COMMON_DOC def __init__(self, *args, **kwargs): # pass in the defining equation so the user doesn't have to later. super().__init__(linear_resonator, *args, **kwargs) self.set_param_hint('Q', min=0) # Enforce Q is positive def guess(self, data, f=None, **kwargs): verbose = kwargs.pop('verbose', None) if f is None: return argmin_s21 = np.abs(data).argmin() fmin = f.min() fmax = f.max() f_0_guess = f[argmin_s21] # guess that the resonance is the lowest point 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. delta_f = np.diff(f) # assume f is sorted min_delta_f = delta_f[delta_f > 0].min() Q_max = f_0_guess/min_delta_f # assume data actually samples the resonance reasonably Q_guess = np.sqrt(Q_min*Q_max) # geometric mean, why not? Q_e_real_guess = Q_guess/(1-np.abs(data[argmin_s21])) if verbose: print("fmin=", fmin, "fmax=", fmax, "f_0_guess=", f_0_guess) print("Qmin=", Q_min, "Q_max=", Q_max, "Q_guess=", Q_guess, "Q_e_real_guess=", Q_e_real_guess) params = self.make_params(Q=Q_guess, Q_e_real=Q_e_real_guess, Q_e_imag=0, f_0=f_0_guess) params['%sQ' % self.prefix].set(min=Q_min, max=Q_max) params['%sf_0' % self.prefix].set(min=fmin, max=fmax) return lmfit.models.update_param_vals(params, self.prefix, **kwargs) Now let's use the model to generate some fake data: .. code-block:: default resonator = ResonatorModel() true_params = resonator.make_params(f_0=100, Q=10000, Q_e_real=9000, Q_e_imag=-9000) f = np.linspace(99.95, 100.05, 100) true_s21 = resonator.eval(params=true_params, f=f) noise_scale = 0.02 np.random.seed(123) measured_s21 = true_s21 + noise_scale*(np.random.randn(100) + 1j*np.random.randn(100)) plt.figure() plt.plot(f, 20*np.log10(np.abs(measured_s21))) plt.ylabel('|S21| (dB)') plt.xlabel('MHz') plt.title('simulated measurement') .. image:: /examples/images/sphx_glr_example_complex_resonator_model_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out Out: .. code-block:: none Text(0.5, 1.0, 'simulated measurement') Try out the guess method we added: .. code-block:: default guess = resonator.guess(measured_s21, f=f, verbose=True) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none fmin= 99.95 fmax= 100.05 f_0_guess= 100.00353535353536 Qmin= 100.00353535354105 Q_max= 99003.50000055433 Q_guess= 3146.537781821432 Q_e_real_guess= 5082.2474265369565 And now fit the data using the guess as a starting point: .. code-block:: default result = resonator.fit(measured_s21, params=guess, f=f, verbose=True) print(result.fit_report() + '\n') result.params.pretty_print() .. rst-class:: sphx-glr-script-out Out: .. code-block:: none [[Model]] Model(linear_resonator) [[Fit Statistics]] # fitting method = leastsq # function evals = 36 # data points = 200 # variables = 4 chi-square = 0.08533642 reduced chi-square = 4.3539e-04 Akaike info crit = -1543.89425 Bayesian info crit = -1530.70099 [[Variables]] f_0: 100.000096 +/- 7.0378e-05 (0.00%) (init = 100.0035) Q: 10059.4926 +/- 142.294761 (1.41%) (init = 3146.538) Q_e_real: 9180.61935 +/- 133.776862 (1.46%) (init = 5082.247) Q_e_imag: -9137.03303 +/- 133.770211 (1.46%) (init = 0) [[Correlations]] (unreported correlations are < 0.100) C(Q, Q_e_real) = 0.518 C(f_0, Q_e_imag) = 0.517 C(f_0, Q_e_real) = 0.515 C(Q, Q_e_imag) = -0.515 Name Value Min Max Stderr Vary Expr Brute_Step Q 1.006e+04 100 9.9e+04 142.3 True None None Q_e_imag -9137 -inf inf 133.8 True None None Q_e_real 9181 -inf inf 133.8 True None None f_0 100 99.95 100 7.038e-05 True None None Now we'll make some plots of the data and fit. Define a convenience function for plotting complex quantities: .. code-block:: default def plot_ri(data, *args, **kwargs): plt.plot(data.real, data.imag, *args, **kwargs) fit_s21 = resonator.eval(params=result.params, f=f) guess_s21 = resonator.eval(params=guess, f=f) plt.figure() plot_ri(measured_s21, '.') plot_ri(fit_s21, 'r.-', label='best fit') plot_ri(guess_s21, 'k--', label='inital fit') plt.legend(loc='best') plt.xlabel('Re(S21)') plt.ylabel('Im(S21)') plt.figure() plt.plot(f, 20*np.log10(np.abs(measured_s21)), '.') plt.plot(f, 20*np.log10(np.abs(fit_s21)), 'r.-', label='best fit') plt.plot(f, 20*np.log10(np.abs(guess_s21)), 'k--', label='initial fit') plt.legend(loc='best') plt.ylabel('|S21| (dB)') plt.xlabel('MHz') .. rst-class:: sphx-glr-horizontal * .. image:: /examples/images/sphx_glr_example_complex_resonator_model_002.png :class: sphx-glr-multi-img * .. image:: /examples/images/sphx_glr_example_complex_resonator_model_003.png :class: sphx-glr-multi-img .. rst-class:: sphx-glr-script-out Out: .. code-block:: none Text(0.5, 0, 'MHz') .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 0.222 seconds) .. _sphx_glr_download_examples_example_complex_resonator_model.py: .. only :: html .. container:: sphx-glr-footer :class: sphx-glr-footer-example .. container:: sphx-glr-download :download:`Download Python source code: example_complex_resonator_model.py ` .. container:: sphx-glr-download :download:`Download Jupyter notebook: example_complex_resonator_model.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_