examples/2d_reflection_matrix_Gaussian_beams/reflection_matrix_Gaussian_beams.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "2bf0dcfe-95cb-4b54-a7fc-83ae9af15dbe",
   "metadata": {},
   "source": [
    "# Reflection Matrix in Gaussian-Beam Basis"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83ba0751-ee17-4a6e-ae57-0bac2a7bb04a",
   "metadata": {},
   "source": [
    "In this example, we compute the reflection matrix of an open system with the input and output bases being Gaussian beams focused to different spatial locations. The magnitude of the diagonal elements of this matrix corresponds to what is measured in confocal microscopy. We first build the list of inputs B and list of outputs C, and then use mesti() to compute the reflection matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a92749f-9678-4e4f-b520-d57817aeeae2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Call necessary packages\n",
    "using MESTI, GeometryPrimitives, Plots"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d1f6311b-09d2-403b-bc85-91b39016c3cc",
   "metadata": {},
   "source": [
    "# Build the system\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe867b1b-4a9e-45fc-b3af-84c401487793",
   "metadata": {},
   "source": [
    "In this example, we consider a test system of a dielectric cylindrical scatterer located at (*y*<sub>0</sub>, *z*<sub>0</sub>) in air."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9cc40aa2-79cd-4907-9485-dc283e447fe5",
   "metadata": {},
   "outputs": [],
   "source": [
    "syst = Syst()\n",
    "syst.length_unit = \"µm\"\n",
    "syst.wavelength = 1.0 # wavelength (µm)\n",
    "syst.dx = syst.wavelength/15 # discretization grid size (µm)\n",
    "pml_npixels = 20    # number of PML pixels\n",
    "W = 20       # width of simulation domain (including PML) (µm)\n",
    "L = 10       # length of simulation domain (including PML) (µm)\n",
    "r_0 = 0.75   # cylinder radius (µm)\n",
    "n_bg   = 1.0 # refractive index of the background\n",
    "n_scat = 1.2 # refractive index of the cylinder\n",
    "y_0 = W/2    # location of the cylinder\n",
    "z_0 = L/2    # location of the cylinder\n",
    "\n",
    "# Build the relative permittivity profile from subpixel smoothing\n",
    "domain = Cuboid([W/2,L/2], [W,L]) # domain for subpixel smoothing cetering at (W/2, L/2) and with witdth W and thickness L\n",
    "domain_epsilon = n_bg^2 # epsilon of the domain for subpixel smoothing\n",
    "object = [Ball([y_0, z_0], r_0)] # object for subpixel smoothing: cylinder locating at (W/2, L/2) with radius r_0\n",
    "object_epsilon = [n_scat^2] # epsilon of the object for subpixel smoothing\n",
    "yBC = \"PEC\"; zBC = \"PEC\" # boundary conditions\n",
    "epsilon_xx = mesti_subpixel_smoothing(syst.dx, domain, domain_epsilon, object, object_epsilon, yBC, zBC) # obtaining epsilon_xx from mesti_subpixel_smoothing()\n",
    "\n",
    "ny_Ex, nz_Ex = size(epsilon_xx)\n",
    "y = syst.dx:syst.dx:(W-syst.dx/2)\n",
    "z = syst.dx:syst.dx:(L-syst.dx/2)\n",
    "\n",
    "# Plot the relative permittivity profile\n",
    "plt1= heatmap(z, y, epsilon_xx, \n",
    "              aspect_ratio=:equal, \n",
    "              xlabel = \"z (µm)\", ylabel = \"y (µm)\", \n",
    "              xlims=(0,10), title=\"εᵣ(y,z)\",   \n",
    "              c = cgrad(:grayC, rev=true))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "472f1246-a086-44c0-8d3b-e96d52fd1f02",
   "metadata": {},
   "source": [
    "# Build the input sources"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "50a506e7-95b1-4e26-8e1e-75bdaba5df7e",
   "metadata": {},
   "source": [
    "We consider inputs being Gaussian beams focused at (*y*<sub>f</sub>, *z*<sub>f</sub>). In this example, we fix the focal depth at  *z*<sub>f</sub> = *z*<sub>0</sub> (i.e., the depth of the scatterer), and scan the transverse coordinate *y*<sub>f</sub> of the focus.\n",
    "\n",
    "Perfect Gaussian beams can be generated with the total-field/scattered-field (TF/SF) method. But since the cross section of the beam decays exponentially in *y*, we can generate Gaussian beams to a high accuracy simply by placing line sources at a cross section on the low side, which is what we do here. We place the line sources at *z* = *z*<sub>s</sub> (location of the source plane), just in front of the PML.\n",
    "\n",
    "To determine the required line sources, we (1) take the field profile of the desired incident Gaussian beam at the focal plane, *E*<sup>in</sup>(*y*,*z*<sub>f</sub>) = *E*<sub>yf</sub>(*y*,*z*<sub>f</sub>) = exp(-(*y* - *y*<sub>f</sub>)<sup>2</sup>/*w*<sup>2</sup>), (2) project it onto the propagating channels (i.e., ignoring evanescent contributions) of free space, (3) back propagate it to the source plane to determine *E*<sub>yf</sub>(*y*,*z*<sub>s</sub>) in the propagating-channel basis, and (4) determine the line source necessary to generate such *E*<sub>yf</sub>(*y*,*z*<sub>s</sub>)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3d2c2b8f-c7c2-45bd-99ab-9f5df14f4301",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameters of the input Gaussian beams\n",
    "NA = 0.5   # numerical aperture\n",
    "z_f = z_0  # location of the focus in z (fixed)\n",
    "y_f_start = 0.3*W # starting location of the focus in y \n",
    "y_f_end   = 0.7*W # ending location of the focus in y\n",
    "dy_f = syst.wavelength/(10*NA)  # spacing of focal spots in y\n",
    "\n",
    "# Parameters of the line sources\n",
    "n_source = pml_npixels + 1 # index of the source plane\n",
    "z_s = z[n_source] # location of the source plane\n",
    "z_d = z_s # location of the detection plane\n",
    "\n",
    "w_0 = syst.wavelength/(pi*NA) # beam radius at z = z_f\n",
    "y_f = y_f_start:dy_f:y_f_end # list of focal positions in y\n",
    "M_in = length(y_f) # number of inputs\n",
    "\n",
    "# Step 1: Generate the list of E_yf(z=z_f, y).\n",
    "# Here, y is an ny_Ex-by-1 column vector, and y_f is a M_in-by-1 column vector.\n",
    "# So, y .- transpose(y_f) is an ny_Ex-by-M_in matrix by implicit expansion.\n",
    "# Then, E_yf is an ny_Ex-by-M_in matrix whose m-th column is the cross section\n",
    "# of the m-th Gaussian beam at z = z_f.\n",
    "E_yf = exp.(-(y .- transpose(y_f)).^2/(w_0^2)) # size(E_yf) = [ny_Ex, M_in]\n",
    "\n",
    "# Get properties of propagating channels in the free space.\n",
    "# We use PEC as the boundary condition for such channels since the default\n",
    "# boundary condition in mesti() is PEC for TM waves, but the choice has\n",
    "# little effect since E_yf should be exponentially small at the boundary of\n",
    "# the simulation domain.\n",
    "channels = mesti_build_channels(ny_Ex, \"PEC\", (2*pi/syst.wavelength)*syst.dx, n_bg^2)\n",
    "\n",
    "# Transverse profiles of the propagating channels. Each column of f_transverse is\n",
    "# one transverse profile. Different columns are orthonormal.\n",
    "f_transverse = channels.f_x_m(channels.kydx_prop) # size(f_transverse) = [ny_Ex, N_prop]\n",
    "\n",
    "# Step 2: Project E_yf(y, z_f) onto the propagating channels and obtain v_f.\n",
    "# sqrt_nu_prop is the square root of nu for the propagating channel.\n",
    "sqrt_nu_prop = reshape(channels.sqrt_nu_prop, :, 1)\n",
    "# amplitude coefficients at z = z_f\n",
    "v_f =  (sqrt_nu_prop.* adjoint(f_transverse))*E_yf # size(v_f) = [N_prop, M_in]\n",
    "\n",
    "# Step 3: Back propagate from z = z_f to z = z_s.\n",
    "# This step assumes a PEC boundary in y, so it is not exact with PML in y,\n",
    "# but it is sufficiently accurate since E_yf decays exponentially in y.\n",
    "# Note we use implicit expansion here.\n",
    "kz = reshape(channels.kzdx_prop/syst.dx, :, 1) # list of wave numbers\n",
    "# amplitude coefficients at z = z_s\n",
    "v_s = exp.(1im*kz*(z_s-z_f)).*v_f # size(v_s) = [N_prop, M_in]\n",
    "\n",
    "# Step 4: Determine the line sources.\n",
    "# In a closed geometry with no PML in y, a line source of\n",
    "# -2i*sqrt(nu[a])*f_transverse[:,a] generates outgoing waves with transverse profile\n",
    "# f_transverse[:,a]/sqrt(nu[a]). With PML in y, this is not strictly true but is sufficiently\n",
    "# accurate since E_yf(y,z=z_s) decays exponentially in y.\n",
    "# Note we use implicit expansion here.\n",
    "B_low = (f_transverse.*transpose(channels.sqrt_nu_prop))*v_s # size(B_low) = [ny_Ex, M_in]\n",
    "\n",
    "# We take the -2i prefactor out, to be multiplied at the end. The reason\n",
    "# will be clear when we handle C below.\n",
    "opts = Opts()\n",
    "opts.prefactor = -2im\n",
    "\n",
    "# In mesti(), Bx.pos = [m1, l1, w, h] specifies the position of a\n",
    "# block source, where (m1, l1) is the index of the smaller-(y,z) corner,\n",
    "# and (w, h) is the width and height of the block. Here, we put line\n",
    "# sources (w=1) at l1 = l_source that spans the whole width of the\n",
    "# simulation domain (m1=1, h=ny_Ex).\n",
    "Bx = Source_struct()\n",
    "Bx.pos = [[1, n_source, ny_Ex, 1]]\n",
    "\n",
    "# Bx.data specifies the source profiles inside such block, with\n",
    "# Bx.data[1][:, a] being the a-th source profile.\n",
    "Bx.data = [B_low]\n",
    "\n",
    "# We check that the input sources are sufficiently localized with little \n",
    "# penetration into the PML; otherwise the Gaussian beams will not be\n",
    "# accurately generated.\n",
    "plt2 = heatmap(1:M_in, collect(y), abs.(B_low), \n",
    "              aspect_ratio=:equal, \n",
    "              xlabel = \"Input index\", ylabel = \"y (µm)\", \n",
    "              title=\"|B_low|\", c =cgrad(:grayC, rev=true))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2f6d22f0-1874-450f-971a-18ba9d99e2ff",
   "metadata": {},
   "source": [
    "# Build the output projections"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "df7d1f03-1061-4626-8feb-4792281d5d24",
   "metadata": {},
   "source": [
    "We consider output projections onto the same set of Gaussian beams *$\\Psi$*<sub>yf</sub>(*y*,*z*<sub>f</sub>) = exp(-(*y* - *y*<sub>f</sub>)<sup>2</sup>/*w*<sup>2</sup>) focused at (*y*<sub>f</sub>, *z*<sub>f</sub>),with the projection done at the detection plane (*z* = *z*<sub>d</sub>) same as the source plane (*z* = *z*<sub>s</sub>).\n",
    "\n",
    "When the system has a closed boundary in y, as is the case in mesti2s(), the set of transverse modes form a complete and orthonormal basis, so it is clear what the output projection should be. But the Gaussian beams here are not orthogonal to each other, are not normalized, and do not form a complete basis. So, it is not obvious how our output projection should be defined in the detection plane.\n",
    "\n",
    "What we do here is to project Gaussian beam at the focal plane onto the complete and orthonormal basis of transverse modes, do the propagating from focal to detection plane, and find the projection profiles. Specifically, we (1) project the Gaussian beam onto the propagating channels (i.e., ignoring evanescent contributions) of free space; (2) propagate such Gaussian beam in transverse mode basis to the detection plane at *z* = *z*<sub>d</sub> (3) reconstruct the Gaussian beams at the detection plane, and (4) construct the projection profile for the Gaussian beams at the detection plane.\n",
    "\n",
    "Above, the incident field *E*<sup>in</sup>(*y*,*z*) was not subtracted. Contribution from the incident field will be subtracted using matrix D in the next step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0183c344-881f-48ec-a36a-9a2cae3f30e4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# We perform the output projection on the same plane as the line source\n",
    "Cx = Source_struct()\n",
    "Cx.pos = Bx.pos\n",
    "\n",
    "# size(Psi_yf) = [ny_Ex, num_of_Gaussian_input]\n",
    "Psi_yf = exp.(-(y .- transpose(y_f)).^2/(w_0^2))\n",
    "# Step 1: projecting Psi_yf onto the propagating channels and obtaining \n",
    "# amplitude coefficients v_f_tilde \n",
    "# size(v_f_tilde) = [num_of_prop_channel, num_of_Gaussian_output]\n",
    "v_f_tilde = (sqrt_nu_prop.*adjoint(f_transverse))*Psi_yf\n",
    "# Step 2: propagating the amplitude coefficient from plane z = z_f to plane z = z_d\n",
    "# size(v_d) = [num_of_prop_channel, num_of_Gaussian_output]\n",
    "v_d = exp.(1im*(-kz)*(z_d-z_f)).*v_f_tilde\n",
    "\n",
    "# Step 3: Reconstruct Psi_yf at the detection plane z = z_d\n",
    "# Psi_yd_zd = (f_transverse./transpose(sqrt_nu_prop))*v_d\n",
    "\n",
    "# Step 4: We obtain the projection profile based on Psi_yd and v_d.\n",
    "C_low = adjoint(v_d)*(sqrt_nu_prop.*adjoint(f_transverse))\n",
    "\n",
    "# Normally, the next step would be\n",
    "# C_struct.data = C_low.';\n",
    "# However, we can see that C_low equals transpose(B_low)\n",
    "println(\"max(|C_low-transpose(B_low)|) = $(maximum(abs.(C_low - transpose(B_low))))\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c183f528-b91a-4a62-82d5-183260f2e927",
   "metadata": {},
   "outputs": [],
   "source": [
    "# That means we will have C = transpose(B). So, we can save some computing\n",
    "# time and memory usage by specifying C = transpose(B).\n",
    "# This is expected by reciprocity -- when the set of inputs equals the set\n",
    "# of outputs, we typically have C = transpose(B) or its permutation.\n",
    "C_x = nothing\n",
    "C = \"transpose(B)\";"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d459dd37-5607-47e8-ad50-c8ed00a071bb",
   "metadata": {},
   "source": [
    "# Compute reflection matrix in Gaussian-beam basis"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "613c146e-1d2d-4bea-8fd2-22645027bed6",
   "metadata": {},
   "source": [
    "The scattering matrix is given by S = C\\*inv(A)\\*B - D, with D = C\\*inv(A<sub>0</sub>)\\*B - S<sub>0</sub> where A<sub>0</sub> is a reference system for which its scattering matrix S<sub>0</sub> is known. We consider A<sub>0</sub> to be a homogeneous space with no scatterers, for which the reflection matrix S<sub>0</sub> is zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2d41219-0ae1-4373-8931-f0e878e13fee",
   "metadata": {},
   "outputs": [],
   "source": [
    "pml = PML(pml_npixels)\n",
    "pml.direction = \"all\"\n",
    "syst.PML = [pml]  # PML on all four sides\n",
    "\n",
    "# For a homogeneous space, the length of the simulation domain doesn't\n",
    "# matter, so we choose a minimal thickness of nz_Ex_temp = n_source + pml_npixels\n",
    "# where n_source = pml_npixels + 1 is the index of the source plane.\n",
    "syst.epsilon_xx = n_bg^2*ones(ny_Ex, n_source + pml_npixels)\n",
    "D, _ = mesti(syst, [Bx], C, opts);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3930ef4d-296d-44e8-829b-b89e1b5d6b16",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute the reflection matrix\n",
    "syst.epsilon_xx = epsilon_xx\n",
    "r, _ = mesti(syst, [Bx], C, D, opts);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4ad2d91-9205-46b1-a009-2503bde105db",
   "metadata": {},
   "source": [
    "# Compute the full field profile"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "96e7c62e-0f99-4cd6-8a4a-db95e1b54389",
   "metadata": {},
   "source": [
    "For most applications, it is not necessary to compute the full field profile, since most experiments measure properties in the far field. Here, we compute the full field profile for the purpose of visualizing the system as the incident Gaussian beams are scanned across *y*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e5daf14-3934-4e76-900e-0c8f200e0141",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Exclude the PML pixels from the returned field profiles\n",
    "opts.exclude_PML_in_field_profiles = true\n",
    "y_interior = y[(pml_npixels+1):(ny_Ex-pml_npixels)]\n",
    "z_interior = z[(pml_npixels+1):(nz_Ex-pml_npixels)]\n",
    "\n",
    "field_profiles, _ = mesti(syst, [Bx], opts);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4463eb75-321e-47c4-860a-1bec1791c130",
   "metadata": {},
   "source": [
    "# Animate the field profiles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e68c449b-3f55-4a98-a66b-f53aba9e3bac",
   "metadata": {},
   "outputs": [],
   "source": [
    "theta = range(0, 2*pi, 100)\n",
    "circ_x = cos.(theta)\n",
    "circ_y = sin.(theta)\n",
    "\n",
    "# Loop through Gaussian beams focused at different locations.\n",
    "anim = @animate for ii ∈ 1:M_in\n",
    "    plt1 = heatmap(collect(z_interior), collect(y_interior), real(field_profiles[:, :, ii]), \n",
    "                   xlims=(z_interior[1], z_interior[end]), aspect_ratio=:equal, \n",
    "                   c = :balance, xticks = false, yticks = false, colorbar = false)\n",
    "    plot!(plt1, z_0.+r_0*circ_x, y_0.+r_0*circ_y, lw=1, color=:black, \n",
    "          xlims=(z_interior[1], z_interior[end]), ylims=(y_interior[1], y_interior[end]), legend=false)\n",
    "\n",
    "    plt2 = heatmap(collect(y_f), collect(y_f), abs.(r).^2, \n",
    "                   xlims=(6,14), aspect_ratio=:equal, \n",
    "                   xlabel = \"Input position\", ylabel = \"Output position\",\n",
    "                   c =cgrad(:grayC, rev=true), xticks = false, yticks = false, colorbar = false)\n",
    "    plot!(plt2, [y_f[ii]], seriestype = :vline, linecolor=:blue, legend=false)    \n",
    "    plot(plt1, plt2, layout = (1, 2))  \n",
    "end\n",
    "gif(anim, \"reflection_matrix_Gaussian_beams.gif\", fps = 10)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.9.3",
   "language": "julia",
   "name": "julia-1.9"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.9.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}