minimal_lensing_sim.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys, os\n",
    "sys.path.append('../')\n",
    "\n",
    "import logging\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import astropy\n",
    "%matplotlib inline\n",
    "import pickle\n",
    "\n",
    "from simulation.wrapper import LensingObservationWithSubhalos\n",
    "from simulation.units import M_s\n",
    "import paper_settings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "paper_settings.setup()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Draw and plot lensed image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "lo = LensingObservationWithSubhalos(f_sub=0.05,  # Population-level substructure fraction\n",
    "                                    beta=-1.9,  # Population-level subhalo mass function slope\n",
    "                                    m_max_sub_div_M_hst_calib=0.01,  # Maximum subhalo mass in units of host halo mass\n",
    "        )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, '$\\\\theta_y$\\\\,[as]')"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQkAAAEFCAYAAADwqhIfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO2d65IdxZWFV7bkO0iNAN/Dg8UbyOIJpnmBCRnPC1j2A4zBfgIPeP7PIOYF8Mj+P4FmHmAA+a8dYdomwtjYRjfw3ahzfnSd0jqLs9epltRdfaT1RRDUUWVlZWWVUvuWe7feO0IIoWJr7gGEEI43WSRCCJYsEiEESxaJEIIli0QIwZJFIoRgOTn3AKbQWut0rOfGY3Xnbm1tlefuchx3dY7vzcfuGmVvb+/A1+kz83U8N9r/VFz/U/78IH0obj7cN1GhffB1PFc6T65/N44TJ06s7PMgc+XGVY1D4et67+/13p/UNhuxSDAnTy4P+ZOf/OR4/Pe//33p3Kc//enx+K9//WvZZ/VBAMuT+KlPfarsg18EfwB677/97W9lf+6D+8Mf/jAe8zMDwIcffjge8/zoM3/sYx8r+/jzn/88HruPj8/xs+i9eQ50PniueOzalu91+/btpXb8bPxc6/pnuH9973w/fk9//OMfy3baB49Lv81HHnlkPOa5//jHP16OV5+F2/K70Ln6xCc+MR7rN/aXv/yFx/H2qvu2TQimOnnyZD916hSA5ckFlidAP1p+aTyh/BcOWH6Z+hf3gw8+GI/5L5b7F1DHUf0l1o+br9PFkP+y6sfIf2F4PnSR4D703oz78Bm+F1D/i6VzpR/xlD50HPyXzn3DfJ3+JeP36f4l1nfBuH98eMFzkgSjY3T/+PA5985EWlg6x+N/99133+y9n9cxxSYRQrDMpm601naGw2d77y/MNY4QgmcWSaK1dg77i8MVAOdaa2fnGEcIYT2zSBK996sArrbWtgHs9t5317QfdTU28gC1oQtYtkPwOTXaMWpgYt3d6Z98neqO3IezZLPuq8/ijHiVMVHtDs5Ay2NxVnn+re+i0p/VduF05KqdM8apDagy6Km9iedbn4XfJ9/7M5/5zFI7/iZ0Trl//a4q47AztOpz8nzzM+t35ewyaohdxdzejfMAbq460Vq7CODicHyUYwohELMaLgd1Y7u1dmHFuUu99/O99/NZJEKYj1kkidbaiwDe6r1fwr4kcca1772P4tqZM8tNWRR0bj1G3V2sfqjYyf1X4h3wUVGQqYKpVITmc86dq9exCMljZB84sCz+6sJbie96L1YjVHTl+XeuUyf+8hir5wJ8/Apfx3PgAqH0XOWidCqWc1/q+6zUU3W3urmqVLWDBJepGrSKuSSJlwHsDh6O7WGxCCEcQ+YyXO4CWBgrr8wxhhDCNBJMFUKwzO3dmERrbbQBTN2sAtTuJNU3+Tp1j7LOxnql9qFuPobtF6zTujBhPce6qeqRrMdWOr3+VptKNVdOV3d7WVwffC91K/Nz8pyqbcS5R3mOXeh15fZVXAg4z7d+O5U9S/vhbQK67aAar/52oePOxjRlW0YkiRCCJYtECMGyEerG1tbWKB47sZ53bAK1eO12FCqVm8y5wtQVW41DRW1up3040ZXbTt0a7VyPLorVqVyV+1VF4SoiEqh3Y7px6DV8P35/qg7wdfoN8Bjd/LpdmlPfGd9Lv29+Tp2D6tvU+XDbz11qgfG+5dUhhIAsEiGENWyEunH79m28//77AD66waZKvAEsi2fOEuySvbAINjW1mIptVcTb3aavU7GZxUnnOWARV0VjFrf5Xk5VcNF7lcgPADdv3tmuo++M54S9J06EVlWhigp1KpbbDMdjUnXAZfHi/t3GMDdGnh83/qlJc1S9qyJLmUgSIQRLFokQgiWLRAjBshE2ia2trUnRfa4N6/Euy7NL0uqS1bhoRob7d7sSVY9nXdjpz3ysejC7bd1zsp6tdo3KVQosuz1d0hlG3xn/5nury5kzobvEOG4HpIvarN61c5fre+f+NbnLn/70p/GYn9m5enUeeSx8b/dudYw8jwu7nxJJIoRgySIRQrBshLrRex9FN5cX0sGimds8peJvFfHm3F1uA5lzW7GY6JLOqMhbRQ7qc7L4rmJn5Y51tUH03KI2CrCsHjjVyW18cuNw1ci4rXP/uehRnkcer4r8rDa44jkKi/luU57b9MdzzONyaqDOr4vQXRBJIoRgySIRQrBkkQghWDbCJrG3tzfqfpqUg3Uq1beqxLjq3mL9Vt1A3KezO7ikp1XILOulei+nq7vkJq6WBJ+bWphW3Wk8p6yPa/8u6a6r/8Fzx/PjEg27ehSMPjPbaJxti+dRv6kqoZCi75rHzLuXdRwuubOzZTBsX9H+XMLisc3aFiGEh5osEiEEy0aoG1tbW2P0oIpLLmquclk6l5lL7MH9O9XGiXAsZmpdDLeTlNUst8t0aq0Kl++Rn9O5BnUOWP3gaE8Vhfne2n9VRlDbTU2k4hLcVK5SHXNVphFYftdOlHfl9FzdDT6nfVRzoH3wezlIwqXxmrUtQggPNVkkQgiWjVA3GGcpVzGLxWa3aYn7VM8BpztnUVMj71xqcieSVug4nCWbRVIeo0s6o/Poog+n3EvHzPOjXpAqvyOw/M5cpXZ+Nldx3CXQcXkyK3TunQfGbS5jeD5U/OfrnGrG8DcLTK84XhFJIoRgmatg8DaAneHnM733F+YYRwhhPXNJEs8BONN7vwwArbWLM40jhLCGuQoGcxXxs9ivMu7aj/qe6pWubkBVa0P1SmcnqJKxurJzU5PJqD7Iv7V/1lVVb+XfVZm8df1PdeHy+HUcPK/s3nX2G2dHchGuLkFxlTxW37sbRxX9OnWHr55zuPoc/C40arPaGewS5qoNSO0Xq5jVcNlaOwvg+lBlXM9dBBAJI4SZmdu7caH3/q1VJwZp4xIAbG1tra9qGkI4FGZbJFprF3rvLw3HO733K6btKDKpOMairKtZ4MR8R7XxR1UbFkOdG5XbTa3sDfhIzaq+iHMNunyMLsrPlTOsys6pC9Ql7+F787HOt4t+rVQAnQ+OClU1gsfoaojwOZ2rqRuw3PfonrPKZek2AKqbVr+lleNb2+IQaK3tAHixtfZma+3NOcYQQpjGXIbLKwCenuPeIYSDkWCqEIJlbsPlJE6cOIFHH30UgN8l53Qxxu3kUx2N9Ux2h6pOz+Nwrk1Xj5P1VrdLU3VT1q15vDoOlwiGdVXn8qvG6/rXdq6GSOUOdK5BTV7MNqCpYfDOnuVwtTqdDYjhb07fmZv/yiXv+nBJlCsiSYQQLFkkQgiWjVA39vb2RheYuplYXHKiq4uIZHFvqiqi0W98L5f8piqFp+1crQo9x6oCi7VnzpxZasdzd/369aVz/DwucQ3johTdjlN1ZzJ8b1YfXTIgdetVblTn9lUx3L2LClV33bfJ34HLr8njcjuDXQlHt8vZvd8FkSRCCJYsEiEESxaJEIJlI2wSW1tbo06ueiXro1ozo8rI49xK6gKtdpI6N6fTn13YNOutavNgnKuXQ6D1WarkrjrGKrOT/nZ1Qrmdq2+p/fP4eUzs5tU+1CbB+rmzAVX30vtVti3F1RBROwz378KyecxT3bJut6j2wW2ruiGRJEIIliwSIQTLRqgbvfdRTFKxkH87McupAyySquhaibwqurpoRha3WZTXdq7MPIvG+pws9qsqwlS7RfW3S4LC6BirqEJVFaYm9WURXcV8jrJ078LtAmVUdZpah8SpItVOUm3Lx6oi8v1cEmgXkcvftH7frh7IgkgSIQRLFokQgmVj1I2FmKTiKf/WqDkWf53F3kW8scjLIq6KrjwOrXzOsDipIi57Z1zJOCfmszjpEp24Gg78zKq+OBGax+/K61XqgI6LRWjdxOW8FgzPo74zvreqoFVldc3XybiENHpvnhMXneqontt5k3Qcbu4WRJIIIViySIQQLFkkQgiWjbFJLPQsdbOxDqe6L8OuHpcgVt11VVl4lzBGo+tu3rw5HrvkIy7JCuMiKZ17kXVmF7Xp3G6su+scsGuWz7ldmmqTqFx5+iwcHajvjK/juVKbgXM5s31ratJd7cNFOlYRrjpXVW0Nhe0aLmGRzhWPuXKHRpIIIViySIQQLBuhbgB3xC63wUbdU5XorSIXoyJ0Vf5OxUIXpVgljFHxtKp9ASyL8ioWViUAdRzO1cZu26kiunMpcrupLmYdoytd5/pndyn3r8/CfaoLlOfHuTKrDYD62yW1cSULXaRw9X07VVjnQDdFriKSRAjBkkUihGDJIhFCsGyETaK1Nupqrnak6mzc1rlKXdIP1iVZb3d6tlLptNoH7/TUMfKzOPdl9cyKS+7K9o+pOw/1HNthVEd24/rggw/GY34WtTe5HZ2s4zsbjXuH1Tn9/rjGh4bjV7t/gWUXrrNd8LvWUH1u61zOVW3UVfdbxWySRGttp7X22lz3DyFMY7ZFwlURDyEcHzZC3ei9jyKTisksWqo4VuVSVJGxqlsBLIuCrA44N5OKtU5VYJyI7tQDFnOdysLPpuPgc1X9DMCXRDx9+vR4zHOgIq1TuapoTBXX3a7YCueGdLlNXRQrw6oH4HeMVhGSmlDI3buK9tR35nZ6umjPcQzl1TPTWrsI4OJwPPNoQnh4ObaLRO/9EoBLAHDixIlpJZRCCPedY7tIMK21UTR0iT2UauOT84I4z4cT+V0ptUpEdVWoXcVuFY0ZFpNVdeLfaomvVB3nEXB5Mt3GIbbsL6rFr+q/8gAAy8/iEhE5NbAqAQD4yFhme3t7PGbPjPYxNRW/tnM5KPk9sYrhyi+6KNxyfGtbHBKttQsAzg//DyEcU2aTJHrvlwFcnuv+IYRpJOIyhGDZCJsEu0BdmT919VS6nktiq1S7O1W343Zq12A37dT6HIrTHfkc73B1LmG1BfBvV4qQx+9KBfKz6L3YVahuyWrH7I0bN5basV3G2aVcgpupCW7ZrezKHuqzTE1w60r5VXYHoN7hqn8PeFzah7NvLYgkEUKwZJEIIVg2Qt3gquKrzlVUuQlVZGRxzG2AcfkpnTpQJTBxLlV9Lh6H1qDguWHXpkbvnTlzZmU7oB6/qkA8py46lZ9TE8s88cQT47G6+N57772V16lY7OqLVPk6nRvVuSi5naoULtlL5c4F6vl27nm9d6VmHSTZkHvucUxrW4QQHmqySIQQLFkkQgiWjbBJ7O3tjfqp6nasf7pkpi7UmHV8p/dVOyUXY5wyDhcmXF2jv1X3ZdcmHz/++ONL7T73uc+Nx7xjE1jW+dnlp/YEbufclzxvTz755FI71pE1lLmyE6jdwblReQ74PWk79ywMj0PtQc5mwN+B2tT4edy2AFcbpLLHObem2qm4/+q6SBIhBEsWiRCCZSPUjdbaKIa6hDEqCnKkX7UzEFhWWVzEm0vewWKtK2tX5YEE/G5AFsO1bsipU6fG4y996Uvj8Re/+MWldjwf+pwsKnNZQlXvXPQhj4PP6Xh5jOyWBZafk1UdVUt4TjXZSxXh6qIlnTua0XljNVDFf3bvuuQ9leoBLH9LLrK0cj/rOf02j/Uu0BDCZpBFIoRg2Rh1YyHWqUjH1llVByqrsUvH7sQvd51TZ6ryd9qfi35zoiaL+Z/97GfHY/UqsGis5d1YxWDRXtvxOedx4KhKVSlYzfr85z+/dI6fhVUgVnMA4J133hmPVc1kdYOvc9+OWvarEgyuDJ+qZm7TGMPfrc4Vz6mrFs7j1294asKbivLLa639EIBV1Hrv31h7hxDCRuMkiV0Ar5rzWSBCeAhwi8SrvfefVCeTwTqEh4NykagWiNbaqd77+24BOQwWupOzGbiIOpfQhd1fqleyvju1VL0mwnVJRRi2O6gNgiPlNFqSXZ1f+MIXxmNO0gosu+TUFca/2Q6hrkf+rf9Q8Nzxe1K7Bs/d008/XfbBdpLf/e53S+34nO4kZVsA2yf0G+DvRd9n5VLUPlxSGOfS5vnm8TrbiNo4KluJK4WpEZfqPl7FJMNla+2/AbwM4FkAj7fW/q/3/m9Trg0hbDZTXaCv9d5/DOC53vtzAG4d4phCCMeIqS7QJ1pr/wLgf1prTwE4e2gjWsHe3t4ohqpIx+qAiuhVnkVXp0HVARbxXA0EVjFUdK2S3+jmKT6nm7NYTGT3IrDsRuRNXJpbkkV5VQE42Qvnk1Txl0V0F33IYqz2weKwJr957LHHxmOOzHz77beX2rH64SJX+ZyK4U71q2qg6Ly5UoQuQpev43Y6V/xNaFRotenPRetOibBUpkoS3wfwi0GK+DqA6we+UwhhI5kkSfTebwH40XD8g0GaCCE8BEw1XD4FgCttfQPAM4cwnhDCMWOqTeIl7AdXnQNwFcAPD21EK+i9j3qh6vusp6mOyb/Zdef0MtXnNBx4gXOFqS5a7TLVRLKudgffT3dVsv2CXaD6nL/61a/GY7WvuMSyjEtoUrkK9Vm4D7aFAMthyZXeruf0Odlu4tzWVXIaxe305P61Dx6XfpvcD78LdZ/zdTr3lT1Ov2G2RWkfahdbxUG8G9+l/99zle/W2oXW2k5r7fl77SuEcHhMXSSebq39OwC01r6PewzJXhQJ7r1fAXCztbZzL/2FEA6PqYbL77bWvtp7/0Vr7ZsALt7jfZ/BnX0hCzXmStWYd4Gq2Dm1NB5f58qxOVVEdxsyt27dCR1R1yOLoSxmqmjJYrmqFFNL9LELUXE1KKqdqqpucTun3rGr8Pr1ZWfYz3/+8/FYn5PVJb63S5Kj8Dy6qEcWtbU/ftfcn6qIPKf6PnkO3I5fPufyWGr/VX/6LK7mDLugNbp2HEN109baP8mNfjH8/5Xe+0/0/AHZlt+Pa4PW2sXW2huttTemZg0KIdx/nCTx7daa82CcA/Dju7zvTQBnXIPe+yUAlwDgxIkTWSVCmAm3SLy85to37uG+r+OONHEWwGuucWutFPVZPXBWXRd55zb68HXczqXlV2t+ZXnWcbA46VLqq0jKkhaLjxqBymqJK9HnJLeqkjVQq2qsigHAu+++Ox5r9Cg/J78/9ca4vJBMlUtS76VqbPWcOqd8b90cyHOs/VcV2F1+Shcp7CI/+Vl0rjSCdBVuF+iP1l59l/TeL7fWnl8YLAcDZgjhGDJb+rre+0tz3TuEMJ0DLRKtte9g38j4HoCrvff/PZRRhRCODQdaJHrvP1gct9a+01o7hyNaLBZuItWDWRdzZefYneZ0R7UFsC7pStW7TF18P9YxXRJSjYxj95Tq55X+rIlwOTJTE67++te/Ho/5mdW+MjUjmZsr3iGq9gq27XBiGX233KdLTst2HtX3uU+1p/A88rfj3Is6DpcIt6r5of1XNVv0fq4eiiaaYdxO2AUHlST+FXeiLa8BeAXA+YP0EULYLA5qk3gVwO6wKxSttdPItvEQHmjsItFa+yr2d3/eAPCG5rUcFosjzXXpqnmruFdtAlLRkkVSzfnH4hiL9Rr9xqKl9lFFKbrNNa7OhkZBXrt2beXxV77ylaV2rH6o65HrXbCY78rwqVhbVRVXcZd/a5Qfqwo8jzwmYFkdU/WRVSRup+I/fweqVlWRmiqe83yo+sjX6fvk385F6ZIqTa274Ta53Y8yf/+I/XDpGwD+s7X2+j1GWoYQNox1i0QD8NYQM/H93vszAG5loQjh4cGqG733V1pr/9FaewzAbmvtLexHSz57JKMLIczOWsNl7/3bg4HyOQD/jP0YicuHPTDm9u3bo37qkruqTYJDlJ2ryp1jfZH1dtVNWSdUvbLahaf2FdbVXR0IdZNxUthf/vKX47G6QHnH5Ze//OWlc5z8hcelLkrWYZ1NhfV/rf/hXHJsh+CQYbXzuPfOVPYJYHk+VN/n/vldu527bkeyXlfZMtRuwt+w9s/Pw+9C59fNzz3VAmUGA+UrU9qGEB4spiadCSE8pMy2d+MgbG1tjWKYikfOLVlFyjkRy7keq9yJwLKI7sq9OdWmKncPLIvbv/3tb5fO8XP/7Gc/G4915yyrGOoCfeqpp1b298477yy1Y1FWXbH8bCw2647T6hpgOZqU3aOq9jh1g98Nz4FLVOOSwrjyizx+95z6zfH7dVGVPGbnWuc5cOq09l/lcF0a+9oWIYSHmiwSIQTLRqgbwB1xTUVGl4yF27r07IzbRFMldwG8+FttMlLxVFUMhsVC7Z+9AJzQ5ac//elSO5eWn71GLGrrnPKzaKRjZUXXPli0VzWi8m7o3LCa6dLtT60E70pEuqhEvrc+p3ufVR8K96ntePysVqnqdNgRlyGEh5wsEiEESxaJEIJlI2wS7AJ1OpTaKyrXlbpKXek61u9YB1S9j/tUXZRtGexC1MhM535l1BbA9gq2SWhEJF83NQry9OnTS795Dtw8ugSxPD+6u7OKsnTvVs9VO201mpHH6L4BRufU1TJxpR8rW5d+E2xPcDU5uH9X9nBKhKUSSSKEYMkiEUKwbIS6sbe3N4qoKiY7V1WlmriNOC5PoRPbnLuLRdmp1audG0vVkt/85jcr+9c+WOznzWrAcnlAfma3qUgjALmtc3OyyK7nOMrSVYJ3eSErNULVDX6HboOXcyHyvVwleI1wrUoRKlPbTXXx6/ftolDH+65tEUJ4qMkiEUKwZJEIIVg2wiaxtbVld9gtUH2/0pGruqKLezGssznbgktEWoXuqs3AhQm7RDCVzqzt2J6ju/+4Lbs23a5EV5uS2+m92NagrsdqV6K2c7tpeb45/NzZDPRdVO5LZ7tQWwDbQHT8fD9nQ+B5dOHn7DrWxEzqqmam1FGZTZJore201myh4BDC/My2SKRIcAibwbFVN1prFwFcBOoaCHpOXVwskrKoqaLlVNHP1fjge7vSae5Zqv4AH3lXub800pHHxTktgWVVhMeloiqP3+3u5LnXKEV+F253p0vo4ly9lYvbvXeF36+rn8HjUtWMn8WVoJxan0O/uSr3pt7LfXPHWt1YR+/9Uu/9fO/9/N2EkoYQ7g+HJkkMkoCyGzUjhM3i0BaJ3vul+9nfQox2VZ3Vq1CJUipmOtG1SuzhktM4D4mL8qu8A9pWN11Vor2K3U69uXHjxso/Z6u53kvntyqbp1XQ3WY49nywSK5qD78nt1HOJYVx46j6OEi+Tjff3LYqR6ljdlXRK/VIf+sYXQTwgjm9GxcAnB/+H0I4psxmuOy9X8YRF/kJIRycWARDCJZj6wJlWmujXqjJUVg3dbsBXZISV1eBdULnkmOdVvtg3Zr1c3W3sq3B2RP0Op6DqtaIjsNFUnJ/+px8ziUA4h2izhWr7lG2E3D/bj6cDu52UU7dkevqofBvZ6/Qb4LHyHOgc8rXufKRVa0RYDna07n4KyJJhBAsWSRCCJaNUDc+/PBDXLt2DcBHRSmXH5BFaBYLVczkPl1SET6nYjLfS2tacFun2jimRmqyuOoSiqjYyXPCIr+K+fwsOgcs1rr55nur+lhFyU5V4Va1nTIOjXCtcpHqRi3+dnRzGs+/qiL8m9u5aElXN4TH6NRp5Vi7QEMIm0EWiRCCJYtECMGyETYJRsOmnWuT9UznxnI7LFkH5eSxqtM7nZB1VdYrVR9nt6ELeVbYhsDjVR2Wn82566pEO4Dflcg4u4mza/A743tpO7632zHLz+LmUHXz6p05u4B+E9xW56qaH7fL2SV65ns7l7DC306V8CeSRAjBkkUihGDZCHWjtTaKxweJGFMRdYGLWFSxsHIHHiQZC6sVzh01NW/G1PoLLkqxmhttp6I8X+dUIlZn3Dh0Pqp5dK66qd+Ey4WpIvnUSEqX95Tvp2obzyuPQ1VhVkGn7hrW+XA7SV0NmgWRJEIIliwSIQRLFokQgmVjbBILXcq53Zxu6txRzm3IOpzuZmRcbQPWVfn4IG5U7n+qXqlj4j71WTi82CWIZbuDS8zKz6auO243NYuXcy/q8/P9nI3G2Qy4D/duq5qhgN+5W+G+P3VR8jkeo7OrVW5OO6YDXxFCeKjIIhFCsGyEugHcEdPVzeSSzlTl/FTUrhKdaFunDqwa66oxsjtNXWsu8pNdj6oOVPVA1M3JIq/ODV/n3GnOjVrV61A3Ko9fxfzKvavPzK5CFa+rnaSuTJ6qFHyOx6iuUhbfnZivrk3GlRF0yXUqdUnfrRuji0JdEEkihGDJIhFCsGyEurG3tzeKtppghMUlFceqqEJVB1x0IJ9jkdEluFGxliMuWdzjaDqgrvWgv1UkZVGTj50Irf1X1cJdEhRXVtElA6ra6W9+t65ivCuN594Zo+eqzXCq9jzyyCPjsaqxfG9NRFQlnXGeGp3vKn+nU220f65zUhFJIoRgySIRQrBkkQghWDbCJtFaG3Uz1U2djly5KTXZC7uIVOdkvZ71eDcO7b9KoKv3cglGqkhE7Yevc3VNVW/l56kSv+j4db4rd53ahpyOz+NwkbAu4pJtAy6ZrksuzPPI34c+i4sKZdTVy7/5vbjdrvpN8PzzuPRZXN3Uqn4rM8si0VrbBrAz/Hym9/7CHOMIIaxnLnXjOQBnhnqgaK1dnGkcIYQ1zCJJ9N4v0c+zAF7WNsPCcXE4HkUyV3bOladnkVTdqNynRqTxb2538+bNpXaVWrLq9wIVH1nUdDkdVXTl52RVR+eKRc1bt24tnXv00UdX3tttinI1SvhYn3PqJj3nsmRcAh12Pbr+dEzcp3Mpuo1VjH4vPC5WG6bWjtF7V1G3es5FdFbMapNorZ0FcL33vqvnhoXkEgCcPHlyfQWREMKhcGiLRKFC7Pber9DvC733bx3WGEII986hLRKiUnyE1tqF3vtLw/GOLB4hhGPCXN6NHQAvtta+N/yR9W5sbW2N4a+6w411R9WRWSfn61zdDe2jKgv/2GOPfWSMC9Q2wi45tjuozYBDZNUVtr29PR5r4hBuy+e0f35ut1OQn1nduYyGlXNbHpOrlcJhzYraXhh+ZzrGauenzmnlQgSW7VY8j/p98LvVOXUJhSs3s9uh7MK+uZ3ON9s/3N+RKqnSXIbLKwCenuPeIYSDkYjLEIJlIyIu9/b2RvFPRXkWJ1VcuhuRzu2Kc/kpWYxTMZldUiwmq9rjksIwKl6zaMx9ul2mGo3KorsX6nIAAAgxSURBVLETT1mUdbsveYxud65zo/I5HW8VZQrUOT9dPQq325XbuaQtLuennqt2froEOndbX8QlIppS6yWSRAjBkkUihGDZCHWDN3i5nHyqRlS5FFUcZRFdrdBV1WgV6bgP5xGYWtlbx8HeB/UI8L35WJ+T1QFVFfh+Lm8oq2Yul6fbTOa8Cgxf58bh1DuX+9Gpd1WyGvUc8G9NLOPK9/F3zM/mVAWXb5RVS7cZbmopyaXrD3xFCOGhIotECMGSRSKEYNkIm0TvfdTpVDdlnU31bNb5nau0iqoE6oQdqu87dx235T5Ul3Yl6Xgc6tqsEvmq65H7V3uCi25kKncuUO+IdK5BZ2Pi96L2m6m6tXP/TXVp871dHwrPz9QSlPoNu5KLjHPFMmrr+v3vf7+270gSIQRLFokQgmUj1I3W2igmqSjP4p5TAapSe9qHSzDCbjLnKnWuQSdq83UqTp8+fXo8VtWgqkGh4+DITB1/NQ63ecqJxixeO9edws/CY9RndpuuqqhNnQ8nllf5UfWaqi4LsDwf7ntxuTZ5HOqm5efRd8GwiuRqzlQVxyNJhBAsWSRCCJYsEiEEy0bYJPb29kZ9ybnMnG7KutdUN6f+dnUgGNUPp+qfPEbnktTrqj7VfXY3uwj1Xjwu1bPZfsF6sN7XuR75940bN8p2rKvrzuCqnQtFV7d4Va/U7YpV+Dq1ZVS7TLU/nn+1k7CNgvtw9VVdrZeKSBIhBEsWiRCCZSPUDd4F6sQjV5PDuTa5nYqdVdSmwmK4U2ecq4rFSacqOBery2NZ7RbVcfE5jTCs8lgCyyKvS67DqDpTuSX1vbskQrxLdmqUqdvV61TcqS5KpapR4tyyOo4qAljnykXJuu9x7G9tixDCQ00WiRCCZSPUjd77KFKrSsGipVp/OYLM5RFkUV7FLxbxuJ1usnIRhix2skiqHoxqUxHgrdyMs9hzn+qZ4Plhb4GKpy5Sk+ek6g/wYrirns44Fc5FhTIuoQvD78xFi7rEMtp/VQVcn5nnW9/F1JKIPFd3s7EvkkQIwZJFIoRgySIRQrBshE0CqCPbWKdy0W8uscfUqDzWCdUdxdepfsh2k6o0m17nks6oHsljnBrNeJBIyqoPte1UEYBTd30C9c5dd52zMVV9azunq7t35ubbueB5/DxvOkbnSq6+OX2WU6dOjcf67UyJwp1tkRjqgQLAs713Wws0hDAfs6gbrbVz2F8crgA411o7O8c4Qgjrmatg8FUAV1tr2wB2e++72qa1dhHAxeF4FMNc/kite+BcbYxzgVZuOCdaqpjMpQOrcnrrxuvcdVWUnivDp2J+FfXncj+q+FsleHFu5ffff3/l2PU65xJ2uTxdO7dJjH+7jXfOxepEeR6Lq/viNuwxU3OF6jfmrlswt03iPICbq0703i8BuAQAJ06cmPa3PYRw3zm0RWKQBJTdQcUAAPTer7TWvt5au9B7v3xYYwkh3D2HtkgMksBKWmsvAnhraHMTwJnDGkcI4d5oU/X2+3rTfUPlwlj59d77t9a0HwfJCWGB6bvf+JyGKzvXoAtRrtqp3srh4dzOhfFOdUMCdYIUHS/r1lN3mbrkrrrrsdqZqd8Y7yRle42Og1H7jdsxy7h3xvOhti63e5Th+VHbkPt2qiQx+iz8TTh72VT3tr6zW7du8Zje7L2fhzCX4XIXwMJYecW1DSHMSyIuQwiWWdSNg9Ja+z2At4vTTwB47wiHU5FxLJNxLLMJ4/iH3vuT+ocbsUg4WmtvrNKjMo6MI+O4P+OIuhFCsGSRCCFYHoRFoozHOGIyjmUyjmU2dhwbb5MIIRwuD4IkEUI4RLJIhBAsWSTCfaW1dqG1ttNae37mcey01l6beQzbw3xcGPYrzTmWneG/A4/jgVok7mUi7uMYjtuHcWR/UVprF4D93b0AblL2sSOHdxvPyHMAzix2OBc7ow+de03y9MAsEsco29Wx+DCAWf6iPIM7e3J2AZw74vsfK3rvl2g39FnMtE+p93619/6CS/LkmDvpzH1jSrarIxoHu5jOAnh5jnHMxLb8fnyWURwzhn+wrs/1TRJlkifHAyNJEHc1EfebY/RhHCXJDbKaC+vSIRwFg2S5vVALp7JRksRxyXY1ZRw4gg9j4jiOktdxR5o4C2BWw+FxYPgOXxqOd+Z4N/ea5OmBCabiiZBJmWMsF8gmMcuHQWN5rff+7BHe73kAVwGcnWv+h3FcAPAKgG/OlRpxMNy+jDuS7QszLRIHSvL0kesfoEXinibiPo7jWHwYw1hm/4sSNp8HZpEIIRwOD6LhMoRwH8kiEUKwZJEIIViySIQQLFkkwn1j2K/y2kGDdVpr51prL7bW/uuwxhbuniwSYRKttYvDf+v2xFw9qLt1CKl/9e5HFw6TLBJhLa21lwH8EMAbAA4kJYTNZ6PCssPRM2yYO9N7vzlIEZP3xSwiXwE8i/2gst2hv+8BuIb9XbtHFg0a7o5IEmEdzwHYHewM38O+NLGWYUHZHkKzd3FHAtkBcG3Yz/Aw7ZDdWLJIhHU8DeDVhZ2h9351SKyzsFGszJcx7H59a9jLcY7+/DKAZ1prN7AvYYRjThaJsI5r2JckdgB8f/izxeatr+FOkpklhsXj8UFiuEp//jz2VY/HAJwZ1I9wjMkiEdZxGfsqwjZ7LYZMYG+iWCSwr5Yscl1eA/AN8ox8a1hEXu+9z577I3hiuAyWQW3QhYANmCoJnKOt8l+jP39J/j8yLDjfuA/DDYdAdoGGECxRN0IIliwSIQRLFokQgiWLRAjBkkUihGDJIhFCsPw/OQi83MHuWaQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(np.log10(lo.image_poiss_psf),\n",
    "        origin=\"lower\",\n",
    "        extent=(-3.2,3.2,-3.2,3.2),\n",
    "        cmap='gray'\n",
    ")\n",
    "\n",
    "plt.xlabel(r\"$\\theta_x$\\,[as]\")\n",
    "plt.ylabel(r\"$\\theta_y$\\,[as]\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.027962647364270506"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Substructure fraction of realization. For reasonable population parameters almost always << 1\n",
    "lo.f_sub_realiz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Latents"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "20048789651023.85"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Host halo mass (in units of solar mass)\n",
    "lo.M_200_hst / M_s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.6255315455943502"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Host halo redshift\n",
    "lo.z_l"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-0.23218257468061865, -0.07666151302764661)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Offset between host halo and source (in arcsec)\n",
    "lo.theta_x_0, lo.theta_y_0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "81"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Total number of subhalos simulated\n",
    "lo.n_sub_roi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([5.16274425e+07, 1.70961667e+07, 1.47169020e+07, 1.32654494e+07,\n",
       "       1.01505565e+07, 1.99775060e+07, 1.44590496e+08, 1.12686441e+07,\n",
       "       1.84481325e+07, 2.88534100e+08, 1.63713408e+07, 1.10539586e+07,\n",
       "       1.06309525e+08, 5.68231924e+08, 1.50569323e+07, 1.06074530e+07,\n",
       "       1.76348788e+07, 4.22143549e+07, 2.58630169e+09, 1.05588652e+08,\n",
       "       1.58864009e+07, 1.80187348e+07, 1.44815492e+08, 1.90040192e+07,\n",
       "       2.17841648e+07, 1.93683632e+07, 1.96155867e+07, 1.01991934e+07,\n",
       "       6.02456731e+07, 1.05548385e+07, 1.83485584e+07, 1.28610950e+08,\n",
       "       1.10583210e+08, 1.04055763e+07, 2.69715783e+07, 1.51373651e+07,\n",
       "       5.50971284e+07, 1.34842165e+07, 1.22148626e+07, 4.34393270e+08,\n",
       "       1.28746275e+07, 1.26912156e+08, 1.21321757e+07, 2.17654997e+07,\n",
       "       1.40268035e+07, 1.00804887e+07, 2.64077868e+07, 1.27132325e+07,\n",
       "       1.24327538e+07, 9.27964847e+07, 1.74704434e+07, 2.79866833e+08,\n",
       "       9.10903937e+07, 1.45551612e+07, 1.42514770e+07, 7.56636584e+07,\n",
       "       1.27017939e+07, 3.26998647e+07, 1.63635408e+07, 1.56301437e+07,\n",
       "       1.77707667e+07, 2.50292428e+07, 1.83448218e+07, 1.60805940e+07,\n",
       "       8.08769491e+07, 4.84785003e+07, 4.14402514e+07, 1.09627757e+07,\n",
       "       1.45175818e+07, 1.16092437e+07, 4.63751005e+07, 1.07999827e+07,\n",
       "       2.11062113e+07, 4.02583359e+08, 1.96089723e+07, 6.09934059e+07,\n",
       "       4.35853326e+07, 3.11685513e+07, 1.17996674e+07, 5.99585840e+08,\n",
       "       3.89668775e+07])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Individual subhalo masses (in units of solar mass)\n",
    "lo.m_subs / M_s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.65108709, -0.8136336 ],\n",
       "       [-0.1246937 , -0.63666673],\n",
       "       [ 0.1869741 ,  1.2226059 ],\n",
       "       [-0.55446315,  0.50298466],\n",
       "       [ 0.56460766, -0.72808265],\n",
       "       [-0.62779132,  0.92975894],\n",
       "       [-0.74174773, -0.39803346],\n",
       "       [-0.49867318,  0.81475312],\n",
       "       [-0.42250352, -1.27951671],\n",
       "       [-0.87960838, -0.74501478],\n",
       "       [-0.18345298, -0.11003733],\n",
       "       [ 0.43654922, -0.62092081],\n",
       "       [-0.58329054, -1.03721378],\n",
       "       [-0.25452263, -0.93019287],\n",
       "       [ 0.85971938,  0.43448988],\n",
       "       [ 0.66514154, -0.36773676],\n",
       "       [-0.88471986, -0.465222  ],\n",
       "       [ 1.2216813 ,  0.39097068],\n",
       "       [ 0.16134747,  1.24755692],\n",
       "       [-1.02873357,  0.66746609],\n",
       "       [-0.59145621, -1.1554476 ],\n",
       "       [ 0.04900974, -0.32003078],\n",
       "       [ 0.28470942,  0.30346029],\n",
       "       [-0.08199624,  0.25187952],\n",
       "       [-0.60143204, -0.18724202],\n",
       "       [ 0.0124062 , -0.62317415],\n",
       "       [-0.97699172,  0.03846383],\n",
       "       [ 0.71162975, -0.37813217],\n",
       "       [-0.22208534, -1.20596586],\n",
       "       [-0.59997878,  0.05548499],\n",
       "       [-1.2703098 , -0.43769786],\n",
       "       [-0.10471441, -0.24586029],\n",
       "       [ 0.86197177,  0.13278475],\n",
       "       [ 0.75403823, -0.20520892],\n",
       "       [-0.31730685, -0.04397015],\n",
       "       [ 1.12815189,  0.6758003 ],\n",
       "       [-0.10285284, -0.01169327],\n",
       "       [-0.05781633, -1.31873144],\n",
       "       [ 1.0961721 ,  0.76128255],\n",
       "       [-0.01422615, -0.0183545 ],\n",
       "       [ 0.46872653, -0.80605826],\n",
       "       [ 0.38284647, -0.55405764],\n",
       "       [-0.74280656, -0.918974  ],\n",
       "       [-0.36445292, -0.44096497],\n",
       "       [ 0.35743113, -1.10630785],\n",
       "       [ 1.08095292, -0.11424464],\n",
       "       [-0.25870556,  0.69029464],\n",
       "       [-0.43534574, -1.21617118],\n",
       "       [ 1.17444362,  0.01383384],\n",
       "       [ 0.05503819,  1.15787217],\n",
       "       [ 1.14479538, -0.2806164 ],\n",
       "       [ 0.80807317,  0.38205517],\n",
       "       [-1.00206651, -0.4952198 ],\n",
       "       [-1.07985756, -0.1440877 ],\n",
       "       [ 0.52621389, -0.45892087],\n",
       "       [-0.40599934, -0.85450753],\n",
       "       [-0.61710558,  0.99227811],\n",
       "       [-0.66295273, -0.76408826],\n",
       "       [-1.05147983, -0.34431318],\n",
       "       [-0.48497144, -0.76587398],\n",
       "       [ 0.17577579, -0.18063283],\n",
       "       [-1.08565308, -0.36384077],\n",
       "       [ 0.76342337, -0.69431405],\n",
       "       [-0.97390699,  0.33116153],\n",
       "       [ 0.54651889, -0.28635678],\n",
       "       [ 1.03239008,  0.03255068],\n",
       "       [-0.54907908,  0.56882513],\n",
       "       [-0.3480698 ,  1.22989026],\n",
       "       [ 0.00151224, -0.74991352],\n",
       "       [ 1.24226201, -0.37773206],\n",
       "       [ 1.14896363, -0.7075067 ],\n",
       "       [ 0.36653309,  0.03237355],\n",
       "       [-1.17538065,  0.26499354],\n",
       "       [ 0.8678504 , -0.37022779],\n",
       "       [ 0.99109003,  0.53224796],\n",
       "       [ 0.57716332, -1.1381609 ],\n",
       "       [-1.25220538, -0.29448906],\n",
       "       [-0.53319796,  0.05737951],\n",
       "       [ 0.94287823,  0.38856831],\n",
       "       [ 0.87484896, -0.84614542],\n",
       "       [ 0.24031837,  0.98246115]])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Positions of individual subhalos (in arcsec)\n",
    "np.transpose([lo.theta_xs, lo.theta_ys])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}