{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "blur angles: [0.00892833 0.01785525 0.03569911 0.07130746]\n",
      "running parallel simulations... \n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "from multiprocessing import Process, Pool\n",
    "\n",
    "# Simulation settings\n",
    "concentrations = np.logspace(-5,0,num=300) # #/mm^3\n",
    "focal_lengths = np.array([25, 50, 100, 200]) # mm\n",
    "focal_distance = 500 # mm\n",
    "pixel_size = 0.0045 # mm\n",
    "particle_area = 0.03**2 # mm^2\n",
    "f_number = 2.8\n",
    "apertures = focal_lengths / f_number # mm\n",
    "blur_angles = np.arctan(apertures/focal_distance/2)\n",
    "print('blur angles: ' + str(blur_angles))\n",
    "particle_contrast = 0.5\n",
    "\n",
    "\n",
    "# Finite iluminated volume case\n",
    "def finite_volume(counter, dx=100, dy=100, dz=100):\n",
    "    # output data\n",
    "    \n",
    "    print(counter)\n",
    "    \n",
    "    results = np.zeros((len(blur_angles), len(concentrations)))\n",
    "\n",
    "    for col, con in enumerate(concentrations):\n",
    "        n_particles = int(dx*dy*dz*con)\n",
    "        z_vals = focal_distance - dz/2 + dz*np.random.random(n_particles)\n",
    "\n",
    "        for row, b in enumerate(blur_angles):\n",
    "\n",
    "            x_pix = focal_lengths[row] *(-dx/2 + dx*np.random.random(n_particles)) / z_vals\n",
    "            y_pix = focal_lengths[row] *(-dy/2 + dy*np.random.random(n_particles)) / z_vals\n",
    "\n",
    "            pix_rad = np.sqrt(x_pix**2 + y_pix**2)\n",
    "\n",
    "            blur_rad = np.abs(z_vals-focal_distance)*np.tan(b)\n",
    "\n",
    "            overlapping_signals = (blur_rad >= pix_rad).astype('uint8')\n",
    "\n",
    "            blur_area = np.pi*(blur_rad)**2\n",
    "            blur_signal = particle_contrast*particle_area/(blur_area + particle_area)\n",
    "            results[row, col] = (particle_contrast - np.sum(blur_signal*overlapping_signals))/particle_contrast\n",
    "\n",
    "    return results\n",
    "\n",
    "\n",
    "\n",
    "results = np.zeros((len(blur_angles), len(concentrations)))\n",
    "\n",
    "iterations = 10\n",
    "print('running parallel simulations... ')\n",
    "output = Pool(10).map(finite_volume, range(0,iterations))\n",
    "print('done.')\n",
    "for out in output:\n",
    "    results += out\n",
    "results /= iterations\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.loglog(concentrations*1000, np.transpose(results))\n",
    "plt.xlabel('Particle Concentration (#/ml)')\n",
    "plt.ylabel('Relative Contrast')\n",
    "plt.legend([str(x) + ' mm' for x in focal_lengths])"
   ]
  },
  {
   "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}