{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import sympy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Notes on `PCi_j` Formalism for Overlappograms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the `PCi_j` formalism, the conversion between intermediate world coordinates and pixel coordinates is given by (from Greisen and Calabretta, 2002),\n",
    "\n",
    "$$\n",
    "x_i = s_iq_i\n",
    "$$\n",
    "$$\n",
    "q_i = m_{ij}(p_j - r_j)\n",
    "$$\n",
    "\n",
    "or in more explicit matrix notation,\n",
    "\n",
    "$$\n",
    "X = sQ\n",
    "$$\n",
    "$$\n",
    "Q = MP\n",
    "$$\n",
    "\n",
    "In other words, we transform our pixel coordinates $P$ to our intermediate world coordinates $Q$ through the transformation $M$. \n",
    "This matrix $M$ is typically referred to as the `PCi_j` matrix, where $i$ indexes the world coordinates and $j$ the pixel coordinates.\n",
    "$s$ is a scalar that represents the `CDELT` conversion factor between pixel units and world units.\n",
    "\n",
    "The resulting `PC_ij` matrix for our overlappogram **assuming that the dispersion and latitude world axes are aligned with the y-like pixel axis, $p_2$**,\n",
    "\n",
    "$$\n",
    "D(\\mu) = \\begin{bmatrix}\n",
    "1 & 0 & 0 \\\\\n",
    "0 & 1 & -\\mu \\\\\n",
    "0 & 0 & 1\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "where $\\mu$ is the order of the dispersion.\n",
    "\n",
    "In general, the dispersion axis may be oriented at some angle $\\gamma$ relative to $p_2$.\n",
    "Thus, before applying $D$, we rotate $P$ by $\\gamma$ into $P^\\prime$ such that $P^\\prime=R(\\gamma)P$ where $R$ is the usual rotation matrix,\n",
    "\n",
    "$$\n",
    "R(\\theta) = \\begin{bmatrix}\n",
    "\\cos\\theta & -\\sin\\theta & 0 \\\\\n",
    "\\sin\\theta & \\cos\\theta & 0 \\\\\n",
    "0 & 0 & 1\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "Lastly, we apply a rotation to align $P^\\prime$ with the intermediate world axes $Q$ such that $Q=R(\\beta)P^\\prime$.\n",
    "Note that this angle $\\beta$ is the relative orientation between $P^\\prime$ and $Q$. \n",
    "However, we will typically be measuring the satellite roll angle relative to the original pixel axes $P$.\n",
    "Thus, $\\beta=\\alpha-\\gamma$ where $\\alpha$ is the orientation of solar north relative to $p_1$.\n",
    "\n",
    "The full `PC` matrix can thus be written as,\n",
    "\n",
    "$$\n",
    "M = R(\\alpha-\\gamma)D(\\mu)R(\\gamma)\n",
    "$$\n",
    "\n",
    "which can be written explicitly as,\n",
    "\n",
    "$$\n",
    "M = \\begin{bmatrix}\n",
    "\\cos{\\left(\\alpha \\right)} & -\\sin{\\left(\\alpha \\right)} & - \\mu\\cos{\\left(\\alpha - \\gamma \\right)}\\\\\n",
    "\\sin{\\left(\\alpha \\right)} & \\cos{\\left(\\alpha \\right)} & - \\mu\\sin{\\left(\\alpha - \\gamma \\right)}\\\\\n",
    "0 & 0 & 1\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "We can evaluate this matrix product with `sympy` and look at a few example cases."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "mu = sympy.symbols('mu')\n",
    "gamma = sympy.symbols('gamma')\n",
    "alpha = sympy.symbols('alpha')\n",
    "theta = sympy.symbols('theta')\n",
    "phi = sympy.symbols('phi')\n",
    "rotation = sympy.Matrix([\n",
    "    [sympy.cos(theta), -sympy.sin(theta), 0],\n",
    "    [sympy.sin(theta), sympy.cos(theta), 0],\n",
    "    [0, 0, 1]\n",
    "])\n",
    "dispersion = sympy.Matrix([\n",
    "    [1, 0, -mu],\n",
    "    [0, 1, 0],\n",
    "    [0 ,0, 1],\n",
    "])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "⎡cos(α)  -sin(α)  -μ⋅cos(α - γ)⎤\n",
      "⎢                              ⎥\n",
      "⎢sin(α)  cos(α)   -μ⋅sin(α - γ)⎥\n",
      "⎢                              ⎥\n",
      "⎣  0        0           1      ⎦\n"
     ]
    }
   ],
   "source": [
    "PC_matrix = rotation.subs(theta, (alpha-gamma)) * dispersion * rotation.subs(theta, gamma)\n",
    "PC_matrix.simplify()\n",
    "sympy.printing.pprint(PC_matrix)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can evaluate the `PCi_j` matrix at a few different combinations of $\\alpha,\\gamma,\\mu$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "alpha=0, gamma=0, mu=0\n",
      "-----------------------------\n",
      "⎡1  0  0⎤\n",
      "⎢       ⎥\n",
      "⎢0  1  0⎥\n",
      "⎢       ⎥\n",
      "⎣0  0  1⎦\n",
      "alpha=0, gamma=0, mu=1\n",
      "-----------------------------\n",
      "⎡1  0  -1⎤\n",
      "⎢        ⎥\n",
      "⎢0  1  0 ⎥\n",
      "⎢        ⎥\n",
      "⎣0  0  1 ⎦\n",
      "alpha=pi/2, gamma=0, mu=1\n",
      "-----------------------------\n",
      "⎡0  -1  0 ⎤\n",
      "⎢         ⎥\n",
      "⎢1  0   -1⎥\n",
      "⎢         ⎥\n",
      "⎣0  0   1 ⎦\n",
      "alpha=0, gamma=pi/2, mu=1\n",
      "-----------------------------\n",
      "⎡1  0  0⎤\n",
      "⎢       ⎥\n",
      "⎢0  1  1⎥\n",
      "⎢       ⎥\n",
      "⎣0  0  1⎦\n",
      "alpha=pi/2, gamma=pi/2, mu=1\n",
      "-----------------------------\n",
      "⎡0  -1  -1⎤\n",
      "⎢         ⎥\n",
      "⎢1  0   0 ⎥\n",
      "⎢         ⎥\n",
      "⎣0  0   1 ⎦\n",
      "alpha=pi/4, gamma=0, mu=1\n",
      "-----------------------------\n",
      "⎡√2  -√2   -√2 ⎤\n",
      "⎢──  ────  ────⎥\n",
      "⎢2    2     2  ⎥\n",
      "⎢              ⎥\n",
      "⎢√2   √2   -√2 ⎥\n",
      "⎢──   ──   ────⎥\n",
      "⎢2    2     2  ⎥\n",
      "⎢              ⎥\n",
      "⎣0    0     1  ⎦\n",
      "alpha=0, gamma=pi/4, mu=1\n",
      "-----------------------------\n",
      "⎡      -√2 ⎤\n",
      "⎢1  0  ────⎥\n",
      "⎢       2  ⎥\n",
      "⎢          ⎥\n",
      "⎢       √2 ⎥\n",
      "⎢0  1   ── ⎥\n",
      "⎢       2  ⎥\n",
      "⎢          ⎥\n",
      "⎣0  0   1  ⎦\n",
      "alpha=pi/4, gamma=pi/4, mu=1\n",
      "-----------------------------\n",
      "⎡√2  -√2     ⎤\n",
      "⎢──  ────  -1⎥\n",
      "⎢2    2      ⎥\n",
      "⎢            ⎥\n",
      "⎢√2   √2     ⎥\n",
      "⎢──   ──   0 ⎥\n",
      "⎢2    2      ⎥\n",
      "⎢            ⎥\n",
      "⎣0    0    1 ⎦\n"
     ]
    }
   ],
   "source": [
    "for a,g,m in [(0,0,0),\n",
    "              (0,0,1),\n",
    "              (sympy.pi/2,0,1),\n",
    "              (0,sympy.pi/2,1),\n",
    "              (sympy.pi/2,sympy.pi/2,1),\n",
    "              (sympy.pi/4,0,1),\n",
    "              (0,sympy.pi/4,1),\n",
    "              (sympy.pi/4,sympy.pi/4,1)]:\n",
    "    print(f'alpha={a}, gamma={g}, mu={m}')\n",
    "    print('-----------------------------')\n",
    "    sympy.printing.pprint(PC_matrix.subs([(alpha, a),(gamma, g),(mu, m)]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can do some simple computations that shows how this transforms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "p = sympy.symarray('p', 3)\n",
    "r = sympy.symarray('r', 3)\n",
    "pix_coord = sympy.Matrix([\n",
    "    p[0] - r[0],\n",
    "    p[1] - r[1],\n",
    "    p[2] - r[2],\n",
    "])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "⎡p₀ - r₀⎤\n",
      "⎢       ⎥\n",
      "⎢p₁ - r₁⎥\n",
      "⎢       ⎥\n",
      "⎣p₂ - r₂⎦\n"
     ]
    }
   ],
   "source": [
    "sympy.printing.pprint(pix_coord)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "world_coord = PC_matrix * pix_coord"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "⎡-μ⋅(p₂ - r₂)⋅cos(α - γ) + (p₀ - r₀)⋅cos(α) - (p₁ - r₁)⋅sin(α)⎤\n",
      "⎢                                                             ⎥\n",
      "⎢-μ⋅(p₂ - r₂)⋅sin(α - γ) + (p₀ - r₀)⋅sin(α) + (p₁ - r₁)⋅cos(α)⎥\n",
      "⎢                                                             ⎥\n",
      "⎣                           p₂ - r₂                           ⎦\n"
     ]
    }
   ],
   "source": [
    "sympy.printing.pprint(world_coord)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "data_shape = (1073, 2000, 750)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "⎡p₀ - p₂ - 374.5⎤\n",
      "⎢               ⎥\n",
      "⎢  p₁ - 1000.5  ⎥\n",
      "⎢               ⎥\n",
      "⎣    p₂ - 1     ⎦\n"
     ]
    }
   ],
   "source": [
    "Q = world_coord.subs([\n",
    "    (alpha, 0),\n",
    "    (gamma, 0),\n",
    "    (mu, 1),\n",
    "    (r[0], (data_shape[2] + 1)/2),\n",
    "    (r[1], (data_shape[1] + 1)/2),\n",
    "    (r[2], 1),#(data_shape[0] + 1)/2),\n",
    "])\n",
    "sympy.printing.pprint(Q)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "⎡-99.0⎤\n",
      "⎢     ⎥\n",
      "⎢  0  ⎥\n",
      "⎢     ⎥\n",
      "⎣ 99  ⎦\n"
     ]
    }
   ],
   "source": [
    "sympy.printing.pprint(Q.subs([\n",
    "    (p[0], (data_shape[2] + 1)/2),\n",
    "    (p[1], (data_shape[1] + 1)/2),\n",
    "    (p[2], 100),\n",
    "]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "483d3e286f47377cd75b0b13b7cda02b28091785fef97b0f4cde78d6f957a4e0"
  },
  "kernelspec": {
   "display_name": "Python [conda env:mocksipipeline]",
   "language": "python",
   "name": "conda-env-mocksipipeline-py"
  },
  "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.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}