{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "

Simple simulations and convolution with instrumental broadening" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "In /Users/holtz/anaconda3/lib/python3.7/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: \n", "The text.latex.preview rcparam was deprecated in Matplotlib 3.3 and will be removed two minor releases later.\n", "In /Users/holtz/anaconda3/lib/python3.7/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: \n", "The mathtext.fallback_to_cm rcparam was deprecated in Matplotlib 3.3 and will be removed two minor releases later.\n", "In /Users/holtz/anaconda3/lib/python3.7/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: Support for setting the 'mathtext.fallback_to_cm' rcParam is deprecated since 3.3 and will be removed two minor releases later; use 'mathtext.fallback : 'cm' instead.\n", "In /Users/holtz/anaconda3/lib/python3.7/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: \n", "The validate_bool_maybe_none function was deprecated in Matplotlib 3.3 and will be removed two minor releases later.\n", "In /Users/holtz/anaconda3/lib/python3.7/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: \n", "The savefig.jpeg_quality rcparam was deprecated in Matplotlib 3.3 and will be removed two minor releases later.\n", "In /Users/holtz/anaconda3/lib/python3.7/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: \n", "The keymap.all_axes rcparam was deprecated in Matplotlib 3.3 and will be removed two minor releases later.\n", "In /Users/holtz/anaconda3/lib/python3.7/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: \n", "The animation.avconv_path rcparam was deprecated in Matplotlib 3.3 and will be removed two minor releases later.\n", "In /Users/holtz/anaconda3/lib/python3.7/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: \n", "The animation.avconv_args rcparam was deprecated in Matplotlib 3.3 and will be removed two minor releases later.\n" ] } ], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

Spectrum simulation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's start with some simple simulations of a spectrum: flux as a function of wavelength. We'll start with simulating a pure emission line spectrum, where we'll assume intrinsic line width to be much less than our instrumental line width." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's further assume that our instrumental broadening profile can be represented by a Gaussian. Start by writing a function to return a Gaussian given a set of input wavelengths, a wavelength center, and the $\\sigma$ corresponding to the broadening." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "def gauss(x,mean=0,sig=1) :\n", " \"\"\" Return normalized Gaussian function here\"\"\"\n", "\n", " return 1/np.sqrt(2*np.pi)/sig*np.exp(-0.5*(x-mean)**2/sig**2) # add code here" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Suppose we have a spectrograph that is chararacterized by a resolution, $R\\equiv \\frac{\\lambda}{\\Delta\\lambda} = 2000$, where $\\Delta\\lambda$ refers to the full-width half-maximum (FWHM) of the broadening profile. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Suppose we are going to simulate a line at $\\lambda=6563$ A. Given the relation between the FWHM and the $\\sigma$ of the Gaussian (which you should easily be able to derive, if you don't already know it), what is the $\\sigma$ corresponding to this resolution? Write a function to return it" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "R: 2000.000000 sig: 1.393418 fwhm: 3.281500\n" ] } ], "source": [ "def get_sigma(R,lam) :\n", " fwhm=lam/R # compute fwhm in wavelength\n", " sig=fwhm/2.355 # compute Gaussian sigma in wavelength\n", " return sig\n", "\n", "R=2000\n", "lam=6563\n", "sig=get_sigma(R,lam)\n", "print('R: {:f} sig: {:f} fwhm: {:f}'.format(R,sig, lam/R))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For our simulation, we'll create a wavelength array and a flux array. Given the FWHM of the instrumental profile, what would be a good choice for a wavelength spacing (the dispersion) of the spectrum, given the sampling theorem, which says you should have at least 2 samples per FWHM?" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "disp= 1 #set an appropriate dispersion (A/pixel)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Create a wavelength array to run between 6000 and 7000 A with your chosen dispersion. Also create an empty flux array:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAO0ElEQVR4nO3cbYxcZ3nG8f9Vu45E2+CkdsDYhnVbg+pWFYSpFURf1IRQO01jPjpqGyuoskBNBagUnEbqy7cAVakiokQWUCUCEQGFxq2MkpC+fGpe1iEJdYybJYV6sUkWpAJqJCKLux/mWEy2s95dz9jrnef/k0ZzzvPcZ85zW+u5fM7OOFWFJKldP7HSC5AkrSyDQJIaZxBIUuMMAklqnEEgSY1bu9ILOBcbNmyoqamplV6GJK0qR44c+U5VbZw/viqDYGpqiunp6ZVehiStKkm+OWzcW0OS1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMYZBJLUOINAkhpnEEhS4wwCSWqcQSBJjTMIJKlxBoEkNc4gkKTGGQSS1LixBEGSXUmOJ5lJcmDIfJLc0c0/neTKefNrknwlyT+NYz2SpKUbOQiSrAHuBHYDO4Abk+yYV7Yb2N499gN3zZt/D3Bs1LVIkpZvHFcEO4GZqnquql4C7gP2zKvZA9xbfY8A65NsAkiyBfgd4ONjWIskaZnGEQSbgRMD+7Pd2FJr/hb4APCjs50kyf4k00mm5+bmRlqwJOnHxhEEGTJWS6lJcj3wQlUdWewkVXWwqnpV1du4ceO5rFOSNMQ4gmAW2DqwvwU4ucSatwI3JPkG/VtKVyf51BjWJElaonEEwePA9iTbkqwD9gKH5tUcAm7qPj10FfC9qjpVVbdW1ZaqmuqO++eq+v0xrEmStERrR32Bqjqd5BbgAWAN8MmqOprkXd383cBh4DpgBngRuHnU80qSxiNV82/nX/x6vV5NT0+v9DIkaVVJcqSqevPH/WaxJDXOIJCkxhkEktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMYZBJLUOINAkhpnEEhS4wwCSWqcQSBJjTMIJKlxBoEkNc4gkKTGGQSS1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQSFLjDAJJatxYgiDJriTHk8wkOTBkPknu6OafTnJlN741yb8kOZbkaJL3jGM9kqSlGzkIkqwB7gR2AzuAG5PsmFe2G9jePfYDd3Xjp4E/qapfBK4C/mjIsZKk82gcVwQ7gZmqeq6qXgLuA/bMq9kD3Ft9jwDrk2yqqlNV9QRAVf0AOAZsHsOaJElLNI4g2AycGNif5f+/mS9ak2QKeBPw6BjWJElaonEEQYaM1XJqkvw08PfAe6vq+0NPkuxPMp1kem5u7pwXK0l6uXEEwSywdWB/C3ByqTVJfpJ+CHy6qr6w0Emq6mBV9aqqt3HjxjEsW5IE4wmCx4HtSbYlWQfsBQ7NqzkE3NR9eugq4HtVdSpJgE8Ax6rqb8awFknSMq0d9QWq6nSSW4AHgDXAJ6vqaJJ3dfN3A4eB64AZ4EXg5u7wtwJ/AHw1yZPd2J9V1eFR1yVJWppUzb+df/Hr9Xo1PT290suQpFUlyZGq6s0f95vFktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMYZBJLUOINAkhpnEEhS4wwCSWqcQSBJjTMIJKlxBoEkNc4gkKTGGQSS1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXEGgSQ1bixBkGRXkuNJZpIcGDKfJHd0808nuXKpx0qSzq+RgyDJGuBOYDewA7gxyY55ZbuB7d1jP3DXMo6VJJ1Ha8fwGjuBmap6DiDJfcAe4JmBmj3AvVVVwCNJ1ifZBEwt4dix+at/PMozJ79/Pl5aki6IHa+5lL/43V8a62uO49bQZuDEwP5sN7aUmqUcC0CS/Ummk0zPzc2NvGhJUt84rggyZKyWWLOUY/uDVQeBgwC9Xm9ozWLGnaKSNAnGEQSzwNaB/S3AySXWrFvCsZKk82gct4YeB7Yn2ZZkHbAXODSv5hBwU/fpoauA71XVqSUeK0k6j0a+Iqiq00luAR4A1gCfrKqjSd7Vzd8NHAauA2aAF4Gbz3bsqGuSJC1d+h/kWV16vV5NT0+v9DIkaVVJcqSqevPH/WaxJDXOIJCkxhkEktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMYZBJLUOINAkhpnEEhS4wwCSWqcQSBJjTMIJKlxBoEkNc4gkKTGGQSS1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQSFLjDAJJatxIQZDk8iQPJXm2e75sgbpdSY4nmUlyYGD8I0m+luTpJF9Msn6U9UiSlm/UK4IDwMNVtR14uNt/mSRrgDuB3cAO4MYkO7rph4BfrqpfAf4TuHXE9UiSlmnUINgD3NNt3wO8Y0jNTmCmqp6rqpeA+7rjqKoHq+p0V/cIsGXE9UiSlmnUIHhVVZ0C6J6vGFKzGTgxsD/bjc33TuBLI65HkrRMaxcrSPJl4NVDpm5b4jkyZKzmneM24DTw6bOsYz+wH+C1r33tEk8tSVrMokFQVW9baC7J80k2VdWpJJuAF4aUzQJbB/a3ACcHXmMfcD1wTVUVC6iqg8BBgF6vt2CdJGl5Rr01dAjY123vA+4fUvM4sD3JtiTrgL3dcSTZBXwQuKGqXhxxLZKkczBqENwOXJvkWeDabp8kr0lyGKD7ZfAtwAPAMeCzVXW0O/5jwM8ADyV5MsndI65HkrRMi94aOpuq+i5wzZDxk8B1A/uHgcND6n5hlPNLkkbnN4slqXEGgSQ1ziCQpMYZBJLUOINAkhpnEEhS4wwCSWqcQSBJjTMIJKlxBoEkNc4gkKTGGQSS1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMYZBJLUOINAkhpnEEhS4wwCSWrcSEGQ5PIkDyV5tnu+bIG6XUmOJ5lJcmDI/PuTVJINo6xHkrR8o14RHAAerqrtwMPd/sskWQPcCewGdgA3JtkxML8VuBb47xHXIkk6B6MGwR7gnm77HuAdQ2p2AjNV9VxVvQTc1x13xkeBDwA14lokSedg1CB4VVWdAuierxhSsxk4MbA/242R5AbgW1X11GInSrI/yXSS6bm5uRGXLUk6Y+1iBUm+DLx6yNRtSzxHhoxVkld0r/H2pbxIVR0EDgL0ej2vHiRpTBYNgqp620JzSZ5PsqmqTiXZBLwwpGwW2DqwvwU4Cfw8sA14KsmZ8SeS7Kyqby+jB0nSCEa9NXQI2Ndt7wPuH1LzOLA9ybYk64C9wKGq+mpVXVFVU1U1RT8wrjQEJOnCGjUIbgeuTfIs/U/+3A6Q5DVJDgNU1WngFuAB4Bjw2ao6OuJ5JUljsuitobOpqu8C1wwZPwlcN7B/GDi8yGtNjbIWSdK58ZvFktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMYZBJLUOINAkhpnEEhS4wwCSWqcQSBJjTMIJKlxBoEkNc4gkKTGGQSS1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxqWqVnoNy5ZkDvjmOR6+AfjOGJezGthzG+y5DaP0/Lqq2jh/cFUGwSiSTFdVb6XXcSHZcxvsuQ3no2dvDUlS4wwCSWpci0FwcKUXsALsuQ323Iax99zc7wgkSS/X4hWBJGmAQSBJjZuIIEiyPsnnk3wtybEkb0lyeZKHkjzbPV82UH9rkpkkx5P89sD4m5N8tZu7I0lWpqPFLdDzR7r9p5N8Mcn6gfqJ7Hlg7v1JKsmGgbGJ7TnJH3d9HU3y4YH6iew5yRuTPJLkySTTSXYO1K/qnpO8oevrzOP7Sd57Qd/DqmrVP4B7gD/sttcB64EPAwe6sQPAh7rtHcBTwCXANuDrwJpu7jHgLUCALwG7V7q3Zfb8dmBtN/ahFnrutrcCD9D/kuGGSe8Z+C3gy8Al3fgVDfT84Jk1A9cB/zpJPQ/0vgb4NvC6C/ketuqvCJJcCvwG8AmAqnqpqv4H2EP/B4ru+R3d9h7gvqr6YVX9FzAD7EyyCbi0qv69+n+i9w4cc1FZqOeqerCqTndljwBbuu2J7bmb/ijwAWDwkw+T3PO7gdur6ofd+AvdIZPccwGXdmWvBE5226u+53muAb5eVd/kAr6HrfogAH4OmAP+LslXknw8yU8Br6qqUwDd8xVd/WbgxMDxs93Y5m57/vjFaKGeB72T/r8IYIJ7TnID8K2qempe/cT2DLwe+PUkjyb5tyS/2tVPcs/vBT6S5ATw18CtXf0k9DxoL/CZbvuCvYdNQhCsBa4E7qqqNwH/S/8yaiHD7pnVWcYvRmftOcltwGng02eGhrzGJPT8l8BtwJ8PqZ/Ung9045cBVwF/Cny2uxc8yT2/G3hfVW0F3kd3xcBk9AxAknXADcDnFisdMjZSz5MQBLPAbFU92u1/nv4P0vPdpRLd8wsD9VsHjt9C/zJzlh/fShkcvxgt1DNJ9gHXA7/XXR6eqZ/UnrcBTyX5Bv31P5Hk1Ux2z7PAF6rvMeBH9P8jsknueR/whW7sc8DOgfrV3vMZu4Enqur5bv+CvYet+iCoqm8DJ5K8oRu6BngGOET/h4fu+f5u+xCwN8klSbYB24HHukuvHyS5qvvX1U0Dx1xUFuo5yS7gg8ANVfXiwCGT2vMTVXVFVU1V1RT9vwhXdrWT2vMzwD8AVwMkeT39X6h+h8nu+STwm93Y1cCz3faq73nAjfz4thBcyPewlfwN+bgewBuBaeBp+n9JLgN+FniY/g/Mw8DlA/W30f9N+3EGfqsO9ID/6OY+RvfN64vxsUDPM/TvHT7ZPe6e9J7nzX+D7lNDk9wz/Tf+T3U9PAFc3UDPvwYcof9pmUeBN09Yz68Avgu8cmDsgr2H+V9MSFLjVv2tIUnSaAwCSWqcQSBJjTMIJKlxBoEkNc4gkKTGGQSS1Lj/A+/3oq3FbdlAAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "waves=np.arange(6000,7000,disp)\n", "flux=np.zeros(len(waves))\n", "plt.plot(waves,flux)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now add in fluxes to your wavelength array given a line centered at 6562.8 Angstroms (H$\\alpha$), using your gaussian routine and your wavelength, and plot your result. You will need to choose some amplitude for the line: start with a total flux of 1000. If you'd like, add in some other lines, e.g. the [NII] lines at 6548.3 and 6583.41 (which come in a fixed flux ratio of 1:3 from atomic physics!). You can also play around with your resolution: what resolution is required to resolve these three lines?" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAXKklEQVR4nO3db4xc13nf8e8zu0uuJFuRZFGMQNGmijJu5aCVE5ax4bpJpSKSo8BSXwig0SZE60BAoAR20DaV6qJpXwhQ7KItitYthPwpA6cRmMSthCRFrTJW0gC2VEq27EiyIKb6x0oWKSeKI9kkd2eevpg73LuzM7vL3XNnZy+/H2Axd+6c2T1HWv54+My550ZmIklql85Wd0CSVJ7hLkktZLhLUgsZ7pLUQoa7JLXQ7FZ3AODqq6/Offv2bXU3JGlbeeKJJ97IzF2jXpuKcN+3bx/Hjx/f6m5I0rYSES+Ne82yjCS1kOEuSS1kuEtSCxnuktRChrsktZDhLkktZLhLUgsZ7tIWe+lbb/O/nz+91d1Qy0zFRUzSxeyHP/MoAC/ef9vWdkSt4sxdklpo3eEeETMR8ZWI+J3q+VUR8UhEPF89Xllre29EnIiI5yLiliY6LrWNd0VTSRcyc/8E8Gzt+T3AsczcDxyrnhMRNwCHgPcBtwKfjYiZMt2V2muha7irnHWFe0RcB9wG/FLt9O3Aker4CHBH7fyDmXk2M18ATgAHi/RWarEzi92t7oJaZL0z938H/DzQq53bnZmvAVSP11Tn9wCv1NqdrM5JWsWZBcNd5awZ7hHx48CpzHxind8zRpxb8e/NiLgrIo5HxPHTp10GpotXp/oTc3aht3pD6QKsZ+b+IeCjEfEi8CBwU0R8Dng9Iq4FqB5PVe1PAntr778OeHX4m2bmA5l5IDMP7No1cq956aKwc7b/kdRZyzIqaM1wz8x7M/O6zNxH/4PS38/Mvw88DByumh0GHqqOHwYORcTOiLge2A88XrznUkvsmO3/MTzjzF0FbeYipvuBoxHxceBl4E6AzHw6Io4CzwCLwN2Z6ZREGmNupl+XOdc13FXOBYV7Zj4KPFodfwu4eUy7+4D7Ntk36aLiMneV5BWq0tQw3VWO4S5NCWfuKslwl6aE2a6SDHdpSjhzV0mGuzQl3DhMJRnu0pQw2lWS4S5NiZ4zdxVkuEtbrtpcxmxXQYa7NCXMdpVkuEtTwqqMSjLcpSmRzt1VkOEuTQln7irJcJemhNmukgx3aUp4EZNKMtylKWG0qyTDXZoWprsKMtylKeFqGZVkuEtTwpK7SjLcpSlhuKskw12aEma7SjLcpSnhUkiVZLhLU8JoV0mGuzQlnLirJMNdmhqmu8ox3KUp4cxdJRnu0pQw21WS4S5NCWfuKslwl6aE2w+oJMNdmhLO3FWS4S5NCbNdJRnu0pbrx7pXqKokw12SWshwl6aEE3eVZLhLU8LVMirJcJemhDN3lWS4S1PCcFdJhrs0Jcx2lWS4S1PCpZAqyXCXpoTRrpIMd2lamO4qaM1wj4j5iHg8Ip6KiKcj4l9V56+KiEci4vnq8crae+6NiBMR8VxE3NLkAKS2cCmkSlrPzP0scFNm/nXgRuDWiPgAcA9wLDP3A8eq50TEDcAh4H3ArcBnI2Kmgb5LrWLJXSWtGe7Z91b1dK76SuB24Eh1/ghwR3V8O/BgZp7NzBeAE8DBkp2W2shsV0nrqrlHxExEfBU4BTySmY8BuzPzNYDq8Zqq+R7gldrbT1bnhr/nXRFxPCKOnz59ehNDkNrBmbtKWle4Z2Y3M28ErgMORsT3r9I8Rn2LEd/zgcw8kJkHdu3ata7OSm3WM91V0AWtlsnMN4FH6dfSX4+IawGqx1NVs5PA3trbrgNe3WxHpfbqz4eMdpW0ntUyuyLiiur4EuDvAN8AHgYOV80OAw9Vxw8DhyJiZ0RcD+wHHi/cb6l9nLmroNl1tLkWOFKteOkARzPzdyLiS8DRiPg48DJwJ0BmPh0RR4FngEXg7szsNtN9qT2MdpW0Zrhn5teA9484/y3g5jHvuQ+4b9O9ky4iTtxVkleoSlPCvWVUkuEuTcg/+NXH+dD9vz/2daNdJa2n5i6pgC8+t/r1HE7cVZIzd2lKmO0qyXCXpoQ1d5VkuEtSCxnu0pRw4q6SDHdpSrifu0oy3KUp4cxdJRnu0pQw21WS4S5NCWfuKslwl6aENXeVZLhLU8KZu0oy3CWphQx3acv1p+xeoaqSDHdpSpjtKslwl6aE2a6SDHdpSjhzV0mGuzQlXAqpkgx3aUo4c1dJhru0xQahbrarJMNd2mLnQ92puwoy3KUtNljfbrSrJMNd2mKDUHfirpIMd2kCVrv6dKnmbrqrHMNdmoDeKrl9vixjtqsgw12agFVn7kOPUgmGuzQB9eBeEfSDsozproIMd2kCerXkHi7RLM3cTXeVY7hLE1CflfeGpujpVUxqgOEuTcCq4T70KJVguEsTUC+5rCi5n6+5G+8qx3CXJmD1mbtLIVWe4S5NQD3Qx83cV1sLL10ow12agHpuj6+5m+4qx3CXJmB5WWb4xZVtpM0y3KUJyGVlmdE1d6kkw12agNVm7q6WURMMd2kC1ldzl8pZM9wjYm9EfDEino2IpyPiE9X5qyLikYh4vnq8svaeeyPiREQ8FxG3NDkAaTtYvv3A6CtUnbirpPXM3BeBf5SZfxX4AHB3RNwA3AMcy8z9wLHqOdVrh4D3AbcCn42ImSY6L20Xy4LbvWU0AWuGe2a+lplPVsd/ATwL7AFuB45UzY4Ad1THtwMPZubZzHwBOAEcLNxvaVvJ1TYOc7WMGnBBNfeI2Ae8H3gM2J2Zr0H/LwDgmqrZHuCV2ttOVueGv9ddEXE8Io6fPn16A12Xto/Vau6j2kibte5wj4h3AL8NfDIzv71a0xHnVvzeZuYDmXkgMw/s2rVrvd2QtqVx2w/kKleuSpuxrnCPiDn6wf7rmfn56vTrEXFt9fq1wKnq/Elgb+3t1wGvlumutD2N235geaCb7ipnPatlAvhl4NnM/De1lx4GDlfHh4GHaucPRcTOiLge2A88Xq7L0vYzriyz/A5NE+uOLgKz62jzIeAngK9HxFerc/8MuB84GhEfB14G7gTIzKcj4ijwDP2VNndnZrd0x6XtZNwHqpZl1JQ1wz0z/4jRdXSAm8e85z7gvk30S2qVHBfo9TaWZVSQV6hKEzBu+4HVtgKWNsNwlyZgPStkzHaVZLhLE7D8A9UxbUx3FWS4SxOQY/aWWT5zN91VjuEuTcDYOjvWZdQMw12aiLUvYjLbVZLhLk3AegLdm3WoJMNdmoDxyx9Hr3mXNstwlyagXlt3+wFNguEuTUCvVzu25q4JMNylCUhGl2KWLZZx6q6CDHdpAsZtP7As9CfYH7Wf4S5NwNiNw1a5t6q0GYa7NAHLP1Ctnx/dRtosw12agPEzd3eFVDMMd2kCermOmbvhroIMd2kCxt5mz43D1BDDXZqA5atlRm8c5sxdJRnu0gSsZ4WM2a6SDHdpAsatirHmrqYY7tIE9Gqfota3Ilge6Ka7yjHcpQkY+4GqNXc1xHCXGvDUK2+y757f5cmX/wxYZfsBa+5qiOEuNeCPTrwBwCPPvA6scuFS7T1uHKaSDHepATtn+3+0zi70C+zLyzJLxznm4iZpswx3qQE752YAOLPYBVZZ525ZRg0x3KUGzFcz9zML/XAfd/elOssyKslwlxowX83cR5Vl3CxMk2C4Sw2Y7QQweubuUkhNguEuNaBbJfW5bnXFUr3mPuYiJjcOU0mGu9SA7tDSl+U363D7ATXPcJcaMAj3iH55Zuxs3fq7GmK4Sw0YnrmPr7lTOzbdVY7hLjWgl8PhvnTsyhlNguEuNaDbW/48x83cvYhJDTHcpQYMVssMQr0+cx83izfdVZLhLjWgt8pqmXGlGGvuKslwlxqwuOID1dpxr/6BqjV3NcNwlxrQG1oKOW73R2vuasqa4R4RvxIRpyLij2vnroqIRyLi+erxytpr90bEiYh4LiJuaarj0jTrrlgts/aFS24cppLWM3P/L8CtQ+fuAY5l5n7gWPWciLgBOAS8r3rPZyNiplhvpW1ixTr3ZRcxjSnLNN4rXUzWDPfM/EPgT4dO3w4cqY6PAHfUzj+YmWcz8wXgBHCwTFel7WP4A9WxFzGNKdFIm7XRmvvuzHwNoHq8pjq/B3il1u5kdW6FiLgrIo5HxPHTp09vsBvSdBouy4y7h+qyNg32Rxef0h+oxohzI39nM/OBzDyQmQd27dpVuBvS1hrM3JfWua89c3fqrpI2Gu6vR8S1ANXjqer8SWBvrd11wKsb7560PQ2WQg5q7+NXyFhzVzM2Gu4PA4er48PAQ7XzhyJiZ0RcD+wHHt9cF6XtZ1CWGYT8si1/a0m/tHukE3eVNbtWg4j4DeBHgKsj4iTwC8D9wNGI+DjwMnAnQGY+HRFHgWeAReDuzOw21Hdpag0CvDdq5l5vVz2ZifAKVRW1Zrhn5sfGvHTzmPb3AfdtplPSdjfYOGx4jxkYrrn3jzudcOauorxCVWrAIMCHZ/D911hxPBOGu8oy3KUGDGrp3RG7Qo6axc90wqKMijLcpQYMQn1Qnhm3FHJw3Am3H1BZhrvUgG53eTlm3EVMg20JZmf8o6iy/I2SGnB+5r7GRUxLM3dr7irLcJcaMPxB6rg7Lg3CfbYTK+67Km2G4S41YPgipvoHp6O2H5idMdxVluEuNaA7tP3AILf74b7UbvnMfaJdVMsZ7lIDekO19kF5Zrj80js/c++s2ANe2gzDXWrA8My9V5u516sv3TGhL22W4S41oDtUax9fc+8fz810VtzgQ9oMw11qwMqaexLRX/I4uiwTK27wIW2G4S41oLqGaVlZJhhcibrUbhD0c52OH6iqKMNdasDwVr9J0okgYvRqmZlOWJZRUYa71IBRH6h2IlbsIXN+KaRlGRVmuEsNWNo4bOkD1ZE192pvGT9QVWmGu9SAQVAv3axjMHP3IiZNhuEuNWDFDbJ7SSf6r43afmDOi5hUmOEuNWB56SWXau4dRm8c5t4yKsxwlxpQn4V3M8fX3JftOWO4qxzDXWrAsnDvJZlJp7Oy5t6trXO3LKOSDHepAcM35BhcxBQxXHOvl2Um3Uu1meEuNWDFzL26iGn4jkuDVTVz1W32XA6pUgx3qQHD90nt9iCqi5hG7i1TLaXxQiaVYrhLDVgcXJ1UHXd7PeZmgmD4A9VBWaaz7Lm0WYa71ICFxeWrZRa6yexMEEMbh+XQzL32d4K0KYa71IDFXm/poqUeLHR7zHU641fLVDN3yzIqxXCXGnBuscf83AzQD+zFaube6YzfOKz+XNosw11qwEI3z4d7r5cs9nrMnp+5r9x+YKksY7irDMNdasBCt8f8bFVq6fVr7nOznZX7ufeWf6DqhUwqxXCXCstMFnvJznpZptdjrhMEo5dCzp0vy0y6t2orw10qbKG6x94lVbgvdpOFxarmPrRaZvAB6kzHmrvKMtylwha6/fWM75ifBeDsYpeFXo+5mQ6zM53zrw9e2zHbYSaqi5icuqsQw10qbLGauV9ehft3z3X7q2U6wfzcDGcWa+G+0K/Nd5y5qzDDXSrsXDUzf+f8HABnFnv9de4zHS6Z63DmXPd827OLXXbOzZxfLTP4i0HarNmt7oDUNgvnw73/x+vMQpfFXlZlmeDM4lK4n1noMT/X4dId/fr8dxe6K7+htAHO3KXCRoZ7t8fsTDA/O8OZhXq4d5mfneHSHf223zm3OPkOq5UMd6mwheGyzEK3v7dMp8P8XIczC/UPVPtXsg5m7m+fdeauMgx3qbA3v7MAwO7LdwL90svZxR47ZoP5HTPLSi9nFrrsnO04c1dxjYV7RNwaEc9FxImIuKepnyNNmzfeOgfAnisuBeCts4v86dtnufodO5mfneHcYu/8lalnFrrMz81w2U5n7iqrkXCPiBngPwIfAW4APhYRNzTxs6Rp88ZbZwHYe9UlvHN+lqdeeZNewu7L58/X4f/sO/2/AP78uwtcumOp5v7WWWfuKqOp1TIHgROZ+X8BIuJB4HbgmZI/5Bvf/DY/81+/ckHvyQ2sI97Q4rQLfNNGfsakxrKRpdd5gT9pQz9jQqsGL/S/86t/fobLdszwrst28le+95184ZnXAdhz5SVcWl21+uFPf5Hdl8/zwhtvc/uNe7jqsh2867Id/MLDT/O5L79UfAyaXj/y3l186rbyc9+mwn0P8Ert+Ungh+oNIuIu4C6Ad7/73Rv6IfOzM7x39zsv/I0xkbcQcWHv2tjP2MB7NvRzLvxdF/yODY1lA/2awH+z2/7ateyY7fDPb7uBI196kb1XXsqH//LVRASfuHk/J069BcAPXX8VP/nB9zDTCe77u9/P7379m3S9Y8dFZffl841839jI7G/NbxpxJ3BLZv5U9fwngIOZ+bOj2h84cCCPHz9evB+S1GYR8URmHhj1WlMfqJ4E9taeXwe82tDPkiQNaSrc/w+wPyKuj4gdwCHg4YZ+liRpSCM198xcjIifAf4nMAP8SmY+3cTPkiSt1NjeMpn5e8DvNfX9JUnjeYWqJLWQ4S5JLWS4S1ILGe6S1EKNXMR0wZ2IOA1s5prrq4E3CnVnO7jYxguO+WLhmC/MezJz16gXpiLcNysijo+7SquNLrbxgmO+WDjmcizLSFILGe6S1EJtCfcHtroDE3axjRcc88XCMRfSipq7JGm5tszcJUk1hrsktdDUhntEXBERvxUR34iIZyPigxFxVUQ8EhHPV49X1trfW92M+7mIuKV2/gcj4uvVa/8+NnJLoQkYM97PVM+/FhH/LSKuqLXf1uOF0WOuvfaPIyIj4uraudaOOSJ+thrX0xHx6Vr7Vo45Im6MiC9HxFcj4nhEHKy139Zjjoj3VuMafH07Ij458fzKzKn8Ao4AP1Ud7wCuAD4N3FOduwf4xer4BuApYCdwPfAnwEz12uPAB+nfKe1/AB/Z6rFdwHh/FJitzv1im8Y7bszV8V7620W/BFzd9jEDfxv4X8DO6vw1F8GYvzDoM/BjwKNtGnNt7DPAN4H3TDq/pnLmHhGXA38L+GWAzDyXmW/Sv8n2karZEeCO6vh24MHMPJuZLwAngIMRcS1weWZ+Kfv/pX6t9p6pMW68mfmFzFysmn2Z/h2tYJuPF1b9fwzwb4GfZ/n9vNs85p8G7s/Ms9X5U9Vb2jzmBC6vmn0PS3dq2/ZjHnIz8CeZ+RITzq+pDHfgLwGngV+NiK9ExC9FxGXA7sx8DaB6vKZqP+qG3Huqr5Mjzk+bceOt+4f0/+aG7T9eGDPmiPgo8P8y86mh9q0dM/B9wIcj4rGI+IOI+BtV+zaP+ZPAZyLiFeBfA/dW7dsw5rpDwG9UxxPNr2kN91ngB4D/lJnvB96m/8+YcUbVoXKV89Nm1fFGxKeAReDXB6dGfI/tNF4YPeZ/CXwK+Bcj2rd1zPdU568EPgD8E+BoVVtt85h/Gvi5zNwL/BzVzJ52jBmA6N9i9KPAb67VdMS5TY95WsP9JHAyMx+rnv8W/V+Q16t/qlA9nqq1H3VD7pMslTLq56fNuPESEYeBHwf+XvVPs0H77TxeGD/m64GnIuJF+v1/MiK+l3aP+STw+ex7HOjR30yqzWM+DHy+OvebwMFa++0+5oGPAE9m5uvV84nm11SGe2Z+E3glIt5bnboZeIb+TbYPV+cOAw9Vxw8DhyJiZ0RcD+wHHq/+6fMXEfGBaib0k7X3TI1x442IW4F/Cnw0M79Te8u2Hi+MHfOTmXlNZu7LzH30f7l/oGrb1jE/A/x34CaAiPg++h86vkG7x/wq8MPVuZuA56vjbT/mmo+xVJKBSefXVn2KvNYXcCNwHPga/V/+K4F3Acfo/yIcA66qtf8U/U+Zn6P2iTJwAPjj6rX/QHVV7rR9jRnvCfq1uK9WX/+5LeMdN+ah11+kWi3T5jHTD/PPVWN4ErjpIhjz3wSeoL9K5DHgB1s25kuBbwHfUzs30fxy+wFJaqGpLMtIkjbHcJekFjLcJamFDHdJaiHDXZJayHCXpBYy3CWphf4/nv0q3RdIUoUAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# repeat these from above to make it convenient for you to play with changing R\n", "R=2000\n", "lam=6563\n", "sig=get_sigma(R,lam)\n", "\n", "# list of lines to add, with associated amplitudes\n", "lines=[6562.8,6548.3,6583.41] # add wavelengths of desired lines\n", "amps=[1000,500,1500] # add amplitudes of desired lines (must be same number as lines)\n", "\n", "# reset the flux in case you're repeating this cell\n", "flux=np.zeros(len(waves))\n", "\n", "# loop over the lines and add i the flux\n", "for line,amp in zip(lines,amps) :\n", " flux+=amp*gauss(waves,line,sig) #use your gaussian routine to add in the flux from this line\n", " \n", "#plot your spectrum\n", "plt.plot(waves,flux)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For the emission spectrum, it was most efficient to just add the flux from each line by adding single Gaussian at the location of each emission line. But, in a more general case, you might want to simulate something that has a continuous flux distribution. \n", "

\n", "To work with this situation, consider an input simulated stellar spectrum that you want to smooth with your instrumental resolution. You can read an input spectrum, from the file linked at star.txt , which you will need to download (and depending where you download it, may need to add a directory to the file location below). This has two columns, one with wavelength and one with flux. The columns are labelled, and I like to read this type of file using astropy ascii.read , which will return an astropy Table with columns 'WAVE' and 'FLUX', as given in the input file, but if you are familiar with reading a file another way, e.g. with pandas.read_csv(), you can do that." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from astropy.io import ascii\n", "a=ascii.read('star.txt')\n", "\n", "plt.plot(a['WAVE'],a['FLUX'])\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Before we can proceed to smooth the spectrum with an instrumental resolution, we need to see what the dispersion of the input spectrum is. We'll do that by calculating the difference in wavelength between adjacent pixels, and check to see if it is constant across the spectrum." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "dispersion=a['WAVE'][1:]-a['WAVE'][0:-1]\n", "print('minimum dispersion: {:f} maximum dispersion: {:f}'.format(dispersion.min(),dispersion.max()))\n", "disp=dispersion.mean()\n", "print('using dispersion: {:f}'.format(disp))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You now want to smear this high resolution input spectrum with the instrumental dispersion by convolving it with the instrumental profile. This is done in the pixel space of the spectrum, so you need to use a kernel which has the correct width in pixels. Most low resolution spectrographs actually have a constant $\\delta\\lambda$, rather than a constant $R$. A typical low resolution spectrum might have $\\delta\\lambda = 2.5$ Angstroms ( $R = 2000$ at $\\lambda=5000$).
\n", "Start by creating a smoothing kernel given your choice of instrumental resolution." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "dlambda = # delta(lambda) in Angtroms\n", "sig = # corresponding sigma in Angstyroms\n", "sig_pixels = # corresponding sigma in pixels, given disp\n", "\n", "# create a gaussian kernel with this sigma. Let this extend from -5*sig to 5*sig\n", "x=np.arange(-5*sig_pixels,5*sig_pixels)\n", "kernel= #use your gauss routine to fill the kernel array\n", "\n", "# plot it\n", "plt.plot(kernel)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "OK, now we are ready to smooth the high resolution spectrum with your kernel. You could loop through each pixel and add the Gaussian corresponding to it to the output image. This is the process of convolution. However, you need to write this routine carefully or it can be quite inefficient.\n", "

\n", "Alternatively, you can use numpy.convolve() or scipy.signal.convolve(). Note that these have a keyword mode=, which can take a value of 'full', 'valid, or 'same'. You can play with all three and should be able to explain the difference.\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# if you want, try to code your own convolution here!\n", "def convolve(spec,kernel):\n", " return\n", " \n", "# or just use np.convolve()\n", "print(a['FLUX'].shape)\n", "c=np.convolve(a['FLUX'],kernel,mode='full')\n", "print(c.shape)\n", "c=np.convolve(a['FLUX'],kernel,mode='valid')\n", "print(c.shape)\n", "c=np.convolve(a['FLUX'],kernel,mode='same')\n", "print(c.shape)\n", "plt.plot(a['WAVE'],a['FLUX'])\n", "\n", "plt.plot(a['WAVE'],c)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

Image simulation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "OK, now let's turn to an imaging simulation, i.e. a 2D simulation of a star field. Start by setting up a blank square array, and define a function to return positions and fluxes for a specified number of artificial stars." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "\n", "def getstars(npix,nstars) :\n", " \"\"\" return positions beween 0 and npix and fluxes for artificial stars\n", " \"\"\"\n", " x= #use random routine to choose pixel positions on the array for nstars \n", " y=\n", " # if you wanted to be weird, you could use one of your random number generators like the triangle routine \n", " # to distribute stars more randomly. \n", "\n", " #choose fluxes\n", " flux= #choose fluxes as you desire\n", " # If you wanted to be more realistic, you could use your Salpeter routine to choose fluxes based on masses!\n", "\n", " return x,y,flux\n", " \n", "# set size of output array and create blank 2D array\n", "n= # specify number of pixels\n", "im=np.zeros([n,n])\n", "\n", "nstars= # set number of stars\n", "x,y,flux=getstars(n,nstars)\n", "\n", "# show the locations of your stars\n", "plt.scatter(x,y)\n", "plt.xlim(0,n)\n", "plt.ylim(0,n)\n", "plt.gca().set_aspect('equal')\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For this simulation, we will assume an aximuthally symmetric smoothing function (the point spread function) that takes a Gaussian form. Define a routine to return the value of a 2D gaussian at distance r. Note that a 2D symmetric normalized gaussian is:\n", "$$ G(r,sig) = {1\\over 2\\pi\\sigma^2} exp{-{r^2\\over 2\\sigma^2}}$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def gauss2d(r,sig=1) :\n", " \"\"\" for input 2D array of distances, return Gaussian at those distances\n", " \"\"\"\n", " return # return value of normalized 2D Gaussian\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Using the sampling theorem, we want to choose pixel size such that there are at least 2 pixels per FWHM. Given the diagonal of a pixel is 1.41 times the pixel size, let's go with a smoothing function of FWHM=3. As above, what gaussian sigma does this correspond to?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fwhm=3\n", "sig= # fill in corresponding Gaussian sigma to this fwhm" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Write a routine to add stars to an input, given input image, list of (x,y) positions, and Gaussian sigma" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def mkstars(im,x,y,sig) :\n", " \"\"\" Add stars to input im array, given x, y positions and Gaussian sigma\n", " \"\"\"\n", " nstars=len(x)\n", " ny,nx = im.shape\n", " for i in range(nstars) :\n", " # add pixels only within 5 sigma, but make sure we don't go off edge of array\n", " xs=np.max([0,int(np.rint(x[i]-5*sig))])\n", " ys=np.max([0,int(np.rint(y[i]-5*sig))])\n", " xe=np.min([nx-1,int(np.rint(x[i]+5*sig))])\n", " ye=np.min([ny-1,int(np.rint(y[i]+5*sig))])\n", " # loop over destination pixels\n", " for iy in range(ys,ye) :\n", " for ix in range(xs,xe) :\n", " r2 = (iy-y[i])**2 + (ix-x[i])**2\n", " im[iy,ix] += # add in flux using your flux array and 2D gaussian array \n", "\n", "#repeat setup from above so you can play with things in a single cell\n", "# set size of output array and create\n", "n= \n", "im=np.zeros([n,n])\n", "\n", "# choose positions for artificial stars and place them at random locations on the array\n", "nstars=\n", "x,y,flux=getstars(n,nstars)\n", "\n", "# note the magic time command, which is very useful\n", "%time mkstars(im,x,y,sig)\n", "\n", "plt.imshow(im)\n", "plt.colorbar()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Adjust some of your parameters to make the task computationally hard, i.e., to take a couple of seconds or so, by adding more stars, or using a larger fwhm." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To demonstrate how coding things can make a difference, let's try a different routine for adding the Gaussians, avoiding for loops in favor of numpy array arithmetic. Here, we will pass the gaussian routine a 2D array with distances from an object (actually, the square of distances, so we avoid taking a square root and then squaring it), and add in the Gaussians based on this. \n", "

\n", "I provide all of the code, but you should read it line by line to understand what it is doing, and ask if you have quesitons!" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# alternative 2D gaussian array\n", "def gauss2(r2,sig=1,siglimit=5):\n", " \"\"\" Give input 2D array of dist**2, return array with a normalized Gaussian at these position\n", " \"\"\"\n", " # only fill in pixels within a distance of siglimit*sig \n", " gd=np.where(r2<(siglimit*sig)**2)\n", " # output array at full size\n", " out=np.zeros(r2.shape)\n", " # populate output with 2D gaussian\n", " out[gd[0],gd[1]]=1/(2*np.pi*sig**2)*np.exp(-r2[gd[0],gd[1]]/2/sig**2)\n", " return out\n", "\n", "\n", "# alternative mkstar routine\n", "def mkstars2(im,x,y,sig) :\n", " \"\"\" Add stars to input im array, given x, y positions and Gaussian sigma\n", " \"\"\"\n", " nstars=len(x)\n", " ny,nx = im.shape\n", " # get 2D arrays that contain the x and y pixel values using np.mgrid()\n", " ypix,xpix = np.mgrid[0:ny,0:nx]\n", " \n", " # loop over stars. For each one, calculate the dist**2 array of every pixel to the star center\n", " for i in range(nstars) :\n", " r2=(xpix-x[i])**2+(ypix-y[i])**2\n", " im += flux[i]*gauss2(r2,sig)\n", " \n", "# make a new array\n", "im=np.zeros([n,n])\n", "\n", "# how fast is this?\n", "%time mkstars2(im,x,y,sig)\n", "\n", "plt.imshow(im)\n", "plt.colorbar()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Alternatively, we could create an image with delta functions at the position of each star and convolve the entire image with our smoothing function. To do this, we might want an input image at a higher sampling, since our initial delta function would have to be full pixels, and thus couldn't account for fractional pixel positions. But let's just see how it would go by putting the stars at the nearest pixel.\n", "

\n", "First, make the smoothing kernel. Again, we'll make use of the numpy.mgrid() routine. Ask if there's anything you don't understand." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import scipy\n", "\n", "ypix,xpix = np.mgrid[-5*sig:5*sig,-5*sig:5*sig]\n", "r2=xpix**2+ypix**2\n", "kernel=1/2./np.pi/sig**2*np.exp(-r2/2/sig**2)\n", "plt.imshow(kernel)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# make a new array and add delta functions at the locations of the artificial stars\n", "im=np.zeros([n,n])\n", "for i in range(nstars) :\n", " im[int(y[i]),int(x[i])] += 1\n", "\n", "plt.imshow(im)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, do the convolution. For an nD array, need to use scipy.signal.convolve()." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%time out=scipy.signal.convolve(im,kernel,mode='same')\n", "plt.imshow(out)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "How does the speed compare? \n", "

Let's see what happens even if we increase the resolution before doing the smoothing, to allow for fractional positions of stars: would it still be faster?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# make a new array, oversampled by factor of 10\n", "im=np.zeros([n*10,n*10])\n", "for i in range(nstars) :\n", " im[int(y[i]*10),int(x[i]*10)] += 1\n", "\n", "plt.imshow(im)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# now make the kernel, oversampled by 10\n", "ypix,xpix = np.mgrid[-5*sig*10:5*sig*10,-5*sig*10:5*sig*10]\n", "r2=xpix**2+ypix**2\n", "kernel=1/2./np.pi/(sig*10)**2*np.exp(-r2/2/(sig*10)**2)\n", "plt.imshow(kernel)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# now convolve\n", "%time out=scipy.signal.convolve(im,kernel,mode='same')\n", "plt.imshow(out)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "How does the speed compare?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since we oversampled the simulation to allow for fractional pixel input, now we could downsample it back to a pixel size appropriate for an obervational setup (with just a few pixels per FWHM)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sim=out[::10,::10]\n", "plt.imshow(sim)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "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.11" } }, "nbformat": 4, "nbformat_minor": 2 }