{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "791727d9",
   "metadata": {},
   "source": [
    "# Calculate the XRD Pattern of Co3O4 \n",
    "COID = 1526734"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "id": "5b6b0d93",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "lattice infos \n",
      "8.065 0.0 4.938388217561701e-16\n",
      "-4.938388217561701e-16 8.065 4.938388217561701e-16\n",
      "0.0 0.0 8.065\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAysAAAGbCAYAAADEAg8AAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABJf0lEQVR4nO3de3RU5b3/8c+ehFxNQshNrcSqS7lYkZtWEAGldtULbam1So8X1vLUarXVtnjs5fSnFK1aldpqtUfXUWldtYItaqttXd6hiBS5iEUo9qiBqiRDQiaBTIbMfn5/YMZcIROS7O9O3q+1WIvszCTPk3cm4cvO7HjOOScAAAAAMCYS9AIAAAAAoCsMKwAAAABMYlgBAAAAYBLDCgAAAACTGFYAAAAAmMSwAgAAAMAkhhUAAAAAJjGsAAAAADCJYQUAAACASQwrADAI/PjHP9bxxx+vE044QZMnT9Y777wjSTr99NP75O3feOONuueee/Z7m4cffljV1dVpvd2HH35Y8+fP73T8iiuu0F133ZV6+T//8z9T7/+Tn/ykxo0bp3HjxmnGjBl67733UrfLzMzUhAkTNHbsWE2aNEkPPPBAl+/3pZde0vDhwzV+/HiNGTNGCxYsSGtvd911lxKJRDpbBQD0AsMKAITcypUr9eKLL2r9+vXauHGjnnjiCQ0fPlyS9OKLLw7YOnozrHTn5ptv1i9+8Qt9+OGHWr16tdauXasrr7wy9fqVK1fqjTfe0GmnnaabbropdXz48OFat26dNm3apGXLlum+++7T/fff3+X7+MxnPqP169drzZo1euSRR7R27doe7603w4rv+2ndHgDAsAIAoffhhx+qtLRUw4YNkyQdccQRKi4uliSVlpZK2ncm4TOf+Yxmz56to446Sj/5yU/0q1/9ShMnTtSnP/1pRaNRSdLMmTP15ptvSpLefPNNzZw5s9P7+9WvfqWTTjpJJ554or761a9q7969WrZsmdasWaMvf/nLmjx5siRpzZo1mjFjhiZNmqTZs2ertrZWkvSnP/1Jxx13nCZPnqzly5d3uaeSkhL913/9l7773e/qm9/8pn7+858rIyOj0+2mTZum7du3d/k2Kisrdeedd+ree+/d78cvPz9fkyZN0ttvv92jvf3yl7/U+++/r6lTp+rzn/+8JOmvf/2rpkyZogkTJuiiiy5KDTIlJSW6+uqrdcIJJ+if//ynTjzxRF166aUaM2aMLrjgAjnn9rs2ABjqGFYAIOTOPPNMbd68WWPHjtU111yjNWvWdHm79evX68EHH9Sbb76pu+66S/F4XGvXrtUZZ5yh3/zmNz1+f1/5ylf097//XRs2bNChhx6qJUuWaM6cOZo8ebIef/xxrVmzRnv37tV3v/tdLVu2TK+//rrmzJmjW265RfF4XFdffbWef/55vfrqq9qyZUu37+fyyy/X66+/rqOPPlqnnXZal7d55plnUgNDVyZOnLjf9yFJO3fu1KpVq3T88cf3aG9XXXWVDj/8cK1cuVJPPfWUotGobr/9dr3wwgtat26djj766NSPn9XW1uqss87Sxo0blZOTo7feekvXX3+9Nm3apB07dmjFihU9+IgDwNCVGfQCAAAHp6CgQOvWrdOLL76o559/XmeeeaaWLFmiM888s93tpkyZorKyMkn7zr6cddZZkqQTTjhBr732Wo/f34YNG/SjH/1I9fX1qq+vV25ubqfbbNmyRRs2bNAZZ5whSWppadHxxx+vzZs367jjjtPIkSMl7Rt8qqqqunw/b7/9tvbs2aN//vOf8n1fkcjH/782depURaNR5ebmtvsxsI72d+biueee04QJExSJRPT9739fxx9/vF588cUD7q2jVatW6Y033tCUKVMkSc3NzTrnnHMkSbm5uam/S9KoUaM0duxYSdKECRP07rvvdjuIAQAYVgBgUMjMzNSZZ56pM888U6WlpXryySc7DSvZ2dmpv0cikdTLkUhEyWQy9XZan1vR3Nzc5fu67LLL9PTTT2vMmDG655579O6773a6je/7mjBhQqfnzKxfv16e5/VoT9/61rd0zz336A9/+IPuv/9+XXHFFanXrVy5UtnZ2Zo7d65uvPFGLVq0qMu3sX79eo0ePbrL133mM5/R448/nvbeOvJ9X+ecc44eeuihTq/Ly8tr93LbBhkZGamPOwCga/wYGACE3JYtW/Svf/1L0r4zCW+++aYqKyt79baOPPJIrV+/XpL0hz/8ocvb7N69WxUVFUokEnr00UdTxwsKCtTQ0CBJGj16tLZt26bXX39d0r7BZ/PmzRo9erT++c9/avv27WppadHSpUu7fB/Lli1TJBLR5z//ed1666366U9/qrq6una3GTZsmH7+859r8eLFnV4nSdu2bdP8+fN19dVX93j/Pdlbx5enTJmiF198MXVVslgslroaGwDg4DCsAEDINTY26qKLLtLxxx+vT33qU/J9X9/85jd79ba+853v6Kc//akmTZrU7dWubrzxRk2ePFnTp0/XuHHjUsfnzZunefPmafLkycrKytJjjz2ma665RieeeKImTZqkDRs2KCcnR7/4xS80a9YsTZkyRccdd1yntx+Px/W9730vdeniQw89VFdddZV+9KMfdbrtJz7xCc2dO1f33XefJGnXrl0aP368xo4dqy9+8Yu64oordNlll/V4/z3ZmyR97Wtf0+mnn67Pf/7zKisr0wMPPKDzzjtP48aN0/Tp09tdThkA0Hue41IkAAAAAAzizAoAAAAAkxhWAAAAAJjEsAIAAADAJIYVAAAAACYxrAAAAAAwiWEFAAAAgEmD5jfY+76v999/XwUFBT3+7cgAAAAABp5zTg0NDTr88MMViXR//mTQDCvvv/++Ro4cGfQyAAAAAPTQtm3bdMQRR3T7+kEzrBQUFEjat+HCwsKAV9MzzjklEgllZWVxNsgImthEF3toYhNd7KGJTXQJXiwW08iRI1P/hu/OoBlWWj/RCgsLQzOs+L6v6upqFRQU7Pf0FwYOTWyiiz00sYku9tDEJrrYcaBhkToAAAAATGJYAQAAAGASw0rAMjMHzU/iDRo0sYku9tDEJrrYQxOb6BIOnnPOBb2IvhCLxVRUVKT6+vrQPGcFAAAAGIp6+m93zqwEyDmnPXv2aJDMi4MCTWyiiz00sYku9tDEJrqEB8NKgJxzisViPFAMoYlNdLGHJjbRxR6a2ESX8GBYAQAAAGASwwoAAAAAkxhWAuR5Hr851Ria2EQXe2hiE13soYlNdAkPrgYGAAAAYECZuxpYfX29Tj75ZB1yyCF68803JUlLly7V1KlTNWvWLG3fvl2StHnzZk2fPl1Tp07V888/P1DLC4RzTg0NDTy5yxCa2EQXe2hiE13soYlNdAmPARtW8vLy9PTTT+vLX/6yJKmlpUWLFi3SSy+9pB//+MdauHChJOkHP/iB/vd//1d/+ctf9P/+3/8bqOUFwjmn3bt380AxhCY20cUemthEF3toYhNdwmPAfnXnsGHDVFZWlnp569atGjNmjLKysnTqqadq/vz5kqT3339fxx57rCRpxIgRikajKi0t7fT2mpub1dzcnHo5FotJknzfl+/7kvb9PKLneXLOtftkPNDx1vv39ngkEun0trs67vt+6u/prtHqnvZ3PAx7ar2Nc65Htw/DntI93t0at23bpp07d6qkpEQjR44MZE9tXzfYPvfCuKe2TQbLnnq7dit7au3S9vtM2Pd0oOPW99T2e/1g2VNvjlvbk9T5e33Y9xS2Th3fd3cGbFjpqK6urt3PpyWTSUntH8hFRUWqra3tcli55ZZbtGDBgk7Ha2pqFI/HJUm5ubkqKipSLBZTU1NT6jb5+fkqKChQXV2dEolE6nhhYaHy8vJUW1urlpaW1PHi4mJlZ2erpqam3Qe9pKREGRkZqq6ubreG8vJyJZNJ7dy5M3XM8zxVVFQokUiorq4utdfGxkZVVFSoqakpNXBJUlZWlkaMGKHGxkbt3r07ddz6niQpMzNTpaWlodyTc06e5ymRSKi+vn5Q7KkvOtXU1OjLXz5f//znFh133Cg9/vjS1H8+DNSedu3aJeecIpHIoPzcC+Oe6uvrU00Gy57C3sn3fdXX16ukpESe5w2KPYW9k+/7qX/jDJY9SeHvlJOTk/oxsEgkMij2FLZODQ0N6okBf4L9vHnzNH/+fEUiEd1xxx168MEHJUlTpkzRq6++qlNOOUWrVq2SJM2ePVsPPfRQj8+sjBw5st0QZH0qdW7fLyQaPnx46uWertHqnvZ3PAx7cs6psbFRBQUF6iise0r3eFdrXLdunU4++WT5/jRFIiu0evVqTZgwYcD21PoPsMLCwtT/iA22z72w7SmZTCoWi6WaDIY9DYZOrd9XioqKUo+dsO/pQMet78m5fc+NKCoq6rSWsO6pN8et7UmSdu3a1e77Stj3FLZOsVhMxcXFB3yCfWBnVo499li99dZbSiQSWrNmjcaNGydJOuyww/Svf/1L5eXl3Z5VkaTs7GxlZ2d3Oh6JRFITcqvWD05H3R3veP/eHO/p+ywuLm73ut6+nQMdH8g9DdTx/tpTUVFRl/fvzRqt7Olgj3/8Ba1Ivu/L87weP876Yk+RSKTdY6Uv9nQwx612Opjj6e4pIyOjU5Ow72mwdGrbZbDs6UDHre+p9T8lu7ptV7cPcu1DqVNX31e6W3t3x63tKUydultTRwM6rJx99tlav369tmzZoq9//eu69tprNXPmTOXk5Gjx4sWSpJtvvlnz5s1TMpns8se8BpPWqbLtVI9g0cQmuthDE5voYg9NbKJLeAzosPLMM890OnbBBRe0e3ns2LFavnz5QC0pUM45NTU1qaCggAeKETSxiS720MQmuthDE5voEh78BnsAAAAAJjGsAAAAADCJYSVAnucpPz+f04+G0MQmuthDE5voYg9NbKJLeAR2NTDse6B0dYlcBIcmNtHFHprYRBd7aGITXcKDMysBcs6ptra20zWrERya2EQXe2hiE13soYlNdAkPhpUAOeeUSCR4oBhCE5voYg9NbKKLPTSxiS7hwbACAAAAwCSGFQAAAAAmMawEyPM8fnOqMTSxiS720MQmuthDE5voEh5cDSxAnucpLy8v6GWgDZrYRBd7aGITXeyhiU10CQ/OrATI931Fo1H5vh/0UvARmthEF3toYhNd7KGJTXQJD4aVgLW0tAS9BHRAE5voYg9NbKKLPTSxiS7hwLACAAAAwCSGFQAAAAAmMawEyPM8FRcXcyUKQ2hiE13soYlNdLGHJjbRJTy4GliAPM9TdnZ20MtAGzSxiS720MQmuthDE5voEh6cWQmQ7/vasWMHV6IwhCY20cUemthEF3toYhNdwoNhJWDOuaCXgA5oYhNd7KGJTXSxhyY20SUcGFYAAAAAmMSwAgAAAMAkhpUAeZ6nkpISrkRhCE1soos9NLGJLvbQxCa6hAfDSoA8z1NGRgYPFENoYhNd7KGJTXSxhyY20SU8GFYC5Pu+qquruRKFITSxiS720MQmuthDE5voEh4MKwAAAABMYlgBAAAAYBLDCgAAAACTGFYCFIlEVF5erkiEDFbQxCa62EMTm+hiD01sokt4UChAzjklk0l+g6ohNLGJLvbQxCa62EMTm+gSHgwrAXLOaefOnTxQDKGJTXSxhyY20cUemthEl/BgWAEAAABgEsMKAAAAAJMYVgLGb061hyY20cUemthEF3toYhNdwiEz6AUMZZFIRBUVFUEvA23QxCa62EMTm+hiD01sokt4cGYlQM45NTc38+QuQ2hiE13soYlNdLGHJjbRJTwYVgLknFNdXR0PFENoYhNd7KGJTXSxhyY20SU8GFYAAAAAmMSwAgAAAMAkhpWAZWZyjQNraGITXeyhiU10sYcmNtElHKgUoEgkotLS0qCXgTZoYhNd7KGJTXSxhyY20SU8OLMSIOec9uzZw5O7DKGJTXSxhyY20cUemthEl/BgWAmQc06xWIwHiiE0sYku9tDEJrrYQxOb6BIeDCsAAAAATGJYAQAAAGASw0qAPM9TVlaWPM8Lein4CE1soos9NLGJLvbQxCa6hAdXAwuQ53kaMWJE0MtAGzSxiS720MQmuthDE5voEh6cWQmQc04NDQ08ucsQmthEF3toYhNd7KGJTXQJD4aVADnntHv3bh4ohtDEJrrYQxOb6GIPTWyiS3gwrAAAAAAwiWEFAAAAgEkMKwHyPE+5ublcicIQmthEF3toYhNd7KGJTXQJD64GFiDP81RUVBT0MtAGTWyiiz00sYku9tDEJrqEB2dWAuScU319PU/uMoQmNtHFHprYRBd7aGITXcKDYSVAzjk1NTXxQDGEJjbRxR6a2EQXe2hiE13Cg2EFAAAAgEkMKwAAAABMYlgJkOd5ys/P50oUhtDEJrrYQxOb6GIPTWyiS3hwNbAAeZ6ngoKCoJeBNmhiE13soYlNdLGHJjbRJTwCPbPi+77mzZun0047TdOmTdPmzZu1YsUKTZ06VdOmTdPGjRuDXF6/c86ptraWJ3cZQhOb6GIPTWyiiz00sYku4RHomZX169erublZy5cv1/Lly7Vo0SJt2bJFTz/9tBoaGnTFFVfomWeeCXKJ/co5p0QiIeccpyGNoIlNdLGHJjbRxR6a2ESX8Ah0WDniiCPknJNzTnV1dcrPz1dGRoaKi4tVXFys2trabu/b3Nys5ubm1MuxWEzSvrM1vu9L2neKz/O81PtodaDjrffv7fFIJNLpbXd13Pf91N/TXaPVPe3veBj21Hob51yPbh+GPaV7vKs1OucUiUTk+x/fp6ePs77cU9vXDbbPvTDuqW2TwbKn3q7dyp5au7T9PhP2PR3ouPU9tf1eP1j21Jvj1vYkdf5eH/Y9ha1Tx/fdnUCHldLSUg0bNkyjR49WPB7X8uXL9a1vfSv1+szMTCUSCWVlZXW67y233KIFCxZ0Ol5TU6N4PC5Jys3NVVFRkWKxmJqamlK3yc/PV0FBgerq6pRIJFLHCwsLlZeXp9raWrW0tKSOFxcXKzs7WzU1Ne0+6CUlJcrIyFB1dXW7NZSXlyuZTGrnzp2pY57nqaKiQolEQnV1dZL2fdFqbGxURUWFmpqaUgOXJGVlZWnEiBFqbGzU7t27U8et70na1620tDSUe3Ju3/+wJBIJ1dfXD4o99UWneDyuyspKvfuuVFlZqXg8nnrfA7WnXbt2ybl9Q9Ng/NwL455af6FaJBIZNHsKeyff91VfX6+SkhJ5njco9hT2Tr7vK5lMStKg2ZMU/k45OTlqaGhIfQ0bDHsKW6eGhgb1hOc6jkgD6C9/+Ysee+wxPfTQQ1qzZo1uu+027dy5Uy+88IIk6ZRTTtGqVau6vG9XZ1ZGjhypuro6FRYWSrI/lTq37xcS5efnp17u6Rqt7ml/x8OwJ+ecmpublZOTo47Cuqd0j3e1xnXr1unkk0+W75+jSORprV69WhMmTBiwPfm+rz179ig3Nzf1P2KD7XMvbHtKJpNqampKNRkMexoMnVq/r+Tl5aUeO2Hf04GOW9+Tc07xeFx5eXmd1hLWPfXmuLU9SdLu3bvbfV8J+57C1ikWi6m4uFj19fWpf7t3JdAzK845lZSUSNp3lqWhoUEtLS3atWuXGhoaNGLEiG7vm52drezs7E7HI5FIakJu1frB6ai74x3v35vjPX2fhxxySLvX9fbtHOj4QO5poI73157y8vK6vH9v1mhlTwd7vO0XNN/35Xlejx9nfbGnSCTS7rGSztr747jVTgdzPN09ZWRkdGoS9j0Nlk5tuwyWPR3ouPU9tf6nZFe37er2Qa59KHXq6vtKd2vv7ri1PYWpU3dr6ijQYeXMM8/Uww8/rBkzZqi5uVmLFi1SS0uLzj77bHmep3vvvTfI5fU73/dVW1urESNG9DgY+hdNbKKLPTSxiS720MQmuoRHoMNKZmamHnvssU7HV65cGcBqgtH25wphA01soos9NLGJLvbQxCa6hAOjJAAAAACTGFYAAAAAmMSwEiDP81RcXNztE+4w8GhiE13soYlNdLGHJjbRJTwCfc7KUOd5XpdXNENwaGITXeyhiU10sYcmNtElPDizEiDf97Vjx44e/wZP9D+a2EQXe2hiE13soYlNdAkPhpWAdfzlOggeTWyiiz00sYku9tDEJrqEA8MKAAAAAJMYVgAAAACYxLASIM/zVFJSwpUoDKGJTXSxhyY20cUemthEl/BgWAmQ53nKyMjggWIITWyiiz00sYku9tDEJrqEB8NKgHzfV3V1NVeiMIQmNtHFHprYRBd7aGITXcKDYQUAAACASQwrAAAAAExiWAEAAABgEsNKgCKRiMrLyxWJkMEKmthEF3toYhNd7KGJTXQJDwoFyDmnZDLJb1A1hCY20cUemthEF3toYhNdwoNhJUDOOe3cuZMHiiE0sYku9tDEJrrYQxOb6BIeDCsAAAAATGJYAQAAAGASw0rA+M2p9tDEJrrYQxOb6GIPTWyiSzhkBr2AoSwSiaiioiLoZaANmthEF3toYhNd7KGJTXQJD86sBMg5p+bmZp7cZQhNbKKLPTSxiS720MQmuoQHw0qAnHOqq6vjgWIITWyiiz00sYku9tDEJrqEB8MKAAAAAJMYVgAAAACYxLASsMxMrnFgDU1soos9NLGJLvbQxCa6hAOVAhSJRFRaWhr0MtAGTWyiiz00sYku9tDEJrqEB2dWAuSc0549e3hylyE0sYku9tDEJrrYQxOb6BIeDCsBcs4pFovxQDGEJjbRxR6a2EQXe2hiE13Cg2EFAAAAgEkMKwAAAABMYlgJkOd5ysrKkud5QS8FH6GJTXSxhyY20cUemthEl/DgamAB8jxPI0aMCHoZaIMmNtHFHprYRBd7aGITXcKDMysBcs6poaGBJ3cZQhOb6GIPTWyiiz00sYku4cGwEiDnnHbv3s0DxRCa2EQXe2hiE13soYlNdAkPhhUAAAAAJjGsAAAAADCJYSVAnucpNzeXK1EYQhOb6GIPTWyiiz00sYku4cHVwALkeZ6KioqCXgbaoIlNdLGHJjbRxR6a2ESX8ODMSoCcc6qvr+fJXYbQxCa62EMTm+hiD01sokt4MKwEyDmnpqYmHiiG0MQmuthDE5voYg9NbKJLeDCsAAAAADCJYQUAAACASQwrAfI8T/n5+VyJwhCa2EQXe2hiE13soYlNdAkPrgYWIM/zVFBQEPQy0AZNbKKLPTSxiS720MQmuoQHZ1YC5JxTbW0tT+4yhCY20cUemthEF3toYhNdwoNhJUDOOSUSCR4ohtDEJrrYQxOb6GIPTWyiS3gwrAAAAAAwiWEFAAAAgEkMKwHyPE+FhYVcicIQmthEF3toYhNd7KGJTXQJD64GFiDP85SXlxf0MtAGTWyiiz00sYku9tDEJrqER4/PrLz22mu65pprdNJJJ2nkyJE67rjj9IUvfEG/+tWvVF9f359rHLR831c0GpXv+0EvBR+hiU10sYcmNtHFHprYRJfw6NGwcu655+qhhx7SZz/7WT311FN65513tHbtWi1YsEDNzc368pe/rKeeeqq/1zootbS0BL0EdEATm+hiD01soos9NLGJLuHQox8De+SRRzR8+PB2xw455BCNHz9e48eP1zXXXKNdu3b1w/IAAAAADFU9OrPScVDp7W0AAAAAoKcO+mpgEydO7It1DEme56m4uJgrURhCE5voYg9NbKKLPTSxiS7hcdBXA1u7dm1frGNI8jxP2dnZQS8DbdDEJrrYQxOb6GIPTWyiS3ikdWalqqqqyz8H46WXXtKsWbN0+umna9myZVqxYoWmTp2qadOmaePGjQf1tq3zfV87duzgShSG0MQmuthDE5voYg9NbKJLeKR1ZuW8886T53lyzqm5uVlbtmzR2LFjtW7dul6986amJt15553685//rKysLEnSjBkz9PTTT6uhoUFXXHGFnnnmmV697bBwzgW9BHRAE5voYg9NbKKLPTSxiS7hkNaw8ve//73dy5s2bdKtt97a63f+6quvKjc3V7Nnz1ZeXp7uvfdeZWRkqLi4WMXFxaqtre32vs3NzWpubk69HIvFJO2blFunZM/zUsNV20/IAx3vOGWnezwSiXR6210d930/9fd012h1T/s7HoY9td7GOdej24dhT+ke72qNzjlFIhH5/sf36enjrC/31PZ1g+1zL4x7attksOypt2u3sqfWLm2/z4R9Twc6bn1Pbb/XD5Y99ea4tT1Jnb/Xh31PYevU07NaB/WclbFjx+r111/v9f137Niht99+W6tWrdJzzz2nG264QYWFhR8vLjNTiUQiddalrVtuuUULFizodLympkbxeFySlJubq6KiIsViMTU1NaVuk5+fr4KCAtXV1SmRSKSOFxYWKi8vT7W1te2uvV1cXKzs7GzV1NS0+6CXlJQoIyND1dXV7dZQXl6uZDKpnTt3po55nqeKigolEgnV1dVJ2vdFq7GxURUVFWpqakoNXJKUlZWlESNGqLGxUbt3704dt74naV+30tLSUO7JOSfP85RIJNr9stMw76kvOsXjcVVWVurdd6XKykrF4/HU+x6oPe3atUvO7RuaBuPnXhj3VF9fn2oyWPYU9k6+76u+vl4lJSXyPG9Q7CnsnXzfVzKZlKRBsycp/J1ycnLU0NCQ+ho2GPYUtk4NDQ3qCc91HJH247rrrktNo77va8OGDSopKdHvfve7nr6Jdv785z/rmWee0d13363m5mZNmzZNBQUFeuGFFyRJp5xyilatWtXlfbs6szJy5EjV1dWlBh7rU6lzTi0tLalhjEk7+D21frPPyMhQR2HdU7rHu1rjunXrdPLJJ8v3z1Ek8rRWr16tCRMmDNiefN/X3r17lZmZmfoaNNg+98K2p2QyqZaWllSTwbCnwdCp9fvKsGHDUo+dsO/pQMet78k5p2QyqWHDhnVaS1j31Jvj1vYkSYlEot33lbDvKWydYrGYiouLVV9f3+5kRUdpnVn51Kc+9fEdMzP1xS9+UdOmTUvnTbRz0kkn6c4775RzTuvXr9fYsWP1zjvvaNeuXWpoaNCIESO6vW92dnaXV3GIRCKpCblV6weno+6Od7x/b4739H22fTmdNVre00Ad7689tf7Dqyth3dPBHm/7Bc33/S4/Rv25p0gkkvrHV1/t6WCOW+10MMfT3VNGRkanz4Ow72mwdGrbZbDs6UDHre+p9XZd3bar2we59qHUqavvK92tvbvj1vYUpk7dramjtIaVSy+9NJ2bH1BpaanmzJmjGTNmyPM8Pfjgg/r3v/+ts88+W57n6d577+3T92eN7/uqrq5WeXl5j4Ohf9HEJrrYQxOb6GIPTWyiS3gc9O9Zufzyy3X//ff3+v5XXXWVrrrqqtTLxxxzjFauXHmwywIAAAAQcgc9Sn7961/vi3UAAAAAQDsHPaxMmjSpL9YBAAAAAO2kNay88cYb+uxnP6uRI0eqvLw89Qe903q5T35W0g6a2EQXe2hiE13soYlNdAmPtJ6zcvnll+uBBx7QJZdcohUrVuihhx5qd71lpKf1cob7u0IIBhZNbKKLPTSxiS720MQmuoRHWuNkS0uLTjjhBCWTSeXn5+vqq6/WsmXL+mttg55zTjt37ux0zWoEhyY20cUemthEF3toYhNdwiOtMyt5eXnau3evxo0bpx/96Ec6/PDDU7+VFQAAAAD6UlpnVh5++GElk0ndc889ysjI0NatW/X73/++v9YGAAAAYAhL68zK0UcfLUnKycnRjTfe2B/rGXL4OUl7aGITXeyhiU10sYcmNtElHHp0ZuX888/Xc889J9/3O73u/fff109/+lMtXry4zxc32EUiEVVUVHAlCkNoYhNd7KGJTXSxhyY20SU8enRm5Wc/+5nuvPNOXXnllTr00ENVVlameDyu//u//9Nhhx2mK6+8Ul/5ylf6e62DjnNOiURCWVlZTPdG0MQmuthDE5voYg9NbKJLePRoWDniiCP0s5/9TD/72c/0zjvv6MMPP1ROTo6OPfZYHXLIIf29xkHLOae6ujqVl5fzQDGCJjbRxR6a2EQXe2hiE13CI63nrEjSUUcdpaOOOqo/1gIAAAAAKfygHgAAAACTGFYClpmZ9skt9DOa2EQXe2hiE13soYlNdAmHtIaVZ555pssrgqF3IpGISktLuRKFITSxiS720MQmuthDE5voEh5pFfrd736nUaNGaf78+dq4cWN/rWnIcM5pz549cs4FvRR8hCY20cUemthEF3toYhNdwiOtYeXXv/611q9fr3Hjxunb3/62Jk+erLvuukvRaLS/1jeoOecUi8V4oBhCE5voYg9NbKKLPTSxiS7hkfa5r/z8fJ166qmaOnWq6uvrtXr1ap1++um6+eab+2N9AAAAAIaotJ5ZdP/99+s3v/mNhg0bpnnz5mnDhg3Ky8tTS0uLjjvuOP3whz/sr3UCAAAAGGLSGla2b9+uxYsX6+ijj27/RjIztWzZsj5d2FDgeR6/OdUYmthEF3toYhNd7KGJTXQJj7R+DCwajXYaVL7xjW9Ikk488cS+W9UQ4XmeRowYwQPFEJrYRBd7aGITXeyhiU10CY+0hpVVq1a1e9n3fb388st9uqChxDmnhoYGntxlCE1soos9NLGJLvbQxCa6hEePhpXbb79dZWVl2rhxo8rLy1VeXq6ysjJ94hOf0OzZs/t7jYOWc067d+/mgWIITWyiiz00sYku9tDEJrqER4+Gleuuu041NTW6/vrrVV1drerqatXU1OiDDz7Qrbfe2t9rBAAAADAE9egJ9u+8846OOuooffWrX9WmTZs6vX7s2LF9vjAAAAAAQ1uPhpVbbrlF999/v6666qpOr/M8Ty+88EKfL2wo8DxPubm5PLnLEJrYRBd7aGITXeyhiU10CY8eDSv333+/JOnFF1/s18UMNZ7nqaioKOhloA2a2EQXe2hiE13soYlNdAmPtK4G9pOf/EQNDQ1KJpM6//zzddxxx+mJJ57op6UNfs451dfX8+QuQ2hiE13soYlNdLGHJjbRJTzSGlaWLl2qgoIC/elPf1JOTo5WrFihBQsW9NfaBj3nnJqamnigGEITm+hiD01soos9NLGJLuGR1rDS1NQkSXryySd14YUXqry8nMgAAAAA+kWPnrPSavbs2TryyCM1fPhw/c///I9qamqUnZ3dX2sDAAAAMISldWbl9ttv17p167R27VoNGzZM+fn5evLJJ/trbYOe53nKz8/nShSG0MQmuthDE5voYg9NbKJLeKR1ZkWS3nrrLb377rtKJpOpY5dcckmfLmqo8DxPBQUFQS8DbdDEJrrYQxOb6GIPTWyiS3ikNazMnTtXH374oSZMmKCMjAxJYiI9CM451dXVqbi4mI+jETSxiS720MQmuthDE5voEh5pDSsbNmzo8jfYo3ecc0okEnLO8UAxgiY20cUemthEF3toYhNdwiOt56ycfPLJ2rJlS3+tBQAAAABS0jqzsn79ep144okaNWqUsrOzU9Po6tWr+2t9AAAAAIaotIYVrvzVtzzPU2FhIacfDaGJTXSxhyY20cUemthEl/BI68fAjjzySG3atElPPfWUjjzySGVnZ2vPnj39tbZBz/M85eXl8UAxhCY20cUemthEF3toYhNdwiOtYWX+/Pn63e9+p1/+8peSpIyMDM2bN68/1jUk+L6vaDQq3/eDXgo+QhOb6GIPTWyiiz00sYku4ZHWj4E9//zzWrdunSZMmCBJKisrUzwe75eFDRUtLS1BLwEd0MQmuthDE5voYg9NbKJLOKR1ZmXYsGHyfT91yqy2tlaRSFpvAgAAAAB6JK1J41vf+pYuuOACRaNRLVy4UNOnT9d1113XX2sDAAAAMISl9WNgF110kSZPnqznn39evu9ryZIlGjt2bH+tbdDzPI/fnGoMTWyiiz00sYku9tDEJrqER1rDyje+8Q3de++9Gj16dKdjSJ/necrOzg56GWiDJjbRxR6a2EQXe2hiE13CI60fA1u1alW7l33f18svv9ynCxpKfN/Xjh07uBKFITSxiS720MQmuthDE5voEh49GlZuv/12lZWVaePGjSovL1dZWZnKysr0iU98Queee25/r3FQc84FvQR0QBOb6GIPTWyiiz00sYku4dCjYeW6665TTU2Nrr/+elVXV6umpkY1NTX64IMPdNttt/X3GgEAAAAMQWk9Z+Wmm27S9u3b9d577ymZTKaOT58+vc8XBgAAAGBoS2tYmT9/vp566il96lOfUkZGhqR9T1BiWOkdz/NUUlLClSgMoYlNdLGHJjbRxR6a2ESX8EhrWPnjH/+of/zjHxo2bFh/rWdI8TxPGRkZPFAMoYlNdLGHJjbRxR6a2ESX8EjramBjx45VfX19f61lyPF9X9XV1VyJwhCa2EQXe2hiE13soYlNdAmPtM6s7N69W2PGjNHUqVPbXZt6yZIlfb4wAACAVlVVVYpGoyotLVVlZWXQywEwQNIaVn74wx/21zoAAAC6VFVVpVGjRiseb1JOTq62bNnMwAIMEWkNKzNmzOivdQAAAHQpGo0qHm+SNF3x+CuKRqMMK8AQ0aNhpaysrMsnIDnn5Hmeqqur+3xhQ0EkElF5ebkikbSeOoR+RBOb6GIPTWwa/F2Kgl5A2gZ/k3CiS3j0aFipqanp73UMSc45JZNJeZ7H1SiMoIlNdLGHJjbRxR6a2ESX8DAxTj766KMqKyuTJC1dulRTp07VrFmztH379oBX1r+cc9q5c6ecc0EvBR+hiU10sYcmNtHFHprYRJfwCHxYSSaTWrp0qUaOHKmWlhYtWrRIL730kn784x9r4cKFQS8PAAAAQEDSeoJ9f3j00Ud1/vnn684779TWrVs1ZswYZWVl6dRTT9X8+fO7vV9zc7Oam5tTL8diMUn7rpvdes3s1lN7zrl2k/OBjne85na6xyORSKe33dVx3/dTf093jVb3tL/jYdhT622ccz26fRj2lO7xrtbonFMkEpHvf3yfnj7O+nJPbV832D73wrintk0Gy556u3Yre2rt0vb7TNj3JKnd/Tp+DbK+p7bf6wfz517Y9iR1/l4f9j2FrVNPf8dNoMNKMpnUkiVL9MQTT+jOO+9UXV2dCgsL272+O7fccosWLFjQ6XhNTY3i8bgkKTc3V0VFRYrFYmpqakrdJj8/XwUFBaqrq1MikUgdLywsVF5enmpra9XS0pI6XlxcrOzsbNXU1LT7oJeUlCgjI6PTBQbKy8uVTCa1c+fO1DHP81RRUaFEIqG6ujpJ+75oNTY2qqKiQk1NTamBS5KysrI0YsQINTY2avfu3anj1vckSZmZmSotLQ3lnpzb94/yRCLR7heghnlPfdEpHo+rsrJS774rVVZWKh6Pp973QO2pvr4+1Wcwfu6FcU+xWCzVZLDsKeydfN9XLBZTSUmJPK/zBXDCuCdp39egnJwcJRKeJk2a1O5rkPU9tf1P1MH8uRe2PeXk5KihoSH1NWww7ClsnRoaGtQTnus4Ig2gxYsXKyMjQxdddJEmT56sX//617rjjjv04IMPSpKmTJmiV199tcv7dnVmZeTIke0GnqEwlbIn9jQQe1q3bp1OPvlk+f45ikSe1urVqzVhwoRQ7+lgjrMn9sSeBnZP69at0+TJkyXN7vQ1KKx72t9x9sSehsKeYrGYiouLVV9f3+5kRUeBnlnZtGmT1q1bp0ceeURbt27V3XffrbfeekuJREJr1qzRuHHjur1vdna2srOzOx2PRCKdLkPX+sHpqLvj3V3GLp3jPXmfzjklEgllZWWlvUarexrI4/2xJ+ecmpublZWVFcjHwGqntl/QfN+X53k9fpz1xZ4kae/evanHSl/s6WCOW+10MMfT3ZPnee2+fgW5djp9fLzt95Xubh+2PbXet1VXX4Ms76ltk8H8udeT45b25Jzr8vtKd2vv7rilPaW79u6OD9Tau1tTR4EOK7fddlvq75MnT9Z9992nxx57TDNnzlROTo4WL14c4Or6n3NOdXV1Ki8v7zImBh5NbKKLPTSxiS720MQmuoRH4E+wb7VmzRpJ0gUXXKALLrgg4NUAAAAACFrgly4GAAAAgK4wrAQsM9PMyS18hCY20cUemthEF3toYhNdwoFKAYpEIiotLQ16GWiDJjbRxR6aBK+qqkrRaFSlpaWqrKyURBeLaGITXcKDMysBcs5pz549nS4Dh+DQxCa62EOTYFVVVWnUqNGaNGmSRo0araqqKkl0sYgmNtElPBhWAtR6jWkeKHbQxCa62EOTYEWjUcXjTZKmKx5vUjQalUQXi2hiE13Cg2EFAIDQKgp6AQDQrxhWAAAAAJjEsBIgz/O6/M2pCA5NbKKLPTSxiS720MQmuoQHVwMLkOd5GjFiRNDLQBs0sYku9tDEJrrYQxOb6BIenFkJkHNODQ0NPLnLEJrYRBd7aGITXeyhiU10CQ+GlQA557R7924eKIbQxCa62EMTm+hiD01sokt4MKwAAAAAMIlhBQAAAIBJDCsB8jxPubm5XInCEJrYRBd7aGITXeyhiU10CQ+uBhYgz/NUVMQv9LKEJjbRxR6a2EQXe2hiE13CgzMrAXLOqb6+nid3GUITm+hiD01soos9NLGJLuHBsBIg55yampp4oBhCE5voYg9NbKKLPTSxiS7hwbACAAAAwCSGFQAAAAAmMawEyPM85efncyUKQ2hiE13soYlNdLGHJjbRJTy4GliAPM9TQUFB0MtAGzSxiS720MQmuthDE5voEh6cWQmQc061tbU8ucsQmthEF3toYhNd7KGJTXQJD4aVADnnlEgkeKAYQhOb6GIPTWyiiz00sYku4cGPgQF9oKqqStFoVKWlpaqsrAx6OQAAAIMCwwpwkKqqqjRq1GjF403KycnVli2bGVgAAAD6AD8GFiDP81RYWMiVKAzpTZNoNKp4vEnSdMXjTYpGo/23wCGKx4o9NLGJLvbQxCa6hAdnVgLkeZ7y8vKCXgbaOLgmRX26FnyMx4o9NLGJLvbQxCa6hAdnVgLk+76i0ah83w96KfgITWyiiz00sYku9tDEJrqEB8NKwFpaWoJeAjqgiU10sYcmNtHFHprYRJdwYFgBAAAAYBLDCgAAAACTGFYC5HmeiouLuRKFITSxiS720MQmuthDE5voEh5cDSxAnucpOzs76GWgDZrYRBd7aGITXeyhiU10CQ/OrATI933t2LGDK1EYQhOb6GIPTWyiiz00sYku4cGwEjDnXNBLQAc0sYku9tDEJrrYQxOb6BIODCsAAAAATGJYAQAAAGASw0qAPM9TSUkJV6IwhCY20cUemthEF3toYhNdwoNhJUCe5ykjI4MHiiE0sYku9tDEJrrYQxOb6BIeDCsB8n1f1dXVXInCEJrYRBd7aGITXeyhiU10CQ+GFQAAAAAmMawAAAAAMIlhBQAAAIBJDCsBikQiKi8vVyRCBitoYhNd7KGJTXSxhyY20SU8KBQg55ySySS/QdUQmthEF3toYhNd7KGJTXQJD4aVADnntHPnTh4ohtDEJrrYQxOb6GIPTWyiS3gwrAAAAAAwiWEFAAAAgEkMKwHjN6faQxOb6GIPTWyiiz00sYku4ZAZ9AKGskgkooqKiqCXgTZoYhNd7KGJTXSxhyY20SU8OLMSIOecmpubeXKXITSxiS720MQmuthDE5voEh4MKwFyzqmuro4HiiE0sYku9tDEJrrYQxOb6BIeDCsAAAAATGJYAQAAAGASw0rAMjO5xoE1NLGJLvbQxCa62EMTm+gSDlQKUCQSUWlpadDLQBs0sYku9tDEJrrYQxOb6BIenFkJkHNOe/bs4cldhtDEJrrYQxOb6GIPTWyiS3gEOqysXr1aU6ZM0fTp0zV37lzt3btXS5cu1dSpUzVr1ixt3749yOX1O+ecYrEYDxRDaGITXeyhiU10sYcmNtElPAIdVkaOHKkXXnhBr7zyij75yU/qySef1KJFi/TSSy/pxz/+sRYuXBjk8gAAAAAEKNDnrBx22GGpv2dlZWnLli0aM2aMsrKydOqpp2r+/Pnd3re5uVnNzc2pl2OxmCTJ9335vi9J8jxPnufJOdducj7Q8db79/Z4JBLp9La7Ou77furv6a7R6p72dzwMe2q9jXOuR7f3PK/d+5A+vo2VPaV7vKs1OucUiUTk+x/fp6ePs77cU9vXDbbPvTDuqW2TwbKn3q59oPfU9jHZ+nLr97+O32fCsqf9HW97v45fg6zvqe33+sHwudfb49b2JHX+Xh/2PYWtU8f33R0TT7B/77339Oyzz+rWW29VTU1N6ngymez2PrfccosWLFjQ6XhNTY3i8bgkKTc3V0VFRYrFYmpqakrdJj8/XwUFBaqrq1MikUgdLywsVF5enmpra9XS0pI6XlxcrOzsbNXU1LT7oJeUlCgjI0PV1dXt1lBeXq5kMqmdO3emjnmep4qKCiUSCdXV1Una9yBpamqS53lqampKDVzSvuFtxIgRamxs1O7du1PHre9J2nd1jdLS0lDuqXWde/fu1a5du3q0J2nf4P2JTxwhaZLi8bhisZiZPfVFp3g8rsrKSr37rlRZWal4PJ563wOxJ9/31djYmNrPYPzcC9ueampq2jUZDHsKU6d4PK5JkyZpy5YcxWJKPSadc2psbFRpaak8zwvVnvbXKR6PKycnR4mEp0mTJrX7GmR9T8651Net1katwvi5t79OYdpTbm6umpqaVFNTkxpewr6nsHVqaGhQT3iu44g0wGKxmM4991w98MADSiaTuuOOO/Tggw9KkqZMmaJXX321y/t1dWZl5MiRqqurU2FhoaShMZWyp+D3tG7dOk2ePFmeN1vSn7R69WpNnDgx1HvquMZ169bp5JNPlu+fo0jkaa1evVoTJkwI9Z4O5jh7Yk9B76ntY1L6o9asWZN6TIZ1T/s73vp1Vprd6WtQWPe0v+PsiT0NhT3FYjEVFxervr4+9W/3rgR6ZqWlpUUXXnihbrjhBo0aNUp79+7VW2+9pUQioTVr1mjcuHHd3jc7O1vZ2dmdjkciEUUi7Z+K0/rB6ai74x3v35vjPXmfzu37H7BDDjkk7TVa3dNAHu+PPTnn1NDQoEMOOSStj8HHD0i/3W0s7Kkvjrf9gub7+/bY08dZX+xJknbv3p16rPTFng7meG/3VFVVpWg0qtLSUlVWVob68eR5XruvX0GuPWyPp4M53rrGjv/IaH1Mtv2+ErY97e942/t19TXI8p7aNhkMn3sHc9zSnpxzXX5f6W7t3R23tKd0197d8YFae3dr6ijQYeXRRx/Va6+9poULF2rhwoW68sorde2112rmzJnKycnR4sWLg1xev2t9oOTn53cZEwOPJjYNhi5VVVUaNWq04vEm5eTkasuWzaqsrAx6Wb02GJoMRnSxhyY20SU8Ah1WLr74Yl188cWdjl9wwQUBrAYA+k80GlU83iRpuuLxVxSNRkM9rAAAMBD4pZAAMKCKgl4AAAChwbASIM/zlJuby+lHQ2hiE13soYlNdLGHJjbRJTxMXLp4qPI8T0VF/C+rJTSxiS720MQmuthDE5voEh6cWQmQc0719fWdLgOH4NDEJrrYQxOb6GIPTWyiS3gwrATIuX2/FJIHih00sYku9tDEJrrYQxOb6BIeDCsAAAAATGJYAQAAAGASw0qAPM/jlxEZQxOb6GIPTWyiiz00sYku4cHVwALkeZ4KCgqCXgbaoIlNdLGHJjbRxR6a2ESX8ODMSoCcc6qtreXJXYbQxCa62EMTm+hiD01sokt4MKwEyDmnRCLBA8UQmthEF3toYhNd7KGJTXQJD4YVAAAAACYxrAAAAAAwiSfY94OqqipFo1GVlpaqsrKy29t5nqfCwkKuRGEITWyiiz00sYku9tDEJrqEB8NKH6uqqtKoUaMVjzcpJydXW7Zs7nZg8TxPeXl5A7xC7A9NbKKLPTSxiS720MQmuoQHPwbWx6LRqOLxJknTFY83KRqNdntb3/cVjUbl+/7ALRD7RROb6GIPTWyiiz00sYku4cGw0m+KenSrlpaWfl4H0kUTm+hiD01soos9NLGJLuHAsAIAAADAJIYVAAAAACYxrATI8zwVFxdzJQpDaGITXeyhiU10sYcmNtElPLgaWIA8z1N2dnbQy0AbNLGJLvbQxCa62EMTm+gSHpxZCZDv+9qxYwdXojCEJjbRxR6a2EQXe2hiE13Cg2ElYM65oJeADmhiE13soYlNdLGHJjbRJRz4MTAAA6aqqkrRaFSlpaXd/rJUAACAVgwrAAZEVVWVRo0arXi8STk5udqyZTMDCwAA2C9+DCxAnueppKSEK1EYQpP+E41GFY83SZqueLxJ0Wi0x/eliz00sYku9tDEJrqEB2dWAuR5njIyMnigGEKTgVCU9j3oYg9NbKKLPTTpnf7+sWG6hAdnVgLk+76qq6u5EoUhNLGJLvbQxCa62EOT9LX+2PCkSZM0atRoVVVV9fn7oEt4cGYFAID9aP0fXklcHAIYAO1/bPgVRaNRHndDGMMKAADdaHthCElcHAIYUOn/2DAGH34MDACAbnz8P7xfkvSltC8OAQA4OJxZCVAkElF5ebkiEWZGK2hiE13sGXpNSoNeQI8MvS720cQmuoQHhQLknFMymeQ3qBpCE5voYg9NbKKLPTSxiS7hwbASIOecdu7cyQPFEJrYRBd7aGITXeyhiU10CQ+GFQAAAAAmMawAAAAAMIlhJWD85lR7aGITXeyhiU10sYcmNtElHLgaWIAikYgqKiqCXgbaoIlNdLGHJjbRxR6a2ESX8ODMSoCcc2pububJXYbQxCa62EMTm+hiD01sokt4MKwEyDmnuro6HiiG0MQmuthDE5voYg9N+ldVVZXWrl2rqqqqtO5Hl/Dgx8AAAAAQOlVVVRo1arTi8Sbl5ORqy5bNqqysDHpZ6GOcWQEAAEDoRKNRxeNNkqYrHm9SNBoNeknoBwwrAcvM5OSWNTSxiS720MQmuthDk/5W1Kt70SUcqBSgSCSi0tLSoJeBNmhiE13soYlNdLGHJjbRJTw4sxIg55z27NnDk7sMoYlNdLGHJjbRxR6a2ESX8GBYCZBzTrFYjAeKITSxiS720MQmuthDE5voEh4MKwAAAABMYlgBAAAAYBJPsA+Q53nKysqS53lBLwUfoYlNdLGHJjbRxZ4wNqmqqlI0GlVpaemg/b0lYewyVDGsBMjzPI0YMSLoZaANmthEF3toYhNd7Albk6HyixbD1mUo48fAAuScU0NDA0/uMoQmNtHFHprYRBd7wtZkqPyixSC6VFVVae3atVq7dq2qqqoG7P2GHWdWAuSc0+7du5Wfn89pSCNoYhNd7KGJTXSxJ7xNeveLFsNioLu0PWMlaVCfteprnFkBAAAA+tHHZ6y+JOlLg/qsVV/jzAoAAEAahsIT0NFfSoNeQOgwrATI8zzl5uaG7LTw4EYTm+hiD01sokv/S/cJ6DSxiS7hwY+BBcjzPBUVFfFAMYQmNtHFHprYRJf+l+4T0Gli01Dt0vok/zA9wZ9hJUDOOdXX14fmCiFDAU1soos9NLGJLgOpZ09Ap4lNQ7FL61nBSZMmadSo0aEZWMwOK9dff71OO+00XXzxxdq7d2/Qy+kXzjk1NTUNqQeKdTSxiS720MQmuthDE5uGYpewXpba5LCyYcMG/fvf/9by5cs1evRoPf7440EvCQCGnDD+uAAA4EDCdVlqk0+wX7lypT772c9Kkj73uc/poYce0ty5c9vdprm5Wc3NzamX6+vrJUm7du2S7/uS9v08oud5cs61m5wPdLz1/r053tDQ8NHfdsrzPDU0NGjXrl2SpEgk0u59+r6vWCymnJyc1Hp6usaB3FNXa+/N8TDsyfd9NTQ0KDs7u9PPsXa3p8bGxo/+XitpX/P6+noze0r3eFdrbGho+OhY58/rnu6p7dtofbn1bRxoT8lkUrFYTFlZWYpEIn2yp4M53ptOrfuXauVc5/1b29P27ds1efJJSiSalZOTq9WrX9MRRxyRun1LS0u7JoP1a8TH3T5svbUaGxu7/PGRgdxTd4+n1u8r2dnZysjIGDRfy9P53trfe9rfWrpa+4cffqj33ntPn/zkJ1VRUXHAvQaxp7bH236tav2eFovFBvT7U3cf445rP9DX1f3t1TnX6fvKwOzpg4/WfHD/XujN46njx6uxsfGgvg91d7yna4/FYpLU6W124gy6+eab3bJly5xzzm3dutXNnTu3021uuOEGJ4k//OEPf/jDH/7whz/84U9I/2zbtm2/c4HJMyvDhw9PTVv19fUaMWJEp9t8//vf13e+853Uy77vq7a2ViUlJZ3+R9yqWCymkSNHatu2bSosLAx6ORBNrKKLPTSxiS720MQmugTPfXR27PDDD9/v7UwOK1OnTtWiRYt0ySWX6K9//atOPfXUTrfJzs5WdnZ2u2PDhw8foBX2rcLCQh4oxtDEJrrYQxOb6GIPTWyiS7CKiooOeBuTT7AfP368KioqdNppp+kf//iHzjvvvKCXBAAAAGCAmTyzIkm333570EsAAAAAECCTZ1aGiuzsbN1www2dfpwNwaGJTXSxhyY20cUemthEl/DwnDvQ9cIAAAAAYOBxZgUAAACASQwrAAAAAExiWAEAAABgEsPKAFq9erWmTJmi6dOna+7cudq7d6+WLl2qqVOnatasWdq+fXvQSxySduzYoalTp2rGjBk644wz9MEHH2jFihWaOnWqpk2bpo0bNwa9xCHr0UcfVVlZmSTxWDHg3XffVVlZmWbOnKmZM2eqpqaGLga89NJLmjVrlk4//XQtW7aMr18GvPrqq6nHyXHHHadvf/vbdAmY7/uaN2+eTjvtNE2bNk2bN2+mSVjs9/fbo0+9//77bs+ePc455773ve+5pUuXulNOOcU1Nze7FStWuMsvvzzgFQ5NLS0tLplMOuece+ihh9zChQvd9OnTXW1trXvvvffcWWedFfAKh6aWlhY3Z84cN2HCBLd3714eKwa888477rzzzku9TJfg7dmzx5177rmuubk5dYyvX7Zceuml7qWXXqJLwF5//XV34YUXOuece+WVV9zXvvY1moQEZ1YG0GGHHabc3FxJUlZWlrZs2aIxY8YoKytLp556qt54442AVzg0ZWRkKBLZ91BoaGjQMccco4yMDBUXF6uyslK1tbUBr3BoevTRR3X++ecrEolo69atPFaM+Nvf/qbTTjtNP/jBD+hiwKuvvqrc3FzNnj1bc+bM0QcffMDXL0MSiYRWr16tyZMn0yVgRxxxhJxzcs6prq5O+fn5NAkJhpUAvPfee3r22Wc1bdo0FRYWpo4nk8kAVzW0rV+/Xp/+9Kd1zz33aOrUqe26ZGZmKpFIBLi6oSeZTGrJkiW64IILJEl1dXU8Vgw47LDD9Pbbb+uVV15RdXW1/vCHP9AlYDt27NDbb7+tP/7xj/ra176mG264ga9fhjz33HOaNWuW6uvr6RKw0tJSDRs2TKNHj9Y3v/lNffvb36ZJSDCsDLBYLKaLL75YDz/8sMrKyhSLxVKvy8jICHBlQ9v48eP12muvaeHChbr55pvbdWlpaVFWVlaAqxt6HnnkEX3lK19JnfEaPnw4jxUDsrOzlZ+fL8/z9KUvfUkbNmygS8CGDx+uU089VVlZWZo1a5bWrVvH1y9Dli5dqvPPP7/T1zC6DLxnn31WmZmZ2rJli37/+9/ru9/9Lk1CgmFlALW0tOjCCy/UDTfcoFGjRunYY4/VW2+9pUQioZUrV2rcuHFBL3FIavs/KUVFRTrkkEPU0tKiXbt2adu2bRoxYkSAqxuaNm3apF//+tf63Oc+p61bt+ruu+/msWJAQ0ND6u/Lly/XOeecQ5eAnXTSSXrrrbfknNP69es1duxYvn4ZsXfvXv3973/XtGnTlJeXR5eAOedUUlIiad9ZloaGBpqERGbQCxhKHn300dT/3i9cuFBXXnmlrr32Ws2cOVM5OTlavHhx0EscktavX6/58+crIyNDOTk5evDBB7V161adffbZ8jxP9957b9BLHHJuu+221N8nT56s++67T4899hiPlYCtWLFC//3f/628vDwdddRRWrhwoXJycugSoNLSUs2ZM0czZsyQ53l68MEH9e9//5uvXwY899xzOuOMM1JniG+66Sa6BOjMM8/Uww8/rBkzZqi5uVmLFi1SS0sLTULAc865oBcBAAAAAB3xY2AAAAAATGJYAQAAAGASwwoAAAAAkxhWAAAAAJjEsAIAAADAJIYVAAAAACYxrAAAAAAwiWEFADBgtm3bppkzZ2rs2LEaN26cli5dmnrdG2+8oQkTJmjUqFG64447un0bc+bM0fvvvy9JeuSRRzRx4kQ98sgjkqQtW7Zo3rx5/boHAMDAYVgBAAyYzMxM3XXXXdq0aZOeffZZXXvttdq9e7ck6Tvf+Y4ef/xxbdq0SX/961/19ttvd7r/+vXrlZubq8MPP1yS9Pvf/16vvfaaHn/8cUnSqFGjtGPHDm3fvn3gNgUA6DcMKwCAAXPYYYdp/PjxkqRDDz1UpaWlqq2tVXV1tXzf1zHHHKOMjAzNnTtXy5Yt63T/3/72t/rCF76Qetk5J8/z5Hle6tg555yjJUuW9PteAAD9j2EFABCI119/XclkUiNHjtS2bdu0du1ajR8/XuPHj9dNN92kqqqqTvdZtWqVJk6cmHp59uzZOumkkzRnzpzUsYkTJ2rlypUDsgcAQP/KDHoBAIChp7a2VpdccokeeOCB1LFp06bpT3/6kyRpyZIlWr58eaf7ffjhhyorK0u9fNlll+myyy5rd5uysjJ98MEH/bRyAMBA4swKAGBANTc364tf/KK+973vaerUqZKkkSNHpp40L+17In5lZWWn++bk5Cgej+/37cfjceXm5vbtogEAgWBYAQAMGOec5s2bpzPOOEMXX3xx6nh5ebkKCwu1bt06JRIJ/fa3v9WXvvSlTvcfPXp0l0+8b+vtt9/WmDFj+nztAICBx7ACABgwf/vb3/TYY4/piSeeSD0/ZePGjZKku+++W5deeqnGjBmj//iP/9AxxxzT6f6f+9zn9PLLL+/3fbz88ss666yz+mX9AICB5TnnXNCLAACgJxobG3XWWWfplVdeaXcFsFaJREJnnHGGXn75ZWVkZASwQgBAX2JYAQCEyjPPPKOJEyfq0EMP7fS6d955R++9955mzpw58AsDAPQ5hhUAAAAAJvGcFQAAAAAmMawAAAAAMIlhBQAAAIBJDCsAAAAATGJYAQAAAGASwwoAAAAAkxhWAAAAAJjEsAIAAADApP8PHg44ONNb5K0AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 800x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "hkl indices: \n",
      "[[{'hkl': (1, 1, 1), 'multiplicity': 8}], [{'hkl': (2, 2, 0), 'multiplicity': 12}], [{'hkl': (3, 1, 1), 'multiplicity': 24}], [{'hkl': (2, 2, 2), 'multiplicity': 8}], [{'hkl': (4, 0, 0), 'multiplicity': 6}], [{'hkl': (3, 3, 1), 'multiplicity': 24}], [{'hkl': (4, 2, 2), 'multiplicity': 24}], [{'hkl': (5, 1, 1), 'multiplicity': 24}, {'hkl': (3, 3, 3), 'multiplicity': 8}], [{'hkl': (4, 4, 0), 'multiplicity': 12}], [{'hkl': (5, 3, 1), 'multiplicity': 48}], [{'hkl': (6, 0, 0), 'multiplicity': 6}, {'hkl': (4, 4, 2), 'multiplicity': 24}], [{'hkl': (6, 2, 0), 'multiplicity': 24}], [{'hkl': (5, 3, 3), 'multiplicity': 24}], [{'hkl': (6, 2, 2), 'multiplicity': 24}], [{'hkl': (4, 4, 4), 'multiplicity': 8}], [{'hkl': (5, 5, 1), 'multiplicity': 24}, {'hkl': (7, 1, 1), 'multiplicity': 24}]]\n"
     ]
    }
   ],
   "source": [
    "from pymatgen.core import Structure\n",
    "from pymatgen.analysis.diffraction.xrd import XRDCalculator\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Load structure from CIF\n",
    "structure = Structure.from_file(\"1526734.cif\")\n",
    "print(f'lattice infos \\n{structure.lattice}')\n",
    "# Calculate XRD pattern\n",
    "xrd = XRDCalculator(wavelength=\"CuKa\")\n",
    "pattern = xrd.get_pattern(structure)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8, 4), constrained_layout = True)\n",
    "ax.bar(pattern.x, pattern.y, width=0.2, color='navy', edgecolor='black')\n",
    "ax.set_xlabel(\"2θ (°)\")\n",
    "ax.set_ylabel(\"Intensity (a.u.)\")\n",
    "ax.set_title(\"Simulated XRD Pattern\")\n",
    "ax.grid(True, linestyle='--', alpha=0.3)\n",
    "plt.show()\n",
    "print(f'hkl indices: \\n{pattern.hkls}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e855b402",
   "metadata": {},
   "source": [
    "# Perform the Analysis of a PXRD Data Set\n",
    "1. Load the data from CPY-131_repro_calc and perform a background correction\n",
    "2. Calculate the XRD pattern for $\\mathrm{Co_3O_4}$\n",
    "3. Fit Psydo-Voigt profiles to the diffraction peaks\n",
    "4. Make a plot of the raw data + fits and save the results to .txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "id": "077f6e53",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "amplitude\t    79.42476\t     0.67533\n",
      "center \t    19.05998\t     0.00000\n",
      "sigma  \t     0.34380\t     0.00253\n",
      "fraction\t     0.41888\t     0.02408\n",
      "fwhm   \t     0.68760\t     0.00507\n",
      "height \t    93.86258\t     0.50521\n",
      "0.01200089133821976\n",
      "amplitude\t    89.16615\t     1.01805\n",
      "center \t    31.37194\t     0.00000\n",
      "sigma  \t     0.39988\t     0.00354\n",
      "fraction\t     0.31401\t     0.03283\n",
      "fwhm   \t     0.79977\t     0.00709\n",
      "height \t    94.13682\t     0.63777\n",
      "0.01395855432340036\n",
      "amplitude\t   404.22082\t     5.14940\n",
      "center \t    36.96689\t     0.00000\n",
      "sigma  \t     0.40192\t     0.00421\n",
      "fraction\t     0.33292\t     0.03731\n",
      "fwhm   \t     0.80385\t     0.00842\n",
      "height \t   421.70951\t     3.26528\n",
      "0.01402974659834673\n",
      "amplitude\t    26.54079\t     1.69446\n",
      "center \t    38.67448\t     0.00000\n",
      "sigma  \t     0.26521\t     0.02468\n",
      "fraction\t     1.00000\t     0.00000\n",
      "fwhm   \t     0.53043\t     0.04937\n",
      "height \t    31.85420\t     1.98764\n",
      "0.009257733331296685\n",
      "amplitude\t    99.51826\t     1.18770\n",
      "center \t    44.95907\t     0.00000\n",
      "sigma  \t     0.38650\t     0.00389\n",
      "fraction\t     0.32507\t     0.03566\n",
      "fwhm   \t     0.77300\t     0.00778\n",
      "height \t   108.27313\t     0.79110\n",
      "0.013491343452559099\n",
      "amplitude\t    23.68307\t     1.77148\n",
      "center \t    55.84610\t     0.00000\n",
      "sigma  \t     0.64389\t     0.07497\n",
      "fraction\t     1.00000\t     0.00000\n",
      "fwhm   \t     1.28777\t     0.14993\n",
      "height \t    11.70788\t     0.86341\n",
      "0.0224759097499496\n",
      "amplitude\t   117.59946\t     2.48989\n",
      "center \t    59.56272\t     0.00000\n",
      "sigma  \t     0.48343\t     0.00779\n",
      "fraction\t     0.22684\t     0.06405\n",
      "fwhm   \t     0.96685\t     0.01559\n",
      "height \t   105.90962\t     1.28945\n",
      "0.016874790182978414\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<>:146: SyntaxWarning: invalid escape sequence '\\m'\n",
      "<>:148: SyntaxWarning: invalid escape sequence '\\T'\n",
      "<>:159: SyntaxWarning: invalid escape sequence '\\m'\n",
      "<>:146: SyntaxWarning: invalid escape sequence '\\m'\n",
      "<>:148: SyntaxWarning: invalid escape sequence '\\T'\n",
      "<>:159: SyntaxWarning: invalid escape sequence '\\m'\n",
      "/var/folders/f2/krn3py8536556jbm969hx5380000gn/T/ipykernel_12269/1062974832.py:146: SyntaxWarning: invalid escape sequence '\\m'\n",
      "  axs.plot(two_theta, counts_smoothed, color = colors_all[1], linestyle = '', marker = 'o', fillstyle = 'none', markersize = 2, label = '$\\mathrm{Co_3O_4}$ nanoparticles')\n",
      "/var/folders/f2/krn3py8536556jbm969hx5380000gn/T/ipykernel_12269/1062974832.py:148: SyntaxWarning: invalid escape sequence '\\T'\n",
      "  axs.set_xlabel('2 $\\Theta$ (°)')\n",
      "/var/folders/f2/krn3py8536556jbm969hx5380000gn/T/ipykernel_12269/1062974832.py:159: SyntaxWarning: invalid escape sequence '\\m'\n",
      "  axs.bar(pattern.x, pattern.y*np.max(counts_smoothed)/200, width=0.2, color='navy', edgecolor='black', label = '$\\mathrm{Co_3O_4}$ COD ID 1526734', zorder=10)\n",
      "/opt/homebrew/anaconda3/lib/python3.12/site-packages/uncertainties/core.py:1024: UserWarning: Using UFloat objects with std_dev==0 may give unexpected results.\n",
      "  warn(\"Using UFloat objects with std_dev==0 may give unexpected results.\")\n",
      "/opt/homebrew/anaconda3/lib/python3.12/site-packages/uncertainties/core.py:1024: UserWarning: Using UFloat objects with std_dev==0 may give unexpected results.\n",
      "  warn(\"Using UFloat objects with std_dev==0 may give unexpected results.\")\n",
      "/opt/homebrew/anaconda3/lib/python3.12/site-packages/uncertainties/core.py:1024: UserWarning: Using UFloat objects with std_dev==0 may give unexpected results.\n",
      "  warn(\"Using UFloat objects with std_dev==0 may give unexpected results.\")\n",
      "/opt/homebrew/anaconda3/lib/python3.12/site-packages/uncertainties/core.py:1024: UserWarning: Using UFloat objects with std_dev==0 may give unexpected results.\n",
      "  warn(\"Using UFloat objects with std_dev==0 may give unexpected results.\")\n",
      "/opt/homebrew/anaconda3/lib/python3.12/site-packages/uncertainties/core.py:1024: UserWarning: Using UFloat objects with std_dev==0 may give unexpected results.\n",
      "  warn(\"Using UFloat objects with std_dev==0 may give unexpected results.\")\n",
      "/opt/homebrew/anaconda3/lib/python3.12/site-packages/uncertainties/core.py:1024: UserWarning: Using UFloat objects with std_dev==0 may give unexpected results.\n",
      "  warn(\"Using UFloat objects with std_dev==0 may give unexpected results.\")\n",
      "/opt/homebrew/anaconda3/lib/python3.12/site-packages/uncertainties/core.py:1024: UserWarning: Using UFloat objects with std_dev==0 may give unexpected results.\n",
      "  warn(\"Using UFloat objects with std_dev==0 may give unexpected results.\")\n",
      "/opt/homebrew/anaconda3/lib/python3.12/site-packages/uncertainties/core.py:1024: UserWarning: Using UFloat objects with std_dev==0 may give unexpected results.\n",
      "  warn(\"Using UFloat objects with std_dev==0 may give unexpected results.\")\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "amplitude\t   154.65847\t     3.00478\n",
      "center \t    65.46660\t     0.00000\n",
      "sigma  \t     0.46648\t     0.00666\n",
      "fraction\t     0.18449\t     0.05907\n",
      "fwhm   \t     0.93296\t     0.01333\n",
      "height \t   146.47068\t     1.60472\n",
      "0.016283287037994153\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWkAAADmCAYAAAAX6nT3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA7MUlEQVR4nO3de1xUdf748ddwGQYvjCCoqJhpppB3RRNEUMpM07S85Lqm3+2+fbPcr9tFLSvWtXV3u2puu2nmmpaYpq3+yjLxkneT1FBDV0FFvCHDKA4ww/n9MTAyw8xwEWYO8H76mIecy5zzPjPDm898zueiURRFQQghhCr5eDsAIYQQrkmSFkIIFZMkLYQQKiZJWgghVEyStBBCqJgkaSGEUDFJ0kIIoWKSpIUQQsUkSQshhIp5NEmnpKSQmJjI4MGDWbt2LTt27CAmJoaBAwdy+PBhALKzsxk6dCixsbEsX77ck+EJIYTqaDzVLfzGjRuMHz+eL7/8Eq1WC0B8fDxfffUVRqORp59+mo0bNzJ9+nRGjBhBQkICcXFxbNmyBZ1O54kQhRBCdfw8daJdu3YRGBjIyJEjadSoER9++CG+vr4EBwcTHBxMTk4OAHv37uXvf/87Pj4+9O3blyNHjtC3b1+7YxUXF3P69Gn8/f3RaDS29QEBAQQEBHjqkoQQokKKomA0GmndujU+PlWvvPBYkr5w4QInTpxg9+7dfP/998yZM4egoKCbgfj5UVhYSFFRke1C9Hq9LXmXlZWVRceOHT0VuhBC3LIzZ87Qtm3bKj/PY0m6WbNmxMbGotVqSUxMZM6cOTRt2tS23Ww2o9Vq8ff3p7i4GB8fHwwGAyEhIeWOVfq8tLQ0u2NISVoIoTZ5eXlERETY5aqq8FiSjo6O5u9//zuKopCamkpUVBSnTp0iNzcXo9FoS8bR0dGkpKQwaNAgDhw4wPz588sdq7SKo02bNnalcSGEUKuyVbNV4bEkHRoaypgxY4iPj0ej0bBkyRLOnTvH8OHD0Wg0fPjhhwC89NJLPProo8yePZunn36awMBAT4UohBCq47HWHTUpLy8PvV6PwWCQkrQQQtVuNV9JZxYhhFAxSdJCCKFikqSFEELFJEkLIYSKSZIWQggVkyQtRB23c+dOEhISiI+PZ8iQIezfv7/Sz01JSSEuLo74+HgeeeQRrl69WouR1qzTp0+zadMmAFJTU1m0aJHT/VJSUpgxY4YnQ6tRHmsnLYSoeTk5OTzzzDN88803hIeHYzAYOHnyZKWfO23aNDZv3kxYWBgrV67kueeeqxOjT5aO37Np0yaGDh1Kz5496dmzp7fDqhVSkhbCW9avh+nTrf9X04YNGxg9ejTh4eGAdbybHj168Nvf/pb4+HhGjBjhsnS8YcMGxowZQ1hYGAATJ05k9+7dFBcXuzxfSkoKw4YNY8yYMfTo0YMjR45w4cIFBg8eTFxcHGPHjsVisTjdD8BisZSLLSUlhaFDhzJy5Eiio6M5fPiwy2OOHDmSMWPGsHTpUhYtWsQXX3xBQkICa9asYcaMGSiKwrPPPktcXByDBw/m0qVLttgVReG5555j8ODB3HPPPZw9e5bdu3fTv39/Bg8ezOuvv17t96E2SZIWwhvWr4cHH4RVq6z/VzNRZ2Vl0bp1a7t1a9eupW3btmzdupVHHnmEDz74oNLPDQsLs0tszhQVFbF27VreeustlixZQnBwMN999x3bt2+nTZs2/PDDD073cxdbfn4+69evZ9myZcyaNcvlMQ0GA2vWrOF3v/sdzzzzDBMmTCAlJcU2rMTXX3+Nj48P27dvZ8uWLTRv3twW94YNGwgODmbLli3MnTuXt956iw0bNjBnzhy2bNnCa6+9VtmX3aPqdJKOjo4mKiqKhQsXejsUIapmyxZo3RrOnrX+n5JSrcO0bt2ac+fO2a07ceIE0dHRgPV3JD09HYBff/2Vd955h9///vdcvnyZ8PBwsrKy7J578eJFQkND+eabbxg1apTTc5ZWK0RERHD16lWuXLnC2LFjiY+PZ+PGjbZjOu7nLrZevXqh0WiIjIzk/PnzLo/Zt29ft2NgHD16lPj4eNty2aFB09LSWLt2LQkJCbz44ovk5uby7LPPsnHjRiZNmsQ333zj8rjeVKeT9L59+0hLS+PZZ5/1dihCVM3gwZCVBW3bWv9PSKjWYUaMGMG6des4f/48YO2C7OPjw969ewHr70inTp0AuPPOO2ndujXZ2dn4+/szYsQI1qxZYys5r1y5krvvvpvjx4+Tl5dHhw4dnJ6zbJJUFIUVK1bwwAMPsHXrVoYNG0bpSBOO+wHccccdTmNLTU1FURSOHz9OeHi4y2OWTbr+/v5YLBa72CIjI9m2bZttuWzVTZcuXRg/fjwpKSls3bqVTz75BL1ez4IFC/jkk0946aWXKveie5jcOBTCG0aNgnXrrCXohATrcjWEhISwaNEiJk6ciKIo+Pr6Mm/ePN577z0GDRpEkyZN7G4ETpgwgaZNm5KZmUm3bt147733eOihh9BoNLRq1YpFixaxbNkyNBoNBw8e5Oeff6ZHjx5uY0hMTGTy5Ml8/fXXFQ6INnr0aNasWWMX26FDh9Dr9YwcOZILFy6wePFiiouLKzxmt27deOWVVxg3bhyTJk0CYOTIkXzzzTcMHDgQf39/Vq1aZdt/5MiR/PDDDwwePBiNRsOkSZMwGo2sWbMGs9nM1KlTK/mqe5YMsCREA/HDDz+wf/9+Tp48yRtvvEGrVq3c7v/CCy/w7rvv1npcKSkp/Oc//+Fvf/tbrZ/LG241X0lJWogGYsiQIQwZMqTS+3siQYuKSUlaCCFqkQxVKoQQ9ZgkaSGEUDFJ0kIIoWKSpIUQQsUkSQshhIp5LEmfPn2asLAwEhISSEhI4NKlSyQnJxMTE0NiYiJnz54F4NixYwwaNIiYmBg2b97sqfCEEEKVPNpOOj4+ntWrVwNgNpt5++232bp1K/v27SMpKYmPPvqImTNnsnjxYlq2bMn9999PYmKiJ0MUQghV8WiS/vHHH4mLiyMuLo7JkycTGRmJVqslNjbWNih3VlaWrT9/SEgIly9fJjQ01Onx8vLy7JYDAgIICAio3YsQopIyMzO5fPnyLR0jNDSUdu3aud1n586dzJw509YtfP78+fTt27dSx09JSeHVV1/Fx8eH8PBwFi1aRHBwcLXO5+xYBoOB6OhounXrhtlsJjo6mqSkJBo1amQ73unTp5kxYwarV6/m9OnTFe5vMBi49957SUtLY/fu3XTt2tX2vLvuuguA5ORkwsLC2Lt3L88//zz+/v60adOGZcuW4e/vT0pKCklJSRQXFzNt2jRatWrFK6+8Alhz0IgRI3jnnXe4cOECY8aMwd/fH19fXz777DPbsLArV65k2rRpFY4aeMsUDzGZTMq1a9eU4uJi5bHHHlP+9Kc/Kc8//7xte3R0tN3/iqIokyZNUo4fP17uWAaDQQHKPebMmVPblyFEpWRkZCg6XaDTz2lVHjpdoJKRkeHyPFeuXFG6d++uZGVlKYqiKLm5ucqBAwcqFeOVK1eUbt26KRcvXlQURVFWrFihTJo0qcLnODufq2OdOnVKefjhhxVFUZTi4mJl9uzZyowZM+yOWXafyuxfWFioXLx4UZkyZYpy+PDhcs8rKysrS8nPz1cURVFefvllJTk5WcnPz1ceeOABpaCgwOk1TpkyRUlJSVEURVHMZrNisVgURVGUTz75RElKSrKtHzNmjNKrVy+3r5ei3MxXBoOhwn2d8ViddEBAAI0bN0aj0fDQQw/x888/25WEfX19AftRrgwGg22cWGfOnDmDwWCwPUr/EgrhbZcvX8ZkugE8BDxZzcdDmEw33JbGPT3ov7Pz9e7du1LH0mg0vPrqq6yv5NjZrvb39/e3naes0m/qpaV8gPDwcNsATVqtFh8fH3bt2kVgYKBtAoHs7GzbMQoLC9m7dy9xcXGANS+V5iSj0Wgrqa9cuZJx48bZ5ava4rEkbTQabT9v376dESNGcPToUQoLC9m5cyfdu3cHrC/qyZMnMRqN5OTkuKzqAAgKCrJ7SFWHUJ9QoHU1H64/+6U8Pei/s+dU5VharZbCwkK311Sd/cPDwzlx4gTbtm3j4sWLrFmzxm57RkYGmzZtso20d+LECb7++mueeOIJuxlZvv/+exITE+2Sb2pqKv3792fBggX07t0bi8XCqlWrmDBhQqWv41Z4LEnv2LGDPn36EBcXx7lz55g0aRIvvPACCQkJzJ49m9mzZwMwd+5cpk6dyn333afa6WyEUAtPD/rv7HyA22OVVVBQUKXCVGX3d/ZNvVReXh6TJ09m6dKl+Pv706xZM2JjY9FqtSQmJvLLL7/Y9k1OTmbcuHF2x+7Zsyd79uwhKSmJefPmsXz5csaPH++RUjR4MEnff//9HDhwgO3bt7Ns2TL8/PyYMGECO3fu5IcffiAiIgKAqKgotm/fzs6dO7n33ns9FZ4QdZKnB/13dr6ffvrJ5bFKqzFLzZs3j9GjR1f6+iq7v+M39TvuuAOwtiJ75JFHmDNnDp07dwasf7iOHj2KoiikpqbarrOoqIh9+/YxcOBA27HKluL1ej2NGjUiLS2NZcuWMWzYMNLT05k2bVqlr6c6ZKhSIWrVrbTuqPi5nh7039n5/vrXv9K8eXOnxzIajWzdupXBgwdjsVjo378/b775pttrqsz+w4cPJzU1lePHj/PUU0/RsmVLZs+eTaNGjbj99ttJSkoCrH8sSkvBSUlJtnkRx4wZQ3x8PBqNxjb/4vfff8+QIUPKVXXMmDEDX19fdDodS5YssdXHg3U6r/fff7/C9+lWyFClQtSCzMxMOnfuUnLzsPp0ukCOHz9WYTO8ylDroP/13a3mK0nSQtQST7WTFuomM7MIoVLt2rWTBCtumQywJIQQKiZJWgghVEyStBBCqJgkaSGEUDFJ0kIIoWJ1OklHR0cTFRXFwoULvR2KEELUijrdBG/fvn3STloIUa/V6ZK0EELUd5KkhRBCxSRJCyGEikmSFkIIFZMkLVQlDbjn0CF6HT5M8t//7u1whPC6Ot26Q9Q/U0+cYF/JVGoTIyPp98wz3LZokZejEsJ7pCQtVCMH2FcyowaAxc+PDTJvpWjgJEkL1Ujftq3cukt6vRciEUI9JEkL1bj6ww/l1uWEhHghEiHUw+NJeuXKlYSFhQHWmXljYmJITEzk7NmzABw7doxBgwYRExPD5s2bPR2e8CJDmUk/S11xmG1aiIbGozcOLRYLycnJREREYDabefvtt9m6dSv79u0jKSmJjz76iJkzZ7J48WJatmzJ/fffT2JioidDFF5kKCoqt85YMou8EA2VR0vSK1euZNy4cfj4+JCenk5kZCRarZbY2FgOHToEQFZWFp06dSIoKIiQkBC3c8Tl5eXZPQoKCjx1KaIWGPz9y60zNmrkhUiEUA+PJWmLxcKqVauYMGECAFevXrUbHMlisQBQXFxsW6fX68nJyXF5zIiICPR6ve0xb968WopeeIKhbdty664VF8P69V6IRgh18Fh1x/Llyxk/fjw+Pta/C82aNSMvL8+23dfXF8C2HcBgMBDi5sbRmTNn7BJ9gDTXqtNye/cut84YFASffw6jRnkhIiG8z2NJOi0tjYMHD7J8+XLS09P54IMPOHr0KIWFhezfv5/uJR0YwsPDOXnyJC1atCAnJ4dQNzeOgoKCZKjSesSYn19+XZMmEBjohWiEUAePJem//OUvtp/79u3LokWL+OKLL0hISECn0/Hpp58CMHfuXKZOnYrFYuGNN97wVHhCBUzO6qSbNoUbN7wQjRDq4JVu4fv37wdgwoQJtjrqUlFRUWzfvt0bYQkvMzn51mTS6aQkLRo06cwiVMPkpMRcGBCAJS3NC9EIoQ6SpIVqFDip7gDn1SBCNBSSpIVqmFyM02EquaksREMkSVqohqmkGaajGyVt6IVoiCRJC9UwabXO1zdp4uFIhFAPSdJCNQp8nH8cb2RkeDgSIdRDkrRQDZclablxKBowSdJCNVzWSXfp4uFIhFAPSdJCNQr8nPetMsmNQ9GA1ekkHR0dTVRUFAsXLvR2KOIWKYDJRZK+ITcORQNWp2cL37dvnwywVE8UAYqrG4fZ2Z4NRggVqdMlaVF/mNxtkxuHogGrdpI+f/48O3fuRFGUmoxHNFDu5tS50b69p8IQQnWqlKT79euH0WgkOzubAQMG8Kc//Yknn3yytmITDYhp0yaX22SgUtGQVSlJFxUV0bRpU5KTk3nyySfZuHEj+/btq63YRANi+vJLl9sKOnXyYCRCqEuVkrTFYmHz5s18+umnjBw5EgCz2VwrgYmGxdUIeOC6/bQQDUGVkvT777/Pe++9x7hx4+jWrRv//e9/GTx4cG3FJhoQU1SU623nznkwEiHUpUpN8EwmE+vLzNzcoUMHRowYUeNBiYbH5OYbmUmn82AkQqhLlUrSM2fOLLdu1qxZNRaMaLhMffu63qbReDASIdSlUiXp7777jk2bNpGVlcWLL75oW5+Xl4ev1BeKGuBqBDwAU0CAByMRQl0qVZJu0aIFXbt2JSAggLvuusv2GDp0KN9++22lTnThwgViYmKIj49nyJAhnD9/nh07dhATE8PAgQM5fPgwANnZ2QwdOpTY2FiWL19e/SsTdYppyxbX26QkLRqwSpWke/ToQY8ePZg8eTI+bko87oSGhrJjxw58fHxYunQpixcv5rvvvmPDhg0YjUaefvppNm7cyF/+8hdefPFFEhISiIuLY+zYseikTrLeczVMKUidtGjYqpRxf/zxR+655x46depEhw4duP322+nQoUOlnuvr62tL8EajkY4dO+Lr60twcDDt2rUjJycHgL179zJkyBD8/Pzo27cvR44ccXnMvLw8u0dBgbt+a0LNCoYMcblNStKiIatS647HH3+cf/zjH/Tp06daddGpqak89dRT5ObmsmnTJr744oubgfj5UVhYSFFRkS2Z6/V6W/J2JiIiwm55zpw5vP7661WOS3ifqVcvl9sK3JSyhajvqpSkQ0JCbqlddM+ePdmzZw+rVq1i7ty55OXl2baZzWa0Wi3+/v4UFxfj4+ODwWAgJCTE5fHOnDljNwpegNxgqrNMR45A167Ot8nNadGAVSlJ9+/fnylTpjB69Gi7hDh8+PAKn1tYWIi2pESk1+tp0qQJZrOZ3NxcjEajLRlHR0eTkpLCoEGDOHDgAPPnz3d5zKCgIBmqtJ4wZWS4TtJSkhYNWJWStMFgwMfHx65Di0ajqVSSTk1NZcaMGfj6+qLT6ViyZAnp6ekMHz4cjUbDhx9+CMBLL73Eo48+yuzZs3n66acJDAys4iWJuqjAzftsqubNaiHqA41SB8cazcvLQ6/XYzAYpCRdT7y0eTPzExOdbrvt6lVOBwd7OCIhasat5qsqlaT/53/+B42TO+1Lliyp8omFKMvkriTtwTiEUJsqJemxY8fafi4oKOCrr74iWEo4oga4azxp0mhg/XoYNcpj8QihFlVK0o6DKT300EPcfffdNRqQaJjcdVgx6XSQkiJJWjRIVUrS+fn5tp+Li4tJTU3l6tWrNR6UaHjcVWkU6HQogYFIlxbREFUpSd91111oNBoURcHPz4/27dvzz3/+s7ZiEw2IKTzc7faCoiKkc7hoiKqUpE+dOlVbcYgGzt0oeACmpk0lSYsGqUpJurCwkA8//JDt27cDMGjQIJ555hlbJxUhqsuUlwctW7reLtO0iQaqSr0EnnzySX799VdmzJjBjBkzSE9Pl9nCRY2oaMxoU5MmHopECHWpUkn64MGD/Pzzz7blAQMG0LNnz5qOSTRAFVZ3WCweikQIdalSSVqn07F//37b8oEDB7w6qFF0dDRRUVEsXLjQazGImlHR+BwmaeopGqgqlaT/8Y9/MGXKFAoLC1EUBZ1Ox9KlS2sptIrt27dPuoXXExWNdCcj4YmGqkpJeuXKlWzdutU23rPFYmH+/Pn0cjMWsBCVYfL3d7/90CGIi/NQNEKoR5WqO7777juCg4PR6/Xo9XpCQkLYtGlTbcUmGpCCikrSjRt7KBIh1KVKSbq4uBij0WhbzsvLo6ioqMaDEg1PRSXpArlxKBqoKlV3PP/88wwcOJAJEyYA8MUXXzB9+vRaCUw0LBWWpHv08FAkQqhLlZL07373O/r168eWLVsAWLFiBXfddVetBCYaDjNgcUjSjU0mrpcZdEluHIqGqkpJGqBr1650dTHNkRDV4WxwpWYGg32SPnEC5Aa1aIBkXiLhdc7GktaXGXER3E8KIER9JklaeJ2rknRF+wjREEiSFl7nLAEHnDsHmZk397njDs8FJISKeCxJ7927lwEDBjBo0CAmTpxIUVERycnJxMTEkJiYyNmzZwE4duwYgwYNIiYmhs2bN3sqPOFF5ao7MjPZ+tDD0LmLLVHLjUPRUHksSUdERPDDDz+wbds22rdvz7p163j77bdJSUnhzTffJCkpCYCZM2eyePFivvnmG1577TVPhSe8yLEk7Xf+PMWFBWC6AZcvW/fJyvJ8YEKogMeSdHh4OIElN3+0Wi3Hjx8nMjISrVZLbGwshw4dAiArK4tOnToRFBRESEgIl0t+SZ3Jy8uzexQUuJvOVKiVY5LWOnkf3c2BKER95vE66YyMDDZt2sTAgQPtBkeylPQoKy4utq3T6/Xk5OS4PFZERISti7per2fevHm1F7ioNY4pWetkgH/TtWueCUYIlalyO+lbkZeXx+TJk1m6dCkWi4W8vDzbNt+SOkefMuMKGwwGQkJCXB7vzJkzdonem8OmiuorV5I2lb+VWFG3cSHqK48labPZzCOPPMKcOXPo3LkzRUVFHD16lMLCQvbv30/37t0Ba7XIyZMnadGiBTk5OYSGhro8ZlBQkAxVWg+Y9u6Ffv1sy1onf2xNbj4HQtRnHkvSK1euZM+ePSQlJZGUlMQzzzzDCy+8QEJCAjqdjk8//RSAuXPnMnXqVCwWC2+88YanwhNeVLB1q12SDrhxo9w+Jo3GkyEJoRoeS9KTJ09m8uTJ5daXDtZUKioqyjbRrWgYHKsytE5GvKtoDkQh6ivpzCK8znT77XbLWidtok0VzIEoRH0ln3zhdQVlWvQAaJ2MUW4KDvZUOEKoiiRp4XWmzp3tlrVNmpTfR0rSooGST77wOscErHUoWQMUFBZ6KhwhVEWStPC6gsOH7ZYDygysVEraSYuGSpK08DqTQx20//Xr5feRJniigZIkLbzO1KyZ3XKAk/pnk1broWiEUBdJ0sLrClq0sFvWNm5cbh+Tvz+KpwISQkUkSQuvc6zKcNYET/HxofxaIW7N+R492Ornx7GQEJR167wdjlMeHWBJCGccexO6qtgwudkmRFWdDQhgaWEhFoCrV7ln9Ghi162DUaO8HZqdOl2Sjo6OJioqioULF3o7FHEL8h16GAY6GbsDZJ5DUYMiI9lWmqBLbAcsTz3lrYhcqtMl6X379skoePVAvsOA/q6G9zd99x3ce2/tByTqPeXYMTIc1hUA2dnZtPFGQG7U6ZK0qB+uOwzyrzt3zul+pv37PRGOaAAMgLPuUZc8HUglSJIWXpfvUCetc3LjEMDkpNWHENVx0c95JYIkaSGcyG/a1G5Z17y50/1MTrqLC1EdBhedowwqnJVekrTwunyHXxhdhmNtoZWUpEVNyXMyZjmAQYUDeakvItHgOFZ3BLqq7jh71hPhiPpu/XqMLr6V5TmZBNnbJEkLr1JwUifdsqXTfaVruKgRW7aQ52KTUVEofuUVj4ZTEUnSwqsKsPYmLEuXn+90X1OHDh6ISNR7jRphdLFJAa59/LEno6mQJGnhVc7Ssc5FvaDpwoXaDaY6Ph8N2zWwTAPtXbXwFq4UvzILY7cYlJmzPHfSw4ddlqQBjC46U3mLR5O0wWCgX79+NGnShCNHjgCQnJxMTEwMiYmJnC2pczx27BiDBg0iJiaGzZs3ezJE4WFOk/SvvzrdV3VjSn80BR5ZB3HAo8DHBRAY6O2o6ozsx2fzzlsFvH3kPj6al0Ne8G0eOW/BuXNO20iXylNZBzmPJulGjRqxYcMGxo4dC4DZbObtt98mJSWFN998k6SkJABmzpzJ4sWL+eabb3jttdc8GaLwMGdJOvD8eaf7FrRrV7vBVFWj5fbL9wCtpPN6ZSiKwtrFOVzD2vzyAq3YlNsPIiNr/dzGy5fdb1fZzUOPJml/f3/CwsJsy+np6URGRqLVaomNjeXQoUMAZGVl0alTJ4KCgggJCeFyBS+qqLsck7RfURF+Lu68m65cqf2AqiLWSZzRng+jLrrYYQAXsb9BfJRICo/9t9bPfa2Cz5HRResib/FqnfTVq1ftxt6wlLRdLC7zS6rX68nJyXH6/Ly8PLtHQUFBtWMpooi5K7sybJMP787wQxk3ttrHquuKiyEp/nt6+fzMUz7/wjRgcK2dy3EOlkYubhoCXHdoBeJdZnD27byf1EtXRsbp8n/givElS1P7I2dcq6BTlFFl1WpeTdLNmjUjL+9mFb5vSW8fnzI3jgwGAyEhIU6fHxERgV6vtz3mzZtX7VjmrOnG7Im/8O1Qhel/s/Bhiy9h3LhqH68uW3T/el7bdg+pSg/+qTzBY7sfh9jYWjnX9V277JYbX78OLnocGiMiaiWG6skCZ53TIpx3khBlzJpFlothjDyRpK9X0CnKaFJXlZVXk3SnTp04evQohYWF7Ny5k+7duwMQHh7OyZMnMRqN5OTkEBoa6vT5Z86cwWAw2B6vVLN9o4LCsn7H7dZ98j/A+vXVOl5d968f7Ju6fc4j5O78pVbOZdi2zW5ZbzDA6NHO91XTXffZU5yvb6Our8qqNH8+OTgveF0udr6+Jl2roKRsVNl8mh4fqnT48OGkpqZy/PhxnnrqKV544QUSEhLQ6XR8+umnAMydO5epU6disVh44403XB4rKCioRoYqzSSTc23t1x3oC4UUNrhB5nNy4GdzV7t1xfjyI7GMqIXz5TVqZLccZLFAfLzzfdX0y3Nmp/P1t2H9Bpac7NFw6hSz2WWSvkJza+GoFgfev+ZQ5+wLXAGalSwbnUyE7E0eT9IbN24st27ChAl2y1FRUWzfvt1TIXH6x8/Aybf5k3doqP17zepy5KMdwMBy63/hrtpJ0g4dVIJ0Oih03kAqT01jd7R1MeNiG+C7DR4Npa4paBzC9etNnG67TCgsXlyrSfp6mfseucBCNCjA/6LQDDBZLBTduIG/SppTSmcWIHOT8x5Gv3ZWUcnNQ07+5Uun649q7qqV8+UZDHbLQRcvutzX0MT5L7ZXtHPxq+MDqhs1XmVyAsJdbsunMTfO1m4rnmtlCgH5QBEKZhS7lkbGrKxajaEqJEkDGYHOE8OZNg1vfuqTBuf1/8eUO2vlfAaHZb2b5lGqKkmHu+kOcac6SmBqlXPNvhIxt+RfqSuXavfma0WtO0CStOpkhju/2XNGZX0nPOEkHZ2uP037Wulo4Jh4g9zUB6oqSbd1sy1MRTc4VShHc7P1Ti65LGAhC1hoS9Q5F2qvM4miKFyvTJJ2MTuQN0iSBjJbOS8VnXX1lbYec5WkswnnxjHn4zzfiqokaYPDTUavctcasK1MTuBOjs/Nb2v55GOmCDNF5JdUOFwx62vt3KYVK6jMu2P87rtai6GqGl4WciLDRYn5TLilwbWVPknZG3mZJY+SJU3Nj61QlSRtbNqUYlU0izRBmJuqsNvUN7uHmlwtcP/HNkfbqtbOXdkR7oxOGjh4S4NP0goKma6SdASgojertuXmQg6lX0UzgS4lD2uizvCt+aFCcx3uoAddcj/L3LUVK2o8hqqrYPIBF9/MBDBrFjnFzdzukqPUXlvpaw6Dd7nqw2pU0RAEDT5JX+EK+S6qOs+1geIG9Aqd+r8FZZYuAzdKHtaxU04rNV+SvtSsmd1yaAV1gTlq6A328Rz326W6w6Wiz1ZhxH3fhiuFTVCU2rlpf+3qVbtlV2V6o4vptbyhAaUg5zLLfJ13VKSF7JYNp4XHqS9/crv9tP8dNX7OSw5dwMOys+2W/Rw6Hria5dmjDqx1v91VG2rB1esVdw8rUALIv+x6DJdbYXAY4a6Riw5Sxlr6I1EdDT5JZ+D+Zpirm4r10al859NWlTpd5K5JQ9UVYK1nLivMobojxKHkc0ENHQxCHUrz6Q7b9Qps/Nxj4dQlObn2KaeR08FqISfd+aBqtyrPIfk2cZOka6s0X1UNPkln/udD99tvazgdWk4r7tscnra0rdHxTJzVPoc5jNPS3GEExIvBwTV2/mpzrPVx1kN82VxPRFLn5Jjtqzr05VrKW11Jr506Ycck3dhFki4CCgzOY/O0Op2ko6OjiYqKYuHChdU+RsZZF2MwlG5vp46/pp5wyuK+zvkkHWFuzSUfxyTtY7EQ7FByDnZI0hfUkKQ7Onxlz2oKjv2hfGRmc2fKtpEGCHIx2+CVf66plfOXK0m7qT5TS4eWOp2k9+3bR1paGs8++2y1j5HRyv0Y1Jlt1HMDobaddjpA8k0XaUnesZr74J7fvdtuOfTyZXxOn7Zb19yhNHOhdesaO3+1RdjfGMw8FUDmTw6/7A2wt2plXNXYt9wI8nPe5DLnwKlaOb/BoSNLkzIz0DuOBG5cvbpWYqiqOp2ka4Jj87uW2Q7bI4BZHpwk00sUpXKtN9ILa64bZsbWrXbLt2VkwG32MYQ4JOms0iR9CdgEbmcUrRUKtLp5MzMzEzovvkLnURYyy96DbqGCVigq5Nj8LkhxXqWQU9jU6fpbUfTii+VqwBuXGcfDsaVH3tKlNR5DdTT4JO1YnTFwt31HhIzbsI7KVc+dPw/XqXgAo1+LO9XYOU87zAre/vTpcn8QWztUf5zo2NGanFuY4D5AfwX8+9dYTBW7QtmX6fJlMBUrmIoU7GZ5a6OuefLUoLBzN3Kxr67S+7qo7lCCa/zG3ZXPPiu3ruyfAseWuMYK2ux7SoNO0kaMXAq1//rTdbOubCc7a5J2SBT1UdrH9nXzOhd33dPN7WvsnKcdJpa9LSOj3BCV7RzqBdNvuw1l2AFufjltDuZ5oPXUyN/2k0O4nHY6vOFUk1XWxV9z7ZY1FBMc63x0xSK0XMu+VqPnv+IwNVsQ4D/w5rC8jRzqp42VGOPDExp0kj7GMbtlzSmY9+EN6KyxJercYLjorgtwPZG22D5J347zOsHjNViSPtytm93yHf8tPwlphEPnFn1uYzRKH4e9hkCRh4YxXfS8/bKrZva30SCqyarigsa+u3col/H7219sy37Yf/uo6WZ4V4z2pfZQgPnzbcuNHUruV9TQcYoGnqTTSLNbjjgIhcXFYFJKO9kBcKi7hwPzgl+y7e+6d8D5rM0H6VUjzfDygeOdO9ut63nkSLn9WgKNyoznMbxML/3Mkn9W8Z4ZZ6XwZ/vl0y72uw348K+1HEzd4jh/YUufy3bLTR1aemT/7HCD6BZddOhF6DiTZnOHknN2cbEq2ko36CR9iEN2yx1czCZ/6K76/9V1d6F96bSjiyR9jC4Ypjx3y+fbOnw4xb436/99LBa6HjpUbj8fk4m7y7QCeeA/1v8zyaQzXehMl5JEfY9n5qTs7lDXXPYLR9lNvkDXejbf4axZ0Ls3mf/7v2Rmuu6p60pmsX1nqPBi+2aKYQ6NMjNn/bPqMbrhOOBAK4c20s0dEnI+6hiytEEn6R3pn9gt3+nYc6zE7mh11E3VFoMBDmM/r2EPUm0/+3Ez2Sj4sDu3y62dMDaWpVPsJ3IdsGuX844FsbEMKpmsNsAE93xvXX2Zy5i4gYkbXOYyMBoKazspFkE/h1VlC/8nHbZ58n5mbZo1i+KWrVj9568YezCM2xd+RMfbOnBYo6v0t5e8ByZymTC7de0dvoaEY3//4bSxOcWWmvndM2ZllZlWwKqNzr7RXRDlB1zK6O/9N7FuJ2m9HjQa8Kn6ZVya+wIHbre/IdjroPN9Nw2FIv/62/Pw60dXoZT5KGgpIJKjtuUuDnX363iwaidYvx6mT4fmzTHo9bzywAOscpjXcvRXX4Gz9u7PPstDX30FwIQvoLHLIR3aAPdWLa6qmt6zfBOAw2V+dpxQfRR1fsZ5ZeYs/vPn1fS8+BDjGM+XtKcYM2YsDGAif1p9ketjxld4nF822A+/EICJVp3tx41u7ZCk82nMqc0101762O9+53B+CHMYx0Pj749jA9Nfq9qhpXlza07y86uxexJ1O0mXUhTrC+Ps4eNjnVEkIAB8fcHXl0sxnZgd+h6WMjdzA/Oh22Hnhzc0gz/NBmNTF+fQaG7OWlL6Jmk0EOtkdttSbdve3M/X1/p/2zJfB0u+Wlb4RpcmQGfJIDLS+mFp2/bmPmX3X7+eS5ow/rbefmqsRDYTUKb0PAj79szL+S3HNF2sMWu19nGWxG2aM4cMYO+995L8738zu3lz7luxghYXL/LWK6/YHa9pXh6PLV3qsjdjt0NHefGtVOa/6P6luKFbQnarh8kPbGyNrWkNtLUtfZ9a+8GD9vcw2AcUlCmN7XB47iBg6RgoHWY+NhYaNXL/uXD3fnrIjRtF7Nt3jnf94+g37xIj+Q2HKT+uy3W0vMoQOn51G7NmbWbLllPkfLbmZvwlnwXjtFnsYoDdc7twDJ9j9q9nENdoxXm7dSn3zaNA42pA0UqIjeWGRsPeb7+1W90J8HGcRGLyZBwniUsDLjpMllxO2d/l0h6yFgv8+c81kqg1ihpqxp146aWX2LlzJ+3bt2fJkiX4+/vbtuXl5aHXW3v932gB/3oCCgKsj0Kt85+vN4YLLSG7FVxxMo3fbz6D//stlNbM9vwXpD5efr/W56DlBQjJAb0BdCbQFkJAgfX/0p+D8uCl+eWf7wnvMY3LhGJCRwEBTv/PIYRztCGbVhRj3zZ8OZOIDN1En5KGv1/RhtFOxlBuxXnachY9BppwjYu/DePozChutAqkINix/5ZrHz/2GI+VaYv+008/0aeP9Z04cM8Ben/fAWhm95yf+Ik+Je/WAQ7Qm962bZebw48xZj566iIdTq1hwXMu6tCnAndgbc3n7BEItMBaSHfWG/0F+Om2t+nzhz9Y4whvTu9frpTftwBrK5DLWDvfGIFrWJvvFZX8X/rzP4EzzsO9FQdpxSruwoQfN/DnRpn/S9flouMcTcnF3SBWWSVBQh+eJJjW+IHdQ4eFRlgIRCEQHwLwQ4P9N9HrHCF61et07Fhw871u3x7z6TD+n8O89L6YaUYuTbiGlkL8KcKfIq5wHAsFJdPIWvChkBAuYcb6UhZhfanzoNxsLI8AnWfO5KeHH755/gMHuLNPH94FHO9C6bE2jw8oefhjve3g+IihfM/FvJLnGwwGgoLcD9PqjArGfSzv559/5ty5c2zfvp25c+eyevVqJk6c6HTfy6Hw6p9u7Xy+ZnhlHhRqtVDSA+nJT/z4/ePlOyRktbE+KhJ6yXtJ+h2mk0H7aj33To4znlUcnjkfSpJPhO4Kk03L+DeP2u2bTTjZlJn5WQtUYRpEH4uFv/7xjzy2ZInrDkPfg2OCBtg2kPIl1xKhV+DBr/14/Y3WhFwt/1ybKUBC5eO1cx74GNgWf3PdH2fBW3+AvzjsG4C16FaZ1ovrqZUknUYYbxF3y8fpwFXbLeU+gPNO+qUpy7mTwL/pyuObTvLMM2W+PX75Jb369GcXMXadXiz4cYVQruBYukoBh5rm81SsFVhLzHPnwk/2w/M28fOjj9nMXofnGCg/abIz4VhrxGpyelRVVnfs3LmToUOHAjBs2DB+/PFHl/tqb3EkUU0xfPh76PoLsGuXbX3/D/bw2hvVP26A+yFBapWO6rXvbEoeK5mIP2aIL5N8fvMb3mE6XcrUUztVyRZT2oICxiYnk9qzJ3945x2IiXG9s5NbAf+9/Qzz/3jzVybdyTDXuXpr08nQy5fLbyxV3flir2Mtilkcykzx8fA2cCuT+Tj5PJ+ncsnHnUBurQfkbeTyb9bwBcm3dBwD8FXJz/7+5RO5f5uWPMxqtFTmF8i/4l0c6IAxWOufnXrxRRKxJvKqygUG4EN3fNyMUl91qkzSV69etX0t0Ov15OS4btRe3WTY6DoM3wDb4+DJfwFduhAaGopOF4hOF0hoaChvvOnD1w9AwhZrabsqbvWPx60IqNQHvOz+Jsawhv30pTcHYexY+9dizhyak8NOYniWBQS5KlM4S9ImIAMSv/2eP86fz6px47gUFkby+PF0O3IEunQBhz/CdufuU7b09BPwFB1O3c62h4ei9fHFz8+fT/+Qyt//kE36HTe/1O4YCMW+5cenLhdbVVwDVgF9gW3AjRvlPjN8uQ4eBP4A/OruYC44+dx8VPK4FYFl7jFUVmcuM5WDbOAzTvA+v9X9SouMDHS+/ujQ0MjlvCblFQGpWCtKSltD+/v7lH/9zp6lLed4nI/pxK9o3E4bW/mKAA3QGXgCaOHvb/vGXO78c+ei7dKFqUA01i+H3qbK6o5mzZqRl2cdOcdgMBAS4nzOMwVoaoQHv7pZJxxQ4Pxnnclal9wqG8LPQ+fj4Fda8dSmDRw9Sjvg+HFrS4Z27drByy/zwJ//zAMb4IYOTtwBp26Hq8GQEwIG/c16b8e68FA3BbjaNpRNdOEYOkwEUGD3f+nPTTHShnOEk0Ukx2jUJgSKimDQWEhOLv9aKArBGg0LeI53eYFjdCGD2zhPONdowjWacCU9hKOPRKLNLkSbXYh/dhEYrF9XP2OStYQOMHMm3LgBCQnluoGXns927gPtIAeY3Qeyb341vcNcMtb+Pz+i3WO/A40GBTA07cqNwPvocNLMulH/5c5f3WTK7VhL0yY3j1ysDWzPAf8FW65bt658rO3aQbt2JYn6QXgHa2+c27F+/9VjHSyiKdYKTv+Sh7bk4Q84GYHgKddXUGkR5PEbDqHDTCBmAikiEHPJsvXnJhTSGiOtMdKWPJqU/YvRpg2cPWv9XPz3BAwZwp6T5/AlEz+K8MOMDxau4ccNfLmOL/n4YsZMIUUUUcjt+DIFX4qeeJqiomL6929r9/q1bNmS119/nVmrVxM2diy/YQUmdGTTEiNNyadRSW209XGaRphogwZfNPgRgIWOZODPzfrxxiUve3NAt25d+WEHHN8/gKNHCRg3juGrVzMU61heV7F+HApKHkVY661LH8Ul/08rucPjrLqjoKCaJUpFhQ4ePKhMnjxZURRFmTt3rrJixQq77QaDQQEUg8GgKDNnKkrv3orSpo2iBAZa/2/Z0vpo00ZRNBrro0kTRfHzs/5sbQ+iKOHhirJunftgZs5UzC1aKFmgmP7v/6zHAUXx9VWUFi0UZezYmzHodNZtfn7WY+v1ihIScvN8Op11P3//m8v9+ilKhw7W/Xx8FCU42HpMvV5RGjWyngesz/Hzu3msso+QEEXRam8u+/oqSpcu1vhiYhRl+nTrdZbG0qSJokyfrhT166dkg1L44IPVe6NiYqyvueP5/f2t19Svn3Wf0tfpVrVpYz1+QICijBrl+r2LiVGKfX0VMyjFpder11tfX2evn6uHRmN9Tni49T0eNcr9eR2NHWu99tL3sOxnz92jbJy+vtbPiY+PojRurJg7dlTSQTFX5lgajfW98PGp/LXrdNbPXbNm1muu6PpK9y39vSh9/zUa6+e1cWPrZyMmpsKXy+73et066+e2S5ebn3vHa/bzs5639NoCAqyfEb3e+nPZ17C6St+70tcGrNdU+lks/Xy1aGFdV/azXvKa5FrLk8qZM2eqFYIqk7SiKMqMGTOUgQMHKr/5zW+UgoICu212b2YVLViwoMrPqe75qnOu6j7Pk9dV3fN58vWoC+9ZdZ7n6fesus+T9+ymM2fO1M8k7c6tfFAjIyM9dr7qnKu6z/PkdVX3fJ58PerCe1ad53n6Pavu8+Q9u+lWk7Qq66QropQ07S6tt64Ki8VS5eeV7l/V51XnXNV9nievq7rn8+TrURfes+o8z9PvWXWfJ+/ZTcaS0fdK81ZVqbYziztnz54lIiLC22EIIUSlnTx5kg4V9V50ok4m6eLiYrKysmjatCkaF7P9CiGEGiiKgtFopHXr1vhUY5yhOpmkhRCioVBlZxYhhBBWkqSFEELFJEk7sXfvXgYMGMCgQYOYOHEiRUVFJCcnExMTQ2JiImfPlh8Rrq64cOECMTExxMfHM2TIEM6fP8+OHTuIiYlh4MCBHD7sYrzWOmLlypWEhVkHl68v79np06cJCwsjISGBhIQELl26VG+uLSUlhcTERAYPHszatWvrzWdx165dtvfrzjvvZPr06dW/tmo13KvnsrKylPz8fEVRFOXll19WkpOTlbvvvlspKChQduzYoTz55JNejrD6zGazYrFYFEVRlE8++URJSkpSBg0apOTk5CgZGRnK/fff7+UIq89sNitjxoxRevXqpRQVFdWb9+zUqVPKww8/bFuuL9eWn5+vPPDAA3ad1erLZ7GsKVOmKCkpKdW+NilJOxEeHk5goHVMXa1Wy/Hjx4mMjESr1RIbG8shJ3Px1RW+vr62O8xGo5GOHTvi6+tLcHAw7dq1czuYldqtXLmScePG4ePjQ3p6er15zwB+/PFH4uLimDlzZr25tl27dhEYGMjIkSMZM2YM58+frzefxVKFhYXs3buXvn37VvvaJEm7kZGRwaZNmxg4cKDdYN0WS92emDY1NZX+/fuzYMECYmJi7K7Nz8+PwkIvDuFXTRaLhVWrVjGhZFqusiMplm6vq8LDwzlx4gTbtm3j4sWLrFmzpl5c24ULFzhx4gRff/01TzzxBHPmzKkXn8Wyvv/+exITE8sN+F+Va5Mk7UJeXh6TJ09m6dKlhIWF2fUw8vV1PaB5XdCzZ0/27NlDUlISc+fOtbs2s9mMVquGARqrZvny5YwfP972LaHsSIpQt9+zgIAAGjdujEaj4aGHHuLnn3+uF9fWrFkzYmNj0Wq1JCYmcvDgwXrxWSwrOTmZcePGlfs8VuXaJEk7YTabeeSRR5gzZw6dO3emU6dOHD16lMLCQnbu3En37t29HWK1lf3rrdfradKkCWazmdzcXM6cOeNyWFi1S0tLY9myZQwbNoz09HQ++OCDevOelXYrBti+fTsjRoyoF9cWHR3N0aNHURSF1NRUoqKi6sVnsVRRURH79u1j4MCBNGrUqPrXVnvV5XXXsmXLlJCQECU+Pl6Jj49XPv/8c+Xzzz9XBgwYoAwePFjJzMz0dojVtmfPHiUuLk5JSEhQhg0bpmRlZSlbt25VBgwYoMTExCipqaneDvGW9enTR1EUpd68Zxs3blR69+6tDBw4UJk8ebJSVFRUb65twYIFSlxcnDJo0CDlxIkT9eqzuHHjRuW5556zLVf32qTHoRBCqJhUdwghhIpJkhZCCBWTJC2EEComSVoIIVRMkrQQQqiYJGkhhFAxSdJCCKFikqSFEELFJEmLeuXMmTMkJCQQFRVF9+7dSU5OdrpfRkYGw4cPp1OnTtxxxx3MmTPH5WzO165d45577kFRFObPn0+fPn3YvHkzYB3IacSIEbV2PUJIkhb1ip+fH++++y5paWls2rSJF154gevXr9vtoygKo0ePZtKkSaSnp3PkyBF++ukn3n//fafH/Pjjjxk/fjzXrl3j0KFD7Nixg48++giA4OBg2rRpw65du2r92kTDJEla1Cvh4eH07NkTgFatWhEaGlpu7N7NmzcTFBTEpEmTANDpdLz//vv89a9/dXrMFStW8OCDD9pK2j4+Pnaz1I8aNYqVK1fWwtUIIUla1GMHDhzAYrEQERFhtz4tLY3evXtz33330bNnT3r27MnGjRu5fv263XCSAAUFBVy4cIGWLVsSFBRE586dGTBgAI8//rhtn969e7Nz506PXJNoePy8HYAQtSEnJ4dHH32Uf/3rXy73+fbbb+2WZ8+eXW6fK1euEBwcbFt+9dVXefXVV+32CQsL4/z587cYsRDOSUla1DsFBQWMHj2al19+mZiYmHLbo6KiOHjwoN26U6dO0bhxY7vZM8BaFWIymdyez2Qy2aZbE6KmSZIW9YqiKEydOpUhQ4YwefJkp/skJiZiNBpt9cgmk4lp06bxxz/+sdy+ISEh3LhxA7PZ7PKcJ06cIDIysmYuQAgHkqRFvfLjjz/yxRdf8NVXX9nqmw8fPmy3j0ajYe3atSxfvpxOnTpx11130atXL6ZNm+b0mPHx8ezZs8flObdu3cr9999fo9chRCkZ9F+ICuzevZulS5fyj3/8w+n2xMREVq9ebVd3LURNkRuHQlTg7rvvJi0tDUVR7JregbUzy3PPPScJWtQaKUkLIYSKSZ20EEKomCRpIYRQMUnSQgihYpKkhRBCxSRJCyGEikmSFkIIFZMkLYQQKiZJWgghVEyStBBCqJgkaSGEULH/DwGQ/BFvp1bIAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 350x218.75 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "from lmfit.models import PseudoVoigtModel\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "from skimage import io\n",
    "from skimage.restoration import rolling_ball\n",
    "from scipy.signal import savgol_filter\n",
    "from matplotlib.ticker import AutoMinorLocator\n",
    "\n",
    "file = 'PXRD_cpy_131_repro_calc.xy'\n",
    "\n",
    "data = pd.read_csv(file, sep= ' ', header = None)\n",
    "counts = data[1]\n",
    "counts_smoothed = savgol_filter(counts, 17,1)\n",
    "two_theta = data[0].to_numpy()\n",
    "background = rolling_ball(counts_smoothed, radius=100)\n",
    "counts_smoothed-=background\n",
    "n_peaks = 8\n",
    "\n",
    "mask = pattern.y > 1\n",
    "\n",
    "\n",
    "\n",
    "peak_positions  = pattern.x[mask][:n_peaks]\n",
    "peak_widths     = np.ones_like(peak_positions)*2\n",
    "peak_amplitudes = pattern.y[mask][:n_peaks]\n",
    "peak_start      = pattern.x[mask][:n_peaks]-3\n",
    "peak_stop       = pattern.x[mask][:n_peaks]+3\n",
    "peak_start[3]   = pattern.x[mask][3]-1\n",
    "peak_stop[3]    = pattern.x[mask][3]+1\n",
    "\n",
    "\n",
    "lambd = 1.54e-10\n",
    "\n",
    "colors_all = [\n",
    "        (0/255,0/255,0/255), (255/255,0/255,0/255), (0/255,255/255,0/255), (0/255,0/255,255/255),\n",
    "        (0/255,255/255,255/255), (255/255,0/255,255/255), (255/255,255/255,0/255),\n",
    "        (0/255,0/255,128/255), (128/255,0/255,128/255), (128/255,0/255,0/255), (0/255,128/255,0/255), (0/255,128/255,128/255),\n",
    "        (0/255,0/255,160/255), (255/255,128/255,0/255), (128/255,0/255,255/255), (255/255,0/255,128/255),\n",
    "        (255/255,255/255,255/255), (0/255,0/255,0/255), (192/255,192/255,192/255), (128/255,128/255,128/255)]\n",
    "\n",
    "\n",
    "\n",
    "    \n",
    "def make_model(num):\n",
    "    #pref = \"f{0}_\".format(num)\n",
    "    model = PseudoVoigtModel()#prefix = pref)\n",
    "    model.set_param_hint('amplitude', value=peak_amplitudes[num], min=0, max=np.inf, vary = True)\n",
    "    model.set_param_hint('center',    value=peak_positions[num], min=peak_positions[num]-1, max=peak_positions[num]+1, vary = False)\n",
    "    model.set_param_hint('fraction',  value=0.1,  vary = True)\n",
    "    model.set_param_hint('sigma',     value=peak_widths[num], min=0, max=5., vary = True)\n",
    "\n",
    "    return model\n",
    "\n",
    "def save_fits(x_exp, y_exp, x_fit, ys_fit, names, params, errs, fname):\n",
    "    fit_labels    = np.zeros(n_peaks+1, dtype = 'U32')\n",
    "    fit_labels[0] = '2Theta'\n",
    "    for i_peak in range(1, n_peaks+1): \n",
    "        fit_labels[i_peak] = 'counts fit {}'.format(i_peak)\n",
    "    file = open(fname, 'w+')\n",
    "    data = np.zeros((n_peaks + 1, len(ys_fit[0,:])))\n",
    "    data[0,:] = x_fit; data[1:,:] =  ys_fit; dataframe = pd.DataFrame(np.transpose(data))\n",
    "    dataframe.to_csv(path_or_buf=file, sep = '\\t', header = fit_labels)\n",
    "    file.close()\n",
    "    file = open(fname, 'a+')\n",
    "    data = np.zeros((2, len(y_exp)))\n",
    "    data[0,:] = x_exp; data[1,:] =  y_exp; dataframe = pd.DataFrame(np.transpose(data))\n",
    "    dataframe.to_csv(path_or_buf=file, sep = '\\t', header = ['angle', 'intensity'])\n",
    "    file.close()\n",
    "    file = open(fname, 'a+')\n",
    "    out_data_pars = np.zeros((len(params), 3), dtype = 'U32')\n",
    "    out_data_pars[:,0] = names\n",
    "    out_data_pars[:,1] = params\n",
    "    out_data_pars[:,2] = errs\n",
    "    dataframe = pd.DataFrame(out_data_pars)\n",
    "    dataframe.to_csv(path_or_buf=fname, sep = '\\t')\n",
    "\n",
    "def find_index(array, value): \n",
    "    index = np.argmin(np.abs(array-value))\n",
    "    return index\n",
    "\n",
    "def calc_size(FWHM_deg, two_theta_deg, lambd = 1.54e-10, K=0.9):\n",
    "    \"\"\"\n",
    "    Scherrer crystallite size.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    FWHM_deg : float\n",
    "        Peak FWHM in degrees (2theta)\n",
    "    two_theta_deg : float\n",
    "        Peak position in degrees (2θ).\n",
    "    lambd : float\n",
    "        X-ray wavelength (same length unit as D).\n",
    "    K : float\n",
    "        Scherrer shape factor (default 0.9).\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    D : float\n",
    "        Crystallite size.\n",
    "    \"\"\"\n",
    "    beta = np.deg2rad(FWHM_deg)\n",
    "    print(beta)\n",
    "    theta = np.deg2rad(two_theta_deg / 2.0)\n",
    "    return K * lambd / (beta * np.cos(theta))\n",
    "\n",
    "\n",
    "\n",
    "names  = np.array([])\n",
    "params = np.array([])\n",
    "errs   = np.array([])\n",
    "fits   = np.ones((n_peaks, len(two_theta)*100))*np.nan \n",
    "angle_fit = np.linspace(start=two_theta[0], stop=two_theta[-1], num = len(two_theta)*100)\n",
    "sizes  = np.array([])\n",
    "mod = None\n",
    "for i_peak in range(n_peaks):\n",
    "    start = find_index(array = two_theta, value = peak_start[i_peak])\n",
    "    stop = find_index(array = two_theta, value = peak_stop[i_peak])\n",
    "    mod = make_model(i_peak)\n",
    "    out=mod.fit(counts_smoothed[start:stop], x=two_theta[start:stop], method ='trf')\n",
    "    for name, param in out.params.items():\n",
    "        print(f'{name:7s}\\t {param.value:11.5f}\\t {param.stderr:11.5f}\\r')\n",
    "        names  = np.append(names, name)\n",
    "        params = np.append(params, param.value)\n",
    "        errs   = np.append(errs, param.stderr)\n",
    "    fits[i_peak, :] = out.eval(x=angle_fit)#mod.eval(x = angle_fit, amplitude = params[-6], center = params[-5], sigma = params[-4], fraction = params[-3] )\n",
    "    sizes = np.append(sizes, calc_size(FWHM_deg = out.params['fwhm'].value, two_theta_deg = out.params['center'].value))\n",
    "\n",
    "\n",
    "\n",
    "save_fits(x_exp = two_theta, y_exp = counts_smoothed, x_fit= angle_fit, ys_fit = fits, names = names, params = params, errs = errs, fname = f'{file[:-3]}_result.txt')\n",
    "\n",
    "Labels = [entry[0]['hkl'] for entry in pattern.hkls]\n",
    "# Legend\n",
    "plt.rcParams['legend.fontsize'] = 6\n",
    "plt.rcParams['legend.frameon'] = False\n",
    "# Axes\n",
    "plt.rcParams['axes.labelsize'] = 7\n",
    "plt.rcParams['axes.titlesize'] = 7\n",
    "# Ticks\n",
    "plt.rcParams['xtick.labelsize'] = 6\n",
    "plt.rcParams['ytick.labelsize'] = 6\n",
    "fig, axs = plt.subplots(nrows=1, ncols=1, sharex = 'col', figsize = [3.5,3.5/1.6], constrained_layout=True)#, gridspec_kw = {'width_ratios' :[2, 1]})\n",
    "axs.tick_params(which = 'both', bottom=True, left=True, direction ='in')\n",
    "axs.plot(two_theta, counts_smoothed, color = colors_all[1], linestyle = '', marker = 'o', fillstyle = 'none', markersize = 2, label = '$\\mathrm{Co_3O_4}$ nanoparticles')\n",
    "axs.set_ylabel('counts')\n",
    "axs.set_xlabel('2 $\\Theta$ (°)')\n",
    "axs.set_ylim(top = 600)\n",
    "axs.set_xlim(left = 15, right = 70)\n",
    "for i_peak in range(n_peaks): \n",
    "    start = find_index(array = angle_fit, value = peak_start[i_peak])\n",
    "    stop = find_index(array = angle_fit, value = peak_stop[i_peak])\n",
    "    if i_peak == 0: \n",
    "        start = start - 2000\n",
    "    axs.plot(angle_fit[start:stop], fits[i_peak, start:stop], color = colors_all[i_peak+2], linewidth = 3)\n",
    "axs.yaxis.set_minor_locator(AutoMinorLocator(n = 2))\n",
    "axs.xaxis.set_minor_locator(AutoMinorLocator(n = 4))\n",
    "axs.bar(pattern.x, pattern.y*np.max(counts_smoothed)/200, width=0.2, color='navy', edgecolor='black', label = '$\\mathrm{Co_3O_4}$ COD ID 1526734', zorder=10)\n",
    "axs.legend()\n",
    "plt.savefig(f'{file[:-3]}_result.tif', dpi = 600, pil_kwargs={\"compression\": \"tiff_lzw\"}, bbox_inches='tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "id": "64e43895",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.17107616e-08, 1.03134901e-08, 1.04163331e-08, 1.58663527e-08,\n",
       "       1.11180458e-08, 6.97913183e-09, 9.46327806e-09, 1.01186744e-08])"
      ]
     },
     "execution_count": 97,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sizes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "id": "972ac497",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([7.94247624e+01, 1.90599847e+01, 3.43800212e-01, 4.18882737e-01,\n",
       "       6.87600424e-01, 9.38625753e+01, 8.91661482e+01, 3.13719356e+01,\n",
       "       3.99883125e-01, 3.14008022e-01, 7.99766251e-01, 9.41368190e+01,\n",
       "       4.04220818e+02, 3.69668946e+01, 4.01922634e-01, 3.32918278e-01,\n",
       "       8.03845268e-01, 4.21709513e+02, 2.65407921e+01, 3.86744779e+01,\n",
       "       2.65214524e-01, 1.00000000e+00, 5.30429048e-01, 3.18542001e+01,\n",
       "       9.95182574e+01, 4.49590655e+01, 3.86498520e-01, 3.25072804e-01,\n",
       "       7.72997040e-01, 1.08273130e+02, 2.36830719e+01, 5.58461034e+01,\n",
       "       6.43887385e-01, 1.00000000e+00, 1.28777477e+00, 1.17078795e+01,\n",
       "       1.17599458e+02, 5.95627179e+01, 4.83427129e-01, 2.26842928e-01,\n",
       "       9.66854258e-01, 1.05909625e+02, 1.54658473e+02, 6.54666032e+01,\n",
       "       4.66481812e-01, 1.84486733e-01, 9.32963624e-01, 1.46470677e+02])"
      ]
     },
     "execution_count": 98,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "id": "96f2c8b3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 19.059984718542005 14.172265569684319 [{'hkl': (1, 1, 1), 'multiplicity': 8}] 1.1710761554684098e-08\n",
      "1 31.371935610991773 31.14589407816431 [{'hkl': (2, 2, 0), 'multiplicity': 12}] 1.0313490140744915e-08\n",
      "2 36.96689464787737 100.0 [{'hkl': (3, 1, 1), 'multiplicity': 24}] 1.0416333083486817e-08\n",
      "3 38.674477930430854 10.885058640504548 [{'hkl': (2, 2, 2), 'multiplicity': 8}] 1.5866352690168263e-08\n",
      "4 44.95906549020098 18.912350668592204 [{'hkl': (4, 0, 0), 'multiplicity': 6}] 1.1118045773626324e-08\n",
      "5 55.84610335426648 9.556509501942594 [{'hkl': (4, 2, 2), 'multiplicity': 24}] 6.9791318336891024e-09\n",
      "6 59.56271785446117 37.38893280332379 [{'hkl': (5, 1, 1), 'multiplicity': 24}, {'hkl': (3, 3, 3), 'multiplicity': 8}] 9.463278061694116e-09\n",
      "7 65.4666031536717 42.22505800536844 [{'hkl': (4, 4, 0), 'multiplicity': 12}] 1.0118674421847211e-08\n"
     ]
    }
   ],
   "source": [
    "hkls_masked = [hkl for hkl, m in zip(pattern.hkls, mask) if m]\n",
    "for i in range(8):\n",
    "    print(i, pattern.x[mask][i], pattern.y[mask][i], hkls_masked[i], sizes[i])\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "id": "0f070f39",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(6.9791318336891024e-09, 1.5866352690168263e-08)"
      ]
     },
     "execution_count": 100,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sizes.min(), sizes.max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "id": "4beed617",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[6.9791318336891024e-09,\n",
       " 9.463278061694116e-09,\n",
       " 1.0118674421847211e-08,\n",
       " 1.0313490140744915e-08,\n",
       " 1.0416333083486817e-08,\n",
       " 1.1118045773626324e-08,\n",
       " 1.1710761554684098e-08,\n",
       " 1.5866352690168263e-08]"
      ]
     },
     "execution_count": 101,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sorted(sizes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "id": "c34ffdbc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0748258444992604e-08"
      ]
     },
     "execution_count": 102,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.average(sizes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "id": "3f159493",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.3401684653773138e-09"
      ]
     },
     "execution_count": 103,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.std(sizes)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "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.12.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
