{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ff5728c5",
   "metadata": {},
   "source": [
    "# 静态反射透射系数矩阵详细表达式\n",
    "\n",
    "+ Author: Zhu Dengda  \n",
    "+ Email:  zhudengda@mail.iggcas.ac.cn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "c62f2eb3",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "np.set_printoptions(linewidth=200)\n",
    "\n",
    "import sympy as sp\n",
    "from sympy.printing.latex import latex\n",
    "from IPython.display import display, Math"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "67b0168a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 定义基本变量\n",
    "k, d = sp.symbols(r'k d')\n",
    "mu1, Delta1 = sp.symbols(r'\\mu_1 \\Delta_1')\n",
    "mu2, Delta2 = sp.symbols(r'\\mu_2 \\Delta_2')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "a892ed51",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{2 \\Delta_{1}}{2 \\Delta_{1} + 2} + \\frac{4 \\mu_{2} k}{4 \\Delta_{1} \\mu_{1} k + 4 \\mu_{1} k} & \\frac{8 \\Delta_{1} \\Delta_{2} \\mu_{2} d k^{2}}{4 \\Delta_{1} \\mu_{1} k + 4 \\mu_{1} k} + \\frac{4 \\Delta_{1} d k}{2 \\Delta_{1} + 2} & \\frac{8 \\Delta_{1} \\mu_{2} d k^{2}}{4 \\Delta_{1} \\mu_{1} k + 4 \\mu_{1} k} - \\frac{4 \\Delta_{1} d k}{2 \\Delta_{1} + 2} & - \\frac{2 \\Delta_{1}}{2 \\Delta_{1} + 2} + \\frac{4 \\Delta_{2} \\mu_{2} k}{4 \\Delta_{1} \\mu_{1} k + 4 \\mu_{1} k}\\\\0 & \\frac{4 \\Delta_{2} \\mu_{2} k}{4 \\Delta_{1} \\mu_{1} k + 4 \\mu_{1} k} + \\frac{2}{2 \\Delta_{1} + 2} & \\frac{4 \\mu_{2} k}{4 \\Delta_{1} \\mu_{1} k + 4 \\mu_{1} k} - \\frac{2}{2 \\Delta_{1} + 2} & 0\\\\- \\frac{8 \\Delta_{1} \\mu_{2} d k^{2}}{4 \\Delta_{1} \\mu_{1} k + 4 \\mu_{1} k} + \\frac{4 \\Delta_{1} d k}{2 \\Delta_{1} + 2} & - \\frac{2 \\Delta_{1}}{2 \\Delta_{1} + 2} + \\frac{4 \\Delta_{2} \\mu_{2} k}{4 \\Delta_{1} \\mu_{1} k + 4 \\mu_{1} k} & \\frac{2 \\Delta_{1}}{2 \\Delta_{1} + 2} + \\frac{4 \\mu_{2} k}{4 \\Delta_{1} \\mu_{1} k + 4 \\mu_{1} k} & - \\frac{8 \\Delta_{1} \\Delta_{2} \\mu_{2} d k^{2}}{4 \\Delta_{1} \\mu_{1} k + 4 \\mu_{1} k} - \\frac{4 \\Delta_{1} d k}{2 \\Delta_{1} + 2}\\\\\\frac{4 \\mu_{2} k}{4 \\Delta_{1} \\mu_{1} k + 4 \\mu_{1} k} - \\frac{2}{2 \\Delta_{1} + 2} & 0 & 0 & \\frac{4 \\Delta_{2} \\mu_{2} k}{4 \\Delta_{1} \\mu_{1} k + 4 \\mu_{1} k} + \\frac{2}{2 \\Delta_{1} + 2}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[                   2*\\Delta_1/(2*\\Delta_1 + 2) + 4*\\mu_2*k/(4*\\Delta_1*\\mu_1*k + 4*\\mu_1*k), 8*\\Delta_1*\\Delta_2*\\mu_2*d*k**2/(4*\\Delta_1*\\mu_1*k + 4*\\mu_1*k) + 4*\\Delta_1*d*k/(2*\\Delta_1 + 2), 8*\\Delta_1*\\mu_2*d*k**2/(4*\\Delta_1*\\mu_1*k + 4*\\mu_1*k) - 4*\\Delta_1*d*k/(2*\\Delta_1 + 2),                   -2*\\Delta_1/(2*\\Delta_1 + 2) + 4*\\Delta_2*\\mu_2*k/(4*\\Delta_1*\\mu_1*k + 4*\\mu_1*k)],\n",
       "[                                                                                          0,                            4*\\Delta_2*\\mu_2*k/(4*\\Delta_1*\\mu_1*k + 4*\\mu_1*k) + 2/(2*\\Delta_1 + 2),                            4*\\mu_2*k/(4*\\Delta_1*\\mu_1*k + 4*\\mu_1*k) - 2/(2*\\Delta_1 + 2),                                                                                                    0],\n",
       "[-8*\\Delta_1*\\mu_2*d*k**2/(4*\\Delta_1*\\mu_1*k + 4*\\mu_1*k) + 4*\\Delta_1*d*k/(2*\\Delta_1 + 2),                  -2*\\Delta_1/(2*\\Delta_1 + 2) + 4*\\Delta_2*\\mu_2*k/(4*\\Delta_1*\\mu_1*k + 4*\\mu_1*k),                   2*\\Delta_1/(2*\\Delta_1 + 2) + 4*\\mu_2*k/(4*\\Delta_1*\\mu_1*k + 4*\\mu_1*k), -8*\\Delta_1*\\Delta_2*\\mu_2*d*k**2/(4*\\Delta_1*\\mu_1*k + 4*\\mu_1*k) - 4*\\Delta_1*d*k/(2*\\Delta_1 + 2)],\n",
       "[                            4*\\mu_2*k/(4*\\Delta_1*\\mu_1*k + 4*\\mu_1*k) - 2/(2*\\Delta_1 + 2),                                                                                                   0,                                                                                          0,                             4*\\Delta_2*\\mu_2*k/(4*\\Delta_1*\\mu_1*k + 4*\\mu_1*k) + 2/(2*\\Delta_1 + 2)]])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 根据 Q 矩阵的定义： Q(z2-, z2+) = D1(z2-)^(-1) * D2(z2+)\n",
    "# 因此在 D1(z2-) 中， (z-zj) 项即为该层厚度\n",
    "# 而在 D2(z2+) 中，(z-zj)项为 0\n",
    "Ds1 = sp.Matrix([\n",
    "    [1,         -(1 + 2*k*Delta1*d),            1,            -(1 - 2*k*Delta1*d)],\n",
    "    [1,          (1 - 2*k*Delta1*d),           -1,            -(1 + 2*k*Delta1*d)],\n",
    "    [2*mu1*k,     2*mu1*Delta1*k*(1 - 2*k*d),   2*mu1*k,      2*mu1*Delta1*k*(1 + 2*k*d)],\n",
    "    [2*mu1*k,    -2*mu1*Delta1*k*(1 + 2*k*d),  -2*mu1*k,      2*mu1*Delta1*k*(1 - 2*k*d)],\n",
    "])\n",
    "\n",
    "z = 0\n",
    "Ds2 = sp.Matrix([\n",
    "    [1,         -(1 + 2*k*Delta2*z),            1,            -(1 - 2*k*Delta2*z)],\n",
    "    [1,          (1 - 2*k*Delta2*z),           -1,            -(1 + 2*k*Delta2*z)],\n",
    "    [2*mu2*k,     2*mu2*Delta2*k*(1 - 2*k*z),   2*mu2*k,      2*mu2*Delta2*k*(1 + 2*k*z)],\n",
    "    [2*mu2*k,    -2*mu2*Delta2*k*(1 + 2*k*z),  -2*mu2*k,      2*mu2*Delta2*k*(1 - 2*k*z)],\n",
    "])\n",
    "\n",
    "Q = sp.expand(Ds1.inv() * Ds2)\n",
    "Q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "b591940e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{2 \\Delta_{1} d k \\left(- \\mu_{1} + \\mu_{2}\\right)}{\\Delta_{1} \\mu_{1} + \\mu_{2}} & \\frac{- 4 \\Delta_{1}^{2} d^{2} k^{2} \\left(\\mu_{1} - \\mu_{2}\\right) \\left(\\Delta_{2} \\mu_{2} + \\mu_{1}\\right) - \\left(\\Delta_{1} \\mu_{1} + \\mu_{2}\\right) \\left(\\Delta_{1} \\mu_{1} - \\Delta_{2} \\mu_{2}\\right)}{\\left(\\Delta_{1} \\mu_{1} + \\mu_{2}\\right) \\left(\\Delta_{2} \\mu_{2} + \\mu_{1}\\right)}\\\\\\frac{- \\mu_{1} + \\mu_{2}}{\\Delta_{1} \\mu_{1} + \\mu_{2}} & \\frac{2 \\Delta_{1} d k \\left(- \\mu_{1} + \\mu_{2}\\right)}{\\Delta_{1} \\mu_{1} + \\mu_{2}}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0 & \\frac{\\Delta_{1} \\mu_{1} - \\Delta_{2} \\mu_{2}}{\\Delta_{1} \\mu_{1} + \\mu_{2}}\\\\\\frac{\\mu_{1} - \\mu_{2}}{\\Delta_{2} \\mu_{2} + \\mu_{1}} & 0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{\\mu_{1} \\left(\\Delta_{1} + 1\\right)}{\\Delta_{1} \\mu_{1} + \\mu_{2}} & \\frac{2 \\Delta_{1} \\mu_{1} d k \\left(\\Delta_{1} + 1\\right)}{\\Delta_{1} \\mu_{1} + \\mu_{2}}\\\\0 & \\frac{\\mu_{1} \\left(\\Delta_{1} + 1\\right)}{\\Delta_{2} \\mu_{2} + \\mu_{1}}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{\\mu_{2} \\left(\\Delta_{2} + 1\\right)}{\\Delta_{2} \\mu_{2} + \\mu_{1}} & \\frac{2 \\Delta_{1} \\mu_{2} d k \\left(\\Delta_{2} + 1\\right)}{\\Delta_{1} \\mu_{1} + \\mu_{2}}\\\\0 & \\frac{\\mu_{2} \\left(\\Delta_{2} + 1\\right)}{\\Delta_{1} \\mu_{1} + \\mu_{2}}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Q11 = Q[:2, :2]\n",
    "Q12 = Q[:2, 2:]\n",
    "Q21 = Q[2:, :2]\n",
    "Q22 = Q[2:, 2:]\n",
    "\n",
    "TD = Q22**(-1)\n",
    "RD = Q12*TD\n",
    "RU = -TD*Q21\n",
    "TU = Q11 - Q12*TD*Q21\n",
    "\n",
    "TD = sp.together(TD).simplify()\n",
    "RD = sp.together(RD).simplify()\n",
    "RU = sp.together(RU).simplify()\n",
    "TU = sp.together(TU).simplify()\n",
    "\n",
    "display(Math(latex(RD)))\n",
    "display(Math(latex(RU)))\n",
    "display(Math(latex(TD)))\n",
    "display(Math(latex(TU)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "18b108d5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{1}{2} + \\frac{\\mu_{2}}{2 \\mu_{1}} & \\frac{1}{2} - \\frac{\\mu_{2}}{2 \\mu_{1}}\\\\\\frac{1}{2} - \\frac{\\mu_{2}}{2 \\mu_{1}} & \\frac{1}{2} + \\frac{\\mu_{2}}{2 \\mu_{1}}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[1/2 + \\mu_2/(2*\\mu_1), 1/2 - \\mu_2/(2*\\mu_1)],\n",
       "[1/2 - \\mu_2/(2*\\mu_1), 1/2 + \\mu_2/(2*\\mu_1)]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Ts1 = sp.Matrix([\n",
    "    [1,          1],\n",
    "    [mu1*k,   -mu1*k]\n",
    "])\n",
    "\n",
    "Ts2 = sp.Matrix([\n",
    "    [1,          1],\n",
    "    [mu2*k,   -mu2*k]\n",
    "])\n",
    "\n",
    "Q = sp.expand(Ts1.inv() * Ts2)\n",
    "Q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "ecf6eba6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{\\mu_{1} - \\mu_{2}}{\\mu_{1} + \\mu_{2}}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{- \\mu_{1} + \\mu_{2}}{\\mu_{1} + \\mu_{2}}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{2 \\mu_{1}}{\\mu_{1} + \\mu_{2}}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{2 \\mu_{2}}{\\mu_{1} + \\mu_{2}}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Q11 = Q[:1, :1]\n",
    "Q12 = Q[:1, 1:]\n",
    "Q21 = Q[1:, :1]\n",
    "Q22 = Q[1:, 1:]\n",
    "\n",
    "TD = Q22**(-1)\n",
    "RD = Q12*TD\n",
    "RU = -TD*Q21\n",
    "TU = Q11 - Q12*TD*Q21\n",
    "\n",
    "TD = sp.together(TD).simplify()\n",
    "RD = sp.together(RD).simplify()\n",
    "RU = sp.together(RU).simplify()\n",
    "TU = sp.together(TU).simplify()\n",
    "\n",
    "display(Math(latex(RD)))\n",
    "display(Math(latex(RU)))\n",
    "display(Math(latex(TD)))\n",
    "display(Math(latex(TU)))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "py310",
   "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": 5
}