.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/calibrate_tsqrtP.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_calibrate_tsqrtP.py: Calculate actual beta as a function of pulse length =================================================== Assuming the data is a series of pulse lengths captured on the GDS oscilloscope (acquired using `FLInst/examples/run_pulse_calibration.py`), here each pulse length capture is frequency filtered and the absolute is integrated and converted to yield :math:`t_{pulse}\sqrt{P_{peak}}`. To convert to rms power this value is divided by :math:`\sqrt{2}` to give the final beta where :math:`\beta = \frac{1}{\sqrt{2}} \int \sqrt{P_{rms}(t)} dt`. From these integrals the programmed pulse length is plotted as a function of the calculated :math:`\beta` and a fit is applied returning the fitting coefficients for the nonlinear and linear regime of the data. The user then copies these coefficients into FLInst/spincore_pp/pulse_length_conv.py in the appropriate section to obtain calibrated pulse lengths for the given amplitude. .. GENERATED FROM PYTHON SOURCE LINES 16-329 .. rst-class:: sphx-glr-horizontal * .. image-sg:: /auto_examples/images/sphx_glr_calibrate_tsqrtP_001.png :alt: amplitude = 0.05 pulse length = 1.15 μs :srcset: /auto_examples/images/sphx_glr_calibrate_tsqrtP_001.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/images/sphx_glr_calibrate_tsqrtP_002.png :alt: amplitude = 0.05 pulse length = 294.55 μs :srcset: /auto_examples/images/sphx_glr_calibrate_tsqrtP_002.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/images/sphx_glr_calibrate_tsqrtP_003.png :alt: amplitude = 0.05 pulse length = 587.95 μs :srcset: /auto_examples/images/sphx_glr_calibrate_tsqrtP_003.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/images/sphx_glr_calibrate_tsqrtP_004.png :alt: amplitude = 0.05 Frequency Domain :srcset: /auto_examples/images/sphx_glr_calibrate_tsqrtP_004.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/images/sphx_glr_calibrate_tsqrtP_005.png :alt: collect filtered analytic :srcset: /auto_examples/images/sphx_glr_calibrate_tsqrtP_005.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/images/sphx_glr_calibrate_tsqrtP_006.png :alt: Measured $\beta$ vs A * $t_{pulse}$ :srcset: /auto_examples/images/sphx_glr_calibrate_tsqrtP_006.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/images/sphx_glr_calibrate_tsqrtP_007.png :alt: Amplitude*$t_{pulse}$ vs $\beta$ :srcset: /auto_examples/images/sphx_glr_calibrate_tsqrtP_007.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/images/sphx_glr_calibrate_tsqrtP_008.png :alt: Amplitude*$t_{pulse}$ vs $\beta$, zoomed :srcset: /auto_examples/images/sphx_glr_calibrate_tsqrtP_008.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/images/sphx_glr_calibrate_tsqrtP_009.png :alt: amplitude = 0.10 pulse length = 0.58 μs :srcset: /auto_examples/images/sphx_glr_calibrate_tsqrtP_009.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/images/sphx_glr_calibrate_tsqrtP_010.png :alt: amplitude = 0.10 pulse length = 184.02 μs :srcset: /auto_examples/images/sphx_glr_calibrate_tsqrtP_010.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/images/sphx_glr_calibrate_tsqrtP_011.png :alt: amplitude = 0.10 pulse length = 367.46 μs :srcset: /auto_examples/images/sphx_glr_calibrate_tsqrtP_011.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/images/sphx_glr_calibrate_tsqrtP_012.png :alt: amplitude = 0.10 Frequency Domain :srcset: /auto_examples/images/sphx_glr_calibrate_tsqrtP_012.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/images/sphx_glr_calibrate_tsqrtP_013.png :alt: amplitude = 0.20 pulse length = 0.29 μs :srcset: /auto_examples/images/sphx_glr_calibrate_tsqrtP_013.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/images/sphx_glr_calibrate_tsqrtP_014.png :alt: amplitude = 0.20 pulse length = 78.89 μs :srcset: /auto_examples/images/sphx_glr_calibrate_tsqrtP_014.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/images/sphx_glr_calibrate_tsqrtP_015.png :alt: amplitude = 0.20 pulse length = 157.49 μs :srcset: /auto_examples/images/sphx_glr_calibrate_tsqrtP_015.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/images/sphx_glr_calibrate_tsqrtP_016.png :alt: amplitude = 0.20 Frequency Domain :srcset: /auto_examples/images/sphx_glr_calibrate_tsqrtP_016.png :class: sphx-glr-multi-img .. rst-class:: sphx-glr-script-out .. code-block:: none ---------- logging output to /home/jmfranck/pyspecdata.0.log ---------- --> simple_functions.py(126):root find_apparent_anal_freq 2024-10-26 15:51:02,441 INFO: Aliasing occurred, but we can still find that frequency! /home/jmfranck/base/lib/python3.11/site-packages/matplotlib/cbook.py:1699: ComplexWarning: Casting complex values to real discards the imaginary part return math.isfinite(val) /home/jmfranck/base/lib/python3.11/site-packages/matplotlib/cbook.py:1345: ComplexWarning: Casting complex values to real discards the imaginary part return np.asarray(x, float) **************** Coefficients for 0.05 **************** Non-linear regime coefficients: [ 3.75442878e+02 2.00599762e+06 -8.17658362e+10 -6.14612522e+15 -2.47655751e+20 -5.88569415e+24 -8.65005568e+28 -7.95171001e+32 -4.45091091e+36 -1.38696656e+40 -1.84433605e+43] Linear regime coefficients: [2.31318373e+00 2.49223410e+06] **************** Coefficients for 0.1 **************** Non-linear regime coefficients: [ 2.93528215e+02 9.59757416e+05 -1.32947207e+10 -5.17620067e+14 -1.05255134e+19 -1.27041564e+23 -9.58050939e+26 -4.56587780e+30 -1.33800953e+34 -2.20250690e+37 -1.55941251e+40] Linear regime coefficients: [4.38010537e+00 1.06971597e+06] **************** Coefficients for 0.2 **************** Non-linear regime coefficients: [ 1.62924179e+02 4.91729568e+05 -1.99215893e+09 -6.38076923e+13 -1.05046642e+18 -1.04830747e+22 -6.88084079e+25 -3.02459439e+29 -8.58707022e+32 -1.41598700e+36 -1.02210088e+39] Linear regime coefficients: [-7.43656432e+00 5.48421364e+05] 1: amplitude = 0.05 pulse length = 1.15 μs |||μs 2: amplitude = 0.05 pulse length = 294.55 μs |||μs 3: amplitude = 0.05 pulse length = 587.95 μs |||μs 4: amplitude = 0.05 Frequency Domain |||('MHz', 'μs') 5: collect filtered analytic |||μs 6: Measured $\beta$ vs A * $t_{pulse}$ |||μs 7: Amplitude*$t_{pulse}$ vs $\beta$ |||μs√W 8: Amplitude*$t_{pulse}$ vs $\beta$, zoomed |||μs√W 9: amplitude = 0.10 pulse length = 0.58 μs |||μs 10: amplitude = 0.10 pulse length = 184.02 μs |||μs 11: amplitude = 0.10 pulse length = 367.46 μs |||μs 12: amplitude = 0.10 Frequency Domain |||('MHz', 'μs') 13: amplitude = 0.20 pulse length = 0.29 μs |||μs 14: amplitude = 0.20 pulse length = 78.89 μs |||μs 15: amplitude = 0.20 pulse length = 157.49 μs |||μs 16: amplitude = 0.20 Frequency Domain |||('MHz', 'μs') | .. code-block:: Python import pyspecdata as psd from pyspecProcScripts import find_apparent_anal_freq import matplotlib.pyplot as plt import numpy as np import logging psd.init_logging() V_atten_ratio = 102.2 # attenutation ratio skip_plots = 45 # diagnostic -- set this to None, and there will be no # diagnostic plots HH_width = 1e6 typical_180 = 40e-6 # typical beta for a 180 -- it's really important to # get pulses in this regime correct # Note: the linear threshold seems to vary from one amplitude to the next with psd.figlist_var() as fl: for filename, nodename, linear_threshold in [ ( "240819_test_amp0p05_calib_pulse_calib.h5", "pulse_calib_3", 150e-6, ), ( "240819_amp0p1_calib_pulse_calib.h5", "pulse_calib_1", 270e-6, ), ( "240819_amp0p2_calib_repeat_pulse_calib.h5", "pulse_calib_8", 310e-6, ), ]: d = psd.find_file( filename, expno=nodename, exp_type="ODNP_NMR_comp/test_equipment" ) assert ( d.get_prop("postproc_type") == "GDS_capture_v1" ), "The wrong postproc_type was set so you most likely used the wrong script for acquisition" amplitude = d.get_prop("acq_params")["amplitude"] fl.basename = f"amplitude = {amplitude:.2f}" # {{{ ensure units are set if d.get_units("t") is None: logging.info( "Units for your t axis weren't set to anything so I am setting them to s" ) d.set_units("t", "s") if d.get_units("t_pulse") is None: logging.info( "Units for your t_pulse axis weren't set to anything so I am setting them to s" ) d.set_units("t_pulse", "s") # }}} d *= V_atten_ratio d /= np.sqrt(50) # V/sqrt(R) = sqrt(P_amp) # {{{ functions that streamline plotting the desired pulse length # datasets def switch_to_plot(d, j): thislen = d["t_pulse"][j] / 1e-6 fl.next(f"pulse length = {thislen:.2f} μs") def indiv_plots(d, thislabel, thiscolor): if skip_plots is None: return for j in range(len(d["t_pulse"])): if j % skip_plots == 0: switch_to_plot(d, j) fl.plot( d["t_pulse", j], alpha=0.2, color=thiscolor, label=thislabel, ) plt.ylabel(r"$\sqrt{P}$ / $\sqrt{\mathrm{W}}$") # }}} # {{{ data is already analytic, and downsampled to below 24 MHz indiv_plots(abs(d), "abs(analytic)", "orange") d, nu_a, _ = find_apparent_anal_freq(d) # find frequency of signal d.ft("t") # {{{ Diagnostic to ensure the frequency properly identified fl.next("Frequency Domain") fl.plot(d) plt.text( x=0.25, y=0.75, s=rf"$\nu_a={nu_a/1e6:0.2f}$ MHz", transform=plt.gca().transAxes, ) assert (0 > nu_a + 0.5 * HH_width) or ( 0 < nu_a - 0.5 * HH_width ), "unfortunately the region I want to filter includes DC -- this is probably not good, and you should pick a different timescale for your scope so this doesn't happen" # }}} # {{{ apply HH frequency filter d["t" : (0, nu_a - 0.5 * HH_width)] *= 0 d["t" : (nu_a + 0.5 * HH_width, None)] *= 0 # }}} d.ift("t") indiv_plots(abs(d), "filtered analytic", "red") # }}} # {{{ set up shape of data to drop the calculated beta values into beta_v_t = d.shape.pop("t").alloc(dtype=np.float64) beta_v_t.copy_axes(d).set_units(r"s√W") # }}} for j in range(len(d["t_pulse"])): s = d["t_pulse", j] int_range = abs(s).contiguous(lambda x: x > 0.05 * s.max())[0] # slightly expand int range to include rising edges int_range[0] -= 5e-6 int_range[-1] += 5e-6 # {{{ plot the integration range of all pulses prior to integrating. # Serves as diagnostic to ensure the beta is consistently # increasing. fl.push_marker() fl.basename = None fl.next("collect filtered analytic") fl.plot(abs(s["t":int_range]), alpha=0.3) fl.pop_marker() # }}} beta_v_t["t_pulse", j] = abs(s["t":int_range]).integrate( "t" ).data.item() / np.sqrt( 2 ) # t * sqrt(P_rms) # {{{ Can't use indiv_plots because we've already indexed the t_pulse # out and we also want to plot the calculated beta on top if skip_plots is not None and j % skip_plots == 0: switch_to_plot(d, j) fl.plot( abs(s["t":int_range]), color="black", label="integrated slice", ) plt.ylabel(r"$\sqrt{P_{pulse}}$ / $\sqrt{\mathrm{W}}$") plt.text( np.mean(int_range) / 1e-6, 0.25, r"$t_{pulse} \sqrt{P_{tx}} = %f \mathrm{μs} \sqrt{\mathrm{W}}$" % (beta_v_t["t_pulse", j].item() / 1e-6), ha="center", ) # }}} # {{{ show what we observe -- how does β vary with the programmed pulse length fl.basename = None # reset so all amplitudes are on same plots beta_v_t.set_plot_color_next() # cycle color list to next color thiscolor = beta_v_t.get_plot_color() # set thiscolor to the above fl.next(r"Measured $\beta$ vs A * $t_{pulse}$") beta_v_t.rename("t_pulse", "$A t_{pulse}$").name(r"$\beta$").set_units( "μs√W" ) beta_v_t[ "$A t_{pulse}$" ] *= amplitude # we only t_pulse to be multiplied by amplitude for plotting purposes only! fl.plot( (beta_v_t / 1e-6), color=thiscolor, label="amplitude = %f" % amplitude, ) plt.axhline( # the linear threshold is the threshold above which beta is linear y=linear_threshold / 1e-6, color=thiscolor, alpha=0.8, label=f"linear threshold for amp={amplitude}", ) psd.gridandtick(plt.gca()) beta_v_t[ "$A t_{pulse}$" ] /= amplitude # need to divide back out for the determination of the coefficients below beta_v_t.rename("$A t_{pulse}$", "t_pulse") # }}} # {{{ Identify captures that don't increase in beta - don't use decreasing_idx = np.nonzero(~(np.diff(beta_v_t.data) > 0))[0] if ( len(decreasing_idx) > 0 ): # beta doesn't always increase with increasing pulse length fl.plot( beta_v_t["t_pulse", : decreasing_idx[-1] + 1], "x", color="r", label="can't use these", ) # throw out betas following a faulty pulse beta_v_t = beta_v_t["t_pulse", decreasing_idx[-1] + 1 :] # }}} # {{{ flip data so beta is on x axis now and t_pulse is the y axis t_v_beta = ( beta_v_t.shape.alloc(dtype=np.float64) .rename("t_pulse", r"$\beta$") .setaxis(r"$\beta$", beta_v_t.data) ) t_v_beta.data[:] = ( beta_v_t["t_pulse"].copy() / 1e-6 ) # because our ppg wants μs t_v_beta.set_units("μs").set_units(r"$\beta$", "s√W") # use as temp for ultimate coeff # }}} # {{{ Determine linear and nonliear coefficients c_nonlinear = t_v_beta[ r"$\beta$":(None, linear_threshold) ].C # make copy to allow the next # line not to cause the data to actually shift c_nonlinear[ r"$\beta$" ] -= linear_threshold # Taylor expand around the linear threshold rather than 0 c_nonlinear = c_nonlinear.polyfit(r"$\beta$", order=10) c_linear = t_v_beta[r"$\beta$":(linear_threshold, None)].polyfit( r"$\beta$", order=1 ) print( f"\n**************** Coefficients for {amplitude} ****************\n" ) print("Non-linear regime coefficients:\n", c_nonlinear) print("Linear regime coefficients:\n", c_linear) # }}} def prog_plen(desired): """If the desired beta is above the linear threshold, the linear fit is applied to return the pulse length required to obtain the desired beta. If the desired beta is below the linear threshold, the nonlinear fit is applied to return the pulse length required to obtain the desired beta. If a series of desired betas is fed, both fits are applied appropriately to return the pulse lengths required to obtain each desired beta in the series. """ def zonefit(desired): if desired > linear_threshold: return np.polyval(c_linear[::-1], desired) else: return np.polyval( c_nonlinear[::-1], desired - linear_threshold ) ret_val = np.vectorize(zonefit)(desired) if ret_val.size > 1: return ret_val else: return ret_val.item() fl.next(r"Amplitude*$t_{pulse}$ vs $\beta$", legend=True) t_v_beta.name("A * $t_{pulse}$").set_units("μs") fl.plot( t_v_beta * amplitude, ".", alpha=0.5, color=thiscolor, label=f"data for {amplitude}", ) fl.next(r"Amplitude*$t_{pulse}$ vs $\beta$, zoomed") fl.plot( t_v_beta[r"$\beta$":(None, typical_180)] * amplitude, ".", alpha=0.5, color=thiscolor, label=f"data for {amplitude}", ) # {{{ we extrapolate past the edges of the data using the linear threshold # to show how the linear is poorly behaved for smaller beta values for_extrap = ( psd.nddata( np.linspace(0.5e-6, t_v_beta[r"$\beta$"].max() + 10e-6, 500), r"$\beta$", ) .set_units("μs") .set_units(r"$\beta$", "s√W") ) fl.plot( for_extrap.eval_poly(c_linear, r"$\beta$")[ r"$\beta$":(None, typical_180) ] * amplitude, "--", alpha=0.25, color=thiscolor, label="linear", ) fl.next(r"Amplitude*$t_{pulse}$ vs $\beta$") fl.plot( for_extrap.eval_poly(c_linear, r"$\beta$") * amplitude, "--", alpha=0.25, color=thiscolor, label="linear", ) # }}} full_fit = for_extrap.fromaxis(r"$\beta$").run(prog_plen) fl.plot(full_fit * amplitude, alpha=0.5, color=thiscolor, label="fit") plt.axvline( x=linear_threshold / 1e-6, # units of μs alpha=0.4, color=thiscolor, ls="dashdot", label="linear threshold", ) fl.next(r"Amplitude*$t_{pulse}$ vs $\beta$, zoomed") fl.plot( full_fit[r"$\beta$":(None, typical_180)] * amplitude, alpha=0.5, color=thiscolor, label="fit", ) # {{{ Add grids to the fit plots for j in ( r"Amplitude*$t_{pulse}$ vs $\beta$", r"Amplitude*$t_{pulse}$ vs $\beta$, zoomed", ): fl.next(j) psd.gridandtick(plt.gca()) # }}} .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 44.446 seconds) .. _sphx_glr_download_auto_examples_calibrate_tsqrtP.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: calibrate_tsqrtP.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: calibrate_tsqrtP.py ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_