{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Stand-Alone Tutorial of S2L2 Representation\n", "This tutorial provide a stand-alone implementation of spin-2 SH with order $l=2$ and relavent S2L2 representations of Stokes vectors. It will be helpful for users about to reproduce our S2L2 representation method." ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "execution": { "iopub.execute_input": "2025-04-09T08:14:49.186238Z", "iopub.status.busy": "2025-04-09T08:14:49.186238Z", "iopub.status.idle": "2025-04-09T08:14:49.387354Z", "shell.execute_reply": "2025-04-09T08:14:49.387354Z" } }, "outputs": [], "source": [ "from typing import Literal\n", "import math\n", "import numpy as np\n", "from scipy.spatial.transform import Rotation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "-----\n", "## 1. Spin-2 SH with Order $l=2$\n", "### 1.1 Implementation of Spin-2 SH Based on $\\theta\\phi$-Frame Field\n", "\n", "$$\n", "\\begin{align}\n", "{}_2Y_{2,-2}\\left(\\theta,\\phi\\right) &= \\frac18 \\sqrt\\frac5\\pi \\left(1+2\\cos\\theta+\\cos^2\\theta\\right)e^{-2i\\phi}\\\\\n", "{}_2Y_{2,-1}\\left(\\theta,\\phi\\right) &= \\frac14\\sqrt\\frac5\\pi\\sin\\theta\\left(-1-\\cos\\theta\\right)e^{-i\\phi}\\\\\n", "{}_2Y_{20}\\left(\\theta,\\phi\\right) &= \\frac14\\sqrt\\frac{15}{2\\pi}\\sin^2\\theta\\\\\n", "{}_2Y_{21}\\left(\\theta,\\phi\\right) &= \\frac14\\sqrt\\frac5\\pi\\sin\\theta\\left(-1+\\cos\\theta\\right)e^{i\\phi}\\\\\n", "{}_2Y_{22}\\left(\\theta,\\phi\\right) &= \\frac18 \\sqrt\\frac5\\pi \\left(1-2\\cos\\theta+\\cos^2\\theta\\right)e^{+2i\\phi}\n", "\\end{align},\n", "$$\n", "where ${}_2Y_{lm}\\coloneqq \\left[\\overset\\leftrightarrow Y_{lm}\\right]^{{\\vec{\\mathbf F}}_{\\theta\\phi}}$." ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "execution": { "iopub.execute_input": "2025-04-09T08:14:49.389358Z", "iopub.status.busy": "2025-04-09T08:14:49.389358Z", "iopub.status.idle": "2025-04-09T08:14:49.394251Z", "shell.execute_reply": "2025-04-09T08:14:49.394251Z" } }, "outputs": [], "source": [ "def S2L2_basis_wrtTP(theta, phi):\n", " \"\"\"\n", " Parameters\n", " ----------\n", " theta: ArrayLike[*]\n", " In the radian units. having a broadcastable shape with `phi`.\n", " phi: ArrayLike[*]\n", " In the radian units. having a broadcastable shape with `theta`.\n", "\n", " Return\n", " ------\n", " s2Y_2m: np.ndarray[*, 5], complex\n", " Complex values for spin-2 SH, with respect to the $\\theta\\phi$-frame field, with $l=2$ at given values of the spherical coordinates.\n", " `s2Y_2m[..., 0]` to `s2Y_2m[..., 4]` indicate $m=-2,\\cdots,2$, respectively.\n", " \"\"\"\n", " theta, phi = np.broadcast_arrays(theta, phi)\n", "\n", " cos = np.cos(theta)\n", " sin = np.sin(theta)\n", " cos_p1 = cos + 1\n", " cos_m1 = cos - 1\n", "\n", " s2Y_2m = np.full((*theta.shape, 5), 1/8*np.sqrt(5/np.pi), dtype=complex)\n", "\n", " s2Y_2m[..., 0] *= (cos_p1)**2\n", " s2Y_2m[..., 1] *= -2*sin*cos_p1\n", " s2Y_2m[..., 2] *= np.sqrt(6) * (sin**2)\n", " s2Y_2m[..., 3] *= 2*sin*cos_m1\n", " s2Y_2m[..., 4] *= (cos_m1)**2\n", " del cos_p1, cos_m1, sin\n", " \n", " m = np.arange(-2, 3)\n", " s2Y_2m[...] *= np.exp(m * 1j * phi[..., None])\n", " return s2Y_2m" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1.2. Implementations of Spin-2 SH at Given Frame\n", "In many cases, we already has our measurement frame (local frames) $\\vec{\\mathbf F}\\in\\vec{\\mathbb F}^3$ and then want Stokes parameters with respect to such frame $\\vec{\\mathbf F}$. Here we introduce two ways.\n", "\n", "#### (a) **Computation using Euler angles**\n", "\n", "The first way is utilizing ${}_2Y_{lm}\\left(\\theta,\\phi\\right)$, which indicates complex-values Stokes parameters $s_1+is_2$ with respect to the $\\theta\\phi$-frame field $\\vec{\\mathbf F}_{\\theta\\phi}$. Then we can convert it with respect to the frame $\\vec{\\mathbf F}$ using the angle between these two frames, i.e.,\n", "\n", "$$\n", "\\begin{align}\n", "{}_2Y_{lm}\\left(\\vec {\\mathbf F}\\right) &\\coloneqq \\left[\\overset\\leftrightarrow Y_{lm}\\left(\\vec{\\mathbf F}\\hat z_g\\right)\\right]^{{\\vec{\\mathbf F}}} = {}_2Y_{lm}\\left(\\theta,\\phi\\right) e^{-2i\\psi}, \\\\\n", "\\text{where } \\vec{\\mathbf F} &= \\vec {\\mathbf F}_{\\theta\\phi} \\begin{bmatrix}\n", " \\cos\\psi & -\\sin\\psi & 0 \\\\\n", " \\sin\\psi & \\cos\\psi & 0 \\\\\n", " 0 & 0 & 1\n", "\\end{bmatrix}.\n", "\\end{align}\n", "$$\n", "\n", "Note that $\\vec{\\mathbf F}_g=\\begin{bmatrix}\\hat x_g & \\hat y_g & \\hat z_g\\end{bmatrix}$ denotes the global (world) frame, so that $\\vec{\\mathbf F}\\hat z_g$ indicates the ray direction (local z-axis) of the measurement frame $\\vec {\\mathbf F}$ and $\\theta$, $\\phi$, and $\\psi$ can be obtained from ZYZ Euler angles of the rotation $\\vec{\\mathbf F}$ with respect to $\\vec{\\mathbf F}_g$, i.e.,\n", "\n", "$$\n", "\\vec{\\mathbf F} = \\vec R_{\\hat z_g \\hat y_g \\hat z_g}\\left(\\phi, \\theta,\\psi \\right) \\vec{\\mathbf F}_g.\n", "$$\n", "\n", "Since $\\vec{\\mathbf F}_g$ is the world frame, we can simply consider it a the 3D identity matrix in the numerical implementation.\n", "\n", "#### (b) **Computation using quaternions**\n", "However, there is more elegant implementation using quaternion representation $\\mathbf q = q_0+q_1i+q_2j+q_3k$ of the rotation $\\vec R$ such that $\\vec{\\mathbf F} = \\vec R \\vec{\\mathbf F}_g$. Simply denoting ${}_2Y_{lm}\\left(\\mathbf q\\right) = {}_2Y_{lm}\\left(\\vec{\\mathbf F}\\right)$, the following formula can be driven through some detailed steps. Here we use two auxiliary complex numbers $\\mathbf q_a = q_0 + q_3i$ and $\\mathbf q_b = q_2 + q_1i$.\n", "\n", "$$\n", "\\begin{align}\n", "\t{}_2Y_{2,-2}\\left({\\mathbf q}\\right) &= \\sqrt\\frac{5}{4\\pi} \\bar{\\mathbf q}_a^4 \\\\\n", "\t{}_2Y_{2,-1}\\left({\\mathbf q}\\right) &= -\\sqrt\\frac{5}{\\pi} \\bar{\\mathbf q}_a^3 \\bar{\\mathbf q}_b \\\\\n", "\t{}_2Y_{20}\\left({\\mathbf q}\\right) &= \\sqrt\\frac{15}{2\\pi}\\bar{\\mathbf q}_a^2 \\bar{\\mathbf q}_b^2 \\\\\n", "\t{}_2Y_{21}\\left({\\mathbf q}\\right) &= -\\sqrt\\frac{5}{\\pi} \\bar{\\mathbf q}_a \\bar{\\mathbf q}_b^3 \\\\\n", "\t{}_2Y_{22}\\left({\\mathbf q}\\right) &= \\sqrt\\frac{5}{4\\pi} \\bar{\\mathbf q}_b^{4} \\\\\n", "\\end{align}\n", "$$" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "execution": { "iopub.execute_input": "2025-04-09T08:14:49.395256Z", "iopub.status.busy": "2025-04-09T08:14:49.395256Z", "iopub.status.idle": "2025-04-09T08:14:49.401254Z", "shell.execute_reply": "2025-04-09T08:14:49.401254Z" } }, "outputs": [], "source": [ "def S2L2_basis_at_frame(F, method: Literal['euler', 'quat'] = 'quat'):\n", " \"\"\"\n", " Parameters\n", " ----------\n", " F: ArrayLike[*, 3, 3]\n", " An ND-array of local (measurement) frames. It can also be considered a local-to-world coordinate conversion matrix, \n", " i.e. `F[...,:,0]`, `F[...,:,1]`, `F[...,:,2]` indicate x-, y-, and z-axis vectors in the global (world) coordinates.\n", " On the other hand, `F[...,0,:]`, `F[...,1,:]`, `F[...,2,:]` indicate global x-, y-, and z-coordinates of each local axis vectors\n", " method: literal 'euler', 'quat'\n", " Implementation to evaluate spin-2 SH basis functions.\n", " * `'euler'`: It first calls the `S2L2_basis_wrtTP()` function, which returns values of the basis functions with respect to the theta-phi-frame field.\n", " Then find `psi` which is the angle between given local frame `F` and the theta-phi-frame, around their common z-axis.\n", " The final value becomes the output of `S2L2_basis_wrtTP()` multiplied by $e^{-2i\\psi}$.\n", " The name 'euler' of this method comes from the fact that `F` is equivalent to the rotation with ZYZ Euler angles phi, theta, psi.\n", " * `'quat'`: The values of the basis functions are evaluated directly from the quaternion formulation of `F`.\n", " \n", " Return\n", " ------\n", " s2Y_2m: np.ndarray[*, 5], complex\n", " Complex values for spin-2 SH, with respect to given local frames `F`, with $l=2$ at given values of the spherical coordinates.\n", " `s2Y_2m[..., 0]` to `s2Y_2m[..., 4]` indicate $m=-2,\\cdots,2$, respectively.\n", " \"\"\"\n", " F = np.asarray(F)\n", " shape_chan = F.shape[:-2]\n", " assert F.shape[-2:] == (3, 3)\n", " R = Rotation.from_matrix(F.reshape(-1, 3, 3))\n", " if method == 'euler':\n", " phi, theta, psi = np.moveaxis(R.as_euler('ZYZ', degrees=False), -1, 0)\n", " s2Y_2m = S2L2_basis_wrtTP(theta, phi)\n", " s2Y_2m *= np.exp(-2j*psi[..., None])\n", " \n", " else:\n", " assert method == 'quat'\n", " # [WARNING] Scipy uses scalar-last order at default.\n", " # https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.transform.Rotation.as_quat.html\n", " q = R.as_quat()\n", " qa_conj = q[..., 3] - 1j*q[..., 2]\n", " qb_conj = q[..., 1] - 1j*q[..., 0]\n", " \n", " s2Y_2m = np.full((math.prod(shape_chan), 5), np.sqrt(5/np.pi), complex)\n", " s2Y_2m[..., 0] *= 0.5 * qa_conj**4\n", " s2Y_2m[..., 1] *= -qa_conj**3 * qb_conj\n", " s2Y_2m[..., 2] *= np.sqrt(3/2) * (qa_conj * qb_conj)**2\n", " s2Y_2m[..., 3] *= -qa_conj * qb_conj**3\n", " s2Y_2m[..., 4] *= 0.5 * qb_conj**4\n", " \n", " return s2Y_2m.reshape(*shape_chan, 5)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1.3. Validation of Consistency of Two Methods" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "execution": { "iopub.execute_input": "2025-04-09T08:14:49.402418Z", "iopub.status.busy": "2025-04-09T08:14:49.402418Z", "iopub.status.idle": "2025-04-09T08:14:49.421572Z", "shell.execute_reply": "2025-04-09T08:14:49.421572Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Clear!\n" ] } ], "source": [ "theta, phi = np.mgrid[:np.pi:20j, :2*np.pi:40j]\n", "psi = np.linspace(0, 2*np.pi, 10)[:, None, None]\n", "angles = np.stack(np.broadcast_arrays(phi, theta, psi), -1).reshape(-1, 3)\n", "F = Rotation.from_euler('ZYZ', angles).as_matrix()\n", "\n", "s2Y_2m_euler = S2L2_basis_at_frame(F, 'euler')\n", "s2Y_2m_quat = S2L2_basis_at_frame(F, 'quat')\n", "assert np.allclose(s2Y_2m_euler, s2Y_2m_quat)\n", "print(\"Clear!\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "-----\n", "## 2. S2L2 Representation of Stokes Vectors\n", "\n", "Given Stokes vector $\\overset\\leftrightarrow s = \\left[\\mathbf s\\right]_{\\vec{\\mathbf F}}$, its 12-parameter S2L2 representation $\\mathbf r = \\mathrm{S2L2}\\left(\\overset\\leftrightarrow s\\right) \\in \\mathbb R^{12}$ is characterized as:\n", "\n", "$$\n", "\\begin{align}\n", " r_0 &= s_0, \\\\\n", " \\begin{bmatrix}\n", " r_1 \\\\ r_2\n", " \\end{bmatrix} &= \\sqrt\\frac{4\\pi}5 \\mathbb R^2 \\left( {}_2Y_{2,-2}\\left(\\vec{\\mathbf F}\\right)^*\\left(s_1+is_2\\right) \\right), \\\\\n", " \\begin{bmatrix}\n", " r_3 \\\\ r_4\n", " \\end{bmatrix} &= \\sqrt\\frac{4\\pi}5 \\mathbb R^2 \\left( {}_2Y_{2,-1}\\left(\\vec{\\mathbf F}\\right)^*\\left(s_1+is_2\\right) \\right), \\\\\n", " \\begin{bmatrix}\n", " r_5 \\\\ r_6\n", " \\end{bmatrix} &= \\sqrt\\frac{4\\pi}5 \\mathbb R^2 \\left( {}_2Y_{20}\\left(\\vec{\\mathbf F}\\right)^*\\left(s_1+is_2\\right) \\right), \\\\\n", " \\begin{bmatrix}\n", " r_7 \\\\ r_8\n", " \\end{bmatrix} &= \\sqrt\\frac{4\\pi}5 \\mathbb R^2 \\left( {}_2Y_{21}\\left(\\vec{\\mathbf F}\\right)^*\\left(s_1+is_2\\right) \\right), \\\\\n", " \\begin{bmatrix}\n", " r_9 \\\\ r_{10}\n", " \\end{bmatrix} &= \\sqrt\\frac{4\\pi}5 \\mathbb R^2 \\left( {}_2Y_{22}\\left(\\vec{\\mathbf F}\\right)^*\\left(s_1+is_2\\right) \\right), \\\\\n", " r_{11} &= s_3,\n", "\\end{align}\n", "$$\n", "\n", "where $\\mathbb R^2$ denotes conversion from complex numbers to 2D numeric real vectors as $\\mathbb R^2\\left(x+yi\\right)=\\left[x,y\\right]^T$.\n" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "execution": { "iopub.execute_input": "2025-04-09T08:14:49.469709Z", "iopub.status.busy": "2025-04-09T08:14:49.469709Z", "iopub.status.idle": "2025-04-09T08:14:49.475083Z", "shell.execute_reply": "2025-04-09T08:14:49.475083Z" } }, "outputs": [], "source": [ "def stk_to_S2L2repr(stkComp, F, method: Literal['euler', 'quat'] = 'quat'):\n", " \"\"\"\n", " Parameters\n", " ----------\n", " stkComp: ArrayLike[*, 4|2]\n", " Numerical vectors of Stokes parameters, which should be considered with the measurement frame `F`.\n", " F: ArrayLike[*, 3, 3]\n", " Measurement frames\n", " method: literal 'euler', 'quat'\n", " Implementation to evaluate spin-2 SH basis functions.\n", " * `'euler'`: It first calls the `S2L2_basis_wrtTP()` function, which returns values of the basis functions with respect to the theta-phi-frame field.\n", " Then find `psi` which is the angle between given local frame `F` and the theta-phi-frame, around their common z-axis.\n", " The final value becomes the output of `S2L2_basis_wrtTP()` multiplied by $e^{-2i\\psi}$.\n", " The name 'euler' of this method comes from the fact that `F` is equivalent to the rotation with ZYZ Euler angles phi, theta, psi.\n", " * `'quat'`: The values of the basis functions are evaluated directly from the quaternion formulation of `F`.\n", " \n", " Return\n", " ------\n", " res: np.ndarray[*, 12|10 or 7|5]\n", " S2L2 representation. It consists of 12 or 7 parameters for each Stokes vector.\n", " If spin-2 Stokes vector (only consists of $s_1$ and $s_$ Stokes parameters) was given,\n", " the number of output parameters changes to 10 for `type==\"COMP\"` and 5 for `type==\"REAL\"`, respectively.\n", " \"\"\"\n", " p = stkComp.shape[-1]\n", " assert p in [2, 4], f\"{stkComp.shape = }\"\n", " assert F.shape[-2:] == (3, 3), f\"{F.shape = }\"\n", " shape_chan = np.broadcast_shapes(stkComp.shape[:-1], F.shape[:-2])\n", " shape = (*shape_chan, p+8) # [*, 12] or [*, 10]\n", " res = np.empty(shape)\n", "\n", " # ---------- Spin-0 ----------\n", " if p == 4:\n", " res[..., 0] = stkComp[..., 0]\n", " res[..., -1] = stkComp[..., 3]\n", " s1, s2 = stkComp[..., 1], stkComp[..., 2]\n", " idx_s12 = np.s_[1:-1]\n", " else:\n", " s1, s2 = stkComp[..., 0], stkComp[..., 1]\n", " idx_s12 = np.s_[:]\n", " # ---------- Spin-2 ----------\n", " s1_is2 = s1 + 1j*s2 # shape [*]\n", " s2Y_2m = S2L2_basis_at_frame(F, method) # shape [*, 5]\n", " prod = s2Y_2m.conj() * s1_is2[..., None] # shape [*, 5]\n", "\n", " res[..., idx_s12] = np.stack([prod.real, prod.imag], -1).reshape(*shape_chan, 10)\n", " res[..., idx_s12] *= 2*np.sqrt(np.pi/5)\n", " return res" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2.1. Properties\n", "S2L2 representation is a frame-free method, i.e., two pairs of `(stkComp, F)` arguments under the relationship of Stokes coordiante conversion produce consistent S2L2 representaion. This representation is also linear for a fixed ray direction.\n", "\n", "\n", "\n", "The constant $\\sqrt\\frac{4\\pi}5$ in the above S2L2 equations are intended to preseve norms. Norms, degrees of polarization, unpolarized intensities can be computed from S2L2 representations." ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "execution": { "iopub.execute_input": "2025-04-09T08:14:49.476130Z", "iopub.status.busy": "2025-04-09T08:14:49.476130Z", "iopub.status.idle": "2025-04-09T08:14:49.565897Z", "shell.execute_reply": "2025-04-09T08:14:49.565897Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Well defined independent of local frames!\n", "Linear!\n", "Norms are preserved!\n", "Degrees of polarization are preserved!\n", "Unpolarized intensity is preserved!\n", "\n" ] } ], "source": [ "# ---------- Utils ----------\n", "def rotation_z(angle):\n", " return np.array([\n", " [np.cos(angle), -np.sin(angle), 0],\n", " [np.sin(angle), np.cos(angle), 0],\n", " [0, 0, 1]\n", " ])\n", "\n", "stk = np.random.rand(*F.shape[:-2], 4)\n", "\n", "# ---------- S2L2 representation ----------\n", "s2l2 = stk_to_S2L2repr(stk, F)\n", "\n", "# ---------- Invariant under Stokes coordinates conversion ----------\n", "for i in range(10):\n", " psi = 2*np.pi*(i+1)/11\n", " F2 = F @ rotation_z(psi)\n", " stk2 = stk.copy()\n", " stk2[..., 1] = np.cos(2*psi) * stk[..., 1] + np.sin(2*psi) * stk[..., 2]\n", " stk2[..., 2] = -np.sin(2*psi) * stk[..., 1] + np.cos(2*psi) * stk[..., 2]\n", " assert np.allclose(s2l2, stk_to_S2L2repr(stk2, F2))\n", "print(\"Well defined independent of local frames!\")\n", "\n", "# ---------- Linearity for Stokes vectors with the same ray direction ----------\n", "stk2 = np.random.rand(*F.shape[:-2], 4)\n", "s2l2_2 = stk_to_S2L2repr(stk2, F)\n", "assert np.allclose(s2l2 + s2l2_2, stk_to_S2L2repr(stk + stk2, F))\n", "print(\"Linear!\")\n", "\n", "# ---------- Other properties related to norms ----------\n", "assert np.allclose(np.linalg.norm(stk, 2, -1), np.linalg.norm(s2l2, 2, -1))\n", "print(\"Norms are preserved!\")\n", "\n", "polar_s0123 = np.sqrt((stk[..., 1:]**2).sum(-1))\n", "polar_s2l2 = np.sqrt((s2l2[..., 1:]**2).sum(-1))\n", "\n", "dop_s0123 = polar_s0123 / stk[..., 0]\n", "dop_s2l2 = polar_s2l2 / s2l2[..., 0]\n", "assert np.allclose(dop_s0123, dop_s2l2)\n", "print(\"Degrees of polarization are preserved!\")\n", "\n", "assert np.allclose(stk[..., 0] - polar_s0123,\n", " s2l2[..., 0] - polar_s2l2)\n", "print(\"Unpolarized intensity is preserved!\\n\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2.2. Reconstructing Stokes Parameters from S2L2 Representations\n", "$\\mathbf s = \\mathrm{S2L2Inv}\\left(\\mathbf r; \\vec{\\mathbf F}\\right)\\in{\\mathbb R}^{4}$\n", "\n", "$$\n", "\\begin{align}\n", " s_0 &= r_0 \\\\\n", " \\begin{bmatrix}\n", " s_1 \\\\ s_2\n", " \\end{bmatrix} &= \\sqrt\\frac{4\\pi}{5}\\mathbb R^2\\left(\\sum_{m=-2}^2 {}_2Y_{2m}\\left(\\vec{\\mathbf F}\\right)\\left(r_{2m+5}+ir_{2m+6}\\right)\\right) \\\\\n", " s_3 &= r_{11} \\\\\n", "\\end{align}\n", "$$\n" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "execution": { "iopub.execute_input": "2025-04-09T08:14:49.566902Z", "iopub.status.busy": "2025-04-09T08:14:49.566902Z", "iopub.status.idle": "2025-04-09T08:14:49.575901Z", "shell.execute_reply": "2025-04-09T08:14:49.575901Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The inversion is correct!\n" ] } ], "source": [ "def S2L2repr_to_stk(s2l2, F, method: Literal['euler', 'quat'] = 'quat'):\n", " \"\"\"\n", " Parameters\n", " ----------\n", " s2l2: ArrayLike[*, 12|10]\n", " S2L2 representation. It consists of 12 parameters for each Stokes vector.\n", " F: ArrayLike[*, 3, 3]\n", " Measurement frames\n", " method: literal 'euler', 'quat'\n", " Implementation to evaluate spin-2 SH basis functions.\n", " * `'euler'`: It first calls the `S2L2_basis_wrtTP()` function, which returns values of the basis functions with respect to the theta-phi-frame field.\n", " Then find `psi` which is the angle between given local frame `F` and the theta-phi-frame, around their common z-axis.\n", " The final value becomes the output of `S2L2_basis_wrtTP()` multiplied by $e^{-2i\\psi}$.\n", " The name 'euler' of this method comes from the fact that `F` is equivalent to the rotation with ZYZ Euler angles phi, theta, psi.\n", " * `'quat'`: The values of the basis functions are evaluated directly from the quaternion formulation of `F`.\n", " \n", " Return\n", " ------\n", " stk: np.ndarray[*, 4]\n", " Numerical vectors of Stokes parameters, which should be considered with the measurement frame `F`.\n", " \"\"\"\n", " F = np.asarray(F)\n", " s2l2 = np.asarray(s2l2)\n", "\n", " assert F.shape[-2:] == (3, 3), f\"{F.shape = }\"\n", " n_param = s2l2.shape[-1]\n", " assert n_param in [10, 12], f\"Invalid shape: {s2l2.shape = }\"\n", " if n_param % 5 == 0:\n", " n_stk = 2 # Spin-2 Stokes vector only\n", " idx_s12 = np.s_[:]\n", " idx_s12_p1 = np.s_[0::2]\n", " idx_s12_p2 = np.s_[1::2]\n", " oidx_s1, oidx_s2 = 0, 1\n", " else:\n", " n_stk = 4 # Full Stokes vector\n", " idx_s12 = np.s_[1:-1]\n", " idx_s12_p1 = np.s_[1:-1:2]\n", " idx_s12_p2 = np.s_[2:-1:2]\n", " oidx_s1, oidx_s2 = 1, 2\n", " \n", " shape_chan = np.broadcast_shapes(s2l2.shape[:-1], F.shape[:-2])\n", " shape_stk = (*shape_chan, n_stk)\n", " stk = np.empty(shape_stk)\n", "\n", " # ---------- Spin-0 ----------\n", " if n_stk == 4:\n", " stk[..., 0] = s2l2[..., 0]\n", " stk[..., -1] = s2l2[..., -1]\n", "\n", " # ---------- Spin-2 ----------\n", " s2l2_5 = s2l2[..., idx_s12_p1] + s2l2[..., idx_s12_p2]*1j # shape [*, 5] complex\n", " factor = np.sqrt(4*np.pi/5)\n", " \n", " s2Y_2m = S2L2_basis_at_frame(F, method) # shape [*, 5] complex\n", " dot = (s2Y_2m * s2l2_5).sum(-1) # shape [*] complex\n", "\n", " stk[..., oidx_s1] = factor * dot.real\n", " stk[..., oidx_s2] = factor * dot.imag\n", " return stk\n", "\n", "assert np.allclose(stk, S2L2repr_to_stk(s2l2, F))\n", "print(\"The inversion is correct!\")" ] } ], "metadata": { "kernelspec": { "display_name": "polarsh", "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.11.7" } }, "nbformat": 4, "nbformat_minor": 2 }