{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6ca6ddd0",
   "metadata": {},
   "source": [
    "# 自由边界和刚性边界的反射系数\n",
    "\n",
    "以 P-SV 为例\n",
    "\n",
    "+ Author: Zhu Dengda  \n",
    "+ Email:  zhudengda@mail.iggcas.ac.cn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "b5a7a686",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy as sp"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5811a951",
   "metadata": {},
   "source": [
    "## 动态解"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "ec9378bd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 定义基本变量\n",
    "k, b, d, a, mu, Omg, kb, ka, w, Delta = sp.symbols(r'k b d, a mu \\Omega k_{b} k_{a} \\omega \\Delta')\n",
    "lam, rho = sp.symbols(r'\\lambda \\rho')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "789ac906",
   "metadata": {},
   "outputs": [],
   "source": [
    "D = sp.Matrix([\n",
    "    [k,            b,             k,            -b],\n",
    "    [a,            k,             -a,           k],\n",
    "    [2*mu*Omg,     2*k*mu*b,      2*mu*Omg,     -2*k*mu*b],\n",
    "    [2*k*mu*a,     2*mu*Omg,      -2*k*mu*a,    2*mu*Omg]\n",
    "])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "239e2dd2",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_R(D, at_bottom:bool):\n",
    "    if not at_bottom:\n",
    "        return - (D[:, 2:])**(-1) * (D[:, :2])\n",
    "    else:\n",
    "        return - (D[:, :2])**(-1) * (D[:, 2:])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8feffbc4",
   "metadata": {},
   "source": [
    "### 自由界面（z平面牵引力为0）"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "eafaff07",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{\\Omega^{2} + a b k^{2}}{- \\Omega^{2} + a b k^{2}} & \\frac{2 \\Omega b k}{- \\Omega^{2} + a b k^{2}}\\\\\\frac{2 \\Omega a k}{- \\Omega^{2} + a b k^{2}} & \\frac{\\Omega^{2} + a b k^{2}}{- \\Omega^{2} + a b k^{2}}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[(\\Omega**2 + a*b*k**2)/(-\\Omega**2 + a*b*k**2),           2*\\Omega*b*k/(-\\Omega**2 + a*b*k**2)],\n",
       "[          2*\\Omega*a*k/(-\\Omega**2 + a*b*k**2), (\\Omega**2 + a*b*k**2)/(-\\Omega**2 + a*b*k**2)]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# z+ 侧， 例如顶层顶界面\n",
    "R = get_R(D[2:, :], False)\n",
    "R.simplify()\n",
    "R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "002464b0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{\\Omega^{2} + a b k^{2}}{- \\Omega^{2} + a b k^{2}} & \\frac{2 \\Omega b k}{\\Omega^{2} - a b k^{2}}\\\\\\frac{2 \\Omega a k}{\\Omega^{2} - a b k^{2}} & \\frac{\\Omega^{2} + a b k^{2}}{- \\Omega^{2} + a b k^{2}}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[(\\Omega**2 + a*b*k**2)/(-\\Omega**2 + a*b*k**2),            2*\\Omega*b*k/(\\Omega**2 - a*b*k**2)],\n",
       "[           2*\\Omega*a*k/(\\Omega**2 - a*b*k**2), (\\Omega**2 + a*b*k**2)/(-\\Omega**2 + a*b*k**2)]])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# z- 侧， 例如底层底界面\n",
    "R = get_R(D[2:, :], True)\n",
    "R.simplify()\n",
    "R"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b171d3f9",
   "metadata": {},
   "source": [
    "### 刚性界面（位移为0）"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "1f167bde",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{a b + k^{2}}{a b - k^{2}} & \\frac{2 b k}{a b - k^{2}}\\\\\\frac{2 a k}{a b - k^{2}} & \\frac{a b + k^{2}}{a b - k^{2}}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[(a*b + k**2)/(a*b - k**2),        2*b*k/(a*b - k**2)],\n",
       "[       2*a*k/(a*b - k**2), (a*b + k**2)/(a*b - k**2)]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# z+ 侧\n",
    "R = get_R(D[:2, :], False)\n",
    "R.simplify()\n",
    "R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "8dd46e32",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{a b + k^{2}}{a b - k^{2}} & - \\frac{2 b k}{a b - k^{2}}\\\\- \\frac{2 a k}{a b - k^{2}} & \\frac{a b + k^{2}}{a b - k^{2}}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[(a*b + k**2)/(a*b - k**2),       -2*b*k/(a*b - k**2)],\n",
       "[      -2*a*k/(a*b - k**2), (a*b + k**2)/(a*b - k**2)]])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# z- 侧\n",
    "R = get_R(D[:2, :], True)\n",
    "R.simplify()\n",
    "R"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "76a678aa",
   "metadata": {},
   "source": [
    "## 静态解"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "2c6d0581",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 在静态解中，矩阵 D 和相对深度 z-zj 有关, 因此 z+ 侧和 z- 侧需使用两种 D 矩阵\n",
    "D1 = sp.Matrix([\n",
    "    [1,         -(1 + 2*k*Delta*d),            1,            -(1 - 2*k*Delta*d)],\n",
    "    [1,          (1 - 2*k*Delta*d),           -1,            -(1 + 2*k*Delta*d)],\n",
    "    [2*mu*k,     2*mu*Delta*k*(1 - 2*k*d),   2*mu*k,      2*mu*Delta*k*(1 + 2*k*d)],\n",
    "    [2*mu*k,    -2*mu*Delta*k*(1 + 2*k*d),  -2*mu*k,      2*mu*Delta*k*(1 - 2*k*d)],\n",
    "])\n",
    "\n",
    "z = 0\n",
    "D2 = sp.Matrix([\n",
    "    [1,         -(1 + 2*k*Delta*z),            1,            -(1 - 2*k*Delta*z)],\n",
    "    [1,          (1 - 2*k*Delta*z),           -1,            -(1 + 2*k*Delta*z)],\n",
    "    [2*mu*k,     2*mu*Delta*k*(1 - 2*k*z),   2*mu*k,      2*mu*Delta*k*(1 + 2*k*z)],\n",
    "    [2*mu*k,    -2*mu*Delta*k*(1 + 2*k*z),  -2*mu*k,      2*mu*Delta*k*(1 - 2*k*z)],\n",
    "])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15d392c1",
   "metadata": {},
   "source": [
    "### 自由界面（z平面牵引力为0）"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "8afb9768",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0 & - \\Delta\\\\- \\frac{1}{\\Delta} & 0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[        0, -\\Delta],\n",
       "[-1/\\Delta,       0]])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# z+ 侧\n",
    "R = get_R(D2[2:, :], False)\n",
    "R.simplify()\n",
    "R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "0faa5c1f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}- 2 d k & \\Delta \\left(- 4 d^{2} k^{2} - 1\\right)\\\\- \\frac{1}{\\Delta} & - 2 d k\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[   -2*d*k, \\Delta*(-4*d**2*k**2 - 1)],\n",
       "[-1/\\Delta,                    -2*d*k]])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# z- 侧\n",
    "R = get_R(D1[2:, :], True)\n",
    "R.simplify()\n",
    "R"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "415c011b",
   "metadata": {},
   "source": [
    "### 刚性界面（位移为0）"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "849e8723",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0 & 1\\\\1 & 0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[0, 1],\n",
       "[1, 0]])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# z+ 侧\n",
    "R = get_R(D2[:2, :], False)\n",
    "R.simplify()\n",
    "R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "ef82bd64",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}2 \\Delta d k & 4 \\Delta^{2} d^{2} k^{2} + 1\\\\1 & 2 \\Delta d k\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[2*\\Delta*d*k, 4*\\Delta**2*d**2*k**2 + 1],\n",
       "[           1,              2*\\Delta*d*k]])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# z- 侧\n",
    "R = get_R(D1[:2, :], True)\n",
    "R.simplify()\n",
    "R"
   ]
  }
 ],
 "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
}