more

docs/theory/force_calibration/diode.rst

@@ -7,20 +7,36 @@ In reality, our measurement is affected by two processes:
 1. The motion of the bead in the trap.
 2. The response of the detector to the incident light.
 
+.. image:: figures/diode_filtering.png
+  :nbattach:
+
 This second factor depends on the type of measurement device being used.
+Unfortunately, silicon has an increased transparency at the near infra-red wavelength of the trapping laser :cite:`berg2003unintended,berg2006power`.
+What this means is that light penetrates into the substrate of the diode where it also results in charge carriers.
+These charge carriers diffuse much more slowly resulting in a fraction of the signal having a slower dynamic response (i.e. a lower bandwidth).
 
-Unfortunately, silicon has an increased transparency at the near infra-red wavelength of the trapping laser.
+In other words, the PSD becomes less sensitive to changes in signal at high frequencies.
+This effect is often referred to as the parasitic filtering effect.
+In literature, it is frequently modelled up to good accuracy with a first order approximation :cite:`berg2003unintended,tolic2006calibration,berg2006power`.
 
+.. math::
 
-The contribution to the power spectrum by the bead in the optical trap is characterized by a diffusion constant and a corner frequency.
-The second contribution refers to a filtering effect where the PSD becomes less sensitive to changes in signal at high frequencies.
-This filtering effect is characterized by a constant that reflects the fraction of light that is transmitted instantaneously and a corner frequency (this is referred to as the diode frequency to be able to differentiate).
+    g(f, f_\mathrm{diode}, \alpha) = \alpha^2 + \frac{1 - \alpha ^ 2}{1 + (f / f_\mathrm{diode})^2} \tag{$-$}
 
-.. image:: figures/diode_filtering.png
-  :nbattach:
+Using this model, the parasitic filtering effect is characterized by a constant (named relaxation factor) that reflects the fraction of light that is transmitted instantaneously and a single corner frequency for the light that is not.
+In the documentation, we generally refer to this second corner frequency as the diode frequency in order to differentiate between the two.
+
+High corner frequencies
+^^^^^^^^^^^^^^^^^^^^^^^
+
+Frequently, these diode parameters are estimated from the same calibration data.
+Unfortunately, there are cases where this can cause issues when calibrating.
 
-.. warning::
+The corner frequency of the physical spectrum can be found in the results as `f_c` and depends on the laser power and bead size (smaller beads resulting in higher corner frequencies).
+The corner frequency of the filtering effect can be found in the results as `f_diode` (which stands for diode frequency) and depends on the incident intensity :cite:`berg2003unintended`.
+When these two frequencies get close, they cannot be estimated from the power spectrum reliably anymore.
+The reason for this is that the effects that these parameters have on the power spectrum becomes very similar.
+When working with small beads or at high laser powers, it is important to verify that the corner frequency `f_c` does not approach the frequency of the filtering effect `f_diode`.
 
-    For high corner frequencies, the parameter estimation procedure can become unreliable.
-    A warning sign for this is when the corner frequency `f_c` approaches or goes beyond the diode frequency `f_diode`.
-    For more information see the section on `High corner frequencies`_.
+Sometimes, the parameters that characterize this diode have been characterized independently.
+In that case, the arguments `fixed_diode` and `fixed_alpha` can be passed to :func:`~lumicks.pylake.calibrate_force()` to fix these parameters to their predetermined values.

docs/theory/force_calibration/fitting.rst

@@ -1,13 +1,15 @@
 Fitting a power spectrum
 ------------------------
 
-In the previous section, the physical origin of the power spectrum was introduced.
-However, there are some practical aspects to consider.
-So far, we have only considered the expectation value of the power spectrum.
+In the previous sections, the physical origin of the power spectrum was introduced.
+However, there are some additional practical aspects to consider.
+
+So far, we have only considered the expected value of the power spectrum.
 In reality, power spectral values follow a distribution.
 
-The real and imaginary part of the frequency spectrum are normally distributed.
+The real and imaginary part of the spectrum are normally distributed.
 As a consequence, the squared magnitude of the power spectrum is exponentially distributed.
+
 This has two consequences:
 
 - Fitting the power spectral values directly using a simple least squares fitting routine, we would
@@ -27,17 +29,32 @@ There are two ways to perform such averaging:
   spectral domain by averaging adjacent bins according to :cite:`berg2004power`. This procedure is
   referred to as *blocking*.
 
-We use the blocking method for spectral averaging, since this allows us to reject noise peaks at high
-resolution prior to averaging. Note however, that the error incurred by this blocking procedure depends
-on :math:`n_b`, the number of points per block, :math:`\Delta f`, the spectral resolution and inversely
-on the corner frequency :cite:`berg2004power`.
+Pylake uses the blocking method for spectral averaging, since this allows us to reject noise peaks
+at high resolution prior to averaging (more on this later).
+Note however, that the error incurred by this blocking procedure depends on :math:`n_b`, the number
+of points per block, :math:`\Delta f`, the spectral resolution and inversely on the corner
+frequency :cite:`berg2004power`.
 
-Setting the number of points per block too low would result in a bias from insufficient averaging
+Setting the number of points per block too low results in a bias from insufficient averaging
 :cite:`berg2004power`. Insufficient averaging would result in an overestimation of the force response
 (:math:`R_f`) and an underestimation of the distance response (:math:`R_d`). In practice, one should
 use a high number of points per block (:math:`n_b \gg 100`), unless a very low corner frequency precludes this.
 In such cases, it is preferable to increase the measurement time.
 
+Noise floor
+^^^^^^^^^^^
+
+When operating at very low powers and corner frequencies, it is important to ensure that the
+upper limit of the fitting range does not include the noise floor.
+
+Noise peaks
+^^^^^^^^^^^
+
+Optical tweezers are precision instruments.
+As such, it is not always possible to exclude all potential sources of noise.
+One of the main advantages of power spectral analysis is that noise sources which are present
+at particular frequencies can easily be excluded prior to analysis.
+
 .. _goodness_of_fit:
 
 Goodness of fit

docs/theory/force_calibration/hyco.rst

@@ -1,6 +1,12 @@
 Hydrodynamically correct model
 ------------------------------
 
+.. only:: html
+
+    :nbexport:`Download this page as a Jupyter notebook <self>`
+
+.. _hydro_model_theory:
+
 While the idealized model of the bead motion is sometimes sufficiently accurate,
 it neglects inertial and hydrodynamical effects of the fluid and the beads.
 
@@ -22,12 +28,54 @@ sizes and higher trap powers the differences can be substantial.
 .. image:: figures/hydro.png
   :nbattach:
 
-The theory discussed above holds for beads far away from any surface (bulk).
-When doing experiments near the surface, the surface starts to have an effect on the friction felt by the bead.
-The hydrodynamically correct model can take this into account up to a certain distance from the surface, after which different models are more suitable.
+.. _fast_sensor_hyco:
+
+Fast sensor measurement
+^^^^^^^^^^^^^^^^^^^^^^^
+
+When fitting a power spectrum, one may ask the question, so why does the fit look good if the model is bad?
+The answer to this lies in the model that is used to capture the parasitic filtering effect.
+When the parameters of this model are estimated, what can happen is that they "hide" the mis-specification of the model.
+
+Fast detectors have the ability to respond much faster to incoming light resulting in no visible filtering effect in the frequency range we are fitting.
+This means that for a fast detector, we do not need to include such a filtering effect in our model and see the power spectrum for what it really is.
+
+We can omit this effect by passing `fast_sensor=True` to the calibration models or to :func:`~lumicks.pylake.calibrate_force()`.
+Note however, that this makes using the hydrodynamically correct model critical, as the simple model doesn't actually capture the data very well.
+The following example data acquired on a fast sensor will illustrate why::
+
+    volts = f.force2y / f.force2y.calibration[0].force_sensitivity
+
+    shared_parameters = {
+        "force_voltage_data": volts.data,
+        "bead_diameter": 4.38,
+        "temperature": 25,
+        "sample_rate": volts.sample_rate,
+        "fit_range": (1e2, 23e3),
+        "num_points_per_block": 200,
+        "excluded_ranges": ([190, 210], [13600, 14600])
+    }
+
+    plt.figure(figsize=(13, 4))
+    plt.subplot(1, 3, 1)
+    fit = lk.calibrate_force(**shared_parameters, hydrodynamically_correct=False, fast_sensor=False)
+    fit.plot()
+    plt.title(f"Simple model + Slow (kappa={fit['kappa'].value:.2f})")
+    plt.subplot(1, 3, 2)
+    fit = lk.calibrate_force(**shared_parameters, hydrodynamically_correct=False, fast_sensor=True)
+    fit.plot()
+    plt.title(f"Simple model + Fast (kappa={fit['kappa'].value:.2f})")
+    plt.subplot(1, 3, 3)
+    fit = lk.calibrate_force(**shared_parameters, hydrodynamically_correct=True, fast_sensor=True)
+    fit.plot()
+    plt.title(f"Hydrodynamically correct + Fast (kappa={fit['kappa'].value:.2f})")
+    plt.tight_layout()
+    plt.show()
+
+.. image:: figures/hydro_fast.png
 
 Mathematical background
------------------------
+^^^^^^^^^^^^^^^^^^^^^^^
 
 The following equation accounts for a frequency dependent drag :cite:`tolic2006calibration`:
 
@@ -63,94 +111,3 @@ Here :math:`f_{\nu}` is the frequency at which the penetration depth equals the
 :math:`4 \nu/(\pi d^2)` with :math:`\nu` the kinematic viscosity.
 
 This approximation is reasonable, when the bead is far from the surface.
-
-When approaching the surface, the drag experienced by the bead depends on the distance between the
-bead and the surface of the flow cell. An approximate expression for the frequency dependent drag is
-then given by :cite:`tolic2006calibration`:
-
-.. math::
-
-    \gamma(f, R/l) = \frac{\gamma_\mathrm{stokes}(f)}{1 - \frac{9}{16}\frac{R}{l}
-    \left(1 - \left((1 - i)/3\right)\sqrt{\frac{f}{f_{\nu}}} + \frac{2}{9}\frac{f}{f_{\nu}}i -
-    \frac{4}{3}(1 - e^{-(1-i)(2l-R)/\delta})\right)} \tag{$\mathrm{kg/s}$}
-
-Where :math:`\delta = R \sqrt{\frac{f_{\nu}}{f}}` represents the aforementioned penetration depth,
-:math:`R` corresponds to the bead radius and :math:`l` to the distance from the bead center to the nearest surface.
-
-While these models may look daunting, they are all available in Pylake and can be used by simply
-providing a few additional arguments to the :class:`~.PassiveCalibrationModel`. It is recommended to
-use these equations when less than 10% systematic error is desired :cite:`tolic2006calibration`.
-No general statement can be made regarding the accuracy that can be achieved with the simple Lorentzian
-model, nor the direction of the systematic error, as it depends on several physical parameters involved
-in calibration :cite:`tolic2006calibration,berg2006power`.
-
-.. note::
-
-    One thing to note is that when considering the surface in the calibration procedure, the drag
-    coefficient returned from the model corresponds to the drag coefficient extrapolated back to its
-    bulk value.
-
-Faxen's law
-^^^^^^^^^^^
-
-The hydrodynamically correct model presented in the previous section works well when the bead center
-is at least 1.5 times the radius above the surface. When moving closer than this limit, we fall back
-to a model that more accurately describes the change in drag at low frequencies, but neglects the
-frequency dependent effects.
-
-To understand why, let's introduce Faxen's approximation for drag on a sphere near a surface under
-creeping flow conditions. This model is used for lateral calibration very close to a surface
-:cite:`schaffer2007surface` and is given by the following equation:
-
-.. math::
-
-    \gamma_\mathrm{faxen}(R/l) = \frac{\gamma_0}{
-        1 - \frac{9R}{16l} + \frac{1R^3}{8l^3} - \frac{45R^4}{256l^4} - \frac{1R^5}{16l^5}
-    } \tag{$\mathrm{kg/s}$}
-
-At frequency zero, the frequency dependent model used in the previous section reproduces this model
-up to and including its second order term in :math:`R/l`. It is, however, a lower order model and the
-accuracy decreases rapidly as the distance between the bead and surface become very small.
-The figure below shows how the model predictions at frequency zero deviate strongly from the higher order model:
-
-.. image:: figures/freq_dependent_drag_zero.png
-  :nbattach:
-
-In addition, the deviation from a Lorentzian due to the frequency dependence of the drag is reduced
-upon approaching a surface :cite:`schaffer2007surface`.
-
-.. image:: figures/freq_dependence_near.png
-  :nbattach:
-
-These two aspects make using Faxen's law in combination with a Lorentzian a more suitable model for
-situations where we have to calibrate extremely close to the surface.
-
-Axial Calibration
-^^^^^^^^^^^^^^^^^
-
-For calibration in the axial direction, no hydrodynamically correct theory exists.
-
-Similarly as for the lateral component, we will fall back to a model that describes the change in
-drag at low frequencies. However, while we had a simple expression for the lateral drag as a function
-of distance, no simple closed-form equation exists for the axial dimension. Brenner et al provide an
-exact infinite series solution :cite:`brenner1961slow`. Based on this solution :cite:`schaffer2007surface`
-derived a simple equation which approximates the distance dependence of the axial drag coefficient.
-
-.. math::
-
-    \gamma_\mathrm{axial}(R/l) = \frac{\gamma_0}{
-        1.0
-        - \frac{9R}{8l}
-        + \frac{1R^3}{2l^3}
-        - \frac{57R^4}{100l^4}
-        + \frac{1R^5}{5l^5}
-        + \frac{7R^{11}}{200l^{11}}
-        - \frac{1R^{12}}{25l^{12}}
-    } \tag{$\mathrm{kg/s}$}
-
-This model deviates less than 0.1% from Brenner's exact formula for :math:`l/R >= 1.1` and less than
-0.3% over the entire range of :math:`l` :cite:`schaffer2007surface`.
-Plotting these reveals that there is a larger effect of the surface in the axial than lateral direction.
-
-.. image:: figures/drag_coefficient.png
-  :nbattach:

docs/theory/force_calibration/index.rst

@@ -8,7 +8,9 @@ Find out about force calibration
     :maxdepth: 1
 
     force_calibration
+    diode
     fitting
     passive
     hyco
+    surface
     active