synphot issue 159

synphot_issue_159.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from astropy import units as u\n",
    "\n",
    "import pysynphot as psyn\n",
    "import synphot as syn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Box centered at 1 micron with width of 1 micron\n",
    "w = 1 * u.micron\n",
    "w_AA = w.to(u.AA).value\n",
    "pbox = psyn.Box(w_AA, w_AA)\n",
    "sbox = syn.SpectralElement(syn.models.Box1D, amplitude=1, x_0=w, width=w)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Flat spectrum with 1 Jy across the board\n",
    "pflat = psyn.FlatSpectrum(1, fluxunits='jy')\n",
    "sflat = syn.SourceSpectrum(syn.models.ConstFlux1D, amplitude=1*u.Jy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Build observation from things above\n",
    "pobs = psyn.Observation(pflat, pbox, binset=pbox.wave)\n",
    "sobs = syn.Observation(sflat, sbox, binset=sbox.waveset)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.999999999999999 Jy\n",
      "0.1657983392178991 PHOTLAM\n"
     ]
    }
   ],
   "source": [
    "# Effstim from pysynphot\n",
    "print(pobs.effstim('jy'), 'Jy')\n",
    "presult = pobs.effstim()\n",
    "print(presult, 'PHOTLAM')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0 Jy\n",
      "0.1509190959275364 PHOTLAM\n"
     ]
    }
   ],
   "source": [
    "# Effstim from synphot\n",
    "print(sobs.effstim('Jy'))\n",
    "sresult = sobs.effstim()\n",
    "print(sresult)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But what does this mean? Let's convert PHOTLAM results from each back to Jy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0986129553256427 Jy\n",
      "1.0000201134350415 Jy\n"
     ]
    }
   ],
   "source": [
    "# Unit conversion in pysynphot\n",
    "print(psyn.units.Photlam().ToJy(w_AA, presult), 'Jy')\n",
    "print(psyn.units.Photlam().ToJy(w_AA, sresult.value), 'Jy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0985914081734778 Jy\n",
      "1.0000004999893348 Jy\n"
     ]
    }
   ],
   "source": [
    "# Unit conversion in synphot\n",
    "print(syn.units.convert_flux(w, presult, 'Jy'))\n",
    "print(syn.units.convert_flux(w, sresult, 'Jy'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So... Low-level unit conversion agree between the two (more or less). But somewhere in `pysynphot`'s `effstim`, conversion between per-Hz and per-Angstrom flux is amiss? From the given `effstim` formula, when spectrum is flat and things are set to 1, things should cancel out and produce result of 1 Jy. Even when the flux is first converted to PHOTLAM, the spectrum is still essentially flat in Jy and should give `effstim` that is equivalent of 1 Jy? This is the case for `synphot`, but not `pysynphot`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# What if we make another observation that is in PHOTLAM but still flat over Jy?\n",
    "pflat2 = psyn.ArraySpectrum(wave=pbox.wave, flux=pflat(pbox.wave))\n",
    "sflat2 = syn.SourceSpectrum(\n",
    "    syn.models.Empirical1D, points=sbox.waveset,\n",
    "    lookup_table=sflat(sbox.waveset, flux_unit='photlam'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "pobs2 = psyn.Observation(pflat2, pbox, binset=pbox.wave)\n",
    "sobs2 = syn.Observation(sflat2, sbox, binset=sbox.waveset)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0000000000000002\n",
      "0.1657983392178993 PHOTLAM (was 0.1657983392178991 above)\n"
     ]
    }
   ],
   "source": [
    "# Effstim from pysynphot\n",
    "print(pobs2.effstim('jy'))\n",
    "print(pobs2.effstim(), 'PHOTLAM', '(was {} above)'.format(presult))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0 Jy\n",
      "0.1509190959275364 PHOTLAM (was 0.1509190959275364 PHOTLAM above)\n"
     ]
    }
   ],
   "source": [
    "# Effstim from synphot\n",
    "print(sobs.effstim('Jy'))\n",
    "print(sobs.effstim(), '(was {} above)'.format(sresult))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This shows that `pysynphot` and `synphot` each is consistent internally, but not with each other."
   ]
  }
 ],
 "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}

Don't leave comments here. I won't get notified. It is GitHub bug.