cx-numerical-G.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "BE5d9wKNPQXk"
   },
   "source": [
    "# From $d\\Sigma/d\\Omega(Q)$ to $G(\\delta)$ using cosine and Hankel transforms\n",
    "In this notebook, two numerical methods for computing $G(\\delta)$ from the $d\\Sigma/d\\Omega(Q)$ of a dilute solid sphere sample are derived.  First, $G(\\delta)$ is computed using a Hankel transform. Next, a more general method is used to compute $G(\\delta)$ by performing the cosine transform once and integrating once. This is equivalent to the Hankel transform for isotropic samples by the projection-slice theorem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "VfZCRGhRt0y4"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from hankel import *\n",
    "from definitions import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "4Tb4uhVMt336"
   },
   "outputs": [],
   "source": [
    "def form_factor(Q, R):\n",
    "    return (3 * (np.sin(Q * R) - Q * R * np.cos(Q * R))/ (Q * R) ** 3)**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 468
    },
    "id": "r09jNCJuurVj",
    "outputId": "94a33e34-20ab-4a41-c648-1156bc1a5326"
   },
   "outputs": [],
   "source": [
    "# Sample volume thickness\n",
    "t = 10e-3 # m\n",
    "# Radius of sphere, equal to 1 um but expressed in AA\n",
    "# Volume ratio\n",
    "phi = 0.015\n",
    "delta_rho = 1.8e10 * (1e2) ** 2 # 1/m^2 (?)\n",
    "wavelength = 4.321e-10 # m\n",
    "\n",
    "R = 50e-9\n",
    "Q_min = 0.01 / R\n",
    "Q_max = 10 / R\n",
    "Q = np.linspace(Q_min, Q_max, 1000)\n",
    "ff = form_factor(Q,R)\n",
    "sigma_an = s_t(R,t,wavelength,phi,delta_rho)\n",
    "plt.plot(Q * 1e-10, ff)\n",
    "print(f\"Sigma * t = {sigma_an} m3\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hankel transformation implementation\n",
    "The implementation closely follows that of the following paper, using straightforward numpy interpretations of the procedure described\n",
    "\n",
    "Bakker, J. H., Washington, A. L., Parnell, S. R., Van Well, A. A., Pappas, C., & Bouwman, W. G. (2020). Analysis of SESANS data by numerical Hankel transform implementation in SasView. Journal of Neutron Research, 22(1), 57-70. https://doi.org/10.3233/jnr-200154\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "y1iTuWnmxCxE"
   },
   "outputs": [],
   "source": [
    "delta_min = R / 10\n",
    "delta_max = 3 * R\n",
    "Z = 250\n",
    "delta, Qn, dQ, N =setup_Q_delta_arrays(delta_min, delta_max, Z)\n",
    "H_0, H_kernel = create_hankel_kernels(Qn, dQ, delta, Z)\n",
    "I_Q = form_factor(Qn,R) * 0.5 * R ** 2\n",
    "G_full, G_full_0 = compute_G_matrices(I_Q, H_0, H_kernel, Z)\n",
    "G_delta_num, sigma_t_num = compute_total_G_sigma(G_full, G_full_0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 490
    },
    "id": "f7YGmTLKz97z",
    "outputId": "74c3a0e8-3289-4705-d0b9-21d11858b9b3"
   },
   "outputs": [],
   "source": [
    "def G_0(xi):\n",
    "    res = np.zeros_like(xi)\n",
    "    res[xi>=2.0] = 0\n",
    "    valid_xi = xi[xi<2.0]\n",
    "    res[xi<2.0] = np.sqrt(1 - (valid_xi / 2) ** 2) * (1 + valid_xi ** 2 / 8)\\\n",
    "         + 1 / 2 * valid_xi ** 2 * (1 - (valid_xi / 4 ) ** 2) * np.log(valid_xi / (2 + np.sqrt(4 - valid_xi ** 2)))\n",
    "    return res\n",
    "G_0_num = G_delta_num / sigma_t_num\n",
    "\n",
    "xi = delta / R\n",
    "G_0_an = G_0(xi)\n",
    "plt.plot(delta * 1e9,G_0_an)\n",
    "plt.plot(delta * 1e9,G_0_num)\n",
    "plt.xlabel(r'$\\delta$ [nm]')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "tkuJssUL72Xd"
   },
   "source": [
    "# Limited $Q$-range effects using Hankel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 493
    },
    "id": "hCYMCsD-0BUU",
    "outputId": "23ea69d1-aad5-4a8c-86e9-d6a0e7e5a322"
   },
   "outputs": [],
   "source": [
    "G_partials, sigma_partials, G_0_partial, P_partial = compute_partial_G_sigma(G_full, G_full_0, sigma_t_num, Z)\n",
    "plt.plot(Qn * 1e-10, sigma_partials)\n",
    "plt.xlabel(r'$Q$ [$\\AA^{-1}$]')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "nZ0qEu9b--i4",
    "outputId": "cad216b8-780e-4268-a980-d6dad527aa2e"
   },
   "outputs": [],
   "source": [
    "Q_max_plot = 2 * np.pi / R\n",
    "Q_indices = indices_within_range(Qn,0, Q_max_plot)\n",
    "Q_max_plot * 1e-10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 429
    },
    "id": "3ms7LHPn8Rm3",
    "outputId": "ea5663c8-2644-460a-bbb5-80511dc554a2"
   },
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "X, Y = np.meshgrid(Qn[Q_indices],delta, indexing='ij')\n",
    "surf = ax.plot_surface(X * 1e-10, Y  * 1e9, G_partials[Q_indices,:], cmap='viridis', edgecolor='none')\n",
    "cbar = fig.colorbar(surf, ax=ax, shrink=0.5, aspect=10, pad=0.12)\n",
    "ax.set_xlabel(r'$Q$ [$\\AA^{-1}$]')\n",
    "ax.set_ylabel(r'$\\delta$ [nm]')\n",
    "ax.set_zlabel(r'$G_{exp}(\\delta)$')\n",
    "ax.view_init(elev=30, azim=130)  # Adjust these values as needed\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 453
    },
    "id": "pF7BtTLm_7Fw",
    "outputId": "f6535578-122e-4b22-8892-977bf190598b"
   },
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "X, Y = np.meshgrid(Qn[Q_indices],delta, indexing='ij')\n",
    "\n",
    "# Supposed to represent\n",
    "\n",
    "\n",
    "surf = ax.plot_surface(X * 1e-10, Y  * 1e9, G_0_partial[Q_indices], cmap='viridis', edgecolor='none')\n",
    "cbar = fig.colorbar(surf, ax=ax, shrink=0.5, aspect=10, pad=0.12)\n",
    "ax.set_xlabel(r'$Q$ [$\\AA^{-1}$]')\n",
    "ax.set_ylabel(r'$\\delta$ [nm]')\n",
    "ax.set_zlabel(r'$G_{0,exp}(\\delta)$')\n",
    "plt.title(r\"$G_{0,exp}(\\delta)$ as a function of the integrated Q-range\")\n",
    "ax.view_init(elev=30, azim=130)  # Adjust these values as needed\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 460
    },
    "id": "n0-1jrsdSqCb",
    "outputId": "fc8b1629-cea9-42fb-be92-5da25ca4c962"
   },
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "X, Y = np.meshgrid(Qn[Q_indices],delta, indexing='ij')\n",
    "\n",
    "surf = ax.plot_surface(X * 1e-10, Y  * 1e9, P_partial[Q_indices], cmap='viridis', edgecolor='none')\n",
    "cbar = fig.colorbar(surf, ax=ax, shrink=0.5, aspect=10, pad=0.1)\n",
    "ax.set_xlabel(r'$Q$ [$\\AA^{-1}$]')\n",
    "ax.set_ylabel(r'$\\delta$ [nm]')\n",
    "plt.title(r\"$P_{exp}(\\delta)$ as a function of the integrated Q-range\")\n",
    "ax.set_zlabel(r'$P_{exp}(\\delta)$')\n",
    "ax.view_init(elev=30, azim=130)  # Adjust these values as needed\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 472
    },
    "id": "INsW-poqP5Q7",
    "outputId": "e2406a18-a1d0-4135-aade-5c4b9d79ce6d"
   },
   "outputs": [],
   "source": [
    "def get_Q_ix(Q):\n",
    "  return int(Q/dQ)\n",
    "\n",
    "Q_lims = [0.001e10, 0.002e10, 0.003e10, 0.004e10, 0.005e10, 0.007e10, 0.012e10]\n",
    "for Q_lim in Q_lims:\n",
    "  plt.plot(delta * 1e9, G_0_partial[get_Q_ix(Q_lim),:], label=r'$Q_{max}$ = ' + str(Q_lim * 1e-10) + r'$\\AA^{-1}$')\n",
    "  plt.xlabel(r'$\\delta$ [nm]')\n",
    "  plt.ylabel(r'$G_{0,exp}(\\delta)$')\n",
    "plt.legend()\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 476
    },
    "id": "oqPCPx1PRADM",
    "outputId": "7bfbfa32-bcf4-4ef9-872a-3cb7fae68116"
   },
   "outputs": [],
   "source": [
    "for Q_lim in Q_lims:\n",
    "  plt.plot(delta * 1e9, P_partial[get_Q_ix(Q_lim),:], label=r'$Q_{max}$ = ' + str(Q_lim * 1e-10) + r'$\\AA^{-1}$')\n",
    "plt.xlabel(r'$\\delta$ [nm]')\n",
    "plt.ylabel(r'$G_{exp}(\\delta)$')\n",
    "plt.legend()\n",
    "plt.grid()\n",
    "plt.xscale('log')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# $G(\\delta)$ from projection and cosine transform"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def setup_Q_delta_arrays(delta_min, delta_max, Z, Q_max = None):\n",
    "    delta = np.linspace(delta_min, delta_max, Z)\n",
    "    Q_max_delta = 2 * np.pi / delta_min\n",
    "    if Q_max is None:\n",
    "        Q_max = Q_max_delta\n",
    "    else:\n",
    "        Q_max = np.min((Q_max_delta, Q_max))\n",
    "    dQ = 0.1 * 2 * np.pi / (Z * (delta_max - delta_min))\n",
    "    N = int(np.ceil(Q_max / dQ))\n",
    "    Qn = np.arange(-N,N+1) * dQ\n",
    "    Qx, Qy = np.meshgrid(Qn, Qn)\n",
    "    Q = np.sqrt(Qx**2 + Qy**2)\n",
    "    return delta, Qn, dQ, N, Qx, Qy, Q\n",
    "\n",
    "R = 50e-9\n",
    "delta_min = R / 10\n",
    "delta_max = 2.5 * R\n",
    "Z = 18\n",
    "\n",
    "delta, Qn, dQ, N, Qx, Qy, Q =setup_Q_delta_arrays(delta_min, delta_max, Z, None)\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "surf = ax.plot_surface(Qx * 1e-10, Qy * 1e-10, Q * 1e-10, cmap='viridis', edgecolor='none')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import integrate\n",
    "\n",
    "def clean_nan(x):\n",
    "    nan_ix = np.isnan(x)\n",
    "    x[nan_ix] = 0\n",
    "    return x\n",
    "\n",
    "I_Q = clean_nan(form_factor(Q,R) * 0.5 * R ** 2)\n",
    "I_Q_x_integ = integrate.trapezoid(I_Q, dx = dQ, axis = 0)\n",
    "sigma_t = integrate.trapezoid(I_Q_x_integ, dx = dQ, axis = 0)\n",
    "plt.plot(Qn, I_Q_x_integ)\n",
    "print(sigma_t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cos_Qy_delta = np.cos(np.outer(Qn,delta))\n",
    "print(cos_Qy_delta.shape)\n",
    "I_Q_mat = np.tile(I_Q_x_integ, (Z,1)).transpose()\n",
    "F_kernel = cos_Qy_delta * I_Q_mat * dQ\n",
    "# print(F_kernel, F_kernel.shape)\n",
    "\n",
    "for i in [0,int(Z/4),int(Z/2),Z-1]:\n",
    "    plt.plot(Qn * 1e-10, F_kernel[:,i])\n",
    "plt.xlim((-0.025,0.025))\n",
    "plt.xlabel(r'$Q$ [Å$^{-1}$]')\n",
    "plt.ylabel(r'$F(\\delta)$')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "G_integ = np.sum(F_kernel,axis=0)\n",
    "G_0_integ = G_integ / sigma_t\n",
    "plt.plot(delta * 1e9, G_0_integ, '.')\n",
    "xi = delta / R\n",
    "G_0_an = G_0(xi)\n",
    "plt.plot(delta * 1e9,G_0_an)\n",
    "plt.xlabel(r'$\\delta$ [nm]')\n",
    "plt.ylabel(r'$G_0(\\delta)$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Limited $Q$-range effects using Fourier and projection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_total_sigma(R, Z):\n",
    "    delta_min = R / 10\n",
    "    delta_max = 3.0 * R\n",
    "    _, Qn, _, _, _, _, _ =setup_Q_delta_arrays(delta_min, delta_max, Z, None)\n",
    "    Qx, Qy = np.meshgrid(Qn, Qn)\n",
    "    Q = np.sqrt(Qx**2 + Qy**2)\n",
    "    I_Q = clean_nan(form_factor(Q,R) * 0.5 * R ** 2)\n",
    "    I_Q_x_integ = integrate.trapezoid(I_Q, dx = dQ, axis = 0)\n",
    "    sigma_t = integrate.trapezoid(I_Q_x_integ, dx = dQ, axis = 0)\n",
    "    return sigma_t\n",
    "\n",
    "delta_min = R / 10\n",
    "delta_max = 3.0 * R\n",
    "sigma_t_absolute = compute_total_sigma(R,Z)\n",
    "Q_max = 0.004039179014103979e10\n",
    "delta, Qn, dQ, N, Qx, Qy, Q = setup_Q_delta_arrays(delta_min, delta_max, Z, Q_max)\n",
    "\n",
    "I_Q = clean_nan(form_factor(Q,R) * 0.5 * R ** 2)\n",
    "\n",
    "# Compute integrals and sigma_t\n",
    "I_Q_x_integ = integrate.trapezoid(I_Q, dx = dQ, axis = 0)\n",
    "sigma_t = integrate.trapezoid(I_Q_x_integ, dx = dQ, axis = 0)\n",
    "\n",
    "# Compute cosine transform\n",
    "cos_Qy_delta = np.cos(np.outer(Qn,delta))\n",
    "I_Q_mat = np.tile(I_Q_x_integ, (Z,1)).transpose()\n",
    "F_kernel = cos_Qy_delta * I_Q_mat * dQ\n",
    "\n",
    "G_integ = np.sum(F_kernel,axis=0)\n",
    "G_0_integ = G_integ / sigma_t\n",
    "# print(sigma_t)\n",
    "plt.plot(delta * 1e9, G_0_integ, '.')\n",
    "xi = delta / R\n",
    "G_0_an = G_0(xi)\n",
    "plt.plot(delta * 1e9,G_0_an)\n",
    "print(sigma_t_absolute, sigma_t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sigma_an = s_t(R,t,wavelength,phi,delta_rho)\n",
    "P_an = np.exp(sigma_an * (G_0_an - 1))\n",
    "P_part = np.exp(sigma_t/sigma_t_absolute * sigma_an *(G_0_integ - 1))\n",
    "plt.plot(delta * 1e9, P_part, '.')\n",
    "plt.plot(delta * 1e9,P_an)"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3",
   "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.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}