Analysis_SOS.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "57560a32",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scattering_2D import get_reflexes, plot_result_gixd, export_results_csv\n",
    "from ase import Atoms, Atom\n",
    "from ase.io import read\n",
    "import numpy as np\n",
    "import itertools"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "b73e54a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "reflections_exp = np.array([15.1, 16.7, 20.8, 25.2, 26.3, 28.1, 30.2])/10  # nm^-1 -> Angstrom^-1"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "bc15ebf6",
   "metadata": {},
   "source": [
    "If we assume that the first reflex is a [01] reflex or equivalent, we know that the lattice planes should be approximately 4.161 Å spaced:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "f56fada8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4.161049872304361\n"
     ]
    }
   ],
   "source": [
    "d = 2*np.pi/reflections_exp[0]\n",
    "print(d)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "758d7357",
   "metadata": {},
   "source": [
    "If we assume for now that the two lattice vectors have the same length and the lattice angle is between 90° and 120°, this means due to geometry the lengths of the lattice vectors should be between 4.161 and 4.805 Å."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "57e251e9",
   "metadata": {},
   "source": [
    "Rough screening with symmetry constraint of |a|=|b|:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "0bd0ddbf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# unit-cell parameters to screen\n",
    "alist = np.linspace(4, 5, 40)\n",
    "#blist = np.linspace(4, 5, 10)\n",
    "gamma = np.linspace(90, 120, 30)\n",
    "\n",
    "# only one basis atom at 0,0,0\n",
    "fractional_coords = [[0,0,0]]\n",
    "n=0\n",
    "best = []\n",
    "for p in itertools.product(alist, gamma):\n",
    "    a, g = p\n",
    "    atoms_dummy = Atoms('C', scaled_positions=fractional_coords, cell=[a, a, 100, 90, 90, g])  # create dummy ase atoms object\n",
    "\n",
    "    result = get_reflexes(atoms_dummy, dummy=False)\n",
    "    \n",
    "    # check if result is close to experimental reflexes\n",
    "    matches = [np.min(np.abs(reflections_exp - Qxy)) for Qxy in result['Qxy_unique'] if np.min(np.abs(reflections_exp - Qxy)) < 0.02]\n",
    "    # sum up the errors between the reflexes and the reflexions as further criterion for a good match\n",
    "    summatches = [np.abs(r-result['Qxy_unique']) for r in reflections_exp]\n",
    "    summatches = np.sum(np.array(summatches))\n",
    "\n",
    "    best.append([a, a, g, summatches, len(matches)])\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "3b5f3044",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[  4.48717949   4.48717949 112.75862069 508.6003416    5.        ]\n",
      " [  4.51282051   4.51282051 112.75862069 503.96308376   5.        ]]\n",
      "[[  4.51282051   4.51282051 112.75862069 503.96308376   5.        ]]\n"
     ]
    }
   ],
   "source": [
    "best = np.array(best)\n",
    "bestmin = best[np.where(best[:,-1]==np.max(best[:,-1]))]  # checking which of the results has the highest number of matching reflexes\n",
    "bestmin2 = bestmin[np.where(bestmin[:,-2]==np.min(bestmin[:,-2]))]  # of these, which have the lowest sum of errors\n",
    "print(bestmin)\n",
    "print(bestmin2)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "e7c7c3f1",
   "metadata": {},
   "source": [
    "Finer screening based on results of coarser screening."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "8f7429cd",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "10000\n",
      "20000\n",
      "30000\n"
     ]
    }
   ],
   "source": [
    "# unit-cell parameters to screen\n",
    "alist = np.linspace(4.4, 4.6, 200)\n",
    "#blist = np.linspace(4.4, 4.6, 20)\n",
    "gamma = np.linspace(112, 113, 200)\n",
    "\n",
    "# only one basis atom at 0,0,0\n",
    "fractional_coords = [[0,0,0]]\n",
    "n=0\n",
    "best = []\n",
    "for p in itertools.product(alist, gamma):\n",
    "    a, g = p\n",
    "    atoms_dummy = Atoms('C', scaled_positions=fractional_coords, cell=[a, a, 100, 90, 90, g])  # create dummy ase atoms object\n",
    "    \n",
    "    if n%10000==0:\n",
    "        print(n)\n",
    "    n += 1\n",
    "\n",
    "    result = get_reflexes(atoms_dummy, dummy=False)\n",
    "    \n",
    "    # check if result is close to experimental reflexes\n",
    "    matches = [np.min(np.abs(reflections_exp - Qxy)) for Qxy in result['Qxy_unique'] if np.min(np.abs(reflections_exp - Qxy)) < 0.01]\n",
    "    # sum up the errors between the reflexes and the reflexions as further criterion for a good match\n",
    "    summatches = [np.abs(r-result['Qxy_unique']) for r in reflections_exp]\n",
    "    summatches = np.sum(np.array(summatches))\n",
    "\n",
    "    best.append([a, a, g, summatches, len(matches)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "37840445",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[  4.49447236   4.49447236 112.63819095 506.6503825    5.        ]\n",
      " [  4.49447236   4.49447236 112.64321608 506.67615781   5.        ]\n",
      " [  4.49447236   4.49447236 112.64824121 506.70194028   5.        ]\n",
      " ...\n",
      " [  4.52060302   4.52060302 112.86432161 503.1288       5.        ]\n",
      " [  4.52060302   4.52060302 112.86934673 503.15423358   5.        ]\n",
      " [  4.52060302   4.52060302 112.87437186 503.17967415   5.        ]]\n",
      "[[  4.52060302   4.52060302 112.75376884 502.57102451   5.        ]]\n"
     ]
    }
   ],
   "source": [
    "best = np.array(best)\n",
    "bestmin = best[np.where(best[:,-1]==np.max(best[:,-1]))]  # checking which of the results has the highest number of matching reflexes\n",
    "bestmin2 = bestmin[np.where(bestmin[:,-2]==np.min(bestmin[:,-2]))]  # of these, which have the lowest sum of errors\n",
    "print(bestmin)\n",
    "print(bestmin2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "cdfbe2a5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAArxUlEQVR4nO3de7xVZbX/8c9gy32DEnAwQNqQFxC5BIooXlAz0UTKn4pkdUiLlymdztFIPZal6flpapq3CMFUfmpoB4XTwSMpB0wUQQooAZWQ4qIIiCLKRWD8/lhr4Wa7LxOYl2fu/X2/XvvlXnPNNcdwrDn3w7w8z2PujoiISGgaZZ2AiIhIddRAiYhIkNRAiYhIkNRAiYhIkNRAiYhIkNRAiYhIkBJroMzsATN7x8z+WsP7ZmZ3mdkyM1tkZv2SykVERPInyTOoB4Ehtbx/JnBY8WcU8KsEcxERkZxJrIFy9+eBd2tZZRjwsBfMAQ4ys88mlY+IiOTLARnG7gSsrPR6VXHZW1VXNLNRFM6yaNmyZf/u3bunkmCDsHkNACve3QZAxWeaQnnHREOuWL26EKtTp0Tj5I3qEk2cdcqq5qF911nnM3/+/PXu3r7q8iwbKKtmWbXjLrn7OGAcwNFHH+2vvPJKknk1LLcXvoZLFg0CYELv2XDlmkRDXnLttYVYN92UaJy8UV2iibNOWdU8tO8663zM7O/VLc+ygVoFHFLpdWcg2b+MUqNvdlqWXqxhw1KLlSeqSzRx1imrmof2XYeWT4klOVismVUAv3f3o6p578vAaOAs4FjgLncfUNc2dQYVs9urOZG9UgMIi0h6zGy+ux9ddXliZ1Bm9hgwGGhnZquAnwCNAdx9LDCNQuO0DPgI+FZSuUjdVnxUDkBFi83Jxwrs+nsoVJdodA8qfqHlU5JYA+XuI+p434HLk4ove+dny/oAxXtQSce6775CrECuv4dCdYkmrjp9/PHHLFi0iLYHHsiWTZviSC2yDe+9B5B63JqklU+zZs3o3LkzjRs3jrR+lvegREQys2rVKg7v1o3mLVvStXPnVGOHdsaSRj7uzoYNG1i1ahVdu3aN9BkNdSQiDdLWrVtp3rIlZtU9UCxxMzPatm3L1q1bI39GDZSINFhqnNK1t/VWAyUiIkHSPSgB4DuHvJ5erAsuSC1Wnqgu0cRZp3Zt2sS2rX2Ne9NNN/Hoo49SVlZGo0aN+PWvf83999/PFVdcwZFHHrnfsSoqKnjllVdo165djev8Zvx4xowZs/v18ccfz4svvrjfsfdXov2gkqB+UDFTPyhpoJYsWUKPHj0yzeGll17iiiuuYObMmTRt2pT169ezfft2OnaMb7ixKA1UeXk5mzcn38UEqq97Tf2gdIlPAFi6uTVLN7dOJ9by5SxdvjyVWHmiukQTZ522bNvGlm2FcSh58aeFf7DF9fPiT+uM+9Zbb9GuXTuaNm0KQLt27ejYsSODBw+m9A/x8vJyrrrqKvr3788Xv/hF5s6dy+DBg+nWrRtTp04F4MEHH2T06NG7t3/22Wczc+bMT8X9yle+Qv/+/enZsyfjxo0D4Oqrr2bLli307tOHiy66aHdMKDx5N2bMGI466ih69erFpEmTAJg5cyaDBw/mvPPOo3v37lx00UUkcbKjBkoAuHV5L25d3iudWBMmcOuECanEyhPVJZo467R2/XrWrl8fy7b2Je6XvvQlVq5cyeGHH85ll13GrFmzPrXuhx9+yODBg5k/fz6tWrXiRz/6EX/4wx948sknue666/Yq7gMPPMD8+fN55ZVXuOuuu9iwYQM333wzzZo1Y+q0aTzyyCN7rD958mQWLFjAwoULefbZZxkzZgxvvVUYz/vPf/4zd955J4sXL2b58uXMnh1/H0o1UCIiGSkvL2f+/PmMGzeO9u3bM3z4cB588ME91mnSpAlDhhSm1uvVqxcnn3wyjRs3plevXqxYsWKv4t1111306dOHgQMHsnLlSt54441a13/hhRcYMWIEZWVldOjQgZNPPpl58+YBMGDAADp37kyjRo3o27fvXucShR6SEBHJUFlZGYMHD2bw4MH06tWLhx56aI/3GzduvPvx7EaNGu2+HNioUSN27NgBwAEHHMCuXbt2f6a6vkYzZ87k2Wef5aWXXqJFixYMHjy4zj5JtV22K+VR+n8o5RInnUGJiGTktdde2+MsZsGCBXzuc5/b6+1UVFSwYMECdu3axcqVK5k7d+6n1nn//fdp06YNLVq0YOnSpcyZM2f3ewc0bszHH3/8qc+cdNJJTJo0iZ07d7Ju3Tqef/55Bgyoc0zv2OgMSkQE4PifFn5StHnzZr73ve/x3nvvccABB3DooYcybtw4zjvvvL3azqBBg+jatSu9evXiqKOOol+/fp9aZ8iQIYwdO5bevXtzxBFHMHDgwN3vjfja1zjz9NM5dsCAPe5DffWrX+Wll16iT58+mBk///nPOfjgg1m6dOm+/0/vBT1m3tAVHzNfsKnQL6Nv642JP2a+oLhz99XMyHtQXaKJq05LlizhcxUVALRo3nx/09orH23ZkkncmqSZz948Zq4zKAGKDVNasfQHuFqqSzRx1imrBiKUhqkktHxKdA9KgMIZVOksKvFYS5fu/lewfEJ1iSbOOn20Zcvus4c0ZRW3JqHlU6IGSgC4e8WR3L1i/4dViRRr4kTunjgxlVh5orpEE2ed3nn3Xd55991YtpWHuDUJLZ8SNVAiIhIkNVAiIhIkNVAiIhIkNVAiIhkpKyujb9++u39uvvnmRONNnTo18RgzZ86MbaoOPWYuAIzp9pf0Yl1ySWqx8kR1iSbOOnWoZQqKJJXiNm/enAULFqQSc8eOHZxzzjmcc845NeYTh5kzZ1JeXs7xxx+/39vSGZQA0L18E93LN6UTq1s3unfrlkqsPFFdoomzTs2bNqV5pTHl0lJb3Pfff58jjjiC1157DYARI0Zw//33A4XBZa+88kr69evHaaedxrp16wD429/+xpAhQ+jfvz8nnnji7pEeRo4cyRVXXMEpp5zCVVddtce0HCNHjuS73/0up5xyCj179GDunDlcfPHF9OjRg5EjR+7OZ/r06Rx33HH069eP888/f/e8URUVFfzkJz+hX79+9OrVi6VLl7JixQrGjh3LHXfcQd++ffnjH/+4X3XSGZQAMGdjewAGtlmXfKyFCwux+vRJPFaeqC7RxFmnzR99BEB5ixZccu21n3r/S4MGMfyss9iybRujb7jhU++fc+qpDDvtNDZu2sQPbrllj/cm3HRTnXG3bNlC3759dy+/5pprGD58OPfccw8jR47k+9//Phs3buQ73/kOUJh6o1+/ftx+++3ccMMNXH/99dxzzz2MGjWKsWPHcthhh/Hyyy9z2WWXMWPGDABef/11nn32WcrKyj41UvrGjRuZMWMGk554gqFDhzJ79mzGjx/PMcccw4IFC+jcuTM33ngjzz77LC1btuSWW27hF7/4xe5pPtq1a8ef/vQn7rvvPm677TbGjx/PpZdeSnl5OT/4wQ/qqH7d1EAJAPevPBxIp4G6//HHC7H0h3gPqks0cdZp/cbCCCrlLVrs97b2JW5Nl/hOP/10nnjiCS6//HIWFhtkKIxgPnz4cAC+/vWvc+6557J582ZefPFFzj///N3rbStNwgicf/75lJWVVZvH0KFDMTM6du5M27Zt6dWrMCdcz549WbFiBatWrWLx4sUMGjQIgO3bt3Pcccft/vy5554LQP/+/Zk8efK+lKJWaqBERKj9jKd506a1vt+mdeta399bu3btYsmSJTRv3px3332Xzp07V7uembFr1y4OOuigGu9ltWzZssY4lafuaFLpkmNpKo+ysjJOP/10HnvssVo/r+k2REQaiDvuuIMePXrw2GOPcfHFF++eCmPXrl387ne/A+DRRx/lhBNOoHXr1nTt2pUnnngCKMzhVPmsa38MHDiQ2bNns2zZMgA++ugjXn/99Vo/06pVKz744INY4quBEhHJSOkeVOnn6quv5vXXX2f8+PHcfvvtnHjiiZx00knceOONQOFs6NVXX6V///7MmDFj972gRx55hAkTJtCnTx969uzJlClTYsmvffv2PPjgg4wYMYLevXszcODAOqfaGDp0KE8++WQsD0louo2GrjjdxiWLCteYJ/Senfh0G6Wb0XFeEqkPVJdo4qrTkiVLaN66NQAVnTrtd157Y8Xq1fsUt7y8fPdTdCHksy803YbstR8fGs8lgUixLrsstVh5orpEE2edPtu+fWzbykPcmoSWT4kaKAGgokX8/yqrMVbK/1rNC9Ulmjjr1LRJk9i2lUbcJM6eILs61EX3oASAWRs6MGtDh3RizZ3LrLlzU4mVJ6pLNHHWadPmzXzw4YexbGtvfPDhh5nErUla+eztLSWdQQkAD68+FICT265NPlbxBu7JAwYkHitPVJdo4qpTs2bNWLV6Nc1btqRVLY9iJ2HDe+8BpB63Jmnk4+5s2LCBZs2aRf6MGigRaZA6d+7Mfz/7LG0PPJCtMT0WHVWpQdiyKZ3hxeqSVj7NmjWrsU9XddRAiUiD1LhxY/7rhReA9J+cDO2JzdDyKdE9KBERCZIaKBERCVKiHXXNbAjwS6AMGO/uN1d5/0Dg/wFdKFxuvM3df1PbNtVRN2bFjrpvbyvcuDy46dbEO+q+XZwi4OBA+15kRXWJJs46ZVXz0L7rrPNJvaOumZUB9wKnA6uAeWY21d0XV1rtcmCxuw81s/bAa2b2iLtvTyovqd7BTbemFyuQgzI0qks0cdYpq5qH9l2Hlk9Jkpf4BgDL3H15scH5LTCsyjoOtDIzA8qBd4H4h8SVOj2zriPPrOuYTqwXXuCZ4s1p+YTqEk2cdcqq5qF916HlU5LkU3ydgJWVXq8Cjq2yzj3AVGAN0AoY7u67qm7IzEYBowC6dOmSSLIN3eNvdQXgjPZrko/19NOFWCeckHisPFFdoomzTlnVPLTvOrR8SpI8g7JqllW9uXEGsADoCPQF7jGz1p/6kPs4dz/a3Y9uH+ipqIiIxCvJBmoVcEil150pnClV9i1gshcsA94EuieYk4iI5ESSDdQ84DAz62pmTYALKVzOq+wfwGkAZtYBOAJYnmBOIiKSE4ndg3L3HWY2GniGwmPmD7j7q2Z2afH9scDPgAfN7C8ULgle5e7rk8pJRETyQxMWNnTFflAbPy4Mt9+m8fbE+0FtLI731ab1p243NmiqSzRx1imrmof2XWedjyYslFq1aZxe17NQDsrQqC7RxFmnrGoe2ncdWj4lGupIAJiy9hCmrD2k7hXjiPXcc0x57rlUYuWJ6hJNnHXKquahfdeh5VOiBkoAmLq2C1PXptPHbOqMGUydMSOVWHmiukQTZ52yqnlo33Vo+ZSogRIRkSCpgRIRkSCpgRIRkSCpgRIRkSCpH1RDV+wHtWVnGQDNy3Ym3g9qy7ZthVhNmyYaJ29Ul2jirFNWNQ/tu846H/WDklo1L9uZXqxADsrQqC7RxFmnrGoe2ncdWj4lusQnAExaU8GkNRXpxJo2jUnTpqUSK09Ul2jirFNWNQ/tuw4tnxI1UALA9PWdmL6+UzqxZs9m+uzZqcTKE9UlmjjrlFXNQ/uuQ8unRA2UiIgESQ2UiIgESQ2UiIgESQ2UiIgESf2gGrpiP6g9JNwPSkSkspr6QekMSkREgqQGSgB4aNXneWjV59OJ9dRTPPTUU6nEyhPVJZo465RVzUP7rkPLp0QNlADw/LsH8/y7B6cTa948np83L5VYeaK6RBNnnbKqeWjfdWj5lKiBEhGRINXZQJnZaDNrk0YyIiIiJVHOoA4G5pnZ42Y2xMyqeexLREQkXnU2UO7+I+AwYAIwEnjDzP7DzNK5oy6paNZoJ80apTOiebOmTWkW6OjJWVJdoomzTlnVPLTvOrR8SiL3gzKzPsC3gCHA/wIDgT+4+w+TS+/T1A8qZuoHJSIZ2+f5oMzsX4B/BtYD44Ex7v6xmTUC3gBSbaBERKRhiDJhYTvgXHf/e+WF7r7LzM5OJi1J27h/HA7AqC6vJx9r0qRCrOHDE4+VJ6pLNHHWKauah/Zdh5ZPSZSHJLpWbZzMbCKAuy9JJCtJ3cvvtefl99qnE2vRIl5etCiVWHmiukQTZ52yqnlo33Vo+ZREaaB6Vn5hZmVA/2TSERERKaixgTKza8zsA6C3mW0q/nwAvANMSS1DERFpkGpsoNz9/7p7K+BWd29d/Gnl7m3d/ZoUcxQRkQaoxockzKy7uy8FnjCzflXfd/c/JZqZpOqgA7anF6tVq9Ri5YnqEk2cdcqq5qF916HlU1JjPygzG+fuo8zsf6t529391GRTq576QcVM/aBEJGN73Q/K3UcV/3tKkomJiIhUJ8pgseebWavi7z8ys8lm9oXkU5M03fVmD+56s0c6sSZO5K6JE1OJlSeqSzRx1imrmof2XYeWT0mUjro/dvcnzOwE4AzgNmAscGyimUmqFn7wmfRiLV2aWqw8UV2iibNOWdU8tO86tHxKovSDKo0g+mXgV+4+BWiSXEoiIiLRGqjVZvZr4AJgmpk1jfg5itNzvGZmy8zs6hrWGWxmC8zsVTObFT11ERGpz6I0NBcAzwBD3P094DPAmLo+VBxx4l7gTOBIYISZHVllnYOA+4Bz3L0ncP7eJC8iIvVXnfeg3P0jM5sCdDCzLsXFUS5YDgCWuftyADP7LTAMWFxpna8Bk939H8VY7+xN8hKfDk22pBerXbvUYuWJ6hJNnHXKquahfdeh5VNS53xQZvY94CfAWmBXcbG7e+86PncehbOubxdffwM41t1HV1rnTqAxhfH+WgG/dPeHq9nWKGAUQJcuXfr//e9/r7qK7Cv1gxKRjO3zfFDA94Ej3H3D3sasZlnVv3wHUBh49jSgOfCSmc1x9z3mfHD3ccA4KHTU3cs8REQkh6I0UCuB9/dh26uAQyq97gysqWad9e7+IfChmT0P9AGSn5RI9vDzvx0FwA8//9fkY40fX4j17W8nHitPVJdo4qxTVjUP7bsOLZ+SKA3UcmCmmf03sK200N1/Ucfn5gGHmVlXYDVwIYV7TpVNAe4xswMoPLp+LHBHxNwlRq99eGB6sd58M7VYeaK6RBNnnbKqeWjfdWj5lERpoP5R/GnCXvR/cvcdZjaawhOAZcAD7v6qmV1afH+suy8xs/8BFlG4vzXe3ZP/J7yIiAQvylN81wOYWcvipbjI3H0aMK3KsrFVXt8K3Lo32xURkfovylh8x5nZYmBJ8XUfM7sv8cxERKRBi3KJ704KY/BNBXD3hWZ2UpJJSfo+13xzerE6dkwtVp6oLtHEWaesah7adx1aPiVR+kG97O7Hmtmf3f0LxWUL3b1PKhlWofmgYqZ+UCKSsf3pB7XSzI4H3MyaAP9C8XKfiIhIUqKMxXcpcDnQiUK/pb7AZQnmJBm44Y0+3PBGOifFN9x7Lzfce28qsfJEdYkmzjplVfPQvuvQ8imJcgZ1hLtfVHmBmQ0CZieTkmTh71vK04u1pmp/bQHVJao465RVzUP7rkPLpyTKGdTdEZeJiIjEpsYzKDM7DjgeaG9mV1R6qzWFjrciIiKJqe0SXxOgvLhOq0rLNwHnJZmUiIhIjQ2Uu88CZpnZg+6u+S3quSNa7st4wPsYq2vX1GLlieoSTZx1yqrmoX3XoeVTEqUf1OHAD4AKKjVo7n5qopnVQP2gYqZ+UCKSsf3pB/UEMBYYD+yMOzEREZHqRGmgdrj7rxLPRDL170v7AfAf3f+UfKw7CjOq/Me//VvisfJEdYkmzjplVfPQvuvQ8imJ0kD9l5ldBjzJnvNBvZtYVpK6tdubpxdr/frUYuWJ6hJNnHXKquahfdeh5VMSpYH65+J/x1Ra5kC3+NMREREpiDIfVJiPd4iISL1WW0fdU919hpmdW9377j45ubRERKShq+0M6mRgBjC0mvccUANVj/Rpld4txT7du6cWK09Ul2jirFNWNQ/tuw4tn5I6+0GFRv2gYqZ+UCKSsZr6QUUZLFZERCR1aqAEgCsXH8OVi49JJ9bNN3PlzTenEitPVJdo4qxTVjUP7bsOLZ+SKI+ZSwPw3o4m6cX64IPUYuWJ6hJNnHXKquahfdeh5VNS5xmUmbUwsx+b2f3F14eZ2dnJpyYiIg1ZlEt8v6EwgsRxxdergBsTy0hERIRoDdTn3f3nwMcA7r4FqObRLxERkfhEuQe13cyaU+j7hJl9nkpj8kn9cOxB69KL1bt3arHyRHWJJs46pV7zF39aiNtmwyevj/9pujlUI9R9L8p8UF8CrgWOBKYDg4CR7j4z8eyqoX5QMUu7H1TxAN1DAAeoSCrU77Ba+zwflLtPN7P5wEAKl/a+7+5hDn0r4Xvp+k8vUwMlItWos4Eys6nAY8BUd/8w+ZQkC5f/dSAA9x41p17FypPLb7gBgHuvuy7jTMIWZ52yqnlox0Co+16UhyRuB04EFpvZE2Z2npk1SzgvSdnWXWVs3VVW72LlydZt29i6Tbd36xJnnbKqeWjHQKj7XpRLfLOAWWZWBpwKfAd4AGidcG4iItKARRpJovgU31BgONAPeCjJpERERKLcg5oEHAv8D3AvMNPddyWdmIiINGxRzqB+A3zN3XcmnYxk56TPvF0vY+XJScekM1hv3sVZp6xqHtoxEOq+V2M/qFBn1FU/qJil3S9D/UCkIdP+X6196QelGXVFRCQzNTZQ7v6T4q83uPubld8zs65RNm5mQ4BfAmXAeHevdsIRMzsGmAMMd/ffRdm2xOuSRYMAmNB7dr2KlSeXXHstABNuuinjTMIWZ52yqnlox0Co+16UflD/Wc2yOhuR4mPp9wJnUhgmaYSZHVnDercAz0TIRUREGogaz6DMrDvQEziwyn2o1kCUjroDgGXuvry4vd8Cw4DFVdb7HoVGMMy7dCIi9VVpbMxNKz55HdDQY7XdgzoCOBs4iD3vQ31AobNuXToBKyu9XkXhcfXdzKwT8FUKHYBrbKDMbBQwCqBLly4RQouISJ1KY2NuGlR8/Wg+Gih3nwJMMbPj3P2lfdh2dXNGVX1c5U7gKnffaVbzFFPuPg4YB4Wn+PYhFxERyZko/aD+bGaXU7jct/vSnrtfXMfnVgGHVHrdGVhTZZ2jgd8WG6d2wFlmtsPdn4qQl8ToS+1W18tYefKlQYOyTiEX4qxTVjUP7RgILZ+SKA3URGApcAZwA3ARsCTC5+YBhxWf+FsNXAh8rfIK7r77aUAzexD4vRqnbAzvuKJexsqT4WedlXUKuRBnnbKqeWjHQGj5lER5iu9Qd/8x8KG7PwR8GehV14fcfQcwmsLTeUuAx939VTO71Mwu3Z+kJX5bdpaxZWc6oyunGStPtmzbxpYAR5QOTZx1yqrmoR0DoeVTEuUM6uPif98zs6OAt4GKKBt392nAtCrLxtaw7sgo25RkjH61MD9NGv0y0oyVJ6OLc/KE1hclNHHWKauah3YMhJZPSZQGapyZtQF+DEwFyoGwZrUSEZF6J8p8UOOLv84CuiWbjoiISEFtHXWvqO2D7v6L+NMREREpqO0MqlVqWYiIiFRRW0fd69NMRLJ1Tod/1MtYeXLOqadmnUIuxFmnrGoe2jEQWj4lUWbUPRz4FdDB3Y8ys97AOe5+Y+LZSWqGdVhZ90o5jJUnw047LesUciHOOmVV89COgdDyKYnSD+p+4BqKj5u7+yIKnW6lHtn4cRM2ftyk3sXKk42bNrFx06as0whenHXKquahHQOh5VMS5THzFu4+t8pYeTsSykcy8oMlhbF60+gHkWasPPnBLbcA6gdVlzjrlFXNQzsGQsunJMoZ1Hoz+zzFgV7N7DzgrUSzEhGRBi/KGdTlFEYS725mq4E3KYzHJyIikpgoHXWXA180s5YUzri2AMOBvyecm4iINGA1XuIzs9Zmdo2Z3WNmpwMfAf8MLAMuSCtBERFpmGo7g5oIbAReojCD7g+BJsBX3H1B8qlJmi747Jv1MlaeXHDmmVmnkAtx1imrmod2DISWT0ltDVQ3d+8FYGbjgfVAF3f/IJXMJFVntK86l2T9iJUnZ5xwQtYp5EKcdcqq5qEdA6HlU1LbU3ylaTZw953Am2qc6q+3tzXj7W3N6l4xZ7Hy5O1163h73bqs0whenHXKquahHQOh5VNS2xlUHzMr9WAzoHnxtQHu7q0Tz05Sc+1r/YF0+kGkGStPrr3zTkD9oOoSZ52yqnlox0Bo+ZTUNhZfeNMriohIgxGlo66IiEjq1ECJiEiQ1ECJiEiQogx1JA3ANzstq5ex8uSbw4ZlnUIuxFmnrGoe2jEQWj4laqAEgJPbrq2XsfLk5AEDsk4hF+KsU1Y1D+0YCC2fEjVQAsCKj8oBqGixuV7FypMVq1cDUNGpU8aZBOrFnwKwYsNWACraNoPjf7pfm8yq5qEdA6HlU6IGSgD42bI+QDr9INKMlSc/u+8+QP2gavTS9QD8bNEgoLj/7GcDlVXNQzsGQsunRA9JiIhIkNRAiYhIkNRAiYhIkNRAiYhIkPSQhADwnUNer5ex8uQ7F2ge0Cji3H+yqnlox0Bo+ZSogRIABrZJb8qBNGPlycA+fbJOIRfi3H+yqnlox0Bo+ZToEp8AsHRza5ZuTmcGlTRj5cnS5ctZunx51mkEL879J6uah3YMhJZPiRooAeDW5b24dXmvehcrT26dMIFbJ0zIOo3gxbn/ZFXz0I6B0PIpUQMlIiJBUgMlIiJBUgMlIiJBUgMlIiJBSvQxczMbAvwSKAPGu/vNVd6/CLiq+HIz8F13X5hkTlK971Usrpex8uR73/hG1inkQpz7T1Y1D+0YCC2fksQaKDMrA+4FTgdWAfPMbKq7V67Em8DJ7r7RzM4ExgHHJpWT1Kxv6431Mlae9O3ePesUciHO/Sermod2DISWT0mSl/gGAMvcfbm7bwd+C+wxfaW7v+jupcrMATonmI/UYsGmNizY1KbexcqTBUuXsmDp0qzTCF6c+09WNQ/tGAgtn5IkG6hOwMpKr1cVl9XkEuDp6t4ws1Fm9oqZvbJuXZg9nvPu7hVHcveKI+tdrDy5e+JE7p44Mes0ghfn/pNVzUM7BkLLpyTJBsqqWebVrmh2CoUG6qrq3nf3ce5+tLsf3b59+xhTFBGRUCX5kMQq4JBKrzsDa6quZGa9gfHAme6+IcF8REQkR5I8g5oHHGZmXc2sCXAhMLXyCmbWBZgMfMPdwxxOV0REMpHYGZS77zCz0cAzFB4zf8DdXzWzS4vvjwWuA9oC95kZwA53PzqpnEREJD8S7Qfl7tOAaVWWja30+7eBbyeZg0Qzpttf6mWsPBlzySVZp5ALce4/WdU8tGMgtHxKNB+UANC9fFO9jJUn3bt1yzqFXIhz/8mq5qEdA6HlU6KhjgSAORvbM2djOk9IphkrT+YsXMichRpIpS5x7j9Z1Ty0YyC0fEp0BiUA3L/ycCCdmTXTjJUn9z/+OKCZdesS5/6TVc1DOwZCy6dEZ1AiIhIkNVAiIhIkNVAiIhIkNVAiIhIkPSQhAPz40PSeZEozVp78+LLLsk4hF+Lcf7KqeWjHQGj5lKiBEgAqWmyul7HypKJTbYP9S0mc+09WNQ/tGAgtnxJd4hMAZm3owKwNHepdrDyZNXcus+bOzTqN4MW5/2RV89COgdDyKdEZlADw8OpDATi57dp6FStPHp4yBYCTBwzIOJOwxbn/ZFXz0I6B0PIp0RmUiIgESQ2UiIgESQ2UiIgESQ2UiIgESQ9JCAA3HTG/XsbKk5v+9V+zTiEX4tx/sqp5aMdAaPmUqIESAA5uurVexsqTg9uHN91BiOLcf7KqeWjHQGj5lOgSnwDwzLqOPLOuY72LlSfPvPACz7zwQtZpBC/O/Sermod2DISWT4nOoASAx9/qCsAZ7dfUq1h58vjTTwNwxgknZJxJ2OLcf7KqeWjHQGj5lOgMSkREgqQGSkREgqQGSkREgqQGSkREgqSHJASA23rMq5ex8uS2q67KOoVciHP/yarmoR0DoeVTogZKAGjTeHu9jJUnbVq3zjqFXIhz/8mq5qEdA6HlU6JLfALAlLWHMGXtIfUuVp5Mee45pjz3XNZpBC/O/Sermod2DISWT4kaKAFg6touTF3bpd7FypOpM2YwdcaMrNMIXpz7T1Y1D+0YCC2fEjVQIiISJDVQIiISJDVQIiISJDVQIiISJD1mLgDc03NOvYyVJ/dcd13WKeRCnPtPVjUP7RgILZ8SNVACQPOynfUyVp40b9o06xRyIc79J6uah3YMhJZPiS7xCQCT1lQwaU1FvYuVJ5OmTWPStGlZpxG8OPefrGoe2jEQWj4laqAEgOnrOzF9fad6FytPps+ezfTZs7NOI3hx7j9Z1Ty0YyC0fErUQImISJASbaDMbIiZvWZmy8zs6mreNzO7q/j+IjPrl2Q+IiKSH4k1UGZWBtwLnAkcCYwwsyOrrHYmcFjxZxTwq6TyERGRfEnyDGoAsMzdl7v7duC3wLAq6wwDHvaCOcBBZvbZBHMSEZGcMHdPZsNm5wFD3P3bxdffAI5199GV1vk9cLO7v1B8/Rxwlbu/UmVboyicYQEcBfw1kaTj1w5Yn3USESjP+OUl17zkCfnJNS95Qji5fs7d21ddmGQ/KKtmWdXWMMo6uPs4YByAmb3i7kfvf3rJy0uuyjN+eck1L3lCfnLNS54Qfq5JXuJbBVSeYKQzsGYf1hERkQYoyQZqHnCYmXU1sybAhcDUKutMBb5ZfJpvIPC+u7+VYE4iIpITiV3ic/cdZjYaeAYoAx5w91fN7NLi+2OBacBZwDLgI+BbETY9LqGUk5CXXJVn/PKSa17yhPzkmpc8IfBcE3tIQkREZH9oJAkREQmSGigREQlSMA2UmT1gZu+YWbV9nGobFqmuIZUyyPWiYo6LzOxFM+tT6b0VZvYXM1tgZq9U9/kU8xxsZu8Xc1lgZtdVei+1mkbIc0ylHP9qZjvN7DPF99Ks5yFm9r9mtsTMXjWz71ezThD7acRcM99PI+YZyn4aJdfM91Uza2Zmc81sYTHP66tZJ4j9tE7uHsQPcBLQD/hrDe+fBTxNoe/UQODl4vIy4G9AN6AJsBA4MuNcjwfaFH8/s5Rr8fUKoF0gNR0M/L6a5anWtK48q6w7FJiRUT0/C/Qr/t4KeL1qXULZTyPmmvl+GjHPUPbTOnMNYV8t7nvlxd8bAy8DA0PcT+v6CeYMyt2fB96tZZWahkWKMqRSqrm6+4vuvrH4cg6F/l2pi1DTmqRa073McwTwWFK51Mbd33L3PxV//wBYAlSdoyCI/TRKriHspxFrWpPgalpFJvtqcd/bXHzZuPhT9Wm4IPbTugTTQEXQCVhZ6fWq4rKalofiEgr/UilxYLqZzbfCEE5ZO654KeBpM+tZXBZkTc2sBTAE+M9KizOpp5lVAF+g8K/TyoLbT2vJtbLM99M68gxqP62rplnvq2ZWZmYLgHeAP7h78PtpdfI05XtNwyJFGi4pC2Z2CoUD/4RKiwe5+xoz+yfgD2a2tHgGkYU/URgDa7OZnQU8RWFk+VBrOhSY7e6Vz7ZSr6eZlVP4w/Ov7r6p6tvVfCSz/bSOXEvrZL6f1pFnUPtplJqS8b7q7juBvmZ2EPCkmR3l7pXv8Qa1n9YkT2dQNQ2LFORwSWbWGxgPDHP3DaXl7r6m+N93gCcpnFJnwt03lS4FuPs0oLGZtSPQmlIYjWSPSyZp19PMGlP44/SIu0+uZpVg9tMIuQaxn9aVZ0j7aZSaFmW+rxZjvQfMpHA2V1kw+2mtsrr5Vd0PUEHNN/S/zJ439eYWlx8ALAe68slNvZ4Z59qFwugYx1dZ3hJoVen3FymM+J5VngfzSWftAcA/ivVNvaa15Vl8/0AK96laZlXPYm0eBu6sZZ0g9tOIuWa+n0bMM4j9NEquIeyrQHvgoOLvzYE/AmeHuJ/W9RPMJT4ze4zC0zrtzGwV8BMKN/fwWoZF8hqGVMo41+uAtsB9ZgawwwsjBnegcLoNhR3hUXf/nwzzPA/4rpntALYAF3phL021phHyBPgqMN3dP6z00VTrCQwCvgH8pXh9H+DfKfyhD20/jZJrCPtplDyD2E8j5grZ76ufBR6ywqSxjYDH3f33FmGYuSz+ntZGQx2JiEiQ8nQPSkREGhA1UCIiEiQ1UCIiEiQ1UCIiEiQ1UCIiEiQ1UCIiEiQ1UCIiEiQ1UCI5YGbfLs4l9K0I63Yzswlm9rs0chNJihookXz4P8CpwPl1reiFqRIuST4lkWSpgRKJgZl1NrMpZvaGmS03s3vMrOk+bqu5mc0qDlVT8jKFqRNerrReLzP7fZWff6phm03M7HkzC2Z4M5G6qIES2U9WGGBtMvCUux9GYSqI5sDP93GTFwOTvTBlQkk5hUE/DywtcPe/uPvZVX7eqW6DXph87jlg+D7mJJI6NVAi++9UYKu7/wZ2z8Xzb8A3i3MH7a2LgCmlF2bWiMIApN8EvlrlzOpTzKytmY0FvmBm11R666nitkVyQaf7IvuvJzC/8gJ332RmK4BDgQVRN2RmTYBu7r6i0uJTgUXuvsLMFhZf/6GmbXhhXqdLq3nrr8AxUXMRyZrOoET2n1H9rKPVzU5al3bAe1WWXcQnk989xj6eBRXP7LabWat9+bxI2tRAiey/V4GjKy8ws9YU5gDqYWb/Ymanmdm/F9+7w8w6mVlTM5tYZVtbgGaVttMcGAbcVjwjuxUYVly+L5oCW/fxsyKpUgMlsv+eA1qY2TcBiveIbgfucffHgMMpXHK7pbj+vcAoCtOC/7byhtx9I1BmZqVG6hzgaXfv4u4V7t4F+C9g6N4maWZtgXXu/vHeflYkC2qgRPZTcXbXrwLnmdkbwAZgl7vfVHzC7yBgF8V7vu6+DOjIJ9NuVzUdOKH4+0XAk1XefxL4+j6kegqFmVRFckEz6orEzMyOp3Cv6FzgROAN4G1guLv/sLjOhcBn3f2Oaj7/BeAKd/9GzHlNBq5x99fi3K5IUtRAiaTMzAYB3wcucfcPaljnYuChKn2h9idmE+BCd384ju2JpEENlIiIBEn3oEREJEhqoEREJEhqoEREJEhqoEREJEhqoEREJEhqoEREJEhqoEREJEhqoEREJEj/H/gvr0tYmEpOAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "a = bestmin2[0][0]\n",
    "b = bestmin2[0][1]\n",
    "gamma = bestmin2[0][2]\n",
    "fractional_coords = [[0,0,0]]\n",
    "atoms_dummy = Atoms(['C'], scaled_positions=fractional_coords, cell=[a, b, 100, 90, 90, gamma])\n",
    "\n",
    "result = get_reflexes(atoms_dummy, dummy=False)\n",
    "\n",
    "fig, ax = plot_result_gixd(result, reflections_exp)\n",
    "fig.tight_layout()\n",
    "fig.savefig('SOS/single_site.png', dpi=300)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "ab4925fc",
   "metadata": {},
   "source": [
    "Best result for a=b=4.52 Å and γ=112.6° for unit cell with one basis particle. However, reflexes are missing compared to experimental result."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "4bbe9a42",
   "metadata": {},
   "source": [
    "This cell is equivalent to a two atom basis with a rectangular unit cell with a=7.529 Å, b=5.007 Å, γ=90°:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "a70f1e15",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(<Figure size 432x288 with 1 Axes>,\n",
       " <AxesSubplot:xlabel='Q$_\\\\mathrm{xy}$ ($\\\\mathrm{\\\\AA}^{-1}$)', ylabel='Relative intensity'>)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAArz0lEQVR4nO3deZgV9ZX/8fehZW9QAgwGkDTEBUSWgCKKC2qMaEQSfyoSkwzRhMcomcxoiDomJhqdn0aNxi0Ewaj8XNAMCpPBEZUBI4ogCZIIqARJWFwAUUQBBc7vj3tv27S9FFDLt7o/r+fpx7731q1zPLeqD3Wrvt8yd0dERCQ0TbJOQEREpCZqUCIiEiQ1KBERCZIalIiIBEkNSkREgqQGJSIiQUqsQZnZPWb2jpn9tZbXzcxuM7PlZrbYzAYklYuIiORPkkdQ9wLD6nj9VOCg4s8Y4DcJ5iIiIjmTWINy92eBd+tYZARwvxfMA/Yzs88nlY+IiOTLPhnG7gKsqvJ4dfG5N6svaGZjKBxl0bp164E9e/ZMJcFGYfNaAFa+uw2Ais81h/LOiYZcuWZNIVaXLonGyRvVJZo465RVzUP7rLPOZ+HChevdvWP157NsUFbDczXOu+TuE4AJAIcffri/9NJLSebVuNxc+BguWDwEgEl958KlaxMNecGVVxZiXXddonHyRnWJJs46ZVXz0D7rrPMxs7/X9HyWDWo1cECVx12BZP8ySq2+3WV5erFGjEgtVp6oLtHEWaesah7aZx1aPiWW5GSxZlYB/MHdD6vhta8CY4HTgCOB29x9UH3r1BFUzG6u4UD2Uk0gLCLpMbOF7n549ecTO4Iys4eAoUAHM1sN/AxoCuDu44EZFJrTcuAj4DtJ5SL1W/lROQAVrTYnHyuw799DobpEo3NQ8Qstn5LEGpS7j6rndQcuTiq+7J5fLO8HFM9BJR3rrrsKsQL5/j0Uqks0cdXpk08+YdHixbTfd1+2bNoUR2qRbXjvPYDU49YmrXxatGhB165dadq0aaTlszwHJSKSmdWrV3Nwjx60bN2a7l27pho7tCOWNPJxdzZs2MDq1avp3r17pPdoqiMRaZS2bt1Ky9atMavpgmKJm5nRvn17tm7dGvk9alAi0mipOaVrd+utBiUiIkHSOSgB4HsHvJZerHPOSS1Wnqgu0cRZpw7t2sW2rj2Ne9111/Hggw9SVlZGkyZN+O1vf8vdd9/NJZdcwqGHHrrXsSoqKnjppZfo0KFDrcv8buJExo0bV/n46KOP5vnnn9/r2Hsr0XFQSdA4qJhpHJQ0UkuXLqVXr16Z5vDCCy9wySWXMHv2bJo3b8769ev5+OOP6dw5vunGojSo8vJyNm9OfogJ1Fz32sZB6Ss+AWDZ5rYs29w2nVgrVrBsxYpUYuWJ6hJNnHXasm0bW7YV5qHk+Z8X/sEW18/zP6837ptvvkmHDh1o3rw5AB06dKBz584MHTqU0j/Ey8vLueyyyxg4cCBf/vKXmT9/PkOHDqVHjx5Mnz4dgHvvvZexY8dWrv/0009n9uzZn4n7ta99jYEDB9K7d28mTJgAwOWXX86WLVvo268f5513XmVMKFx5N27cOA477DD69OnDlClTAJg9ezZDhw7lrLPOomfPnpx33nkkcbCjBiUA3LiiDzeu6JNOrEmTuHHSpFRi5YnqEk2cdXp7/XreXr8+lnXtSdyvfOUrrFq1ioMPPpiLLrqIOXPmfGbZDz/8kKFDh7Jw4ULatGnDT37yE5566ikee+wxrrrqqt2Ke88997Bw4UJeeuklbrvtNjZs2MD1119PixYtmD5jBg888MAuy0+dOpVFixbx8ssv8/TTTzNu3DjefLMwn/ef//xnbr31VpYsWcKKFSuYOzf+MZRqUCIiGSkvL2fhwoVMmDCBjh07MnLkSO69995dlmnWrBnDhhVurdenTx+OP/54mjZtSp8+fVi5cuVuxbvtttvo168fgwcPZtWqVbz++ut1Lv/cc88xatQoysrK6NSpE8cffzwLFiwAYNCgQXTt2pUmTZrQv3//3c4lCl0kISKSobKyMoYOHcrQoUPp06cP99133y6vN23atPLy7CZNmlR+HdikSRO2b98OwD777MPOnTsr31PTWKPZs2fz9NNP88ILL9CqVSuGDh1a75ikur62K+VR+n8o5RInHUGJiGTk1Vdf3eUoZtGiRXzhC1/Y7fVUVFSwaNEidu7cyapVq5g/f/5nlnn//fdp164drVq1YtmyZcybN6/ytX2aNuWTTz75zHuOO+44pkyZwo4dO1i3bh3PPvssgwbVO6d3bHQEJSICcPTPCz8p2rx5Mz/4wQ9477332GeffTjwwAOZMGECZ5111m6tZ8iQIXTv3p0+ffpw2GGHMWDAgM8sM2zYMMaPH0/fvn055JBDGDx4cOVro77xDU49+WSOHDRol/NQX//613nhhRfo168fZsYvf/lL9t9/f5YtW7bn/9O7QZeZN3bFy8wXbSqMy+jfdmPil5kvKm7c/XVn5F2oLtHEVaelS5fyhYoKAFq1bLm3ae2Wj7ZsySRubdLMZ3cuM9cRlADFxpRWLP0BrpHqEk2cdcqqQYTSmEpCy6dE56AEKBxBlY6iEo+1bFnlv4LlU6pLNHHW6aMtWyqPHtKUVdzahJZPiRqUAHD7ykO5feXeT6sSKdbkydw+eXIqsfJEdYkmzjq98+67vPPuu7GsKw9xaxNaPiVqUCIiEiQ1KBERCZIalIiIBEkNSkQkI2VlZfTv37/y5/rrr0803vTp0xOPMXv27Nhu1aHLzAWAcT3+kl6sCy5ILVaeqC7RxFmnTnXcgiJJpbgtW7Zk0aJFqcTcvn07Z5xxBmeccUat+cRh9uzZlJeXc/TRR+/1unQEJQD0LN9Ez/JN6cTq0YOePXqkEitPVJdo4qxTy+bNaVllTrm01BX3/fff55BDDuHVV18FYNSoUdx9991AYXLZSy+9lAEDBnDSSSexbt06AP72t78xbNgwBg4cyLHHHls508Po0aO55JJLOOGEE7jssst2uS3H6NGj+f73v88JJ5xA7169mD9vHueffz69evVi9OjRlfnMnDmTo446igEDBnD22WdX3jeqoqKCn/3sZwwYMIA+ffqwbNkyVq5cyfjx47nlllvo378/f/zjH/eqTjqCEgDmbewIwOB265KP9fLLhVj9+iUeK09Ul2jirNPmjz4CoLxVKy648srPvP6VIUMYedppbNm2jbHXXPOZ18848URGnHQSGzdt4kc33LDLa5Ouu67euFu2bKF///6Vz19xxRWMHDmSO+64g9GjR/PDH/6QjRs38r3vfQ8o3HpjwIAB3HzzzVxzzTVcffXV3HHHHYwZM4bx48dz0EEH8eKLL3LRRRcxa9YsAF577TWefvppysrKPjNT+saNG5k1axZTHn2U4cOHM3fuXCZOnMgRRxzBokWL6Nq1K9deey1PP/00rVu35oYbbuBXv/pV5W0+OnTowJ/+9CfuuusubrrpJiZOnMiFF15IeXk5P/rRj+qpfv3UoASAu1cdDKTToO5+5JFCLP0h3oXqEk2cdVq/sTCDSnmrVnu9rj2JW9tXfCeffDKPPvooF198MS8XGzIUZjAfOXIkAN/85jc588wz2bx5M88//zxnn3125XLbSjdhBM4++2zKyspqzGP48OGYGZ27dqV9+/b06VO4J1zv3r1ZuXIlq1evZsmSJQwZMgSAjz/+mKOOOqry/WeeeSYAAwcOZOrUqXtSijqpQYmIUPcRT8vmzet8vV3btnW+vrt27tzJ0qVLadmyJe+++y5du3atcTkzY+fOney33361nstq3bp1rXGq3rqjWZWvHEu38igrK+Pkk0/moYceqvP9ut2GiEgjccstt9CrVy8eeughzj///MpbYezcuZPf//73ADz44IMcc8wxtG3blu7du/Poo48ChXs4VT3q2huDBw9m7ty5LF++HICPPvqI1157rc73tGnThg8++CCW+GpQIiIZKZ2DKv1cfvnlvPbaa0ycOJGbb76ZY489luOOO45rr70WKBwNvfLKKwwcOJBZs2ZVngt64IEHmDRpEv369aN3795MmzYtlvw6duzIvffey6hRo+jbty+DBw+u91Ybw4cP57HHHovlIgndbqOxK95u44LFhe+YJ/Wdm/jtNkono+P8SqQhUF2iiatOS5cupWXbtgBUdOmy13ntjpVr1uxR3PLy8sqr6ELIZ0/odhuy2356YDxfCUSKddFFqcXKE9Ulmjjr9PmOHWNbVx7i1ia0fErUoASAilbx/6us1lgp/2s1L1SXaOKsU/NmzWJbVxpxkzh6guzqUB+dgxIA5mzoxJwNndKJNX8+c+bPTyVWnqgu0cRZp02bN/PBhx/Gsq7d8cGHH2YStzZp5bO7p5R0BCUA3L/mQACOb/928rGKJ3CPHzQo8Vh5orpEE1edWrRoweo1a2jZujVt6rgUOwkb3nsPIPW4tUkjH3dnw4YNtGjRIvJ71KBEpFHq2rUr//3007Tfd1+2xnRZdFSlhrBlUzrTi9UnrXxatGhR65iumqhBiUij1LRpU/7rueeA9K+cDO2KzdDyKdE5KBERCZIalIiIBCnRgbpmNgz4NVAGTHT366u9vi/w/4BuFL5uvMndf1fXOjVQN2bFgbpvbSucuNy/+dbEB+q+VbxFwP6Bjr3IiuoSTZx1yqrmoX3WWeeT+kBdMysD7gROBlYDC8xsursvqbLYxcASdx9uZh2BV83sAXf/OKm8pGb7N9+aXqxAdsrQqC7RxFmnrGoe2mcdWj4lSX7FNwhY7u4rig3nYWBEtWUcaGNmBpQD7wLxT4kr9XpyXWeeXNc5nVjPPceTxZPT8inVJZo465RVzUP7rEPLpyTJq/i6AKuqPF4NHFltmTuA6cBaoA0w0t13Vl+RmY0BxgB069YtkWQbu0fe7A7AKR3XJh/riScKsY45JvFYeaK6RBNnnbKqeWifdWj5lCR5BGU1PFf95MYpwCKgM9AfuMPM2n7mTe4T3P1wdz+8Y6CHoiIiEq8kG9Rq4IAqj7tSOFKq6jvAVC9YDrwB9EwwJxERyYkkG9QC4CAz625mzYBzKXydV9U/gJMAzKwTcAiwIsGcREQkJxI7B+Xu281sLPAkhcvM73H3V8zswuLr44FfAPea2V8ofCV4mbuvTyonERHJD92wsLErjoPa+Elhuv12TT9OfBzUxuJ8X+3afuZ0Y6OmukQTZ52yqnlon3XW+eiGhVKndk3TG3oWyk4ZGtUlmjjrlFXNQ/usQ8unRFMdCQDT3j6AaW8fUP+CccR65hmmPfNMKrHyRHWJJs46ZVXz0D7r0PIpUYMSAKa/3Y3pb6czxmz6rFlMnzUrlVh5orpEE2edsqp5aJ91aPmUqEGJiEiQ1KBERCRIalAiIhIkNSgREQmSxkE1dsVxUFt2lAHQsmxH4uOgtmzbVojVvHmicfJGdYkmzjplVfPQPuus89E4KKlTy7Id6cUKZKcMjeoSTZx1yqrmoX3WoeVToq/4BIApayuYsrYinVgzZjBlxoxUYuWJ6hJNnHXKquahfdah5VOiBiUAzFzfhZnru6QTa+5cZs6dm0qsPFFdoomzTlnVPLTPOrR8StSgREQkSGpQIiISJDUoEREJkhqUiIgESeOgGrviOKhdJDwOSkSkqtrGQekISkREgqQGJQDct/qL3Lf6i+nEevxx7nv88VRi5YnqEk2cdcqq5qF91qHlU6IGJQA8++7+PPvu/unEWrCAZxcsSCVWnqgu0cRZp6xqHtpnHVo+JWpQIiISpHoblJmNNbN2aSQjIiJSEuUIan9ggZk9YmbDzKyGy75ERETiVW+DcvefAAcBk4DRwOtm9h9mls4ZdUlFiyY7aNEknRnNWzRvTotAZ0/OkuoSTZx1yqrmoX3WoeVTEnkclJn1A74DDAP+FxgMPOXuP04uvc/SOKiYaRyUiGRsj+8HZWb/AvwzsB6YCIxz90/MrAnwOpBqgxIRkcYhyg0LOwBnuvvfqz7p7jvN7PRk0pK0TfjHwQCM6fZa8rGmTCnEGjky8Vh5orpEE2edsqp5aJ91aPmURLlIonv15mRmkwHcfWkiWUnqXnyvIy++1zGdWIsX8+LixanEyhPVJZo465RVzUP7rEPLpyRKg+pd9YGZlQEDk0lHRESkoNYGZWZXmNkHQF8z21T8+QB4B5iWWoYiItIo1dqg3P3/unsb4EZ3b1v8aePu7d39ihRzFBGRRqjWiyTMrKe7LwMeNbMB1V939z8lmpmkar99Pk4vVps2qcXKE9UlmjjrlFXNQ/usQ8unpNZxUGY2wd3HmNn/1vCyu/uJyaZWM42DipnGQYlIxnZ7HJS7jyn+94QkExMREalJlMlizzazNsXff2JmU83sS8mnJmm67Y1e3PZGr3RiTZ7MbZMnpxIrT1SXaOKsU1Y1D+2zDi2fkigDdX/q7o+a2THAKcBNwHjgyEQzk1S9/MHn0ou1bFlqsfJEdYkmzjplVfPQPuvQ8imJMg6qNIPoV4HfuPs0oFlyKYmIiERrUGvM7LfAOcAMM2se8X0Ub8/xqpktN7PLa1lmqJktMrNXzGxO9NRFRKQhi9JozgGeBIa5+3vA54Bx9b2pOOPEncCpwKHAKDM7tNoy+wF3AWe4e2/g7N1JXkREGq56z0G5+0dmNg3oZGbdik9H+cJyELDc3VcAmNnDwAhgSZVlvgFMdfd/FGO9szvJS3w6NduSXqwOHVKLlSeqSzRx1imrmof2WYeWT0m994Mysx8APwPeBnYWn3Z371vP+86icNT13eLjbwFHuvvYKsvcCjSlMN9fG+DX7n5/DesaA4wB6Nat28C///3v1ReRPaVxUCKSsT2+HxTwQ+AQd9+wuzFreK76X759KEw8exLQEnjBzOa5+y73fHD3CcAEKAzU3c08REQkh6I0qFXA+3uw7tXAAVUedwXW1rDMenf/EPjQzJ4F+gHJ35RIdvHLvx0GwI+/+NfkY02cWIj13e8mHitPVJdo4qxTVjUP7bMOLZ+SKA1qBTDbzP4b2FZ60t1/Vc/7FgAHmVl3YA1wLoVzTlVNA+4ws30oXLp+JHBLxNwlRq9+uG96sd54I7VYeaK6RBNnnbKqeWifdWj5lERpUP8o/jRjN8Y/uft2MxtL4QrAMuAed3/FzC4svj7e3Zea2f8Aiymc35ro7sn/E15ERIIX5Sq+qwHMrHXxq7jI3H0GMKPac+OrPb4RuHF31isiIg1flLn4jjKzJcDS4uN+ZnZX4pmJiEijFuUrvlspzME3HcDdXzaz45JMStL3hZab04vVuXNqsfJEdYkmzjplVfPQPuvQ8imJMg7qRXc/0sz+7O5fKj73srv3SyXDanQ/qJhpHJSIZGxvxkGtMrOjATezZsC/UPy6T0REJClR5uK7ELgY6EJh3FJ/4KIEc5IMXPN6P655PZ2D4mvuvJNr7rwzlVh5orpEE2edsqp5aJ91aPmURDmCOsTdz6v6hJkNAeYmk5Jk4e9bytOLtbb6eG0B1SWqOOuUVc1D+6xDy6ckyhHU7RGfExERiU2tR1BmdhRwNNDRzC6p8lJbCgNvRUREElPXV3zNgPLiMm2qPL8JOCvJpERERGptUO4+B5hjZve6u+5v0cAd0npP5gPew1jdu6cWK09Ul2jirFNWNQ/tsw4tn5Io46AOBn4EVFClobn7iYlmVguNg4qZxkGJSMb2ZhzUo8B4YCKwI+7EREREahKlQW13998knolk6t+XDQDgP3r+KflYtxTuqPIf//ZvicfKE9UlmjjrlFXNQ/usQ8unJEqD+i8zuwh4jF3vB/VuYllJ6t7+uGV6sdavTy1Wnqgu0cRZp6xqHtpnHVo+JVEa1D8X/zuuynMO9Ig/HRERkYIo94MK8/IOERFp0OoaqHuiu88yszNret3dpyaXloiINHZ1HUEdD8wChtfwmgNqUA1IvzbpnVLs17NnarHyRHWJJs46ZVXz0D7r0PIpqXccVGg0DipmGgclIhmrbRxUlMliRUREUqcGJQBcuuQILl1yRDqxrr+eS6+/PpVYeaK6RBNnnbKqeWifdWj5lES5zFwagfe2N0sv1gcfpBYrT1SXaOKsU1Y1D+2zDi2fknqPoMyslZn91MzuLj4+yMxOTz41ERFpzKJ8xfc7CjNIHFV8vBq4NrGMREREiNagvujuvwQ+AXD3LUANl36JiIjEJ8o5qI/NrCWFsU+Y2RepMiefNAxH7rcuvVh9+6YWK09Ul2jirFPqNX/+54W47TZ8+vjon6ebQw1C3fai3A/qK8CVwKHATGAIMNrdZyeeXQ00DipmaY+DKu6guwhgBxVJhcYd1miP7wfl7jPNbCEwmMJXez909zCnvpXwvXD1Z59TgxKRGtTboMxsOvAQMN3dP0w+JcnCxX8dDMCdh81rULHy5OJrrgHgzquuyjiTsMVZp6xqHto+EOq2F+UiiZuBY4ElZvaomZ1lZi0SzktStnVnGVt3ljW4WHmydds2tm7T6d36xFmnrGoe2j4Q6rYX5Su+OcAcMysDTgS+B9wDtE04NxERacQizSRRvIpvODASGADcl2RSIiIiUc5BTQGOBP4HuBOY7e47k05MREQatyhHUL8DvuHuO5JORrJz3OfeapCx8uS4I9KZrDfv4qxTVjUPbR8IddurdRxUqHfU1TiomKU9LkPjQKQx0/Zfoz0ZB6U76oqISGZqbVDu/rPir9e4+xtVXzOz7lFWbmbDgF8DZcBEd6/xhiNmdgQwDxjp7r+Psm6J1wWLhwAwqe/cBhUrTy648koAJl13XcaZhC3OOmVV89D2gVC3vSjjoP6zhufqbSLFy9LvBE6lME3SKDM7tJblbgCejJCLiIg0ErUeQZlZT6A3sG+181BtgSgDdQcBy919RXF9DwMjgCXVlvsBhSYY5lk6EZGGqjQ35qaVnz4OaOqxus5BHQKcDuzHruehPqAwWLc+XYBVVR6vpnC5eiUz6wJ8ncIA4FoblJmNAcYAdOvWLUJoERGpV2luzE1Dio8fzEeDcvdpwDQzO8rdX9iDddd0z6jql6vcClzm7jvMar/FlLtPACZA4Sq+PchFRERyJso4qD+b2cUUvu6r/GrP3c+v532rgQOqPO4KrK22zOHAw8Xm1AE4zcy2u/vjEfKSGH2lw5oGGStPvjJkSNYp5EKcdcqq5qHtA6HlUxKlQU0GlgGnANcA5wFLI7xvAXBQ8Yq/NcC5wDeqLuDulVcDmtm9wB/UnLIxsvPKBhkrT0aedlrWKeRCnHXKquah7QOh5VMS5Sq+A939p8CH7n4f8FWgT31vcvftwFgKV+ctBR5x91fM7EIzu3Bvkpb4bdlRxpYd6cyunGasPNmybRtbApxROjRx1imrmoe2D4SWT0mUI6hPiv99z8wOA94CKqKs3N1nADOqPTe+lmVHR1mnJGPsK4X706QxLiPNWHkytnhPntDGooQmzjplVfPQ9oHQ8imJ0qAmmFk74KfAdKAcCOuuViIi0uBEuR/UxOKvc4AeyaYjIiJSUNdA3UvqeqO7/yr+dERERArqOoJqk1oWIiIi1dQ1UPfqNBORbJ3R6R8NMlaenHHiiVmnkAtx1imrmoe2D4SWT0mUO+oeDPwG6OTuh5lZX+AMd7828ewkNSM6rap/oRzGypMRJ52UdQq5EGedsqp5aPtAaPmURBkHdTdwBcXLzd19MYVBt9KAbPykGRs/adbgYuXJxk2b2LhpU9ZpBC/OOmVV89D2gdDyKYlymXkrd59fba687QnlIxn50dLCXL1pjINIM1ae/OiGGwCNg6pPnHXKquah7QOh5VMS5QhqvZl9keJEr2Z2FvBmolmJiEijF+UI6mIKM4n3NLM1wBsU5uMTERFJTJSBuiuAL5tZawpHXFuAkcDfE85NREQasVq/4jOztmZ2hZndYWYnAx8B/wwsB85JK0EREWmc6jqCmgxsBF6gcAfdHwPNgK+5+6LkU5M0nfP5NxpkrDw559RTs04hF+KsU1Y1D20fCC2fkroaVA937wNgZhOB9UA3d/8glcwkVad0rH4vyYYRK09OOeaYrFPIhTjrlFXNQ9sHQsunpK6r+Eq32cDddwBvqDk1XG9ta8Fb21rUv2DOYuXJW+vW8da6dVmnEbw465RVzUPbB0LLp6SuI6h+ZlYawWZAy+JjA9zd2yaenaTmylcHAumMg0gzVp5ceeutgMZB1SfOOmVV89D2gdDyKalrLr7wbq8oIiKNRpSBuiIiIqlTgxIRkSCpQYmISJCiTHUkjcC3uyxvkLHy5NsjRmSdQi7EWaesah7aPhBaPiVqUALA8e3fbpCx8uT4QYOyTiEX4qxTVjUPbR8ILZ8SNSgBYOVH5QBUtNrcoGLlyco1awCo6NIl40wC9PzPK39duWErABXDr9/r1WZV89D2gdDyKVGDEgB+sbwfkM44iDRj5ckv7roL0DioGr1wdeWvv1g8BIBJw/d+tVnVPLR9ILR8SnSRhIiIBEkNSkREgqQGJSIiQVKDEhGRIOkiCQHgewe81iBj5cn3ztF9QKOIc/vJquah7QOh5VOiBiUADG6X3i0H0oyVJ4P79cs6hVyIc/vJquah7QOh5VOiBiUALNtcuHtKz/JN9SyZr1h5smzFCgB69uiRcSZhq9x+4lhXRjUPbR8ILZ8SNSgB4MYVfYB0xkGkGStPbpw0CdA4qPpUbj9xrCujmoe2D4SWT4kukhARkSCpQYmISJDUoEREJEhqUCIiEqREL5Iws2HAr4EyYKK7X1/t9fOAy4oPNwPfd/eXk8xJavaDiiUNMlae/OBb38o6hVyIc/vJquah7QOh5VOSWIMyszLgTuBkYDWwwMymu3vVSrwBHO/uG83sVGACcGRSOUnt+rfd2CBj5Un/nnFcON3wxbn9ZFXz0PaB0PIpSfIIahCw3N1XAJjZw8AIoLJBufvzVZafB3RNMB+pw6JN7YB0NtQ0Y+XJomXLADWq+lRuP3GsK6Oah7YPhJZPSZLnoLoAq6o8Xl18rjYXAE/U9IKZjTGzl8zspXXrwhzxnHe3rzyU21ce2uBi5cntkydz++TJWacRvDi3n6xqHto+EFo+JUk2KKvhOa9xQbMTKDSoy2p63d0nuPvh7n54x44dY0xRRERCleRXfKuBA6o87gqsrb6QmfUFJgKnuvuGBPMREZEcSfIIagFwkJl1N7NmwLnA9KoLmFk3YCrwLXcPczpdERHJRGJHUO6+3czGAk9SuMz8Hnd/xcwuLL4+HrgKaA/cZWYA29398KRyEhGR/Eh0HJS7zwBmVHtufJXfvwt8N8kcJJpxPf7SIGPlybgLLsg6hVyIc/vJquah7QOh5VOi2cwFSHea/dCm9A+FbrMRTZzbT1Y1D20fCC2fEjUoAWDexsLVkWncuCzNWHky7+XCJCq6cWHdKrefONaVUc1D2wdCy6dEDUoAuHvVwUA6G2iasfLk7kceAdSg6lO5/cSxroxqHto+EFo+JZosVkREgqQGJSIiQVKDEhGRIKlBiYhIkHSRhADw0wPTuw1XmrHy5KcXXZR1CrkQ5/aTVc1D2wdCy6dEDUoAqGi1uUHGypOKLnVN9i8lcW4/WdU8tH0gtHxK1KAEgDkbOgFwfPu3G1SsPJkzfz4Axw8alHEmYavcfuJYV0Y1D20fCC2fEjUoAeD+NQcC6WygacbKk/unTQPUoOpTuf3Esa6Mah7aPhBaPiW6SEJERIKkBiUiIkFSgxIRkSCpQYmISJB0kYQAcN0hCxtkrDy57l//NesUciHO7Sermoe2D4SWT4kalACwf/OtDTJWnuzfsWPWKeRCnNtPVjUPbR8ILZ8SNSgB4Ml1nQE4pePaBhUrT5587jkATjnmmIwzCVvl9hPHujKqeWj7QGj5lKhBCQCPvNkdSGcDTTNWnjzyxBOAGlR9KrefONaVUc1D2wdCy6dEF0mIiEiQ1KBERCRIalAiIhIkNSgREQmSLpIQAG7qtaBBxsqTmy67LOsUciHO7Sermoe2D4SWT4kalADQrunHDTJWnrRr2zbrFHIhzu0nq5qHtg+Elk+JGpQAMO3tAwAY0WlVg4qVJ9OeeQaAESedlHEmYavcfuJYV0Y1D20fCC2fEjUoAWD6292AdDbQNGPlyfRZswA1qPpUbj9xrCujmoe2D4SWT4kukhARkSCpQYmISJDUoEREJEhqUCIiEiRdJCEA3NF7XoOMlSd3XHVV1inkQpzbT1Y1D20fCC2fEjUoAaBl2Y4GGStPWjZvnnUKuRDn9pNVzUPbB0LLp0QNSgCYsrYCgJGdVzaoWHkyZcYMAEaedlrGmYStcvuJY10Z1Ty0fSC0fEp0DkoAmLm+CzPXd2lwsfJk5ty5zJw7N+s0ghfn9pNVzUPbB0LLp0QNSkREgpRogzKzYWb2qpktN7PLa3jdzOy24uuLzWxAkvmIiEh+JNagzKwMuBM4FTgUGGVmh1Zb7FTgoOLPGOA3SeUjIiL5kuQR1CBgubuvcPePgYf57PRZI4D7vWAesJ+ZfT7BnEREJCfM3ZNZsdlZwDB3/27x8beAI919bJVl/gBc7+7PFR8/A1zm7i9VW9cYCkdYAIcBf00k6fh1ANZnnUQEyjN+eck1L3lCfnLNS54QTq5fcPeO1Z9M8jJzq+G56t0wyjK4+wRgAoCZveTuh+99esnLS67KM355yTUveUJ+cs1LnhB+rkl+xbcaOKDK467A2j1YRkREGqEkG9QC4CAz625mzYBzgenVlpkOfLt4Nd9g4H13fzPBnEREJCcS+4rP3beb2VjgSaAMuMfdXzGzC4uvjwdmAKcBy4GPgO9EWPWEhFJOQl5yVZ7xy0uueckT8pNrXvKEwHNN7CIJERGRvaGZJEREJEhqUCIiEqRgGpSZ3WNm75hZjWOc6poWqb4plTLI9bxijovN7Hkz61fltZVm9hczW2RmL9X0/hTzHGpm7xdzWWRmV1V5LbWaRshzXJUc/2pmO8zsc8XX0qznAWb2v2a21MxeMbMf1rBMENtpxFwz304j5hnKdhol18y3VTNrYWbzzezlYp5X17BMENtpvdw9iB/gOGAA8NdaXj8NeILC2KnBwIvF58uAvwE9gGbAy8ChGed6NNCu+PuppVyLj1cCHQKp6VDgDzU8n2pN68uz2rLDgVkZ1fPzwIDi722A16rXJZTtNGKumW+nEfMMZTutN9cQttXitlde/L0p8CIwOMTttL6fYI6g3P1Z4N06FqltWqQoUyqlmqu7P+/uG4sP51EY35W6CDWtTao13c08RwEPJZVLXdz9TXf/U/H3D4ClQPV7FASxnUbJNYTtNGJNaxNcTavJZFstbnubiw+bFn+qXw0XxHZan2AaVARdgFVVHq8uPlfb86G4gMK/VEocmGlmC60whVPWjip+FfCEmfUuPhdkTc2sFTAM+M8qT2dSTzOrAL5E4V+nVQW3ndaRa1WZb6f15BnUdlpfTbPeVs2szMwWAe8AT7l78NtpTfJ0R93apkWKNF1SFszsBAo7/jFVnh7i7mvN7J+Ap8xsWfEIIgt/ojAH1mYzOw14nMLM8qHWdDgw192rHm2lXk8zK6fwh+df3X1T9ZdreEtm22k9uZaWyXw7rSfPoLbTKDUl423V3XcA/c1sP+AxMzvM3aue4w1qO61Nno6gapsWKcjpksysLzARGOHuG0rPu/va4n/fAR6jcEidCXffVPoqwN1nAE3NrAOB1pTCbCS7fGWSdj3NrCmFP04PuPvUGhYJZjuNkGsQ22l9eYa0nUapaVHm22ox1nvAbApHc1UFs53WKauTXzX9ABXUfkL/q+x6Um9+8fl9gBVAdz49qdc741y7UZgd4+hqz7cG2lT5/XkKM75nlef+fDpYexDwj2J9U69pXXkWX9+Xwnmq1lnVs1ib+4Fb61gmiO00Yq6Zb6cR8wxiO42SawjbKtAR2K/4e0vgj8DpIW6n9f0E8xWfmT1E4WqdDma2GvgZhZN7eB3TInktUyplnOtVQHvgLjMD2O6FGYM7UTjchsKG8KC7/0+GeZ4FfN/MtgNbgHO9sJWmWtMIeQJ8HZjp7h9WeWuq9QSGAN8C/lL8fh/g3yn8oQ9tO42SawjbaZQ8g9hOI+YK2W+rnwfus8JNY5sAj7j7HyzCNHNZ/D2ti6Y6EhGRIOXpHJSIiDQialAiIhIkNSgREQmSGpSIiARJDUpERIKkBiUiIkFSgxIRkSCpQYnkgJl9t3gvoe9EWLaHmU0ys9+nkZtIUtSgRPLh/wAnAmfXt6AXbpVwQfIpiSRLDUokBmbW1cymmdnrZrbCzO4ws+Z7uK6WZjanOFVNyYsUbp3wYpXl+pjZH6r9/FMt62xmZs+aWTDTm4nURw1KZC9ZYYK1qcDj7n4QhVtBtAR+uYerPB+Y6oVbJpSUU5j0c9/SE+7+F3c/vdrPOzWt0As3n3sGGLmHOYmkTg1KZO+dCGx1999B5b14/g34dvHeQbvrPGBa6YGZNaEwAem3ga9XO7L6DDNrb2bjgS+Z2RVVXnq8uG6RXNDhvsje6w0srPqEu28ys5XAgcCiqCsys2ZAD3dfWeXpE4HF7r7SzF4uPn6qtnV44b5OF9bw0l+BI6LmIpI1HUGJ7D2j5ruO1nR30vp0AN6r9tx5fHrzu4fYw6Og4pHdx2bWZk/eL5I2NSiRvfcKcHjVJ8ysLYV7APUys38xs5PM7N+Lr91iZl3MrLmZTa62ri1AiyrraQmMAG4qHpHdCIwoPr8nmgNb9/C9IqlSgxLZe88Arczs2wDFc0Q3A3e4+0PAwRS+cruhuPydwBgKtwV/uOqK3H0jUGZmpSZ1BvCEu3dz9wp37wb8FzB8d5M0s/bAOnf/ZHffK5IFNSiRvVS8u+vXgbPM7HVgA7DT3a8rXuG3H7CT4jlfd18OdObT225XNxM4pvj7ecBj1V5/DPjmHqR6AoU7qYrkgu6oKxIzMzuawrmiM4FjgdeBt4CR7v7j4jLnAp9391tqeP+XgEvc/Vsx5zUVuMLdX41zvSJJUYMSSZmZDQF+CFzg7h/Ussz5wH3VxkLtTcxmwLnufn8c6xNJgxqUiIgESeegREQkSGpQIiISJDUoEREJkhqUiIgESQ1KRESCpAYlIiJBUoMSEZEgqUGJiEiQ/j+SEmvReS1jngAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "a = 7.529\n",
    "b = 5.007\n",
    "gamma = 90\n",
    "fractional_coords = [[0,0,0 ], [0.5, 0.5, 0]]\n",
    "\n",
    "\n",
    "atoms_dummy = Atoms(['C', 'C'], scaled_positions=fractional_coords, cell=[a, b, 100, 90, 90, gamma])\n",
    "\n",
    "result = get_reflexes(atoms_dummy)\n",
    "\n",
    "plot_result_gixd(result, reflections_exp)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "d072719e",
   "metadata": {},
   "source": [
    "The missing reflexes can be obtained by making the two sites distinguishable, for example by using two different basis atoms (for example C and S). This means the actual unit cell has two distinguishable sites."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "bed5170d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAsHklEQVR4nO3de7wWZbn/8c/FkvMCJWBrgATkAUUOgSKKB9RMNJXyp6JZbdPiZUq7vTVSt2Vqun+amuYpQjSNnxraVmG3cWvKBpODIIWWgEqIcTAERBEFFLh+f8zzLBeLdRhgDves9X2/XuvFeuaZZ66La2a4mWfmvm9zd0RERELTLO8EREREaqMGSkREgqQGSkREgqQGSkREgqQGSkREgqQGSkREgpRaA2Vm95vZO2b21zreNzO7w8wWm9krZjYwrVxERKR40ryCegAYXs/7JwP7l35GAb9MMRcRESmY1Bood38eeLeeVUYAv/HIbGAvM/tsWvmIiEix7JFj7K7Asmqvl5eWvV1zRTMbRXSVRdu2bQf17t07kwSbhA0rAVj67mYAenymJVR2STXk0hUrolhdu6Yap2hUl3iSrFNeNQ9tX+edz7x589a4e+eay/NsoKyWZbWOu+Tu44BxAIceeqi/9NJLaebVtNwa7YYLXxkKwH39ZsBlK1MNeeFVV0Wxbrgh1ThFo7rEk2Sd8qp5aPs673zM7K3alufZQC0H9q32uhuQ7r+MUqdvdl2cXawRIzKLVSSqSzxJ1imvmoe2r0PLp8zSHCzWzHoAv3f3Q2p578vAaOAU4HDgDncf3NA2dQWVsFtruZC9TAMIi0h2zGyeux9ac3lqV1Bm9ggwDOhkZsuBnwDNAdx9LDCFqHFaDHwEfCutXKRhSz+qBKBHmw3pxwrs+/dQqC7x6B5U8kLLpyy1Bsrdz23gfQcuSSu+7JyfLu4PlO5BpR3rnnuiWIF8/x4K1SWepOr0ySefMP+VV+i4555sXL8+idRiW/veewCZx61LVvm0atWKbt260bx581jr53kPSkQkN8uXL+eAXr1o3bYtPbt1yzR2aFcsWeTj7qxdu5bly5fTs2fPWJ/RUEci0iRt2rSJ1m3bYlbbA8WSNDOjY8eObNq0KfZn1ECJSJOlxilbO1tvNVAiIhIk3YMSAL6z7+vZxTr77MxiFYnqEk+SderUoUNi29rVuDfccAMPP/wwFRUVNGvWjF/96lfce++9XHrppRx88MG7HatHjx689NJLdOrUqc51fj1+PGPGjKl6feSRRzJz5szdjr27Uu0HlQb1g0qY+kFJE7Vw4UIOOuigXHOYNWsWl156KdOmTaNly5asWbOGjz/+mC5dkhtuLE4DVVlZyYYN6XcxgdrrXlc/KH3FJwAs2tCeRRvaZxNryRIWLVmSSawiUV3iSbJOGzdvZuPmaBxKZl4T/YctqZ+Z1zQY9+2336ZTp060bNkSgE6dOtGlSxeGDRtG+T/ilZWVXH755QwaNIgvfvGLzJkzh2HDhtGrVy8mT54MwAMPPMDo0aOrtn/qqacybdq0HeJ+5StfYdCgQfTp04dx48YBcMUVV7Bx40b69e/PeeedVxUToifvxowZwyGHHELfvn2ZOHEiANOmTWPYsGGceeaZ9O7dm/POO480LnbUQAkANy/py81L+mYT6777uPm++zKJVSSqSzxJ1mnVmjWsWrMmkW3tStwvfelLLFu2jAMOOICLL76Y6dOn77Duhx9+yLBhw5g3bx7t2rXjRz/6EX/4wx944oknuPrqq3cq7v3338+8efN46aWXuOOOO1i7di033ngjrVq1YvKUKTz00EPbrf/4448zf/58Xn75ZZ599lnGjBnD229H43n/+c9/5vbbb2fBggUsWbKEGTOS70OpBkpEJCeVlZXMmzePcePG0blzZ0aOHMkDDzyw3TotWrRg+PBoar2+ffty7LHH0rx5c/r27cvSpUt3Kt4dd9xB//79GTJkCMuWLeONN96od/0XXniBc889l4qKCvbee2+OPfZY5s6dC8DgwYPp1q0bzZo1Y8CAATudSxx6SEJEJEcVFRUMGzaMYcOG0bdvXx588MHt3m/evHnV49nNmjWr+jqwWbNmbNmyBYA99tiDbdu2VX2mtr5G06ZN49lnn2XWrFm0adOGYcOGNdgnqb6v7cp5lP8O5VySpCsoEZGcvPbaa9tdxcyfP5/Pfe5zO72dHj16MH/+fLZt28ayZcuYM2fODuu8//77dOjQgTZt2rBo0SJmz55d9d4ezZvzySef7PCZY445hokTJ7J161ZWr17N888/z+DBDY7pnRhdQYmIABx5TfSToQ0bNvC9732P9957jz322IP99tuPcePGceaZZ+7UdoYOHUrPnj3p27cvhxxyCAMHDtxhneHDhzN27Fj69evHgQceyJAhQ6reO/drX+PkE0/k8MGDt7sP9dWvfpVZs2bRv39/zIyf/exn7LPPPixatGjX/9I7QY+ZN3Wlx8znr4/6ZQxovy71x8znlw7uAZoZeTuqSzxJ1WnhwoV8rkcPANq0br27ae2UjzZuzCVuXbLMZ2ceM9cVlAClhimrWPoHuFaqSzxJ1imvBiKUhqkstHzKdA9KgOgKqnwVlXqsRYuq/hcsn1Jd4kmyTh9t3Fh19ZClvOLWJbR8ytRACQB3Lj2YO5fu/rAqsWJNmMCdEyZkEqtIVJd4kqzTO+++yzvvvpvItooQty6h5VOmBkpERIKkBkpERIKkBkpERIKkBkpEJCcVFRUMGDCg6ufGG29MNd7kyZNTjzFt2rTEpurQY+YCwJhef8ku1oUXZharSFSXeJKs0971TEGRpnLc1q1bM3/+/ExibtmyhdNPP53TTz+9znySMG3aNCorKznyyCN3e1u6ghIAeleup3fl+mxi9epF7169MolVJKpLPEnWqXXLlrSuNqZcVuqL+/7773PggQfy2muvAXDuuedy7733AtHgspdddhkDBw7khBNOYPXq1QD87W9/Y/jw4QwaNIijjz66aqSH888/n0svvZTjjjuOyy+/fLtpOc4//3y++93vctxxx9HnoIOYM3s2F1xwAQcddBDnn39+VT7PPPMMRxxxBAMHDuSss86qmjeqR48e/OQnP2HgwIH07duXRYsWsXTpUsaOHcttt93GgAED+OMf/7hbddIVlAAwe11nAIZ0WJ1+rJdfjmL17596rCJRXeJJsk4bPvoIgMo2bbjwqqt2eP9LQ4cy8pRT2Lh5M6Ovu26H908//nhGnHAC69av5wc33bTde/fdcEODcTdu3MiAAQOqll955ZWMHDmSu+66i/PPP5/vf//7rFu3ju985ztANPXGwIEDufXWW7nuuuu49tprueuuuxg1ahRjx45l//3358UXX+Tiiy9m6tSpALz++us8++yzVFRU7DBS+rp165g6dSoTH3uM0047jRkzZjB+/HgOO+ww5s+fT7du3bj++ut59tlnadu2LTfddBM///nPq6b56NSpE3/605+45557uOWWWxg/fjwXXXQRlZWV/OAHP2ig+g1TAyUA3LvsACCbBureRx+NYukf4u2oLvEkWac166IRVCrbtNntbe1K3Lq+4jvxxBN57LHHuOSSS3i51CBDNIL5yJEjAfj617/OGWecwYYNG5g5cyZnnXVW1Xqby5MwAmeddRYVFRW15nHaaadhZnTp1o2OHTvSt280J1yfPn1YunQpy5cvZ8GCBQwdOhSAjz/+mCOOOKLq82eccQYAgwYN4vHHH9+VUtRLDZSICPVf8bRu2bLe9zu0b1/v+ztr27ZtLFy4kNatW/Puu+/SrVu3WtczM7Zt28Zee+1V572stm3b1hmn+tQdLap95VieyqOiooITTzyRRx55pN7Pa7oNEZEm4rbbbuOggw7ikUce4YILLqiaCmPbtm387ne/A+Dhhx/mqKOOon379vTs2ZPHHnsMiOZwqn7VtTuGDBnCjBkzWLx4MQAfffQRr7/+er2fadeuHR988EEi8dVAiYjkpHwPqvxzxRVX8PrrrzN+/HhuvfVWjj76aI455hiuv/56ILoaevXVVxk0aBBTp06tuhf00EMPcd9999G/f3/69OnDpEmTEsmvc+fOPPDAA5x77rn069ePIUOGNDjVxmmnncYTTzyRyEMSmm6jqStNt3HhK9F3zPf1m5H6dBvlm9FJfiXSGKgu8SRVp4ULF9K6fXsAenTtutt57YylK1bsUtzKysqqp+hCyGdXaLoN2Wk/3i+ZrwRixbr44sxiFYnqEk+Sdfps586JbasIcesSWj5laqAEgB5tkv9fWZ2xMv7falGoLvEkWaeWLVoktq0s4qZx9QT51aEhugclAExfuzfT1+6dTaw5c5g+Z04msYpEdYknyTqt37CBDz78MJFt7YwPPvwwl7h1ySqfnb2lpCsoAeA3K/YD4NiOq9KPVbqBe+zgwanHKhLVJZ6k6tSqVSuWr1hB67ZtaVfPo9hpWPveewCZx61LFvm4O2vXrqVVq1axP6MGSkSapG7duvHfzz5Lxz33ZFNCj0XHVW4QNq7PZnixhmSVT6tWrers01UbNVAi0iQ1b96c/3rhBSD7JydDe2IztHzKdA9KRESCpAZKRESClGpHXTMbDvwCqADGu/uNNd7fE/h/QHeirxtvcfdf17dNddRNWKmj7j82Rzcu92m5KfWOuv8oTRGwT6B9L/KiusSTZJ3yqnlo+zrvfDLvqGtmFcDdwInAcmCumU129wXVVrsEWODup5lZZ+A1M3vI3T9OKy+p3T4tN2UXK5CTMjSqSzxJ1imvmoe2r0PLpyzNr/gGA4vdfUmpwfktMKLGOg60MzMDKoF3geSHxJUGPb26C0+v7pJNrBde4OnSzWn5lOoST5J1yqvmoe3r0PIpS/Mpvq7AsmqvlwOH11jnLmAysBJoB4x09201N2Rmo4BRAN27d08l2abu0bd7AnBS55Xpx3rqqSjWUUelHqtIVJd4kqxTXjUPbV+Hlk9ZmldQVsuymjc3TgLmA12AAcBdZtZ+hw+5j3P3Q9390M6BXoqKiEiy0myglgP7VnvdjehKqbpvAY97ZDHwJtA7xZxERKQg0myg5gL7m1lPM2sBnEP0dV51fwdOADCzvYEDgSUp5iQiIgWR2j0od99iZqOBp4keM7/f3V81s4tK748Ffgo8YGZ/IfpK8HJ3X5NWTiIiUhyasLCpK/WDWvdJNNx+h+Yfp94Pal1pvK8O7Xe43dikqS7xJFmnvGoe2r7OOx9NWCj16tA8u65noZyUoVFd4kmyTnnVPLR9HVo+ZRrqSACYtGpfJq3at+EVk4j13HNMeu65TGIVieoST5J1yqvmoe3r0PIpUwMlAExe1Z3Jq7LpYzZ56lQmT52aSawiUV3iSbJOedU8tH0dWj5laqBERCRIaqBERCRIaqBERCRIaqBERCRI6gfV1JX6QW3cWgFA64qtqfeD2rh5cxSrZctU4xSN6hJPknXKq+ah7eu881E/KKlX64qt2cUK5KQMjeoST5J1yqvmoe3r0PIp01d8AsDElT2YuLJHNrGmTGHilCmZxCoS1SWeJOuUV81D29eh5VOmBkoAeGZNV55Z0zWbWDNm8MyMGZnEKhLVJZ4k65RXzUPb16HlU6YGSkREgqQGSkREgqQGSkREgqQGSkREgqR+UE1dqR/UdlLuByUiUl1d/aB0BSUiIkFSAyUAPLj88zy4/PPZxHrySR588slMYhWJ6hJPknXKq+ah7evQ8ilTAyUAPP/uPjz/7j7ZxJo7l+fnzs0kVpGoLvEkWae8ah7avg4tnzI1UCIiEqQGGygzG21mHbJIRkREpCzOFdQ+wFwze9TMhptZLY99iYiIJKvBBsrdfwTsD9wHnA+8YWb/YWbZ3FGXTLRqtpVWzbIZ0bxVy5a0CnT05DypLvEkWae8ah7avg4tn7LY/aDMrD/wLWA48L/AEOAP7v7D9NLbkfpBJUz9oEQkZ7s8H5SZ/Qvwz8AaYDwwxt0/MbNmwBtApg2UiIg0DXEmLOwEnOHub1Vf6O7bzOzUdNKSrI37+wEAjOr+evqxJk6MYo0cmXqsIlFd4kmyTnnVPLR9HVo+ZXEekuhZs3EyswkA7r4wlawkcy++15kX3+ucTaxXXuHFV17JJFaRqC7xJFmnvGoe2r4OLZ+yOA1Un+ovzKwCGJROOiIiIpE6Gygzu9LMPgD6mdn60s8HwDvApMwyFBGRJqnOBsrd/6+7twNudvf2pZ927t7R3a/MMEcREWmC6nxIwsx6u/si4DEzG1jzfXf/U6qZSab22uPj7GK1a5dZrCJRXeJJsk551Ty0fR1aPmV19oMys3HuPsrM/reWt93dj083tdqpH1TC1A9KRHK20/2g3H1U6c/j0kxMRESkNnEGiz3LzNqVfv+RmT1uZl9IPzXJ0h1vHsQdbx6UTawJE7hjwoRMYhWJ6hJPknXKq+ah7evQ8imL01H3x+7+mJkdBZwE3AKMBQ5PNTPJ1MsffCa7WIsWZRarSFSXeJKsU141D21fh5ZPWZx+UOURRL8M/NLdJwEt0ktJREQkXgO1wsx+BZwNTDGzljE/R2l6jtfMbLGZXVHHOsPMbL6ZvWpm0+OnLiIijVmchuZs4GlguLu/B3wGGNPQh0ojTtwNnAwcDJxrZgfXWGcv4B7gdHfvA5y1M8mLiEjj1eA9KHf/yMwmAXubWffS4jhfWA4GFrv7EgAz+y0wAlhQbZ2vAY+7+99Lsd7ZmeQlOXu32JhdrE6dMotVJKpLPEnWKa+ah7avQ8unrMH5oMzse8BPgFXAttJid/d+DXzuTKKrrm+XXn8DONzdR1db53agOdF4f+2AX7j7b2rZ1ihgFED37t0HvfXWWzVXkV2lflAikrNdng8K+D5woLuv3dmYtSyr+S/fHkQDz54AtAZmmdlsd99uzgd3HweMg6ij7k7mISIiBRSngVoGvL8L214O7FvtdTdgZS3rrHH3D4EPzex5oD+Q/qREsp2f/e0QAH74+b+mH2v8+CjWt7+deqwiUV3iSbJOedU8tH0dWj5lcRqoJcA0M/tvYHN5obv/vIHPzQX2N7OewArgHKJ7TtVNAu4ysz2IHl0/HLgtZu6SoNc+3DO7WG++mVmsIlFd4kmyTnnVPLR9HVo+ZXEaqL+XflqwE/2f3H2LmY0megKwArjf3V81s4tK749194Vm9j/AK0T3t8a7e/r/hRcRkeDFeYrvWgAza1v6Ki42d58CTKmxbGyN1zcDN+/MdkVEpPGLMxbfEWa2AFhYet3fzO5JPTMREWnS4nzFdzvRGHyTAdz9ZTM7Js2kJHufa70hu1hdumQWq0hUl3iSrFNeNQ9tX4eWT1mcflAvuvvhZvZnd/9CadnL7t4/kwxr0HxQCVM/KJHszLxmx2VH1rKsidmdflDLzOxIwM2sBfAvlL7uExGRnTDr2h2XqYGqU5yx+C4CLgG6EvVbGgBcnGJOkoPr3ujPdW9kc1F83d13c93dd2cSq0hUl3iSrFNeNc/yfIsj1GMvzhXUge5+XvUFZjYUmJFOSpKHtzZWZhdrZc3+2gKqS1xJ1imvmmd5vsUR6rEX5wrqzpjLREREElPnFZSZHQEcCXQ2s0urvdWeqOOtiIhIaur7iq8FUFlap1215euBM9NMSkREpM4Gyt2nA9PN7AF31/wWjdyBbXdlPOBdjNWzZ2axikR1iSfJOuVV8yzPtzhCPfbi9IM6APgB0INqDZq7H59qZnVQP6iEqR+USHZ0vtVqd/pBPQaMBcYDW5NOTEREpDZxGqgt7v7L1DORXP37ooEA/EfvP6Uf67ZoRpX/+Ld/Sz1Wkagu8SRZp7xqnuX5Fkeox16cBuq/zOxi4Am2nw/q3dSyksyt+rh1drHWrMksVpGoLvEkWae8ap7l+Vav0tBLq5a88enrgEa2iNNA/XPpzzHVljnQK/l0REQkM+Whl9YPLb1+uFgNlLuH+XiHiIg0avV11D3e3aea2Rm1ve/uj6eXloiINHX1XUEdC0wFTqvlPQfUQDUi/dtld0uxf+/emcUqEtUlniTrlFfNszzf4ggtn7IG+0GFRv2gEqZ+GSLZCe18CySfuvpBxRksVkREJHNqoASAyxYcxmULDssm1o03ctmNN2YSq0hUl3iSrFNeNc/yfIsjtHzK4jxmLk3Ae1taZBfrgw8yi1Ukqks8SdYpr5pneb7FEVo+ZQ1eQZlZGzP7sZndW3q9v5mdmn5qIiLSlMX5iu/XRCNIHFF6vRy4PrWMREREiNdAfd7dfwZ8AuDuG4FaHv0QiWHmNbB+afQz85qqoVZERGqKcw/qYzNrTdT3CTP7PNXG5JPG4fC9VmcTaNa1HN7sgNLvD0d/BjS0Sp4O79cv7xQKIck65VXzzM63mELLpyzOfFBfAq4CDgaeAYYC57v7tNSzq4X6QSUs634QgfS7EMlFaMd/IPns8nxQ7v6Mmc0DhhB9tfd9d9ewyyIikqoGGygzmww8Akx29w/TT0nycMlfhwBw9yGzG1WsIrnkuusAuPvqq3POJGxJ1imvmod2DoSWT1mchyRuBY4GFpjZY2Z2ppm1SjkvydimbRVs2lbR6GIVyabNm9m0Wbd3G5JknfKqeWjnQGj5lMX5im86MN3MKoDjge8A9wPtU85NRESasFgjSZSe4jsNGAkMBB5MMykREZE496AmAocD/wPcDUxz921pJyYiIk1bnCuoXwNfc/etaScj+TnmM/9olLGK5JjDwhusM0RJ1imvmod2DoSWT1md/aBCnVFX/aASpn5QItkJ7fgPJJ9d6QelGXVFRCQ3dTZQ7v6T0q/Xufub1d8zs55xNm5mw4FfABXAeHevdeIVMzsMmA2MdPffxdm2JOvCV4YCcF+/GY0qVpFceNVVANx3ww05ZxK2JOuUV81DOwdCy6csTj+o/6xlWYONSOmx9LuBk4mGSTrXzA6uY72bgKdj5CIiIk1EnVdQZtYb6APsWeM+VHsgTkfdwcBid19S2t5vgRHAghrrfY+oEdQdYhERqVLfPagDgVOBvdj+PtQHRJ11G9IVWFbt9XKix9WrmFlX4KtEHYDrbKDMbBQwCqB79+4xQouISNHVdw9qEjDJzI5w91m7sO3a5oyq+XjI7cDl7r7VrO4pptx9HDAOoqf4diEXEREpmDj9oP5sZpcQfd1X9dWeu1/QwOeWA/tWe90NWFljnUOB35Yap07AKWa2xd2fjJGXJOhLnVY0ylhF8qWhQ/NOoRCSrFNeNQ/tHAgtn7I4DdQEYBFwEnAdcB6wMMbn5gL7l574WwGcA3yt+gruXvU0oJk9APxejVM+RnZZ2ihjFcnIU07JO4VCSLJOedU8tHMgtHzK4jzFt5+7/xj40N0fBL4M9G3oQ+6+BRhN9HTeQuBRd3/VzC4ys4t2J2lJ3satFWzcms1oxlnGKpKNmzezUaOZNyjJOuVV89DOgdDyKYtzBfVJ6c/3zOwQ4B9Ajzgbd/cpwJQay8bWse75cbYp6Rj9ajQfTBb9ILKMVSSjS3MTqR9U/ZKsU141D+0cCC2fsjgN1Dgz6wD8GJgMVAKaUU1ERFIVZz6o8aVfpwO90k1HREQkUl9H3Uvr+6C7/zz5dERERCL1XUG1yywLERGRGurrqHttlolIvk7f+++NMlaRnH788XmnUAhJ1imvmod2DoSWT1mcGXUPAH4J7O3uh5hZP+B0d78+9ewkMyP2XtbwSgWMVSQjTjgh7xQKIck65VXz0M6B0PIpi9MP6l7gSkqPm7v7K0SdbqURWfdJC9Z90qLRxSqSdevXs279+rzTCF6Sdcqr5qGdA6HlUxbnMfM27j6nxlh5W1LKR3Lyg4XRWL1Z9IPIMlaR/OCmmwD1g2pIknXKq+ahnQOh5VMW5wpqjZl9ntJAr2Z2JvB2qlmJiEiTF+cK6hKikcR7m9kK4E2i8fhERERSE6ej7hLgi2bWluiKayMwEngr5dxERKQJq/MrPjNrb2ZXmtldZnYi8BHwz8Bi4OysEhQRkaapviuoCcA6YBbRDLo/BFoAX3H3+emnJlk6+7NvNspYRXL2ySfnnUK4Zl5T9evZvTYkttm8ah7aORBaPmX1NVC93L0vgJmNB9YA3d39g0wyk0yd1LnmXJKNI1aRnHTUUXmnEK5Zn44bcFLVb7/Y7c3mVfPQzoHQ8imrr4EqT7NBaUr2N9U4NV7/2BxNlrxPy02NKlaR/GP1agD26dw550zCVnX8JLGtnGoe2jkQWj5l9TVQ/c2s3IPNgNal1wa4u7dPPTvJzFWvDQKy6QeRZawiuer22wH1g2pI1fGTxLZyqnlo50Bo+ZTVNxZfeNMriohIkxGno66IiEjm1ECJiEiQ1ECJiEiQ4gx1JE3AN7subpSxiuSbI0bknUIhJHn85FXz0M6B0PIpUwMlABzbcVWjjFUkxw4enHcKhZDk8ZNXzUM7B0LLp0wNlACw9KNKAHq0Sa6XfgiximTpihUA9OjaNedMwlZ1/CSxrZxqHto5EFo+ZWqgBICfLu4PZNMPIstYRfLTe+4B1A+qIVXHTxLbyqnmoZ0DoeVTpockREQkSGqgREQkSGqgREQkSGqgREQkSHpIQgD4zr6vN8pYRfKdszUPaBxJHj951Ty0cyC0fMrUQAkAQzqsbpSximRI//55p1AISR4/edU8tHMgtHzK1EAJAIs2RLOn9K5c38CaxYpVJIuWLAGgd69eOWcStqrjJ4lt5VTz0M6B0PIpUwMlANy8pC+QTT+ILGMVyc33RT171A+qflXHTxLbyqnmoZ0DoeVTpockREQkSGqgREQkSGqgREQkSGqgREQkSKk+JGFmw4FfABXAeHe/scb75wGXl15uAL7r7i+nmZPU7ns9FjTKWEXyvW98I+8UCiHJ4yevmod2DoSWT1lqDZSZVQB3AycCy4G5ZjbZ3atX4k3gWHdfZ2YnA+OAw9PKSeo2oP26RhmrSAb0TuLB6cYvyeMnr5qHdg6Elk9ZmldQg4HF7r4EwMx+C4wAqhood59Zbf3ZQLcU85F6zF/fAcjmQM0yVpHMX7QIUEPVkKrjJ4lt5VTz0M6B0PIpS/MeVFdgWbXXy0vL6nIh8FRtb5jZKDN7ycxeWr06zB7PRXfn0oO5c+nBjS5Wkdw5YQJ3TpiQdxrBS/L4yavmoZ0DoeVTlmYDZbUs81pXNDuOqIG6vLb33X2cux/q7od27tw5wRRFRCRUaX7FtxzYt9rrbsDKmiuZWT9gPHCyu69NMR+RsK1fGv0585pPlx15TS0rijQNaTZQc4H9zawnsAI4B/ha9RXMrDvwOPANdw9zOF2RrKx/K/pz1sOfLlMDJU1Yag2Uu28xs9HA00SPmd/v7q+a2UWl98cCVwMdgXvMDGCLux+aVk4iIlIcqfaDcvcpwJQay8ZW+/3bwLfTzEHiGdPrL40yVpGoLvEkWacxF16Y2LZ2Km5g+zq0fMo0mrkA2Q6zH9qQ/qFQXeJJsk55TW0S2r4OLZ8yNVACwOx10dORWUxclmWsIlFd4qmqUxLbejkauCbriQtD29eh5VOmBkoAuHfZAUA2B2iWsYpEdYmnqk5JbOvRR6NtZdxAhbavQ8unTA2UiDRt65du/2g/6OnJQKiBEpGmbf1b2z/aD2qgAqHpNkREJEhqoEREJEj6ik8A+PF+2U3DlWWsIlFd4kmyTj+++GK4/47Ethc7bmD7OrR8ytRACQA92mxolLGKRHWJJ8k69ejaFXKoe2j7OrR8yvQVnwAwfe3eTF+7d6OLVSSqSzxJ1mn6nDm51Dy0fR1aPmW6ghIAfrNiPwCO7biqUcUqEtUlnqo6JbGtSZNgxX6Z1zy0fR1aPmW6ghIRkSDpCioNNTv9gfpViIjsJDVQaZh17Y7L1ECJiOwUfcUnIiJB0hWUAHDDgfMaZawiUV3iSbJON/zrv8K9P09se7HjBravQ8unTA2UALBPy02NMlaRqC7xJFmnfTp3hhzqHtq+Di2fMn3FJwA8vboLT6/u0uhiFYnqEk+SdXr6hRdyqXlo+zq0fMp0BSUAPPp2TwBO6ryyUcUqEtUlnqo6JbGtp56Ct3tmXvPQ9nVo+ZTpCkpERIKkBkpEwlZbv0JpEtRAiUjYautXKE2CGigREQmSHpIQAG45aG6jjFUkqks8Sdbplssvh3tuSWx7seMGtq9Dy6dMDZQA0KH5x40yVpGoLvEkWacO7dtDDnUPbV+Hlk+ZvuITACat2pdJq/ZtdLGKRHWJJ8k6TXruuVxqHtq+Di2fMjVQAsDkVd2ZvKp7o4tVJKpLPEnWafLUqbnUPLR9HVo+ZWqgREQkSGqgREQkSGqgREQkSGqgREQkSHrMXAC4q8/sRhmrSFSXeJKs011XXw2/aJvY9mLHDWxfh5ZPmRooAaB1xdZGGatIVJd4kqxT65YtIYe6h7avQ8unTF/xCQATV/Zg4soejS5Wkagu8SRZp4lTpuRS89D2dWj5lKmBEgCeWdOVZ9Z0bXSxikR1iSfJOj0zY0YuNQ9tX4eWT5kaKBERCVKqDZSZDTez18xssZldUcv7ZmZ3lN5/xcwGppmPiIgUR2oNlJlVAHcDJwMHA+ea2cE1VjsZ2L/0Mwr4ZVr5iIhIsaR5BTUYWOzuS9z9Y+C3wIga64wAfuOR2cBeZvbZFHMSEZGCMHdPZ8NmZwLD3f3bpdffAA5399HV1vk9cKO7v1B6/Rxwubu/VGNbo4iusAAOAf6aStLJ6wSsyTuJGJRn8oqSa1HyhOLkWpQ8IZxcP+funWsuTLMflNWyrGZrGGcd3H0cMA7AzF5y90N3P730FSVX5Zm8ouRalDyhOLkWJU8IP9c0v+JbDlSfYKQbsHIX1hERkSYozQZqLrC/mfU0sxbAOcDkGutMBr5ZeppvCPC+u7+dYk4iIlIQqX3F5+5bzGw08DRQAdzv7q+a2UWl98cCU4BTgMXAR8C3Ymx6XEopp6EouSrP5BUl16LkCcXJtSh5QuC5pvaQhIiIyO7QSBIiIhIkNVAiIhKkYBooM7vfzN4xs1r7ONU3LFJDQyrlkOt5pRxfMbOZZta/2ntLzewvZjbfzF6q7fMZ5jnMzN4v5TLfzK6u9l5mNY2R55hqOf7VzLaa2WdK72VZz33N7H/NbKGZvWpm369lnSCO05i55n6cxswzlOM0Tq65H6tm1srM5pjZy6U8r61lnSCO0wa5exA/wDHAQOCvdbx/CvAUUd+pIcCLpeUVwN+AXkAL4GXg4JxzPRLoUPr95HKupddLgU6B1HQY8Ptalmda04byrLHuacDUnOr5WWBg6fd2wOs16xLKcRoz19yP05h5hnKcNphrCMdq6dirLP3eHHgRGBLicdrQTzBXUO7+PPBuPavUNSxSnCGVMs3V3We6+7rSy9lE/bsyF6Omdcm0pjuZ57nAI2nlUh93f9vd/1T6/QNgIVBzjoIgjtM4uYZwnMasaV2Cq2kNuRyrpWNvQ+ll89JPzafhgjhOGxJMAxVDV2BZtdfLS8vqWh6KC4n+p1LmwDNmNs+iIZzydkTpq4CnzKxPaVmQNTWzNsBw4D+rLc6lnmbWA/gC0f9OqwvuOK0n1+pyP04byDOo47ShmuZ9rJpZhZnNB94B/uDuwR+ntSnSlO91DYsUa7ikPJjZcUQn/lHVFg9195Vm9k/AH8xsUekKIg9/IhoDa4OZnQI8STSyfKg1PQ2Y4e7Vr7Yyr6eZVRL9w/Ov7r6+5tu1fCS347SBXMvr5H6cNpBnUMdpnJqS87Hq7luBAWa2F/CEmR3i7tXv8QZ1nNalSFdQdQ2LFORwSWbWDxgPjHD3teXl7r6y9Oc7wBNEl9S5cPf15a8C3H0K0NzMOhFoTYlGI9nuK5Os62lmzYn+cXrI3R+vZZVgjtMYuQZxnDaUZ0jHaZyaluR+rJZivQdMI7qaqy6Y47Reed38qu0H6EHdN/S/zPY39eaUlu8BLAF68ulNvT4559qdaHSMI2ssbwu0q/b7TKIR3/PKcx8+7aw9GPh7qb6Z17S+PEvv70l0n6ptXvUs1eY3wO31rBPEcRoz19yP05h5BnGcxsk1hGMV6AzsVfq9NfBH4NQQj9OGfoL5is/MHiF6WqeTmS0HfkJ0cw+vZ1gkr2NIpZxzvRroCNxjZgBbPBoxeG+iy22IDoSH3f1/cszzTOC7ZrYF2Aic49FRmmlNY+QJ8FXgGXf/sNpHM60nMBT4BvCX0vf7AP9O9A99aMdpnFxDOE7j5BnEcRozV8j/WP0s8KBFk8Y2Ax51999bjGHm8vj3tD4a6khERIJUpHtQIiLShKiBEhGRIKmBEhGRIKmBEhGRIKmBEhGRIKmBEhGRIKmBEhGRIKmBEikAM/t2aS6hb8VYt5eZ3Wdmv8siN5G0qIESKYb/AxwPnNXQih5NlXBh+imJpEsNlEgCzKybmU0yszfMbImZ3WVmLXdxW63NbHppqJqyF4mmTnix2np9zez3NX7+qY5ttjCz580smOHNRBqiBkpkN1k0wNrjwJPuvj/RVBCtgZ/t4iYvAB73aMqEskqiQT/3LC9w97+4+6k1ft6pbYMeTT73HDByF3MSyZwaKJHddzywyd1/DVVz8fwb8M3S3EE76zxgUvmFmTUjGoD0m8BXa1xZ7cDMOprZWOALZnZltbeeLG1bpBB0uS+y+/oA86ovcPf1ZrYU2A+YH3dDZtYC6OXuS6stPh54xd2XmtnLpdd/qGsbHs3rdFEtb/0VOCxuLiJ50xWUyO4zap91tLbZSRvSCXivxrLz+HTyu0fYxaug0pXdx2bWblc+L5I1NVAiu+9V4NDqC8ysPdEcQAeZ2b+Y2Qlm9u+l924zs65m1tLMJtTY1kagVbXttAZGALeUrshuBkaUlu+KlsCmXfysSKbUQInsvueANmb2TYDSPaJbgbvc/RHgAKKv3G4qrX83MIpoWvDfVt+Qu68DKsys3EidDjzl7t3dvYe7dwf+CzhtZ5M0s47Aanf/ZGc/K5IHNVAiu6k0u+tXgTPN7A1gLbDN3W8oPeG3F7CN0j1fd18MdOHTabdregY4qvT7ecATNd5/Avj6LqR6HNFMqiKFoBl1RRJmZkcS3Ss6AzgaeAP4BzDS3X9YWucc4LPuflstn/8CcKm7fyPhvB4HrnT315Lcrkha1ECJZMzMhgLfBy509w/qWOcC4MEafaF2J2YL4Bx3/00S2xPJghooEREJku5BiYhIkNRAiYhIkNRAiYhIkNRAiYhIkNRAiYhIkNRAiYhIkNRAiYhIkNRAiYhIkP4/4ouqN17bAfoAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "a = 7.529\n",
    "b = 5.007\n",
    "gamma = 90\n",
    "fractional_coords = [[0,0,0 ], [0.5, 0.5, 0]]\n",
    "\n",
    "\n",
    "atoms_dummy = Atoms(['C', 'S'], scaled_positions=fractional_coords, cell=[a, b, 100, 90, 90, gamma])\n",
    "\n",
    "result = get_reflexes(atoms_dummy)\n",
    "\n",
    "fig, ax = plot_result_gixd(result, reflections_exp)\n",
    "fig.tight_layout()\n",
    "fig.savefig('SOS/two_sites.png', dpi=300)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "973646d7",
   "metadata": {},
   "source": [
    "Distinguishable sites for the packing of SOS are the case for a herringbone-packing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "4459ca1d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_herringbone(a, b, c, theta, sos):\n",
    "    sos2 = sos.copy()\n",
    "    sos2.rotate(theta, 'z', sos2.get_center_of_mass())\n",
    "    sos2_pos = sos2.get_positions()\n",
    "    sos2_pos[:,0] += a/2\n",
    "    sos2_pos[:,1] += b/2\n",
    "    names_merged = np.append(sos.get_chemical_symbols(), sos2.get_chemical_symbols())\n",
    "    pos_merged = np.append(sos.get_positions(), sos2_pos, 0)\n",
    "    merged = Atoms(names_merged, positions=pos_merged, cell=[a, b, c, 90, 90, 90], pbc=True)\n",
    "    \n",
    "    return merged\n",
    "\n",
    "sos = read('SOS/HSOS.xyz')  # reading structure of HSOS\n",
    "\n",
    "a = 7.529\n",
    "b = 5.007\n",
    "\n",
    "s = get_herringbone(a, b, 1000, 90, sos)  # generating HB packing of HSOS\n",
    "s.write('SOS/HSOS_HB_from_GIXD.cif')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "b2ddddcf",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAArc0lEQVR4nO3de7xVdbnv8c/DkvsCJeBggATkBUUugSKCF9QsNJHyaEhWm7R4meL2bI3UbZmadjQzzVuEYCrHDG2jsNu4NWWDiSKXAkxAJMS4eAFEEUUQeM4fc05aLNZcawDj8hus7/v1Wq/WHGPM8Tw9Y7h+jMvv9zN3R0REJDQNsk5ARESkJmqgREQkSGqgREQkSGqgREQkSGqgREQkSGqgREQkSIk1UGb2gJm9a2Z/K7PezOwuM1tmZgvNrE9SuYiISP4keQX1IDC4lvVnAIcVf0YCv04wFxERyZnEGih3fx54r5ZNhgIPe8Es4CAz+2xS+YiISL4ckGHsDsDKKp9XFZe9VX1DMxtJ4SqL5s2b9+3WrVsqCdYLm9YAsOK9LQB0/kxjqGyfaMgVq1cXYnXokGicvFFdoomzTlnVPLRjnXU+8+bNW+fubasvz7KBshqW1TjukruPBcYCHHPMMT537twk86pfbi8chosWDgRgfM+ZcOWaRENedO21hVg335xonLxRXaKJs05Z1Ty0Y511Pmb2Zk3Ls2ygVgGHVPncEUj2L6OU9e0Oy9KLNXRoarHyRHWJJs46ZVXz0I51aPmUWJKDxZpZZ+CP7n50Deu+AowCzgSOA+5y93517VNXUDG7vYYL2Ss1gLCIpMfM5rn7MdWXJ3YFZWaPAoOANma2CvgJ0BDA3ccAUyk0TsuAj4HvJJWL1G3Fx5UAdG62KflYgd1/D4XqEo2eQcUvtHxKEmug3H14HesduDSp+LJnfrqsF1B8BpV0rPvuK8QK5P57KFSXaOKq06effsr8hQtpfeCBbN64MY7UIlv//vsAqcctJ618mjRpQseOHWnYsGGk7bN8BiUikplVq1ZxeNeuNG3enC4dO6YaO7QrljTycXfWr1/PqlWr6NKlS6TvaKgjEamXPvnkE5o2b45ZTS8US9zMjNatW/PJJ59E/o4aKBGpt9Q4pWtP660GSkREgqRnUALA9w5Zml6sr389tVh5orpEE2ed2rRqFdu+9jbuzTffzO9+9zsqKipo0KABv/nNb7j//vu54oorOOqoo/Y5VufOnZk7dy5t2rQpu81vx41j9OjROz8PGDCAF198cZ9j76tE+0ElQf2gYqZ+UFJPLV68mCOPPDLTHF566SWuuOIKpk+fTuPGjVm3bh1bt26lffv4hhuL0kBVVlayaVPyXUyg5rqX6welW3wCwJJNLVmyqWU6sZYvZ8ny5anEyhPVJZo467R5yxY2bymMQ8mL1xf+wRbXz4vX1xn3rbfeok2bNjRu3BiANm3a0L59ewYNGkTpH+KVlZVcddVV9O3bly9+8YvMnj2bQYMG0bVrV6ZMmQLAgw8+yKhRo3bu/6yzzmL69Om7xf3qV79K37596d69O2PHjgXg6quvZvPmzfTs1YsLLrhgZ0wovHk3evRojj76aHr06MHEiRMBmD59OoMGDeLcc8+lW7duXHDBBSRxsaMGSgC4bXkPblveI51Y48dz2/jxqcTKE9Ulmjjr9M66dbyzbl0s+9qbuF/60pdYuXIlhx9+OJdccgkzZszYbduPPvqIQYMGMW/ePFq0aMGPfvQj/vSnP/HEE09w3XXX7VHcBx54gHnz5jF37lzuuusu1q9fzy233EKTJk2YMnUqjzzyyC7bT5o0ifnz57NgwQKeffZZRo8ezVtvFcbz/utf/8qdd97JokWLWL58OTNnxt+HUg2UiEhGKisrmTdvHmPHjqVt27YMGzaMBx98cJdtGjVqxODBhan1evTowcknn0zDhg3p0aMHK1as2KN4d911F7169aJ///6sXLmS119/vdbtX3jhBYYPH05FRQXt2rXj5JNPZs6cOQD069ePjh070qBBA3r37r3HuUShlyRERDJUUVHBoEGDGDRoED169OChhx7aZX3Dhg13vp7doEGDnbcDGzRowLZt2wA44IAD2LFjx87v1NTXaPr06Tz77LO89NJLNGvWjEGDBtXZJ6m223alPEr/H0q5xElXUCIiGXnttdd2uYqZP38+n/vc5/Z4P507d2b+/Pns2LGDlStXMnv27N22+eCDD2jVqhXNmjVjyZIlzJo1a+e6Axo25NNPP93tOyeddBITJ05k+/btrF27lueff55+/eoc0zs2uoISEQEYcH3hJ0WbNm3isssu4/333+eAAw7g0EMPZezYsZx77rl7tJ+BAwfSpUsXevTowdFHH02fPn1222bw4MGMGTOGnj17csQRR9C/f/+d64Z/4xuccfrpHNev3y7Pob72ta/x0ksv0atXL8yMn//85xx88MEsWbJk7/9P7wG9Zl7fFV8zn7+x0C+jd8sNib9mPr94cvfWzMi7UF2iiatOixcv5nOdOwPQrGnTfU1rj3y8eXMmcctJM589ec1cV1ACFBumtGLpD3CNVJdo4qxTVg1EKA1TSWj5lOgZlACFK6jSVVTisZYs2fmvYPkn1SWaOOv08ebNO68e0pRV3HJCy6dEDZQAcPeKo7h7xb4PqxIp1oQJ3D1hQiqx8kR1iSbOOr373nu8+957sewrD3HLCS2fEjVQIiISJDVQIiISJDVQIiISJDVQIiIZqaiooHfv3jt/brnllkTjTZkyJfEY06dPj22qDr1mLgCM7vpKerEuuii1WHmiukQTZ53a1TIFRZJKcZs2bcr8+fNTiblt2zbOPvtszj777LL5xGH69OlUVlYyYMCAfd6XrqAEgG6VG+lWuTGdWF270q1r11Ri5YnqEk2cdWrauDFNq4wpl5ba4n7wwQccccQRvPbaawAMHz6c+++/HygMLnvllVfSp08fTjvtNNauXQvA3//+dwYPHkzfvn058cQTd470MGLECK644gpOOeUUrrrqql2m5RgxYgTf//73OeWUU+h+5JHMnjWLCy+8kCOPPJIRI0bszOeZZ57h+OOPp0+fPpx33nk7543q3LkzP/nJT+jTpw89evRgyZIlrFixgjFjxnDHHXfQu3dv/vznP+9TnXQFJQDM2tAWgP6t1iYfa8GCQqxevRKPlSeqSzRx1mnTxx8DUNmsGRdde+1u6780cCDDzjyTzVu2MOrGG3dbf/appzL0tNPYsHEjP7j11l3Wjb/55jrjbt68md69e+9cfs011zBs2DDuueceRowYweWXX86GDRv43ve+BxSm3ujTpw+33347N954IzfccAP33HMPI0eOZMyYMRx22GG8/PLLXHLJJUybNg2ApUuX8uyzz1JRUbHbSOkbNmxg2rRpTHz8cYYMGcLMmTMZN24cxx57LPPnz6djx47cdNNNPPvsszRv3pxbb72VX/7ylzun+WjTpg1/+ctfuO+++/jFL37BuHHjuPjii6msrOQHP/hBHdWvmxooAeD+lYcD6TRQ9z/2WCGW/hDvQnWJJs46rdtQGEGlslmzfd7X3sQtd4vv9NNP5/HHH+fSSy9lQbFBhsII5sOGDQPgm9/8Jueccw6bNm3ixRdf5Lzzztu53ZbSJIzAeeedR0VFRY15DBkyBDOjfceOtG7dmh49CnPCde/enRUrVrBq1SoWLVrEwIEDAdi6dSvHH3/8zu+fc845APTt25dJkybtTSlqpQZKRITar3iaNm5c6/pWLVvWun5P7dixg8WLF9O0aVPee+89OnbsWON2ZsaOHTs46KCDyj7Lat68edk4VafuaFTllmNpKo+KigpOP/10Hn300Vq/r+k2RETqiTvuuIMjjzySRx99lAsvvHDnVBg7duzgD3/4AwC/+93vOOGEE2jZsiVdunTh8ccfBwpzOFW96toX/fv3Z+bMmSxbtgyAjz/+mKVLl9b6nRYtWvDhhx/GEl8NlIhIRkrPoEo/V199NUuXLmXcuHHcfvvtnHjiiZx00kncdNNNQOFq6NVXX6Vv375MmzZt57OgRx55hPHjx9OrVy+6d+/O5MmTY8mvbdu2PPjggwwfPpyePXvSv3//OqfaGDJkCE888UQsL0louo36rjjdxkULC/eYx/ecmfh0G6WH0XHeEtkfqC7RxFWnxYsX07RlSwA6d+iwz3ntiRWrV+9V3MrKyp1v0YWQz97QdBuyx358aDy3BCLFuuSS1GLlieoSTZx1+mzbtrHtKw9xywktnxI1UAJA52bx/6usbKyU/7WaF6pLNHHWqXGjRrHtK424SVw9QXZ1qIueQQkAM9a3Y8b6dunEmj2bGbNnpxIrT1SXaOKs08ZNm/jwo49i2dee+PCjjzKJW05a+ezpIyVdQQkAD68+FICTW7+TfKziA9yT+/VLPFaeqC7RxFWnJk2asGr1apo2b06LWl7FTsL6998HSD1uOWnk4+6sX7+eJk2aRP6OGigRqZc6duzIfz37LK0PPJBPYnotOqpSg7B5YzrDi9UlrXyaNGlStk9XTdRAiUi91LBhQ/7zhReA9N+cDO2NzdDyKdEzKBERCZIaKBERCVKiHXXNbDDwK6ACGOfut1RbfyDw/4BOFG43/sLdf1vbPtVRN2bFjrpvbyk8uDy48SeJd9R9uzhFwMGB9r3IiuoSTZx1yqrmoR3rrPNJvaOumVUA9wKnA6uAOWY2xd0XVdnsUmCRuw8xs7bAa2b2iLtvTSovqdnBjT9JL1Yg/1GGRnWJJs46ZVXz0I51aPmUJHmLrx+wzN2XFxuc3wNDq23jQAszM6ASeA+If0hcqdPTa9vz9Nr26cR64QWeLj6cln9SXaKJs05Z1Ty0Yx1aPiVJvsXXAVhZ5fMq4Lhq29wDTAHWAC2AYe6+o/qOzGwkMBKgU6dOiSRb3z32VhcAvtx2TfKxnnqqEOuEExKPlSeqSzRx1imrmod2rEPLpyTJKyirYVn1hxtfBuYD7YHewD1m1nK3L7mPdfdj3P2YtoFeioqISLySbKBWAYdU+dyRwpVSVd8BJnnBMuANoFuCOYmISE4k2UDNAQ4zsy5m1gg4n8LtvKr+AZwGYGbtgCOA5QnmJCIiOZHYMyh332Zmo4CnKbxm/oC7v2pmFxfXjwF+CjxoZq9QuCV4lbuvSyonERHJD01YWN8V+0Ft+LQw3H6rhlsT7we1oTjeV6uWuz1urNdUl2jirFNWNQ/tWGedjyYslFq1aphe17NQ/qMMjeoSTZx1yqrmoR3r0PIp0VBHAsDkdw5h8juH1L1hHLGee47Jzz2XSqw8UV2iibNOWdU8tGMdWj4laqAEgCnvdGLKO+n0MZsybRpTpk1LJVaeqC7RxFmnrGoe2rEOLZ8SNVAiIhIkNVAiIhIkNVAiIhIkNVAiIhIk9YOq74r9oDZvrwCgacX2xPtBbd6ypRCrceNE4+SN6hJNnHXKquahHeus81E/KKlV04rt6cUK5D/K0Kgu0cRZp6xqHtqxDi2fEt3iEwAmrunMxDWd04k1dSoTp05NJVaeqC7RxFmnrGoe2rEOLZ8SNVACwDPrOvDMug7pxJo5k2dmzkwlVp6oLtHEWaesah7asQ4tnxI1UCIiEiQ1UCIiEiQ1UCIiEiQ1UCIiEiT1g6rviv2gdpFwPygRkarK9YPSFZSIiARJDZQA8NCqz/PQqs+nE+vJJ3noySdTiZUnqks0cdYpq5qHdqxDy6dEDZQA8Px7B/P8ewenE2vOHJ6fMyeVWHmiukQTZ52yqnloxzq0fErUQImISJDqbKDMbJSZtUojGRERkZIoV1AHA3PM7DEzG2xmNbz2JSIiEq86Gyh3/xFwGDAeGAG8bmY/M7N0nqhLKpo02E6TBumMaN6kcWOaBDp6cpZUl2jirFNWNQ/tWIeWT0nkflBm1gv4DjAY+B+gP/And/9hcuntTv2gYqZ+UCKSsb2eD8rM/hX4F2AdMA4Y7e6fmlkD4HUg1QZKRETqhygTFrYBznH3N6sudPcdZnZWMmlJ2sb+43AARnZamnysiRMLsYYNSzxWnqgu0cRZp6xqHtqxDi2fkigvSXSp3jiZ2QQAd1+cSFaSupffb8vL77dNJ9bChby8cGEqsfJEdYkmzjplVfPQjnVo+ZREaaC6V/1gZhVA32TSERERKSjbQJnZNWb2IdDTzDYWfz4E3gUmp5ahiIjUS2UbKHf/v+7eArjN3VsWf1q4e2t3vybFHEVEpB4q+5KEmXVz9yXA42bWp/p6d/9LoplJqg46YGt6sVq0SC1Wnqgu0cRZp6xqHtqxDi2fkrL9oMxsrLuPNLP/qWG1u/upyaZWM/WDipn6QYlIxva4H5S7jyz+7ylJJiYiIlKTKIPFnmdmLYq//8jMJpnZF5JPTdJ01xtHctcbR6YTa8IE7powIZVYeaK6RBNnnbKqeWjHOrR8SqJ01P2xuz9uZicAXwZ+AYwBjks0M0nVgg8/k16sJUtSi5Unqks0cdYpq5qHdqxDy6ckSj+o0giiXwF+7e6TgUbJpSQiIhKtgVptZr8Bvg5MNbPGEb9HcXqO18xsmZldXWabQWY238xeNbMZ0VMXEZH9WZSG5uvA08Bgd38f+Awwuq4vFUecuBc4AzgKGG5mR1Xb5iDgPuBsd+8OnLcnyYuIyP6rzmdQ7v6xmU0G2plZp+LiKDcs+wHL3H05gJn9HhgKLKqyzTeASe7+j2Ksd/ckeYlPu0ab04vVpk1qsfJEdYkmzjplVfPQjnVo+ZTUOR+UmV0G/AR4B9hRXOzu3rOO751L4arru8XP3wKOc/dRVba5E2hIYby/FsCv3P3hGvY1EhgJ0KlTp75vvvlm9U1kb6kflIhkbK/ngwIuB45w9/V7GrOGZdX/8h1AYeDZ04CmwEtmNsvdd5nzwd3HAmOh0FF3D/MQEZEcitJArQQ+2It9rwIOqfK5I7Cmhm3WuftHwEdm9jzQC0h+UiLZxc//fjQAP/z835KPNW5cIdZ3v5t4rDxRXaKJs05Z1Ty0Yx1aPiVRGqjlwHQz+y9gS2mhu/+yju/NAQ4zsy7AauB8Cs+cqpoM3GNmB1B4df044I6IuUuMXvvowPRivfFGarHyRHWJJs46ZVXz0I51aPmURGmg/lH8acQe9H9y921mNorCG4AVwAPu/qqZXVxcP8bdF5vZfwMLKTzfGufuyf8TXkREghflLb4bAMysefFWXGTuPhWYWm3ZmGqfbwNu25P9iojI/i/KWHzHm9kiYHHxcy8zuy/xzEREpF6LcovvTgpj8E0BcPcFZnZSkklJ+j7XdFN6sdq3Ty1Wnqgu0cRZp6xqHtqxDi2fkij9oF529+PM7K/u/oXisgXu3iuVDKvRfFAxUz8oEcnYvvSDWmlmAwA3s0bAv1K83SciIpKUKGPxXQxcCnSg0G+pN3BJgjlJBm58vRc3vp7ORfGN997Ljffem0qsPFFdoomzTlnVPLRjHVo+JVGuoI5w9wuqLjCzgcDMZFKSLLy5uTK9WGuq99cWUF2iirNOWdU8tGMdWj4lUa6g7o64TEREJDZlr6DM7HhgANDWzK6osqolhY63IiIiiantFl8joLK4TYsqyzcC5yaZlIiISNkGyt1nADPM7EF31/wW+7kjmu/NeMB7GatLl9Ri5YnqEk2cdcqq5qEd69DyKYnSD+pw4AdAZ6o0aO5+aqKZlaF+UDFTPygRydi+9IN6HBgDjAO2x52YiIhITaI0UNvc/deJZyKZ+vclfQD4Wbe/JB/rjsKMKj/7t39LPFaeqC7RxFmnrGoe2rEOLZ+SKA3Uf5rZJcAT7Dof1HuJZSWpe2dr0/RirVuXWqw8UV2iibNOWdU8tGMdWj4lURqofyn+7+gqyxzoGn86IiIiBVHmgwrz9Q4REdmv1dZR91R3n2Zm59S03t0nJZeWiIjUd7VdQZ0MTAOG1LDOATVQ+5FeLdJ7pNirW7fUYuWJ6hJNnHXKquahHevQ8impsx9UaNQPKmbqByUiGSvXDyrKYLEiIiKpUwMlAFy56FiuXHRsOrFuuYUrb7kllVh5orpEE2edsqp5aMc6tHxKorxmLvXA+9sapRfrww9Ti5Unqks0cdYpq5qHdqxDy6ekzisoM2tmZj82s/uLnw8zs7OST01EROqzKLf4fkthBInji59XATcllpGIiAjRGqjPu/vPgU8B3H0zUMOrXyIiIvGJ8gxqq5k1pdD3CTP7PFXG5JP9w3EHrU0vVs+eqcXKE9UlmjjrlFXNQzvWoeVTEmU+qC8B1wJHAc8AA4ER7j498exqoH5QMVM/KBHJ2F7PB+Xuz5jZPKA/hVt7l7t7mEPfiojIfqPOBsrMpgCPAlPc/aPkU5IsXPq3/gDce/Ss5GPdeGMh1nXXJR4rT1SXaOKsU1Y1D+1Yh5ZPSZSXJG4HTgQWmdnjZnaumTVJOC9J2Sc7KvhkR0U6sbZs4ZMteoxZneoSTZx1yqrmoR3r0PIpiXKLbwYww8wqgFOB7wEPAC0Tzk1EROqxSCNJFN/iGwIMA/oADyWZlIiISJRnUBOB44D/Bu4Fprv7jqQTExGR+i3KFdRvgW+4+/akk5HsnPSZt9OLdWw6g9LmjeoSTZx1yqrmoR3r0PIpKdsPKtQZddUPKmbqByUiGdubflCaUVdERDJTtoFy958Uf73R3d+ous7MukTZuZkNBn4FVADj3L3GCUfM7FhgFjDM3f8QZd8Sr4sWDgRgfM+Zyce69tpCrJtvTjxWnqgu0cRZp6xqHtqxDi2fkij9oP6jhmV1NiLF19LvBc6gMEzScDM7qsx2twJPR8hFRETqibJXUGbWDegOHFjtOVRLIEpH3X7AMndfXtzf74GhwKJq211GoREM8ymdiIhkorZnUEcAZwEHsetzqA8pdNatSwdgZZXPqyi8rr6TmXUAvkahA3DZBsrMRgIjATp16hQhtIiI5F1tz6AmA5PN7Hh3f2kv9l3TnFHVXw+7E7jK3beblZ9iyt3HAmOh8BbfXuQiIiI5E6Uf1F/N7FIKt/t23tpz9wvr+N4q4JAqnzsCa6ptcwzw+2Lj1AY408y2ufuTEfKSGH2pzer0Yg0cmFqsPFFdoomzTlnVPLRjHVo+JVHmg3ocWAJ8A7gRuABY7O6X1/G9A4ClwGnAamAOhQ6/r5bZ/kHgj3W9xad+UDFTPygRyVi5flBR3uI71N1/DHzk7g8BXwF61PUld98GjKLwdt5i4DF3f9XMLjazi/csfUna5u0VbN6ezmjmm7dsYXOAIydnTXWJJs46ZVXz0I51aPmURLnF92nxf983s6OBt4HOUXbu7lOBqdWWjSmz7Ygo+5RkjHq1MB9UGv2gRhXnngmtz0XWVJdo4qxTVjUP7ViHlk9JlAZqrJm1An4MTAEqgbBmtRIRkf1OlPmgxhV/nQF0TTYdERGRgto66l5R2xfd/ZfxpyMiIlJQ2xVUi9SyEBERqaa2jro3pJmIZOvsdv9IL9app6YWK09Ul2jirFNWNQ/tWIeWT0mUflCHA78G2rn70WbWEzjb3W9KI8Hq1A8qZuoHJSIZ25d+UPcD11B83dzdFwLnx5ueZG3Dp43Y8GmjdGJt3MiGjRtTiZUnqks0cdYpq5qHdqxDy6ckymvmzdx9drWx8rYllI9k5AeLC2P1ptEP6ge33lqIFVifi6ypLtHEWaesah7asQ4tn5IoV1DrzOzzFAd6NbNzgbcSzUpEROq9KFdQl1IYSbybma0G3qAwHp+IiEhionTUXQ580cyaU7ji2gwMA95MODcREanHyt7iM7OWZnaNmd1jZqcDHwP/AiwDvp5WgiIiUj/VdgU1AdgAvERhBt0fAo2Ar7r7/ORTkzR9/bNvpBfrjDNSi5Unqks0cdYpq5qHdqxDy6ekbD8oM3vF3XsUf68A1gGd3P3DFPPbjfpBxUz9oEQkY3vTD6o0zQbuvh14I+vGSZLz9pYmvL2lSd0bxhFr7VreXrs2lVh5orpEE2edsqp5aMc6tHxKarvF18vMSj23DGha/GyAu3vLxLOT1Fz7Wl8gnX5Q1955ZyFWYH0usqa6RBNnnbKqeWjHOrR8Smobiy+d6VVFRERqEKWjroiISOrUQImISJDUQImISJCiDHUk9cC3OyxLL9bQoanFyhPVJZo465RVzUM71qHlU1LnfFChUT+omKkflIhkbF/mg5J6YMXHlaz4uDKdWKtXs2L16lRi5YnqEk2cdcqq5qEd69DyKVEDJQD8dFkvfrqsVzqx7ruPn953Xyqx8kR1iSbOOmVV89COdWj5lKiBEhGRIKmBEhGRIKmBEhGRIKmBEhGRIKkflADwvUOWphfr65rvsiaqSzRx1imrmod2rEPLp0T9oOo79YMSkYypH5TUasmmlizZlM4MKkuWL2fJ8uWpxMoT1SWaOOuUVc1DO9ah5VOiBkoAuG15D25b3iOdWOPHc9v48anEyhPVJZo465RVzUM71qHlU6IGSkREgqQGSkREgqQGSkREgqQGSkREgpRoPygzGwz8CqgAxrn7LdXWXwBcVfy4Cfi+uy9IMiep2WWdF6UT6MXruaz31p2/AzDg+nRiB263uoBqU4PLvvWtIPeVh7jlhJZPSWL9oMysAlgKnA6sAuYAw919UZVtBgCL3X2DmZ0BXO/ux9W2X/WDilna/aDU76o81UbqqSz6QfUDlrn7cnffCvwe2GXaRnd/0d03FD/OAjommI/UYv7GVszf2Gq/i5Unqks085csYf6SJcHtKw9xywktn5IkG6gOwMoqn1cVl5VzEfBUTSvMbKSZzTWzuWvXro0xRSm5e8VR3L3iqP0uVp6oLtHcPWECd0+YENy+8hC3nNDyKUmygarhfgU13q8ws1MoNFBX1bTe3ce6+zHufkzbtm1jTFFEREKV5EsSq4BDqnzuCKypvpGZ9QTGAWe4+/oE8xERkRxJ8gpqDnCYmXUxs0bA+cCUqhuYWSdgEvAtd09vOG0REQleYldQ7r7NzEYBT1N4zfwBd3/VzC4urh8DXAe0Bu4zM4BtNb3JISIi9U+i/aDcfSowtdqyMVV+/y7w3SRzkGhGd31lv4yVJ6pLNKMvuijIfeUhbjmh5VOiCQsFgG6VG/fLWHmiukTTrWvXIPeVh7jlhJZPiYY6EgBmbWjLrA3pvCGZZqw8UV2imbVgAbMWxDPgTJz7ykPcckLLp0RXUALA/SsPB6B/q+T7maUZK09Ul2juf+wxAPr36hXUvvIQt5zQ8inRFZSIiARJDZSIiARJDZSIiARJDZSIiARJL0kIAD8+NL03eNKMlSeqSzQ/vuSSIPeVh7jlhJZPiRooAaBzs037Zaw8UV2i6dyhtkkRsttXHuKWE1o+JbrFJwDMWN+OGevb7Xex8kR1iWbG7NnMmD07uH3lIW45oeVToisoAeDh1YcCcHLrd/arWHmiukTz8OTJAJzcr19Q+8pD3HJCy6dEV1AiIhIkNVAiIhIkNVAiIhIkNVAiIhIkc/esc9gjxxxzjM+dOzfrNPYftxsAb29pAsDBjT+BKxM8J263XWNBsvFy5O2fNQWq1AVUmxq8vbYwmO7Bbfd95Pc495WHuOVknY+Zzatpslq9xSdAtT+K+1GsPFFdoonzj2hWf5BDaZh48XoADi59fh0YcH02udRADZQA8PTa9gB8ue2a/SpWnqgu0Tz9wgsAfPmEE4LaVx7i7ualG4Bq554aKAnNY291AdL545hmrDxRXaJ57KmngHj+uMe5rzzELSfUc08vSYiISJDUQImISJDUQImISJDUQImISJDUD6q+K/aD2vBpIwBaNdyaeD+oXWKB+voUbbilMVClLqDa1GDDxo0AtGrZMqh95SHubtL+778M9YOSWu3yR3E/ipUnqks0cf5Rz6qByLxhqibUc08NlAAw+Z1DABjabuV+FStPVJdoJj/3HABDTzstqH1FUuwYO/mV9YW4PVoH0e8o1HNPDZQAMOWdTkA6J2iasfJEdYlmyrRpQDyNSpz7iqTYMXbKwoGFuJtmBtFAhXru6SUJEREJkhooEREJkhooEREJkhooEREJkl6SEADu6T5rv4yVJ6pLNPdcd12Q+9qjuIEd69DyKVEDJQA0rdi+X8bKE9UlmqaNG8e7r+Kr37tI+M260I51aPmUqIESACau6QzAsPYr9qtYeaK6lFGtAZn413XQ5UyGnXnmPu964tSp8NxDu9c84QYqtGMdWj4laqAEgGfWdQDSOUHTjJUnqksZxb5DJc8sHAhLNsTSQD0zcyas65B6zUM71qHlU6KXJEQkfza+mXUGkoJEGygzG2xmr5nZMjO7uob1ZmZ3FdcvNLM+SeYjIiL5kVgDZWYVwL3AGcBRwHAzO6raZmcAhxV/RgK/TiofERHJlySvoPoBy9x9ubtvBX4PDK22zVDgYS+YBRxkZp9NMCcREcmJxOaDMrNzgcHu/t3i528Bx7n7qCrb/BG4xd1fKH5+DrjK3edW29dICldYAEcDf0sk6fi1AdZlnUQEyjN+eck1L3lCfnLNS54QTq6fc/e21Rcm+Raf1bCsemsYZRvcfSwwFsDM5tY0sVWI8pKr8oxfXnLNS56Qn1zzkieEn2uSt/hWAYdU+dwRWLMX24iISD2UZAM1BzjMzLqYWSPgfGBKtW2mAN8uvs3XH/jA3d9KMCcREcmJxG7xufs2MxsFPA1UAA+4+6tmdnFx/RhgKnAmsAz4GPhOhF2PTSjlJOQlV+UZv7zkmpc8IT+55iVPCDzXxF6SEBER2RcaSUJERIKkBkpERIIUTANlZg+Y2btmVmMfp9qGRaprSKUMcr2gmONCM3vRzHpVWbfCzF4xs/lmNrem76eY5yAz+6CYy3wzu67KutRqGiHP0VVy/JuZbTezzxTXpVnPQ8zsf8xssZm9amaX17BNEOdpxFwzP08j5hnKeRol18zPVTNrYmazzWxBMc8batgmiPO0Tu4exA9wEtAH+FuZ9WcCT1HoO9UfeLm4vAL4O9AVaAQsAI7KONcBQKvi72eUci1+XgG0CaSmg4A/1rA81ZrWlWe1bYcA0zKq52eBPsXfWwBLq9cllPM0Yq6Zn6cR8wzlPK0z1xDO1eK5V1n8vSHwMtA/xPO0rp9grqDc/XngvVo2KTcsUpQhlVLN1d1fdPcNxY+zKPTvSl2EmpaTak33MM/hwKNJ5VIbd3/L3f9S/P1DYDHQodpmQZynUXIN4TyNWNNygqtpNZmcq8Vzb1PxY8PiT/W34YI4T+sSTAMVQQdgZZXPq4rLyi0PxUUU/qVS4sAzZjbPCkM4Ze344q2Ap8yse3FZkDU1s2bAYOA/qizOpJ5m1hn4AoV/nVYV3HlaS65VZX6e1pFnUOdpXTXN+lw1swozmw+8C/zJ3YM/T2uSpwkLyw2LFGm4pCyY2SkU/sM/ocrige6+xsz+F/AnM1tSvILIwl8ojIG1yczOBJ6kMLJ8qDUdAsx096pXW6nX08wqKfzh+T/uvrH66hq+ktl5WkeupW0yP0/ryDOo8zRKTcn4XHX37UBvMzsIeMLMjnb3qs94gzpPy8nTFVS5YZGCHC7JzHoC44Ch7r6+tNzd1xT/913gCQqX1Jlw942lWwHuPhVoaGZtCLSmFEYj2eWWSdr1NLOGFP44PeLuk2rYJJjzNEKuQZyndeUZ0nkapaZFmZ+rxVjvA9MpXM1VFcx5WqusHn7V9AN0pvwD/a+w60O92cXlBwDLgS7886Fe94xz7URhdIwB1ZY3B1pU+f1FCiO+Z5Xnwfyzs3Y/4B/F+qZe09ryLK4/kMJzquZZ1bNYm4eBO2vZJojzNGKumZ+nEfMM4jyNkmsI5yrQFjio+HtT4M/AWSGep3X9BHOLz8wepfC2ThszWwX8hMLDPbyWYZG8zJBKGed6HdAauM/MALZ5YcTgdhQut6FwIvzO3f87wzzPBb5vZtuAzcD5XjhLU61phDwBvgY84+4fVflqqvUEBgLfAl4p3t8H+HcKf+hDO0+j5BrCeRolzyDO04i5Qvbn6meBh6wwaWwD4DF3/6NFGGYui7+ntdFQRyIiEqQ8PYMSEZF6RA2UiIgESQ2UiIgESQ2UiIgESQ2UiIgESQ2UiIgESQ2UiIgESQ2USA6Y2XeLcwl9J8K2Xc1svJn9IY3cRJKiBkokH/43cCpwXl0bemGqhIuST0kkWWqgRGJgZh3NbLKZvW5my83sHjNrvJf7ampmM4pD1ZS8TGHqhJerbNfDzP5Y7ed/ldlnIzN73syCGd5MpC5qoET2kRUGWJsEPOnuh1GYCqIp8PO93OWFwCQvTJlQUklh0M8DSwvc/RV3P6vaz7s17dALk889Bwzby5xEUqcGSmTfnQp84u6/hZ1z8fwb8O3i3EF76gJgcumDmTWgMADpt4GvVbuy2o2ZtTazMcAXzOyaKqueLO5bJBd0uS+y77oD86oucPeNZrYCOBSYH3VHZtYI6OruK6osPhVY6O4rzGxB8fOfyu3DC/M6XVzDqr8Bx0bNRSRruoIS2XdGzbOO1jQ7aV3aAO9XW3YB/5z87lH28iqoeGW31cxa7M33RdKmBkpk370KHFN1gZm1pDAH0JFm9q9mdpqZ/Xtx3R1m1sHMGpvZhGr72gw0qbKfpsBQ4BfFK7LbgKHF5XujMfDJXn5XJFVqoET23XNAMzP7NkDxGdHtwD3u/ihwOIVbbrcWt78XGElhWvDfV92Ru28AKsys1EidDTzl7p3cvbO7dwL+Exiyp0maWWtgrbt/uqffFcmCGiiRfVSc3fVrwLlm9jqwHtjh7jcX3/A7CNhB8Zmvuy8D2vPPaberewY4ofj7BcAT1dY/AXxzL1I9hcJMqiK5oBl1RWJmZgMoPCs6BzgReB14Gxjm7j8sbnM+8Fl3v6OG738BuMLdvxVzXpOAa9z9tTj3K5IUNVAiKTOzgcDlwEXu/mGZbS4EHqrWF2pfYjYCznf3h+PYn0ga1ECJiEiQ9AxKRESCpAZKRESCpAZKRESCpAZKRESCpAZKRESCpAZKRESCpAZKRESCpAZKRESC9P8BQ8Agz0yRLFUAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sos_hb = read('/home/dbodesheim/zih_home_drive/PhD/SMAIS/GIXD/HSOS_HB_from_GIXD.cif')\n",
    "del sos_hb[[atom.index for atom in sos_hb if atom.symbol=='H']]  # deleting H-atoms\n",
    "result = get_reflexes(sos_hb)\n",
    "\n",
    "fig, ax = plot_result_gixd(result, reflections_exp)\n",
    "fig.tight_layout()\n",
    "fig.savefig('SOS/SOS_HB.png', dpi=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "24849aaa",
   "metadata": {},
   "outputs": [],
   "source": [
    "# exporting results\n",
    "export_results_csv(result, filename='SOS/results_full.csv', full=True)\n",
    "export_results_csv(result, filename='SOS/results.csv', full=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "992c0cf7",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8.5"
  },
  "vscode": {
   "interpreter": {
    "hash": "f73dddfd4a60ef34f659ef713707f286d0fb6cc1ea54e0d810b2463eb16cf75e"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}