doc/tutorial_matching_automatic.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Automatic Matching\n",
    "This notebooks explains how the automatic matching is implemented on WEST ICRH antenna."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import skrf as rf\n",
    "import numpy as np\n",
    "import sys\n",
    "from tqdm.notebook import tqdm\n",
    "\n",
    "sys.path.append('..')\n",
    "from west_ic_antenna import WestIcrhAntenna"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "The WEST automatic matching of the capacitors of an ICRH antenna is made with a negative feedback loop. ICRH Operator defines the desired input impedance at the T-junction setpoint $Z_{T,SP}$. This setpoint is compared to the actual impedance at the T-junction $Z_T$ and capacitors are adjusted to minimize this difference in realtime. The input impedance at the T-junction is determined from the input impedance at the bidirective coupler located behind the antennas.\n",
    "\n",
    "<img src=\"automatic_matching_control_loop.png\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The different steps are described below:\n",
    "\n",
    "- RF Measurements: the forward $P_{fwd}$ and reflected powers $P_{ref}$ are measured at the bidirective coupler located behind an antenna. The phase difference between these two $\\Phi_\\Gamma$ is also measured. The impedance at the bidirective coupler $Z_{coupler}$ is deduced from the complex reflection coefficient $\\Gamma$ :\n",
    "\n",
    "$$ \\Gamma = \\sqrt{\\frac{P_{ref}}{P_{fwd}}} \\exp⁡ \\left[ j(\\Phi_\\Gamma + \\Phi_c ) \\right]$$\n",
    "\n",
    "$$ Z_{coupler} = Z_0  \\frac{1 + \\Gamma}{1 - \\Gamma} $$\n",
    "\n",
    "where $Z_0$=30 Ω is the characteristic impedance of the transmission line and $\\Phi_c$ a possible phase correction that the IC operator can tune. \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- $Z_T$ Calculation: \n",
    "\n",
    "Let $A$ the ABCD matrix of the assembly of the piece of transmission line between the bidirective coupler and the antenna elements up to the T-junction. This matrix has been deduced from the CAD model of the antenna ($S_{SWIT}$) and from the measurement of the piece of transmission lines ($S_{ciseau}$) for each antennas (which can differ in length). \n",
    "\n",
    "The impedance at the T-junction is deduced from the impedance $Z_{coupler}$ (here asssuming port 1 is located at the T junction):\n",
    "\n",
    "$$ Z_T = \\frac{A_{11} Z_{coupler} + A_{12}}{A_{21} Z_{coupler} + A_{22} } $$\n",
    "\n",
    "NB: if port 1 is instead the 30 Ohm side, the equation would have been:\n",
    "\n",
    "$$ Z_T = \\frac{A_{12}  - Z_{coupler} A_{22}}{A_{21} Z_{coupler} - A_{11} } $$\n",
    "\n",
    "\n",
    "<img src=\"automatic_matching_transmission_line_definition.png\"/>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Normalized error signal: the relative change of the measured error $\\varepsilon$ is :\n",
    "\n",
    "$$\n",
    "\\mathbf{\\varepsilon}\n",
    "=\n",
    "\\left(\n",
    "\\begin{array}{c}\n",
    "\\Re(\\varepsilon) \\\\\n",
    "\\Im(\\varepsilon)\n",
    "\\end{array}\n",
    "\\right)\n",
    "=\n",
    "\\left(\n",
    "\\begin{array}{c}\n",
    "\\Re\\left( \\frac{Z_{T,sp}}{Z_T} \\right) - 1 \\\\\n",
    "\\Im\\left( \\frac{Z_{T,sp}}{Z_T} \\right) \n",
    "\\end{array}\n",
    "\\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Controller: gives a velocity request $V_r=(V_{r,top}, V_{r,bot})$ to the capacitors actuator:\n",
    "\n",
    "$$\n",
    "\\left(\n",
    "\\begin{array}{c}\n",
    "V_{r,top} \\\\\n",
    "V_{r,bot}\n",
    "\\end{array}\n",
    "\\right)\n",
    "=\n",
    "\\mathbf{T} \n",
    "\\times\n",
    "\\mathbf{\\varepsilon}\n",
    "$$\n",
    "\n",
    "The derivation of the matrix $\\mathbf{T}$ is detailed in the next section"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Controler Development\n",
    "### Introduction\n",
    "The impedance at the T-junction $Z_T$ is a function of the capacitances set $\\mathbf{C}=( C_{top}, C_{bot})^t$:\n",
    "\n",
    "$$\n",
    "Z_T = Z_T(\\mathbf{C})\n",
    "$$\n",
    "\n",
    "Close to the desired capacitances solution $\\mathbf{C}_{SP}=(C_{top,SP}, C_{bot,SP})^t$, the previous equation can be approximated thanks to the Taylor theorem:\n",
    "\n",
    "$$\n",
    "Z_T(\\mathbf{C}) \n",
    "\\approx \n",
    "Z_T(\\mathbf{C}_{SP})\n",
    "+ \n",
    "\\nabla Z_T(\\mathbf{C}_{SP})\n",
    "\\cdot\n",
    "(\\mathbf{C} - \\mathbf{C}_{SP})\n",
    "$$\n",
    "with \n",
    "$$\n",
    "\\nabla Z_T(\\mathbf{C}_{SP})\n",
    "=\n",
    "\\left.\n",
    "\\left(\n",
    "\\begin{array}{c}\n",
    "\\frac{\\partial Z_T}{\\partial C_{top}} \\\\\n",
    "\\frac{\\partial Z_T}{\\partial C_{bot}}\n",
    "\\end{array}\n",
    "\\right)\n",
    "\\right|_{\\mathbf{C}=\\mathbf{C}_{SP}}\n",
    "$$\n",
    "We define later $\\delta\\mathbf{C}=\\mathbf{C} - \\mathbf{C}_{SP}$ and recalls that the first term on the right hand side $Z_T(\\mathbf{C}_{SP})$ corresponds to the desired Set Point $Z_{T,SP}$. \n",
    "\n",
    "Splitting $Z_T$ into real and imaginary parts leads to:\n",
    "\n",
    "$$\n",
    "\\left( \n",
    "\\begin{array}{c} \n",
    "    \\Re(Z_T)(\\mathbf{C}) \\\\ \n",
    "    \\Im(Z_T)(\\mathbf{C})\n",
    "\\end{array}\n",
    "\\right)\n",
    "\\approx\n",
    "\\left( \n",
    "\\begin{array}{c} \n",
    "    \\Re(Z_{T,SP})\\\\ \n",
    "    \\Im(Z_{T,SP})\n",
    "\\end{array}\n",
    "\\right)\n",
    "+ \n",
    "\\mathbf{D}\n",
    "\\cdot\n",
    "\\delta\\mathbf{C}\n",
    "$$\n",
    "where\n",
    "$$\n",
    "\\mathbf{D}\n",
    "=\n",
    "\\left.\n",
    "\\left(\n",
    "\\begin{array}{cc}\n",
    "\\frac{\\partial \\Re(Z_T)}{\\partial C_{top}} & \\frac{\\partial \\Re(Z_T)}{\\partial C_{bot}} \\\\\n",
    "\\frac{\\partial \\Im(Z_T)}{\\partial C_{top}} & \\frac{\\partial \\Im(Z_T)}{\\partial C_{bot}}\n",
    "\\end{array}\n",
    "\\right)\n",
    "\\right|_{\\mathbf{C}=\\mathbf{C}_{opt}}\n",
    "$$\n",
    "\n",
    "The previous formula can be inversed to given the required $\\delta \\mathbf{C}$ increment to reach the Set Point:\n",
    "\n",
    "$$\n",
    "\\delta\\mathbf{C} \n",
    "\\approx \n",
    "\\mathbf{D}^{-1} \n",
    "\\left(\n",
    "\\left( \n",
    "\\begin{array}{c} \n",
    "    \\Re(Z_T)(\\mathbf{C}) \\\\ \n",
    "    \\Im(Z_T)(\\mathbf{C})\n",
    "\\end{array}\n",
    "\\right)\n",
    "-\n",
    "\\left( \n",
    "\\begin{array}{c} \n",
    "    \\Re(Z_{T,SP})\\\\ \n",
    "    \\Im(Z_{T,SP})\n",
    "\\end{array}\n",
    "\\right)\n",
    "\\right)\n",
    "$$\n",
    "\n",
    "The goal of the this section is to determine the elements of the Jacobian matrix $\\mathbf{D}$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Numerical Values of the Jacobian\n",
    "Let's define a WEST ICRH antenna for a given frequency:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f0 = 55  # MHz\n",
    "freq0 = rf.Frequency(f0, f0, unit='MHz', npoints=1)\n",
    "antenna = WestIcrhAntenna(frequency=freq0, \n",
    "                          front_face='../west_ic_antenna/data/Sparameters/front_faces/TOPICA/S_TSproto12_55MHz_Hmode_LAD6-2.5cm.s4p')\n",
    "antenna.front_face_Rc()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we calculate the ideal matching point (for an half side only):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "C_match = antenna.match_one_side(f_match=f0*1e6, \n",
    "                                 side='left', solution_number=1)\n",
    "print(C_match)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "power = [1,1e-6]\n",
    "phase = [0,0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's calculate the derivate of $Z_T$ around the match point:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "delta_C_tops = np.linspace(-10, 10, 51)\n",
    "delta_C_bots = np.linspace(-10, 10, 51)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# derivatives vs Ctop\n",
    "z_Ts = []\n",
    "z_T_SP = (2.9-0.3j)\n",
    "\n",
    "for delta_C_top in delta_C_tops:\n",
    "    # z_T\n",
    "    _z_T = antenna.z_T(power, phase, Cs=[C_match[0]+delta_C_top, C_match[1], 150, 150])\n",
    "    z_Ts.append(_z_T[:,0]) # left side\n",
    "z_Ts = np.array(z_Ts).squeeze()\n",
    "dz_Ts_dCtop = np.diff(z_Ts)\n",
    "\n",
    "# normalized error and derivative\n",
    "eps_Ctop = (z_T_SP - z_Ts)/z_Ts\n",
    "deps_dCtop = np.diff(eps_Ctop)\n",
    "\n",
    "# diff vs Cbot\n",
    "z_Ts = []\n",
    "for delta_C_bot in delta_C_bots:\n",
    "    # z_T\n",
    "    _z_T = antenna.z_T(power, phase, Cs=[C_match[0], C_match[1]+delta_C_bot, 150, 150])\n",
    "    z_Ts.append(_z_T[:,0]) # left side\n",
    "z_Ts = np.array(z_Ts).squeeze()/z_T_SP\n",
    "dz_Ts_dCbot = np.diff(z_Ts)\n",
    "\n",
    "# normalized error and derivative\n",
    "eps_Cbot = (z_T_SP - z_Ts)/z_Ts\n",
    "deps_dCbot = np.diff(eps_Cbot)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below are the traces of the derivatives of $Z_T$ around the Set Point:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(2, 2, sharex=True, sharey=True)\n",
    "axes[0,0].plot(delta_C_tops[:-1], dz_Ts_dCtop.real, label='$\\partial \\Re(z_T)/\\partial C_{top}$')\n",
    "axes[0,1].plot(delta_C_bots[:-1], dz_Ts_dCbot.real, label='$\\partial \\Re(z_T)/\\partial C_{bot}$')\n",
    "axes[1,0].plot(delta_C_tops[:-1], dz_Ts_dCtop.imag, label='$\\partial \\Im(z_T)/\\partial C_{top}$')\n",
    "axes[1,1].plot(delta_C_bots[:-1], dz_Ts_dCbot.imag, label='$\\partial \\Im(z_T)/\\partial C_{bot}$')\n",
    "axes[1,0].set_xlabel('$\\delta C$')\n",
    "axes[1,1].set_xlabel('$\\delta C$')\n",
    "[ax.legend() for ax in axes.ravel()]\n",
    "[ax.axvline(0, color='gray', ls='--', alpha=0.8) for ax in axes.ravel()]\n",
    "fig.subplots_adjust(wspace=0.0, hspace=0.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below are the plot of the derivatives of the error signal $\\varepsilon$ around the Set Point:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(2, 2, sharex=True, sharey=True)\n",
    "axes[0,0].plot(delta_C_tops[:-1], deps_dCtop.real, label=r'$\\partial \\Re(\\varepsilon)/\\partial C_{top}$')\n",
    "axes[0,1].plot(delta_C_bots[:-1], deps_dCbot.real, label=r'$\\partial \\Re(\\varepsilon)/\\partial C_{bot}$')\n",
    "axes[1,0].plot(delta_C_tops[:-1], deps_dCtop.imag, label=r'$\\partial \\Im(\\varepsilon)/\\partial C_{top}$')\n",
    "axes[1,1].plot(delta_C_bots[:-1], deps_dCbot.imag, label=r'$\\partial \\Im(\\varepsilon)/\\partial C_{bot}$')\n",
    "axes[1,0].set_xlabel('$\\delta C$')\n",
    "axes[1,1].set_xlabel('$\\delta C$')\n",
    "[ax.legend() for ax in axes.ravel()]\n",
    "[ax.axvline(0, color='gray', ls='--', alpha=0.8) for ax in axes.ravel()]\n",
    "fig.subplots_adjust(wspace=0.0, hspace=0.0)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below are the values of the derivative at the Set Point:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "idx = np.argmin(np.abs(delta_C_bots))\n",
    "print('dz_T/dCtop:', dz_Ts_dCtop[idx])\n",
    "print('dz_T/dCbot:', dz_Ts_dCbot[idx])\n",
    "\n",
    "print('de/dCtop:', deps_dCtop[idx])\n",
    "print('de/dCbot:', deps_dCbot[idx])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As the $\\mathbf{D}$ matrix depends of the conditions (coupling, excitation), it is not possible to determine it exactly in advance in real experiments. \n",
    "\n",
    "An approximative solution is to take only the sign of the derivatives: hence few iterations are needed until converging to the solution. \n",
    "\n",
    "Note that this algorithm will diverge if the start point is not close enough from the optimal solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "D = np.array([[dz_Ts_dCtop[idx].real, dz_Ts_dCbot[idx].real],\n",
    "              [dz_Ts_dCtop[idx].imag, dz_Ts_dCbot[idx].imag]])\n",
    "print(np.sign(D))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, for solution 1, one would obtain:\n",
    "$$\n",
    "\\left(\n",
    "\\begin{array}{cc}\n",
    "-1 & -1 \\\\\n",
    "-1 & +1\n",
    "\\end{array}\n",
    "\\right)\n",
    "$$\n",
    "and for solution 2:\n",
    "$$\n",
    "\\left(\n",
    "\\begin{array}{cc}\n",
    "-1 & -1 \\\\\n",
    "+1 & -1\n",
    "\\end{array}\n",
    "\\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note the sign change if one uses the Jacobian relative to the error signal $\\mathbf{\\varepsilon}$ instead of the Jacobian relative to the $z_T$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerical Examples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the example below, we setup and antenna and we look forward a solution by iterating on the capacitor predictor detailled below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f0 = 55  # MHz\n",
    "freq0 = rf.Frequency(f0, f0, unit='MHz', npoints=1)\n",
    "antenna = WestIcrhAntenna(frequency=freq0, \n",
    "                          front_face='../west_ic_antenna/data/Sparameters/front_faces/TOPICA/S_TSproto12_55MHz_Hmode_LAD6-2.5cm.s4p')\n",
    "C_match_sol1 = antenna.match_one_side(f_match=f0*1e6, side='left', solution_number=1)\n",
    "C_match_sol2 = antenna.match_one_side(f_match=f0*1e6, side='left', solution_number=2)\n",
    "print(C_match_sol1)\n",
    "print(C_match_sol2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Matching a single side"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to illustrate the convergence, we first calculate a capacitor map $\\mathrm{VSWR}(C_{top}, C_{bot})$ for the left side around the ideal match point:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "delta_Cs = np.linspace(-10, +10, num=31, endpoint=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "C_top_lefts, C_bot_lefts = np.meshgrid(C_match_sol1[0] + delta_Cs, \n",
    "                                       C_match_sol1[1] + delta_Cs)\n",
    "print(C_top_lefts.size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "power = [1, 1e-12]  # small value on right side to avoid division by zero\n",
    "phase = [0, 0]\n",
    "\n",
    "vswrs = []\n",
    "for (C_top, C_bot) in tqdm(np.nditer([C_top_lefts, C_bot_lefts]), total=C_top_lefts.size):\n",
    "    _vswr = antenna.vswr_act(power, phase, Cs=[C_top, C_bot, 120, 120])\n",
    "    vswrs.append(_vswr)\n",
    "vswrs = np.array(vswrs).squeeze()\n",
    "vswrs_left = vswrs[:,0].reshape(C_top_lefts.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Which looks like that:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "cs=ax.contourf(C_bot_lefts, C_top_lefts, vswrs_left, np.linspace(1, 6, 41))\n",
    "cs2=ax.contour(C_bot_lefts, C_top_lefts, vswrs_left, np.linspace(1, 6, 11), colors='k', alpha=0.6)\n",
    "ax.clabel(cs2, inline=1, fontsize=10)\n",
    "ax.set_xlabel('$C_{bot}$ [pF]')\n",
    "ax.set_ylabel('$C_{top}$ [pF]')\n",
    "ax.set_title('SWR at feeding line - left side')\n",
    "ax.grid(True, alpha=0.2)\n",
    "ax.plot(C_match_sol1[1], C_match_sol1[0], 'o', color='C1', label=\"solution 1 ($C_{top} > C_{bot}$)\")\n",
    "ax.plot(C_match_sol2[1], C_match_sol2[0], 'o', color='C2', label=\"solution 2 ($C_{top} > C_{bot}$)\")\n",
    "ax.legend()\n",
    "fig.colorbar(cs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's define a starting point and iterate on the capacitor predictor. Playing with the code below, you'll see rapidly that the convergence is heavily linked to the starting point. Using a starting point far from the expected solution or with capacitor inversed with respected to the desired solution will lead to the divergence of the algorithm:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "C0_start = [60,40,120,120]\n",
    "sol_num = 1\n",
    "\n",
    "C0 = list(C0_start)  # new list to avoid reference\n",
    "C0s = []\n",
    "cont = True\n",
    "iterations = 0 \n",
    "while cont:\n",
    "    C_next_left, C_next_right, eps = antenna.capacitor_predictor(power, phase, Cs=C0, solution_number=sol_num)\n",
    "    C0 = [C_next_left.squeeze()[0], C_next_left.squeeze()[1], 120, 120]\n",
    "\n",
    "    iterations += 1\n",
    "    if iterations > 1 and (np.abs(C0s[-1][0] - C0[0]) < 0.01) and (np.abs(C0s[-1][1] - C0[1]) < 0.01):\n",
    "        cont = False\n",
    "    if iterations > 60:\n",
    "        cont = False\n",
    "    C0s.append([C0[0], C0[1]])\n",
    "print(f'Stopped after {iterations} iterations')    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's illustrate the trace of the matching convergence:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "cs=ax.contourf(C_bot_lefts, C_top_lefts, vswrs_left, np.linspace(1, 6, 41))\n",
    "cs2=ax.contour(C_bot_lefts, C_top_lefts, vswrs_left, np.linspace(1, 6, 11), colors='k', alpha=0.6)\n",
    "ax.clabel(cs2, inline=1, fontsize=10)\n",
    "ax.set_xlabel('$C_{bot}$ [pF]')\n",
    "ax.set_ylabel('$C_{top}$ [pF]')\n",
    "ax.set_title('SWR at feeding line - left side')\n",
    "ax.grid(True, alpha=0.2)\n",
    "ax.plot(C_match_sol1[1], C_match_sol1[0], 'o', color='C1', label=\"solution 1 ($C_{top} > C_{bot}$)\")\n",
    "ax.plot(C_match_sol2[1], C_match_sol2[0], 'o', color='C2', label=\"solution 2 ($C_{top} > C_{bot}$)\")\n",
    "ax.legend()\n",
    "fig.colorbar(cs)\n",
    "\n",
    "\n",
    "C0s = np.array(C0s)\n",
    "ax.plot(C0_start[1], C0_start[0], 'x', color='r')\n",
    "ax.plot(C0s[:,1], C0s[:,0], '-x', color='r', alpha=0.8)\n",
    "print(C0s[-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below is another representation of the convergence:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.plot(np.arange(len(C0s)), C0s, lw=2)\n",
    "ax.set_xlabel('# Steps')\n",
    "ax.set_ylabel('Capacitance [pF]')\n",
    "ax.grid(True, alpha=.4)\n",
    "ax.set_title('Capacitor Matching Alg Convergence')\n",
    "ax.axhline(eval(f'C_match_sol{sol_num}')[0], ls='--', color='gray')\n",
    "ax.axhline(eval(f'C_match_sol{sol_num}')[1], ls='--', color='gray')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Matching both sides at the same time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "power = [1, 1]\n",
    "phase = [0, np.pi]\n",
    "# reference solution we wish to obtain\n",
    "C_opt_dipole = antenna.match_both_sides(f_match=55e6, power=power, phase=phase)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "C0_start = [60, 40, 60, 40] # note that the start point matches solution 1 (Ctop>Cbot)\n",
    "sol_num = 1\n",
    "\n",
    "C0 = list(C0_start)  # new list to avoid reference\n",
    "C0s = []\n",
    "cont = True\n",
    "iterations = 0 \n",
    "while cont:\n",
    "    C_next_left, C_next_right, eps = antenna.capacitor_predictor(power, phase, Cs=C0, solution_number=sol_num)\n",
    "    C0 = [C_next_left.squeeze()[0], C_next_left.squeeze()[1], \n",
    "          C_next_right.squeeze()[0], C_next_right.squeeze()[1]]\n",
    "    print(C0)\n",
    "    iterations += 1\n",
    "    if iterations > 1 and (np.abs(C0s[-1][0] - C0[0]) < 0.01) and (np.abs(C0s[-1][1] - C0[1]) < 0.01):\n",
    "        cont = False\n",
    "    if iterations > 60:\n",
    "        cont = False\n",
    "    C0s.append([C0[0], C0[1], C0[2], C0[3]])\n",
    "print(f'Stopped after {iterations} iterations')    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.plot(np.arange(len(C0s)), C0s, lw=2)\n",
    "ax.set_xlabel('# Steps')\n",
    "ax.set_ylabel('Capacitance [pF]')\n",
    "ax.grid(True, alpha=.4)\n",
    "ax.set_title('Capacitor Matching Alg Convergence')\n",
    "[ax.axhline(C, ls='--', color=f'C{idx}') for idx, C in enumerate(C_opt_dipole)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The antenna is matched:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "antenna.vswr_act(power, phase, Cs=C0s[-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The automatic matching algorithm is implemented in the  `matching_both_sides_iterative` method of the antenna object:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "C_opt_2 = antenna.match_both_sides_iterative(f_match=55e6, power=power, phase=phase, Cs=[50, 50, 50, 50])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "antenna.vswr_act(power, phase, Cs=C_opt_2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Correspondance Between Error Signal Amplitude and Capacitance\n",
    "In the control room, the IC operator job is often to improve the antenna tuning for the next shot. Because turning on the matching feedback is not always desirable, another option is to look to the error signals $\\mathbf{\\varepsilon}$ generated in the latest pulse and to deduce from them the correction to apply to the capacitors.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, let's build a half-matched antenna:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f0 = 55  # MHz\n",
    "freq0 = rf.Frequency(f0, f0, unit='MHz', npoints=1)\n",
    "antenna = WestIcrhAntenna(frequency=freq0, \n",
    "                          front_face='../west_ic_antenna/data/Sparameters/front_faces/TOPICA/S_TSproto12_55MHz_Hmode_LAD6-2.5cm.s4p')\n",
    "C_match = antenna.match_one_side(f_match=f0*1e6, \n",
    "                                 side='left', solution_number=1)\n",
    "print(C_match)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's depart one capacitor from the ideal point and look to the evolution of the error signal:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "delta_Cs = np.linspace(-10, +10, 11)\n",
    "\n",
    "_Cs = []\n",
    "for delta_C in tqdm(delta_Cs):\n",
    "    Cs = [C_match[0]+delta_C, *C_match[1:]]\n",
    "    _C = antenna.capacitor_predictor(power, phase, Cs)[0].squeeze()\n",
    "    _Cs.append(_C)\n",
    "_Cs = np.array(_Cs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.plot(delta_Cs, _Cs-C_match[:2])\n",
    "ax.plot([-10, 10], [-10, 10])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Velocity Controler\n",
    "The velocity controler matrix $\\mathbf{T}$ is in fact $\\mathbf{D}^{-1}$, since the velocity $\\mathbf{V}$ is proportional to the change $\\delta C$:\n",
    "\n",
    "$$\n",
    "\\delta C = \\mathbf{D}^{-1} (z_T - z_{T,SP})\n",
    "$$\n",
    "so that\n",
    "$$\n",
    "\\mathbf{V} = k \\times \\delta C \n",
    "= k \\mathbf{D}^{-1} (z_T - z_{T,SP})\n",
    "= \\mathbf{T} (z_T - z_{T,SP})\n",
    "$$\n",
    "Following the previous result on $\\mathbf{D}$, the matrix $\\mathbf{T}$ can be parametrized as:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "$$\n",
    "\\mathbf{T}\n",
    "=\n",
    "k\\times\n",
    "\\left(\n",
    "\\begin{array}{cc}\n",
    "-S \\cos \\alpha & S K_{ri} \\sin\\alpha \\\\\n",
    "-\\sin\\alpha & - K_{ri} \\cos\\alpha\n",
    "\\end{array}\n",
    "\\right)\n",
    "$$\n",
    "\n",
    "Where $k$ is gain, $S$ and $\\alpha$ two parameters used to define the choice of the matching solution (\"1\" for $C_{top}>C_{bot}$, \"2\" for the opposite). $K_{ri}$ is used to define the relative weight between real and imaginary part of the error signal $\\mathbf{\\varepsilon}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.core.display import HTML\n",
    "def _set_css_style(css_file_path):\n",
    "    \"\"\"\n",
    "    Read the custom CSS file and load it into Jupyter\n",
    "    Pass the file path to the CSS file\n",
    "    \"\"\"\n",
    "    styles = open(css_file_path, \"r\").read()\n",
    "    s = '<style>%s</style>' % styles\n",
    "    return HTML(s)\n",
    "\n",
    "_set_css_style('custom.css')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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": 4
}