{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Interpretation of Mueller Matrices based on the Polar decomposition, or the Lu-Chipman decomposition\n",
    "\n",
    "_written by Jaren N. Ashcraft_\n",
    "\n",
    "`Katsu` has features to simulate polarimetric instrumentation and reduce polarimetric data, resulting in a Mueller matrix. For further insights, we can look into the Polar decomposition by Lu and Chipman [1]. This technique decomposes a Mueller matrix $\\mathbf{M}$ into its constituent depolarizer $\\mathbf{M_{\\Delta}}$, diattenuator $\\mathbf{M_{D}}$, and retarder $\\mathbf{M_{R}}$, as shown in the following Equation,\n",
    "\n",
    "$$\\mathbf{M} = \\mathbf{M_{\\Delta}}\\mathbf{M_{R}}\\mathbf{M_{D}}. $$\n",
    "\n",
    "This function is critical for separating depolarization from the Mueller matrix. In this tutorial, we review the methods available in `Katsu` to perform this decomposition. Below, we initialize a random diattenuator, retarder, and depolarizer.\n",
    "\n",
    "**References**\n",
    "- [1] Shih-Yau Lu and Russell A. Chipman, \"Interpretation of Mueller matrices based on polar decomposition,\" J. Opt. Soc. Am. A 13, 1106-1113 (1996) https://doi.org/10.1364/JOSAA.13.001106"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Diattenuator\n",
      "[[0.85802071 0.10370314 0.09697308 0.        ]\n",
      " [0.10370314 0.85250275 0.00590091 0.        ]\n",
      " [0.09697308 0.00590091 0.85171026 0.        ]\n",
      " [0.         0.         0.         0.8461923 ]]\n",
      "------------------------------\n",
      "Retarder\n",
      "[[ 1.          0.          0.          0.        ]\n",
      " [ 0.          0.91262605 -0.0027427  -0.40878621]\n",
      " [ 0.         -0.0027427   0.99991391 -0.01283196]\n",
      " [ 0.          0.40878621  0.01283196  0.91253996]]\n",
      "------------------------------\n",
      "Depolarizer\n",
      "[[1.         0.         0.         0.        ]\n",
      " [0.         0.11585728 0.12764087 0.        ]\n",
      " [0.         0.12764087 0.82168709 0.        ]\n",
      " [0.         0.         0.         0.27848385]]\n",
      "------------------------------\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from katsu.mueller import (\n",
    "    linear_diattenuator,\n",
    "    linear_retarder,\n",
    "    depolarizer\n",
    ")\n",
    "\n",
    "# create a Mueller matrix \n",
    "M_diattenuator = linear_diattenuator(np.random.random(), np.random.random())\n",
    "M_retarder = linear_retarder(np.random.random(), np.random.random())\n",
    "M_depolarizer = depolarizer(np.random.random(), np.random.random(), np.random.random(), np.random.random())\n",
    "\n",
    "print('Diattenuator')\n",
    "print(M_diattenuator)\n",
    "print('-'*30)\n",
    "print('Retarder')\n",
    "print(M_retarder)\n",
    "print('-'*30)\n",
    "print('Depolarizer')\n",
    "print(M_depolarizer)\n",
    "print('-'*30)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In reality, optics we measure will be some combination of these three. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "As-measured Mueller Matrix\n",
      "[[ 0.85802071  0.10370314  0.09697308  0.        ]\n",
      " [ 0.02327453  0.09059165  0.1090549  -0.04146236]\n",
      " [ 0.09148722  0.10223165  0.700155   -0.05307461]\n",
      " [ 0.01215214  0.0970703   0.00371534  0.21504085]]\n"
     ]
    }
   ],
   "source": [
    "M_under_test = M_depolarizer @ M_retarder @ M_diattenuator\n",
    "print('As-measured Mueller Matrix')\n",
    "print(M_under_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we wanted to know what component of this is a diattenuator, `Katsu` has the routines in [1] built into `katsu.mueller` to do so. As shown below, the decomposition is capable of extracting the diattenuator from the total Mueller matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Diattenuattor from Polar Decomposition\n",
      "[[0.85802071 0.10370314 0.09697308 0.        ]\n",
      " [0.10370314 0.85250275 0.00590091 0.        ]\n",
      " [0.09697308 0.00590091 0.85171026 0.        ]\n",
      " [0.         0.         0.         0.8461923 ]]\n",
      "------------------------------\n",
      "Diattenuattor we Specified\n",
      "[[0.85802071 0.10370314 0.09697308 0.        ]\n",
      " [0.10370314 0.85250275 0.00590091 0.        ]\n",
      " [0.09697308 0.00590091 0.85171026 0.        ]\n",
      " [0.         0.         0.         0.8461923 ]]\n",
      "------------------------------\n"
     ]
    }
   ],
   "source": [
    "from katsu.mueller import decompose_diattenuator\n",
    "\n",
    "M_d = decompose_diattenuator(M_under_test)\n",
    "print('Diattenuattor from Polar Decomposition')\n",
    "print(M_d)\n",
    "print('-'*30)\n",
    "print('Diattenuattor we Specified')\n",
    "print(M_diattenuator)\n",
    "print('-'*30)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`Katsu` can perform a simmilar operation to parse the retarder from the total Mueller matrix. The `decompose_retarder` method has the option to return both the diattenuator and retarder, which we show below. The diattenuator that returned is the same as what we gave it, but the retarder is not. Why could this be? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Diattenuattor from Polar Decomposition\n",
      "[[0.85802071 0.10370314 0.09697308 0.        ]\n",
      " [0.10370314 0.85250275 0.00590091 0.        ]\n",
      " [0.09697308 0.00590091 0.85171026 0.        ]\n",
      " [0.         0.         0.         0.8461923 ]]\n",
      "------------------------------\n",
      "Diattenuattor we Specified\n",
      "[[0.85802071 0.10370314 0.09697308 0.        ]\n",
      " [0.10370314 0.85250275 0.00590091 0.        ]\n",
      " [0.09697308 0.00590091 0.85171026 0.        ]\n",
      " [0.         0.         0.         0.8461923 ]]\n",
      "------------------------------\n",
      "Retarder from Polar Decomposition\n",
      "[[ 1.00000000e+00  0.00000000e+00 -7.04731412e-18  0.00000000e+00]\n",
      " [ 0.00000000e+00  1.05384290e-01  1.27312114e-01 -4.89987408e-02]\n",
      " [ 1.38777878e-17  1.14234735e-01  8.21266266e-01 -6.27216840e-02]\n",
      " [ 1.73472348e-18  1.13840357e-01  3.57349440e-03  2.54127639e-01]]\n",
      "------------------------------\n",
      "Retarder we Specified\n",
      "[[ 1.          0.          0.          0.        ]\n",
      " [ 0.          0.91262605 -0.0027427  -0.40878621]\n",
      " [ 0.         -0.0027427   0.99991391 -0.01283196]\n",
      " [ 0.          0.40878621  0.01283196  0.91253996]]\n",
      "------------------------------\n"
     ]
    }
   ],
   "source": [
    "from katsu.mueller import decompose_retarder\n",
    "M_r, M_d = decompose_retarder(M_under_test, return_all=True) # if false, just returns retarder\n",
    "\n",
    "print('Diattenuattor from Polar Decomposition')\n",
    "print(M_d)\n",
    "print('-'*30)\n",
    "print('Diattenuattor we Specified')\n",
    "print(M_diattenuator)\n",
    "print('-'*30)\n",
    "\n",
    "print('Retarder from Polar Decomposition')\n",
    "print(M_r)\n",
    "print('-'*30)\n",
    "print('Retarder we Specified')\n",
    "print(M_retarder)\n",
    "print('-'*30)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's because we haven't separated depolarization yet. The structure of a depolarizer looks very much like a kind of retarder. Indeed this is because depolarization can arize from rapidly varying polarization, like scatter. To perform the final step in the decomposition, we call `decompose_depolarizer` with the `return_all=True` keyword argument to get our original Mueller matrix separated into a depolarizer, retarder, and diattenuator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Diattenuattor from Polar Decomposition\n",
      "[[0.85802071 0.10370314 0.09697308 0.        ]\n",
      " [0.10370314 0.85250275 0.00590091 0.        ]\n",
      " [0.09697308 0.00590091 0.85171026 0.        ]\n",
      " [0.         0.         0.         0.8461923 ]]\n",
      "------------------------------\n",
      "Diattenuattor we Specified\n",
      "[[0.85802071 0.10370314 0.09697308 0.        ]\n",
      " [0.10370314 0.85250275 0.00590091 0.        ]\n",
      " [0.09697308 0.00590091 0.85171026 0.        ]\n",
      " [0.         0.         0.         0.8461923 ]]\n",
      "------------------------------\n"
     ]
    }
   ],
   "source": [
    "from katsu.mueller import decompose_depolarizer\n",
    "\n",
    "M_a, M_r, M_d = decompose_depolarizer(M_under_test, return_all=True) # if false, just returns depolarizer\n",
    "\n",
    "print('Diattenuattor from Polar Decomposition')\n",
    "print(M_d)\n",
    "print('-'*30)\n",
    "print('Diattenuattor we Specified')\n",
    "print(M_diattenuator)\n",
    "print('-'*30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Retarder from Polar Decomposition\n",
      "[[ 1.00000000e+00  0.00000000e+00 -7.04731412e-18  0.00000000e+00]\n",
      " [-2.69840694e-33  9.12626052e-01 -2.74270322e-03 -4.08786212e-01]\n",
      " [ 3.06626967e-33 -2.74270322e-03  9.99913905e-01 -1.28319629e-02]\n",
      " [ 7.70371978e-34  4.08786212e-01  1.28319629e-02  9.12539957e-01]]\n",
      "------------------------------\n",
      "Retarder we Specified\n",
      "[[ 1.          0.          0.          0.        ]\n",
      " [ 0.          0.91262605 -0.0027427  -0.40878621]\n",
      " [ 0.         -0.0027427   0.99991391 -0.01283196]\n",
      " [ 0.          0.40878621  0.01283196  0.91253996]]\n",
      "------------------------------\n"
     ]
    }
   ],
   "source": [
    "print('Retarder from Polar Decomposition')\n",
    "print(M_r)\n",
    "print('-'*30)\n",
    "print('Retarder we Specified')\n",
    "print(M_retarder)\n",
    "print('-'*30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Depolarizer from Polar Decomposition\n",
      "[[ 1.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00]\n",
      " [ 0.00000000e+00  1.15857278e-01  1.27640865e-01 -6.81314585e-18]\n",
      " [ 1.38777878e-17  1.27640865e-01  8.21687089e-01 -5.84221487e-18]\n",
      " [ 1.73472348e-18 -6.81314585e-18 -5.84221487e-18  2.78483848e-01]]\n",
      "------------------------------\n",
      "Depolarizer we Specified\n",
      "[[1.         0.         0.         0.        ]\n",
      " [0.         0.11585728 0.12764087 0.        ]\n",
      " [0.         0.12764087 0.82168709 0.        ]\n",
      " [0.         0.         0.         0.27848385]]\n",
      "------------------------------\n"
     ]
    }
   ],
   "source": [
    "print('Depolarizer from Polar Decomposition')\n",
    "print(M_a)\n",
    "print('-'*30)\n",
    "print('Depolarizer we Specified')\n",
    "print(M_depolarizer)\n",
    "print('-'*30)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And with the full depolarizing decomposition, we successfully return the depolarizer, retarder, and diattenuator that we put into the initial Mueller matrix. This tool is a powerful method of decomposing a Mueller matrix into components that are easier to digest. "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "katsudev",
   "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.10.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}