c4-form-factor-modelling.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analytical form factor\n",
    "Dilute SANS samples have $Q$-ranges corresponding to their form factor. The above simulates a very naive scattering model in which $Q_y = Q$ is used. Next, proper radial scattering is considered."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from definitions import *\n",
    "from scipy import integrate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def form_factor(Q, R):\n",
    "    return (3 * (np.sin(Q * R) - Q * R * np.cos(Q * R))/ (Q * R) ** 3)**2\n",
    "\n",
    "Rs = [50e-9, 300e-9, 2000e-9]\n",
    "Q = np.linspace(0, 0.015e10,100000)\n",
    "# for R in Rs:\n",
    "R = 100e-9\n",
    "ff = form_factor(Q,R)\n",
    "plt.plot(Q * 1e-10,ff)\n",
    "# plt.xscale('log')\n",
    "plt.yscale('log')\n",
    "plt.xlabel(r'$Q$ [Å$^{-1}$]')\n",
    "plt.ylabel(r'$P(Q)$')\n",
    "plt.xlim((0.0,0.015))\n",
    "\n",
    "plt.ylim((1e-7,1.1))\n",
    "plt.axvline(0.1 / R * 1e-10, linestyle='--', color='red', label=r'$Q_\\text{min} = 0.0001$ Å$^{-1}$')\n",
    "plt.axvline(10 / R * 1e-10, linestyle='--', color='green', label=r'$Q_\\text{max} = 0.01$ Å$^{-1}$')\n",
    "plt.legend()\n",
    "plt.grid()\n",
    "plt.savefig(\"docs/analytical-P-log.eps\",bbox_inches=\"tight\", pad_inches=0, format='eps')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1D scattering model with $Q = Q_y$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "L_s = 1.8\n",
    "\n",
    "def form_factor_y_naive(Q, R):\n",
    "    Q = np.abs(Q)\n",
    "    return (3 * (np.sin(Q * R) - Q * R * np.cos(Q * R))/ (Q * R) ** 3)**2\n",
    "\n",
    "def Qy_to_y(Q_y, s_y, k):\n",
    "    y = L_s * np.tan(2 * np.arcsin(Q_y / (2 * k)))\n",
    "    return y + s_y\n",
    "\n",
    "def y_to_Qy(y, s_y, k):\n",
    "    Qy = 2 * k * np.sin(np.arctan((y - s_y) / L_s)/2)\n",
    "    return Qy\n",
    "\n",
    "def form_factor_y(y,s_y, k, R):\n",
    "    Qy = y_to_Qy(y,s_y, k)\n",
    "    return form_factor_y_naive(Qy,R)\n",
    "    \n",
    "def ff_area(R, k):\n",
    "    y_lim = 0.025e-7 / R\n",
    "    y_oversized = np.linspace(-y_lim,y_lim, 10000)\n",
    "    dy = y_oversized[1] - y_oversized[0]\n",
    "    Q = y_to_Qy(y_oversized,0,k)\n",
    "    ff = form_factor_y_naive(Q,R)\n",
    "    ff_area = integrate.trapezoid(ff) * dy\n",
    "    # print(ff_area)\n",
    "    return ff_area\n",
    "\n",
    "R = 50e-9\n",
    "Q_max = 0.015e10\n",
    "y_lim = 0.025e-7 / R\n",
    "y_oversized = np.linspace(-y_lim,y_lim, 10000)\n",
    "# dQ = Q[1] - Q[0]\n",
    "norm_map = {}\n",
    "for lambda_0 in [4.321e-10, 8e-10]:\n",
    "    k = 2 * np.pi / lambda_0\n",
    "    Q = y_to_Qy(y_oversized,0,k)\n",
    "    ff = form_factor_y_naive(Q,R)\n",
    "    # plt.plot(Q * 1e-10,ff)\n",
    "    ff_norm_factor = ff_area(R, k)\n",
    "    ff_norm = ff / ff_norm_factor\n",
    "    # norm_map[lambda_0] = ff_area\n",
    "    # print(ff_area)\n",
    "    # ff_norm = ff / ff_area\n",
    "    y = Qy_to_y(Q, 0, k=k)\n",
    "    plt.plot(y * 1e3,ff_norm, label=rf'$\\lambda_0 = {lambda_0 * 1e10}Å$')\n",
    "    # plt.yscale('log')\n",
    "    plt.legend()\n",
    "    plt.xlabel(r'$y$ [mm]')\n",
    "    plt.ylabel(r'$P(y)$')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Q_max = 0.015e10\n",
    "# Q = np.linspace(-Q_max, Q_max,100000)\n",
    "h_d = detector_size\n",
    "y = np.linspace(-h_d / 2, h_d/2, 4000)\n",
    "# dQ = Q[1] - Q[0]\n",
    "dy = y[1] - y[0]\n",
    "print(R)\n",
    "for (lambda_0, delta_limit) in [(4.321e-10, 70.71e-9), (8e-10,130.91)]:\n",
    "    for R in [50e-9, 100e-9, 200e-9, 2000e-9]:\n",
    "        k = 2 * np.pi / lambda_0\n",
    "        N_sample_points = 4000\n",
    "        s_y_range = np.linspace(-0.005,0.005,N_sample_points)\n",
    "        I_t = np.zeros_like(y)\n",
    "        norm_factor = ff_area(R, k)\n",
    "        s_y_grid, y_grid = np.meshgrid(s_y_range, y, indexing='ij')\n",
    "        ff = form_factor_y(y_grid, s_y_grid, k, R) / (N_sample_points * norm_factor)\n",
    "        # Sum along the s_y dimension and scale by 0.01\n",
    "        I_t = np.sum(ff, axis=0) * 0.01\n",
    "        \n",
    "        # The slightly more readable, non-vectorized form is this:\n",
    "        #\n",
    "        # for s_y in s_y_range:\n",
    "        #     ff = form_factor_y(y, s_y, k, R) / (N_sample_points * norm_factor)\n",
    "        #     # plt.plot(y, ff)\n",
    "        #     # plt.plot(y[indices] * 1e3,ff[indices])\n",
    "        #     I_t += ff * 0.01\n",
    "\n",
    "        plt.plot(y * 1e3, I_t, label=rf'$R = {round(R * 1e9)}$nm')\n",
    "        plt.xlabel(r'$y$ [mm]')\n",
    "        plt.ylabel(r'$I(y)$')\n",
    "    plt.legend()\n",
    "    plt.xlim((-h_d/2 * 1e3,h_d/2 * 1e3))\n",
    "    plt.grid()\n",
    "    plt.savefig(f\"docs/simplified-I-scattering-{lambda_0 * 1e10}.eps\",bbox_inches=\"tight\", pad_inches=0, format='eps')\n",
    "\n",
    "    plt.show()\n",
    "# y  = np.linspace(-10e-3, 10e-3)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2D scattering model\n",
    "The above model simply equates $Q_y = Q$, whereas in fact $P(Q)$ says something about a radial distribution. In the actual detector, for a given $Q_y$ or $y$, all $Q_x$/$x$ will be integrated over the long detector element. This corresponds also to the $G_\\text{exp}(\\delta)$ expression which performs a cosine transform in one dimension and simply integrates in the other. This can be accounted for by consider the true radially symmetric scattering due to the sphere form factor $P(Q)$, which is done here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "PLOT_3D = True\n",
    "\n",
    "N_x = 1000\n",
    "N_y = 1000\n",
    "\n",
    "def form_factor_2D(x, s_x, y, s_y, k, R):\n",
    "    Qx = y_to_Qy(x,s_x, k)\n",
    "    Qy = y_to_Qy(y,s_y, k)\n",
    "    Q = np.sqrt(Qx ** 2 + Qy ** 2)\n",
    "    return form_factor(Q,R)\n",
    "\n",
    "    \n",
    "def ff_norm_factor_2D(R, k):\n",
    "    x_lim = 0.025e-7 / R\n",
    "    y_lim = x_lim\n",
    "    x = np.linspace(-x_lim,x_lim, N_x)\n",
    "    y = np.linspace(-y_lim,y_lim, N_y)\n",
    "    X, Y = np.meshgrid(x, y)\n",
    "    Z = form_factor_2D(X, 0, Y, 0, k, R)\n",
    "    dx = x[1] - x[0]\n",
    "    dy = y[1] - y[0]\n",
    "    Z_y = integrate.trapezoid(Z, axis = 1, dx = dx)\n",
    "    scale_factor = integrate.trapezoid(Z_y, dx = dy)\n",
    "    return scale_factor\n",
    "\n",
    "def plot_3D_intensity(X, Y, I):\n",
    "    if PLOT_3D:\n",
    "        fig = plt.figure()\n",
    "        ax = fig.add_subplot(111, projection='3d')\n",
    "        surf = ax.plot_surface(X * 1e3, Y * 1e3, I * 1e-6, cmap='viridis', edgecolor='none')\n",
    "        cbar = fig.colorbar(surf, ax=ax, shrink=0.5, aspect=10, pad=0.1)\n",
    "        cbar.set_label(r'$I_{(s_x,s_y) = (0,0)}$ [a.u.]', fontsize=12)\n",
    "        ax.set_xlabel(r'$x$ [mm]')\n",
    "        ax.set_ylabel(r'$y$ [mm]')\n",
    "        plt.savefig(f\"docs/simplified-I-scattering-3D-plot-point-{lambda_0 * 1e10}.eps\",bbox_inches=\"tight\", pad_inches=0, format='eps')\n",
    "        plt.show()\n",
    "    else:\n",
    "        x_lim, y_lim = np.min(X[0,:]), np.min(Y[0,:])\n",
    "        plt.imshow(I, extent=[-x_lim, x_lim, -y_lim, y_lim], origin='lower', cmap='viridis')\n",
    "        plt.colorbar(label='Function Value')\n",
    "        plt.xlabel('x')\n",
    "        plt.ylabel('y')\n",
    "        plt.show()\n",
    "for (lambda_0, delta_limit) in [(4.321e-10, 70.71e-9), (8e-10,130.91)]:\n",
    "    k = 2 * np.pi / lambda_0\n",
    "    R = 200e-9\n",
    "    x_lim = h_d / 2\n",
    "    y_lim = x_lim\n",
    "    x = np.linspace(-x_lim,x_lim, N_x)\n",
    "    y = np.linspace(-y_lim,y_lim, N_y)\n",
    "    X, Y = np.meshgrid(x, y, indexing='ij')\n",
    "    Z = form_factor_2D(X, 0, Y, 0, k, R)\n",
    "    scale_factor = ff_norm_factor_2D(R,k)\n",
    "    print(scale_factor)\n",
    "    Z_unit = Z / scale_factor\n",
    "    plot_3D_intensity(X,Y,Z_unit)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Integrate intensity curves from full sample 2D intensity pattern"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h_d = detector_size\n",
    "h_s = 0.01\n",
    "N_x = 200\n",
    "N_y = 200\n",
    "N_sample_points = 100\n",
    "dy = y[1] - y[0]\n",
    "print(R)\n",
    "I_t_map = {}\n",
    "for (lambda_0, delta_limit) in [(4.321e-10, 70.71e-9), (8e-10,130.91)]:\n",
    "    for R in [50e-9, 100e-9, 200e-9, 300e-9, 2000e-9]:\n",
    "        k = 2 * np.pi / lambda_0\n",
    "        x_lim = h_d / 2\n",
    "        y_lim = x_lim\n",
    "        x = np.linspace(-x_lim,x_lim, N_x)\n",
    "        y = np.linspace(-y_lim,y_lim, N_y)\n",
    "        X, Y = np.meshgrid(x, y, indexing='ij')\n",
    "        scale_factor = ff_norm_factor_2D(R,k)\n",
    "        s_x_range = np.linspace(-h_s/2,h_s/2,N_sample_points)\n",
    "        s_y_range = np.linspace(-h_s/2,h_s/2,N_sample_points)\n",
    "        I_t = np.zeros_like(X)\n",
    "        for s_x in s_x_range:\n",
    "            s_y_grid, X_grid, Y_grid = np.meshgrid(s_y_range, x, y, indexing='ij')\n",
    "            ff = form_factor_2D(X_grid, s_x, Y_grid, s_y_grid, k, R) / (N_sample_points * N_sample_points * scale_factor)\n",
    "            I_t += np.sum(ff, axis=0)\n",
    "        I_t *= 0.01\n",
    "        I_t_map[(lambda_0, R)] = I_t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for (lambda_0, delta_limit) in [(4.321e-10, 70.71e-9), (8e-10,130.91)]:\n",
    "    for R in [50e-9, 100e-9, 200e-9, 300e-9, 2000e-9]:\n",
    "        x_lim = h_d / 2\n",
    "        y_lim = x_lim\n",
    "        y = np.linspace(-y_lim,y_lim, N_y)\n",
    "        I_t = I_t_map[(lambda_0, R)]\n",
    "        I_t_y = integrate.trapezoid(I_t, x=y, axis=1)\n",
    "        plt.plot(y * 1e3, I_t_y, label=rf'$R = {round(R * 1e9)}$nm')\n",
    "        plt.xlabel(r'$y$ [mm]')\n",
    "        plt.ylabel(r'$I(y)$')\n",
    "    plt.legend()\n",
    "    plt.xlim((-h_d/2 * 1e3,h_d/2 * 1e3))\n",
    "    plt.grid()\n",
    "    plt.savefig(f\"docs/simplified-I-scattering-{lambda_0 * 1e10}.eps\",bbox_inches=\"tight\", pad_inches=0, format='eps')\n",
    "    plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import pickle\n",
    "# with open('objects/I_t.pkl', 'wb') as file:\n",
    "#     pickle.dump(I_t, file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pickle\n",
    "with open('objects/I_t.pkl', 'rb') as file:\n",
    "    I_t = pickle.load(file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for (lambda_0, delta_limit) in [(4.321e-10, 70.71e-9), (8e-10,130.91)]:\n",
    "    # for R in [50e-9, 100e-9, 200e-9, 300e-9, 2000e-9]:\n",
    "    for R in [200e-9]:\n",
    "        I_t = I_t_map[(lambda_0, R)]\n",
    "        fig = plt.figure()\n",
    "        ax = fig.add_subplot(111, projection='3d')\n",
    "        surf = ax.plot_surface(X * 1e3, Y * 1e3, I_t, cmap='viridis', edgecolor='none')\n",
    "        cbar = fig.colorbar(surf, ax=ax, shrink=0.5, aspect=10, pad=0.1)\n",
    "        cbar.set_label(r'$I$ [a.u.]', fontsize=12)\n",
    "        ax.set_xlabel(r'$x$ [mm]')\n",
    "        ax.set_ylabel(r'$y$ [mm]')\n",
    "        plt.savefig(f\"docs/simplified-I-scattering-3D-plot-{lambda_0 * 1e10}.eps\",bbox_inches=\"tight\", pad_inches=0, format='eps')\n",
    "\n",
    "        plt.show()\n",
    "        plt.imshow(I_t, extent=[-x_lim * 1e3, x_lim * 1e3, -y_lim * 1e3, y_lim * 1e3], origin='lower', cmap='viridis')\n",
    "        plt.colorbar(label=r'$I$ [a.u.]')\n",
    "        plt.xlabel(r'$x$ [mm]')\n",
    "        plt.ylabel(r'$y$ [mm]')\n",
    "        plt.show()\n",
    "        plt.savefig(f\"docs/simplified-I-scattering-2D-image-{lambda_0 * 1e10}.eps\",bbox_inches=\"tight\", pad_inches=0, format='eps')\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "bep",
   "language": "python",
   "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": 2
}