cigreTB880_case01ecl.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",
    "\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": 1,
   "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. Introductory Case study 0\n",
    "\n",
    "<div>\n",
    "<img src=\"./img/img1_cs0.png\" width=\"150\"/>\n",
    "</div>\n",
    "\n",
    "**Input data:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "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": 3,
   "id": "735eb1b1-ac83-46fc-b39c-da838c8f339a",
   "metadata": {},
   "outputs": [],
   "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] Diammeter after the layer\n",
    "dia[0] = 30.3 #[mm] Conductor diameter (630 mm2)\n",
    "\n",
    "for i in range(1,len(dia)):\n",
    "    dia[i] = dia[i-1] + 2*thk[i]\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": 4,
   "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.6 CASE #0-1: Possible variations\n",
    "\n",
    "### 4.6.4 Variant with non-neglected eddy-current losses \n",
    "\n",
    "This option was adopted by the CIGRE WG B1.56, not to neglect the eddy-current losses when the metal sheaths are bonded at both ends of an electrical section (Guidance Point 6)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1feb518f-360f-4cdc-ab3f-dc7103b9db03",
   "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": 5,
   "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": 6,
   "id": "7fab0e9f-c7a2-4273-8192-1e2b1914dd64",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "X = 5.0403313985E-05 Ω/m\n"
     ]
    }
   ],
   "source": [
    "# 4.4.2 Reactance\n",
    "\n",
    "omega = 2*mt.pi*f #[rad/s] Angular frequency\n",
    "s = 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",
    "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": 7,
   "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": 8,
   "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": 9,
   "id": "1fa6cb36-8d13-4b68-bb91-8cf052a9de24",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "T3 = 0.0867193748 Km/W\n"
     ]
    }
   ],
   "source": [
    "# 4.4.5 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 = 1.6*T3 #[Km/W] From IEC 60287-2-1 Section 4.2.4.3.2\n",
    "print(\"T3 = %.10f Km/W\" % T3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "2bd97694-089a-49fe-b90a-925ea00a3eb7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "T4 = 1.5946928925 Km/W\n"
     ]
    }
   ],
   "source": [
    "# 4.4.6 External thermal resistance\n",
    "\n",
    "L = 1000 #[mm]\n",
    "u = 2*L/dia[5]\n",
    "rho_T_ex = 1 #[Km/W]\n",
    "\n",
    "T4 = 1.5/mt.pi*rho_T_ex*(mt.log(2*u)-0.630) #[Km/W] From IEC 60287-2-1 Section 4.2.4.3.2\n",
    "print(\"T4 = %.10f Km/W\" % T4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "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**2*tan_delta #[W/m] From IEC 60287-1-1 Section 2.2\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": 12,
   "id": "871b93a1-733c-4503-abc3-37f9555caa4c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration:  1\n",
      "R = 3.9521526380E-05 Ohm/m\n",
      "Rs = 2.0727239574E-04 Ohm/m\n",
      "λ1 = 0.3654736398\n",
      "I = 803.3635402401 A\n",
      "Wc = 25.5069155968 W/m\n",
      "Ws = 9.3221052841 W/m\n",
      "θj = 76.1557692286 °C\n",
      "θs = 79.2095190893 °C\n",
      "θc = 90.0000000000 °C\n",
      "Iteration:  2\n",
      "R = 3.9521526380E-05 Ohm/m\n",
      "Rs = 2.0674067176E-04 Ohm/m\n",
      "λ1 = 0.3662997154\n",
      "I = 803.1581818270 A\n",
      "Wc = 25.4938769413 W/m\n",
      "Ws = 9.3383998691 W/m\n",
      "θj = 76.1609614361 °C\n",
      "θs = 79.2149936490 °C\n",
      "θc = 90.0000000000 °C\n",
      "Iteration:  3\n",
      "R = 3.9521526380E-05 Ohm/m\n",
      "Rs = 2.0674435427E-04 Ohm/m\n",
      "λ1 = 0.3662939825\n",
      "I = 803.1596064582 A\n",
      "Wc = 25.4939673828 W/m\n",
      "Ws = 9.3382868431 W/m\n",
      "θj = 76.1609254209 °C\n",
      "θs = 79.2149556752 °C\n",
      "θc = 90.0000000000 °C\n",
      "Iteration:  4\n",
      "R = 3.9521526380E-05 Ohm/m\n",
      "Rs = 2.0674432873E-04 Ohm/m\n",
      "λ1 = 0.3662940223\n",
      "I = 803.1595965765 A\n",
      "Wc = 25.4939667554 W/m\n",
      "Ws = 9.3382876271 W/m\n",
      "θj = 76.1609256707 °C\n",
      "θs = 79.2149559386 °C\n",
      "θc = 90.0000000000 °C\n",
      "Iteration:  5\n",
      "R = 3.9521526380E-05 Ohm/m\n",
      "Rs = 2.0674432891E-04 Ohm/m\n",
      "λ1 = 0.3662940220\n",
      "I = 803.1595966450 A\n",
      "Wc = 25.4939667598 W/m\n",
      "Ws = 9.3382876217 W/m\n",
      "θj = 76.1609256689 °C\n",
      "θs = 79.2149559368 °C\n",
      "θc = 90.0000000000 °C\n",
      "Iteration:  6\n",
      "R = 3.9521526380E-05 Ohm/m\n",
      "Rs = 2.0674432891E-04 Ohm/m\n",
      "λ1 = 0.3662940220\n",
      "I = 803.1595966445 A\n",
      "Wc = 25.4939667597 W/m\n",
      "Ws = 9.3382876217 W/m\n",
      "θj = 76.1609256690 °C\n",
      "θs = 79.2149559368 °C\n",
      "θc = 90.0000000000 °C\n",
      "Iteration:  7\n",
      "R = 3.9521526380E-05 Ohm/m\n",
      "Rs = 2.0674432891E-04 Ohm/m\n",
      "λ1 = 0.3662940220\n",
      "I = 803.1595966445 A\n",
      "Wc = 25.4939667597 W/m\n",
      "Ws = 9.3382876217 W/m\n",
      "θj = 76.1609256690 °C\n",
      "θs = 79.2149559368 °C\n",
      "θc = 90.0000000000 °C\n",
      "Convergence reached after  7 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",
    "    # 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",
    "    \n",
    "    # 4.6.3.1 Losses caused by eddy currents (Dielectric losses taken into account)\n",
    "    rho_s = rho_s_20*(1+alpha_s_20*(theta_s[i]-20)) #[Ohm/m]\n",
    "    \n",
    "    beta_1 = mt.sqrt(4*mt.pi*omega/(1e7*rho_s)) # From IEC 60287-1-1 Section 2.3.6.1\n",
    "    m = (omega/R_s[i])*1e-7\n",
    "    D_s = dia[4] # Diameter after metal sheath\n",
    "    g_s = 1 + (t_s/D_s)**1.74*(beta_1*D_s*1e-3-1.6) # From IEC 60287-1-1 Section 2.3.6.1\n",
    "    # For cables in trefoil formation\n",
    "    lambda_0 = 3*(m**2/(1+m**2))*(d/(2*s))**2 # From IEC 60287-1-1 Section 2.3.6.1\n",
    "    Delta_1 = (1.14*m**2.45+0.33)*(d/(2*s))**(0.92*m+1.66)\n",
    "    Delta_2 = 0\n",
    "    \n",
    "    lambda_1_prime2 = R_s[i]/R[i]*(g_s*lambda_0*(1+Delta_1+Delta_2)+((beta_1*t_s)**4/12e12))\n",
    "    \n",
    "    # 4.6.4.1 Reduction of the eddy currents losses by the circulating currents \n",
    "    M = R_s[i]/X # From IEC 60287-1-1 Section 2.3.5\n",
    "    N = M # For cables in trefoil formation\n",
    "    F = (4*M**2*N**2+(M+N)**2)/(4*(M**2+1)*(N**2+1))\n",
    "    \n",
    "    lambda_1_prime2 = F*lambda_1_prime2\n",
    "    \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",
    "        theta_s[i+1] = theta_s[i] # Update for next iteration\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": 13,
   "id": "5957dc03-af4f-4841-870a-9018df3a14ad",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlIAAAGwCAYAAABiu4tnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABJQUlEQVR4nO3deXxU9b3/8fdkJSErktUsBIwiECDBVrYLKBIFpbhUBFHhcm1RwEqtgloXqkIElYuSCtVaRJQr91eFa7UsESGKiAJhl0UEkgCJIZiNQBKSnN8fMFNSEGfChDPL6/l4zMPMmTPnfDL6aN79ns98jsUwDEMAAABwmI/ZBQAAALgrghQAAEAzEaQAAACaiSAFAADQTAQpAACAZiJIAQAANBNBCgAAoJn8zC7AkzU2NurIkSMKDQ2VxWIxuxwAAGAHwzBUVVWl+Ph4+fhceM2JINWCjhw5osTERLPLAAAAzVBYWKiEhIQL7kOQakGhoaGSTv+LCAsLM7kaAABgj8rKSiUmJtr+jl8IQaoFWS/nhYWFEaQAAHAz9rTl0GwOAADQTAQpAACAZiJIAQAANBNBCgAAoJkIUgAAAM1EkAIAAGgmghQAAEAzEaQAAACaiSAFAADQTAQpAACAZiJIAQAANBNBCgAAoJkIUm7IMAwVHDuhw+UnzS4FAACvRpByQ9M+2aV+L63W218eMLsUAAC8GkHKDXWMC5Mk5RWUm1sIAABejiDlhnokR0qSth+uUG19g8nVAADgvQhSbqjdZcFq0zpAdfWN2nmk0uxyAADwWgQpN2SxWJSeGCFJyssvM7cYAAC8GEHKTWWcuby3mT4pAABMQ5ByUxlJp4NUXgErUgAAmIUg5aa6JYbL18eioooaHWGeFAAApiBIuangAD91jA2VxKoUAABmIUi5MdvlvfxycwsBAMBLEaTcmHWeFCtSAACYgyDlxqwrUjuPVKjmFIM5AQC41AhSbiyxTZDahgToVIOhHYcrzC4HAACvQ5ByYxaLhTEIAACYiCDl5qyDOWk4BwDg0iNIuTnritSmgjIZhmFyNQAAeBeClJvrmhAuPx+LjlbV6lAZgzkBALiUCFJurpW/rzrHh0miTwoAgEuNIOUB0pO4gTEAAGYwNUjV19frqaeeUkpKioKCgtS+fXs999xzamxstO1jGIamTp2q+Ph4BQUFacCAAdq5c2eT44wbN04dOnRQUFCQoqKiNGzYMO3evfuC5547d666du2qsLAwhYWFqVevXlq2bFmTfcaMGSOLxdLk0bNnT+d9AE5ibTjflM+KFAAAl5KpQWrGjBmaN2+esrOztWvXLs2cOVMvvfSS5syZY9tn5syZmjVrlrKzs7VhwwbFxsZq0KBBqqqqsu3To0cPzZ8/X7t27dKKFStkGIYyMzPV0PDTQyoTEhL04osvauPGjdq4caOuv/56DRs27JyQdtNNN6moqMj2+Oc//+n8D+IiZSRFSJJ2FVXqZB2DOQEAuFQsholf9brlllsUExOjt956y7btjjvuUHBwsBYuXCjDMBQfH69JkyZpypQpkqTa2lrFxMRoxowZGjdu3HmPu23bNnXr1k379u1Thw4d7K6nTZs2eumll/Rf//Vfkk6vSJWXl2vp0qXN+v0qKysVHh6uiooKhYWFNesY9jAMQz2zVumHylot/m1PXdv+shY7FwAAns6Rv9+mrkj17dtXq1at0t69eyVJW7du1dq1azVkyBBJ0oEDB1RcXKzMzEzbewIDA9W/f3+tW7fuvMesrq7W/PnzlZKSosTERLvqaGho0Pvvv6/q6mr16tWryWtr1qxRdHS0rrzySv3mN79RSUnJTx6ntrZWlZWVTR6XwtmDOTfRcA4AwCVjapCaMmWKRo4cqY4dO8rf31/p6emaNGmSRo4cKUkqLi6WJMXExDR5X0xMjO01q9dff10hISEKCQnR8uXLlZOTo4CAgAuef/v27QoJCVFgYKAeeOABLVmyRJ06dbK9PnjwYL333nv67LPP9Morr2jDhg26/vrrVVtbe97jZWVlKTw83PawN8g5g23COYM5AQC4ZEwNUosXL9a7776rRYsWKS8vTwsWLNDLL7+sBQsWNNnPYrE0eW4YxjnbRo0apc2bNys3N1epqakaPny4ampqLnj+q666Slu2bNH69ev14IMPavTo0fr2229tr9911126+eab1aVLFw0dOlTLli3T3r179cknn5z3eE888YQqKipsj8LCQkc+jouSkRwhSdrMYE4AAC4ZPzNP/thjj+nxxx/XiBEjJElpaWnKz89XVlaWRo8erdjYWEmnV6bi4uJs7yspKTlnlcq6CpSamqqePXsqMjJSS5Yssa1unU9AQICuuOIKSdI111yjDRs26NVXX9Vf/vKX8+4fFxen5ORkfffdd+d9PTAwUIGBgfZ/AE7UOT5cAb4+OlZdp4IfTyj5stam1AEAgDcxdUXqxIkT8vFpWoKvr69t/EFKSopiY2OVk5Nje72urk65ubnq3bv3BY9tGMZPXoJr7nuOHTumwsLCJqHOVbTy91Xny083xDEGAQCAS8PUIDV06FBNmzZNn3zyiQ4ePKglS5Zo1qxZuu222ySdvqQ3adIkTZ8+XUuWLNGOHTs0ZswYBQcH6+6775Yk7d+/X1lZWdq0aZMKCgr01Vdfafjw4QoKCrI1rUvSwIEDlZ2dbXv+5JNP6osvvtDBgwe1fft2/fGPf9SaNWs0atQoSdLx48f16KOP6quvvtLBgwe1Zs0aDR06VG3btrXV52psfVI0nAMAcEmYemlvzpw5evrppzV+/HiVlJQoPj5e48aN0zPPPGPbZ/LkyTp58qTGjx+vsrIyXXvttVq5cqVCQ0MlSa1atdIXX3yh2bNnq6ysTDExMerXr5/WrVun6Oho23G+//57lZaW2p7/8MMPuvfee1VUVKTw8HB17dpVy5cv16BBgySdXhnbvn273nnnHZWXlysuLk7XXXedFi9ebDu3q8lIitRbOkDDOQAAl4ipc6Q83aWaI2VVVHFSvbI+k49F2j71RrUONDUnAwDgltxmjhScKy48SPHhrdRoSFsLy80uBwAAj0eQ8jDpyfRJAQBwqRCkPMy/Gs7LzS0EAAAvQJDyMD3OWpGi/Q0AgJZFkPIwneLCFOjno/ITp7S/tNrscgAA8GgEKQ8T4OejtMvDJUl5DOYEAKBFEaQ8UEYyfVIAAFwKBCkPZG0438w39wAAaFEEKQ+UkRwhSdrzQ5Uqa06ZWwwAAB6MIOWBokNbKSEySAaDOQEAaFEEKQ9lmyfFffcAAGgxBCkP1YMJ5wAAtDiClIf614TzMjU2MpgTAICWQJDyUB3jQtXK30dVNfX6/uhxs8sBAMAjEaQ8lL+vj7olREji8h4AAC2FIOXBbIM5aTgHAKBFEKQ8mLVPahMrUgAAtAiClAdLT4qQJO0rOa6KEwzmBADA2QhSHqxtSKDaXRYsSdpcyKoUAADORpDycP8azEmQAgDA2QhSHi7dNpiz3NxCAADwQAQpD5dxpk9qS2G5GhjMCQCAUxGkPNxVMaFqHeCr47X1+q6kyuxyAADwKAQpD+fn66NuiRGSpE30SQEA4FQEKS/wr4bzcnMLAQDAwxCkvEBGcoQkaTODOQEAcCqClBdITzy9IrW/tFo/VteZXA0AAJ6DIOUFIlsHqH1Ua0msSgEA4EwEKS9h65MiSAEA4DQEKS9BwzkAAM5HkPISPc5MON96qFz1DY0mVwMAgGcgSHmJ1OgQhQb66URdg3YXM5gTAABnIEh5CR8fi7qfuV0MDecAADgHQcqLpCdxA2MAAJyJIOVFrH1SfHMPAADnIEh5ke5n7rmXf+yESo/XmlsMAAAegCDlRcKD/JUaHSJJyuMGxgAAXDSClJfJoE8KAACnIUh5GfqkAABwHoKUl8lIjpAkbTtUrlMM5gQA4KIQpLxM+7YhCmvlp5pTjdpVVGl2OQAAuDWClJfx8bEow3p5j4ZzAAAuCkHKC9FwDgCAcxCkvJA1SG1iRQoAgItCkPJC3RLDZbFIh8tPqqSyxuxyAABwWwQpLxTayl9XxYRKYgwCAAAXgyDlpawN51zeAwCg+QhSXoqGcwAALh5ByktlJEVIkrYfrlBdPYM5AQBoDoKUl0pp21qRwf6qq2/UziMVZpcDAIBbIkh5KYvFwhgEAAAuEkHKi1kbzjfTJwUAQLMQpLxY+pk+KUYgAADQPAQpL9YtIUK+PhYVVdToSPlJs8sBAMDtEKS8WOtAP3WMZTAnAADNRZDycrZ5Uvnl5hYCAIAbIkh5uYzkCEmsSAEA0BwEKS/XI6mNJGnnkQrVnGowuRoAANwLQcrLJbYJUtuQAJ1qMLTjMIM5AQBwBEHKy1ksFqXb7rvH5T0AABxBkAIN5wAANBNBCupxZsL5poIyGYZhcjUAALgPghTUNSFcfj4WHa2q1aEyBnMCAGAvP3t2atOmjUMHtVgsysvLU3JycrOKwqXVyt9XneLDtO1QhfIKypTYJtjskgAAcAt2Bany8nLNnj1b4eHhP7uvYRgaP368Ghr4Kr07yUiK1LZDFdpcUK5h3S83uxwAANyCXUFKkkaMGKHo6Gi79n3ooYeaXRDMkZEcqbfXHeSbewAAOMCuINXY2OjQQauqqppVDMyTkRQhSfr2SKVO1jUoKMDX3IIAAHADTms2b2ho0NKlS511OFxil0cEKTo0UPWNhrYdKje7HAAA3MJFB6ndu3dr8uTJio+P1/Dhw51RE0xgsVj+NU+qoNzcYgAAcBPNClLV1dX629/+pj59+qhz587Ky8vTtGnTdOTIEWfXh0vIOk+KPikAAOxjd7O5JH311Vf661//qv/93/9VamqqRo0apa+//lqvvfaaOnXq1FI14hLJSI6QJOXlnx7MabFYzC0IAAAXZ/eKVKdOnTRy5EjFxMTo66+/Vl5env7whz9c1B/b+vp6PfXUU0pJSVFQUJDat2+v5557rklzu2EYmjp1quLj4xUUFKQBAwZo586dTY4zbtw4dejQQUFBQYqKitKwYcO0e/fuC5577ty56tq1q8LCwhQWFqZevXpp2bJlTfax59yepHN8uPx9LTpWXaeCH0+YXQ4AAC7P7iC1b98+9evXT9ddd52uvvpqp5x8xowZmjdvnrKzs7Vr1y7NnDlTL730kubMmWPbZ+bMmZo1a5ays7O1YcMGxcbGatCgQU2+GdijRw/Nnz9fu3bt0ooVK2QYhjIzMy84yyohIUEvvviiNm7cqI0bN+r666/XsGHDmgQle87tSVr5+6rL5adnhXF5DwAAOxh2OnTokPHCCy8YHTp0MOLj440//OEPRl5enuHv72/s3LnT3sM0cfPNNxtjx45tsu3222837rnnHsMwDKOxsdGIjY01XnzxRdvrNTU1Rnh4uDFv3ryfPO7WrVsNSca+ffscqicyMtL461//elHnPltFRYUhyaioqHCoDjM994+dRvKUj42nlmw3uxQAAEzhyN9vu1ekLr/8cv3xj3/Uvn37tHDhQhUXF6tPnz6qr6/X22+/rb179zoc4vr27atVq1bZ3rt161atXbtWQ4YMkSQdOHBAxcXFyszMtL0nMDBQ/fv317p16857zOrqas2fP18pKSlKTEy0q46Ghga9//77qq6uVq9evZp97traWlVWVjZ5uBvrN/c25bMiBQDAz2nWt/auv/56vfvuuyoqKlJ2drY+++wzdezYUV27dnXoOFOmTNHIkSPVsWNH+fv7Kz09XZMmTdLIkSMlScXFxZKkmJiYJu+LiYmxvWb1+uuvKyQkRCEhIVq+fLlycnIUEBBwwfNv375dISEhCgwM1AMPPKAlS5bYmuYdObdVVlaWwsPDbQ97g5wrsTac7y6uVHVtvbnFAADg4i5qjlR4eLjGjx+vjRs3Ki8vTwMGDHDo/YsXL9a7776rRYsWKS8vTwsWLNDLL7+sBQsWNNnv3xvajfN8o2zUqFHavHmzcnNzlZqaquHDh6umpuaC57/qqqu0ZcsWrV+/Xg8++KBGjx6tb7/91uFzWz3xxBOqqKiwPQoLCy94flcUFx6k+PBWajSkrQzmBADgghwaf3Ah3bt312uvvebQex577DE9/vjjGjFihCQpLS1N+fn5ysrK0ujRoxUbGyvp9OpQXFyc7X0lJSXnrBRZV4FSU1PVs2dPRUZGasmSJbbVrfMJCAjQFVdcIUm65pprtGHDBr366qv6y1/+4tC5rQIDAxUYGOjQZ+CK0pMjdWRbkfLyy9S7Q1uzywEAwGXZtSKVkZGhsjL7e2b69u2rw4cP/+x+J06ckI9P0xJ8fX1t4w9SUlIUGxurnJwc2+t1dXXKzc1V7969L3hswzBUW1trd83//p6LObe7Y8I5AAD2sWtFasuWLdq6davatGlj10G3bNliV4gZOnSopk2bpqSkJHXu3FmbN2/WrFmzNHbsWEmnL6tNmjRJ06dPV2pqqlJTUzV9+nQFBwfr7rvvliTt379fixcvVmZmpqKionT48GHNmDFDQUFBtqZ1SRo4cKBuu+02TZw4UZL05JNPavDgwUpMTFRVVZXef/99rVmzRsuXL7f73J7KegPjzQUM5gQA4ELsvrQ3cOBAGYZh1772/uGdM2eOnn76aY0fP14lJSWKj4/XuHHj9Mwzz9j2mTx5sk6ePKnx48errKxM1157rVauXKnQ0FBJUqtWrfTFF19o9uzZKisrU0xMjPr166d169YpOjradpzvv/9epaWltuc//PCD7r33XhUVFSk8PFxdu3bV8uXLNWjQILvP7ak6x4crwM9HZSdO6UBptdpHhZhdEgAALsli2JGO8vPzHT5wQkKCfH19m1WUp6isrFR4eLgqKioUFhZmdjkO+fXcddqYX6aXft1Vd17jft8+BACguRz5+23XilRycrJTCoP7yEiO1Mb8MuUVlBOkAAD4CRc1/gCe6+w+KQAAcH4EKZyX9Zt7e36oUlXNKZOrAQDANRGkcF7RYa2UEBkkw5C2FJabXQ4AAC6JIIWfZJsnlV9ubiEAALgoh4NU+/btdezYsXO2l5eXq3379k4pCq7B2ieVR58UAADn5XCQOnjwoBoaGs7ZXltba9c0c7iPHsmnB7BuLihTY6N9M8QAAPAmdg/k/Oijj2w/r1ixQuHh4bbnDQ0NWrVqldq1a+fU4mCujnGhauXvo8qaen1/9LhSYzx7ECkAAI6yO0jdeuutkk5PLR89enST1/z9/dWuXTu98sorTi0O5vL39VHXhAh9c+BH5RWUEaQAAPg3dl/aa2xsVGNjo5KSklRSUmJ73tjYqNraWu3Zs0e33HJLS9YKE9BwDgDAT7N7RcrqwIEDLVEHXFSP5DNBioZzAADO4XCQkqRVq1Zp1apVtpWps/3tb39zSmFwDelnvrn3XclxVZw4pfBgf3MLAgDAhTj8rb0//elPyszM1KpVq1RaWqqysrImD3iWtiGBSr4sWJK0uZB/vwAAnM3hFal58+bp7bff1r333tsS9cAFZSRFKv/YCeUVlGvAVdFmlwMAgMtweEWqrq5OvXv3bola4KIyzvRJcQNjAACacjhI3X///Vq0aFFL1AIXZZ1wvrmgXA0M5gQAwMbhS3s1NTV644039Omnn6pr167y92/afDxr1iynFQfXcFVMqIIDfHW8tl7flVSpY2yY2SUBAOASHA5S27ZtU/fu3SVJO3bsaPKaxWJxSlFwLX6+PuqWEKGv9h9TXn45QQoAgDMcDlKrV69uiTrg4nokR54OUgVluvvaJLPLAQDAJTjcI2W1b98+rVixQidPnpQkGQa9M54sIzlCkpSXT8M5AABWDgepY8eOaeDAgbryyis1ZMgQFRUVSTrdhP6HP/zB6QXCNaQnnv7m3v7SapVV15lcDQAArsHhIPX73/9e/v7+KigoUHBwsG37XXfdpeXLlzu1OLiOyNYBah/VWhKDOQEAsHI4SK1cuVIzZsxQQkJCk+2pqanKz893WmFwPdYbGG/i8h4AAJKaEaSqq6ubrERZlZaWKjAw0ClFwTVZg1Refrm5hQAA4CIcDlL9+vXTO++8Y3tusVjU2Niol156Sdddd51Ti4NrsTacbz1UrvqGxgvvDACAF3B4/MFLL72kAQMGaOPGjaqrq9PkyZO1c+dO/fjjj/ryyy9boka4iNToUIUG+qmqtl57fqhS5/hws0sCAMBUDq9IderUSdu2bdMvf/lLDRo0SNXV1br99tu1efNmdejQoSVqhIvw9bGo+5nbxTAGAQAAB1ekTp06pczMTP3lL3/Rn/70p5aqCS4sPSlSX3xXqryCct3by+xqAAAwl0MrUv7+/tqxYwe3gvFi1hsY5xWwIgUAgMOX9u677z699dZbLVEL3ED6mW/u5R87odLjtSZXAwCAuRxuNq+rq9Nf//pX5eTk6JprrlHr1q2bvD5r1iynFQfXEx7kr9ToEH1Xclx5+WXK7BxrdkkAAJjG4SC1Y8cOZWRkSJL27t3b5DUu+XmHjKTI00GqoJwgBQDwag4FqYaGBk2dOlVpaWlq06ZNS9UEF5eRHKHFGwvpkwIAeD2HeqR8fX114403qqKioqXqgRvokXy6T2rboXKdYjAnAMCLOdxsnpaWpv3797dELXAT7duGKKyVn2pONWpXUaXZ5QAAYBqHg9S0adP06KOP6uOPP1ZRUZEqKyubPOD5fHwstm/vMZgTAODNHG42v+mmmyRJv/rVr5o0lxuGIYvFooaGBudVB5eVkRSp3L1HlVdQrjF9zK4GAABzOBykVq9e3RJ1wM1Y+6RoOAcAeDOHg1T//v1bog64mW6J4bJYpENlJ1VSWaPosFZmlwQAwCXncJD6/PPPL/h6v379ml0M3EdoK39dFROq3cVVyiso001d4swuCQCAS87hIDVgwIBztp3dK0WPlPdIT4o8E6TKCVIAAK/k8Lf2ysrKmjxKSkq0fPly/eIXv9DKlStboka4KFufFN/cAwB4KYdXpMLDw8/ZNmjQIAUGBur3v/+9Nm3a5JTC4PoykiIkSdsOV6iuvlEBfg7ncgAA3JrT/vJFRUVpz549zjoc3EBK29aKDPZXXX2jdh5h2j0AwPs4vCK1bdu2Js8Nw1BRUZFefPFFdevWzWmFwfVZLKcHc362u0R5BeW2IZ0AAHgLh4NU9+7dZbFYZBhGk+09e/bU3/72N6cVBvfQI9kapMr0X0oxuxwAAC4ph4PUgQMHmjz38fFRVFSUWrVijpA3Sj/TJ0XDOQDAGzkcpJKTk1uiDripbgkR8rFIRRU1Kqo4qbjwILNLAgDgknG42fx3v/udXnvttXO2Z2dna9KkSc6oCW6kdaCfro4LkyTl5ZebWwwAAJeYw0Hqgw8+UJ8+596ltnfv3vr73//ulKLgXjKSuO8eAMA7ORykjh07dt5ZUmFhYSotLXVKUXAvGckRkqRN9EkBALyMw0Hqiiuu0PLly8/ZvmzZMrVv394pRcG9WFekdh6pUM0pbhEEAPAeDjebP/LII5o4caKOHj2q66+/XpK0atUqvfLKK5o9e7az64MbSGoTrLYhASo9XqedRyrUI7mN2SUBAHBJOBykxo4dq9raWk2bNk3PP/+8JKldu3aaO3eu7rvvPqcXCNdnHcyZ8+0P2pRfRpACAHgNh4OUJD344IN68MEHdfToUQUFBSkkJMTZdcHNZJwJUnxzDwDgTZoVpKyioqKcVQfcnPUGxnkFZTIMQxaLxdyCAAC4BJx202J4t64JEfLzsaikqlaHy0+aXQ4AAJcEQQpOERTgq07xpwdzMgYBAOAtCFJwGusYhM0F5eYWAgDAJXJRQaqmpsZZdcADpJ/VJwUAgDdwOEg1Njbq+eef1+WXX66QkBDt379fkvT000/rrbfecnqBcB89kk+vSH17pFIn6xjMCQDwfA4HqRdeeEFvv/22Zs6cqYCAANv2tLQ0/fWvf3VqcXAvl0cEKTo0UPWNhrYdKje7HAAAWpzDQeqdd97RG2+8oVGjRsnX19e2vWvXrtq9e7dTi4N7sVgsZ93AuNzcYgAAuAQcDlKHDx/WFVdccc72xsZGnTp1yilFwX1Zb2BMnxQAwBs4HKQ6d+6sL7744pzt/+///T+lp6c7pSi4L2uf1OYzgzkBAPBkDk82f/bZZ3Xvvffq8OHDamxs1Icffqg9e/bonXfe0ccff9wSNcKNdI4Pl7+vRaXH61Tw4wklX9ba7JIAAGgxDq9IDR06VIsXL9Y///lPWSwWPfPMM9q1a5f+8Y9/aNCgQS1RI9xIK39fdY4Pl8TlPQCA52vWvfZuvPFG3Xjjjc6uBR4iIylSWwrLlZdfrtvSE8wuBwCAFsNkczidtU+KFSkAgKeza0UqMjJSFovFrgP++OOPF1UQ3J/1m3u7iipVXVuv1oHNWvgEAMDl2fUXbvbs2S1cBjxJXHiQ4sJbqaiiRlsPlat3h7ZmlwQAQIuwK0iNHj26RU5eX1+vqVOn6r333lNxcbHi4uI0ZswYPfXUU/LxOX3V0TAM/elPf9Ibb7yhsrIyXXvttfrzn/+szp07244zbtw4ffrppzpy5IhCQkLUu3dvzZgxQx07dvzJc2dlZenDDz/U7t27FRQUZHvPVVddZdtnzJgxWrBgQZP3XXvttVq/fr2TPwnPk5EUqU+2F2lzAUEKAOC5mtUj1dDQoL///e96/vnn9cILL+iDDz5QfX29w8eZMWOG5s2bp+zsbO3atUszZ87USy+9pDlz5tj2mTlzpmbNmqXs7Gxt2LBBsbGxGjRokKqqqmz79OjRQ/Pnz9euXbu0YsUKGYahzMxMNTT89P3ecnNzNWHCBK1fv145OTmqr69XZmamqqurm+x30003qaioyPb45z//6fDv6Y0yrH1S+fRJAQA8l8VwcGrijh07NGzYMBUXF9tWb/bu3auoqCh99NFHSktLs/tYt9xyi2JiYprc7PiOO+5QcHCwFi5cKMMwFB8fr0mTJmnKlCmSpNraWsXExGjGjBkaN27ceY+7bds2devWTfv27VOHDh3squXo0aOKjo5Wbm6u+vXrJ+n0ilR5ebmWLl1q9+90tsrKSoWHh6uiokJhYWHNOoa72lxQptteX6fIYH/lPT3I7h47AADM5sjfb4dXpO6//3517txZhw4dUl5envLy8lRYWKiuXbvqt7/9rUPH6tu3r1atWqW9e/dKkrZu3aq1a9dqyJAhkqQDBw6ouLhYmZmZtvcEBgaqf//+Wrdu3XmPWV1drfnz5yslJUWJiYl211JRUSFJatOmTZPta9asUXR0tK688kr95je/UUlJyU8eo7a2VpWVlU0e3qpzfLgC/HxUduKUDpRW//wbAABwQw4Hqa1btyorK0uRkZG2bZGRkZo2bZq2bNni0LGmTJmikSNHqmPHjvL391d6eromTZqkkSNHSpKKi4slSTExMU3eFxMTY3vN6vXXX1dISIhCQkK0fPly5eTkKCAgwK46DMPQI488or59+6pLly627YMHD9Z7772nzz77TK+88oo2bNig66+/XrW1tec9TlZWlsLDw20PR4Kcpwnw81Ha5dbBnOXmFgMAQAtxOEhdddVV+uGHH87ZXlJSct6bGV/I4sWL9e6772rRokXKy8vTggUL9PLLL5/T4P3vl4UMwzhn26hRo7R582bl5uYqNTVVw4cPV01NjV11TJw4Udu2bdP//M//NNl+11136eabb1aXLl00dOhQLVu2THv37tUnn3xy3uM88cQTqqiosD0KCwvtOr+nYp4UAMDT2fWtvbMvUU2fPl2/+93vNHXqVPXs2VOStH79ej333HOaMWOGQyd/7LHH9Pjjj2vEiBGSpLS0NOXn5ysrK0ujR49WbGysJNm+0WdVUlJyziqVdRUoNTVVPXv2VGRkpJYsWWJb3fopDz30kD766CN9/vnnSki48BTuuLg4JScn67vvvjvv64GBgQoMDPzZ39tbZCRFSKLhHADguewKUhEREU1WgAzD0PDhw23brP3qQ4cOveA35f7diRMnbGMOrHx9fdXY2ChJSklJUWxsrHJycpSeni5JqqurU25u7s+GNsMwfvISnPX1hx56SEuWLNGaNWuUkpLys/UeO3ZMhYWFTUIdflpG0ukVqT0/VKmq5pRCW/mbXBEAAM5lV5BavXp1i5x86NChmjZtmpKSktS5c2dt3rxZs2bN0tixYyWdvqQ3adIkTZ8+XampqUpNTdX06dMVHBysu+++W5K0f/9+LV68WJmZmYqKitLhw4c1Y8YMBQUF2ZrWJWngwIG67bbbNHHiREnShAkTtGjRIv3f//2fQkNDbT1X4eHhCgoK0vHjxzV16lTdcccdiouL08GDB/Xkk0+qbdu2uu2221rk8/A00WGtlBAZpENlJ7W1sEJ9U5knBQDwLHYFqf79+7fIyefMmaOnn35a48ePV0lJieLj4zVu3Dg988wztn0mT56skydPavz48baBnCtXrlRoaKgkqVWrVvriiy80e/ZslZWVKSYmRv369dO6desUHR1tO87333+v0tJS2/O5c+dKkgYMGNCkpvnz52vMmDHy9fXV9u3b9c4776i8vFxxcXG67rrrtHjxYtu58fMykiJ1qOykNuWXEaQAAB7HrjlS27ZtU5cuXeTj46Nt27ZdcN+uXbs6rTh3581zpKze/vKApv7jW/W/MkoLxv7S7HIAAPhZjvz9tmtFqnv37iouLlZ0dLS6d+8ui8Wi8+Uvi8XiUI8UPJ91wvnmgjI1Nhry8WEwJwDAc9gVpA4cOKCoqCjbz4C9ro4LUyt/H1XW1Gt/6XFdEc1lUQCA57ArSCUnJ5/3Z+Dn+Pv6qGtChL458KM25ZcRpAAAHsXhgZwLFixoMpBy8uTJioiIUO/evZWfn+/U4uAZrGMQ8vLLzS0EAAAnczhITZ8+XUFBQZKkr776StnZ2Zo5c6batm2r3//+904vEO7PNpiTCecAAA9j16W9sxUWFtpuBbN06VL9+te/1m9/+1v16dPnnFECgPSvhvPvSo6r4uQphQcxmBMA4BkcXpEKCQnRsWPHJEkrV67UDTfcIOn0PKeTJ086tzp4hLYhgUq+LFjS6W/vAQDgKRwOUoMGDdL999+v+++/X3v37tXNN98sSdq5c6fatWvn7PrgIWx9UgXl5hYCAIATORyk/vznP6tXr146evSoPvjgA1122WWSpE2bNv3sDYLhvax9UqxIAQA8iV2TzdE8TDb/l51HKnTza2sVGuinLc9mypfBnAAAF+XI32+HV6SWL1+utWvX2p7/+c9/Vvfu3XX33XerrIzVBpzfVTGhCg7wVVVtvb4rqTK7HAAAnMLhIPXYY4+psrJSkrR9+3b94Q9/0JAhQ7R//3498sgjTi8QnsHP10fdEiIkMU8KAOA5HA5SBw4cUKdOnSRJH3zwgW655RZNnz5dr7/+upYtW+b0AuE5MpIjJDFPCgDgORwOUgEBATpx4oQk6dNPP1VmZqYkqU2bNraVKuB8eiRbv7lHkAIAeAaHB3L27dtXjzzyiPr06aNvvvlGixcvliTt3btXCQkJTi8QniM98XSQ2n+0WmXVdYpsHWByRQAAXByHV6Sys7Pl5+env//975o7d64uv/xySdKyZct00003Ob1AeI7I1gFq37a1JGlzIatSAAD35/CKVFJSkj7++ONztv/3f/+3UwqCZ0tPitT+0mrl5Zfr+o4xZpcDAMBFsStIVVZW2uYo/FwflLfPS8KF9UiO1Ad5h+iTAgB4BLuCVGRkpIqKihQdHa2IiAhZLOcOUzQMQxaLRQ0NDU4vEp7D+s29LYXlqm9olJ+vw1eXAQBwGXYFqc8++0xt2rSRJK1evbpFC4JnS40OVUign47X1mvPD1XqHB9udkkAADSbXUGqf//+5/0ZcJSvj0XdEyO0dl+p8grKCVIAALfmcLO5JNXU1Gjbtm0qKSlRY2Njk9d+9atfOaUweK6M5Eit3VeqzfllurdnstnlAADQbA4HqeXLl+u+++5TaWnpOa/RIwV7ZCRFSJI20XAOAHBzDnf6Tpw4UXfeeaeKiorU2NjY5EGIgj2sgznzj51Q6fFak6sBAKD5HA5SJSUleuSRRxQTwwwgNE94sL9So0MkSZsLys0tBgCAi+BwkPr1r3+tNWvWtEAp8CYZSadXpTblc3kPAOC+HO6Rys7O1p133qkvvvhCaWlp8vf3b/L67373O6cVB8+VkRyhxRsLGcwJAHBrDgepRYsWacWKFQoKCtKaNWuaDOe0WCwEKdjFuiK17VC5TjU0yp/BnAAAN+RwkHrqqaf03HPP6fHHH5ePD3/80DwdokIU1spPlTX12l1UpbQE5kkBANyPw0morq5Od911FyEKF8XHx6J0W5/UjyZXAwBA8zichkaPHq3Fixe3RC3wMtbLe3l8cw8A4KYcvrTX0NCgmTNnasWKFerates5zeazZs1yWnHwbNYbGNNwDgBwVw4Hqe3btys9PV2StGPHjiavnd14Dvyc7okRslikQ2UnVVJZo+iwVmaXBACAQxwKUg0NDZo6darS0tLUpk2blqoJXiK0lb+uignV7uIq5RWU6aYucWaXBACAQxzqkfL19dWNN96oioqKlqoHXiadPikAgBtzuNk8LS1N+/fvb4la4IWsNzDOY8I5AMANORykpk2bpkcffVQff/yxioqKVFlZ2eQBOKJH8pnBnIcrVFffaHI1AAA4xuFm85tuukmS9Ktf/apJc7lhGLJYLGpoaHBedfB4KW1bKzLYX2UnTmnnkQrbpT4AANyBw0Fq9erVLVEHvJTFcnow52e7S5RXUE6QAgC4FYeDVP/+/VuiDnixjKSIM0GqTP+lFLPLAQDAbs26z8sXX3yhe+65R71799bhw4clSQsXLtTatWudWhy8Q8aZPqnNNJwDANyMw0Hqgw8+0I033qigoCDl5eWptrZWklRVVaXp06c7vUB4vm4JEfKxSEcqalRUcdLscgAAsJvDQeqFF17QvHnz9Oabbza5PUzv3r2Vl5fn1OLgHVoH+qljbJgkKS+/3NxiAABwgMNBas+ePerXr98528PCwlReXu6MmuCFuO8eAMAdORyk4uLitG/fvnO2r127Vu3bt3dKUfA+1nlSBCkAgDtxOEiNGzdODz/8sL7++mtZLBYdOXJE7733nh599FGNHz++JWqEF8g4M/Zgx+EK1ZxiFhkAwD04PP5g8uTJqqio0HXXXaeamhr169dPgYGBevTRRzVx4sSWqBFeIKlNsC5rHaBj1XXaeaRCPZK5KTYAwPU1a/zBtGnTVFpaqm+++Ubr16/X0aNH9fzzzzu7NngR62BOiYZzAID7sDtInThxQhMmTNDll1+u6Oho3X///WrXrp1++ctfKiQkpCVrhJegTwoA4G7sDlLPPvus3n77bd18880aMWKEcnJy9OCDD7ZkbfAyGUkRkqRN+WUyDMPcYgAAsIPdPVIffvih3nrrLY0YMUKSdM8996hPnz5qaGiQr69vixUI79E1IUJ+PhaVVNXqcPlJJUQGm10SAAAXZPeKVGFhof7jP/7D9vyXv/yl/Pz8dOTIkRYpDN4nKMBXV8edGcxZUG5uMQAA2MHuINXQ0KCAgIAm2/z8/FRfX+/0ouC9bH1S3HcPAOAG7L60ZxiGxowZo8DAQNu2mpoaPfDAA2rdurVt24cffujcCuFV0pMi9PY6Gs4BAO7B7iA1evToc7bdc889Ti0GsA7m/PZIpWpONaiVP/13AADXZXeQmj9/fkvWAUiSEiKDFB0aqJKqWm07VKFfpjCYEwDgupo1kBNoKRaLxbYqtYk+KQCAiyNIweVkJEdIok8KAOD6CFJwOdYVqc0FDOYEALg2ghRcTpfLw+Xva1Hp8ToV/njS7HIAAPhJBCm4nFb+vuocHy5J2lTwo8nVAADw0whScEnWy3t5+eXmFgIAwAUQpOCSaDgHALgDghRckvVWMbuLq1Rdy22IAACuiSAFlxQXHqS48FZqaDS09VC52eUAAHBeBCm4rH+NQSg3txAAAH4CQQouKz0pQpKUx4RzAICLIkjBZVn7pPIYzAkAcFEEKbiszvHhCvDzUdmJUzpQWm12OQAAnIMgBZcV4OejtMtPD+bMo08KAOCCCFJwaRnWPinmSQEAXJCpQaq+vl5PPfWUUlJSFBQUpPbt2+u5555TY2OjbR/DMDR16lTFx8crKChIAwYM0M6dO5scZ9y4cerQoYOCgoIUFRWlYcOGaffu3Rc8d1ZWln7xi18oNDRU0dHRuvXWW7Vnz54m+9hzbrQsW58UDecAABdkapCaMWOG5s2bp+zsbO3atUszZ87USy+9pDlz5tj2mTlzpmbNmqXs7Gxt2LBBsbGxGjRokKqqqmz79OjRQ/Pnz9euXbu0YsUKGYahzMxMNTQ0/OS5c3NzNWHCBK1fv145OTmqr69XZmamqqv/1Ytjz7nRsqwjEPb8UKWqmlMmVwMAQFMWw8SvQ91yyy2KiYnRW2+9Zdt2xx13KDg4WAsXLpRhGIqPj9ekSZM0ZcoUSVJtba1iYmI0Y8YMjRs37rzH3bZtm7p166Z9+/apQ4cOdtVy9OhRRUdHKzc3V/369Wv2uc9WWVmp8PBwVVRUKCwszK46cK4+L36mw+Un9e5/Xau+qW3NLgcA4OEc+ftt6opU3759tWrVKu3du1eStHXrVq1du1ZDhgyRJB04cEDFxcXKzMy0vScwMFD9+/fXunXrznvM6upqzZ8/XykpKUpMTLS7loqKCklSmzZtmn3u2tpaVVZWNnng4mWcNQYBAABXYmqQmjJlikaOHKmOHTvK399f6enpmjRpkkaOHClJKi4uliTFxMQ0eV9MTIztNavXX39dISEhCgkJ0fLly5WTk6OAgAC76jAMQ4888oj69u2rLl26OHxuq6ysLIWHh9sejgQ5/LQeNJwDAFyUqUFq8eLFevfdd7Vo0SLl5eVpwYIFevnll7VgwYIm+1kslibPDcM4Z9uoUaO0efNm5ebmKjU1VcOHD1dNTY1ddUycOFHbtm3T//zP/5zzmj3ntnriiSdUUVFhexQWFtp1flxYxlkN542NDOYEALgOPzNP/thjj+nxxx/XiBEjJElpaWnKz89XVlaWRo8erdjYWEmnV4fi4uJs7yspKTlnpci6CpSamqqePXsqMjJSS5Yssa1u/ZSHHnpIH330kT7//HMlJCTYtjtybqvAwEAFBgY68AnAHlfHhamVv48qa+q1v/S4rogONbskAAAkmbwideLECfn4NC3B19fXNv4gJSVFsbGxysnJsb1eV1en3Nxc9e7d+4LHNgxDtbW1F3x94sSJ+vDDD/XZZ58pJSWlyesXc244l7+vj7peHiFJyssvN7UWAADOZmqQGjp0qKZNm6ZPPvlEBw8e1JIlSzRr1izddtttkk5fVps0aZKmT5+uJUuWaMeOHRozZoyCg4N19913S5L279+vrKwsbdq0SQUFBfrqq680fPhwBQUF2ZrWJWngwIHKzs62PZ8wYYLtsmJoaKiKi4tVXFyskydP2n1uXDo0nAMAXJGpl/bmzJmjp59+WuPHj1dJSYni4+M1btw4PfPMM7Z9Jk+erJMnT2r8+PEqKyvTtddeq5UrVyo09PTlnVatWumLL77Q7NmzVVZWppiYGPXr10/r1q1TdHS07Tjff/+9SktLbc/nzp0rSRowYECTmubPn68xY8bYdW5cOtYJ55sYzAkAcCGmzpHydMyRcp7S47W65oVPJUlbn81UeJC/yRUBADyV28yRAuzVNiRQSW2CJUlbCsvNLQYAgDMIUnAb1vvucXkPAOAqCFJwG9Y+qc00nAMAXARBCm4j/cwNjLcUlDOYEwDgEghScBsdY0MVHOCrqtp6fVdy3OxyAAAgSMF9+Pn6qFtChCT6pAAAroEgBbeSkRwhicGcAADXQJCCW8lIYsI5AMB1EKTgVqwN5/uPVqusus7kagAA3o4gBbfSpnWA2rdtLUnaXMiqFADAXAQpuB3rqlRefrm5hQAAvB5BCm6HhnMAgKsgSMHtWG8Vs7WwXPUNjSZXAwDwZgQpuJ3U6FCFBPqpuq5Be36oMrscAIAXI0jB7fj6WNQ9MUKSlFdQbmotAADvRpCCW7LdwJgJ5wAAExGk4JYykhnMCQAwH0EKbik98XSQOnjshEqP15pcDQDAWxGk4JbCg/11RXSIJGkzfVIAAJMQpOC2rH1SXN4DAJiFIAW3ZZ0nlUfDOQDAJAQpuK2MM7eK2XqoXKcYzAkAMAFBCm6rQ1SIwlr5qeZUo3YXMZgTAHDpEaTgtnx8LOqexBgEAIB5CFJwaz0IUgAAExGk4NYykiMkSZtoOAcAmIAgBbfWPTFCFot0qOykSqpqzC4HAOBlCFJwa6Gt/HVldKgkKS+/3NxiAABehyAFt2e9795m+qQAAJcYQQpuzzrhnD4pAMClRpCC27OuSG07XKG6egZzAgAuHYIU3F77tq0VEeyvuvpGfVtUaXY5AAAvQpCC27NYLLbbxXB5DwBwKRGk4BGsfVIM5gQAXEoEKXgE64rUZlakAACXEEEKHqFbYoR8LNKRihoVVZw0uxwAgJcgSMEjtA70U8fYMEkM5gQAXDoEKXgM63336JMCAFwqBCl4DGufFEEKAHCpEKTgMXqcGcy583Clak41mFwNAMAbEKTgMZLaBOuy1gGqa2jUziMVZpcDAPACBCl4DIvFonTr5T0azgEAlwBBCh6FhnMAwKVEkIJH6XFWw7lhGCZXAwDwdAQpeJSuCRHy87Hoh8paHS5nMCcAoGURpOBRggJ8dXXcmcGcBeXmFgMA8HgEKXgc2w2Mue8eAKCFEaTgcTLOzJPaTMM5AKCFEaTgcawTznceYTAnAKBlEaTgcRIigxQVGqj6RkPbDjGYEwDQcghS8DgWi+VffVJc3gMAtCCCFDyS9b57NJwDAFoSQQoeKYPBnACAS4AgBY/U5fJw+ftaVHq8ToU/MpgTANAyCFLwSK38fdUpPlwSfVIAgJZDkILHOvu+ewAAtASCFDxWRnKEJGkTDecAgBbiZ3YBQEuxNpzvLq7Sibp6BQd4z3/ujY2G6hoaVXuqUbX1DaqtP/3PmlONtp9r6//99UbVnmqQM3vzLRbnHcuZLE4szEV/RcBrdE0I1zXt2ph2fu/5ywKvEx8RpNiwViqurNHWwgr16nDZJTt3fUPjWUHmp4JLw0++XmcNNj+z39kB6Oz9TzXwTUUA3mH8gA4EKaCl9EiO1Cfbi7R6T4kiW/s3DSTnhJPT/6yr/5nQYkewaWh0nSDjYzndfB/o56NAP18F+vv862c/HwX6+yjA9/TzAD8f+fk4Z43FmZ+AM0dYOLcuJx4MQLNcHRdm6vkJUvBo6UkR+mR7kd74fL/e+Hy/KTX4+1r+FVr8fBTo73s6uJwn0Pz7fraffzIAnX/72cfx86UVEgBaCkEKHm1wWpwWfHVQ5dWnmgSMADuDSpNgY9d+TV8P8PORr5NWeAAArocgBY92eUSQvph8vdllAAA8FGv+AAAAzUSQAgAAaCaCFAAAQDMRpAAAAJqJIAUAANBMBCkAAIBmIkgBAAA0E0EKAACgmUwNUvX19XrqqaeUkpKioKAgtW/fXs8995waGxtt+xiGoalTpyo+Pl5BQUEaMGCAdu7c2eQ448aNU4cOHRQUFKSoqCgNGzZMu3fvvuC5P//8cw0dOlTx8fGyWCxaunTpOfuMGTNGFoulyaNnz55O+d0BAID7MzVIzZgxQ/PmzVN2drZ27dqlmTNn6qWXXtKcOXNs+8ycOVOzZs1Sdna2NmzYoNjYWA0aNEhVVVW2fXr06KH58+dr165dWrFihQzDUGZmphoaGn7y3NXV1erWrZuys7MvWONNN92koqIi2+Of//znxf/iAADAI1gMZ95W3UG33HKLYmJi9NZbb9m23XHHHQoODtbChQtlGIbi4+M1adIkTZkyRZJUW1urmJgYzZgxQ+PGjTvvcbdt26Zu3bpp37596tChw8/WYbFYtGTJEt16661Nto8ZM0bl5eXnXa06n9raWtXW1tqeV1ZWKjExURUVFQoLM/fu1AAAwD6VlZUKDw+36++3qStSffv21apVq7R3715J0tatW7V27VoNGTJEknTgwAEVFxcrMzPT9p7AwED1799f69atO+8xq6urNX/+fKWkpCgxMfGia1yzZo2io6N15ZVX6je/+Y1KSkp+ct+srCyFh4fbHs44PwAAcF2mBqkpU6Zo5MiR6tixo/z9/ZWenq5JkyZp5MiRkqTi4mJJUkxMTJP3xcTE2F6zev311xUSEqKQkBAtX75cOTk5CggIuKj6Bg8erPfee0+fffaZXnnlFW3YsEHXX399k1Wnsz3xxBOqqKiwPQoLCy/q/AAAwLX5mXnyxYsX691339WiRYvUuXNnbdmyRZMmTVJ8fLxGjx5t289isTR5n2EY52wbNWqUBg0apKKiIr388ssaPny4vvzyS7Vq1arZ9d111122n7t06aJrrrlGycnJ+uSTT3T77befs39gYKACAwObfT4AAOBeTA1Sjz32mB5//HGNGDFCkpSWlqb8/HxlZWVp9OjRio2NlXR6ZSouLs72vpKSknNWqayX01JTU9WzZ09FRkZqyZIlttUtZ4iLi1NycrK+++47u/a3tp9VVlY6rQYAANCyrH+37WkjNzVInThxQj4+Ta8u+vr62sYfpKSkKDY2Vjk5OUpPT5ck1dXVKTc3VzNmzLjgsQ3D+MlLcM117NgxFRYWNgl1F2L9ZiG9UgAAuJ+qqiqFh4dfcB9Tg9TQoUM1bdo0JSUlqXPnztq8ebNmzZqlsWPHSjp9SW/SpEmaPn26UlNTlZqaqunTpys4OFh33323JGn//v1avHixMjMzFRUVpcOHD2vGjBkKCgqyNa1L0sCBA3Xbbbdp4sSJkqTjx49r3759ttcPHDigLVu2qE2bNkpKStLx48c1depU3XHHHYqLi9PBgwf15JNPqm3btrrtttvs+v3i4+NVWFio0NDQcy5FXizrNwILCwv5RuDP4LOyH5+V/fis7MdnZT8+K/u15GdlGIaqqqoUHx9v186mqaysNB5++GEjKSnJaNWqldG+fXvjj3/8o1FbW2vbp7Gx0Xj22WeN2NhYIzAw0OjXr5+xfft22+uHDx82Bg8ebERHRxv+/v5GQkKCcffddxu7d+9ucq7k5GTj2WeftT1fvXq1Iemcx+jRow3DMIwTJ04YmZmZRlRUlOHv728kJSUZo0ePNgoKClr0M7FXRUWFIcmoqKgwuxSXx2dlPz4r+/FZ2Y/Pyn58VvZzlc/K1DlSaD5HZlx4Oz4r+/FZ2Y/Pyn58Vvbjs7Kfq3xW3GsPAACgmQhSbiowMFDPPvss4xbswGdlPz4r+/FZ2Y/Pyn58VvZzlc+KS3sAAADNxIoUAABAMxGkAAAAmokgBQAA0EwEKQAAgGYiSLmZzz//XEOHDlV8fLwsFouWLl1qdkkuKSsrS7/4xS8UGhqq6Oho3XrrrdqzZ4/ZZbmsuXPnqmvXrgoLC1NYWJh69eqlZcuWmV2Wy8vKyrLdgQHnmjp1qiwWS5OH9R6qONfhw4d1zz336LLLLlNwcLC6d++uTZs2mV2Wy2nXrt05/11ZLBZNmDDBlHoIUm6murpa3bp1U3Z2ttmluLTc3FxNmDBB69evV05Ojurr65WZmanq6mqzS3NJCQkJevHFF7Vx40Zt3LhR119/vYYNG6adO3eaXZrL2rBhg9544w117drV7FJcWufOnVVUVGR7bN++3eySXFJZWZn69Okjf39/LVu2TN9++61eeeUVRUREmF2ay9mwYUOT/6ZycnIkSXfeeacp9Zh6rz04bvDgwRo8eLDZZbi85cuXN3k+f/58RUdHa9OmTerXr59JVbmuoUOHNnk+bdo0zZ07V+vXr1fnzp1Nqsp1HT9+XKNGjdKbb76pF154wexyXJqfnx+rUHaYMWOGEhMTNX/+fNu2du3amVeQC4uKimry/MUXX1SHDh3Uv39/U+phRQpeoaKiQpLUpk0bkytxfQ0NDXr//fdVXV2tXr16mV2OS5owYYJuvvlm3XDDDWaX4vK+++47xcfHKyUlRSNGjND+/fvNLsklffTRR7rmmmt05513Kjo6Wunp6XrzzTfNLsvl1dXV6d1339XYsWNlsVhMqYEgBY9nGIYeeeQR9e3bV126dDG7HJe1fft2hYSEKDAwUA888ICWLFmiTp06mV2Wy3n//feVl5enrKwss0txeddee63eeecdrVixQm+++aaKi4vVu3dvHTt2zOzSXM7+/fs1d+5cpaamasWKFXrggQf0u9/9Tu+8847Zpbm0pUuXqry8XGPGjDGtBi7tweNNnDhR27Zt09q1a80uxaVdddVV2rJli8rLy/XBBx9o9OjRys3NJUydpbCwUA8//LBWrlypVq1amV2Oyzu7DSEtLU29evVShw4dtGDBAj3yyCMmVuZ6Ghsbdc0112j69OmSpPT0dO3cuVNz587VfffdZ3J1ruutt97S4MGDFR8fb1oNrEjBoz300EP66KOPtHr1aiUkJJhdjksLCAjQFVdcoWuuuUZZWVnq1q2bXn31VbPLcimbNm1SSUmJevToIT8/P/n5+Sk3N1evvfaa/Pz81NDQYHaJLq1169ZKS0vTd999Z3YpLicuLu6c/9Ny9dVXq6CgwKSKXF9+fr4+/fRT3X///abWwYoUPJJhGHrooYe0ZMkSrVmzRikpKWaX5HYMw1Btba3ZZbiUgQMHnvOts//8z/9Ux44dNWXKFPn6+ppUmXuora3Vrl279B//8R9ml+Jy+vTpc86Ilr179yo5Odmkilyf9UtEN998s6l1EKTczPHjx7Vv3z7b8wMHDmjLli1q06aNkpKSTKzMtUyYMEGLFi3S//3f/yk0NFTFxcWSpPDwcAUFBZlcnet58sknNXjwYCUmJqqqqkrvv/++1qxZc863H71daGjoOX12rVu31mWXXUb/3Xk8+uijGjp0qJKSklRSUqIXXnhBlZWVGj16tNmluZzf//736t27t6ZPn67hw4frm2++0RtvvKE33njD7NJcUmNjo+bPn6/Ro0fLz8/kKGPAraxevdqQdM5j9OjRZpfmUs73GUky5s+fb3ZpLmns2LFGcnKyERAQYERFRRkDBw40Vq5caXZZbqF///7Gww8/bHYZLumuu+4y4uLiDH9/fyM+Pt64/fbbjZ07d5pdlsv6xz/+YXTp0sUIDAw0OnbsaLzxxhtml+SyVqxYYUgy9uzZY3YphsUwDMOcCAcAAODeaDYHAABoJoIUAABAMxGkAAAAmokgBQAA0EwEKQAAgGYiSAEAADQTQQoAAKCZCFIAAADNRJACACdq166dZs+ebXYZAC4RghQAtzVmzBjdeuutkqQBAwZo0qRJl+zcb7/9tiIiIs7ZvmHDBv32t7+9ZHUAMBc3LQaAs9TV1SkgIKDZ74+KinJiNQBcHStSANzemDFjlJubq1dffVUWi0UWi0UHDx6UJH377bcaMmSIQkJCFBMTo3vvvVelpaW29w4YMEATJ07UI488orZt22rQoEGSpFmzZiktLU2tW7dWYmKixo8fr+PHj0uS1qxZo//8z/9URUWF7XxTp06VdO6lvYKCAg0bNkwhISEKCwvT8OHD9cMPP9henzp1qrp3766FCxeqXbt2Cg8P14gRI1RVVWXb5+9//7vS0tIUFBSkyy67TDfccIOqq6tb6NME4AiCFAC39+qrr6pXr176zW9+o6KiIhUVFSkxMVFFRUXq37+/unfvro0bN2r58uX64YcfNHz48CbvX7Bggfz8/PTll1/qL3/5iyTJx8dHr732mnbs2KEFCxbos88+0+TJkyVJvXv31uzZsxUWFmY736OPPnpOXYZh6NZbb9WPP/6o3Nxc5eTk6Pvvv9ddd93VZL/vv/9eS5cu1ccff6yPP/5Yubm5evHFFyVJRUVFGjlypMaOHatdu3ZpzZo1uv3228X95gHXwKU9AG4vPDxcAQEBCg4OVmxsrG373LlzlZGRoenTp9u2/e1vf1NiYqL27t2rK6+8UpJ0xRVXaObMmU2OeXa/VUpKip5//nk9+OCDev311xUQEKDw8HBZLJYm5/t3n376qbZt26YDBw4oMTFRkrRw4UJ17txZGzZs0C9+8QtJUmNjo95++22FhoZKku69916tWrVK06ZNU1FRkerr63X77bcrOTlZkpSWlnYRnxYAZ2JFCoDH2rRpk1avXq2QkBDbo2PHjpJOrwJZXXPNNee8d/Xq1Ro0aJAuv/xyhYaG6r777tOxY8ccuqS2a9cuJSYm2kKUJHXq1EkRERHatWuXbVu7du1sIUqS4uLiVFJSIknq1q2bBg4cqLS0NN1555168803VVZWZv+HAKBFEaQAeKzGxkYNHTpUW7ZsafL47rvv1K9fP9t+rVu3bvK+/Px8DRkyRF26dNEHH3ygTZs26c9//rMk6dSpU3af3zAMWSyWn93u7+/f5HWLxaLGxkZJkq+vr3JycrRs2TJ16tRJc+bM0VVXXaUDBw7YXQeAlkOQAuARAgIC1NDQ0GRbRkaGdu7cqXbt2umKK65o8vj38HS2jRs3qr6+Xq+88op69uypK6+8UkeOHPnZ8/27Tp06qaCgQIWFhbZt3377rSoqKnT11Vfb/btZLBb16dNHf/rTn7R582YFBARoyZIldr8fQMshSAHwCO3atdPXX3+tgwcPqrS0VI2NjZowYYJ+/PFHjRw5Ut98843279+vlStXauzYsRcMQR06dFB9fb3mzJmj/fv3a+HChZo3b9455zt+/LhWrVql0tJSnThx4pzj3HDDDeratatGjRqlvLw8ffPNN7rvvvvUv3//815OPJ+vv/5a06dP18aNG1VQUKAPP/xQR48edSiIAWg5BCkAHuHRRx+Vr6+vOnXqpKioKBUUFCg+Pl5ffvmlGhoadOONN6pLly56+OGHFR4eLh+fn/6fv+7du2vWrFmaMWOGunTpovfee09ZWVlN9undu7ceeOAB3XXXXYqKijqnWV06vZK0dOlSRUZGql+/frrhhhvUvn17LV682O7fKywsTJ9//rmGDBmiK6+8Uk899ZReeeUVDR482P4PB0CLsRh8hxYAAKBZWJECAABoJoIUAABAMxGkAAAAmokgBQAA0EwEKQAAgGYiSAEAADQTQQoAAKCZCFIAAADNRJACAABoJoIUAABAMxGkAAAAmun/A4Ww4Y6llvQ4AAAAAElFTkSuQmCC",
      "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()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "b5519496-3d0e-462c-aa1b-7c3ac8e4d9bf",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The permissible current rating of one conductor is equal to   800 A\n"
     ]
    }
   ],
   "source": [
    "# Permissible current rounding\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.10.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}