cigreTB880_case02_trefoil_touching HDPE ducts.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "647901f1-6d85-4372-92dc-ff8d2354df61",
   "metadata": {
    "tags": []
   },
   "source": [
    "# CIGRE TB 880 - 2022\n",
    "## Power cable rating examples for calculation tool verification\n",
    "## Case #0-2: Sub-case study with touching HDPE ducts\n",
    "\n",
    "Reference IEC standards:\n",
    "- \"IEC 60287-1-1, Electric Cables - Calculation of the Current Rating - Part 1-1: Current Rating Equations (100 % load factor) and Calculation of Losses – General\", IEC,  2014-11.  \n",
    "- \"IEC 60287-2-1, Electric Cables – Calculation of the Current Rating – Part 2-1: Thermal Resistance – Calculation of Thermal Resistance, Edition 2.0\", IEC, 2015-2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 717,
   "id": "f27d4e7f-0172-4478-a595-63ec4a5f2497",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import libraries\n",
    "import math as mt\n",
    "#import numpy as np\n",
    "#import pandas as pd\n",
    "import matplotlib\n",
    "#matplotlib.use('TkAgg')\n",
    "import matplotlib.pyplot as plt "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3604b992-3112-4825-b45a-617695c13105",
   "metadata": {},
   "source": [
    "# 4.7 Case #0-2: Sub-case study with touching HDPE ducts\n",
    "\n",
    "<div>\n",
    "<img src=\"./img/4_7_trefoil_HDPE_pipes.png\" width=\"150\"/> \n",
    "</div>\n",
    "\n",
    "**Input data:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 718,
   "id": "82ce3f24-f0e2-4ac5-b9b8-057ad0175386",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 4.2.1 Grid parameters\n",
    "\n",
    "U = 132 #[kv] Voltage level\n",
    "f = 50 #[Hz] Frequency"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e1f051e-fd4b-490a-ab4e-24c94a5ffcbe",
   "metadata": {},
   "source": [
    "<div>\n",
    "<img src=\"./img/img2_xlpe132kV.png\" width=\"300\"/> <img src=\"./img/img3_cross.png\" width=\"450\"/>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 719,
   "id": "735eb1b1-ac83-46fc-b39c-da838c8f339a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[30.3, 33.3, 64.3, 66.89999999999999, 68.49999999999999, 75.49999999999999]\n"
     ]
    }
   ],
   "source": [
    "# 4.2.2 Cable parameters\n",
    "\n",
    "# List of the cable composition: [Conductor, Conductor_screen, XLPE_insulation, Insulation_screen, Sheath, Oversheath]\n",
    "thk = [0, 1.5, 15.5, 1.3, 0.8, 3.5] #[mm] Thickness of the layer\n",
    "dia = [0]*6 #[mm] Diameter over the layer. Initialize the size of list \n",
    "dia[0] = 30.3 #[mm] Conductor diameter (630 mm2)\n",
    "\n",
    "\n",
    "for i in range(1,len(dia)):\n",
    "    dia[i] = dia[i-1] + 2*thk[i]\n",
    "print(dia)\n",
    "\n",
    "# Conductor\n",
    "R_0 = 28.3e-6 #[Ohm/m] DC resistance at 20°C. From IEC 60228 Table 2\n",
    "alpha_20 = 3.93e-3 #[K^-1] Temperature coefficient. From IEC 60287-1-1 Table 1 (Copper)\n",
    "k_s = 1 # Skin effect coefficient (round stranded conductor and extruded insulation system). From IEC 60287-1-1 Table 2\n",
    "k_p = 1 # Proximity effect coefficient (round stranded conductor and extruded insulation system). From IEC 60287-1-1 Table 2\n",
    "\n",
    "# Insulation layer\n",
    "epsilon = 2.5 # Relative permittivity. From IEC 60287-1-1 Table 3\n",
    "tan_delta = 0.001 # Loss factor. From IEC 60287-1-1 Table 3\n",
    "rho_i = 3.5 #[Km/W] Thermal resistivity. From CIGRE TB 640 Table 1\n",
    "rho_sc = 2.5 #[Km/W] Thermal resistivity screen (semi-conducting layer). CIGRE WG B1.56, IEC 60949\n",
    "\n",
    "# Metal sheath\n",
    "t_s = thk[4] #[mm] From Manufacturer datasheet\n",
    "rho_s_20 = 2.84e-8 #[Ohm/m] Resistivity at 20°c. From IEC 60287-1-1 Table 1\n",
    "alpha_s_20 = 4.03e-3 #[K^-1] Temperature coefficient. From IEC 60287-1-1 Table 1 (Aluminium)\n",
    "\n",
    "# Oversheath\n",
    "rho_Tos = 3.5 #[Km/W] Thermal resistivity. From IEC 60287-2-1 Table 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 720,
   "id": "04f470c7-a635-4761-8a4e-2f27c63076a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 4.2.5 Operating conditions\n",
    "\n",
    "theta_c_max = 90 #[°C] Maximum operating temperature of conductor. User\n",
    "theta_a = 20 #[°C] Ambient ground temperature. From IEC 60287-3-1 Section 3.2.1\n",
    "rho_Tsoil = 1 #[Km/W] Thermal resistivity. From IEC 60287-3-1 Section 3.2.2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c97cea0-4b59-48f4-8c7e-d0072f004986",
   "metadata": {},
   "source": [
    "##  4.7 Case #0-2: Sub-case study with touching HDPE ducts\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "616e12f1-ee33-4de9-a0f5-016be0c5a1c6",
   "metadata": {},
   "source": [
    "### 4.4.1 Capacitance\n",
    "\n",
    "$ C = \\frac{\\epsilon}{18\\ln(\\frac{D_i}{d_c})} 10^{-9} $ (F/m). From IEC 60287-1-1 Section 2.2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 721,
   "id": "993418d5-0a09-43e5-a8d5-1b2211d7999e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "C = 2.1107662202E-10 F/m\n"
     ]
    }
   ],
   "source": [
    "# 4.4.1 Capacitance\n",
    "\n",
    "d_c = dia[1] #[mm] diameter of conductor, including screen, if any\n",
    "D_i = dia[2] #[mm] external diameter of the insulation (excluding screen)\n",
    "\n",
    "C = (epsilon/(18*mt.log(D_i/d_c)))*1e-9 #[F/m] From IEC 60287-1-1 Section 2.2\n",
    "print(\"C = %.10E F/m\" % C)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 722,
   "id": "7fab0e9f-c7a2-4273-8192-1e2b1914dd64",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "X = 8.9202596895E-05 Ω/m\n"
     ]
    }
   ],
   "source": [
    "# 4.7.1.1 Reactance\n",
    "\n",
    "omega = 2*mt.pi*f #[rad/s] Angular frequency\n",
    "D0 = 140 #[mm] -  diameter of duct (external)\n",
    "Dd = 119.4 #[mm] -  diameter of duct (internal)\n",
    "Dpipe=(D0-Dd)/2 #[mm] - thickness of pipe\n",
    "\n",
    "s = D0#dia[5] #[mm] distance between conductor axes in the electrical section being considered\n",
    "d = D_i + 2*thk[3] + thk[4] #[mm] mean diameter of the sheath\n",
    "\n",
    "\n",
    "X = 2*omega*1e-7*mt.log(2*s/d) #[Ohm/m] From IEC 60287-1-1 Section 2.3.1\n",
    "print(\"X = %.10E \\u03A9/m\" % X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 723,
   "id": "8bd37659-8cb9-4329-a1f1-d5b8d57e87be",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Rs0 = 1.6691286499E-04 Ω/m\n"
     ]
    }
   ],
   "source": [
    "# 4.4.3 Electrical resistance of the metal sheath at 20°C\n",
    "\n",
    "A_s = mt.pi*d*t_s #[mm2] cross section\n",
    "\n",
    "R_s0 = rho_s_20/(A_s*1e-6) #[Ohm/m]\n",
    "print(\"Rs0 = %.10E \\u03A9/m\" % R_s0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 724,
   "id": "8e009965-6a4a-422c-baf0-9f2a64ed77e1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "T1 = 0.4198714890 Km/W\n"
     ]
    }
   ],
   "source": [
    "# 4.4.4 Thermal resistance between conductor and sheath\n",
    "\n",
    "t_1 = thk[1:4] #[mm] insulation thickness\n",
    "d_c_arr = dia[0:3] #[mm] diameter\n",
    "rho_T = [rho_sc, rho_i, rho_sc] #[Km/W] Thermal resistivity\n",
    "T1_layer = [0]*len(rho_T) #[Km/W] Thermal resistance\n",
    "\n",
    "for i in range(len(T1_layer)):\n",
    "    T1_layer[i] = (rho_T[i]/(2*mt.pi))*mt.log(1+2*t_1[i]/d_c_arr[i]) #[Km/W] From IEC 60287-2-1 Section 4.1.2.1\n",
    "\n",
    "T1 = sum(T1_layer) #[Km/W]\n",
    "print(\"T1 = %.10f Km/W\" % T1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 725,
   "id": "1fa6cb36-8d13-4b68-bb91-8cf052a9de24",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "T3 = 0.0541996092 Km/W\n"
     ]
    }
   ],
   "source": [
    "# 4.7.1.2 Thermal resistance of the oversheath\n",
    "\n",
    "t_3 = thk[5] #[mm] oversheath thickness\n",
    "D_a = dia[4] #[mm] reference diameter\n",
    "\n",
    "T3 = (1/(2*mt.pi))*rho_Tos*mt.log(1+2*t_3/D_a) #[Km/W] From IEC 60287-2-1 Section 4.1.4.1\n",
    "T3 = T3 #[Km/W] From IEC 60287-2-1 Section 4.1.4.1; without multiplication by a factor 1.6\n",
    "print(\"T3 = %.10f Km/W\" % T3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 726,
   "id": "2bd97694-089a-49fe-b90a-925ea00a3eb7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "T4_1 = 0.3520961015 Km/W\n",
      "T4_2 = 0.08866064719 Km/W\n",
      "T4_3 = 1.38002093961 Km/W\n",
      "T4 = 1.82077768833 Km/W\n"
     ]
    }
   ],
   "source": [
    "#4.7.2 External thermal resistance\n",
    "#The external thermal resistance of a cable in a duct consists of three parts:\n",
    "# a) The thermal resistance of the air space between the cable surface and duct internal surface T′4,\n",
    "# b) The thermal resistance of the duct itself, T’′4,\n",
    "# c) The external thermal resistance of the duct T’’′4.\n",
    "\n",
    "\n",
    "#4.7.2.1 Thermal resistance between cable and duct (or pipe) T4’\n",
    "Uconstant=1.87 #constant for plastic duct; table 4\n",
    "V=0.312 #constant for plastic duct; table 4\n",
    "Y=0.0037 #constant for plastic duct; table 4\n",
    "De = dia[5] #[mm], diameter\n",
    "Om=70 #[celsius], maximum operating temperature\n",
    "\n",
    "T4_1 = Uconstant/(1+0.1*(V+Y*Om)*De)\n",
    "print(\"T4_1 = %.10f Km/W\" % T4_1)\n",
    "\n",
    "\n",
    "#4.7.2.2 Thermal resistance of the duct (or pipe) itself T4′′\n",
    "T4_2 = (1/(2*mt.pi))*rho_Tos*mt.log(D0/Dd)\n",
    "print(\"T4_2 = %.11f Km/W\" % T4_2)\n",
    "\n",
    "#4.7.2.3 External thermal resistance of the duct (or pipe) T4′′′\n",
    "L = 1000 #[mm]\n",
    "u = 2*L/D0\n",
    "rho_T_ex = 1 #[Km/W]\n",
    "T4_3 = 1/(2*mt.pi)*(mt.log(2*u)+2*mt.log(u))\n",
    "print(\"T4_3 = %.11f Km/W\" % T4_3)\n",
    "\n",
    "\n",
    "T4 = T4_1+T4_2+T4_3 #[Km/W]\n",
    "print(\"T4 = %.11f Km/W\" % T4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 727,
   "id": "2dd3e5ff-c2f1-45ae-952c-0c3d6958b721",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wd = 0.3851382172 W/m\n"
     ]
    }
   ],
   "source": [
    "# 4.4.7 Dielectric losses\n",
    "\n",
    "U_0 = U/mt.sqrt(3)*1e3 #[V]\n",
    "\n",
    "W_d = omega*C*U_0*U_0*tan_delta #[W/m] From IEC 60287-1-1 Section 2.2\n",
    "\n",
    "print(\"Wd = %.10f W/m\" % W_d)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6346190d-c4c9-48f4-a3a5-a9d4d7f0e054",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 4.5 case #0-1: Results dependent on the temperature\n",
    "### 4.5.1  Results linked to the maximum operating temperature\n",
    "### 4.5.2 First iteration & 4.5.3 Next iterations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 728,
   "id": "871b93a1-733c-4503-abc3-37f9555caa4c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration:  1\n",
      "λ1 = 0.8387035149\n",
      "I = 680.6595349744 A\n",
      "Wc = 17.8923931138 W/m\n",
      "Ws = 15.0064129937 W/m\n",
      "θj = 80.6026632061 °C\n",
      "θs = 82.4066399825 °C\n",
      "θc = 90.0000000000 °C\n",
      "θs = 82.4066399825 °C\n",
      "Iteration:  2\n",
      "λ1 = 0.8342165243\n",
      "I = 682.8094655127 A\n",
      "Wc = 18.0056014123 W/m\n",
      "Ws = 15.0205702283 W/m\n",
      "θj = 80.5482271071 °C\n",
      "θs = 82.3591070457 °C\n",
      "θc = 90.0000000000 °C\n",
      "θs = 82.3591070457 °C\n",
      "Iteration:  3\n",
      "λ1 = 0.8343048317\n",
      "I = 682.8140883734 A\n",
      "Wc = 18.0058452217 W/m\n",
      "Ws = 15.0223636681 W/m\n",
      "θj = 80.5480143204 °C\n",
      "θs = 82.3590046771 °C\n",
      "θc = 90.0000000000 °C\n",
      "θs = 82.3590046771 °C\n",
      "Iteration:  4\n",
      "λ1 = 0.8343050219\n",
      "I = 682.8144604072 A\n",
      "Wc = 18.0058648428 W/m\n",
      "Ws = 15.0223834628 W/m\n",
      "θj = 80.5480039457 °C\n",
      "θs = 82.3589964387 °C\n",
      "θc = 90.0000000000 °C\n",
      "θs = 82.3589964387 °C\n",
      "Iteration:  5\n",
      "λ1 = 0.8343050372\n",
      "I = 682.8144647689 A\n",
      "Wc = 18.0058650729 W/m\n",
      "Ws = 15.0223839303 W/m\n",
      "θj = 80.5480038113 °C\n",
      "θs = 82.3589963421 °C\n",
      "θc = 90.0000000000 °C\n",
      "θs = 82.3589963421 °C\n",
      "Iteration:  6\n",
      "λ1 = 0.8343050374\n",
      "I = 682.8144648684 A\n",
      "Wc = 18.0058650781 W/m\n",
      "Ws = 15.0223839379 W/m\n",
      "θj = 80.5480038084 °C\n",
      "θs = 82.3589963399 °C\n",
      "θc = 90.0000000000 °C\n",
      "θs = 82.3589963399 °C\n",
      "Iteration:  7\n",
      "λ1 = 0.8343050374\n",
      "I = 682.8144648700 A\n",
      "Wc = 18.0058650782 W/m\n",
      "Ws = 15.0223839380 W/m\n",
      "θj = 80.5480038084 °C\n",
      "θs = 82.3589963399 °C\n",
      "θc = 90.0000000000 °C\n",
      "θs = 82.3589963399 °C\n",
      "Iteration:  8\n",
      "λ1 = 0.8343050374\n",
      "I = 682.8144648700 A\n",
      "Wc = 18.0058650782 W/m\n",
      "Ws = 15.0223839380 W/m\n",
      "θj = 80.5480038084 °C\n",
      "θs = 82.3589963399 °C\n",
      "θc = 90.0000000000 °C\n",
      "θs = 82.3589963399 °C\n",
      "Iteration:  9\n",
      "λ1 = 0.8343050374\n",
      "I = 682.8144648700 A\n",
      "Wc = 18.0058650782 W/m\n",
      "Ws = 15.0223839380 W/m\n",
      "θj = 80.5480038084 °C\n",
      "θs = 82.3589963399 °C\n",
      "θc = 90.0000000000 °C\n",
      "Convergence reached after  9 iterations\n"
     ]
    }
   ],
   "source": [
    "# Criteria of convergence\n",
    "tempTol = 1e-11 #[K] Temperature tolerance threshold for convergence\n",
    "currTol = 1e-11 #[A] Current tolerance threshold for convergence\n",
    "maxIte = 10 # Maximum number of iterations\n",
    "\n",
    "# Other parameters\n",
    "n = 1 # Single-core cable\n",
    "d_c = dia[0] # Conductor diameter\n",
    "lambda_2 = 0 # The cable is unarmoured\n",
    "T2 = 0 # The cable is unarmoured\n",
    "# NOTE: IEC ArmoUred Cable is not the same that NEC(NFPA 70) Armored Cable (Type AC - Article 320), which is similar to Metal-Clad Cable (Type MC Article 330)\n",
    "\n",
    "# Variable Preallocation\n",
    "R = [0]*maxIte #[Ohm/m] AC resistance of conductor\n",
    "R_s = [0]*maxIte #[Ohm/m] Resistance of sheath\n",
    "lambda_1 = [0]*maxIte # Loss factor\n",
    "I = [0]*maxIte #[A] Current rating\n",
    "W_c = [0]*maxIte #[W/m] Conductor losses\n",
    "W_s = [0]*maxIte #[W/m] Sheath losses\n",
    "theta_c = [0]*maxIte #[°C] Conductor operation temperature\n",
    "theta_s = [0]*maxIte #[°C] Sheath temperature\n",
    "\n",
    "# Initial values\n",
    "theta_c[0] = theta_c_max # Maximum allowable conductor temperature\n",
    "theta_s[0] = theta_c_max-10 # Assumption to initializate!\n",
    "\n",
    "for i in range(maxIte):\n",
    "    print(\"Iteration: %2d\" % (i+1))\n",
    "\n",
    "\n",
    "    T4_1 = Uconstant/(1+0.1*(V+Y*Om)*De)\n",
    "    T4 = T4_1+T4_2+T4_3 #[Km/W]\n",
    "    #print(\"T4 iteration = %.10f Km/W\" % T4)\n",
    "\n",
    "    # 4.5.1  Results linked to the maximum operating temperature\n",
    "    # 4.5.1.1 DC resistance of conductor\n",
    "    R_prime = R_0*(1+alpha_20*(theta_c[i]-20)) #[Ohm/m] From IEC 60287-1-1 Section 2.1.1\n",
    "    #print(\"R' = %.10E W/m\" % R_prime)\n",
    "    \n",
    "    # 4.5.1.2 Skin effect factor\n",
    "    x_s = mt.sqrt(8*mt.pi*f/R_prime*1e-7*k_s) # From IEC 60287-1-1 Section 2.1.2\n",
    "    # For 0 < x_s <= 2,8\n",
    "    y_s = x_s**4/(192+0.8*x_s**4)\n",
    "    #print(\"xs = %.10f\\nys = %.10f\" % (x_s, y_s))\n",
    "    \n",
    "    # 4.5.1.3 Proximity effect factor \n",
    "    x_p = mt.sqrt(8*mt.pi*f/R_prime*1e-7*k_p) # From IEC 60287-1-1 Section 2.1.4.1\n",
    "    y_p = x_p**4/(192+0.8*x_p**4)*(d_c/s)**2*(0.312*(d_c/s)**2+1.18/(x_p**4/(192+0.8*x_p**4)+0.27))\n",
    "    #print(\"xp = %.10f\\nyp = %.10f\" % (x_p, y_p))\n",
    "    \n",
    "    # 4.5.1.4 AC resistance of conductor \n",
    "    R[i] = R_prime*(1+y_s+y_p) #[Ohm/m] From IEC 60287-1-1 Section 2.1\n",
    "    #print(\"R = %.10E Ohm/m\" % R[i])\n",
    "    \n",
    "    # 4.5.2 First iteration & 4.5.3 Next iterations\n",
    "    # 4.5.2.1 Resistance of sheath at operating temperature\n",
    "    R_s[i] = R_s0*(1+alpha_s_20*(theta_s[i]-20)) #[Ohm/m] From IEC 60287-1-1 Section 2.3\n",
    "    #print(\"Rs = %.10E Ohm/m\" % R_s[i])\n",
    "    \n",
    "    # 4.5.2.2 Losses caused by circulating currents\n",
    "    lambda_1_prime = (R_s[i]/R[i])*(1/(1+(R_s[i]/X)**2)) # From IEC 60287-1-1 Section 2.3.1\n",
    "    lambda_1_prime2 = 0\n",
    "    lambda_1[i] = lambda_1_prime + lambda_1_prime2\n",
    "    print(\"\\u03BB1 = %.10f\" % lambda_1[i])\n",
    "    \n",
    "    # 4.5.2.3 Permissible current rating\n",
    "    Delta_theta = theta_c[i] - theta_a\n",
    "    I[i] = mt.sqrt((Delta_theta-W_d*(0.5*T1+n*(T2+T3+T4)))/(R[i]*T1+n*R[i]*(1+lambda_1[i])*T2+n*R[i]*(1+lambda_1[i]+lambda_2)*(T3+T4))) #[A] From IEC 60287-1-1 Section 1.4.1.1\n",
    "    print(\"I = %.10f A\" % I[i])\n",
    "    \n",
    "    # 4.5.2.4 Losses in conductor and sheath\n",
    "    W_c[i] = R[i]*I[i]**2 #[W/m]\n",
    "    W_s[i] = lambda_1[i]*W_c[i] #[W/m]\n",
    "    print(\"Wc = %.10f W/m\\nWs = %.10f W/m\" % (W_c[i], W_s[i]))\n",
    "    \n",
    "    # 4.5.2.5 Temperature on cable components\n",
    "    theta_j = theta_a + n*(W_c[i]+W_s[i]+W_d)*T4 #[°C] Jacket. From IEC 60287-1-1 Section 1.4.1.1\n",
    "    theta_s[i] = theta_j + n*(W_c[i]+W_s[i]+W_d)*T3 #[°C] Metal sheath \n",
    "    theta_c[i] = theta_s[i] + n*(W_c[i]+W_d/2)*T1 #[°C] Conductor\n",
    "    print(\"\\u03B8j = %.10f °C\\n\\u03B8s = %.10f °C\\n\\u03B8c = %.10f °C\" % (theta_j, theta_s[i], theta_c[i]))\n",
    "    \n",
    "    # Convergence evaluation\n",
    "    Delta_theta_s = abs(theta_s[i] - theta_s[i-1])\n",
    "    Delta_theta_c = abs(theta_c[i] - theta_c[i-1])\n",
    "    Delta_theta1 = max(Delta_theta_s, Delta_theta_c)\n",
    "    Delta_I = abs(I[i] - I[i-1])\n",
    "    \n",
    "    if (Delta_theta1 < tempTol) and (Delta_I < currTol):\n",
    "        print(\"Convergence reached after %2d iterations\" % (i+1))\n",
    "        break\n",
    "    elif i+1 == maxIte:\n",
    "        print(\"Convergence criteria not reached after %2d iterations\" % maxIte)\n",
    "    else:        \n",
    "        Om = theta_j-0.5*T4_1*(W_c[i]+W_d+W_s[i]) #mean air gap temperature value change\n",
    "        theta_s[i+1] = theta_s[i] # Update for next iteration \n",
    "        print(\"\\u03B8s = %.10f °C\" % (theta_s[i+1]))\n",
    "        theta_c[i+1] = theta_c[i]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "637311d6-9f48-47e4-8413-62b6267f9d5e",
   "metadata": {},
   "source": [
    "### Results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 729,
   "id": "5957dc03-af4f-4841-870a-9018df3a14ad",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkkAAAGwCAYAAAC99fF4AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAABAe0lEQVR4nO3deXhU5f338c9kmTCQhS2BhISwhSUQFglFwIIKRZC6Xg+iDTVIrUVZRJQqKOLGYkWLRQuClkWr+KsYjP4eREWWiiIgoiA8JEAUJYQAgUwWyDbn+QNnMCRCJpnJLHm/ritXyZkz53wHvZxP7/O979tkGIYhAAAAVBLg6QIAAAC8ESEJAACgGoQkAACAahCSAAAAqkFIAgAAqAYhCQAAoBqEJAAAgGoEeboAX2Wz2ZSdna2wsDCZTCZPlwMAAGrAMAwVFBQoJiZGAQGXHisiJNVSdna24uLiPF0GAACohR9//FGxsbGXPIeQVEthYWGSzv8lh4eHe7gaAABQE1arVXFxcY7v8UshJNWS/RFbeHg4IQkAAB9Tk1YZGrcBAACqQUgCAACoBiEJAACgGoQkAACAahCSAAAAqkFIAgAAqAYhCQAAoBqEJAAAgGoQkgAAAKpBSAIAAKgGIQkAAKAahCQAAIBqsMEtUEeGYajCZqjcZsj2859ttovOkXHRey5+veo1L/36pa9/mV8v+35n6wMAdwgLCVZE42CP3Z+QhHp3IKdA/y/HqgqbceHHMGSzXRw2pAqb7fz//vz6L4PIxe8tt/18DeMXf/75/PJf/Lnyey/cw/7eSq//sjbDUEXFRecYRpVAAQBwjfuu7qi/jujqsfsTklCvzhSX6saXPlNJue3yJ6NGTKaLfq/yuukSr1383l+/2MXvBQB3Cwrw7H95CEmoV/uyrSoptyk0JEh92jZVUIBJgQEmBZjO/6/jx2RSQIBJQQHn/zfwotcDTBe/JgUGBCgwQI5r/fJ1+7Wq3Kua61atQwq66NqXqzPAZLpkeDFd9OJlw8vFBwAAbkdIQr3KOF4gSbqyQwu9mprs4WoAAPh1zG5DvcrILZQkdW4V6uFKAAC4NEIS6lXmzyNJnVuFebgSAAAujZCEemMYhjKOnx9JSmAkCQDg5QhJqDcnCkqUf7ZMASapYyQhCQDg3QhJqDf2UaT4Fk3UKDjQw9UAAHBphCTUG/vMtoQoRpEAAN6PkIR6k5lL0zYAwHcQklBvaNoGAPgSQhLqxfmZbedHkrq0ZiQJAOD9CEmoF8etJSo4V67AAJPat2zi6XIAALgsQhLqhX0UqV2LxgoJYmYbAMD7EZJQLzJYaRsA4GMISagXmY6mbUISAMA3EJJQLzIc0/+Z2QYA8A2EJLidYRiOkSQetwEAfAUhCW6XnX9OhSXlCgowqV0LZrYBAHwDIQluZ2/abt+yicxB/CsHAPANfGPB7TKZ2QYA8EGEJLgd25EAAHwRIQlux0gSAMAXEZLgVjabocxc+8w2RpIAAL6DkAS3OnrmrIpLKxQcaFI8M9sAAD6EkAS3yvx5EckOLUMVHMi/bgAA38G3FtyKpm0AgK8iJMGt2NgWAOCrCElwqwvbkTCSBADwLYQkuI3NZuhgrv1xGyNJAADfQkiC2/x0+qzOllXIHBig+OaNPV0OAABOISTBbez9SB0imyiImW0AAB/DNxfcJiOXpm0AgO8iJMFtaNoGAPgyQhLcxv64jaZtAIAvIiTBLSp+MbOtCyEJAOCDCElwix/zilVSblNIUIDimNkGAPBBhCS4hf1RW6eoUAUGmDxcDQAAziMkwS0yc+1N2zxqAwD4JkIS3OJC0zYz2wAAvomQBLfIsE//j2IkCQDgmwhJcLnyCpsO8bgNAODjCElwuR/yilVaYZMlOFCxzSyeLgcAgFohJMHlMn8xsy2AmW0AAB9FSILL2fuRaNoGAPgyj4eko0ePauzYsWrRooUsFouSkpK0c+dOx+uFhYWaNGmSYmNjZbFYlJiYqCVLljhez8vL0+TJk9WlSxdZLBa1bdtWU6ZMUX5+/iXvO27cOJlMpko/I0aMcNvnbEjsM9voRwIA+LIgT9789OnTGjRokK655hqtW7dOkZGRyszMVLNmzRznTJs2TZ9++qneeOMNtWvXTh999JHuu+8+xcTE6MYbb1R2drays7O1YMECJSYm6ocfftCECROUnZ2td95555L3HzFihJYvX+74PSQkxG2ftSFhY1sAgD/waEh69tlnFRcXVymotG/fvtI5n3/+uVJTU3X11VdLku655x698sor2r59u2688Ub16NFDa9ascZzfsWNHzZkzR2PHjlV5ebmCgn79I4aEhKh169au/VANXFmFTYdP/vy4jen/AAAf5tHHbenp6UpOTtbo0aMVFRWlPn36aNmyZZXOGThwoNLT03X06FEZhqGNGzcqIyNDw4cP/9Xr5ufnKzw8/JIBSZI2bdqkqKgodenSRffee69OnTr1q+eWlJTIarVW+kFVP5wqUlmFocbmQLVpysw2AIDv8mhIOnz4sBYvXqyEhAStX79e9957r6ZMmaKVK1c6zlm0aJESExMVGxsrs9msESNG6OWXX9bgwYOrvebJkyf19NNP65577rnkvUeMGKFVq1Zpw4YNevbZZ7V582aNHDlSFRUV1Z4/b948RUREOH7i4uJq/8H9mKNpm5ltAAAfZzIMw/DUzc1ms5KTk/X55587jk2ZMkU7duzQF198IUlasGCBli1bpgULFig+Pl5btmzRjBkzlJaWpmHDhlW6ntVq1e9+9zs1b95c6enpCg4OrnEthw8fVseOHfXJJ59o6NChVV4vKSlRSUlJpXvFxcU5Rq1w3sJPMrTwk0z9n76xWjC6l6fLAQCgEqvVqoiIiBp9f3u0Jyk6OlqJiYmVjnXr1s3RY3T27FnNnDlTaWlpGjVqlCSpZ8+e2r17txYsWFApJBUUFGjEiBEKCwtTWlqaUwFJkjp06KCWLVvq4MGD1YakkJAQGrtrgKZtAIC/8OjjtkGDBunAgQOVjmVkZCg+Pl6SVFZWprKyMgUEVC4zMDBQNpvN8bvVatXw4cNlNpuVnp6uRo0aOV3LTz/9pFOnTik6OroWnwR2Fza2pWkbAODbPBqSHnjgAW3btk1z587VwYMH9eabb2rp0qWaOHGiJCk8PFxDhgzR9OnTtWnTJmVlZWnFihVatWqVbrnlFkkXAlJRUZFee+01Wa1W5eTkKCcnp1J/UdeuXZWWlibp/NpL06dP17Zt2/T9999rw4YNuummm9SpUyddd9119f8X4SdKy23KOlkkiTWSAAC+z6OP2/r166e0tDTNmDFDTz31lNq3b6+FCxcqJSXFcc7q1as1Y8YMpaSkKC8vT/Hx8ZozZ44mTJggSdq1a5e+/PJLSVKnTp0qXT8rK0vt2rWTJB04cMCxwGRgYKC+/fZbrVy5UmfOnFFMTIyGDx+up59+mkdqdfD9qSKV2wyFhgQpJsL50TwAALyJRxu3fZkzjV8NxQffZmvSm1+rd1xTrZ04yNPlAABQhTPf3x7flgT+I4OmbQCAHyEkwWUy2bMNAOBHCElwGTa2BQD4E0ISXKKkvELfnyqWREgCAPgHQhJcIutkkSpshsIaBalVODMEAQC+j5AEl7jQtB0mk4k92wAAvo+QBJe40LTNzDYAgH8gJMElHNuRRNGPBADwD4QkuMQvH7cBAOAPCEmos3NlFfrhlH3PNh63AQD8AyEJdXboRKFshhRhCVZkGDPbAAD+gZCEOsv8xXYkzGwDAPgLQhLqzNG0TT8SAMCPEJJQZ46m7Sj6kQAA/oOQhDrLzGXPNgCA/yEkoU7OllboSN75Pdt43AYA8CeEJNTJoROFMgypWeNgtQw1e7ocAABchpCEOvll0zYz2wAA/oSQhDrJ+MX0fwAA/AkhCXVyYWNb+pEAAP6FkIQ6ychlY1sAgH8iJKHWikvL9WPeWUk8bgMA+B9CEmrtYO75fqQWTcxqEcqebQAA/0JIQq3Zm7YTGEUCAPghQhJqjaZtAIA/IySh1tjYFgDgzwhJqDX747YuhCQAgB8iJKFWikrKdfQMM9sAAP6LkIRayfx5ZltkWIiaNmbPNgCA/yEkoVYyHE3bjCIBAPwTIQm1Yp/ZxkrbAAB/RUhCrVzY2JaQBADwT4Qk1AqP2wAA/o6QBKdZz5XpWP45SayRBADwX4QkOC3z50dtrcJDFGEJ9nA1AAC4ByEJTmM7EgBAQxBUk5OuuOIKpy5qMpmUnp6uNm3a1KooeDfHxrbMbAMA+LEahaTdu3frwQcfVGjo5Zt0DcPQ/PnzVVJSUufi4J0yc2naBgD4vxqFJEmaPn26oqKianTu888/X+uC4P3Y2BYA0BDUKCRlZWUpMjKyxhfdt2+fYmJial0UvFf+2TIdt54fJUxgJAkA4Mdq1LgdHx8vk8l02fP27t0rSYqLi1NgYGDdKoNXsjdtR0c0UngjZrYBAPxXnWe3FRQUaOnSpfrNb36jXr16uaImeDFH0zaP2gAAfq7WIWnLli1KTU1VdHS0FixYoGuvvVbbtm1zZW3wQo6VtqN41AYA8G81btyWpJycHK1YsUKvvfaarFarbrvtNpWUlGjt2rVKTEx0V43wIhdmtjGSBADwbzUeSbrhhhvUpUsXffvtt1q4cKGys7O1aNEid9YGL3ThcRsjSQAA/1bjkaR169ZpypQpuvfee5WQkODOmuClzhSX6kSBfWYbI0kAAP9W45Gkzz77TAUFBerbt6/69++vl156SSdPnnRnbfAy9lGkNk0tCg1x6kktAAA+p8Yh6corr9SyZct07Ngx/eUvf9Hq1asVExMjm82mjz/+WAUFBe6sE17gwiKSPGoDAPg/p2e3NWnSROPHj9dnn32mPXv26MEHH9T8+fMVFRWlG2+80R01wkuwsS0AoCGp0zpJXbp00d/+9jf99NNPeuutt1xVE7yU/XEbIQkA0BDUeTFJSQoMDNTNN9+s9PR0V1wOXoqNbQEADUmNQtKtt94qq9Va44umpKQoNze31kXB++QVlepkYakkqRMLSQIAGoAaTVF67733dOLEiRpd0DAMvf/++3r66acVFRVVp+LgPexN23HNLWpsZmYbAMD/1ejbzjAMde7c2d21wIs5mraj6EcCADQMNQpJGzdudPrCbdq0cfo98F5sbAsAaGhqFJKGDBni7jrg5Q4cp2kbANCwuGR2G/ybYRiskQQAaHAISbisk4WlOl1cJpNJ6hjJSBIAoGEgJOGy7KNIbZs3lsUc6OFqAACoH4QkXJZjzzZmtgEAGhCnQ9K1116rM2fOVDlutVp17bXXuqImeJmMXPt2JDxqAwA0HE6HpE2bNqm0tLTK8XPnzum///2vS4qCd6FpGwDQENV46eRvv/3W8ed9+/YpJyfH8XtFRYU+/PBD1kbyQ4Zh/GKNJEaSAAANR41DUu/evWUymWQymap9rGaxWLRo0SKXFgfPO1FQovyzZQpgZhsAoIGpcUjKysqSYRjq0KGDtm/frsjISMdrZrNZUVFRCgxk5pO/sY8ixbdookbB/PMFADQcNQ5J8fHxkiSbzea2YuB9LsxsYxQJANCw1Go798zMTG3cuFG5ublVQtPjjz/uksLgHTJzadoGADRMTs9uW7Zsmbp166bHH39c77zzjtLS0hw/a9eudbqAo0ePauzYsWrRooUsFouSkpK0c+dOx+uFhYWaNGmSYmNjZbFYlJiYqCVLljhez8vL0+TJk9WlSxdZLBa1bdtWU6ZMUX5+/iXvaxiGHn/8cUVHR8tisWjYsGHKzMx0un5/R9M2AKChcnok6ZlnntGcOXP08MMP1/nmp0+f1qBBg3TNNddo3bp1ioyMVGZmppo1a+Y4Z9q0afr000/1xhtvqF27dvroo4903333KSYmRjfeeKOys7OVnZ2tBQsWKDExUT/88IMmTJig7OxsvfPOO79677/97W/6xz/+oZUrV6p9+/aaNWuWrrvuOu3bt0+NGjWq82fzB+dntjGSBABomEyGYRjOvCE8PFy7d+9Whw4d6nzzRx55RFu3br3k+ko9evTQmDFjNGvWLMexvn37auTIkXrmmWeqfc9//vMfjR07VkVFRQoKqpoDDcNQTEyMHnzwQT300EOSpPz8fLVq1UorVqzQ7bffftnarVarIiIilJ+fr/Dw8Mue74ty8s/pynkbFBhg0r6nrlNIEI3bAADf5sz3t9OP20aPHq2PPvqo1sX9Unp6upKTkzV69GhFRUWpT58+WrZsWaVzBg4cqPT0dB09elSGYWjjxo3KyMjQ8OHDf/W69g9eXUCSzs/Uy8nJ0bBhwxzHIiIi1L9/f33xxRfVvqekpERWq7XSj7+zjyLFt2hMQAIANDhOP27r1KmTZs2apW3btikpKUnBwcGVXp8yZUqNr3X48GEtXrxY06ZN08yZM7Vjxw5NmTJFZrNZqampkqRFixbpnnvuUWxsrIKCghQQEKBly5Zp8ODB1V7z5MmTevrpp3XPPff86n3tC2G2atWq0vFWrVpVWiTzl+bNm6cnn3yyxp/NHzgetbFnGwCgAXI6JC1dulShoaHavHmzNm/eXOk1k8nkVEiy2WxKTk7W3LlzJUl9+vTR3r17tWTJkkohadu2bUpPT1d8fLy2bNmiiRMnKiYmptJIkHR+CG3UqFFKTEzUE0884exHu6QZM2Zo2rRple4VFxfn0nt4m8zj7NkGAGi4nA5JWVlZLrt5dHS0EhMTKx3r1q2b1qxZI0k6e/asZs6cqbS0NI0aNUqS1LNnT+3evVsLFiyoFJIKCgo0YsQIhYWFKS0trcoI1y+1bt1aknT8+HFFR0c7jh8/fly9e/eu9j0hISEKCQmp1ef0VRn26f+tGUkCADQ8Tvck2ZWWlurAgQMqLy+v9c0HDRqkAwcOVDqWkZHhWLiyrKxMZWVlCgioXGZgYGCl9ZmsVquGDx8us9ms9PT0y85Oa9++vVq3bq0NGzZUusaXX36pAQMG1Prz+BPDMHTQMZJESAIANDxOh6Ti4mL96U9/UuPGjdW9e3cdOXJEkjR58mTNnz/fqWs98MAD2rZtm+bOnauDBw/qzTff1NKlSzVx4kRJ52fSDRkyRNOnT9emTZuUlZWlFStWaNWqVbrlllskXQhIRUVFeu2112S1WpWTk6OcnBxVVFQ47tW1a1elpaVJOv9YcOrUqXrmmWeUnp6uPXv26M4771RMTIxuvvlmZ/9K/NKx/HMqKClXUIBJ7Vo08XQ5AADUP8NJU6ZMMfr27Wv897//NZo0aWIcOnTIMAzDWLt2rdG7d29nL2e8//77Ro8ePYyQkBCja9euxtKlSyu9fuzYMWPcuHFGTEyM0ahRI6NLly7G888/b9hsNsMwDGPjxo2GpGp/srKyHNeRZCxfvtzxu81mM2bNmmW0atXKCAkJMYYOHWocOHCgxnXn5+cbkoz8/HynP7Mv2Pj/jhvxD39gDHt+k6dLAQDAZZz5/nZ6naT4+Hi9/fbbuvLKKxUWFqZvvvlGHTp00MGDB3XFFVc0iKnxkv+vk7Rsy2HN+b/7NSopWi+nXOHpcgAAcAm3rpN04sQJRUVFVTleVFQkk8nk7OXgpRwb2zKzDQDQQDkdkpKTk/W///u/jt/twejVV1+l6dmPsB0JAKChc3oJgLlz52rkyJHat2+fysvL9eKLL2rfvn36/PPPq6ybBN9ksxnKzGWNJABAw+b0SNJVV12lb775RuXl5UpKStJHH32kqKgoffHFF+rbt687akQ9O3rmrIpLKxQcaFI8M9sAAA2UUyNJZWVl+stf/qJZs2ZV2WMN/iPz50UkO7QMVXBgrZfSAgDApzn1DRgcHOxYDRv+K+PnRSRp2gYANGRODxPcfPPNWrt2rRtKgbegaRsAgFo0bickJOipp57S1q1b1bdvXzVpUrlnxZkNbuGd2NgWAADJ6cUk27dv/+sXM5l0+PDhOhflC/x1MUmbzVD32et1tqxCGx4coo6RBCUAgP9w5vvbqZEkwzC0adMmRUVFyWKx1KlIeKefTp/V2bIKmQMDFN+8safLAQDAY5zqSTIMQwkJCfrpp5/cVQ88zN6P1CGyiYKY2QYAaMCc+hYMCAhQQkKCTp065a564GEZuTRtAwAg1WJ22/z58zV9+nTt3bvXHfXAw2jaBgDgPKdnt915550qLi5Wr169ZDabq/Qm5eXluaw41L8LG9sykgQAaNicDkkLFy50QxnwBhU2Qwcde7YRkgAADZvTISk1NdUddcAL/JhXrJJym0KCAtSWmW0AgAbO6ZB05MiRS77etm3bWhcDz7I/ausYGarAAJOHqwEAwLOcDknt2rWTyfTrX6AVFRV1Kgiek5lL0zYAAHZOh6Svv/660u9lZWX6+uuv9cILL2jOnDkuKwz1j6ZtAAAucDok9erVq8qx5ORkxcTE6LnnntOtt97qksJQ/zJ+nv7fhZAEAIDz6yT9mi5dumjHjh2uuhzqWYXN0KETzGwDAMDO6ZEkq9Va6XfDMHTs2DE98cQTSkhIcFlhqF8/nCpSablNluBAxTZjXz4AAJwOSU2bNq3SuG0YhuLi4rR69WqXFYb6ZX/U1ikqVAHMbAMAwPmQ9Omnn1YKSQEBAYqMjFSnTp0UFOT05eAlMh1N28xsAwBAqkVIuvrqq91QBjwtg5W2AQCoxOnG7Xnz5ulf//pXleP/+te/9Oyzz7qkKNS/jJzzI0mskQQAwHlOh6RXXnlFXbt2rXK8e/fuWrJkiUuKQv0qq7Dp8MnzI0kJUYwkAQAg1SIk5eTkKDo6usrxyMhIHTt2zCVFoX79cKpIZRWGGpsD1aYpM9sAAJBqEZLi4uK0devWKse3bt2qmJgYlxSF+mWf2ZbAzDYAABycbtz+85//rKlTp6qsrEzXXnutJGnDhg3661//qgcffNDlBcL92I4EAICqnA5J06dP16lTp3TfffeptLRUktSoUSM9/PDDmjFjhssLhPtlHmdjWwAALuZ0SDKZTHr22Wc1a9Ys7d+/XxaLRQkJCQoJCXFHfagHjCQBAFBVrVd/DA0NVb9+/VxZCzygtNymrJNFklgjCQCAX3LZBrfwTd+fKlK5zVBoSJBiIhp5uhwAALwGIamBsz9q6xQVWmVPPgAAGjJCUgOXQdM2AADVIiQ1cPaNbelHAgCgslqFpNdff12DBg1STEyMfvjhB0nSwoUL9d5777m0OLgfM9sAAKie0yFp8eLFmjZtmq6//nqdOXNGFRUVkqSmTZtq4cKFrq4PblRSXqHvTxVL4nEbAAAXczokLVq0SMuWLdOjjz6qwMBAx/Hk5GTt2bPHpcXBvbJOFqnCZigsJEitw5nZBgDALzkdkrKystSnT58qx0NCQlRUVOSSolA/HHu2tWJmGwAAF3M6JLVv3167d++ucvzDDz9Ut27dXFET6glN2wAA/DqnV9yeNm2aJk6cqHPnzskwDG3fvl1vvfWW5s2bp1dffdUdNcJNMghJAAD8KqdD0t133y2LxaLHHntMxcXF+sMf/qCYmBi9+OKLuv32291RI9zkwsa2hCQAAC5Wq73bUlJSlJKSouLiYhUWFioqKsrVdcHNzpVV6PtT9j3bmNkGAMDFar3BrSQ1btxYjRs3dlUtqEeHTxTJZkgRlmBFhoV4uhwAALxOjUJSnz59ajz7adeuXXUqCPUjM9fej8TMNgAAqlOjkHTzzTe7uQzUN1baBgDg0moUkmbPnu3uOlDPDuT83LQdRT8SAADVqXVP0s6dO7V//35JUmJiovr27euyouB+Fx63MZIEAEB1nA5JP/30k+644w5t3bpVTZs2lSSdOXNGAwcO1OrVqxUbG+vqGuFiZ0srdCTv/J5tPG4DAKB6Tq+4fffdd6usrEz79+9XXl6e8vLytH//ftlsNt19993uqBEuduhEoQxDatY4WC1DzZ4uBwAAr+T0SNLmzZv1+eefq0uXLo5jXbp00aJFi/Tb3/7WpcXBPX7ZtM3MNgAAquf0SFJcXJzKysqqHK+oqFBMTIxLioJ7ZThW2qZpGwCAX+N0SHruuec0efJk7dy503Fs586duv/++7VgwQKXFgf3YGNbAAAur0aP25o1a1bpsUxRUZH69++voKDzby8vL1dQUJDGjx/Pmko+IOPnmW0JUYQkAAB+TY1C0sKFC91cBupLcWm5fsw7K4nHbQAAXEqNQlJqaqq760A9OZh7vh+pRROzWoSyZxsAAL+mRiHJarUqPDzc8edLsZ8H72Rv2k5gFAkAgEuqcU/SsWPHFBUVpaZNm1Y7bdwwDJlMJlVUVLi8SLgOTdsAANRMjULSp59+qubNm0uSNm7c6NaC4F5sbAsAQM3UKCQNGTKk2j/D9zjWSGJjWwAALsnpdZI+/PBDffbZZ47fX375ZfXu3Vt/+MMfdPr0aZcWB9cqKinX0TP2mW2MJAEAcClOh6Tp06c7mrf37NmjadOm6frrr1dWVpamTZvm8gLhOpk/z2xrGRqiZk3Ysw0AgEtxeu+2rKwsJSYmSpLWrFmjG264QXPnztWuXbt0/fXXu7xAuE6Go2mbR20AAFyO0yNJZrNZxcXFkqRPPvlEw4cPlyQ1b978sssDwLOY2QYAQM05PZJ01VVXadq0aRo0aJC2b9+ut99+W5KUkZGh2NhYlxcI17mwsS0hCQCAy3F6JOmll15SUFCQ3nnnHS1evFht2rSRJK1bt04jRoxwuoCjR49q7NixatGihSwWi5KSkiptnltYWKhJkyYpNjZWFotFiYmJWrJkSaVrLF26VFdffbXCw8NlMpl05syZy973iSeekMlkqvTTtWtXp+v3JZk8bgMAoMacHklq27atPvjggyrH//73vzt989OnT2vQoEG65pprtG7dOkVGRiozM1PNmjVznDNt2jR9+umneuONN9SuXTt99NFHuu+++xQTE6Mbb7xRklRcXKwRI0ZoxIgRmjFjRo3v3717d33yySeO3+0b9vqjgnNlys4/J4k1kgAAqAmnU8GuXbsUHByspKQkSdJ7772n5cuXKzExUU888YTM5prPmnr22WcVFxen5cuXO461b9++0jmff/65UlNTdfXVV0uS7rnnHr3yyivavn27IyRNnTpVkrRp0yanPktQUJBat27t1Ht8lX1mW6vwEEVYgj1cDQAA3s/px21/+ctflJGRIUk6fPiwbr/9djVu3Fj/+c9/9Ne//tWpa6Wnpys5OVmjR49WVFSU+vTpo2XLllU6Z+DAgUpPT9fRo0dlGIY2btyojIwMR8N4XWRmZiomJkYdOnRQSkqKjhw58qvnlpSUyGq1VvrxJTRtAwDgHKdDUkZGhnr37i1J+s9//qPBgwfrzTff1IoVK7RmzRqnrnX48GEtXrxYCQkJWr9+ve69915NmTJFK1eudJyzaNEiJSYmKjY2VmazWSNGjNDLL7+swYMHO1t6Jf3799eKFSv04YcfavHixcrKytJvf/tbFRQUVHv+vHnzFBER4fiJi4ur0/3rm2Nj2yhCEgAANeH04zbDMGSz2SSdXwLg97//vSQpLi5OJ0+edOpaNptNycnJmjt3riSpT58+2rt3r5YsWaLU1FRJ50PStm3blJ6ervj4eG3ZskUTJ05UTEyMhg0b5mz5DiNHjnT8uWfPnurfv7/i4+P1P//zP/rTn/5U5fwZM2ZUWizTarX6VFBijSQAAJzjdEhKTk7WM888o2HDhmnz5s1avHixpPOLTLZq1cqpa0VHRzsWprTr1q2bY0Tq7NmzmjlzptLS0jRq1ChJ5wPN7t27tWDBgjqFpIs1bdpUnTt31sGDB6t9PSQkRCEhIS67X31jY1sAAJzj9OO2hQsXateuXZo0aZIeffRRderUSZL0zjvvaODAgU5da9CgQTpw4EClYxkZGYqPj5cklZWVqaysTAEBlcsMDAx0jGa5SmFhoQ4dOqTo6GiXXtcb5J8t03FriSQpgZEkAABqxOmRpJ49e2rPnj1Vjj/33HMKDAx06loPPPCABg4cqLlz5+q2227T9u3btXTpUi1dulSSFB4eriFDhmj69OmyWCyKj4/X5s2btWrVKr3wwguO6+Tk5CgnJ8cxCrRnzx6FhYWpbdu2at68uSRp6NChuuWWWzRp0iRJ0kMPPaQbbrhB8fHxys7O1uzZsxUYGKg77rjD2b8Sr2dv2o6OaKTwRsxsAwCgJly2MFCjRo2cfk+/fv2UlpamGTNm6KmnnlL79u21cOFCpaSkOM5ZvXq1ZsyYoZSUFOXl5Sk+Pl5z5szRhAkTHOcsWbJETz75pON3e1P38uXLNW7cOEnSoUOHKvVM/fTTT7rjjjt06tQpRUZG6qqrrtK2bdsUGRnp9Ofwdo6mbR61AQBQYybDMIzLndS8eXNlZGSoZcuWatasmUwm06+em5eX59ICvZXValVERITy8/MVHh7u6XIu6Yn077Ti8+9191Xt9djvEy//BgAA/JQz3981Gkn6+9//rrCw86MQCxcurHOBqF+ZuayRBACAs2oUkuzT8S/+M3zDhcdtNG0DAFBTte5Jys3NVW5ubpVZZj179qxzUXCdM8WlOlFgn9nGSBIAADXldEj66quvlJqaqv379+vidiaTyaSKigqXFYe6s48itWlqUWiI/27gCwCAqzn9rTl+/Hh17txZr732mlq1anXJJm543oVFJHnUBgCAM5wOSYcPH9aaNWsci0jCu7GxLQAAteP0ittDhw7VN998445a4AYXNrZlJAkAAGc4PZL06quvKjU1VXv37lWPHj0UHFx5Becbb7zRZcWh7pj+DwBA7Tgdkr744gtt3bpV69atq/IajdveJa+oVCcLSyVJnRhJAgDAKU4/bps8ebLGjh2rY8eOyWazVfohIHkXe9N2bDOLmjCzDQAApzgdkk6dOqUHHnhArVq1ckc9cCGatgEAqD2nQ9Ktt96qjRs3uqMWuJi9aZuQBACA85x+BtO5c2fNmDFDn332mZKSkqo0bk+ZMsVlxaFuMhwjSfQjAQDgLJNx8bLZl9G+fftfv5jJpMOHD9e5KF/gzC7CnnLF0x8rr6hUH0y+Sj3aRHi6HAAAPM6Z72+nRpIMw9CmTZsUFRUli8VSpyLhXicLS5RXVCqTSeoYyUgSAADOcqonyTAMJSQk6KeffnJXPXAR+6O2ts0by2IO9HA1AAD4HqdCUkBAgBISEnTq1Cl31QMXyXSstE3TNgAAteH07Lb58+dr+vTp2rt3rzvqgYscoGkbAIA6cXp225133qni4mL16tVLZrO5Sm9SXl6ey4pD7bFGEgAAdeN0SFq4cKEbyoArGYZxYWNbRpIAAKgVp0NSamqqO+qAC50oKFH+2TIFMLMNAIBac7onSZIOHTqkxx57THfccYdyc3MlSevWrdN3333n0uJQO/ZRpPgWTdQomJltAADUhtMhafPmzUpKStKXX36pd999V4WF57+Qv/nmG82ePdvlBcJ59un/CVGMIgEAUFtOh6RHHnlEzzzzjD7++GOZzWbH8WuvvVbbtm1zaXGoncxcmrYBAKgrp0PSnj17dMstt1Q5HhUVpZMnT7qkKNQNTdsAANSd0yGpadOmOnbsWJXjX3/9tdq0aeOSolB752e2MZIEAEBdOR2Sbr/9dj388MPKycmRyWSSzWbT1q1b9dBDD+nOO+90R41wwnFriQrOlSswwKQOkU08XQ4AAD7L6ZA0d+5cde3aVXFxcSosLFRiYqIGDx6sgQMH6rHHHnNHjXCCfRQpvkVjhQQxsw0AgNpyep0ks9msZcuW6fHHH9eePXtUWFioPn36KCEhwR31wUmOR23s2QYAQJ3UOCTZbDY999xzSk9PV2lpqYYOHarZs2dX2ZYEnmXf2JY92wAAqJsaP26bM2eOZs6cqdDQULVp00YvvviiJk6c6M7aUAsZP0//T6BpGwCAOqlxSFq1apX++c9/av369Vq7dq3ef/99/fvf/5bNZnNnfXCCYRg66BhJIiQBAFAXNQ5JR44c0fXXX+/4fdiwYTKZTMrOznZLYXDesfxzKigpV1CASe1bMrMNAIC6qHFIKi8vV6NGjSodCw4OVllZmcuLQu3Ym7bbtWwic1CttuUDAAA/q3HjtmEYGjdunEJCQhzHzp07pwkTJqhJkwujFu+++65rK0SN0bQNAIDr1DgkpaamVjk2duxYlxaDumGlbQAAXKfGIWn58uXurAMukJFL0zYAAK5C44qfOD+zzT6SxOM2AADqipDkJ46eOaui0goFB5oU34KZbQAA1BUhyU/Ym7Y7tAxVcCD/WAEAqCu+Tf2EvWk7gUdtAAC4BCHJTxxgZhsAAC5FSPITrJEEAIBrEZL8gM1m6ODP0//Z2BYAANcgJPmBn06f1dmyCpkDAxTfvLGnywEAwC8QkvyAvWm7Q2QTBTGzDQAAl+Ab1Q9k5NK0DQCAqxGS/ABN2wAAuB4hyQ9cWCOJkSQAAFyFkOTjKn4xs43HbQAAuA4hycf9mFesknKbQoIC1JaZbQAAuAwhycfZH7V1jAxVYIDJw9UAAOA/CEk+LjOXpm0AANyBkOTjaNoGAMA9CEk+LuM4TdsAALgDIcmHVdgMHTrB4zYAANyBkOTDfjhVpNJymxoFByiuGTPbAABwJUKSD7M/ausUFaoAZrYBAOBShCQflvlz03bnKPqRAABwNUKSD8uwT/9vTUgCAMDVCEk+zDGSRNM2AAAuR0jyUeUVNh0+USRJSuBxGwAALkdI8lHfnypWaYVNjc2BatPU4ulyAADwO4QkH2V/1JbAzDYAANyCkOSj7NP/2Y4EAAD3ICT5qAyatgEAcCtCko9iY1sAANzL4yHp6NGjGjt2rFq0aCGLxaKkpCTt3LnT8XphYaEmTZqk2NhYWSwWJSYmasmSJZWusXTpUl199dUKDw+XyWTSmTNnanTvl19+We3atVOjRo3Uv39/bd++3ZUfzW1Ky23KOnl+Zhsb2wIA4B4eDUmnT5/WoEGDFBwcrHXr1mnfvn16/vnn1axZM8c506ZN04cffqg33nhD+/fv19SpUzVp0iSlp6c7zikuLtaIESM0c+bMGt/77bff1rRp0zR79mzt2rVLvXr10nXXXafc3FyXfkZ3+P5UkcpthkJDghQT0cjT5QAA4JdMhmEYnrr5I488oq1bt+q///3vr57To0cPjRkzRrNmzXIc69u3r0aOHKlnnnmm0rmbNm3SNddco9OnT6tp06aXvHf//v3Vr18/vfTSS5Ikm82muLg4TZ48WY888shla7darYqIiFB+fr7Cw8Mve74rffBttia9+bV6xzXV2omD6vXeAAD4Mme+vz06kpSenq7k5GSNHj1aUVFR6tOnj5YtW1bpnIEDByo9PV1Hjx6VYRjauHGjMjIyNHz48Frft7S0VF999ZWGDRvmOBYQEKBhw4bpiy++qPY9JSUlslqtlX48xT6zjaZtAADcx6Mh6fDhw1q8eLESEhK0fv163XvvvZoyZYpWrlzpOGfRokVKTExUbGyszGazRowYoZdfflmDBw+u9X1PnjypiooKtWrVqtLxVq1aKScnp9r3zJs3TxEREY6fuLi4Wt+/ri5sR0I/EgAA7hLkyZvbbDYlJydr7ty5kqQ+ffpo7969WrJkiVJTUyWdD0nbtm1Tenq64uPjtWXLFk2cOFExMTGVRoLcbcaMGZo2bZrjd6vV6rGgxMw2AADcz6MhKTo6WomJiZWOdevWTWvWrJEknT17VjNnzlRaWppGjRolSerZs6d2796tBQsW1DoktWzZUoGBgTp+/Hil48ePH1fr1q2rfU9ISIhCQkJqdT9XKimv0PeniiXxuA0AAHfy6OO2QYMG6cCBA5WOZWRkKD4+XpJUVlamsrIyBQRULjMwMFA2m63W9zWbzerbt682bNjgOGaz2bRhwwYNGDCg1tetD1kni1RhMxQWEqTW4cxsAwDAXTw6kvTAAw9o4MCBmjt3rm677TZt375dS5cu1dKlSyVJ4eHhGjJkiKZPny6LxaL4+Hht3rxZq1at0gsvvOC4Tk5OjnJycnTw4EFJ0p49exQWFqa2bduqefPmkqShQ4fqlltu0aRJkySdX1ogNTVVycnJ+s1vfqOFCxeqqKhId911Vz3/LTjnwnYkoTKZ2LMNAAB38WhI6tevn9LS0jRjxgw99dRTat++vRYuXKiUlBTHOatXr9aMGTOUkpKivLw8xcfHa86cOZowYYLjnCVLlujJJ590/G5v6l6+fLnGjRsnSTp06JBOnjzpOGfMmDE6ceKEHn/8ceXk5Kh379768MMPqzRzexuatgEAqB8eXSfJl3lqnaS/vL5T6787rlm/T9Sfrmpfb/cFAMAf+Mw6SXBeJmskAQBQLwhJPuRcWYW+P8WebQAA1AdCkg85fKJINkMKbxSkqDDPL0cAAIA/IyT5kMzcC03bzGwDAMC9CEk+xL7SdufWPGoDAMDdCEk+xLGxbRRN2wAAuBshyYewRhIAAPWHkOQjzpVV6Ie883u2sbEtAADuR0jyEQdzC2UYUrPGwWoZavZ0OQAA+D1Cko+wz2xLYGYbAAD1gpDkIw7ksNI2AAD1iZDkI2jaBgCgfhGSfESG/XFbFCEJAID6QEjyAcWl5fox76wkHrcBAFBfCEk+4GDu+X6kFk3MahHKnm0AANQHQpIPsK+0ncAoEgAA9YaQ5ANo2gYAoP4RknyAfWNbVtoGAKD+EJJ8ABvbAgBQ/whJXq6opFxHz9hntjGSBABAfSEkebnMn2e2tQwNUbMm7NkGAEB9ISR5uQxH0zaP2gAAqE+EJC/HzDYAADyDkOTlWCMJAADPICR5OUaSAADwDEKSFys4V6bs/HOSpM5sbAsAQL0iJHkx+8y2qLAQRTQO9nA1AAA0LIQkL8ajNgAAPIeQ5MUcK20TkgAAqHeEJC/GGkkAAHgOIcmLZTqm/zOSBABAfSMkean8s2XKsZ6f2cYaSQAA1D9Ckpc6mHv+UVt0RCOFN2JmGwAA9Y2Q5KUyeNQGAIBHEZK81IGcn5u2o3jUBgCAJxCSvFRmLmskAQDgSYQkL8XGtgAAeBYhyQudKS7ViYISSfQkAQDgKYQkL2QfRWrT1KLQkCAPVwMAQMNESPJC9pW2edQGAIDnEJK8EBvbAgDgeYQkL+Ro2mb6PwAAHkNI8kJM/wcAwPMISV4mr6hUJwtLJUmdGEkCAMBjCElext60HdvMoibMbAMAwGMISV6Gpm0AALwDIcnLFJSUq1FwANP/AQDwMJNhGIani/BFVqtVERERys/PV3h4uEuvbbMZKim3yWIOdOl1AQBo6Jz5/mYkyQsFBJgISAAAeBghCQAAoBqEJAAAgGoQkgAAAKpBSAIAAKgGIQkAAKAahCQAAIBqEJIAAACqQUgCAACoBiEJAACgGoQkAACAahCSAAAAqkFIAgAAqAYhCQAAoBpBni7AVxmGIUmyWq0ergQAANSU/Xvb/j1+KYSkWiooKJAkxcXFebgSAADgrIKCAkVERFzyHJNRkyiFKmw2m7KzsxUWFiaTyeTSa1utVsXFxenHH39UeHi4S6/tDfh8vs/fPyOfz/f5+2fk89WeYRgqKChQTEyMAgIu3XXESFItBQQEKDY21q33CA8P98t/+e34fL7P3z8jn8/3+ftn5PPVzuVGkOxo3AYAAKgGIQkAAKAahCQvFBISotmzZyskJMTTpbgFn8/3+ftn5PP5Pn//jHy++kHjNgAAQDUYSQIAAKgGIQkAAKAahCQAAIBqEJIAAACqQUjyIlu2bNENN9ygmJgYmUwmrV271tMludS8efPUr18/hYWFKSoqSjfffLMOHDjg6bJcZvHixerZs6dj8bMBAwZo3bp1ni7LbebPny+TyaSpU6d6uhSXeeKJJ2QymSr9dO3a1dNludTRo0c1duxYtWjRQhaLRUlJSdq5c6eny3KZdu3aVflnaDKZNHHiRE+X5hIVFRWaNWuW2rdvL4vFoo4dO+rpp5+u0T5kvqKgoEBTp05VfHy8LBaLBg4cqB07dnikFlbc9iJFRUXq1auXxo8fr1tvvdXT5bjc5s2bNXHiRPXr10/l5eWaOXOmhg8frn379qlJkyaeLq/OYmNjNX/+fCUkJMgwDK1cuVI33XSTvv76a3Xv3t3T5bnUjh079Morr6hnz56eLsXlunfvrk8++cTxe1CQ//xn8vTp0xo0aJCuueYarVu3TpGRkcrMzFSzZs08XZrL7NixQxUVFY7f9+7dq9/97ncaPXq0B6tynWeffVaLFy/WypUr1b17d+3cuVN33XWXIiIiNGXKFE+X5xJ333239u7dq9dff10xMTF64403NGzYMO3bt09t2rSp32IMeCVJRlpamqfLcKvc3FxDkrF582ZPl+I2zZo1M1599VVPl+FSBQUFRkJCgvHxxx8bQ4YMMe6//35Pl+Qys2fPNnr16uXpMtzm4YcfNq666ipPl1Gv7r//fqNjx46GzWbzdCkuMWrUKGP8+PGVjt16661GSkqKhypyreLiYiMwMND44IMPKh2/4oorjEcffbTe6+FxGzwmPz9fktS8eXMPV+J6FRUVWr16tYqKijRgwABPl+NSEydO1KhRozRs2DBPl+IWmZmZiomJUYcOHZSSkqIjR454uiSXSU9PV3JyskaPHq2oqCj16dNHy5Yt83RZblNaWqo33nhD48ePd/lG5J4ycOBAbdiwQRkZGZKkb775Rp999plGjhzp4cpco7y8XBUVFWrUqFGl4xaLRZ999lm91+M/48jwKTabTVOnTtWgQYPUo0cPT5fjMnv27NGAAQN07tw5hYaGKi0tTYmJiZ4uy2VWr16tXbt2eaw/wN369++vFStWqEuXLjp27JiefPJJ/fa3v9XevXsVFhbm6fLq7PDhw1q8eLGmTZummTNnaseOHZoyZYrMZrNSU1M9XZ7LrV27VmfOnNG4ceM8XYrLPPLII7JareratasCAwNVUVGhOXPmKCUlxdOluURYWJgGDBigp59+Wt26dVOrVq301ltv6YsvvlCnTp3qv6B6H7tCjcjPH7dNmDDBiI+PN3788UdPl+JSJSUlRmZmprFz507jkUceMVq2bGl89913ni7LJY4cOWJERUUZ33zzjeOYvz1uu9jp06eN8PBwv3lkGhwcbAwYMKDSscmTJxtXXnmlhypyr+HDhxu///3vPV2GS7311ltGbGys8dZbbxnffvutsWrVKqN58+bGihUrPF2ayxw8eNAYPHiwIckIDAw0+vXrZ6SkpBhdu3at91oYSUK9mzRpkj744ANt2bJFsbGxni7Hpcxms+P/7fTt21c7duzQiy++qFdeecXDldXdV199pdzcXF1xxRWOYxUVFdqyZYteeukllZSUKDAw0IMVul7Tpk3VuXNnHTx40NOluER0dHSVkc1u3bppzZo1HqrIfX744Qd98sknevfddz1diktNnz5djzzyiG6//XZJUlJSkn744QfNmzfPb0YDO3bsqM2bN6uoqEhWq1XR0dEaM2aMOnToUO+10JOEemMYhiZNmqS0tDR9+umnat++vadLcjubzaaSkhJPl+ESQ4cO1Z49e7R7927HT3JyslJSUrR7926/C0iSVFhYqEOHDik6OtrTpbjEoEGDqiy7kZGRofj4eA9V5D7Lly9XVFSURo0a5elSXKq4uFgBAZW/ugMDA2Wz2TxUkfs0adJE0dHROn36tNavX6+bbrqp3mtgJMmLFBYWVvp/rFlZWdq9e7eaN2+utm3berAy15g4caLefPNNvffeewoLC1NOTo4kKSIiQhaLxcPV1d2MGTM0cuRItW3bVgUFBXrzzTe1adMmrV+/3tOluURYWFiV/rEmTZqoRYsWftNX9tBDD+mGG25QfHy8srOzNXv2bAUGBuqOO+7wdGku8cADD2jgwIGaO3eubrvtNm3fvl1Lly7V0qVLPV2aS9lsNi1fvlypqal+tYSDJN1www2aM2eO2rZtq+7du+vrr7/WCy+8oPHjx3u6NJdZv369DMNQly5ddPDgQU2fPl1du3bVXXfdVf/F1PsDPvyqjRs3GpKq/KSmpnq6NJeo7rNJMpYvX+7p0lxi/PjxRnx8vGE2m43IyEhj6NChxkcffeTpstzK33qSxowZY0RHRxtms9lo06aNMWbMGOPgwYOeLsul3n//faNHjx5GSEiI0bVrV2Pp0qWeLsnl1q9fb0gyDhw44OlSXM5qtRr333+/0bZtW6NRo0ZGhw4djEcffdQoKSnxdGku8/bbbxsdOnQwzGaz0bp1a2PixInGmTNnPFKLyTD8aJlOAAAAF6EnCQAAoBqEJAAAgGoQkgAAAKpBSAIAAKgGIQkAAKAahCQAAIBqEJIAAACqQUgCAACoBiEJAGqoXbt2WrhwoafLAFBPCEkAvNK4ceN08803S5KuvvpqTZ06td7uvWLFCjVt2rTK8R07duiee+6ptzoAeJZ/7fwHAJdQWloqs9lc6/dHRka6sBoA3o6RJABebdy4cdq8ebNefPFFmUwmmUwmff/995KkvXv3auTIkQoNDVWrVq30xz/+USdPnnS89+qrr9akSZM0depUtWzZUtddd50k6YUXXlBSUpKaNGmiuLg43XfffSosLJQkbdq0SXfddZfy8/Md93viiSckVX3cduTIEd10000KDQ1VeHi4brvtNh0/ftzx+hNPPKHevXvr9ddfV7t27RQREaHbb79dBQUFjnPeeecdJSUlyWKxqEWLFho2bJiKiorc9LcJwBmEJABe7cUXX9SAAQP05z//WceOHdOxY8cUFxenM2fO6Nprr1WfPn20c+dOffjhhzp+/Lhuu+22Su9fuXKlzGaztm7dqiVLlkiSAgIC9I9//EPfffedVq5cqU8//VR//etfJUkDBw7UwoULFR4e7rjfQw89VKUum82mm266SXl5edq8ebM+/vhjHT58WGPGjKl03qFDh7R27Vp98MEH+uCDD7R582bNnz9fknTs2DHdcccdGj9+vPbv369Nmzbp1ltvFfuOA96Bx20AvFpERITMZrMaN26s1q1bO46/9NJL6tOnj+bOnes49q9//UtxcXHKyMhQ586dJUkJCQn629/+Vumav+xvateunZ555hlNmDBB//znP2U2mxURESGTyVTpfhfbsGGD9uzZo6ysLMXFxUmSVq1ape7du2vHjh3q16+fpPNhasWKFQoLC5Mk/fGPf9SGDRs0Z84cHTt2TOXl5br11lsVHx8vSUpKSqrD3xYAV2IkCYBP+uabb7Rx40aFhoY6frp27Srp/OiNXd++fau895NPPtHQoUPVpk0bhYWF6Y9//KNOnTql4uLiGt9///79iouLcwQkSUpMTFTTpk21f/9+x7F27do5ApIkRUdHKzc3V5LUq1cvDR06VElJSRo9erSWLVum06dP1/wvAYBbEZIA+KTCwkLdcMMN2r17d6WfzMxMDR482HFekyZNKr3v+++/1+9//3v17NlTa9as0VdffaWXX35Z0vnGblcLDg6u9LvJZJLNZpMkBQYG6uOPP9a6deuUmJioRYsWqUuXLsrKynJ5HQCcR0gC4PXMZrMqKioqHbviiiv03XffqV27durUqVOln4uD0S999dVXstlsev7553XllVeqc+fOys7Ovuz9LtatWzf9+OOP+vHHHx3H9u3bpzNnzigxMbHGn81kMmnQoEF68skn9fXXX8tsNistLa3G7wfgPoQkAF6vXbt2+vLLL/X999/r5MmTstlsmjhxovLy8nTHHXdox44dOnTokNavX6+77rrrkgGnU6dOKisr06JFi3T48GG9/vrrjobuX96vsLBQGzZs0MmTJ6t9DDds2DAlJSUpJSVFu3bt0vbt23XnnXdqyJAhSk5OrtHn+vLLLzV37lzt3LlTR44c0bvvvqsTJ06oW7duzv0FAXALQhIAr/fQQw8pMDBQiYmJioyM1JEjRxQTE6OtW7eqoqJCw4cPV1JSkqZOnaqmTZsqIODX/9PWq1cvvfDCC3r22WfVo0cP/fvf/9a8efMqnTNw4EBNmDBBY8aMUWRkZJXGb+n8CNB7772nZs2aafDgwRo2bJg6dOigt99+u8afKzw8XFu2bNH111+vzp0767HHHtPzzz+vkSNH1vwvB4DbmAzmmgIAAFTBSBIAAEA1CEkAAADVICQBAABUg5AEAABQDUISAABANQhJAAAA1SAkAQAAVIOQBAAAUA1CEgAAQDUISQAAANUgJAEAAFTj/wMT8iCO58O+dAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Permissible current convergence plot\n",
    "I = I[:i+1]\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(list(range(1,i+2)),I)\n",
    "plt.xlabel('Iterations')\n",
    "plt.ylabel('Permissible current [A]')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 730,
   "id": "b5519496-3d0e-462c-aa1b-7c3ac8e4d9bf",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The permissible current rating of one conductor is equal to 682.814465 A\n",
      "The permissible current rating of one conductor is equal to   680 A\n"
     ]
    }
   ],
   "source": [
    "# Permissible current rounding\n",
    "print(\"The permissible current rating of one conductor is equal to %10f A\" % I[i]) \n",
    "\n",
    "if I[i] < 200:\n",
    "    I_rating = int(I[i])\n",
    "elif I[i] < 500:\n",
    "    I_rating = int(I[i]/5)*5\n",
    "else:\n",
    "    I_rating = int(I[i]/10)*10\n",
    "     \n",
    "print(\"The permissible current rating of one conductor is equal to %5d A\" % I_rating)"
   ]
  }
 ],
 "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.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}