sim.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Stabilizer current sense AFE noise/stability model\n",
    "\n",
    "All values are referred to the coil current (RTC)\n",
    "\n",
    "See *Merkel et al* in prep"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import constants\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def johnson(res, temp=300):\n",
    "    \"\"\" Johnson Voltage noise density caluclator\n",
    "    :param res: resistance (Ohm)\n",
    "    :param temp: temperature (K)\n",
    "    :return: RMS noise voltage density (V/rtHz)\n",
    "    \"\"\"\n",
    "    return np.sqrt(4*constants.k*res*temp) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# simulation parameters\n",
    "I_coil = 137  # coil current (A)\n",
    "field = 28.8e-3  # Magnetic field (Tesla)\n",
    "loop_bw = 1e3  # feedback loop noise bandwidth (Hz)\n",
    "\n",
    "G_coil = field/I_coil  # Coil gain (T/A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "power dissipated in sense resistor: 0.19W\n"
     ]
    }
   ],
   "source": [
    "# pick one of ..\n",
    "\n",
    "# Sensor option I: LEM-IT-400-S Sensor\n",
    "# data sheet: https://www.lem.com/sites/default/files/products_datasheets/it_400-s_ultrastab.pdf\n",
    "# NB a newer model is now available: https://www.lem.com/sites/default/files/products_datasheets/in_1000-s.pdf\n",
    "\n",
    "n_turns = 2  # number of turns of wire around the sensor\n",
    "sensor_Ifs = 400  # sensor full-scale current (A)\n",
    "G_sensor = 1/2000 * n_turns  # A/A\n",
    "R_sense = 10  # transimpedance resistor in analog front-end (Ohm)\n",
    "\n",
    "if (I_coil * n_turns) > sensor_Ifs:\n",
    "    raise ValueError(\"Too much flux through the sensor\")\n",
    "\n",
    "# Broad-band noise:\n",
    "# - We measured the sensor noise floor to be <~5ppb/rtHz (RMS of full-scale)\n",
    "#   this measurement was limited by the noise floor of our analog front end\n",
    "# - The sensor data sheet specifies 8ppm RMS MAX integrated over 50kHz\n",
    "#   eqivalent to 35ppb/rtHz\n",
    "# - so the specification appears to be rather cautious\n",
    "# - The newer IN-1000-S model specifies both typical and max noise,\n",
    "#   with a factor of three difference between the two\n",
    "in_sensor = 5e-9*(sensor_Ifs/n_turns) # sensor current noise (A RTC)\n",
    "\n",
    "# sensor low frequency noise:\n",
    "# 0 to 10Hz integrated noise (pk-pk normalized to full-scale coil current)\n",
    "# taken from datasheet, probably quite pessimistic\n",
    "In_sensor = 0.05e-6*(sensor_Ifs/n_turns)/I_coil\n",
    "\n",
    "print(\"power dissipated in sense resistor: {:.2f}W\".format((I_coil*G_sensor)**2*R_sense))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# or...\n",
    "\n",
    "# Sensor option II: resistor\n",
    "G_sensor = 1\n",
    "R_sense = 1/2000*6*10  # chosen to give same output voltage as the flux-gate\n",
    "in_sensor = 0  # no broard-band noise on top of the sense resistor Johnson noise\n",
    "In_sensor = 0  # no low-frequency noise\n",
    "\n",
    "print(\"power dissipated in {:.1f} mOhm sense resistor: {:.1f}W\".format(R_sense*1e3, I_coil**2*R_sense))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Broadband noise\n",
    "\n",
    "- Our measurements indicate that the flux-gate sensor noise is far lower than the MAX specification on the data sheet. We weren't able to measure it accurately in the previous tests we performed. For now we assume that it can be neglected.\n",
    "- With that assumption, the broadband noise is dominated by the 1k resistors we use in the InAmp and subtractor.\n",
    "- So, to lower the circuit noise we would need to lower the value of these resistors. Johnson noise scales like the root of the resistance, so we'd really need to lower it by a factor of 5-10 to make this worthwhile. Note, however, that when using resistor ladders we don't have a huge amount of choice in the resistor values.\n",
    "- If we lower the resistance to 100R then, assuming a 2V sensor voltage, the current through the resistors is 20mA. That's within the drive capability of many OpAmps, but still quite a bit.\n",
    "- NB 100R puts the noise contribution from the resistors about equal to using OpAmps with 1nV/rtHz (which is about as low as one can get). So there isn't much point pushing it lower than this.\n",
    "- The power dissipated in the resistor would be 40mW (160mW for a quad resistor package) which is OK.\n",
    "- Potentially worse would be the power dissipated in the OpAmp. Assuming we keep the +-12V rails, this would be 240mW per OpAmp. Depending on the OpAmp and package layout, that could gi~ve quite a large temperature increase and impact stability.\n",
    "- So, it looks like using 100R resistors and an OpAmp with 1nV/rtHz noise we can probably get the AFE broadband noise down by a factor of 2-3 at the cost of a significant increase in power dissipation. However, in practice, we may well not see that improvement since, for example, the sensor noise could start to dominate. Note also, that the low-frequency noise is harder to reduce without doing something daft like parallelling large numbers of voltage references.\n",
    "- My current feeling is that it's not worth trying to squeeze out that last iota of performance, but the option is there if a user has an application where they do need it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def print_white_noise(name, i_n):\n",
    "    \"\"\":param i_n: noise density (A/rtHz)\"\"\"\n",
    "    print(\"{}: {:.2f} uA/rtHz ({:.2f} nT/rtHz), {:.2f} uA RMS ({:.2f} nT RMS)\".format(\n",
    "        name,\n",
    "        i_n*1e6,\n",
    "        i_n*G_coil*1e9,\n",
    "        i_n*np.sqrt(loop_bw)*1e6,\n",
    "        i_n*np.sqrt(loop_bw)*G_coil*1e9))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "sense resistor: 0.04 uA/rtHz (0.01 nT/rtHz), 1.29 uA RMS (0.27 nT RMS)\n"
     ]
    }
   ],
   "source": [
    "# transimpedance stage\n",
    "G_sense = G_sensor*R_sense  # V/A RTC\n",
    "\n",
    "vn_sense = johnson(R_sense)\n",
    "in_sense = vn_sense/G_sense  # A/rtHz RTC\n",
    "\n",
    "print_white_noise(\"sense resistor\", in_sense)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "InAmp resistors: 0.81 uA/rtHz (0.17 nT/rtHz), 25.74 uA RMS (5.41 nT RMS)\n",
      "InAmp OpAmp vn: 0.61 uA/rtHz (0.13 nT/rtHz), 19.36 uA RMS (4.07 nT RMS)\n",
      "InAmp OpAmp in: 0.06 uA/rtHz (0.01 nT/rtHz), 1.79 uA RMS (0.38 nT RMS)\n",
      "InAmp total: 1.02 uA/rtHz (0.22 nT/rtHz), 32.35 uA RMS (6.80 nT RMS)\n"
     ]
    }
   ],
   "source": [
    "# InAmp\n",
    "# NB stage bandwidth is large compared with loop BW so assume white noise\n",
    "\n",
    "# ADA4004 (used in previous version of the design)\n",
    "# opamp_vn = 1.8e-9  # voltage noise density (V/rtHz)\n",
    "# opamp_in = 1.2e-12  # current noise density (A/rtHz)\n",
    "\n",
    "# LT1007\n",
    "opamp_vn = 2.5e-9  # voltage noise density (V/rtHz)\n",
    "opamp_in = 400e-15  # current noise density (A/rtHz)\n",
    "\n",
    "R_in = 16  # RFI filter resistors\n",
    "R_inamp = 1e3  # InAmp resistors\n",
    "\n",
    "noise_sources = [\n",
    "    johnson(R_in)*np.sqrt(2),\n",
    "    johnson(R_inamp)*np.sqrt(4),\n",
    "    opamp_vn*np.sqrt(1**2 + 1**2 + 2**2),  # noise gain of final OpAmp is 2\n",
    "    opamp_in*R_inamp*np.sqrt(2)\n",
    "]\n",
    "\n",
    "vn_inamp = np.sqrt(sum([noise**2 for noise in noise_sources]))\n",
    "in_inamp = vn_inamp/G_sense\n",
    "\n",
    "print_white_noise(\"InAmp resistors\", noise_sources[1]/G_sense)\n",
    "print_white_noise(\"InAmp OpAmp vn\", noise_sources[2]/G_sense)\n",
    "print_white_noise(\"InAmp OpAmp in\", noise_sources[3]/G_sense)\n",
    "print_white_noise(\"InAmp total\", in_inamp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Subtractor resistors: 0.41 uA/rtHz (0.09 nT/rtHz), 12.81 uA RMS (2.69 nT RMS)\n",
      "Subtractor OpAmp vn: 0.25 uA/rtHz (0.05 nT/rtHz), 7.91 uA RMS (1.66 nT RMS)\n",
      "Subtractor OpAmp in: 0.04 uA/rtHz (0.01 nT/rtHz), 1.25 uA RMS (0.26 nT RMS)\n",
      "Subtractor reference: 0.25 uA/rtHz (0.05 nT/rtHz), 7.83 uA RMS (1.65 nT RMS)\n",
      "Subtractor total: 0.54 uA/rtHz (0.11 nT/rtHz), 17.01 uA RMS (3.58 nT RMS)\n"
     ]
    }
   ],
   "source": [
    "# error subtractor\n",
    "R_fb = 100e3\n",
    "R_g = 1e3\n",
    "G_subtractor = 1 + R_fb/R_g\n",
    "\n",
    "# all numbers referred to subtractor input\n",
    "vn_r_g = johnson(R_g)*(R_fb/R_g)/G_subtractor\n",
    "vn_r_fb = johnson(R_fb)/G_subtractor\n",
    "vn_res = np.sqrt(vn_r_g**2+vn_r_fb**2)\n",
    "\n",
    "vn_opamp_in = opamp_in*R_fb/G_subtractor\n",
    "vn_opamp = np.sqrt(opamp_vn**2 + vn_opamp_in**2)\n",
    "\n",
    "vn_ref = opamp_vn*(R_fb/R_g)/G_subtractor  # only OpAmp vn as all else filtered!\n",
    "\n",
    "vn_sub = np.sqrt(vn_res**2 + vn_opamp**2 + vn_ref**2)\n",
    "\n",
    "in_sub = vn_sub/G_sense\n",
    "\n",
    "print_white_noise(\"Subtractor resistors\", vn_res/G_sense)\n",
    "print_white_noise(\"Subtractor OpAmp vn\", opamp_vn/G_sense)\n",
    "print_white_noise(\"Subtractor OpAmp in\", vn_opamp_in/G_sense)\n",
    "print_white_noise(\"Subtractor reference\", vn_ref/G_sense)\n",
    "print_white_noise(\"Subtractor total\", in_sub)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total: 1.53 uA/rtHz (0.32 nT/rtHz), 48.33 uA RMS (10.16 nT RMS)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWQAAADuCAYAAAAOR30qAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl8XGW9x/HPM0uSmexb9yVdhkLpdIFCKFCosiOURQREvQVXuMDVq6IVvFLxAvcqiFdAQUCMQFlEhECVTRZZUwpdQoE0pfvepGkyk9nOmfPcP850obS0TSY5s/zer1dIOplz5pch880zz3kWpbVGCCGE81xOFyCEEMImgSyEEBlCAlkIITKEBLIQQmQICWQhhMgQEshCCJEhJJCFECJDSCALIUSGkEAWQogMIYEshBAZQgJZCCEyhASyEEJkCAlkIYTIEBLIQgiRISSQhRAiQ0ggCyFEhpBAFkKIDCGBLIQQGcLjdAFCHKi62fNcQBVQA5QCXsANuEoPm62AJGCmPrqBrUB786xmy5mKhTg4SvbUE5mgbvY8N1AHHAqMAw4BBmOHbw1QDVSy13d1llV62LX7erdnAe3Y4bwl9XkrsBb4AFgKrJTQFplAAln0u7rZ86qB44CjgMOwQ3gsUNiT82md1GXjr1O9KCkKfIQdzkuxg/rd5lnN63txTiEOmgSy6HN1s+eNBo7f8aG1PlQp1ZsA/QStTV02/qdpO99ulgMv7/hontW8qQ8eQ4idJJBF2tXNnucDTgPO01qfrJQa0pePl4YW8oH6CHgFO6Cfb57VvL0fHlPkEQlkkRZ1s+dVAGdrrc8DTlNK+fvrsfsxkHeXAJ4DHgEam2c1h/v58UUOkkAWPVY3e14xcLHW+iJghlLK60QdDgXy7qLAPOBRYF7zrOaog7WILCaBLA5a3ex5E7W2rgC+qpSrxOl6MiCQdxcGngTuaJ7V3OR0MSK7SCCLA5LqF75IW8krlcs91el6dpdhgby7JuA3wOPNs5pNp4sRmU8CWXymutnzBmjLugb4tnK5ypyuZ28yOJB3WAfcCfyheVbztv54QKXUdcAl2JNlLOA7Wuu9ttiVUnOAsNb6lj1urwOO1VrPTVNN5wLLtNYfpON8uUimTou9qps9b9CIHzxxp7asNcrl+mGmhnGWGAbcDKwNNgTvCjYER/XlgymlpgFnAUdorScCJ2NPhDlYddihvrfH6Mks33OB8QdzQA8fJ2tJC1l8Qt3seUMsI369cnsvVS5XgdP1HIgsaCHvKQHcDfyieVbz1nSfXCl1PnCZ1vrsPW5fBUzVWrcppaYCt2itZ6RayGOAocBw4Jda63uUUm9jT9xZCTQAHcAXgCKgGJgJPIU9g9IL/FRr/VTqsf4N+CGggSXA74FngM7Uxxexp7/fBfiBj4Gva607lFKvAG9iTx5q1Frfmu7nKFPl1V8fsW91s+cNsIz4jcrtneXyFjoyWiKPFABXA5cGG4K3ALc2z2ruTuP5nwd+ppRaBrwIPKq1fnU/x0wEjsEO2oVKqXnAbOCHWuuzAJRSlwLTgIla622p1ut5WusupVQN8LZSqhG7FXwdcFwq/KtS928EntFaP5463xLgaq31q0qpG4Drge+l6qnQWp+YrickW0gg57m62fM8yVjox64C/7Uub2G/jR0WgN1C/Dnw78GG4A3APc2zmo3enlRrHVZKHQlMBz4HPKqUmr2fw57SWkeBqFLqZeBoYG8TX17QWu/oB1fATUqpE7D7qYcCA4HPA49rrdtS9Xyq31wpVY4dujv+UDQAf9ntLo8ewI+ac6QPOY8Nv/qh061E9GN3Uel/K5dbwtg5A7Ev+n0QbAielY4Taq2TWutXtNbXA1dhdxGY7HrNF+15yH7+vcPuLfmvALXAkVrrycDm1HnVZxx/oNL5jiFrSCDnoRHfe3T48O8+/Ly7uOIfrgLfCKfrETuNBZ4ONgTnBhuCNT09iVJqnFIqsNtNk4HVwCrgyNRtX9zjsHOUUkVKqWpgBvAOEMJuxe9LObBFa20opT4HjEzd/k/gwtS5UEpVpW7feT6tdSfQoZSanvre14D9davkPOmyyCN1s+e5kuGOG1y+0mtcbk9WXLDLU18GTgk2BP+jeVbzwz04vgS4XSlVgd0qXg58G/sC3X1KqWuxx0jvbj72bMMRwC+01huUUlsBUym1GPgT9kW93T0EPK2UWgAswl7rA631UqXUjcCrSqkksBC4FHua+T1Kqf8ALgBmAXelptmvAC7rwc+aU2SURZ4Y8vXbJ7pLa//i9pUe4nQt6ZaFoywOxtPAFbIUaH6QLosc5w/UqyHf/P2N3uoRC3IxjPPA2dh9y992uhDR96SFnMMGXnjDCG/NyCc9ZTVTnK6lL+V4C3l3TwGXyrKfuUtayDlq0Ndu+VrhsPFLcz2M88w5wMJgQ/AopwsRfUMCOcf4A/XewbNu+2PhkHENrgKf4yuxibSrA14PNgSvdLoQkX7SZZFDyo46d1DJxFOfLqgdmVGrsfW1POqy2NP92Bf84k4XItJDWsg5ourk70wtO+q8d/MtjPPcZcBrwYbgMKcLEekhgZzl/IF6VXPWD75aMuGklzxlNX26d53ISEcBbwcbgoc7XYjoPQnkLOYP1LtLJp5yU/Fh0//oKir+rBlVIrcNBf4VbAge43QhonckkLOUP1BfXBw8+T7fmPofK7dXVmcTVcCLwYbgaU4XInpOAjkL+QP1VcUTTnrYH5j2b8rlyseLWWLvirHXwviy04WInpFAzjL+QH1NycRTH/YfcuzZSikJY7EnL/BQsCF4ldOFiIMngZxF/IH6gSWTz/iLP3DMqZLF4jMo4PZgQ/BapwsRB0cCOUv4A/VDSo84+6/+MUfNcLoWkTVulAkk2UUCOQv4A/XDSyae9rBv1JTjnK5FZJ3bpU85e0ggZzh/oH6EL3DMPb6xR0/f/72F+BQFNAQbgqc7XYjYPwnkDOYP1NcWDp/w65IJJ50sF/BEL3iBvwYbgsc6XYj4bBLIGcofqC/zDhh9Y9kRZ52lXG630/WIrOcHngk2BINOFyL2TQI5A/kD9UWe8kE/KT/6/C8rT0Gh0/WInFEJPBdsCA51uhCxdxLIGcYfqPe4Ckv+vXzahd92Ffpl+UyRboOBvwQbgjK7MwNJIGcQf6BeAReVHX3ele7iiqr9HiBEz0wDbnO6CPFpEsiZ5aTi8TOuKhgwarTThYicd2WwIfgVp4sQnySBnCH8gfqxBQPHXO0fd5xszyP6yx/kIl9mkUDOAP5AfZmrqPT7ZUede6KMqBD9yA88EWwIljtdiLBJIDvMH6h3A98sP+ZLZ7oKi+WFIfrbWODPThchbBLIzjujJHjKxd7qYSOdLkTkrZnBhuA3nC5CSCA7yh+oP8xTMfgy39ijpzhdi8h7twYbgrIFmMMkkB3iD9SXAleUHXXuUcrl9jhdj8h75cBdTheR7ySQHZAab3xh8fgZR3jKaoc7XY8QKWfLUDhnSSA7Y7y7pPp0/yHTpjpdiBB7+L9gQ3CA00XkKwnkfuYP1PuAr5cdfd4U5fbKOhUi01QDdzpdRL6Svsv+d44vcMwUb+WQMU4XIj5boj3B+nvWY3aaoKByRiU1p9bQOb+TLU9uIb4xzpifjcE3yrfX40NLQmycuxEsqDyhktqzagFYe9daYutilE4uZdAFgwDY8tQWioYXUXZEWb/9fJ/hgmBD8JzmWc1POV1IvpEWcj/yB+rH4vaeWXzo9COcrkXsn3IrBl08iMDNAUb/12i2/XMbsfUxCocVMuLqEfgP8e/zWG1pNjywgbrv1zH2prF0NnUSWx8jtjYGQOC/A0SWRUhGkhjbDaIropkSxjv8KtgQlAZbP5NA7if+QL0X+EbpxFNGuQp8MgEkC3grvPjq7Nav2+emcEghZodJ0ZAiCgd/dm9TdEWUwoGFFAwowOVxUV5fTmhhCNygDY22NNrU4IItT2xhwPkZ120bAL7tdBH5RgK5/0xzFZWMLBoxSS7kZaHE1gSx1TF8Y/bePbEno8PAW7VrhUtPpQejw6BoSBHeKi8fX/8x5UeVk9icAMA38sDO28+uDzYEZQnYfiRvSfqBP1DvBy4snXzmGOXxFjldjzg4yViSNXesYdAlg3D7DnCpEb3vbw3+yuCdX6++bTVDLh3ClsYtxNbGKDm8hKoZGbPy6gDgGuB6pwvJF9JC7h+fd5fV1hYMDkjfcZbRpmbtHWupmFZB+dQD72nyVnkxthk7/212mHgrP7kmfNd7XfhG+bDiFvH1cUZcOYLtb27Hiltpqz8NfhBsCA5yuoh8IYHcx/yB+grgnNIpZx4mM/Kyi9aa9X9cT+HgQmpOrzmoY32jfMQ3x0lsTWCZFp1NnZROKd11blPT/kI7NWfUYCUse29oAG1/L4MUAz93uoh8obTOqP/5OccfqP+Kt2bEeRUnzDpfdo7uG1onddn469L+3HYv62blTSspHFbIjv91Ay8YiDY1Gx7cQDKUxOV34Rvho+6HdRgdBuvvX0/d9+sACC22h71pS1M5vZIBM3dduGt7rg13sZvK4yvRWrPurnXE1sconVjKoAszrkGaBMY3z2pe5nQhuU4CuQ/5A/WDgZsqTpg1uaB25GSn68lVfRXI4hPubp7VfLnTReQ66bLoW2e5iyvd3urhsiuDyHb/FmwIVjtdRK6TQO4j/kB9DTCteMJJI5XLJbuAiGznA77jdBG5TgL5ICjbgT5nJypPoSocNFbGHYtccVWwIVjgdBG5LOsCWSlVrJSap5RarJR6Xyl1kVLqSKXUq0qpd5VSzymlBqfu+4pS6n+VUvOVUsuUUtNTtx+eum2RUmqJUiqQuv37qXO+r5T6Xuq2OqXUh0qp3wHvAftdLjM17viU4sNnDFKegn3PrxUiuwwGLnK6iFyWdYEMnA5s0FpP0lpPAJ4Fbgcu0FofCfwRuHG3+3u01kcD32PXAPfLgf/TWk8GpgLrlFJHApcB9cAxwLeUUjt28hgH/FlrPUVrvfoAajwGKCgaNkF2kBa55j+dLiCXZWMgNwMnp1q+07FbrBOAF5RSi4CfAsN2u/8Tqc/vAnWpr98CrlVK/RgYqbWOAscDf9Nad2utw6njpqfuv1pr/faBFOcP1HuAs32jp5a6iorlIojINVOCDcEZTheRq7JuooLWelmqNXsmcDPwArBUaz1tH4fEU5+TpH5erfVcpVQT8AXgOaXUN9k1NH9vug+ixCBQUVQ3+bCDOEaIbPIt4BWni8hFWddCVkoNASJa6weBW7C7GGqVUtNS3/cqpQ7fzzlGAyu01r8FGoGJwL+Ac5VSfqVUMXAe8FoPSjzTVVQa95QPOrQHxwqRDc4JNgTl2kgfyLoWMnYL9FdKKQswgCsAE/itUqoc+2f6DbD0M85xEfBVpZQBbAJu0FpvU0r9CZifus+9WuuFSqm6Ay3MH6gfCIz1jzuuVoa6iRxWDJwFPOZ0IblGZuqlkT9Q/wXg/OrTrjrZXVI1wul68oXM1HPEk82zms9zuohck3VdFpnKH6h3ASd5ygfGXcWVEsYi150RbAjKRgtpJoGcPqOASt/Y+jGyhpDIA4XA+U4XkWskkNNnKmAW1NZ95gVFIXLIxU4XkGskkNPAH6h3A9Nd/oqwy18+xOl6hOgnJwUbgrVOF5FLJJDTYxTg89VNGiZrHos84gZOdbqIXCKBnB6HAdpbWzfG6UKE6GcSyGkkgZweRwGdnrKBEsgi35zsdAG5RAK5l/yB+jJgWMGA0UWugqIyp+sRor9orXWRZYXv+M0IWSYgTbJxpl6mGQNQOGz86L5+IG0m2DT3x2jTAMvCP+44KqZ/hejqxWx/+Y/opEHBoLFUn/Fd1F4mCm5+7GfEN7RQNGw8Ay7YtbP71qd/hbF1Nb4xR1F54iwAtr/xMAUDRuEPHNPXP5bIIu6I2TEqHO86IRz1XaKihQPRhwCfBz50urZcIIHce5OAhLdqaJ8HMm4vAy++CVeBD5002fTQjygadQTt825j4MU34q0ayvbXHiTc/E9KJ326a6/s6PPRZpzwomd33pbYshKAIV+/g00P/Qgr3o1lxElsXEbFcV/u8x9JZDYVS3aNCMW3nxiJqvPNaNUYl1UJVO6xFNcJwJ2OFJhjJJB7wR+oV8ARQIe7uHLY/u7fW0opVIEPAG2ZYCVRLhfK7cVbNRSAorrJdL79l70Gsq9uMrE1Sz55TpcHbSbQ2kInTVAuOl97kIrpX+3rH0dkorgVHhqKdxwfierzjUjFeJdVBthdcfvu4Jy+z++IgyKB3DuDgBJPxWCzv3YG0VaSjQ3fw+zYSOkRX6Bg8CFoyyS+sZXCwQEiLW+Q7Go74PN5a4bjKa1l45++S8nhn8Ps2AhAgVyfzA8JKzIwlNg2rTuaPNeMlB+pkhVACXAwV5gGM6d8NHM6V/RRlXlDArl3hgIUDBg1uL8eULncDLnsdqxYmC1/uxGjbTW1M39Ex0v3oJMGRXVHwEEuNFd18rd3fr3l8Z9TddpVdL75KIktKymqm0zp5NPT/WMIh2jTitd0JdqO7o6Z5xiRkmMwqt1K2Y2J3o2gnwhIIPeSBHLv1AFJT+XgQf39wK6iEoqGB4mueI/y+vMZ9JVfAhBd+R7mtvU9Omek9W0KBgXQRoxE22pqz53Npod+TPHhM3B5i9JZvugnOqkTFaFE25HhWGKmESk5USeqPUrZ/Vtq53/SYQLwZLpOlq8kkHvnECDsLqnplxZyMtKJcrlxFZVgGXFiqxdRVn8Bye7tuIsr0KZBV9PjlE07+H0oddKka0EjAy74GWbHBna+ULWGpAne9P4som/opDZLw8bWyeFY/Kx4xH8y8ZpCe1MHW99NJA321YnziQRyD6XWr6gDNrv95f0TyOFttM27DbQF2sJ/6HT8Y4+m4+U/Elk+H9CUTj4T38hJAMQ3thJe9A+qz/gPADY99COM9nVoI8a6O2dRfcZ/4Bt9JACh9+ZRMuEkXN4ivLWjAM2G+67EN2YqrqKS/vjxRA9oS1vFYWPrhHAsekY86jvDitUUu+xd19PcAt6fCf31QLlMFqjvodTuIDe5iyvbqk+/+kdO15PP8mmBeq21LgqbbYeFY5HTYtGCs3WsplxlxPsXEyhmTmfC6UKymbSQe24QgLdmpKx2JfqM1pqCSLI9EIqFTolFveda0Zoaxa7fucz5M+QBDgWW7O+OYt8kkHtuBIC7pEp2TRBp5Y6YHaND8a6TolH3uclo9VCXrgaqgUwK4L0JIoHcKxLIPVcHRNy+spFOFyKymyuW7BzRFe88MRpVXzSjVaN2zIaDbFttRgav95IEcs8NAGIuX6m0kMVBUfFkeEgo0XFCJMp5RqT8MJdVDti/R9kVwHvq9+GfuUYCueeqgW2uQr+s8CY+W8KKDAzFtx0biSXPMyLlU3o2Gy4bSCD3kgRyD/gD9UWADzBVgU9ayOITtGnFaroS7fX2bLjSY5VZBaRjNlymk0DuJQnknikDLACX1yct5DynTZ2oDCXapnZHjbMS0eIT2XM2XN6QQO4lCeSeKQe08hZ5lMfrc7oY0b90UptlYWPr5HA08YVE1Hey3m02XP9Oxsg0A50uINtJIPdMOaBcRcWFThci+p62dLI4bLQF7dlwRWdasVrfjtlw0JfTkbONnznlZczp7HK6kGwlgdwzZYBLeX2ZMENKpFlqNtzW8eFY5PRYtPAL9my4Xa0/lwTwZxgISCD3kARyzxQDSZenQJ6/HLBjNtwhoVj4lFjUe44Vra5RDNh5B8nfgyFdeL0ggdIzhYClvIXSQs5S7ojZMcaeDec5LxmtGpw9s+Ey3cEtxi0+QQK5Z4oAS3kKCpwuRBwYVyzZObIr3jkjGlXnZfdsuEwnmdIL8uT1TBGQVB5pIWcqFU+GhnbFO6ZHY5xvRCoPzZ3ZcJlOMqUX5MnrGbvLwuOV5y9TJKzIoFC8/djuqD7PiJZPdiXLgVJAArh/yWuiF+TJ65kiICnjnTLDpYkH1v58/boR7JgNJwHsJMmUXpBf3Z4pACydNEynC8l3R3b/M35d6fPDna5D7CSB3AsSyD2jALSZMJwuJJ8Ni39s3F96n6vAJeMiRG6QQO6ZBODSZlwC2SE+M2Q97LkhXlaQEdsXiV06nS4gm0kg90wMcGtDAtkR2uJe4yfh4cWG7L6aebY7XUA2k0DumRjgsiSQHfHT0M2h4yq3ySp7mUkCuRckkHsm1UKOSiD3s3O7Hol8o/b9UqfrEPskgdwLEsg9Y7eQ4xHZ8rwfBSPz47+sbixQMtowU0WZ0ymviV6QQO4Zu4WciBo6acadLiYf1CQ2mH8u/q1V4JZhVRlMWse9JIHcMxFSz51lROWqch/zWjH9oPqvaGWhJSuJZTYJ5F6SQO6ZEKABdFwCua/9NnJd6NDSqPQbZ76NTheQ7SSQe6aT1J56VrxbFuPuQ1d2/TZ8Rs1GGVGRHVqdLiDbSSD3TCep2XpWLCwt5D7yudDfYz+ofavY6TrEAVvudAHZTgK5ZzpJPXfJSKcEch8YFf0g8bvKB10uGVKRTaSF3EsSyD0QaW2KAVHAk+zukEBOsxKzw5pbdLPh8yAbAGQXaSH3kgRyz7UBhca29e1OF5JLlDb5k/mT8GBfUroqsosGPna6iGwngdxzW4CiZNeWsGXEw04XkytuDN8QmlrRJRfxss865nTGnC4i20kg99waUjvsWpHOTQ7XkhMu6bo/ckntchnelp2WOV1ALpBA7rk1pJ4/M9wu4y97aWr3a/Ebal4odLoO0WMLnC4gF0gg99xGUpNDzO2bJJB7YXB8lfnH0ru0xyVbyGext5wuIBdIIPdcG5AE3IktKyWQe6gw2a0fcs+JlRXoIqdrEb0igZwGEsg9FGltSgKrgGJz27rt2jSiDpeUle6K/yQ0uiQhC81ntxXM6dzidBG5QAK5d1pIbTWfjHRKK/kgXdP1y9DnqtpkREX2k9Zxmkgg984qdlzY275xlZOFZJszQ3+NXlG7UFrGuUECOU0kkHtnZ6s4vql1hZOFZJPDIgsTt1X91SPTonOGBHKaSCD3zibsxeoL4us+2KCThixWvx8ViS3JB/y3moVu2S06F2itO4DFTteRKySQeyHS2mQBS4AKtKXNUNtKp2vKZG4roR/kukhNkeV3uhaRHkqpeczpTDpdR66QQO69RUARgLF1tax29RlujfwsNKGsW2bi5ZZGpwvIJRLIvbdzQZXY6sUSyPvwza7fd59bs0ZGVOQQrXUCeNbpOnKJBHLvtWMvNFRsdm4OJaNdsq7FHo4PPx/7Sc2/ZD+83PMyczpDTheRSySQeynS2qSBt4EKAGPrqqXOVpRZRsRbjbvL/6TcLiW/azlGKfWU0zXkGnmRpMdSUs9lZPn8Zq21w+VkBr/Zac31/ne82IssGpRjtP1LLv3HaSaBnB4rsXcQKTQ7NnQmuztWO12Q47TFfcZPwsP8hkz+yE0LmNO53ukico0EchpEWptM4BWgBiCxcdkSRwvKANeHbgxNq9wuF/FylFLqz07XkIskkNPnHbCXj4y0vrVUW1bejs38YtdDkcsGfCjD23JUanTFQ07XkYskkNNnDfZoixIrGoqbnZtanC7ICZO7347/T/U82Zw0h2l4ijmdHU7XkYskkNMkNdriJaASIL7ug7zrtqhNrDfvL7nd8rrxOF2L6Dsupe5xuoZcJYGcXguxn1MVWd7UaiViXU4X1F+8VkzPVf8VrSzUMt44h5mWXgW86HQduUoCOY0irU1bsTd7rMBKWvENH73tdE395Y7otaFAaUz6jXOcW3EHczplXGcfkUBOv+eBcoDupS+/p5Nmzq8A992u28KnVW+SERU5ztI6rpS63+k6cpkEcvotAbYDfisWiie2rHzP6YL60smhp6PfrZ1f7HQdou8lLe5jTuc2p+vIZRLIaRZpbTKAp0mNSe7+4JW3tbYsZ6vqG4Ho+8btlQ+7ZaH53GdpbXjd6ian68h1Esh9423AALzm9o1dZseGnFvfotxoTz5Y9L8JnwcZ4pYHYiYPycy8vieB3AcirU3dwAvAIIBIy5tvOltRerkskz9Z13YP9CWlqyIPWFon/V41x+k68oEEct95BVCAK77ho03G9o0fOVxP2tzcPadrSnlILuLliZjJI8zplPVZ+oEEch+JtDa1YXddDAQIL3nhxVzoS/5a173dF9WukDDOE5bWlt+rrne6jnwhgdy3ngEKALexdVW7sXV1Vo+4qA+/Gr++5qUip+voKzFTc/Q9YSbdFebw34W5/uUYAC+tNDni7jATfhdm1pNRTGvvw3B//EKMCb+z7/fo+8bO27/yRISJvw9z7T9jO2/7xatxnvrI2NtpMkrM5K/M6fx4//cU6SCB3IcirU0bgJeBwQChRf94RSfNhLNV9czQ2Arj3rI/aI/LXkApFxW64aVZxSy+vIRF3ynm2Y9N3lxrMuvJKI9c4OP9fy9hZLmiYdGng3TeMoP3NiVZdHkxTd8s5ldvxumKa5ZstteYWnJFCa+tSdIZ02wMWczfkOScQzN7423T0nG/V/2n03XkEwnkvvc0oAFvMtTWHd/YknUX+IqSIesh7w3x0gKds61jAKUUJQX2CD7DAiMJbmUH9SHV9t+hU0Z7+OuH5qeO/WCrxYkjPXhciuICxaSBbp5dbuJ1QdQAS2sSSY3bBT97Oc4NMzJ/zf7OmL5NRlb0LwnkPhZpberA3lkh1Up+9i3LiIedrerg/CF+bbiuOJEXC80nLc3ku8IM+FWIU0Z7OHqoG8OCBRvslu7jH5is7fr0pYBJg9z8Y7lJxNC0RSxeXmWyttPisFo3I8pdHHF3NxeO97J8m4UGpgzO7Dca3Qm9udrvkr7jfiarcvWPfwKnAT4d747GVi18yR84ZqbTRR2I2V03h04Y0J43F/HcLsWiy0vYHtOc92iEpVstHvmij/98Lkbc1Jw6xoNnL82YU8d4eGd9kmPv66a2WDFtuHvn/X5z+q43Fmc/HOHus4q48V9xFm9OcspoD986MvOGcocT+qrim7qysnstm0kLuR9EWpsiwGPsHHHx/EIz1L7K0aIOwMzQY5Hv1C7JywWDKooUM0Z6eHa5ybThHl67rJg4JrpfAAAOeElEQVT53yrhhJEeAlV7f9lcd0Ihiy4v4YWvFaM1BKo/eb+nPjKYOthNd0Lz/tYkj33JzwNLDCJGZq3Vsy2q3xp4S+hxp+vIRxLI/ectYDOp9ZJD7z7dqK1kxl5mHx95N3FL1ZPefJoVvbXbYnvMDseooXlxpcmhNS62dNtdFHFT879vxLl86qdbtElL0x6x77dkc5Ilmy1OHbPrDaiR1PxfU4JrjisgYtgD1AEsDYkM2lsmaemkS3Gp03XkK+my6CeR1ibDH6i/F/gp0GW0r+mIrVnysq9uyqlO17anKmNT8gH/bWaBG7/TtfSnjWHNrCcjJC07KC883MtZh3i55vkYz7SaWBqumOrl86Psl82CDUnuWpDg3pk+DAum3x8BoKxQ8eD5PjyuXX/M7nwnwaxJXvxexcSBLjQQ/H2YM8d6qCjKnD96WyP6zkG3hJY5XUe+UrJlff/yB+ovAU4C1qKUqj7t6m+4iyuGOl3XDh4rrhvNK8PjyyJ52VWRz9oj1spqv+sQ5nR+ehiJ6BfSZdH/ngS6gDK01qGFf38qkzZEvS3yXyEJ4/yTSGqzdZt1oYSxsySQ+1nqAt992MtzuhKbl2+Nr//gXw6XBcDlXXeEz65ZlzcjKsQuLW3Wr465t3uB03XkOwlkZywFXgWGAHS98+RrTo+6ODH8XOya2jfyqs9Y2NZ1WYuveyl+ndN1CAlkR6R2qH4MiAKlaEt3vvXo405NGKmLtRi/L29QbqXk9yHPRAwdaW23zm1sybCxd3lKXoAOibQ2hYDfY3ddeJKhtu7wkucf1/18lbXY7LDmFtyY8HvJ/Lm8Iq201ny41frPzzV0r3K6FmGTQHZQpLXpA+AJYDhAbNXC1fE1zS/11+MrbXK/eW14iN+Uhebz0OLN1iNH/iH8B6frELtIIDvvGeyNUe3+5AVPvm52bu6XcaA3hP87dHRFp1zEy0MtbcmlP3s5fqnTdYhPkkB2WKS1KQncC0SAcoDtbz7yNyse6ejLx72w64Hur9Uuk+FteWhjyNry58XGFxpbjLjTtYhPkkDOAJHWpk7gDuxp1QVWpDPW2fT4Q9pMRPri8aZ0vxG/sfof0mech0JxHXnyI/NLN74Wly2ZMpAEcoaItDYtB+YCwwCXsXVVe9d7z8xN93oXAxNrzD+V/s7yumXafL4xkjo5r9X8wRXzohkx7l18mgRyZnkRe6nOkQDxte+v7176ctpGXhQmu/VD6vpoeYH2peN8IntYWusXV5h3z2027na6FrFvEsgZJNLaZAEPAe8CIwAiy95cFl3e9Ew6zn9n7KehsaVx6TfOM1prnl1uPvH7Bcb3ZLxxZpNAzjCR1iYTuAdYQWqXkfCS59+LrVv6Sm/O+/2uW0MnV2+WERV56O+t5j/vWmB8vbHFyNjlXoVNAjkDRVqbosBvge1ALUBX019fjW9a3tST850eejJ6Ve2CvNiCSXzSs8vN1+5+17ikscXocroWsX8SyBkq0trUBdwKWKQWte98Y+6z8Q0trx/MeQ6JLjH+r/IxjyufVpoXALy00mz63TuJixpbjC1O1yIOjARyBou0Nm0BbgEK2BHKbz36z9ja918+kOMrjK3Jh3y/Mgo9ZPZ+8yLt/rXafO83bycuaGwxNjpdizhwEsgZLtLatBr4H+zdXaoBuuY/8a/oqkXPf9ZxLsvkz9Z13bVFSVnBLc+8uMJ855Y3E+c3thjrnK5FHBwJ5CyQCuWbAY29GBGhdxvfinz8zrx9jYj7ZffPuiaWh+UiXh6xtLYebjZe/W1T4uLGFkMmfmQhCeQsEWltWgfcBBjAAIDwon8siCx780mtLWv3+17WdXf3BbWrJIzzSCKpjd+9k/jHw+8b32xsMVY4XY/oGdlTL8v4A/W1wDVAKfYu1hSNOrKudOKpFymPt+jY7pdiD1TeU+B2ydrG+aI7oSO3vZ14Yv765I8bW4wNTtcjek4COQv5A/VVwA+xh8StB/BXDx55+KRJF86r+Y27RNY2zhvtEWv7/7yeaGhpt37e2GL06YJUou9JIGcpf6C+FLgCGO8muXGsWn/a2KKuTbee5AqMKHeNdbo+0feWtSfX3fpm4u6NYf3rxhajTxaiEv1LAjmL+QP1XuCSEWz+9hDVbhYTe9ntQs0+vuDzxwzzHO90faJvWFpbf281F93zrnGnhgdkBl7ukEDOcv5AvZqkll9cQ9fpStEGhAAuPNxz2JfGe88u9ChZSCiHhBM6dHtT4vW31iXvBp5ubDGs/R4ksoYEco6YOc57KHA19siZzQB1Far0mmMLzx1e7hrtaHEiLT7eZq29+fX4C1u69a8bW4ylTtcj0k8COYfMHOetBa4E6oB1gKmAK47yHn3yaM8pHpeSNZCzkKW19dxyc9FdC4zHNdzd2GJsc7om0TckkHPMzHHeAuAsYCb24kTbAaYMctVcXV/wxRq/a5CT9YmD0xaxtt7elGhauMm6H2hsbDFMp2sSfUcCOUfNHOcNAJcDFdhD46wCN64fHlvwuaOHuo+TxYYyW9LSyZdWJhfetSDxrmFxR2OL8b7TNYm+J4Gcw2aO8xYDFwEzsPuVuwGmDnEN+MaUgjOGlrnqnKtO7Mu6Lmvt7U2JhR+2WS8BDza2GO1O1yT6hwRyjps5zquAScC3sFeN24i9pCcXHe4Zf86h3lNLClS5gyWKlKihI098aLzz6FLzA6ABmC87fOQXCeQ8MXOctwL4InAC0AW0A5QU4LnyqILj64e5j5OLfs4wLW28tTa58N73Eh93xHgVeExm3eUnCeQ8M3OcdwzwNWAUdjdGBGBctaviO1O9p4ypdI2X7uX+kbS0tWSzteiuBYnWjWG9Frgf+FBaxflLAjkPzRzndQPHAJcARcAGIAlw1BDXgIsneE8YU+UaLxf++obWmpZ26/27FyQ+/LhDdwCNwIuNLUb8QI5XSoW11vvdkkspVYv9//YqrbXsNp0FJJDz2Mxx3hLgC8BpgAlsItW/PGmgq/qSoHf6uBpX0KVk5bh0sLTWq7brZX9alFi6aJO1HXgWeL6xxeg8mPMcRCD/O/BlIKm1ntGjokW/kkAWzBznHQycARyP3VLelPrMoTWuiq9N9B4/vtY12e1SbgfLzFqJpI43b7YWzm02VrVus2LAK8C8xhajrSfn2xHISqkZwBygDZgAvAt8Vade1Eqp14AfAHOBE7XW63ccD9wJnAx0ANcCvwRGAN/TWjcqpS4FzgMKsbu35mqtf96TesWBk0AWO80c5x2A3Vqegb07ySbsljMjy1XJF8d7pxwx2D2lrFBVOldl9uiK622vrzHfmdtsbOyyOyMWAH/r7ZrFewTyU8Dh2F0TbwDXaK1fV0oNB17SWgeUUjcBbVrrX6eO18CZWut/KKX+BhRjv1MaDzRorSenAvlm7KCPAO8Al2qtF/SmdvHZ5Kq62Cm1O/EDM8d5/w6cgt2CUsCW1Z06/Ou3Eq8peO2MgGf0SaPcR4yudB0qreZPsrS21nXpFc9/bC5+usXs0vY7jdeAlxtbjLV98JDztdbrAJRSi7Cnzb8OXAw8lrrPI8B9wK9T/05gd5cANANxrbWhlGpOHb/DC1rr9tS5n8B+ByWB3IckkMWnpCYiPDJznPc54ETsYB4EdGto/3urueLvreaKQSXK96Xx3klTh7iPqPSpWkeLdpDWmk1hvea9jcnmp1rMDZvC2os9tHAe8HZji9HVhw+/+4XAJLte018GBiqlvpL69xClVEBr3QoYetdbY2vHObTWllKfGPq459tneTvdxySQxT6lxsI+OXOcdx7229mTsN/CAmzdFNbR2+cn3gbenjTQVf25UZ5Dx9e6xg0oVsPyYYTG1m5rw6JN1vtPLzOWr9quC7FX2lsHPAM0O7VOsVJqHFCstR66220/x241/+IgTnWKUqoKiALnAl9Pa6HiUySQxX6lgmUxsHjmOG8NcDRwKvYWUibQvniz1b54c+IN4I3BJcp/+ljPIZMGuccNL1NjvG7lda769ImZOrKuy1r1UZu18tVVyTUt7ZYXcGN3AbwAvAeszYBxxF8G/rbHbX/F7ro4mEB+HXgAGIt9UU+6K/qYXNQTPZIay3wIMAWYBvhT3+oAwjvu5/fiOX6Ee+jhte7hdRWuYYNK1DCfVxX3f8UHL5HU8Q0hvbq13Vo5f31y5fz1yS4NldgNmW7svuEFwKpcWyg+dVFvqtb6KqdryScSyKLXZo7zuoDhQBA7nAenvhUBOoFPvHU/rMZVeeQQ97BAlWvY0DI1tKJI1RS4laMbs0YMHWqL6M2bwnrLmk5rc0ubtXnBhmRnUlOOPXlGY/+xeRc7hD9ubDGSTtbclySQnSGBLNIqtZhRDXAo9qJGhwE7tpFKYl/sCrPHBaLBJcp/SLWramSFq2pQiaqq9auqiiJVWVqoygvd+Hs7miNm6mjUINxt6FAorsOhhA61R3TXxx3W5sWbrC2bu3UMu5VfnKrXwv6DsghYAqwE2jOgO0LkMAlk0ad2C+hhQAB7zOxw7EDeceEvutvHXlud5YUUDCxx+SqLVFFxAQXFXlXg91KgFCppYZn2R9KwsExLJ00Ly0iS3NKtI2u7rHDM3Hlehb3qXXHqQ2GHrwt7LO9K4MPU500SwKI/SSCLfjdznLcIu1ujGhiAHdBDgYGAl10BuUMS++Lh7p/BDtPdP1y7fe3d7Vw7fsl3fL8DWAV8jL14/2agTXZvFk6TQBYZI9WaLgWqgHLsroPi1NdlqX/7sLsWXNjhbGKH7o6vk9h91tth5y7c4d0+d+dy36/IbhLIQgiRIWQVLyGEyBASyEIIkSEkkIUQIkNIIAshRIaQQBZCiAwhgSyEEBlCAlkIITKEBLIQQmQICWQhhMgQEshCCJEhJJCFECJDSCALIUSGkEAWQogMIYEshBAZQgJZCCEyhASyEEJkCAlkIYTIEBLIQgiRIf4fWg+eLoNzK8sAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "in_t = np.sqrt(in_sub**2 + in_inamp**2 + in_sensor**2)\n",
    "\n",
    "labels = 'sensor', 'InAmp', 'Subtractor'\n",
    "sizes = np.array([in_sensor, in_inamp, in_sub])\n",
    "sizes /= sum(sizes)\n",
    "\n",
    "fig1, ax1 = plt.subplots()\n",
    "ax1.pie(sizes, labels=labels, autopct='%1.1f%%',\n",
    "        shadow=True, startangle=90)\n",
    "ax1.axis('equal')  # Equal aspect ratio ensures that pie is drawn as a circle.\n",
    "\n",
    "print_white_noise(\"Total\", in_t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Low-frequency noise\n",
    "\n",
    "Dominated by the DAC and voltage reference"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [],
   "source": [
    "def print_low_freq_noise(name, I_n):\n",
    "    \"\"\":param I_n: peak-peak current noise as a fraction of the full-scale coil current\"\"\"\n",
    "    print(\"{}: {:.2f} ppm pk-pk ({:.2f} nT pk-pk)\".format(\n",
    "        name,\n",
    "        I_n*1e6,\n",
    "        I_n*field*1e9))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "sensor: 0.07 ppm pk-pk (2.10 nT pk-pk)\n"
     ]
    }
   ],
   "source": [
    "# sensor noise\n",
    "print_low_freq_noise(\"sensor\", In_sensor)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "reference: 0.14 ppm pk-pk (3.95 nT pk-pk)\n"
     ]
    }
   ],
   "source": [
    "# reference LTC6655\n",
    "# NB there is a newer -LN part number with lower noise\n",
    "# the reference is filtered with a pole at ~16Hz, so we don't need to worry about the integrated white noise\n",
    "V_ref = 2.5\n",
    "i_dac_fs = V_ref/G_sense  # equivalent full-scale coil current\n",
    "G_dac = (I_coil/i_dac_fs)  # conversion between DAC full-scale value and the coil current\n",
    "\n",
    "Vn_ref = 0.25e-6  # 0..10Hz integrated noise (pk-pk as a fraction of full-scale)\n",
    "In_ref = Vn_ref*G_dac  # pk-pk noise as a fraction of coil current after scaling by dac\n",
    "\n",
    "print_low_freq_noise(\"reference\", In_ref)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dac: 0.10 ppm pk-pk (2.82 nT pk-pk)\n"
     ]
    }
   ],
   "source": [
    "# DAC (AD5541)\n",
    "# ignore DAC white noise density as small compared with the reference (and is filtered out anyway)\n",
    "\n",
    "Vn_dac = 0.134e-6  # 0..10Hz pk-pk DAC noise (V)\n",
    "In_dac = (Vn_dac/V_ref)*(i_dac_fs/I_coil)  # pk-pk as a fraction of the coil current\n",
    "\n",
    "print_low_freq_noise(\"dac\", In_dac)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "OpAmp voltage noise: 0.11 ppm pk-pk (3.09 nT pk-pk)\n"
     ]
    }
   ],
   "source": [
    "# OpAmps\n",
    "\n",
    "# ADA4004-4 voltage noise is 0.1uV pk-pk 0..10Hz\n",
    "# vn_opamp = 0.15e-6  # 0..10Hz pk-pk DAC noise (V)\n",
    "# in_opamp = (vn_opamp/G_sense)/I_coil  # pk-pk as a fraction of the coil current\n",
    "\n",
    "# LT1007 voltage noise is 0.1uV pk-pk 0..10Hz\n",
    "Vn_opamp = 60e-9  # 0..10Hz pk-pk DAC noise (V)\n",
    "In_opamp = (Vn_opamp/G_sense)/I_coil  # pk-pk as a fraction of the coil current\n",
    "In_opamp_T = In_opamp*np.sqrt(1**2 + 1**2 + 2**2)  # contribution from all OpAmops (noise gain of final InAmp OpAmp is 2)\n",
    "\n",
    "print_low_freq_noise(\"OpAmp voltage noise\", In_opamp_T)\n",
    "# NB this is the noise for each OpAmp, so the total contribution will be several times this...\n",
    "\n",
    "# ADA4004-4 low-frequency current noise is not well specified\n",
    "# ballpark as follows:\n",
    "# - current noise converted to voltage by resistors, which are all ~1k\n",
    "# - typical bias current is 40nA, so gives offsets around 40uV\n",
    "# - assume 0.1% stability for this (should be less thank that in a well temp controlled environment)\n",
    "# - gives around 40nV fluctuations, which are small comapred with other contributions\n",
    "# - so let's ignore this\n",
    "\n",
    "# LT1007 current noise is still only 1.5pA/rtHz at 10Hz so we can ignore this safely!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total : 0.21 ppm pk-pk (6.12 nT pk-pk)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWQAAADuCAYAAAAOR30qAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl8XFX9//HXuXcmk8xM1iZpmm5pISRDC4WyRARkRwRsQQG/gloQN3BBcKsiLig/3L+IIKB8QREVZCuRImuRfZW2hjJJ9zZNmyZplsk26z2/P+50gbbQJpO5k5nP8/GYR5qZe+98Jpm+c+bcc89RWmuEEEI4z3C6ACGEEDYJZCGEyBASyEIIkSEkkIUQIkNIIAshRIaQQBZCiAwhgSyEEBlCAlkIITKEBLIQQmQICWQhhMgQEshCCJEhJJCFECJDSCALIUSGkEAWQogMIYEshBAZQgJZCCEyhASyEEJkCAlkIYTIEC6nCxBif9QsXFwMlAGlgAf7PWze9tQv4tMGOhQQA+LJr4PA5kBzcMipeoXYHxLIIiPULFxsAFOBA4EDdvk6mZ0BXMJe3rMVw73NQP2eHgvWB0JAG7B5L7dVgeZgZwpfjhAjIoEs0q5m4eICYC5wdPJ2GDADu8U7IoNWoqxg7w8XJW+BvW0QrA9sBN4AXk/e3gg0B/tGWo8QIyGBLMZczcLF1cDpwDHA0Vrr2UqplL73LK3dozzEtOTtY8nvdbA+sJpkOCe/vhZoDkZH+TxC7JUEcpZTSilAaa2tdD1nzcLFbuBY4Ayt9UeUUoe+q6aUP6dCJVJ+SKhN3i5M3hcK1gceBRYBjwaag/0pfk6R4ySQ00Ap5QP+AUwBTOAnwGrgN4Af6AIu1lpvUUr9G3gVOAm7z/RSrfXzSqlZwJ1AHvbomI9rrVcppa4CPpt8qtu11jcopWqAfwHPYLdKzwE2jOVrrFm42Aucm6zrFKVUUfK1j+XTplsR8D/JWzRYH3gaO5wfDjQHtzpamcgKEsjpcQawWWt9FoBSqhg7MOdrrTuVUp8ArmNnsLq01kcrpc4EfgicCnwJ+K3W+q9KqTzAVEodAVwCNGC36F5VSj0L9AB1wCVa68vH6kXVLFxsAqdprS8Czk3+4cm2EN6bPOAjydstwfrAK9jh/ECgObjW0crEuCWBnB5NwK+UUj8HHsEOzNnAk8nwMoEtu2z/YPLrf4Ca5L9fBq5WSk0BHky2jo8DHtJaDwIopR4EjgcagQ1a61fG4sXULFx8qNb6EtAXKWVU5EgAvxcD+GDy9vNgfeAp4FagMdAcjDtamRhXJJDTQGu9MtmaPRO4HngSWKG1PmYvu0SSXxMkf0da678ppV4FzgIeV0p9DrtVvDeDKSk+KTks7aPaSnxTGeZxya7pVD5FtlDAaclbW7A+cDtwW6A5uOW9dxNCrtRLC6VUNTCktb4b+BV2F0OFUuqY5OPuZB/xex1jJrBWa30jdgv4UOA54ByllDfZXXAu8Hwqa69ZuNg//Tv/vEIn4muBRcowj0vl8bPcZOwup/XB+sCdwfrAIU4XJDKbtJDT4xDgl0opC/sKssuwrya7Mdmf7AJuAFa8xzE+AXxKKRUD2oFrtdbdSqk/Aa8lt7lda700eVJvVGoWLi7XVuK7oL6gDMOPKX+7RyEPuBi4OFgfeAL4ZaA5+JSzJYlMpLTWTtcgMkjNwsVFVix8tTLdX1WG+R7XWmSWuxZ9p6uCRLnTdeyHp4FvBJqDy50uRGQOaSELAGoWLs63osPfVKb724Y7v9DpenLAKcCbwfrAX4CrA83BNqcLEs6TQM5xNQsXG1Y0fJkyzR8beQUTnK4nxxjAAuD8YH3gN8DPA83BAYdrEg6SjsEcNvWKe463osPNRl7+Tcp0Sxg7xwt8H1gdrA98IVgfMJ0uSDhDAjkHTb3inrKpX/3rIiPf/5yRV1DrdD1ih4nAbcDyYH3gVKeLEekngZxjplx255eV27PB9JXMlws6MtYs4MlgfeCmYH1g3JxYFaMnfcg5ovqzN00xfaUPuoorj3K6FrHPvgycGqwPfCrQHHzD6WLE2JMWcg6YtOCGS12l1c2mr0TCePypA14K1geukb7l7CeBnMUq5n3bN/kLtz2SV3Xg7Ybb43O6HjFibuBa4IVgfeBAp4sRY0cCOUtVXfTzD+XPOHy1u2zKWdJXnDU+ACwL1gc+73QhYmxIIGcZb22DMenTv77WU13/tFlQVOV0PSLlfMAfgvWBh4L1AX+qDqqUmqKUelgptUoptUYp9dvkNK/vt9+5SimtlNrjeoZi/0ggZxHfwScWFh974aOeyfXXKNMlJ2yz2znA88H6wOTRHii5qsyDwCKtdS1wEPbCCdftw+6fBF7AnrRfjJIEcpYoOWFBbcnxF/3HU3Xgh52uRaTNYcCrwfrAnFEe52QgrLW+E0BrnQCuBD6rlLo82XJ+TCnVopT64fadlFJ+7KW6LmWXQFZKnaiUelYp9Q+l1Eql1M+UUhcppV5TSjUppQ5IbvcnpdStSqnnk9udPcrXMe5JIGeBCWde8eHCOR9+xV1aLRd55J7J2Cf7zhzFMWZhL4awg9Y6BGzEHhp7NHAR9h+A85VSRyY3Owd4TGu9EuhWSs3d5RBzgCuwZzr8NHCQ1vpo4Hbgq7tsVwOcgD3P961KqfxRvI5xTwJ5HPPWNqiKed++whf40CLTW1zmdD3CMX6gMVgfGOlyXQrY07SP2+9/Umu9TWs9jN21sX1O7E8C9yT/fU/y++1e11pv0VpHgDXAE8n7m9i5Cg7AP7TWltZ6FbAWyOm+aOlnHKe8tQ0eb92xP/HWHXul9BcL7GXAbk4Oi/tmoDm4P6uMrwA+vusdyUVqp2KvWvPusNZKqQnYXR2zlVI6+fxaKfXt5DaRXba3dvne4p25s9ux96PurCMt5HHIW9vg89Yde6vv4BOvkjAW73Il8ECwPuDdj32eBrxKqc8AKKVM4NfAn4Ah4DSlVJlSqgC7m+JF4DzgLq31dK11jdZ6KrCOna3nfXW+UspI9ivPBFr2c/+sIoE8znhrG4p8B59wh+/gkz6jDFOu3BJ7cg7wyL7Og6HtVSrOxQ7HVcBKIAx8L7nJC8BfgGXAA1rrN7C7Jx5616EeAC7cz1pbgGexV2H/ktY6vJ/7ZxVZMWQc8dY2lPoPOfWugtpjzlJytcc7jMMVQ9LhSWBeoDk44pBTSl0MHKm1/krKqtp57D8Bj2it70/1sccraSGPE97ahnLf7FP+WlB7zNkSxmIfnQY8GKwPeJwuROwbaSGPA97ahqqC2g/c6j/k1HlKGRLGeyAt5Pf0T+BjgeZg3OlCxHuTFnKG89Y2TMifduiN/tmnnC1hLEboo8AdwfqAvH8ynARyBvPWNhTlTTro54Vzz5ovJ/DEKH0ae+SEyGASyBnKW9tQ4C6f/qOio865UJnu953kRYh9cGWwPrDQ6SLE3kkgZyBvbYPb9JVeWdTw8c8Z7nxZwkek0vXB+oBMBJShJJAzjLe2wUQZny1qOO/LZr6/0Ol6RFa6PVgfCDhdhNidBHLmOavoiHlfcZdOqna6EJG1fNhX86VsPmWRGhLIGcRb2zCn4MCjr8qffuhsp2sRWS+APfOayCASyBnCW9tQ7S6f/n3/7FM/6HQtImd8Ilgf+Or7bybSRQI5A3hrG3yGx/fN4oaPn6ZMl9vpekRO+XWwPvABp4sQNglkh3lrGwzgs4VHzj/LyPcXO12PyDlu4L5gfaDC6UKEBHImOK1g5hHneKoOPMjpQkTOmgL8LVgfkDxwmPwCHOStbZhqeIsX+GafKh8ZhdNOBa5yuohcJ4HsEG9tQx7wheKG8xoMt2d/JhMXYqz8KFgfmOZ0EblMVptwzke99cef4C6bPNPpQnbV9egNDK95HdNbTPWlvweg8+GfE+veBIAVHsTI91F9ye9229cKD7DtXzcS7doIQPmZV+CZHKDn33cyvPY/5FXOoPzsbwAw8NYSrHA/RUfOT9MrE/vAB9wEzHO6kFwlgewAb21DreEtOd9Xf9zRTtfybv5DTqVw7tlsW/ybHfdVzP/Ojn93L7kdw+Pb477dT/+B/JlHUHHu99CJGDoWwYoMEmkLUv3Zm+j85y+Jdq7HVTKJwbeeovL8a8f89Yj99tFgfeDcQHPw3auBiDSQLos089Y2eIEvFh3x0UOU6c64icPzp87GLNjzFdtaa4aaX8AX+NBuj1mRIcKtK/AfejoAynRj5PsBhU7E0Vqj41GUYRJ67UEKj5iHMqU9kKFuDNYH5LJ9B0ggp9+ZeZPq6twVNQc7Xcj+imxagekrwV02ebfH4r3tmN4itj16A5vv/Brb/nUjVjSM4fHirfsgW/70NVzFE1EeH9EtK/HWynnMDDYFkI8vDpBATiNvbUM1Sp1ZeNgZR4zHVZgG3352j61jAG0liLavofDwM6m+5EaU20PolfsAKG44j+pLfkfZyZ+j7/m7KTn+U/Qvf5zORT+j96V70vkSxL77arA+cLjTReQaCeQ08dY2KOBC36yTDzC9xZOcrmd/aSvB0MqX8dbvOZBdheWYheV4qusA8NYdS3Trmndss/17V+lkBt9aQsU5C4l1biDW3Ta2xYuRMIHbZGxyeskPO33mGPmFR3gPOOoopwsZifD6ZbgnTMFVtOdl60x/Ka6icmLb7NEY4Q3LcZe/cwRV7/N3U3zcRWDFQVv2ncpAxyNjWrsYsaOALzpdRC6RQE4Db21DPrCg8PCP1CpXXkZPON/Z+Ava//JNYt1tbLp5Af3LnwBgMPjcbt0V8f5tbL3vhzu+Lzv1S3Q98is23/EVoh3rKDrmgh2PDa18mbyqWlyFEzDy/Xiq69n8f18GBXmVGTXyT7zT1bJqdfrIqtNp4K1t+Kjpn3Bh2WlfOk8ZMrRgLMiq02Pq8kBz8Bani8gFEg5jzFvbUALMKzzsjFoJ4+yxJRbju1u20JWIo4ALSkr4dGkZN3V1cn9fH6XJNWm/Xl7BCf7d54H/c3c39/f1ooCDPB6uq5qExzD41ubNrIpEOMHv58oKe76fW7q6OMjj4ZRCx0aiLQzWB24PNAdjThWQK6TLYux9OFAc8T9b8fMpHw4/LmevsoRLKb5dWckjM2Zyz/Tp/K2nh9URuy/8M6WlPFQzg4dqZuwxjLfGYtzd28N902tonDGTBPBof4iWcBiARTNm8J/hIfoTCTrjcZrCw06GMcA04DNOFpArJJDHUG3tzMkl9H/yZ3O31tZ6BybeVvLnyUv05ztPDz+x2enaxOhUuFwcnJ8PgM8wmenx0BGP7/P+Ca0Ja01ca8KWRaXLjUspItrC0pqY1hhK8buuTr5SnhEzY343WB8wnS4i20kgj6Fa1Xby8b7WGSdMHJq7/b6ZBYMVfyj5U/US/fnOU8NPSTBngbZYlGA4zKHJgP5bTw/nrFvH1Vu20JdI7Lb9RLebS8rKOGXNak5Ysxq/YXKsz8cBHg+T3G4+vmE9ZxQWsjEaRcOO4HfYAcCFTheR7eSk3hiZV+cuA37xnWPz6o6d5jpub9utHvJ3XB/9ROLp/FPG3djkTOLUSb1By2LBxg18cUI5pxUW0hWPU2qaKODGri4643Gum/TOX21fIsHXN7fx60nVFJomV25u43R/IfOK37k+weWbWvlRVRUP9vXREonwQa+P80tK0vjqdtMCHBxoDlpOFpHNpIU8dk7xmBhzJ5nvebXTgd6Byv8r+b9JT1pf6DgpvERazONITGu+3tbG2UXFnJbs4y13uTCVwlCK80uKaQoP77bfy0ODTHa7KXO5cCvFaf5Clr1ru6f7+5mVX8CQpVkdifC/1ZNpDPUxbDmahXXABe+7lRgxCeQxMK/O7QVOu2CWe0KBW+15arR3qfUOVN5Zcnv1k9YXO06M/HvLGJcoRklrzTXtW5jpyePisrId93fu0o/8VP8AtZ7dh/BOcrlZPjzMsGWhteaVoUFm5uXteDymNXf39vDZsjLCloXCvsxeJx9z2PecLiCbSSCPjSMB9wk15hH7u2Ott7/yT8V/mPS49aWtH4o8K8Gcod4cHqYxFOLVwSHOXb+Oc9ev49mBAX7V2cH8des4Z906XhsaYmHlRAA64jG+uKkVgDkFBZxeWMh5G9Yzf/06LOCC4p1dEX/v6WF+UTEFhkGdx4NGM3/dOg4vKKDIdPy82iHB+oCsjD5GpA85xebVuQ3gZx+cak5eeJxn1CdBmoeKtv40epF+If/4qhSUl7XkwpC0uiPQHLzU6SKykbSQU68OqJhf55qdioPVe0MT7y65pepficvajw2/0J6KYwoxShcE6wO7D7AWoyaBnHqnFeYRq51gBFJ50ICvr+qvJb+vetS6bOsxkRclmIWT/MAnnC4iG0kgp9C8OncFcPg59e5yl6HcY/EcB3v7Jv69+OaqxYnL2xsiL28di+cQYh9c4nQB2UgCObUOB/RRk82UdFe8l1m+3qp7i3838REJZuGMD8oK1akngZwi8+rcCjip0qeGphapA9L1vLOTwfzPxJe3HBV5VYJZpIsCPul0EdlGAjl1JgMTz6l3TTcNlfaf6yG+nkn3Ff92YmPiK1uOjLzWke7nFzlJAjnFJJBTZy5gzZ1kHuJkEYf6uifdX3xD5aLEVzbPjbzR6WQtIuvNCdYHxt1ivZlMAjkFkmOPT5xQoAaq/Coj+tUO83VXP1j8m4qH4l/dfHj0P9JiFmPlbKcLyCYSyKlRA5ScOtOsNjJsOenD/duqHyr6deWD8a9tPiz6prSYRaqd7HQB2UQCOTUOBqw5VeaBTheyN3P9XdWLin5V8UD8a5vnRJZKMItUOS5YHxiTIZ65SAI5NY4G+maUGGkbXTFSR/i7qh8u/mXF/fEr2g6JLO9yuh4x7vmABqeLyBYSyKM0r85dDEw5vMrI9+WpIqfr2VdH+jsn/7P45+X/iH+9bXb0vxLMYjROcrqAbCGBPHoHABw3zZXxreM9OdrfMfmRop+V3xv/etusaJMEsxgJ6UdOEQnk0ZsDRA8oM6Y7XchoNPg7Jj9SeH35PbErNwWiK7Y5XY8YV44J1gcyYp2p8U4CeRSSw93mAt1VfjXZ6XpGSyn4QOHWKY8WXjfh77GrNtVHV3Q7XZMYFzzAsU4XkQ0kkEenAvDOKFEer1s5uk57KikFxxS2T/lX4XVlf419Y1Nd9G1pMYv3I90WKSCBPDrVgDqi2hz3reM9UQqOLdwy5bHCn064O/bNTQdFm6XFLPbmRKcLyAYSyKMzE0gcNMHIykDeTik4rnDzlMcKry39S+ybrQdGW3qcrklknFlOF5ANJJBH52Cgf0pRdgfydoZCHV+4eeoThT8uuSv2LQlmsaviYH1gotNFjHcSyCM0r87tBqYDA+VelVPr3RkK9aHCtqlPFP645E+xb7ceEF0pwSzAXr5MjIIE8shVAarSp/LyXarA6WKcYCjUiYWbpj5Z+KOSO6PfaZ0RWy3BnNskkEdJAnnkJgKqboJR5nQhTjMU6qSi1qlP+39Qckf0O601sTW9TtckHHGQ0wWMdxLII1cJML3EmOB0IZnCUKiTi1qnLvFfU3x7dOHG6bG1Esy5RVrIoySBPHLTgOHqQpXzLeR3MxTq1KKN057xf7/4j9Hvtk6LretzuiaRFtJCHiUJ5JGbDAyVeyWQ98ZQqNOKNkx9xn914R+i39s4VYI5280M1gdcThcxnkkgj0ByQdMqIFySr0qdrifTmQrj9KL10/7tv7rw1ujVG6fENoScrkmMCTf22HwxQhLII1MMmEDC61Z+p4sZL0yFcUbRumnP+r7rvyXy/Y3VsY0SzNlnXM56mCkkkEemFLAA8l14Ha5l3DENjI8Ur532vG+h/+boNRLM2aXE6QLGMwnkkfECym1g5JnK43Qx45VpYJxVtGbac76FvpsiP9hQHWuVYB7/smaSLSdIII+MF1BV/ty8ICTVXAbm2cWrpz/n+47vd5EfbJgUa+t3uiYxYtKFNwoSyCPjxb5KT7orUshlYH60ePX0533f8v428qMNVXEJ5nFIWsijIIE8MiVAoqxAAnksuAzM+cUrp7/g/Zb3hsiPNkyMbx5wuiaxz6SFPAoSyCNTAsQK3EqWPx9DLgPznOKV01/0frPgN5FrN1TGt0gwZz4J5FGQQB6ZYiBmKJTTheQCl4H5seLm6S96v5H/68i1Gyri7YNO1yT2SrosRkECeWQ8gGUq+fmlk9vA9fHi5ukvea/y/DLy0w3lsS3DTtckdiMt5FGQQBkZA9CmIT8/J7gNXOcXvz39Je838q7vuSZUEWsfcromsYO0kEdBrjsfGQPQ0mXhrDwT85OT1hSdl7gq9o0PTuzzRBj2DWP5wmhvRCtvBApiqIIYZn5cGfkJXHmWcuVp3G6t8lyQp5T8ClNMhoKOggTyyJiANpWSFnIGcJu4/zNbFYcm5JXv146WRcEQUe+QjvqGdMw3SNw3pOO+IW15h3XCH0b7wlp7w+CNogpiGPkxZRTEMT0JXHlaufIs3G5UnkvJCd6kqNMFjGcSyCNjAtrpIsRORjjh2/+dDIb95A37yds22uePWzp/SEd8Q8S8Qzrmt4M97h3C8g9ryzuM9kV2hLuRbLkb+Qnl8iQw8yzlztN2uJtKjef/l9KvPwrj+RfvJAMgHNcxpwsRNrNAdWNPieoIy2WooSLyh4rIH+2xzJhlFQzpiG9Ix72DxHx2wCe8QyR8Ya19w8mWewRVEAXv9i6ZOKbHwrVLuHuM9H+Kk0AeBQnkkdEAQzH5eJYp3Ja2nK4hVRJuwxgopmCgePTHckWshHdIR7xDOuYfJO61W+4J/zAJ77C2fGHwhaEgond0yxTElelJ7OiWcedp5XZB3j6Gu5xgHQUJ5JGJAEWDMS2BnCE8WksX0h7EPYYZ8uANpWDW7rywFS8Y0lH/oI567f72hG9Yx33DWAUDlrvScK0zXGptYPRPlbMkkEdmGDBDER1xuhBh81gSyGMtmm+4ovm4+sr2OuVsFfDCV9JZVJaRUQIjEwbMbUNa+ssyRL4leZwhup0uYDyTQB6ZfsDVOaTDThcibAVaEjlD9DhdwHgmgTwy/YArHCcRjms5iZEBCqTHIlNIC3kUJJBHJoQ9FpmBqJaVlDOA19JyyV1m2OJ0AeOZBPLI9JNcUy8UQQI5A3i1JYGcGVY6XcB4JoE8Mr0kxyL3hqWFnAl8WlrIGaCvaUFTh9NFjGcSyCPTB/bEQt3DEsiZwIcEcgZY5XQB450E8siESP7stg5YEsgZwK8teS87T7orRknexCPQ2BKLAgOAe12vHu28NCIFfFqbTtcgpIU8WhLII9cJ5C9rT3QmLJ1wuphcVygt5EwgLeRRkjfxyG0EvNEEVm9YdzldTK7zIS3kDCCBPEoSyCO3BnttPToGdbvDteS8QrTMy+I86bIYJQnkkWsnOfRtU0gC2WkSyI7b2LSgSU5wj5IE8si1kxz6tqrb2upwLTnPjyVdFs56zukCsoEE8sgNJG95b25JbLFkPl5H+dGypp2znnW6gGwggTxCjS0xDawF/B2DOtwzrKWV7KB8pQydRauGjEPSQk4BCeTRaQZ8AOt79XpnSxFYyBqHDtBatzctaJIRFikggTw6a7b/462OxHoH6xCAkvHgjlBKSes4RSSQR2dj8qvx7IbEBulHdpaydNzpGnKUBHKKSCCPQmNLLAKsBoq6hnR425AMf3OSIS1kp8gJvRSRQB69N4FCgHW91lqHa8lpprSQ005rvQ1Y4XQd2UICefR29CO/sDHR7EQB4bjm6D8OMOfWAWb9foAfPmMv9Xfpw8PMuXWAQ28Z4Lx/DDEQ3XOPyvXPRzjwxn7qbhrg8dV2pnUOWhx3xyCzfz/Aouad58rm3zPE5v7MHMzgklEWaaeUeqxpQZN01aWIBPLobcRePcR8dn1i02BU96e7AI8JSxb4WP4lP8u+6OOxNXFe2RTnf8/IZ/mX/Pz3Mj/Tig1uei26275vdya4Z0WMFZf7eewiL5c/OkzC0vz9rRgL5rh5+VIfv3zJ3u+fLTHmVplUF2bm28YlXRZO+LvTBWSTzPyfNY4kp+JcBkzQwKpuK+2tZKUU/jx7fvaYBbGEfQlhkce+T2vNcEyzpxncH26O8z+z3HhcihmlBgeWGbzWlsBtKIbjmkhCYyiIW5obXo3yrWPz0vfC9pPbkpOq6aS17gGecLqObCKBnBovA/kAL7Um3naigISlOezWASp/2c9pM100TLGndrjk4WGqfj1A8zaLrzbsHqZt/RZTi3dG9ZRCg7Z+zYWHuHl8TYIz7h7iRyd4+P3rUT5zqBuvO3MX5siTLot0e6BpQZOM/U4hCeTUaMGeaMh4am18Qziuh9JdgGkoln3Jz6arCnltc4K3OuxP73fOL2DzVX4C5Qb3vrX7/509tSkVUJyvWHyhlze+4GfuJJNHVsb5+MFuPt84zHn/GOLl1sw7f+aRYYdppZS6x+kaso0Ecgo0tsQGgbeAsriFXrUt/d0W25XkK06c7uKx1TsD0zQUn5jl5oHg7iE6pcigtW9njm3qt6gufGcr+NpnI1x9vIe/N8U4otrkjvkFfG9JZOxexAh5pMsibbTWW4FnnK4j20ggp86LgBfgiTXxpel84s5Bi96wnUXDMc1T6+LUTTBY3W1/gtda88+VcerLd/91z6tzcc+KGJG4Zl2PxaptFkdP3jlx2qptCTYPWJxQ42IoZvcnKyCceQ1k8i1LAjlNlFL3NS1oki6iFJM5ZFOnmWS3xbMbEpsuOVx3lBWoynQ88ZYBzYJFQyQssDRcMMvNWQe5OP7OIUIRjdYwp8rglrMKAGhsifHG5gTXnpTPrEqTCw52c/DvB3AZipvPzMc0draQr14S4bqTPQB88hA359wzzG9fjXLtiZ50vLT9UiAN5HSS7ooxoKTbLXXm1bkvBw4Btl5+VF7DGQe6znC6plzytbziDc9MLp7udB3ZTlt6ozJUjYw/Tj3pskitf5McbXH/27HlcblyLK28lnyCTgdlqN9KGI8N6bJIrRagF/B2DOqhNd1WsK7cPMTponKFV+u0jMmLbovS9sc24n1xUFB6Yinlp5ez9YGthJaGUEphFplM+dwU3KW7z5sf3Ral7Y4g7QYqAAATSklEQVQ24t32/tOvnE5eRR6tt7YS3hSm8LBCqs6rAqDj4Q7yp+ZTNLcoHS/tfWlLDypD3e50HdlKAjmFGltiiXl17seB84Ghx1bH35BATh9fmgJZmYqq/6mioKaAxHCCNT9ag3+Wn/Izy5n48YkAbHtyGx0PdzD54sm77b/pD5uo/Ggl/tl+EuEESinCrfbl7rU/rWXt/1tLYiiBFbUYXjtM5fy0nIrYN5o/Ni1oCjldRraSLovUew17IILx9LrExvYBq9XpgnKFT1tpeT+7S9wU1NgnSM0CE0+1h3hPHLNg5+gUK2Kh1O5/H8JtYbDAP9tv759vYngMMEHHNNrS6LgGAzoe7KDyY5kTxlrrhDLVb5yuI5tJIKdYY0usG/tS6gqAf62Kv+BsRbnDT3payLuKdkYJbwhTcIAd0Fvv30rzVc30vtxL5bm7h2m0PYrpNdn4u42s/sFq2u9pR1ua/Op83GVu1vxwDcVHFRPdas8fUjC9IK2v573ohH6oaUGTNDDGkATy2HgUKABY1Bxf2T2sOxyuJyf4tU7rytOJcIKNN22k6sKqHa3jiedNpP439ZQcU8K2p7ftto+2NIMrB6n6RBUH/PAAop1Rep7vAWDSRZM48CcHUv6Rcrt1fG4lHY0dbLx5I93/7k7nS9sjw2Vc73QN2U4CeWysAVaSnHDoqbXSSk4HP1baAlnHNa03tVJyTAnFRxbv9njxB4oJvbF7V6u71E3+tHzyKvNQpqJwbiHhDeF3bBN6M0TBjAKsiEWkLcK0L0+j96VerIhzo0isuPVi04KmNx0rIEdIII+B5IrUD5OcuP6et2JvhSK6x9mqsl+6Wshaa9ruaMMzyUP5GeU77o+077ycvH9pP55Ju188UzCzAGvIIh6yR0QOBgfxVO/cTsc1257cRvlHyrGiFjum6NP2Y06R1nF6yCiLsRMEWoHiuEXfknXx58+pd89zuqhsVkh6Anlo1RC9L/XimeJh9TWrAburoue5HjuUFeRNyKP64moAhtcN0/1MN5M/OxllKKo+UcW6X6wD7D7i0hNLdxx729PbKDm2BMNjkD81HzSs+v4qCg8txPSltUdmBytqvbni8ysWO/LkOUau1BtD8+rcc4GvAesNhbpjfsGX0nU5dS7aoNXQ2TOnep2uI5torUFz9FuXvPW607XkAumyGFv/BbYCxZZGPxiMPel0QdmsUGv5xJdi1pD1iIRx+kggj6HGllgc+CtQZn8fX93aZ615773ESPmUBHIq6YSOGR7jMqfryCUSyGOvCXsmuAqAPy+PPWlJP9GY8ChlaFk1JGUSw4mb37r0rU1O15FLJJDHWHLExT2ADzBea0tsDXZayxwuK3tZyJJCKWDFrF6X3/V9p+vINRLIadDYElsHvARUAfzutehTkbgedraq7KRkhr2UsCLW95oWNA06XUeukUBOn4cAE3Bv7tdDj6+JP+V0QdlIWTrhdA3jXSKcWOXyu251uo5cJIGcJo0tsU5gEVAN8H9vxt7c0m9tcLaq7GNKC3lUtKUTOqb/R+Y7doYEcno9AXQAJRq48dVoo0xin1qmtJBHJdoV/V3wK0G5RNohWRPISqnjlVIrlFLLlFKZM0XWLhpbYhHgTqAUMFZ0Wt3PrEsscbisrGLKKIsRi/fFVw+vHf6m03XksnEVyMq2t5ovAn6ltT5M6/c/YfY+xxozjS2xZuBpYDLAza9HX5Gui9RxSyCPiBW3otHO6Lmtt7TKJwwHZXwgK6VqlFJBpdTvgTeBTyulXlZKvamUuk8p5VdKfQ64APiBUuqvyf2+pZR6XSn1X6XUj/dyrKlKqdPffbzktuuVUj9O3t+klKpP3u9XSt2ZvO+/SqmPJ+/f43H24gGgHyiyNPpnL0QeCMf10Bj9CHOKW7osRiTaHv3Jmp+secvpOnJdxgdyUh1wF3AacClwqtZ6LvAGcJXW+nagEfiW1voipdTpQC1wNHAYcIRS6kO7HktrfTgwCHz/3cfb5Xm7kvffAmz/KHcN0Ke1PkRrfSiwRClV/j7HeYfGltgg8EdgAmCu69X9f2uKLZLrRUbPbckPcX/FumNvbL1/63VO1yHGTyBv0Fq/AnwAOBh4USm1DFgA7GnZ99OTt6XYLeF67IDe9Vjsw/EeTH79D1CT/PepwM3bN9Ba9+xHXTs0tsRWAI8AUwEWNcdX/WeL9fJ77SPen0cCeb9YEWsw2hH9WGhpSH5uGWC8XPu/fYC6Ap7UWn/yfbZXwPVa69vecadSNbsca1+Ot32C2wQ7f1YKePebd1/rerdF2H8sqoGtP3sh8tQtZ+VPr/AZ1ft5HJGULx8z9pm2tA63hS9fe/1aWZYpQ4yXFvJ2rwDHKqUOBFBKeZVSB+1hu8eBz+7SHzxZqT1Oe7mvx9vVE8BXtn+jlCod4XFobInFgNuwfw/eaALrly9F7wvLVXwjlm/JOb19NbRm6A9rfrzmLqfrEDuNq0DWWncCFwN/V0r9FzsI6/ew3RPA34CXlVJNwP0kV+8YyfHe5adAqVLqLaXUcuCkER4HgMaWWAfwB+zLqo3mLqv3j/+J3puQk1MjUiBdFvtkeP3wK12Lu77mdB3inWSC+gwwr86tgAuxT1quB7jkMPdh5wbc852sazy6wl28fsmU4hqn68hk0Y7o+q2Lth7V+2Jvl9O1iHcaVy3kbJWcEe4fwAqS45PvXBZb9lpb/EVHCxuHvNpS779V7oqH4j29L/eeLWGcmSSQM0SyP/lWYBtQDnDdc9Gn1vZYzY4WNs54tZZA3gsrYoX73uj71NYHt65wuhaxZxLIGaSxJdYP3IA9K1yhBq5ZEn6wfcCSs+D7yCeBvEc6oROhN0NXb/7z5kedrkXsnQRyhmlsibVjh3IZ4OmPEvvuU5G/dg5aWxwubVzwaUve0++iLa37l/ff1vdq3w1O1yLem7x5M1BjS6wF+0q+yUDetmEd+d7Tkb90D+sOh0vLeH6t5T29C21pHXojdE/Pcz1XhpaGZExghpM3b4ZqbIm9BPwZmAK4tw7q4WuWhO/qDettDpeW0fxYptM1ZAqtte57re+R3pd6LwstDUWdrke8PwnkzLYEez2+qYCrNaQHf/hM+M+hiO5xuK6MJS1km7a07nut79G+V/ouDi0N9Tldj9g38ubNYMnhcP/CnlNjGsmJiL6/JHyndF/smR+d8y1kbWmr9+Xef/a93LcgtDTU7XQ9Yt9JIGe4ZCg/jD2b3XTAXN+r+7/9ZPhOGX2xu0Kd210WOqETPc/1LAq9Hro0tDQk3VvjjATyOJAM5Qexg3k64O4Y1OGrHg/ftb7XWulsdZmlUOnxMmFWyllxK9b9TPd9/cv6vxhaGnrfCz+UUonkCjsrlFLLlVJXvXvRBqXUb5VSbXu4/yNKqTeS84s3K6V+lerXk4skkMeJxpaYhR3Kd2P3KecPRIl/4/HwvW93JpY7W13m8OvcDOTEYCLU2dh5x8BbA5ftSxgnDSdX2JmFfdn+mcAPtz+YDOFzgVbgQ7vcPxu4CfiU1joAzAbWpuil5DSZy2IcmlfnPgb4ItAJDCrgW8fmnXDsVPNEpXL7uoiY1nruzOk59UOIdkTbOh7p+EsilLg+tDQU2tf9lFIDWmv/Lt/PBF4HyrXWWil1CvAN4F7gg1rrLya3uwv4t9b6jtS+EiEt5HGosSX2MvBr7ItHijXwixejz961PHZPNKEj7713dnMrpXQOzZQ3uHIwuOXeLT9PhBI/2Z8w3hOt9VrsTNg+Ve0ngb8DDwFnK6XcyftnYy/aIFJMAnmcamyJNQHXY0+cPxHggWC85afPRW7P+bHKlo47XcJY05ZO9L7U+2LXo13fI8HvQ0tDqVqTUQEopfKwuzAWaa1DwKvYq/CIMSSBPI41tsRWAz8G2rGHxall7VbX1x8L/zGXT/Ypi6wOZCtiDXYu7lzc91rfFaGloUWhpaGUfCJIdlkkgA7gDKAYaFJKrQeOw24xgz0r4RGpeE7xThLI41xjS6wL+BnwIjADyOse1pEr/hX++5J18adycaJ7I4tbyLFtsfb2f7TfPbxm+KuhpaGUdRsopSqwZxu8Sdsnlj4JfE5rXaO1rsF+b52ulPICvwS+t31VHKWUoZTa66K+Yt/JSb0skZzk/mTg00A3EAL4wBSz6rIj8z5WWqAqnKwvnQ6fWN0d97rKnK4jlXRCx/uX9y/tea7nAeDWVFx9p5RKAE2AG4gDfwF+A+QDm4CaZHfF9u0fBO7VWt+rlDob+9OZF3uNycVa62+NtqZcJ4GcZebVueuAywEf0AZorxvXt4/1nHp4ldGQC6MwjqyY1Bnxu7PmD1CsN7a5619dr0W3Ru8H/hFaGoo5XZMYGxLIWWhenbsQ+BRwDLAFGAY4s9Y189OHuuf78lSRk/WNtYYJVe1DRXlVTtcxWjqh4/1N/W/2PNvzKprbgabQ0pD8h81iEshZKtmF8QHsxVctYCtAaT55X2vwnHT4JONo411XX2WL40onbu4r8VQ7XcdoxPpiW7Y9tu21yJbIYuCe0NJQv9M1ibEngZzl5tW5K4DPYa+CvQUIAxxVbVR+/oi8s6v8xlQn6xsLJxVXbuoqy5/idB0joeM6OrBiYFn3s91vYPFHYLm0inOHBHIOmFfnNoETgE9gjzPdAmgFXHyY+7CP1LpOy3cpr5M1ptKHCytaN5cXjKs/NNrSVnhjeHn3ku7V8VB8CfC30V7oIcYfCeQcMq/OXQqchz2mtAfoBaj0qfwvHJF3/NxJxtEuQ437uSDm+co3rKv0Tne6jn0V2RoJdj/T/Xa0PboVuAN4U1rFuUkCOcck+5brgAXAJOzWcgTggFJVdOncvBMPrjDmjOf+5fO9EzY0T/RlfCDHemMbe57vWTa8ZrgLeBR4PLQ0NOB0XcI5Esg5al6dOw84CXs2Lzd2MMcBZlUYZZ+Z4z6hrtw4xBiH4+Q+nV+6ftmkwhqn69ibxGCiq+/1vjf6l/W3Ay8AD+/HDG0ii0kg57h5de4i7DkKzsAe4N+OffkssyqMsgtmuRtmVxqHu80dE8tkvC94Sja8XF2UcS3kWHdsbf/y/hX9/+3vQfNf4L7Q0tBGp+sSmUMCWQAwr849ATuUT2bnMLk4QLlX5V90iHtuwxSzwT8OxjBf4S5ev2RKcY3TdQDohI5FNkf+2/tqbzCyKWIBG7BnUGuWfmLxbhLI4h2Sw+ROB04ETOw5l4cB3AbG+bPcgRNrzIZMHi73XVfhhkemljraQk6EE33Da4bf6H2ptzUxmABYD/wTexhbzs0vIvaNBLLYo2RXxgeBswA/0EdyVAbAwRVG6UcPcs05dKI5p9CjShwqc49+avo33jutbFq6n1cndDzWHVsz2DL4dujNUC8WGngFeApYLy1i8X4kkMV7Sp78mwPMA6Zgd2N0ADGwBzWfMtOcdlKNa85BE4yDPS6V71ixSf9r+DbdMX1CWi4M0Qkdi3ZFVw2vHw72L+vfag1bXmAQeAx4MbQ01JOOOkR2kEAW+yQ5XG4G9uXYx2PPCDYEbMPuc8ZtYJw8wzXt6Mlmbe0E46CSfFXuRK1/VAVtN9ZUTB6r4+u4jka7oiuH1w6/3b+8v92KWIXYf5tagUeAZaGloehYPb/IXhLIYr8lW80B7Kv/DsMOo0Hsi0129I8Gyo3Sk2aYtYFy88Aqv5qartbz38jfcv2MykmpOp5O6Hi8P7451h1rDW8KbxhoGujSMV2YfHgL8Bz2NJZbpFtCjIYEshiVZF/zLOBo7LXWDOxQ7iY5bwbYiT13klFx+CRz6gGlxtTJRcbUYg8TxmKYcyN5W6+eUTVxpPtbUWsg1htrjXXGWsOt4dahNUPtOqb92CtoAGzEDuG3gA4JYZEqEsgiZebVuT3AAcCh2F0b2z/KR7EnzH/Hum+VPpV/SKVROaPUKK8uNMorvKqirECV+/MoGU1QL9HuritmTnrP7hKttbYiVsgatnoSg4mexECiJ9Yb6x5eP9wWbY8OAUXsnHwd7GXuXwRWhJaGOkdcnBDvQQJZjIlkn3MVMBV7prlZ2KsZW8lNhpK3MDtDDwB/Hq6ZpUbxRJ/yl3sNf2mB8pfk4y/MU36vWxWYBi6XgWkqXKaBaSrlMhRmf4KCBITfVq6uK6sr3TquIzquo1bMiuioDsf7473x3nhPtDPaE+2I9um4trCH9hVgB/D2eTwGgWbsFvBGoE36hEU6SCCLtElOnD8VqME+QTgNqMAOacXO1nQs+XX7bZ/epGtdrlPaXWZr1DB2XeBVAR7s0M0H8pLPp7G7V2LYF8G8DazEPjHXJd0QwgkSyMJR8+rcbmACUI7dgp6Y/L4MKMVuuSp2tqz3qtU0aze6XX0xw+hgZ8Br7P7sduwlrTZjj6fuSd4GJXxFppBAFhltXp3bwG7derG7FMw9fLWARL+hjGUez1BCqTB2yzoCRENLQ+8b5kJkAglkIYTIEON2zlshhMg2EshCCJEhJJCFECJDSCALIUSGkEAWQogMIYEshBAZQgJZCCEyhASyEEJkCAlkIYTIEBLIQgiRISSQhRAiQ0ggCyFEhpBAFkKIDCGBLIQQGUICWQghMoQEshBCZAgJZCGEyBD/H+jpGfC4uKGKAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "In_t = np.sqrt(In_sensor**2 + In_ref**2 + In_dac**2 + In_opamp_T**2)\n",
    "\n",
    "labels = 'sensor', 'reference', 'DAC', 'OpAmp'\n",
    "sizes = np.array([In_sensor, In_ref, In_dac, In_opamp_T])\n",
    "sizes /= sum(sizes)\n",
    "\n",
    "fig1, ax1 = plt.subplots()\n",
    "ax1.pie(sizes, labels=labels, autopct='%1.1f%%',\n",
    "        shadow=True, startangle=90)\n",
    "ax1.axis('equal')  # Equal aspect ratio ensures that pie is drawn as a circle.\n",
    "\n",
    "print_low_freq_noise(\"Total \", In_t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Long-term/Thermal stability"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## DAC\n",
    "- offset (zero-code) temp co is: +-0.05 ppm/C\n",
    "- gain temp co is +-0.1 ppm/C\n",
    "- negligible\n",
    "\n",
    "## OpAmps\n",
    "- Offset voltage: 0.3uV/C typical, assuming 2V as the sensor volage, this corresponds to approx 0.15ppm/C\n",
    "- offset/bias current: 20nA typical, corresponding to 20uV total offsets.\n",
    "  Assuming this is stable to even 1%/C, it's negligible\n",
    "\n",
    "## sense resistor\n",
    "- ± 0.05 ppm/°C\n",
    "- negligible\n",
    "\n",
    "## InAmp resistors:\n",
    "- LT5400-A-4 0.2ppm/C\n",
    "- negligible\n",
    "\n",
    "## fluxgate sensor:\n",
    "- 1 ppm/K (MAX, probably pessimistic)\n",
    "\n",
    "## reference:\n",
    "- 20ppm/rt(k Hour)\n",
    "- 1ppm/C typical assuming B-grade (else 2.5ppm/C)\n",
    "\n",
    "## subtractor resistors:\n",
    " - no temp co specified. Assume these are 5ppm/K\n",
    " - $V_0 = V_{\\mathrm{sensor}} + \\frac{R_{\\mathrm{f}}}{R_{\\mathrm{g}}}\\left(V_{\\mathrm{ref}}-V_{\\mathrm{sensor}}\\right) $ \n",
    " - for closed-loop operation, $V_0=0$ so $V_{\\mathrm{sensor}}\\left(1- \\frac{R_{\\mathrm{g}}}{R_{\\mathrm{f}}}\\right)= V_{\\mathrm{ref}}$\n",
    " - so, all temperature coefficients are suppressed by a factor of order $\\frac{R_{\\mathrm{f}}}{R_{\\mathrm{g}}} = 100 $\n",
    " - so, we do want something better than standard 100ppm/K, but don't need to be ultra-low temp co. (25ppm/K or lower is fine)\n",
    " - essentially, this is non-critical since the gain for the sensor and reference is almost identical, so drifts in the resistor values common-mode out"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Digital control loop\n",
    "\n",
    "In one experiment we measured the power supply noise to be ~1mA. This is very dependent on the PSU and coil inductance, but it gives a ball-park to aim for in terms of dynamic range."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total circuit gain: 212.1 V/A (1.0 V/uT)\n",
      "1 LSB at the ADC is: 1.4 uA (0.30 nT)\n",
      "ADC dynamic range is: 94.3 mA (19.8 uT)\n"
     ]
    }
   ],
   "source": [
    "# Feedback ADC\n",
    "\n",
    "G_error_amp = 1 + 20e3/1e3\n",
    "G_AFE = G_sense*G_subtractor*G_error_amp # Volts per amp of coil current measured after the error amplifier\n",
    "G_InAmp = 10**1  # InAmp gain\n",
    "G_T = G_AFE*G_InAmp  # Total gain (V/A) at the ADC input\n",
    "ADC_res = 20/(2**16)  # ADC resolution (V)\n",
    "\n",
    "print(\"Total circuit gain: {:.1f} V/A ({:.1f} V/uT)\".format(G_T, G_T/G_coil*1e-6))\n",
    "print(\"1 LSB at the ADC is: {:.1f} uA ({:.2f} nT)\".format(ADC_res/G_T*1e6, 1e9*ADC_res/G_T*G_coil))\n",
    "print(\"ADC dynamic range is: {:.1f} mA ({:.1f} uT)\".format(20./G_T*1e3, 1e6*20./G_T*G_coil))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 LSB at the DAC is: 305.2 nA (0.06 nT)\n",
      "DAC dynamic range is: 20.0 mA (4.20 uT)\n"
     ]
    }
   ],
   "source": [
    "# Feedback Iset DAC\n",
    "DAC_res = 20/2**16  # 1 DAC LSB in V\n",
    "G_shunt = 1e-3 # gain of current shunt (A/V)\n",
    "\n",
    "print(\"1 LSB at the DAC is: {:.1f} nA ({:.2f} nT)\".format(DAC_res*G_shunt*1e9, DAC_res*G_shunt*G_coil*1e9))\n",
    "print(\"DAC dynamic range is: {:.1f} mA ({:.2f} uT)\".format(20.*G_shunt*1e3, 20.*G_shunt*G_coil*1e6))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:artiq-env]",
   "language": "python",
   "name": "conda-env-artiq-env-py"
  },
  "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.5.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}