{ "cells": [ { "cell_type": "markdown", "id": "dba61810-609c-4b0c-82cd-cf2f822f9b37", "metadata": {}, "source": [ "# Image-based analysis (experimental)" ] }, { "cell_type": "markdown", "id": "8a5adb47-4130-47ae-b741-951ff94ab4e4", "metadata": {}, "source": [ "Although FRAP primarily operates on the visibility plane, there are situations where visibility analysis is not suitable. For instance, if a disk exhibits significant azimuthal asymmetry, azimuthal averaging is not adequate. In such cases, FRAP can also process radial intensity profiles. While this function has not been extensively tested, we demonstrate its usage in this notebook." ] }, { "cell_type": "markdown", "id": "39501131-e135-4e7a-b4b1-33ecc9310706", "metadata": {}, "source": [ "The basic usage is the same as the [normal procedure](tutorial_1.ipynb)." ] }, { "cell_type": "code", "execution_count": 2, "id": "e69039c3-1759-40e0-a114-195055e34b85", "metadata": {}, "outputs": [], "source": [ "import frap\n", "import numpyro\n", "import pickle\n", "import dsharp_opac\n", "import jax.numpy as jnp" ] }, { "cell_type": "code", "execution_count": 3, "id": "f62cba31-87b0-4962-881e-15725e5dfe58", "metadata": {}, "outputs": [], "source": [ "numpyro.set_host_device_count(10)" ] }, { "cell_type": "markdown", "id": "8ece8314-ccc3-4e61-ae81-3e743a759b4d", "metadata": {}, "source": [ "In this example, we use a mock data set, and assume a cut along the major axis. " ] }, { "cell_type": "code", "execution_count": 17, "id": "1db1197f-2b98-47df-a1d5-ef6958166791", "metadata": {}, "outputs": [], "source": [ "bands = [ 'band3', 'band6' ]\n", "data = {\n", " 'band3' : 'data_band3.pkl',\n", " 'band6' : 'data_band6.pkl'\n", "}" ] }, { "cell_type": "code", "execution_count": 18, "id": "181c6c3f-9b39-4827-bda0-bfe4097fa3d3", "metadata": {}, "outputs": [], "source": [ "obs = {}\n", "\n", "for band in bands:\n", " with open(data[band], 'rb') as f:\n", " obs[band] = pickle.load(f)\n", " " ] }, { "cell_type": "markdown", "id": "09b745dc-40cd-4350-817a-692d3180e3a0", "metadata": {}, "source": [ "he structure of the data should be like this." ] }, { "cell_type": "code", "execution_count": 19, "id": "99c1de53-93dd-43ec-ab3b-40801b98b7f3", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "dict_keys(['r', 'Tb', 'Cov', 'nu'])" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "obs['band6'].keys()" ] }, { "cell_type": "markdown", "id": "037d68e4-f164-4659-b4ac-9090ac0de5d8", "metadata": {}, "source": [ "The unit of Tb should be brightness temprature. Cov is the covariance matrix of the noise." ] }, { "cell_type": "code", "execution_count": 20, "id": "108a0cf5-7ca2-4414-9b05-c7b51dfb07f7", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAGdCAYAAACyzRGfAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjksIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvJkbTWQAAAAlwSFlzAAAPYQAAD2EBqD+naQAAR6pJREFUeJzt3Ql8lOW5//8r+0YWEghJIEEWWWRVVETU4lIRWyvVVtQelx5ba6v+fsrpsdJ6TrUb1vbf2gWx7fEntsfdiru4YIFWFgVFQPY9CAlr9j0z/9d1zzxjAgEyk2eW55nP+/UaZpJMkpkkzPN9rvu67zvB6/V6BQAAIEISI/WNAAAACB8AACDiqHwAAICIInwAAICIInwAAICIInwAAICIInwAAICIInwAAICISpYY4/F4ZO/evZKdnS0JCQnRfjgAAKAbdM3S2tpaKSkpkcTERGeFDw0epaWl0X4YAAAgBOXl5TJgwABnhQ+teFgPPicnJ9oPBwAAdENNTY0pHljHcUeFD2uoRYMH4QMAAGfpTssEDacAACCiCB8AACCiCB8AACCiCB8AACCiCB8AACCiCB8AACCiCB8AACCiCB8AACCiCB8AACCiCB8AACCiCB8AACB2w8fcuXNl7NixgX1XJk2aJG+++Wbg401NTXL77bdLQUGB9OrVS66++mqprKwMx+MGAADxED50i9wHH3xQVq1aJStXrpSLLrpIrrzySvn000/Nx++++2559dVX5fnnn5fFixfL3r175aqrrgrXY0eIGlra5PH3d8iDb26UT/dW83MEAERUgtfr9fbkC+Tn58uvfvUr+drXviZ9+/aVp556ytxWGzdulJEjR8qyZcvknHPO6faWvLm5uVJdXc2utmEy+40N8qcl283t0vwMWThziqQmMwIHAAhdMMfvkI847e3t8swzz0h9fb0ZftFqSGtrq1xyySWB+4wYMULKyspM+Die5uZm84A7XhA+1Y2t8uSK3YG3yw83yrMry/mRAwAiJujwsXbtWtPPkZaWJrfddpvMnz9fTjvtNKmoqJDU1FTJy8vrdP9+/fqZjx3P7NmzTVKyLqWlpaE9E3TL0x/slrrmNhneL1se+Moo877/+aevCgIAQEyGj+HDh8vq1atlxYoV8t3vflduuukmWb9+fcgPYNasWaZEY13KyzkLD6dFm/ab6387p0y+NmGAJCcmyK5DDVJ+uCGs3xcAAEuyBEmrG0OHDjW3J0yYIB9++KH87ne/kxkzZkhLS4tUVVV1qn7obJeioqLjfj2toOgF4dfa7pHV5VXm9jmDCyQrLVnGl+bJyl1HZOm2gzIjv4xfAwAg7HrcZejxeEzfhgaRlJQUWbhwYeBjmzZtkt27d5ueEETf+r010tTqkdyMFBnSt5d537lD+5jrpdsORfnRAQDiRXKwQyTTpk0zTaS1tbVmZsuiRYvkrbfeMv0at9xyi8ycOdPMgNFO1zvvvNMEj+7OdEF4aYVDTRjYWxITE8ztc4cUyO8XbjHhQyc+JST43g8AQEyEj/3798uNN94o+/btM2FDFxzT4PHFL37RfPy3v/2tJCYmmsXFtBoydepUeeSRR8L12BGkjzqED4sOuyQlJsiB2mbZX9ss/XLS+bkCAGInfDz22GMn/Hh6errMmTPHXBB7NlXWmusx/XMD70tPSZJTCjJl24F62VRRS/gAAIQdK0vFibZ2j+w6VG9uD+6b1eljw4uyzfVmfzgBACCcCB9x4rOqRmlt90pacqKU5GZ0+tiwfr7woZUPAADCjfARJ7Yf8FU9BvXJCjSbWnTBMUXlAwAQCYSPOLH9YNdDLmpYYNilTjyeHm31AwDASRE+4sT2A3WBysfRBuZnSmpSojS2tpvhGQAAwonwEWfDLoP7+BYX6yg5KVEG5Pv6QFhmHQAQboSPOGHNdDmli8qHGtA701yXH2GPFwBAeBE+4kC7xyuVtc3mdv+8zjNdLKW9rcoHwy4AgPAifMSBg3XNJoDoSqZ9s7vexK8031f52EPlAwAQZoSPOLDX30TaLzvNBJCulAaGXah8AADCi/ARByqqm8x1Ue7x920ZEBh2oecDABBehI84sM8fPoqPWtm0q2EX3VyuqbU9Yo8NABB/CB9xYF+1byil+ASVj96ZKZKVmmRus9YHACCcCB9xVPk40bBLQkJCYLrtHvo+AABhRPiIo56PEw27qH7+cLK/xnd/AADCgfARTz0fecevfFizYay+DwAAwoXw4XK6UVylv5JRlHOS8OH/uHV/AADCgfDhclWNrdLm36m2T6+uFxiz9MvxfZzwAQAIJ8KHyx2q8w2h5GakSGryiX/dhYHKB8MuAIDwIXy43MG6FnNd0Cv1pPe1hl1oOAUAhBPhIw72denOkEvHYRdtONVeEQAAwoHwESfDLn26UfnQgJKQIKZH5HCDr2ICAIDdCB8ud6jeP+ySdfLKR0pSYuB+NJ0CAMKF8OFywQy7dBp6oekUABAmhA+XC6bhVBX6Fxqj8gEACBfCh8sF0/Phu19ap+EaAADsRviIk8pHd4dd8v0h5ZD/8wAAsBvhI04qHwXdDB99/A2nh+tZaAwAEB6EDxdrbGmX+pb2oHo+8rP8lQ+GXQAAYUL4cLFD/upFalKiZKcld+tzrJDCsAsAIFwIHy5W1dBqrntnpUiCrh7WDdY6H1ZwAQDAboQPFzviX6U0L6N7Qy4dKx+H61vE62WJdQCA/QgfcVD5yMtM6fbnWD0fre1eqWlqC9tjAwDEL8KHi1VZlY8gwkd6SpL08veHWDNlAACwE+EjHno+Mrs/7NKx+qFDLwAA2I3w4WJH/OEjN4jKR8e+D2uBMgAA7ET4cLGqxpaQKh8FVD4AAGFE+IiHhtOMICsf1nRbej4AAGFA+IiLhtPgKh+9/ZUPa9gGAAA7ET5cLJSpth3vb4UXAADsRPhwsarG0Ga79LbCh//zAQCwE+HDpTweb6ByYYWJ7sr1r4hK5QMAEA6ED5eqbW4Tj3919GCn2gYqH/R8AADCgPDhUlbVIjM1SdKSk4L6XKtBlWEXAEDUw8fs2bPlrLPOkuzsbCksLJTp06fLpk2bOt1nypQpZgfVjpfbbrvN7seNME2zPbrhVIdvAACIWvhYvHix3H777bJ8+XJ55513pLW1VS699FKpr6/vdL9vf/vbsm/fvsDloYcesvVB4+Sq/c2iOSGEj1z/52ju0OEbAADs5NtBrJsWLFjQ6e158+aZCsiqVavkggsuCLw/MzNTioqK7HuUCFpNk39p9RDCh24ul5GSJI2t7VLd0BrS1wAAICw9H9XV1eY6Pz+/0/uffPJJ6dOnj4wePVpmzZolDQ0Nx/0azc3NUlNT0+mCnqtpbAu58tFx6OUIa30AAKJZ+ejI4/HIXXfdJZMnTzYhw3L99dfLwIEDpaSkRNasWSM/+MEPTF/Iiy++eNw+kgceeCDUh4GTDbukhxo+UmVfdRNNpwCA2Akf2vuxbt06+de//tXp/bfeemvg9pgxY6S4uFguvvhi2bZtmwwZMuSYr6OVkZkzZwbe1spHaWlpqA8LRw275GSE9iu2GlVZ6wMAYLeQjkx33HGHvPbaa7JkyRIZMGDACe87ceJEc71169Yuw0daWpq5wF41Pax89M5irQ8AQAyED6/XK3feeafMnz9fFi1aJIMGDTrp56xevdpcawUEkVPT1LOej89XOWWJdQBAFMOHDrU89dRT8vLLL5u1PioqKsz7c3NzJSMjwwyt6Mcvv/xyKSgoMD0fd999t5kJM3bsWJsfOrpX+Qhx2CWwvwubywEA7BXUkWnu3LmBhcQ6evzxx+Xmm2+W1NRUeffdd+Xhhx82a39o78bVV18t9913n72PGkH0fIQ47MIS6wCAWBl2ORENG7oQGZzf82Gt7WF9HQAA7MLeLq7v+Qht2MUKLVYFBQAAuxA+XKqnlQ9ruMZarAwAALsQPlyoqbVdmts8Per5sEKLtVgZAAB2IXy4UK1/yCUhQSQ7LcRhF/9wDcMuAAC7ET5cyAoMvdKSJTExoUcNpw0t7dLa7quiAABgB8KHC/W038MKLkdXUgAAsAPhw4V6urqpSk5KDAQQptsCAOxE+HChnq5uarE+n74PAICdCB8uZA2TZPdg2KVj5YQZLwAAOxE+XKi2ya7KB2t9AADsR/hwobpmX+WjV0/Dh7XQGKucAgBsRPhw8bBLxxkrPVrrg4XGAAA2Iny4kG09H+zvAgAIA8KHC9U1t9o77ML+LgAAGxE+XNzzEerS6harYZXZLgAAOxE+XKguMOzSs/BhLbFOwykAwE6EDxeyr+HUGnZhZ1sAgH0IHy5Ua9NUW2vYxhrGAQDADoQPFw+79GRjuY7hxfp6AADYgfDhMm3tHmlsbbdl2MWaqmtVUgAAsAPhw2U6DpH0dNjFCi/6NT0eb48fGwAAivDh0mbT9JRESUnq2a/Xmi3j9Yo0+KspAAD0FOHDrfu6pPWs30OlJSdKcmKC7+vS9wEAsAnhw7VLq/dsyEUlJCR83nTqXzUVAICeIny4jBUS7AgfHfs+rFADAEBPET5cxq4FxrpqOgUAwA6ED9f2fNgTPqy1Quj5AADYhfDh2p6PnjecKqvng7U+AAB2IXy4jFWh6JWWZMvXo+cDAGA3wodbh13sajhliXUAgM0IHy7T0OILH5mp9oSPzzeXY6otAMAehA+XqW+2Z18XC7NdAAB2I3y4dNgly67wYTWcss4HAMAmhA+XDrtkpdrTcGrNmmGdDwCAXQgfLlPnH3axrfJh9XxQ+QAA2ITw4TL1Ng+7WMu0M+wCALAL4cOtwy42r/PBsAsAwC6ED7c2nKba3XDKVFsAgD0IHy7S1u6RplaPrVNtP1/no028Xq8tXxMAEN8IHy7S0OprNlWZdg27+CsfHq9IY4evDwBAqAgfLmw2TUlKkLRke8JHRkqSJCUmmNvMeAEA2IHw4SJ2z3RRCQkJn28u5//6AAD0BOHDjWt82NRsamFnWwBA1MLH7Nmz5ayzzpLs7GwpLCyU6dOny6ZNmzrdp6mpSW6//XYpKCiQXr16ydVXXy2VlZW2Pmh0rSFQ+bBnyOXotT4YdgEARDx8LF682ASL5cuXyzvvvCOtra1y6aWXSn19feA+d999t7z66qvy/PPPm/vv3btXrrrqKlseLCK7r8uxa30w3RYA0HNBHaUWLFjQ6e158+aZCsiqVavkggsukOrqannsscfkqaeekosuusjc5/HHH5eRI0eawHLOOefY8JBxPPWBfV1sDh+scgoAiJWeDw0bKj8/31xrCNFqyCWXXBK4z4gRI6SsrEyWLVvW5ddobm6WmpqaTheEpj6wr4u9wy6scgoAiInw4fF45K677pLJkyfL6NGjzfsqKiokNTVV8vLyOt23X79+5mPH6yPJzc0NXEpLS0N9SHEvHLNdOu1sy+ZyAIBohg/t/Vi3bp0888wzPXoAs2bNMhUU61JeXt6jrxfPAuHD5mGXQMMpU20BADYI6Sh1xx13yGuvvSZLliyRAQMGBN5fVFQkLS0tUlVV1an6obNd9GNdSUtLMxf0XH1Le1gbTmuofAAAIl350L09NHjMnz9f3nvvPRk0aFCnj0+YMEFSUlJk4cKFgffpVNzdu3fLpEmT7Hi86Ebloxc9HwCAGJYc7FCLzmR5+eWXzVofVh+H9mpkZGSY61tuuUVmzpxpmlBzcnLkzjvvNMGDmS7hZw2LZIZptksdO9sCAGwQ1FFq7ty55nrKlCmd3q/TaW+++WZz+7e//a0kJiaaxcV0JsvUqVPlkUceseOx4iQa/MMudu1o29XOtgAA9FRQR6nubKmenp4uc+bMMRe4ZJEx1vkAANiIvV1c2PORafvy6v6ptlQ+AAA2IHy4SLiGXVhkDABgJ8KHG4ddwrXOR1Nbt4beAAA4EcKHK1c4Dc/y6m0erzS1emz92gCA+EP4cAmPxxsYdrG74TQzNUkSEny36fsAAPQU4cMlGlp9wSMcPR8JCQmBoRyrugIAQKgIHy5hhYLEBJG0ZPt/rdZQDpUPAEBPET5cuKOtVirsZg3lED4AAD1F+HCJ+ubwTLO1WF+XYRcAQE8RPlwiXKubWljrAwBgF8KHSzS0WGt82DvN1mKFGqvCAgBAqAgfLhGpygfDLgCAniJ8uIRVkQhX+GC2CwDALoQPl4jUsAuzXQAAPUX4cImwD7uwyBgAwCaED5ewejHCNdWWygcAwC6ED5eo9+/rkmnzjraWXv6dbWk4BQD0FOHDJcK1o+2xs12YagsA6BnCh0sw7AIAcArCh0tYFYnMsK3z4auo1Ptn1QAAECrCh0tYocAKCWFrOG0ifAAAeobw4baptmFqOLW+Lut8AAB6ivDhEg1hXuHUajhtbvNIW7snLN8DABAfCB+um+0S3nU+fN+LGS8AgNARPlzA6/UGej7CNdU2NTnRXFQdTacAgB4gfLhAU6tHPF7f7aww9XwodrYFANiB8OECVhNoQoKucBqeyodiZ1sAgB0IH27q90hNlgRNIGGiX18x3RYA0BOEDxew+j3CWfVQDLsAAOxA+HABa/ZJuHa0tbCzLQDADoQPFwj3NFsLO9sCAOxA+HCBiA27+Hs+6ltY5wMAEDrChwuEe0dbC8MuAAA7ED5coC7MS6tbrE3rmO0CAOgJwocLNAR6PsI77GKFG6vSAgBAKAgfLmAtd54VxtVNO4YPdrYFAPQE4cMFIjbbxap8sLcLAKAHCB8u0BDo+YjMImNWjwkAAKEgfLhAXYQqH/R8AADsQPhwAWsYJNxTbVleHQBgB8KHi5ZXzwx7wylTbQEAPUf4cFXDaYQ2lmtpE6/XG9bvBQBwL8KHC0R6hVOPV6SxlaZTAECEwseSJUvkiiuukJKSEklISJCXXnqp08dvvvlm8/6Ol8suuyzEh4fusPZaCfewi+4dk5Dgu81aHwCAiIWP+vp6GTdunMyZM+e499GwsW/fvsDl6aefDvkB4sR0+CNSlQ8NkoHN5ZhuCwAIUdBHq2nTppnLiaSlpUlRUVGojwlBaG7zSJuOg0Sg58P3PZKltrmNJdYBALHV87Fo0SIpLCyU4cOHy3e/+105dOhQOL4NdIGxDtvbh3vYpWPAqW1ifxcAQGhsP1rpkMtVV10lgwYNkm3btskPf/hDUylZtmyZJCUde2be3NxsLpaamhq7H5KrWUMuGSlJkpTob8gII9b6AADEXPi49tprA7fHjBkjY8eOlSFDhphqyMUXX3zM/WfPni0PPPCA3Q8jbkRqddNjVjllfxcAQKxOtR08eLD06dNHtm7d2uXHZ82aJdXV1YFLeXl5uB+SqzRYO9pGoN/D933Y2RYA0DNhP13es2eP6fkoLi4+bnOqXhAaa5O3rAj0eyiGXQAAPRX0Eauurq5TFWPHjh2yevVqyc/PNxcdQrn66qvNbBft+bjnnntk6NChMnXq1B4/WBwrUtNsLexsCwDoqaCPWCtXrpQLL7ww8PbMmTPN9U033SRz586VNWvWyBNPPCFVVVVmIbJLL71UfvrTn1LdCHP4yIzwsIv1fQEACHv4mDJlygn39XjrrbeCfhCwY1+XSFU+2FwOANAz7O3ikqXVrZVHI9ZwymwXAECICB8Ox7ALAMBpCB8OF62GU3o+AAChIny4ZapthBcZs74vAADBIny4ZZGx1MjMdqHyAQDoKcKHw0V6efXP1/lgqi0AIDSED4eL9FRbaxl3wgcAIFSED4dr8E+1zYrw8uotbR5pbfdE5HsCANyF8OGaYZfIrnCqmPECAAgF4cPhIj3VNiUpUVKTfX82DL0AAEJB+HDJCqeZEQofnWe8MN0WABA8woeDac+F9l5Ecnl1872Y8QIA6AHCh4N17LmIVM+H73sx3RYAEDrCh4PVNvnCR3pKoiQnRe5Xae1sS8MpACAUhA8Hq2+JbLOphcoHAKAnCB8OVtcU3fBB5QMAEArCh4NFeml1i9XcSvgAAISC8OFg0Qof7GwLAOgJwoeDWZWH7EhXPtKt2S6tEf2+AAB3IHw4WJ1/ka+ID7sEZruwyBgAIHiEDzc0nPorEZHCbBcAQE8QPhwsWlNtP19e/fNFzgAA6C7ChwsWGYv4VFtmuwAAeoDw4WD1UZ/tQuUDABA8woeDWQd/qwE0UtjVFgDQE4QPV4SPlChNtaXyAQAIHuHDBbNdIrmjbcfvpw2vXq83ot8bAOB8hA8XzHbJjvBUW2vYRXNHQwtrfQAAgkP4cEXlI7LhIyMlSRITfLeZbgsACBbhwxU9H5ENHwkJCYHptvR9AACCRfhwqNZ2jzS3eaISPjpWW1hiHQAQLMKHQ3Uc7oj0sEvHplMqHwCAYBE+HMo66KclJ0pKUuR/jb3SfdN7CR8AgGARPhwqWv0ex+5sy1ofAIDgED4cyjroR3pHW0sWDacAgBARPhy+qVyWPwREGjvbAgBCRfhwKGuWSdQqH4HZLgy7AACCQ/hwqLrm1qj2fHy+sy0rnAIAgkP4cCjroB/thlMrBAEA0F2ED4eK1tLqx/Z8UPkAAASH8OHwTeWsCkT0hl3o+QAABIfw4fDZLr3SfIt9RRqzXQAAoSJ8OJQ1y8Ra5jzSqHwAAEJF+HAoa7gjO9pTbf3DPwAAhC18LFmyRK644gopKSkxW6u/9NJLnT7u9Xrlv//7v6W4uFgyMjLkkksukS1btgT7bdDN8EHDKQDA9eGjvr5exo0bJ3PmzOny4w899JD8/ve/l0cffVRWrFghWVlZMnXqVGlqarLj8cL6PUR5b5fArrb+3hMAALor6CPXtGnTzKUrWvV4+OGH5b777pMrr7zSvO+vf/2r9OvXz1RIrr322mC/HWJ0Y7lsf6NrS7tHWto8kprMCB4AoHtsPWLs2LFDKioqzFCLJTc3VyZOnCjLli3r8nOam5ulpqam0wXBNJxGt/LR8bGg5zweLz9GAK5na/jQ4KG00tGRvm197GizZ882AcW6lJaW2vmQ4mCqbXTCR3JSoqT5qx2s9dFz+jO89a8rZfT9b8lv3t4kbe0eG74qAMSmqNfKZ82aJdXV1YFLeXl5tB9SzGtt90hzmyeq4aPj92bGS8/ocOW/z/tQ3l5fKQ0t7fL797bK3EXbbPkdAYDrw0dRUZG5rqys7PR+fdv62NHS0tIkJyen0wUn1nGYI1rDLh2/N8MuPfPexv3ywY7DkpGSJLdeMNi87/+9v0MamMYMwKVsDR+DBg0yIWPhwoWB92kPh856mTRpkp3fKq5Zwxza5BnNRk8rfFhDQAieadJ+1zcV/cZJA+WeqcNlYEGmHGlolec+pAoIwJ2CPnLV1dXJ6tWrzcVqMtXbu3fvNut+3HXXXfKzn/1MXnnlFVm7dq3ceOONZk2Q6dOnh+Pxx/cCY1GsenTcV4bN5USaWttldXmVNLcFt9Hep3trZO1n1SZEatVDe2lunHSK+dhbn3auIAKAWwR99Fq5cqVceOGFgbdnzpxprm+66SaZN2+e3HPPPWYtkFtvvVWqqqrkvPPOkwULFkh6erq9jzyORXumi4X9XXwaW9rlur8sN+FDV5z94/VnyBeG9e3Wz/DVNXvN9cUjCqWgV5q5feHwvvLT10RW7jpshl4yU6P7ewaAqFc+pkyZYkrFR180eCitfvzkJz8xs1t0YbF3331Xhg0bZvsDj2fWMEe0wwf7u/jc99I6Ezys383MZ1fLwbrmbk2rfe2Tfeb2V8aVBN4/qE+W9M/LkNZ2r6zYfjhMvz0AiOPZLnDevi4WKh8iB2qb5aXVn5mfx1///WwZUZQth+pb5Jdvbjzpz+/j8iPyWVWjZKUmyYUjCgPv1wB/wbA+5vaSLQfC+BsEgOggfDhQTaMvfOSkx0jlI45nZbz6yV5p93hlXGmeXDCsr/z8q2PM+zWQ7K858ZYCr/qrHpeOKpL0lM67E587xBc+Vu06ErbHDgDRQvhwoNqmVnOdk+5b4jxamGor8uLHe8zP4qrT+5vrCQN7y5kDe5shk3lLdx73Z6eLiL225tghF8vYAbnmeuO+WrOuCwC4CeHDwT0f0R92ie/N5bSvY91nvu0Avjy2OPD+b53vW6vjmQ/LzSyYrqzYcdh8fl5mikwe6qtydFSWn2l+v7p3zubK2rA9BwCIBsKHgysf2VGufPTyby5X1xzc9FK3WLnTNyQyvF92YKaKumRkoRTnpsvh+hZ5c52vutHVcI2aNrq4y7VatO9jdImv+vGpP+AAgFsQPhyoxl9pyMmIcuXDX3mxwlC8WbnTNxNlwim9O71f1+q4/uwyc/tvy3Yd83m6C/Cb63x7HV0x7vOKydHG+IdedB0QAHATwocDxUrlw2p4jdcVTlf6m0HPOip8qBlnl0pyYoJ8tLtKPt3bOTws3nxAqhtbpTA7TSYOKjju1x9V4ttqgPABwG0IHw6ufES758MKP7XNrXG5sNg6f0XizIH5x3y8MDtdpo727Wf0v8t3d/rY35b7qiFXji+RpMSE436P04p94WPr/jqzlg4AuAXhw4FqGmNjtotV+bCm/saTTZW10ubxSp9eqTKgd0aX97nhnIHmev7He6SiuikQJJZsPiAJCfpx3zLqxzOwIMuEE13XpeIk03YBwEkIHw5UG2OVDz04xtuZ+aYKXxPoiKIc0xzalYmD8uWMsjxpavXILxdsND+jB9/cYD52ych+UlaQecLvoY2ousmcFVoAwC0IHw4UMz0f/oZXXWSroSW+ZrxsrPBNfx1elH3c+2gouf8ro0yVY/7Hn8nXHl0m727YL6lJiTLzi93bcmBo317mmvABwE0IHw6j+4FYy6tHe4XTjJSkQM9CvDWd6uJfJwsfauyAPPnPqcM7rVY66/IRMtLfz3EyQwsJHwDch+0yHaa+pU08/hGOnIzoVj70zF6HfqoaWk01pig3PnYu1uET7flQI4tOHiK+N2WojCrJlWXbDpkdaycOPv4Ml6MRPgC4EeHDYawKQ0pSgqR1sThVpGnTq4aPmjha6+NAXbNZQEyLPqf281UmTuYLw/qaS7Cs8LHtAD0fANwj+kcvhNhsmnLcRsdIspperem/8WBLZV1gNsrRG8LZbbC/5+NgXUvcLuYGwH0IHw5jVRii3e9xdPiIp56P7QfrzfWQvllh/1690pKlICvV3N51qCHs3w8AIoHw4TCxMtPlmIXG4uisfMcBX/gY1Cf84UNZU3J3HyZ8AHAHwofDxMoaHxZrobN4Wmhsx0HfsMugPt3r9+ipgfm+8EHlA4BbED4cJlaWVj922CWOKh8HI1358H2f3Yd93xcAnI7w4dCl1XOjPM02XjeX0x1py480mtuDI9DzoU7xD7tQ+QDgFoQPh9HdUGMpfMRbz0f5kQazomtmapLZlTYSrCXWCR8A3ILw4TBVDS3mOi/TNwMi2qwl1uNlqm3HZtNITXUuy/dVWPZVN0pzW3wtYw/AnQgfDq18RHt103itfOzyzzixqhGRoDvnaqVFV7bd4x/yAQAnI3w4TKwNu8TbbJdyf/iwqhGRoBWWMv+Ml92s9QHABQgfDlPtP8jHTPgIDLu0xlX4KM3PiOj3/bzvgxkvAJyP8OEwsTbbxXocVkXG7ayFvqxKRKToUu4dh30AwMkIH05tOI2x8NHQ0i6t7R5x+262OtslGuGDYRcAbkL4cBA9uNe3tMdU5aPjMu9ur37obrZNrR6zm21JXpSGXah8AHABwocDh1xiabZLUmJCYJVTt4cPq9+jODdDUpIi+19noL/BVYd9PDrtBQAcjPDhINbBPTst2Rz0Y0W89H2UH26MSrOpKslLl+TEBLPCamVtU8S/PwDYifDhIFUxtsZHvIWPaDWbquSkRBnQ2xd6WOkUgNMRPhwk1tb4sFiPp+OwkBtFM3yoUmutD/o+ADgc4cNBrIN7XmZshg+3Vz4+X+MjOuHDCj3W4wAApyJ8OEisVz6qGwgf4RSYbkv4AOBwhA8HsQ7uMRs+XFz50A3d9tU0RXXYhfABwC0IHw5sOI218JETB+Fjb1WTeL0iGSlJUpAVnR2FreEehl0AOB3hw0GO+Fc37R2lg188Vz46NpvqRm/RUOZfaOxgXYvUN8fHRn4A3Inw4SBH6n3hIz+T8BEvG8odvYNwb3+zsbXMOwA4EeHDQQ77ez6ofMTfTBcLe7wAcAPChxMrH1mx1fMRD+t8RHuNDwtrfQBwA8KHA8NH7xgddrEaYt3IGuYo7R0jlQ+m2wJwMMKHQ+ieHrX+JsP8GGs4tcJQQ0u7mZLqRrsPNXRq+owWwgcANyB8OESVf6aL7ienjYexRHe1tTa6q3LhQmO6vkpNky/4UfkAgBgMH/fff7+ZitjxMmLECLu/Tdw5bE2zzUyVxBja0Vbp48nzD71Y04HdOOTSp1eaZKQmxUTPx57DjeLxeKP6WAAgVMkSBqNGjZJ3333382+SHJZvE1eO1MfmTBeLPq5D9S1y2N+X4s5m0+hNs7UU56ZLcmKCtLR7pLK2SYpzo/+YACBYYUkFGjaKiorC8aXjllVRiLU1PizW+hNuHHaJlWm2KjkpUfr3zpBdhxpMHwrhA4AThaXnY8uWLVJSUiKDBw+Wb3zjG7J79+7j3re5uVlqamo6XXAsq6LQO8am2R7ddOruykf0w4ei6RSA09kePiZOnCjz5s2TBQsWyNy5c2XHjh1y/vnnS21tbZf3nz17tuTm5gYupaWldj8kl63xEauVj9ROjbFuDB/Rbja1WBUYrX4AgBPZHj6mTZsmX//612Xs2LEydepUeeONN6Sqqkqee+65Lu8/a9Ysqa6uDlzKy8vtfkiuajjNi9VhF38oOuzvTXET6yA/MMrTbC2DCrLM9Y5D9dF+KAAQkrB3gubl5cmwYcNk69atXX48LS3NXODMfV2O7flocd36Knv8s10G9fEd9KPtFP/j2HmQ8AHAmcK+zkddXZ1s27ZNiouLw/2tXE1nksT0sItV+XBZ+NBptjqjNTM1Sfpmx0ZIHtQnMxA+vF6m2wJwHtvDx/e//31ZvHix7Ny5U5YuXSpf/epXJSkpSa677jq7v1VcOVDbbK5j5QB4vJ6PIy6b7WJVFwYWZJk1a2Kl50OXeqlvaZcDdb6/CwCI62GXPXv2mKBx6NAh6du3r5x33nmyfPlycxuhO1gX2+HD2uzOGh5yix3+8GFVG2JBWnKSlORlyJ4jjbLjQL0UZqdH+yEBQHTDxzPPPGP3l4x7be2ewLCLrrIZi6xGWLetcLrT39R5ir/JM1Zo/4mGD318EwcXRPvhAEBQ2NvFAbSPQof2tdQeqz0fViNsbVObtLZ7xC12Hmzo1OQZK6zm1x3+xwcATkL4cFC/R0GvtMAGbrEmJyPFhCO3VT8+H3aJzfCx7UBdtB8KAASN8OGg8BGrQy5KQ5FVlTlY647w0dTaLnurG2Ny2OXUwmxzvW0/4QOA8xA+HOBgXUtMN5tarHBkNce6YU8XHe7qlZYsfXrF1nDXqf16mWvt+Whua4/2wwGAoBA+nDTNNoYrH24MH9aQyyl9MmNmmq2lMDtNstOSzRokVl8KADgF4cNJwy7ZsXX2fbQCf3XALeEjVme6KA1DQ/3Vjy37u943CQBiFeHDSWt8OKTyccg/TOR01kySWGs2tZxa6A8flfR9AHAWwocDxPrqpkeHD7esummtbhqLlY+OTadbaToF4DCEDwfYX9vkkPBhDbu0uGvYJYZWN+2q6XRjRU20HwoABIXwEeN047B91b7wUZybIbGsjz8cHfRXapysoaUt8HOP1crHaSU55nr7wXrzeAHAKQgfMa6msU0aWnxTKYtzY3sPD6snxQ0Np9ZQhlZzdHG3WKR7umg1TKcDb6yg6RSAcxA+Ypy1yJUu4JWekiSOaDitbxGPzgF1MOtgPqyfr68iVo3yVz8+3cvQCwDnIHzEuAp/6b8oJ7arHspa4bTd45WqxlZxss0OCR+nFfvCx/q91dF+KADQbYQPh1Q+SvJiP3ykJidKbkaKK4ZeNlX6wseIotgOH6NKcs01lQ8ATkL4iHH7qpzRbNpx5U21v8bZ4WOzP3wMi/HwMaa/L3xs2Fdj9qIBACcgfDik8lHsgMqHKvI3xe7zP24nqmpokUp/eIr1YZfS/AzTdNra7pU1exh6AeAMhA/HVD6cET6sx2n1qjiRNYShB3bdVC6W6TLrZw7sbW5/uPNwSFO5F23aL79csFH+smS72UwPAMIttl9ZEaggOGXYxWqMrahxbviwKghjB+SJE5x5Sr68ua5CVu06EtTnHa5vkbufXS2LNx8IvO8372yWB68eI1eO7x+GRwoAPlQ+YphOV93rryCUOCV8+B+nkysfaz+rMtdj/f0Usc6qfKzcedjMNOqOmqZWueGxFSZ4pCUnyjVnDpDTy/KksbVd7np2tSzpEEgAwG6EjxhWWdskLW0eSU5McMRsl47DLtbqoE6ufIwZ4IzwoSudZqcnS01Tm6wuP9KtoZa7n1lthpd0EbWX75gsD31tnLxw27nytQkDzKJlGkCcPmMJQOwifMSwnf5dVUvzMyU5yRm/qn7+YZdKhw676FDEniO+oa7RDql8pCQlypThheb2wg37T3r/J5bulIUb95up0fO+ebaMKPKtFZKUmCA/mz7aTC/Wn8Pv3t0S9scOID4544gWp3b5NzYbWBCbG5udqPKhq5w6cernJ3t8Qy6D+mRJTrpvzRInuGRk98KHTsn9xZsbze0fThtxTMDSVXR/fMUoc/upD3bLtgO+ZeYBwE6Ejxi285Cv8jEw3znhIy8zxfQQOHWtjxXbD3fqo3CKKcMKTeVCF0c7XmBobGmX//P0x2Yo76IRhXLTuad0eb9JQwrk4hGFpn/kz4u3h/mRA4hHhA9HVD5ic1fV4039tNb6cOKMlxU7DpnriYMLxElyM1PkwuF9ze2/LdvV5X1+9vp62bK/zuzB89DXxprf1fF8d8oQcz3/48/kgAt2KQYQWwgfDqh8nNLHOZWPjtNt91Y5a6Gx+ua2QLPpxEH54jQ3TvJVMl5YtUfqmts6fWzBugp5csVuc/u3M8YFNgE8ngkDe5vZLy3tHnnK/3kAYBfCR4zSGQlOrHxYDbLKaQtW6ToZOtTQPy8j8Byc5LyhfWRw3ywTPH61wNfXYS0V/4O/rzG3v3PBYDn/VF+F5ES0KnLjpIHm9gsflTt+l2IAsYXwEaP21zZLQ0u7JCaIDOjtjDU+LFaPyi6HhQ9rsa1zHDbkYklMTJD7/c2iTyzbJY8u3iZ/X7VHrvvzcqlubJXxpXnyH5cO7/bXu2xUsVnhtfxwo3wQwuqpAHA8hI8YtdG/pbvOukhLThInKfPPztntoPChlaZ3N1Sa2188zTdzxIkuGNZXvjnZN/zy4Jsb5T+e/8TMPBpVkiPzvnmWmV7bXRmpSfLlscXm9vMr94TtMQOIP4SPGLVxn29/kRHFvjUYnKTMgcMuW/fXya5DDebg3J1hiVj2318+TWZfNcas1zGsXy+5+5JhZgGxvMzUoL/W188cYK7fWLvvmD4SAAgVe7vEeOVjZIxv6X6i8KGzXXStD107Ita99WmFuZ48pECyYnwzue70a1x3dpm59NQZZb1lcJ8s2X6w3gSQa84steUxAohvVD5ilC4GpazVJ50kPyvV9AroMt3WaqGxPuTy948+M7enjfENM+DzIHP1hAGBWTQAYAfCRwzSRaCshaJGluQ48oDlpBkvH+w4LDsO1ktWapJ8ifBxjK+e3l90SRD9OTlt+jSA2ET4iNH+g9Z2r9ksrMS/YJfTlOX7ZuhY04Vj2f/617G4YlyJ44dcwqEkL0POPsW37skrn+yN9sMB4AKEjxi0ardvZ9JxA/JOuAplLBvct5e53hrje4Nohen1Nb4D6r+d41vXAse6cnx/c/3Sx77hKQDoCcJHDPpwh29NhbP8Z5tONLyfr1F2c0Vsh48/LNwiun6WbszmlF1so+HyMUWSkpRgGqE3+ZuhASBUhI8YbH780L+g01mnOGtzs46G+cOHbnSmzykWLd9+SF5a7at6/J+LT432w4lpOk13ynDf+icvr6b6AaBnCB8x5rOqRtlX3STJiQkyvixPnEqX+dbVWXVlzVjcmKy2qVVmvbjW3P7GxDIZO8C5P+tIme4fenl59V6WWwfQI4SPGLN0q29X1VH9cyUz1bnNj7q2xyl9sgLVj1ii+7f8x3OfmBkuxbnpcs9lI6L9kBzh4pGFZgq1BmSrLwkAQkH4iDFvrttnri8e4dwlvo/u+4ilHoG2do/ZZO3t9ZWmh2Huv02Q3IyUaD8sxwTKqaOKzG0aTwH0hHNPrV2opqlV/rX1oLk9bbTvRd7JtO/jzXUVst6/YFq0VTW0yF3PrpZFmw5IUmKC/O7a081ma+i+6aeXyN8/2iOvr90nP75iVFB7xbiNVtA+3VttluWvrGmS5jaPmR6vi+yNKsmVUwoyHTtbDQg3woeNdCnxZdsOmembKUmJMrwo28xY0QNdd7y7vtKs7zGkb5ac6q8aOJnVs/Lx7qqoPg5teH1tzT756WvrzW7B6SmJ8vCM8XLZaFYzDda5Q/pIYXaa+Tm+vb5Cvjy2ROKJ/i3p7sfPr9oj/9py0PQ0HY9W1C4c3leunzhQzh7k3JlrQDgQPmxqXvzzku0y7/2dUnvU5lvaU3DnRafKjLNKTxhC9EXt8fd3dlpTwenOKPXN1tHeioN1zdKnV5ptX7u13WNWTz3S0CppyYlmB9a8jBQzK0N/lg2t7bL7UIOpJD2/sly2HagPNML+/trTmVYbIv0bvvasUvn9e1vlb8t2xVX4WF1eJfe/8qm5tuSkJ5stEPrlppu/w/rmNtMwrtU+DSY6m0ovGj7+44vDZOLggqg+B0SPvi5ptUz/DyVQEYvP8KEHro92HZGVu46YFTi1XNo7M1VO7ddLLji1b2Bp8O5YuKFS7ntpnXnBUboi6RkDe5vvsXz7YfP+H85fa0rVv/ra2MDiW0dbuu2QrP2s2pyVu2Wxq9zMFLOr6ubKOvPzvtTfL9ATerb512U75Z9bDkpja3u3P08bJW+9YLB8+/zBJqggdNdNLJM5i7bJih2HZXNlbWBatZsrmlo1e9K/Em5mapI5mdDgNW5AriQnJXa5RcKaPVXm/73uiaNL08/483L598mD5J7Lhjtis0X0rLdMT3ze27hf1n1Wbfre6lt8r1cZKUkyoHeGDCzIktNKcmSU/9I/LyOuQkmCN8YWYaipqZHc3Fyprq6WnBx79zXRM/BnPtwtf1+1Rw7WtRz3frq+xlVnDJCvnGC5bT3r/sUbG0xPg9Lx3XunjZRLT+snif4KR3Nbuzy5fLf85p3NZjtyPTPSmRXfPPeUwH1UY0u7XPHHf5ll1W8+9xS5/yujxC1mvbhGnv6gXL7zhcEya9rIkL9OhT/E6X9mix4EdHxdX+j1Z3h01UnPSk8v6y2XjupnfpfZ6TSW2uU7f1spb31aKTdOGig/uXK0uJVWz7775Cr5dK+vb+mqM/rLvZeNkMKc7m97sK+6UR5+Z4s8u7I80Ij9lxvPlLKC7p/kwDkbgr740R5T7Qp2iYHcjBQ5rdgXRDTQD8jPkNLemaZ63lXAjUXBHL/DFj7mzJkjv/rVr6SiokLGjRsnf/jDH+Tss8+OWvhYtGm/3Pz4h4G39aB17pAC80vWg9iBumZZvbtKPth52OzGap0tf2V8iVw+uljG9M+V9NRE2VJZJ/M//syUnFvaPWYti29fMFjuvmTYcc9mdGriD15YE2gm1XDzs+ljTE+IlmlnPrfavJDrWPqb//d8KbBxeCLa9Kzv+89/ImMH5Mord5wX0tfQPprvPbnKDLHoDJVvTBwoXz9zgIwsyukU4vRso6qxVRITEiQrLUlSkxLj6kwikrQC9W+PrTD/R1b88GJX7omzcudhueWJlWb4RF8vdLjuvFP7hPz13ttYKfe8sNYMQRZkpcpjN59la8OzvpTz9x4cPUH8pFwrEzXm5O9gfYs5mdGXleLcDFONGFrYS0YUZ0tRTnqXP9/tB+rk3Q2VMv/jvYHdyFXvzBT50thi0/enoUKHhHX9Jn2NKj/cYHoDNdTqZUtlrbTpUstd0GEaDSBl+Zly5sDecv6wvnJ6aV5MBpKoh49nn31WbrzxRnn00Udl4sSJ8vDDD8vzzz8vmzZtksLCwqiEDy2dTpq9UMaV5sm1Z5WZNQu0KbSrM2xdwfHZD8tl+8ETb4p23tA+ct+XR3Zr23v9MT/1wW75xesbAuU3/WM6VNds3tY/ynnfPLtHL26xaH9tk0z8xUIT6JbPuliKgtwo7/U1++SuZz82jbh6RqAzVPTFANHl8Xjlkt8sNv9HfnzFafLNyYMi+v117HzVriPywY5DJtyrAb0zZdKQAvPC3NODsA6nfu/Jj8yQrL5mzP3GGWaDvZ7SWTH/Pu9Dc8DRIdY5158hF4/sF9LX2nOkQeZ/9JkZgty8v9aEpMyUJDNsrD0mulHihLLenQJ6rPXK6c83Jz0lYrOm9O9W+3H0RPD9rQfNatJNrZ5uVyZGFGXLwIJMc+A/Ut9ivpbOdrLoydHFI/rJ1RMGyBeG9e3282puazcntutNGKk2/68+O9Ioe440mpPco2n/3NUT+suMM0uPO5Qfl+FDA8dZZ50lf/zjH83bHo9HSktL5c4775R77703asMuOvShZ2rdoT+WZdsPycsf75V/bjkge/09HdlpyXLOkAK5fmKZTBnWN+gXOU28P399g7yzodK8gKrS/Az57TXj5UwH7+VyIlc98r58tLtKfjp9tNwQRD+LhsCZz31ifk66t8hvrhnPWHkM+d/lu0y/k1bsltxzYUR+N3oS8cTSnaY5u6LG93/yaDpbTJfL16G2UEKIDsve8/c15u/uohGFJiDY2Sekr0O3P/mRmTWjZ7WzvzpGrjmrtNufrycsD7650fSTHOdkOUB7C277whC55szSqE+L1p/nki0HTGBauu1gYOhbD9haWdbf14yzymzvydIlDP652dd/sXjz/mOG3Pv0SjWbeOpJjZ4c6c9JH6v26+nrtfY1bT9Qf9zKhD5+DXuXjSoygU8rHHaGpf21zVJ+pEG27a+T97cdkn9tOWCqwBbdcfras0tl2ujiqPezRTV8tLS0SGZmprzwwgsyffr0wPtvuukmqaqqkpdffrnT/Zubm82l44PXoBKO8NETDS1t0trmlZyMZFtKm1p61T9oHfLRklysnp3YYe6ibfLLBRvl/FP7yN9umditz5n/8R6zCqn+f//6hAHy4NVjuz1lGZGhZ2sX/mqRCeb3fWmkfOv8wWH9ftq0ee+La8z/G6un54JhfeXUwmzxeL3mILFk84FAZVFflH8yfVS3KpNKXwofXbzd/K1a/R2/vHpslxXSntKGdF3eX4cl1fcvHSa3Xzj0hK8tekB8+oPd8qu3NgWm+E4eWuBvfM0zB1ENNvpzeGf9fnnr0wrzttLhg/+4dJhZIj/SrzXa1K/PUy9WY/7xaJD9v5ecasJST37uWsF+bc1eWbhhv6ludAwOegJ6zuB8M2188tA+pin+ZK/p+re+bX+9bKyoMc9Bf39asRnUN8sMq3T3pNYOre0eE6S0Oq/tBNZT0zVm9Pd72egimTCwd5cnA9p/tHLnERN8dXbgfV8+zT3hY+/evdK/f39ZunSpTJo0KfD+e+65RxYvXiwrVqzodP/7779fHnjggWO+TqyFD4Ru58F6mfLrRaL/vxfO/MJJy4Q6NVbPPPUvU6d1/uKrY1wdzpzs2Q93yw/+vtZUBN/7/hTpm50WlvL8Qws2yd+W7zJv6/f4z6nD5crxJZKW3PkFVg+2897fIX/8x1ZTTtfhTD3zv+OioSeszOgL+n3z1wWaQr99/iDTIB3Ovzt96f3125tkzj+2mbevO7vULNzW1ePU/pOfvLZe1uypNm/rCYtWEvUgc6Iq0TMf7DYzk6zmR52d819fPi3kKqt+TV23R4PgzkP15m3fVPdkyctMMQc0HZ7QoQKd4aHTkq1mXaX3+erp/c1Zus70yEpNMkML2jPxP//cERhC00CgjfcaEII5QdTA9eJHn5lhlY5HNp1ir6tGXzii0ISFcATKaNhX3Wgqdfp3W37Y97NTWr0Z0reXCaT6+zlU3yJ7qxqlsubzE33tO/rwR5fY+jfuqPDhlMoHeuaWeR/Kwo37zSZuP//qmOPe7y9LtsvP39hgbut9f3rlaIJHDNOz8a8+8r45KH5pTLH88frTbW161CbNH83/fCq7htFZl4886ZL4ehB74JVPzTL6anCfLJl91Zgu19nQcKxN0Tr1Xl+H//vLp8nNEexh0bD0wGvrzcFSZ83dcdGpZrhHDxr6c9Wp5dasOg15WsHQ6fjdbTjUgPD/3t8hc97bGqgKaSOkztrp7rICOm1Ym+xfXbO32z0SFv1zOP/UvqaCqTPPjg6MHasLT6/YLb9buCUwrKB/U7MuH2H6ebqih6+Pdh+R51fuMQsJWpUeq7FfQ47+LK19ptzK4/G1CfiC14FOIaMjrR7rbCudbKG7VGuPlJ0VZUcNu0Sy5wPRozNWrvvLcjMD5ZU7Jx9TCteD2EMLNsqflmw3b99y3iBTyqd7P/bpgemqR5aa0rYe4K87u6zHX1P7Gh54db288sneQHP2g1eNkXOHBteQvWDdPvmvlz8NnPnr0J+Oy2tfyJH6VnlnfaW8+PEe09CsZ+F/uP50uWhEaA2gPaF9ZTrMqOP7XdHjg/ZDzPzisJCrS9r8/Zu3N5uzZH3V17Njnfav/Wu65sTRqhta5Y11++SZD8vlkw4Lq+nQiIY4rb7orLLmVo/Ut7SZoSD9HL1OTkowX1Nn8+gMjWCmJus2CP/f25vlyRW7zJCCVq90TyGdJKBn83qw1F4MPdjq8EHHhk/9O9HhMq2udPWc4oHX6zXLSuw+3GD6W3T2js7W0r+bkcU6uzPZ3Q2nOq1Wp9daDadlZWVyxx13RLXhFNGjf2bfemKlqX5oSfXJb50TeBHV/yg/mr/WLLSm7p02Qr5zwWCChwP7evRA8acbJoQ8g0P/Tl5a/Zn85NX15uxXD7raS6JT2UNtptOD4YNvbjBj5Mdr0NTekZ9PHx3UAoN207N2babVdSKsFXm1NK4H3X8/b1C3e1dORmdU/Oz19YH/b+qMsjw5o6y35PdKlZrGNjOLY/m2Q4GZFnrSoE3fN0waaO4XiZMCfZw/f2O9vO/f6ft4dNGuy8cUm+n32ufDEG30xMRUW610/OlPfzIhRKfaPvfcc7Jx40bp1+/EL0qED/fSs89pv/tnYKl1HYPVszGdKqhnzfoiomfO0093x/Ly8UTLvrpejS6upAFEq1Y3Tuq8mN7JaH+AriSqU2iVTmt86GtjZeyAPNsWDHvhoz2ydOtBM1NGF50b2z/XnCnH2rLnOlTS4m9qDAd92ddmzCeW7TRTTo8XyrREr/8f9cBu5/YIwdAVQnW4Z9XOI2Y4TZuL9cRFQ5A2jOolkg2fiOHwoXSarbXI2Pjx4+X3v/+9qYicDOHD3XRhnVv/ujJwZmeZMryvGWuPpTnrCI42bWrvxMurfUMlWnLX5ex17Zrj9WjoyrQLN1bKUyt2B87ENYRqg6guh++WxsBYpmuP6OwJXWdCp6XqgXxQnyyZODjftmoL4kNNLISPUBE+3E/P6nRb+/V7qyUnI8UsxuOGXXzx+QaJOoujwd/cqMWPkcU5Zs2J/Kw0cx8dYtAps5sqawPr3WglX8fq75k6IujF6ABEH+EDQNTPprV/QWdpaE/PiegaFFef0V9mnF1mbgNwJsIHgJihCz6tLj9ien50vYEks/dOstmvYnxZntlDA0B8hQ+6dACElQ6hXJZbzE8ZQADdXAAAIKIIHwAAIKIIHwAAIKIIHwAAIKIIHwAAIKIIHwAAIKIIHwAAIKIIHwAAIKIIHwAAIKIIHwAAIKIIHwAAIKIIHwAAgPABAADcK+Z2tfV6vYGteQEAgDNYx23rOO6o8FFbW2uuS0tLo/1QAABACMfx3NzcE94nwdudiBJBHo9H9u7dK9nZ2ZKQkGB7KtNQU15eLjk5OeJ2PF934/frbvH2+43H51zjsuercUKDR0lJiSQmJjqr8qEPeMCAAWH9HvpLdsMvurt4vu7G79fd4u33G4/POcdFz/dkFQ8Ls10AAEBEET4AAEBExVX4SEtLkx//+MfmOh7wfN2N36+7xdvvNx6fc1qcPd+YbjgFAADuFleVDwAAEH2EDwAAEFGEDwAAEFGEDwAAEFGuCx9z5syRU045RdLT02XixInywQcfnPD+zz//vIwYMcLcf8yYMfLGG2+IW5/vX/7yFzn//POld+/e5nLJJZec9Ofj9N+v5ZlnnjEr5k6fPl3c/Hyrqqrk9ttvl+LiYtNBP2zYMEf9TQf7fB9++GEZPny4ZGRkmJUi7777bmlqahInWLJkiVxxxRVmNUj923zppZdO+jmLFi2SM844w/xuhw4dKvPmzROnCPb5vvjii/LFL35R+vbtaxbgmjRpkrz11lviFKH8fi3vv/++JCcny/jx48WtXBU+nn32WZk5c6aZuvTRRx/JuHHjZOrUqbJ///4u77906VK57rrr5JZbbpGPP/7YHJj0sm7dOnHj89UXLn2+//jHP2TZsmXmxfrSSy+Vzz77TNz4fC07d+6U73//+yZ4OUmwz7elpcW8WOvzfeGFF2TTpk0mcPbv31/c+Hyfeuopuffee839N2zYII899pj5Gj/84Q/FCerr681z1MDVHTt27JAvfelLcuGFF8rq1avlrrvukm9961uOOSAH+3z14K1/zxqeV61aZZ63Hsz1tdqNz7fjCcSNN94oF198sbia10XOPvts7+233x54u7293VtSUuKdPXt2l/e/5pprvF/60pc6vW/ixIne73znO143Pt+jtbW1ebOzs71PPPGE163PV5/jueee6/2f//kf70033eS98sorvU4R7POdO3eud/Dgwd6WlhavEwX7fPW+F110Uaf3zZw50zt58mSv0+hL8fz58094n3vuucc7atSoTu+bMWOGd+rUqV43Pt+unHbaad4HHnjA6+bnO2PGDO99993n/fGPf+wdN26c161cU/nQsz5NxzqU0HGfGH1bz/K7ou/veH+lZ1rHu7/Tn+/RGhoapLW1VfLz88Wtz/cnP/mJFBYWmuqWk4TyfF955RVTmtZhl379+sno0aPlF7/4hbS3t4sbn++5555rPscamtm+fbs5S7788svFjZz8emXXpqO6aZkTXq9C9fjjj5u/Y63muV3MbSwXqoMHD5oXWX3R7Ujf3rhxY5efU1FR0eX99f1ufL5H+8EPfmDGI49+QXPL8/3Xv/5lSvFaonaaUJ6vvmi999578o1vfMMchLdu3Srf+973TMCM9RezUJ7v9ddfbz7vvPPOM7tptrW1yW233eaYYZdgHe/1SndGbWxsNH0vbvbrX/9a6urq5JprrhE32rJlixlG/Oc//2n6PdzONZUPBOfBBx80TZjz5883zX1uo2dIN9xwg+l56NOnj8TLmaFWef785z/LhAkTZMaMGfKjH/1IHn30UXEj7WHSys4jjzxiekS0QfH111+Xn/70p9F+aLCZ9vc88MAD8txzz5m/cbdpb283YVqfozaJxwPXxCs9wCQlJUllZWWn9+vbRUVFXX6Ovj+Y+zv9+XY8g9Dw8e6778rYsWPFCYJ9vtu2bTONl9qg1vHgrPSsQpsxhwwZIm76/eoMl5SUFPN5lpEjR5ozZh3WSE1NFTc93//6r/8yAVObLpXOVtMmv1tvvdWELh22cZPjvV7pTBA3Vz30JEl/xzoz0QlV2lBPllauXGmaae+4447A65VW9PT16u2335aLLrpI3MQ1/zv1hVXP9hYuXBh4n/7y9G0dB++Kvr/j/dU777xz3Ps7/fmqhx56yJwZLliwQM4880xximCfr06fXrt2rRlysS5f+cpXAjMFdKaP236/kydPNkMtVshSmzdvNqEkloNHqM9Xe5aODhhW8HLjllVOfr0K1dNPPy3f/OY3zbXO9HGrnJycY16vdAhRp5HrbZ127jpeF3nmmWe8aWlp3nnz5nnXr1/vvfXWW715eXneiooK8/EbbrjBe++99wbu//7773uTk5O9v/71r70bNmww3cUpKSnetWvXet34fB988EFvamqq94UXXvDu27cvcKmtrfW68fkezWmzXYJ9vrt37zazl+644w7vpk2bvK+99pq3sLDQ+7Of/czrxuer/1/1+T799NPe7du3e99++23vkCFDzCw2J9D/dx9//LG56Evxb37zG3N7165d5uP6XPU5W/Q5ZmZmev/zP//TvF7NmTPHm5SU5F2wYIHXjc/3ySefNK/P+jw7vl5VVVV53fh8j+b22S6uCh/qD3/4g7esrMwcZHXq3vLlywMf+8IXvmAOQB0999xz3mHDhpn76zS2119/3evW5ztw4EDzn+Doi/6Ru/X36+TwEcrzXbp0qZkurgdxnXb785//3Ew3duPzbW1t9d5///0mcKSnp3tLS0u93/ve97xHjhzxOsE//vGPLv8/Ws9Rr/U5H/0548ePNz8f/f0+/vjjXqcI9vnq7RPd342/33gKHwn6T7SrLwAAIH64pucDAAA4A+EDAABEFOEDAABEFOEDAABEFOEDAABEFOEDAABEFOEDAABEFOEDAABEFOEDAABEFOEDAABEFOEDAABEFOEDAABIJP3/Xj+mAKUuAlYAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import matplotlib.pyplot as plt\n", "plt.plot( obs['band6']['r'], obs['band6']['Tb'] )" ] }, { "cell_type": "markdown", "id": "5319f71c-2766-42bd-ad40-5eeabe087f4a", "metadata": {}, "source": [ "Those settings are identical to the visibility-based analysis" ] }, { "cell_type": "code", "execution_count": 22, "id": "e55218c2-565e-4a6f-9a51-f313420c837e", "metadata": {}, "outputs": [], "source": [ "dust_opacity = jnp.load(dsharp_opac.get_datafile('default_opacities_smooth.npz'))\n", "\n", "disk = frap.model( incl= 6.0,\n", " r_out = 1.3,\n", " N_GP = 130,\n", " flux_uncert = True )\n", "\n", "D = 114.9 # pc\n", "\n", "disk.set_parameter('log10_T', \n", " free = False, \n", " profile = lambda r: jnp.log10(77.0*(r*D/10)**(-0.5)) )\n", "\n", "\n", "disk.set_parameter('q', \n", " free = False, \n", " profile = lambda r: 3.5 )\n", "\n", "\n", "lengthscale = 0.03\n", "\n", "disk.set_parameter('log10_Sigma_d', free = True,\n", " bounds = ( -6.0, 3.0 ), \n", " lengthscale = lengthscale)\n", "\n", "disk.set_parameter('log10_a_max', free = True, \n", " bounds = ( -2.0, 2.0 ),\n", " lengthscale = lengthscale)\n" ] }, { "cell_type": "markdown", "id": "baad9efd-7c3e-4b02-9164-801c7a6aba33", "metadata": {}, "source": [ "Then you can use `set_radialprofile()` to set your radial profile. Here you also have to give the covariance matrix." ] }, { "cell_type": "code", "execution_count": 24, "id": "03f1ff66-a3e9-49e9-b7d2-e2a2c84f11b9", "metadata": {}, "outputs": [], "source": [ "f_s = { 'band3': 0.025, 'band6': 0.05 }\n", "f_mean = { 'band3': 1.0, 'band6': 1.0 }\n", " \n", "for band in bands:\n", " disk.set_radialprofile( band = band,\n", " r = [obs[band]['r']],\n", " Tb = [obs[band]['Tb']],\n", " nu = [obs[band]['nu']],\n", " Nch = 1,\n", " f_s = f_s[band],\n", " f_mean = f_mean[band],\n", " FWHM = 0.05, \n", " Cov = [obs[band]['Cov']] )" ] }, { "cell_type": "markdown", "id": "a94be697-d804-4452-8830-dc75388608ad", "metadata": {}, "source": [ "Remaining part is the same as before!" ] }, { "cell_type": "code", "execution_count": 25, "id": "2dc9525f-3984-4bd3-aa23-8de77785fabe", "metadata": {}, "outputs": [], "source": [ "disk.set_opacity( opac_dict = dust_opacity )" ] }, { "cell_type": "code", "execution_count": 26, "id": "6eed7ca1-0533-45d6-8a65-c14391e99801", "metadata": {}, "outputs": [], "source": [ "inference = frap.inference(disk)" ] }, { "cell_type": "code", "execution_count": 27, "id": "14ec02bb-1bf0-4ca3-b59d-4c3df4384319", "metadata": { "scrolled": true }, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "6d114fa969774f56a85620e103207d65", "version_major": 2, "version_minor": 0 }, "text/plain": [ " 0%| | 0/500 [00:00" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plotter.sample_paths( 'posterior_f', nskip = 10, plot_kwargs = { 'alpha' : 0.01, 'color' : 'b' } )\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "4fb9a6e6-af54-4346-baae-068b29398767", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.13.12" } }, "nbformat": 4, "nbformat_minor": 5 }