#!/usr/bin/env python3
"""
optical-tweezer-v1.py -- the numbers behind "Held Without Touching" on
solvetheuniverse.com, computed from stated inputs rather than quoted.

  1. Transverse gradient force and trap stiffness for a Rayleigh
     (point-dipole) particle in a focused Gaussian beam
  2. The same calculation run on Neuman & Block's own stated experimental
     case (0.5 um polystyrene bead, water, 1064 nm, 1.2 NA), compared
     against their stated stiffness for it (a worked example, not itself
     tied to a measurement) and, separately, a real measured stiffness
     from a different trapped bead elsewhere in the same paper (Fig. 10)
  3. Why the computed value disagrees with both by about two orders of
     magnitude, and why that is expected rather than a bug (the bead's
     radius is comparable to the wavelength in the medium, the regime the
     point-dipole approximation is not supposed to cover)

KNOWN DISCREPANCY, disclosed rather than hidden: Neuman & Block's own
Eqs. (3)-(4), read directly off a rendered page image (not OCR), are
F = 2 pi alpha / (c n_m^2) grad(I0) with alpha = n_m^2 a^3 (m^2-1)/(m^2+2).
Substituting one into the other cancels every factor of n_m, leaving no
medium-index dependence at all. This script instead uses the more
commonly cited form below, linear in n_m, which matches both an
independent hand derivation from the induced-dipole force and the form
given on Wikipedia's Optical tweezers page (sourced there to Harada &
Asakura, Opt. Commun. 124, 529, 1996). Which reading is intended by this
particular paper's own printed equations is not settled here. The choice
changes every number this script prints by a factor of n_m = 1.33, not
by an order of magnitude, and changes nothing about which size regime
the point-dipole approximation is valid in.

Run with no arguments to print the numbers used in the article.
Run with --selftest to check the scaling behaviour every treatment of
this model agrees on (proportional to power, to 1/waist^4, to particle
volume, and zero for an index-matched particle) -- the parts that do not
depend on which regime the point-dipole approximation is valid in.

EVERY INPUT BELOW IS AN ASSUMPTION, chosen either to be typical of
published optical-tweezers experiments, or copied directly from Neuman &
Block, "Optical trapping," Rev. Sci. Instrum. 75, 2787 (2004), Sec. III.B
"Trapping laser" (the "0.16 pN/nm per W" case) and Fig. 10 (the "0.08 pN/nm"
measured example). The gradient-force formula, F = (2 pi n_m a^3 / c)
((m^2-1)/(m^2+2)) grad(I), is the standard Rayleigh point-dipole result
used throughout the optical-trapping literature -- see the KNOWN
DISCREPANCY note above for exactly which printed source it does and does
not match. (Ashkin, Biophys. J. 61, 569 (1992) is the ray-optics, large-
particle regime, a different approximation than this script's point-
dipole treatment, and is not the source for this formula.) The beam-waist
approximation w0 = lambda_vac / (pi * NA) is a standard diffraction-limit
estimate, not a substitute for a full aperture-filling calculation.
"""
import argparse
import math
import sys

# ---- inputs, all assumed or quoted (see docstring) -------------------------
C = 2.998e8                # m/s, speed of light in vacuum
LAMBDA_VAC = 1064e-9       # m, Nd:YAG trapping laser, the field standard wavelength
N_WATER = 1.33             # refractive index of water (medium)
N_POLYSTYRENE = 1.57       # Neuman & Block's stated bead index (Sec. III.B)
BEAD_RADIUS = 0.25e-6      # m, half of their stated "0.5 um polystyrene sphere"
NA_OBJECTIVE = 1.2         # their stated objective numerical aperture
REF_POWER = 1.0            # W, power in the specimen plane (their stiffness is quoted "per W")
STATED_STIFFNESS_PER_W = 0.16     # pN/nm per W, Neuman & Block Sec. III.B, a worked example, not itself tied to a citation
FIG10_MEASURED_STIFFNESS = 0.08   # pN/nm, their Fig. 10, a real trapped bead's measured rolloff


def relative_index(n_particle, n_medium):
    return n_particle / n_medium


def polarizability_factor(m):
    """The Clausius-Mossotti factor (m^2-1)/(m^2+2). Zero when m=1 (index-matched, invisible)."""
    return (m ** 2 - 1.0) / (m ** 2 + 2.0)


def beam_waist_from_na(lambda_vac, na):
    """Diffraction-limit estimate, not an aperture-filling calculation."""
    return lambda_vac / (math.pi * na)


def gaussian_intensity(r, i0, w0):
    return i0 * math.exp(-2.0 * r ** 2 / w0 ** 2)


def intensity_log_slope(r, w0):
    """Exact d(ln I)/dr for a Gaussian profile: -4r/w0^2."""
    return -4.0 * r / w0 ** 2


def gradient_force(r, i0, w0, a, n_m, m):
    """Transverse gradient force, Rayleigh point-dipole approximation."""
    mfac = polarizability_factor(m)
    i_r = gaussian_intensity(r, i0, w0)
    di_dr = i_r * intensity_log_slope(r, w0)
    return (2.0 * math.pi * n_m * a ** 3 / C) * mfac * di_dr


def trap_stiffness(power, w0, a, n_m, m):
    """k = -dF/dr at r=0, closed form: k = 16 n_m a^3 (m^2-1)/(m^2+2) P / (c w0^4)."""
    mfac = polarizability_factor(m)
    return (16.0 * n_m * a ** 3 * mfac * power) / (C * w0 ** 4)


def main():
    parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    parser.add_argument("--selftest", action="store_true")
    args = parser.parse_args()
    if args.selftest:
        return selftest()

    print("1. Transverse trap stiffness, Rayleigh point-dipole model")
    w0 = beam_waist_from_na(LAMBDA_VAC, NA_OBJECTIVE)
    m = relative_index(N_POLYSTYRENE, N_WATER)
    mfac = polarizability_factor(m)
    print(f"   beam waist estimate from NA {NA_OBJECTIVE}: {w0*1e9:.0f} nm")
    print(f"   relative index m = {N_POLYSTYRENE}/{N_WATER} = {m:.4f}, (m^2-1)/(m^2+2) = {mfac:.4f}")

    k = trap_stiffness(REF_POWER, w0, BEAD_RADIUS, N_WATER, m)
    k_pn_nm = k * 1000.0  # 1 N/m = 1000 pN/nm
    print(f"   computed stiffness at {REF_POWER} W: {k:.3e} N/m = {k_pn_nm:.2f} pN/nm per W")

    print("\n2. Compared against Neuman & Block's own stated case")
    print(f"   their worked-example stiffness for this exact bead/laser/objective (not itself a cited measurement): {STATED_STIFFNESS_PER_W} pN/nm per W")
    ratio = k_pn_nm / STATED_STIFFNESS_PER_W
    print(f"   ratio, computed / their example: {ratio:.1f}x")
    print(f"   a separate, actually measured trapped bead (their Fig. 10, rolloff-frequency method): {FIG10_MEASURED_STIFFNESS} pN/nm")

    print("\n3. Why the point-dipole estimate overshoots by roughly two orders of magnitude")
    lambda_medium = LAMBDA_VAC / N_WATER
    print(f"   bead radius: {BEAD_RADIUS*1e9:.0f} nm; wavelength in water: {lambda_medium*1e9:.0f} nm; ratio a/lambda_medium = {BEAD_RADIUS/lambda_medium:.2f}")
    print("   Neuman & Block, p. 2789: 'When the dimensions of the trapped particle are")
    print("   comparable to the wavelength of the trapping laser (a ~ lambda), neither the")
    print("   ray optic nor the point-dipole approach is valid.' This bead sits exactly there.")
    print("   The point-dipole model treats the bead as sampling the field gradient at a single")
    print("   point; a bead this size instead averages the gradient over its own volume, and the")
    print("   real focus is not a simple paraxial Gaussian at this NA -- both effects cut the")
    print("   real restoring force well below the naive estimate. Full treatments use Lorenz-Mie")
    print("   or T-matrix electromagnetic theory instead; real instruments calibrate stiffness")
    print("   directly from a trapped bead's own Brownian motion (their Fig. 10), not from this formula.")
    return 0


def selftest():
    ok = True
    w0 = beam_waist_from_na(LAMBDA_VAC, NA_OBJECTIVE)
    m = relative_index(N_POLYSTYRENE, N_WATER)

    # (a) Force is zero exactly at the beam center (symmetry: di_dr=0 at r=0).
    f0 = gradient_force(0.0, 1.0, w0, BEAD_RADIUS, N_WATER, m)
    print(f"[a] force at r=0: {f0:.3e} N (expect 0)")
    ok &= abs(f0) < 1e-30

    # (b) Force is restoring: positive displacement gives a force back toward r=0 (negative).
    f_pos = gradient_force(50e-9, 1.0, w0, BEAD_RADIUS, N_WATER, m)
    print(f"[b] force at r=+50nm: {f_pos:.3e} N (expect < 0, restoring)")
    ok &= f_pos < 0.0

    # (c) Stiffness scales linearly with power.
    k1 = trap_stiffness(1.0, w0, BEAD_RADIUS, N_WATER, m)
    k2 = trap_stiffness(2.0, w0, BEAD_RADIUS, N_WATER, m)
    print(f"[c] k(2W)/k(1W) = {k2/k1:.6f} (expect 2)")
    ok &= abs(k2 / k1 - 2.0) < 1e-9

    # (d) Stiffness scales as 1/w0^4 (halving the waist should give 16x the stiffness).
    k_w0 = trap_stiffness(1.0, w0, BEAD_RADIUS, N_WATER, m)
    k_half_w0 = trap_stiffness(1.0, w0 / 2.0, BEAD_RADIUS, N_WATER, m)
    print(f"[d] k(w0/2)/k(w0) = {k_half_w0/k_w0:.4f} (expect 16)")
    ok &= abs(k_half_w0 / k_w0 - 16.0) < 1e-6

    # (e) Stiffness scales as particle volume, a^3 (doubling radius -> 8x).
    k_a = trap_stiffness(1.0, w0, BEAD_RADIUS, N_WATER, m)
    k_2a = trap_stiffness(1.0, w0, BEAD_RADIUS * 2.0, N_WATER, m)
    print(f"[e] k(2a)/k(a) = {k_2a/k_a:.4f} (expect 8)")
    ok &= abs(k_2a / k_a - 8.0) < 1e-6

    # (f) Index-matched particle (m=1) is invisible: zero polarizability factor, zero stiffness.
    k_matched = trap_stiffness(1.0, w0, BEAD_RADIUS, N_WATER, 1.0)
    print(f"[f] k at m=1 (index-matched): {k_matched:.3e} (expect 0)")
    ok &= abs(k_matched) < 1e-30

    # (g) Hand check of the closed-form stiffness against the two-step derivation
    # (I0 from power and waist, then k from the log-slope formula), not the same
    # code path as trap_stiffness() itself.
    i0 = 2.0 * REF_POWER / (math.pi * w0 ** 2)
    mfac = polarizability_factor(m)
    k_hand = (2.0 * math.pi * N_WATER * BEAD_RADIUS ** 3 / C) * mfac * i0 * (4.0 / w0 ** 2)
    k_fn = trap_stiffness(REF_POWER, w0, BEAD_RADIUS, N_WATER, m)
    err = abs(k_hand - k_fn) / k_fn
    print(f"[g] closed-form vs. independent hand derivation: rel err {err:.1e}")
    ok &= err < 1e-9

    print("SELFTEST", "PASS" if ok else "FAIL")
    return 0 if ok else 1


if __name__ == "__main__":
    sys.exit(main())