view_output.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from microlensing.Util.util import *\n",
    "\n",
    "import mpl_style\n",
    "plt.style.use(mpl_style.style1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "is_ipm = True\n",
    "is_double = False\n",
    "prefix = ''\n",
    "if is_ipm:\n",
    "    prefix = f'{prefix}ipm_'\n",
    "else:\n",
    "    prefix = f'{prefix}irs_'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "params = read_params(f'{prefix}parameter_info.txt')\n",
    "params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stars, _, _, _ = read_stars(f'{prefix}stars.bin', np.float32)\n",
    "print(stars.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.scatter(stars[:, 0], stars[:, 1], s=stars[:,2])\n",
    "\n",
    "ax.set_xlabel('$x_1 / \\\\theta_★$')\n",
    "ax.set_ylabel('$x_2 / \\\\theta_★$')\n",
    "\n",
    "ax.set_aspect(1)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dat = read_array(f'{prefix}magnifications.bin', np.float32)\n",
    "# magnification\n",
    "if is_ipm:\n",
    "    mu = dat\n",
    "else:\n",
    "    mu = dat / params['num_rays_y']\n",
    "print(np.min(mu))\n",
    "print(np.max(mu))\n",
    "print(np.mean(mu))\n",
    "# astronomical magnitudes\n",
    "mags = -2.5*np.log10(mu / np.abs(params['mu_ave']))\n",
    "print(np.max(mags))\n",
    "print(np.min(mags))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.imshow(mags[:,:], extent = [params['center_y1'] - params['half_length_y1'], \n",
    "                               params['center_y1'] + params['half_length_y1'], \n",
    "                               params['center_y2'] - params['half_length_y2'], \n",
    "                               params['center_y2'] + params['half_length_y2']],\n",
    "                               cmap = 'viridis_r')\n",
    "\n",
    "ax.set_xlabel('$y_1 / \\\\theta_★$')\n",
    "ax.set_ylabel('$y_2 / \\\\theta_★$')\n",
    "\n",
    "ax.set_aspect(params['half_length_y1'] / params['half_length_y2'])\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "y1, step = np.linspace(params['center_y1'] - params['half_length_y1'],\n",
    "                       params['center_y1'] + params['half_length_y1'],\n",
    "                       params['num_pixels_y1'],\n",
    "                       endpoint = False, retstep=True)\n",
    "y1 += step/2\n",
    "\n",
    "y2, step = np.linspace(params['center_y2'] - params['half_length_y2'],\n",
    "                       params['center_y2'] + params['half_length_y2'],\n",
    "                       params['num_pixels_y2'],\n",
    "                       endpoint = False, retstep=True)\n",
    "y2 += step/2\n",
    "\n",
    "where = slice(0,1000)\n",
    "\n",
    "ax.plot(y1[where], mags[100, where], label = '$y_1$')\n",
    "ax.plot(y2[where], mags[where, 100], label = '$y_2$')\n",
    "\n",
    "ax.legend()\n",
    "\n",
    "ax.set_xlabel('$y_i / \\\\theta_★$')\n",
    "ax.set_ylabel('microlensing $\\\\Delta m$ (magnitudes)')\n",
    "\n",
    "ax.invert_yaxis()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# magnification\n",
    "if is_ipm:\n",
    "    dat, num_pixels = read_hist(f'{prefix}mags_numpixels.txt').T\n",
    "    mu = dat / 1000\n",
    "else:\n",
    "    dat, num_pixels = read_hist(f'{prefix}numrays_numpixels.txt').T\n",
    "    mu = dat / params['num_rays_y']\n",
    "# astronomical magnitudes\n",
    "mags = -2.5*np.log10(mu / np.abs(params['mu_ave']))\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.hist(mags, weights=num_pixels, density=True, bins = 1000)\n",
    "\n",
    "ax.set_xlabel(\"microlensing $\\\\Delta m$ (magnitudes)\")\n",
    "ax.set_ylabel(\"$p(\\\\Delta m)$\")\n",
    "\n",
    "ax.invert_xaxis()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "microlensing",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.5"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}