diff --git a/docs/examples/EmLine_examples.ipynb b/docs/examples/EmLine_examples.ipynb index 0305bdc4..83c8c9a8 100644 --- a/docs/examples/EmLine_examples.ipynb +++ b/docs/examples/EmLine_examples.ipynb @@ -315,21 +315,21 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 2", + "display_name": "Python 3", "language": "python", - "name": "python2" + "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", - "version": 2 + "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", - "pygments_lexer": "ipython2", - "version": "2.7.11" + "pygments_lexer": "ipython3", + "version": "3.6.1" } }, "nbformat": 4, diff --git a/docs/examples/EmSystem_examples.ipynb b/docs/examples/EmSystem_examples.ipynb new file mode 100644 index 00000000..63a6b79f --- /dev/null +++ b/docs/examples/EmSystem_examples.ipynb @@ -0,0 +1,384 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Examples with EmSystem" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# imports\n", + "from imp import reload \n", + "from linetools.isgm import emsystem as liems\n", + "from linetools import io as lio\n", + "from linetools.analysis import emline as lta_eml\n", + "import linetools" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## With ALIS" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "alisb_file = linetools.__path__[0]+'/isgm/tests/files/spec1d_J0018p2345_KASTb_coadd.mod.out'\n", + "alisr_file = linetools.__path__[0]+'/isgm/tests/files/spec1d_J0018p2345_KASTr_coadd.mod.out'" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "linetools.lists.parse: Reading linelist --- \n", + " /home/xavier/local/Python/linetools/linetools/data/lines/galaxy_forbidden.ascii\n", + "linetools.lists.parse: Reading linelist --- \n", + " /home/xavier/local/Python/linetools/linetools/data/lines/galaxy_recomb.ascii\n", + "linetools.lists.parse: Reading linelist --- \n", + " /home/xavier/local/Python/linetools/linetools/data/lines/galaxy_abs.ascii\n", + "read_sets: Using set file -- \n", + " /home/xavier/local/Python/linetools/linetools/lists/sets/llist_v1.2.ascii\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/xavier/.pyenv/versions/anaconda3-4.3.1/lib/python3.6/site-packages/astropy/table/column.py:1096: MaskedArrayFutureWarning: setting an item on a masked array which has a shared mask will not copy the mask and also change the original mask array in the future.\n", + "Check the NumPy 1.11 release notes for more information.\n", + " ma.MaskedArray.__setitem__(self, index, value)\n", + "/home/xavier/local/Python/linetools/linetools/lists/linelist.py:376: UserWarning: Not implemented: will not set relative strength for LineList: Galaxy.\n", + " warnings.warn('Not implemented: will not set relative strength for LineList: {}.'.format(self.list))\n" + ] + }, + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "reload(liems)\n", + "emsys = liems.EmSystem.from_alis(alisb_file, 'J001859.32+234540.32')\n", + "emsys" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "0.01540425" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "emsys.zem" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "[,\n", + " ,\n", + " ,\n", + " ,\n", + " ,\n", + " ,\n", + " ,\n", + " ]" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "emsys._emlines" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "{'EW': ,\n", + " 'coord': ,\n", + " 'flag_EW': 0,\n", + " 'flag_flux': 1,\n", + " 'flux': ,\n", + " 'sig_EW': ,\n", + " 'sig_flux': }" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "emsys._emlines[0].attrib" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### More lines" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/xavier/local/Python/linetools/linetools/lists/linelist.py:376: UserWarning: Not implemented: will not set relative strength for LineList: Galaxy.\n", + " warnings.warn('Not implemented: will not set relative strength for LineList: {}.'.format(self.list))\n" + ] + }, + { + "data": { + "text/plain": [ + "[,\n", + " ,\n", + " ,\n", + " ,\n", + " ,\n", + " ,\n", + " ,\n", + " ,\n", + " ,\n", + " ,\n", + " ,\n", + " ,\n", + " ]" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "emsys.add_emlines_from_alis(alisr_file, chk_z=False)\n", + "emsys._emlines" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Utilities" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "Hbeta = emsys.get_emline('Hbeta')\n", + "Hbeta" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Metallicity" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "(7.160611600976936, 7.159297218972787)" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "reload(lta_eml)\n", + "r_val, s_val = lta_eml.metallicity('PG16', emsys)\n", + "r_val, s_val" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "collapsed": true + }, + "source": [ + "## Write to disk" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Wrote J001859.32+234540.32_z0.015 system to J001859.32+234540.32_emsys.json file\n" + ] + } + ], + "source": [ + "tdict = emsys.write_json('J001859.32+234540.32_emsys.json')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Read from disk" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/xavier/local/Python/linetools/linetools/lists/linelist.py:376: UserWarning: Not implemented: will not set relative strength for LineList: Galaxy.\n", + " warnings.warn('Not implemented: will not set relative strength for LineList: {}.'.format(self.list))\n", + "/home/xavier/local/Python/linetools/linetools/analysis/linelimits.py:165: UserWarning: Redshift=0. If this is unexpected, set _z and reset limits\n", + " warnings.warn(\"Redshift=0. If this is unexpected, set _z and reset limits\")\n" + ] + }, + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "reload(liems)\n", + "liems.EmSystem.from_json('J001859.32+234540.32_emsys.json')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "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.6.1" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/linetools/analysis/emline.py b/linetools/analysis/emline.py new file mode 100644 index 00000000..5929c354 --- /dev/null +++ b/linetools/analysis/emline.py @@ -0,0 +1,63 @@ +""" Utlities for the analysis of emission lines +""" +from __future__ import print_function, absolute_import, division, unicode_literals + +import numpy as np +import pdb +import warnings + +from astropy import units as u +from astropy import constants as const +from astropy.io import ascii +from astropy.utils import isiterable +from linetools.lists.linelist import LineList + +# Atomic constant +atom_cst = (const.m_e.cgs*const.c.cgs / (np.pi * + (const.e.esu**2).cgs)).to(u.AA*u.s/(u.km*u.cm**2)) + +# e2/(me*c) in CGS +e2_me_c_cgs = (const.e.esu**2 / (const.c.to('cm/s') * const.m_e.to('g'))) + + +# Estimate metallicity from a tried-and-true method +def metallicity(method, emsystem): + """ Calculate a nebular abundance from an input set of Emission lines + + Parameters + ---------- + method : str + Name of the method to be applied + PG16 -- http://adsabs.harvard.edu/cgi-bin/bib_query?arXiv:1601.08217 + Requires Hbeta, [OII], [OIII], [NII], [SII] + Return r_val and s_val + emsystem : EmSystem + + Returns + ------- + """ + if method == 'PG16': + # Requires Hbeta, [OII], [OIII], [NII], [SII] + R2 = (emsystem.get_emline('[OII] 3726').attrib['flux'] + + emsystem.get_emline('[OII] 3729').attrib['flux']) / emsystem.get_emline('Hbeta').attrib['flux'] + R3 = (emsystem.get_emline('[OIII] 4959').attrib['flux'] + + emsystem.get_emline('[OIII] 5007').attrib['flux']) / emsystem.get_emline('Hbeta').attrib['flux'] + N2 = (emsystem.get_emline('[NII] 6548').attrib['flux'] + + emsystem.get_emline('[NII] 6584').attrib['flux']) / emsystem.get_emline('Hbeta').attrib['flux'] + S2 = (emsystem.get_emline('[SII] 6716').attrib['flux'] + + emsystem.get_emline('[SII] 6731').attrib['flux']) / emsystem.get_emline('Hbeta').attrib['flux'] + # Proceed + if np.log10(N2) < -0.6: + r_val = 7.932 + 0.944*np.log10(R3/R2) + 0.695*np.log10(N2) + \ + ((0.97 - 0.291*np.log10(R3/R2)) - 0.019*np.log10(N2))*np.log10(R2) + + s_val = 8.072 + 0.789*np.log10(R3/S2) + 0.726*np.log10(N2) + \ + (1.069 - 0.170*np.log10(R3/S2) +0.022*np.log10(N2))*np.log10(S2) + else: + r_val = 8.589 + 0.022*np.log10(R3/R2) + 0.399*np.log10(N2) + \ + (-0.137 + 0.164*np.log10(R3/R2) + 0.589*np.log10(N2))*np.log10(R2) + + s_val = 8.424 + 0.030*np.log10(R3/S2) + 0.751*np.log10(N2) + \ + (-0.349 + 0.182*np.log10(R3/S2) +0.508*np.log10(N2))*np.log10(S2) + return r_val.decompose().value, s_val.decompose().value + diff --git a/linetools/analysis/setup_package.py b/linetools/analysis/setup_package.py new file mode 100644 index 00000000..fe716ce2 --- /dev/null +++ b/linetools/analysis/setup_package.py @@ -0,0 +1,4 @@ +def get_package_data(): + # Installs the testing data files. Unable to get package_data + # to deal with a directory hierarchy of files, so just explicitly list. + return {'linetools.analysis.tests': ['files/*.fits', 'files/*.all', 'files/*.out']} diff --git a/linetools/analysis/tests/files/spec1d_J0018p2345_KASTb_coadd.mod.out b/linetools/analysis/tests/files/spec1d_J0018p2345_KASTb_coadd.mod.out new file mode 100644 index 00000000..4bfbda3b --- /dev/null +++ b/linetools/analysis/tests/files/spec1d_J0018p2345_KASTb_coadd.mod.out @@ -0,0 +1,140 @@ +# +# Generated by ALIS on 18/10/16 at 15:19:12 +# +# Running Time (hrs) = 0.001082 +# Initial Chi-Squared = 25517.886155 +# Bestfit Chi-Squared = 1077.242298 +# Degrees-of-Freedom = 231 +# Num. of Iterations = 15 +# Convergence Reason = The relative reduction in the sum of squares is less than atol + +run blind False +run nsubpix 20 +chisq ftol 1.0E-10 +chisq atol 0.001 +chisq miniter 10 +chisq maxiter 1000 +out fits False +out plots spec1d_J0018p2345_KASTb_coadd.pdf +plot dims 3x3 +plot labels True +plot fits False + +data read + spec1d_J0018p2345_KASTb_coadd.dat specid=0 fitrange=[4917.0,4957.0] loadrange=[4907.0,4967.0] resolution=vfwhm(282.946vtie) columns=[wave,flux,error] label=Hb4863 + spec1d_J0018p2345_KASTb_coadd.dat specid=1 fitrange=[4388.0,4428.0] loadrange=[4378.0,4438.0] resolution=vfwhm(282.946vtie) columns=[wave,flux,error] label=Hg4342 + spec1d_J0018p2345_KASTb_coadd.dat specid=2 fitrange=[4146.0,4186.0] loadrange=[4136.0,4196.0] resolution=vfwhm(282.946vtie) columns=[wave,flux,error] label=Hd4103 + spec1d_J0018p2345_KASTb_coadd.dat specid=3 fitrange=[4012.0,4052.0] loadrange=[4002.0,4062.0] resolution=vfwhm(282.946vtie) columns=[wave,flux,error] label=He3970 + spec1d_J0018p2345_KASTb_coadd.dat specid=4 fitrange=[3772.0,3802.0] loadrange=[3767.0,3807.0] resolution=vfwhm(282.946vtie) columns=[wave,flux,error] label=[OII]3727+3729 + spec1d_J0018p2345_KASTb_coadd.dat specid=5 fitrange=[5021.0,5051.0] loadrange=[5016.0,5056.0] resolution=vfwhm(282.946vtie) columns=[wave,flux,error] label=[OIII]4960 + spec1d_J0018p2345_KASTb_coadd.dat specid=6 fitrange=[5065.0,5105.0] loadrange=[5060.0,5110.0] resolution=vfwhm(282.946vtie) columns=[wave,flux,error] label=[OIII]5008 +data end + +model read +fix vfwhm value False +fix gaussian dispersion True + emission + legendre 136.86663193 -8.85011622 scale=1.0,1.0 specid=0 continuum=True + legendre 195.06532795 11.12511048 scale=1.0,1.0 specid=1 continuum=True + legendre 205.78720694 -7.67128769 scale=1.0,1.0 specid=2 continuum=True + legendre 238.51636060 4.99858582 scale=1.0,1.0 specid=3 continuum=True + legendre 254.20178241 -61.15994519 scale=1.0,1.0 specid=4 continuum=True + legendre 154.08236287 31.95230380 scale=1.0,1.0 specid=5 continuum=True + legendre 152.53327300 35.17673753 scale=1.0,1.0 specid=6 continuum=True + gaussian 7.85749594E+03 0.01540425ra 30.00000000dh specid=0 IntFlux=True wave=4862.721 + gaussian 3.59948935E+03 0.01540425ra 30.00000000dh specid=1 IntFlux=True wave=4341.684 + gaussian 1.65562745E+03 0.01540425ra 30.00000000dh specid=2 IntFlux=True wave=4102.891 + gaussian 1.56279300E+03 0.01540425ra 30.00000000dh specid=3 IntFlux=True wave=3971.195 + gaussian 3.46048372E+03 0.01540425ra 30.00000000do specid=4 IntFlux=True wave=3727.092 + gaussian 4.88103542E+03 0.01540425ra 30.00000000do specid=4 IntFlux=True wave=3729.874 + gaussian 4.79330012E+03orata 0.01540425ra 30.00000000dox specid=5 IntFlux=True wave=4960.295 + gaussian 1.43799004E+04oratb 0.01540425ra 30.00000000dox specid=6 IntFlux=True wave=5008.239 +model end + +link read +oratb(orata) = orata*3.0 +link end + +# +# Errors: +# +#emission +# legendre 4.83938676 15.35821602 scale=1.0,1.0 specid=0 continuum=True +# legendre 7.08376397 21.93070237 scale=1.0,1.0 specid=1 continuum=True +# legendre 7.53275066 22.28386888 scale=1.0,1.0 specid=2 continuum=True +# legendre 10.32401860 33.21308807 scale=1.0,1.0 specid=3 continuum=True +# legendre 14.91070787 41.98922735 scale=1.0,1.0 specid=4 continuum=True +# legendre 5.34940707 16.10995876 scale=1.0,1.0 specid=5 continuum=True +# legendre 4.65726318 13.63971131 scale=1.0,1.0 specid=6 continuum=True +# gaussian 1.07203070E+02 0.00000346ra 0.00000000dh specid=0 IntFlux=True wave=4862.721 +# gaussian 1.24076227E+02 0.00000346ra 0.00000000dh specid=1 IntFlux=True wave=4341.684 +# gaussian 1.34774396E+02 0.00000346ra 0.00000000dh specid=2 IntFlux=True wave=4102.891 +# gaussian 1.51170291E+02 0.00000346ra 0.00000000dh specid=3 IntFlux=True wave=3971.195 +# gaussian 2.13668992E+02 0.00000346ra 0.00000000do specid=4 IntFlux=True wave=3727.092 +# gaussian 2.14767652E+02 0.00000346ra 0.00000000do specid=4 IntFlux=True wave=3729.874 +# gaussian 9.65686646E+01orata 0.00000346ra 0.00000000dox specid=5 IntFlux=True wave=4960.295 +# gaussian 0.00000000E+00oratb 0.00000346ra 0.00000000dox specid=6 IntFlux=True wave=5008.239 + +# Convolution Models: +# vfwhm 282.946vtie +# +# Errors: +# +# vfwhm 2.315vtie + +# Shift Models: +# Ashift 0.0000000 +# +# Errors: +# +# Ashift 0.0000000 + + +################################################### +# # +# HERE IS A COPY OF THE INPUT MODEL # +# # +################################################### +# +# run blind False +# run nsubpix 20 +# chisq ftol 1.0E-10 +# chisq atol 0.001 +# chisq miniter 10 +# chisq maxiter 1000 +# out fits False +# out plots spec1d_J0018p2345_KASTb_coadd.pdf +# plot dims 3x3 +# plot labels True +# plot fits False +# data read +# spec1d_J0018p2345_KASTb_coadd.dat specid=0 fitrange=[4917.0,4957.0] loadrange=[4907.0,4967.0] resolution=vfwhm(150.0vtie) columns=[wave,flux,error] label=Hb4863 +# spec1d_J0018p2345_KASTb_coadd.dat specid=1 fitrange=[4388.0,4428.0] loadrange=[4378.0,4438.0] resolution=vfwhm(150.0vtie) columns=[wave,flux,error] label=Hg4342 +# spec1d_J0018p2345_KASTb_coadd.dat specid=2 fitrange=[4146.0,4186.0] loadrange=[4136.0,4196.0] resolution=vfwhm(150.0vtie) columns=[wave,flux,error] label=Hd4103 +# spec1d_J0018p2345_KASTb_coadd.dat specid=3 fitrange=[4012.0,4052.0] loadrange=[4002.0,4062.0] resolution=vfwhm(150.0vtie) columns=[wave,flux,error] label=He3970 +# spec1d_J0018p2345_KASTb_coadd.dat specid=4 fitrange=[3772.0,3802.0] loadrange=[3767.0,3807.0] resolution=vfwhm(150.0vtie) columns=[wave,flux,error] label=[OII]3727+3729 +# spec1d_J0018p2345_KASTb_coadd.dat specid=5 fitrange=[5021.0,5051.0] loadrange=[5016.0,5056.0] resolution=vfwhm(150.0vtie) columns=[wave,flux,error] label=[OIII]4960 +# spec1d_J0018p2345_KASTb_coadd.dat specid=6 fitrange=[5065.0,5105.0] loadrange=[5060.0,5110.0] resolution=vfwhm(150.0vtie) columns=[wave,flux,error] label=[OIII]5008 +# data end +# model read +# fix vfwhm value False +# fix gaussian dispersion True +# emission +# legendre 146.7779 0.01 scale=[1.0,1.0] specid=0 continuum=True +# legendre 222.0621 0.01 scale=[1.0,1.0] specid=1 continuum=True +# legendre 220.9873 0.01 scale=[1.0,1.0] specid=2 continuum=True +# legendre 278.4923 0.01 scale=[1.0,1.0] specid=3 continuum=True +# legendre 289.6434 0.01 scale=[1.0,1.0] specid=4 continuum=True +# legendre 169.9477 0.01 scale=[1.0,1.0] specid=5 continuum=True +# legendre 161.8335 0.01 scale=[1.0,1.0] specid=6 continuum=True +# gaussian 6547.75804488 0.01541ra 30.0dh wave=4862.721 specid=0 IntFlux=True +# gaussian 4038.49307181 0.01541ra 30.0dh wave=4341.684 specid=1 IntFlux=True +# gaussian 2312.54298859 0.01541ra 30.0dh wave=4102.891 specid=2 IntFlux=True +# gaussian 2603.38462088 0.01541ra 30.0dh wave=3971.195 specid=3 IntFlux=True +# gaussian 22165.9457491 0.01541ra 30.0do wave=3727.092 specid=4 IntFlux=True +# gaussian 7758.0810122 0.01541ra 30.0do wave=3729.874 specid=4 IntFlux=True +# gaussian 3654.59671255orata 0.01541ra 30.0dox wave=4960.295 specid=5 IntFlux=True +# gaussian 10963.7901376oratb 0.01541ra 30.0dox wave=5008.239 specid=6 IntFlux=True +# model end +# + diff --git a/linetools/analysis/tests/test_io.py b/linetools/analysis/tests/test_io.py new file mode 100644 index 00000000..76799d0f --- /dev/null +++ b/linetools/analysis/tests/test_io.py @@ -0,0 +1,25 @@ +# Tests of linetools.utils + +from __future__ import print_function, absolute_import, division, unicode_literals + +# TEST_UNICODE_LITERALS + +import pytest +import numpy as np +import os + +import linetools.io as lio + + +def data_path(filename): + data_dir = os.path.join(os.path.dirname(__file__), 'files') + return os.path.join(data_dir, filename) + + +def test_emline_from_alis(): + J0018_emlines = lio.emlines_from_alis_output(data_path('spec1d_J0018p2345_KASTb_coadd.mod.out')) + + assert len(J0018_emlines) == 8 + np.testing.assert_allclose(J0018_emlines[0].z, 0.01540425) + np.testing.assert_allclose(J0018_emlines[1].wrest.value, 4341.684) #Hgamma + diff --git a/linetools/io.py b/linetools/io.py new file mode 100644 index 00000000..9500b281 --- /dev/null +++ b/linetools/io.py @@ -0,0 +1,148 @@ +""" I/O routines useful to linetools +""" +from __future__ import print_function, absolute_import, division, unicode_literals + +import pdb + +import numpy as np +from astropy import units as u + +from linetools.spectralline import EmLine +from linetools.lists.linelist import LineList + +def emlines_from_alis_output(alis_file): + """ Generate a list of emission lines from a standard ALIS output file + Currently loads + + Parameters + ---------- + alis_file : str + + Returns + ------- + emlines : list + List of EmLine objects + """ + linelist = LineList('Galaxy') + + def get_value(istr): + eqpos = istr.find('=') + return istr[eqpos+1:] + + def strip_end(istr): + for kk in range(-1, -1*len(istr),-1): + try: + _ = int(istr[kk]) + except ValueError: + pass + else: + # Cut + if kk == -1: + cutstr = istr + else: + cutstr = istr[:kk+1] + return float(cutstr) + return + + flg_model = False + flg_error = False + flg_emiss = False + flg_sigemiss = False + line_dict = dict(continuum={}, lines=[], errors=[]) + with open(alis_file) as file_to_search: + for line in file_to_search: + # model + if 'model read' in line: + flg_model = True + # Errors + if 'Errors:' in line: + flg_error = True + if flg_model: + # End + if 'model end' in line: + flg_model = False + flg_emiss = False + if 'emission' in line: + flg_emiss = True + if flg_error: + # End + if line[0] != '#': + flg_error = False + flg_sigemiss = False + flg_emiss = False + if 'emission' in line: + flg_sigemiss = True + flg_emiss = True + # lines + if flg_emiss: + flg_intflux = True + if flg_sigemiss: + conti_sp = line.split()[1:] + else: + conti_sp = line.split() + if len(conti_sp) == 0: + continue + # Find specid + for istr in conti_sp: + if 'specid' in istr: + specid = int(get_value(istr)) + # Add to dict + if 'continuum=True' in line: # Continuum + pass + elif conti_sp[0] in ['gaussian']: # Line :: Could be others + for istr in conti_sp: + if 'IntFlux=True' in istr: + flg_intflux = True + elif 'wave=' in istr: + wrest = float(get_value(istr)) + # Values + if flg_sigemiss is False: + gauss_fit = {} + if flg_intflux: + gauss_fit['flux'] = strip_end(conti_sp[1]) + else: + gauss_fit['amp'] = strip_end(conti_sp[1]) + gauss_fit['z'] = strip_end(conti_sp[2]) + gauss_fit['sigma'] = strip_end(conti_sp[3]) + gauss_fit['type'] = 'gaussian' + gauss_fit['wrest'] = wrest*u.AA + # Init? + line_dict['lines'].append(gauss_fit) + else: + waves = np.array([idict['wrest'].value for idict in line_dict['lines']]) + mt = np.where(np.abs(waves-wrest) < 1e-3)[0] + if len(mt) != 1: + raise ValueError("Bad wrest") + # Fill in + if flg_intflux: + line_dict['lines'][mt[0]]['sig_flux'] = strip_end(conti_sp[1]) + else: + line_dict['lines'][mt[0]]['sig_amp'] = strip_end(conti_sp[1]) + + # Generate EmLines + emlines = [] + for iline in line_dict['lines']: + obj = EmLine(iline['wrest'], z=iline['z'], linelist=linelist, closest=True) + # Fill + try: + obj.attrib['flux'] = iline['flux']*u.erg/u.s + except KeyError: + # Amplitude + pdb.set_trace() + else: + obj.attrib['sig_flux'] = iline['sig_flux']*u.erg/u.s + obj.attrib['flag_flux'] = 1 # Measured + # Append + emlines.append(obj) + # Return + return emlines + + + + + + + + + + diff --git a/linetools/isgm/emsystem.py b/linetools/isgm/emsystem.py new file mode 100644 index 00000000..674ab0fb --- /dev/null +++ b/linetools/isgm/emsystem.py @@ -0,0 +1,533 @@ +""" Class for emission line systems +""" +from __future__ import print_function, absolute_import, division, unicode_literals + +# Python 2 & 3 compatibility +try: + basestring +except NameError: + basestring = str + +import numpy as np +import pdb +import warnings +from abc import ABCMeta + +from astropy import units as u +from astropy.units import Quantity +from astropy.table import Table +from astropy import constants as const +from astropy.coordinates import SkyCoord + +from linetools.isgm import utils as ltiu +from linetools import utils as ltu +from linetools import line_utils as ltlu +from linetools.spectralline import EmLine +from linetools import io as lio + +# Globals to speed things up +c_mks = const.c.to('km/s').value + +class EmSystem(object): + """ + Class for an emission line system, e.g. Planetary Nebula, Galaxy, SN, quasar + + Parameters + ---------- + radec : tuple or coordinate + RA/Dec of the sightline or astropy.coordinate + zem : float + Emission redshift + vlim : Quantity array (2) + Velocity limits of the system + zem : float, optional + Emission redshift of background source + em_type : str or unicode + Type of emission system, e.g. galaxy, quasar, LAE + name : str, optional + Name for the system + + Attributes + ---------- + em_type : str + coord : SkyCoord + RA/Dec of the sightline + zem : float + Emission redshift + zem : float + Emission redshift of background source + vlim : Quantity array (2) + Velocity limits of the system + ZH : float + Metallicity (log10) + name : str + Name of the system + """ + + __metaclass__ = ABCMeta + + @classmethod + def from_emlines(cls, emlines, vlim=None, **kwargs): + """Instantiate from a list of EmLines + + Parameters + ---------- + emlines : list + List of EmLine objects + vlim : list, optional + Velocity limits for the system + If not set, the first components sets vlim + """ + # Check + #assert ltiu.chk_components(components) + # Instantiate with the first component + init_eml = emlines[0] + # + slf = cls(init_eml.attrib['coord'], init_eml.z, vlim=vlim) + if slf.chk_emline(init_eml): + slf._emlines.append(init_eml) + else: + raise IOError("Bad component input") + # Append with component checking + for emline in emlines[1:]: + slf.add_emline(emline) + # Return + return slf + + @classmethod + def from_alis(cls, alis_file, radec, **kwargs): + """ Load from an ALIS output file + + Parameters + ---------- + alis_file : .mod ALIS output file + radec : SkyCoord or tuple of RA,DEC + kwargs + + Returns + ------- + EmSystem + + """ + emlines = lio.emlines_from_alis_output(alis_file) + # Add coordinates + coord = ltu.radec_to_coord(radec) + for emline in emlines: + emline.attrib['coord'] = coord + # + return cls.from_emlines(emlines, **kwargs) + + @classmethod + def from_dict(cls, idict, use_coord=False, **kwargs): + """ Instantiate from a dict. Usually read from the hard-drive + + Parameters + ---------- + idict : dict + + Returns + ------- + + """ + slf = cls(SkyCoord(ra=idict['RA']*u.deg, dec=idict['DEC']*u.deg), + idict['zem'], vlim=idict['vlim']*u.km/u.s, + name=idict['Name'], em_type=idict['em_type']) + # Emission lines + em_dict = idict['emlines'] + for key in em_dict: + iline = em_dict[key] + obj = EmLine.from_dict(iline) + slf.add_emline(obj) + # Return + return slf + + @classmethod + def from_json(cls, json_file, **kwargs): + """ Load from a JSON file (via from_dict) + + Parameters + ---------- + json_file + kwargs + + Returns + ------- + EmSystem + + """ + idict = ltu.loadjson(json_file) + slf = cls.from_dict(idict, **kwargs) + return slf + + + def __init__(self, radec, zem, vlim=None, em_type=None, name=None): + + self.zem = zem + if vlim is None: + self.vlim = [-300., 300.]*u.km/u.s + else: + self.vlim = vlim + self.coord = ltu.radec_to_coord(radec) + if name is None: + self.name = 'J{:s}{:s}_z{:.3f}'.format( + self.coord.ra.to_string(unit=u.hour,sep='',pad=True), + self.coord.dec.to_string(sep='',pad=True,alwayssign=True), + self.zem) + else: + self.name = name + + # Em type + if em_type is None: + self.em_type = 'NONE' + else: + self.em_type = em_type + + # Components + self._emlines = [] # List of EmLine objects + + # Components + #self._components = [] # List of EmComponent objects + + # Kinematics + self.kin = {} + + # Metallicity + self.ZH = 0. + self.sig_ZH = 0. + + # Abundances and Tables + self._EW = Table() + self._fluxes = None # Needs to be None for fill_ion + #self._trans = Table() + self._abund = Table() + + # Refs (list of references) + self.Refs = [] + + def add_emlines_from_alis(self, alis_file, **kwargs): + ''' + Code assumes the coordinates of the class (i.e. no checking) + Parameters + ---------- + alis_file : str + + + + Returns + ------- + + ''' + dum_sys = EmSystem.from_alis(alis_file, self.coord) + # Add in lines + for emline in dum_sys._emlines: + self.add_emline(emline, **kwargs) + + def add_emline(self, emline, tol=0.2*u.arcsec, + chk_sep=True, chk_z=True, overlap_only=False, + vtoler=1., debug=False, **kwargs): + """Add an EmLine object if it satisfies all of the rules. + + For velocities, we demand that the new line has a velocity + range that is fully encompassed by the EmSystem. + We allow a small tolerance for round-off error + + Should check for duplicates.. + + Parameters + ---------- + emline : EmLine + tol : Angle, optional + Tolerance on matching coordinates + Only used if chk_sep=True + chk_sep : bool, optional + Perform coordinate check (expensive) + chk_z : bool, optional + Perform standard velocity range test + overlap_only : bool, optional + Only require that the components overlap in redshift + vtoler : float, optional + Tolerance for velocity in km/s + + Returns + ------- + test : bool + True if successful + """ + # Coordinates + if chk_sep: + testcoord = bool(self.coord.separation(emline.attrib['coord']) < tol) + else: + testcoord = True + # Now redshift/velocity + if chk_z: + # Will avoid Quantity for speed + sys_vlim_mks = self.vlim.to('km/s').value + dz_toler = (1 + self.zem) * vtoler / c_mks + zlim_sys = self.zem + (1 + self.zem) * (sys_vlim_mks / c_mks) + testz = (emline.z >= (zlim_sys[0]-dz_toler)) & ( + emline.z <= (zlim_sys[1]+dz_toler)) + else: + testz = True + + # Additional checks (specific to EmSystem type) + testcomp = self.chk_emline(emline) + test = testcoord & testcomp & testz + + # Append? + if test: + self._emlines.append(emline) + else: + warnings.warn('Input EmLine with wrest={} does not match EmSystem rules. Not appending'.format(emline.wrest)) + if not testcoord: + warnings.warn('Failed coordinate match') + if not testcomp: + warnings.warn('Failed component check') + if not testz: + warnings.warn('Failed velocity overlap') + # + return test + + def chk_emline(self, component): + """Additional checks on the component + + Parameters + ---------- + component : EmComponent + + """ + return True + + def fill_ionN(self, **kwargs): + """ Fills the ionN Table from the list of components + """ + self._ionN = ltiu.iontable_from_components(self._components, **kwargs) + +# def fill_trans(self, **kwargs): +# """ Fills the ionN Table from the list of emission lines +# """ +# self._trans = ltlu.transtable_from_speclines(self.list_of_abslines()) + + def get_emline(self, inp): + """ Returns an EmLine from the EmSystem + + Parameters + ---------- + inp : str or Quantity + str -- Name of the transition, e.g. 'Halpha' + Quantity -- Rest wavelength of the transition, e.g. 6564.61*u.AA + to 0.01 precision + + Returns + ------- + emline -- EmLine object or list of EmLines + More than one will be returned if this line exists in + multiple components. The returned quantity will then + be a list instead of a single object + """ + # Generate the lines + if isinstance(inp,basestring): + names = np.array([emline.name for emline in self._emlines]) + mt = np.where(names == inp)[0] + elif isinstance(inp,Quantity): + wrest = Quantity([emline.wrest for emline in self._emlines]) + mt = np.where(np.abs(wrest-inp) < 0.01*u.AA)[0] + else: + raise IOError("Bad input to emline") + # Finish + if len(mt) == 0: + warnings.warn("No emline with input={}".format(inp)) + return None + elif len(mt) == 1: + return self._emlines[mt[0]] + else: + return [self._emlines[ii] for ii in mt] + +# def get_component(self, inp): +# """ Returns the component related to the given input +# Parameters +# ---------- +# inp : tuple or AbsLine +# tuple -- (Z,ion) integers +# AbsLine -- actual absorption line object + +# Returns +# ------- +# component +# """ +# if isinstance(inp, tuple): +# # Assume Zion for now, e.g. (26,2) +# tuples = [comp.Zion for comp in self._components] +# try: +# idx = tuples.index(inp) +# except ValueError: +# warnings.warn("Input Zion {} is not in any component".format(inp)) +# return None +# else: +# return self._components[idx] +# elif isinstance(inp, AbsLine): +# return self.get_comp_from_absline(inp) +# else: +# raise IOError("Bad input to get_component method") + +# def get_comp_from_absline(self, aline): +# """ Returns the component that holds the input AbsLine + +# Parameters +# ---------- +# aline : AbsLine + +# Returns +# ------- +# comp -- AbsComponent object that holds this AbsLine +# """ +# # Loop on components +# for comp in self._components: +# # Is the line present? +# try: +# idx = comp._abslines.index(aline) +# except ValueError: +# pass +# else: +# return comp +# # Raise error? +# warnings.warn("Input absorption line is not in any component") +# return None + + def measure_restew(self, spec=None, **kwargs): + """ Measure rest-frame EWs for lines in the EmSystem + Parameters + ---------- + spec : XSpectrum1D, optional + kwargs + + Returns + ------- + + """ + # Grab Lines + abs_lines = self.list_of_abslines() + # Loop + for iline in abs_lines: + # Fill in spec? + if spec is not None: + iline.analy['spec'] = spec + # Measure + iline.measure_restew(**kwargs) + + def measure_aodm(self, spec=None, **kwargs): + """ Measure ADOM columns for the list of lines + Note: Components are *not* updated by default + + Parameters + ---------- + spec : XSpectrum1D, optional + kwargs + + Returns + ------- + + """ + # Grab Lines + abs_lines = self.list_of_abslines() + # Loop + for iline in abs_lines: + # Fill in spec? + if spec is not None: + iline.analy['spec'] = spec + # Measure + iline.measure_aodm(**kwargs) + # + print("You may now wish to update the component column densities with update_component_colm()") + + def update_component_colm(self, **kwargs): + """ Synthesize/update column density measurements for components + + Parameters + ---------- + kwargs + + Returns + ------- + + """ + for comp in self._components: + comp.synthesize_colm(**kwargs) + + def stack_plot(self, pvlim=None, **kwargs): + """Show a stack plot of the system, if spec are loaded + Assumes the data are normalized. + + Parameters + ---------- + pvlim : Quantities, optional + Over-ride system vlim for plotting + """ + from linetools.analysis import plots as ltap + if pvlim is not None: + vlim = pvlim + else: + vlim = self.vlim + ltap.stack_plot(self.list_of_abslines(), vlim=vlim, **kwargs) + + def to_dict(self): + """ Write EmSystem data to a dict that can be written with JSON + """ + import datetime + import getpass + date = str(datetime.date.today().strftime('%Y-%b-%d')) + user = getpass.getuser() + # Generate the dict + outdict = dict(Name=self.name, em_type=self.em_type, + vlim=self.vlim.to('km/s').value, zem=self.zem, + RA=self.coord.ra.value, DEC=self.coord.dec.value, + kin=self.kin, Refs=self.Refs, CreationDate=date, + ZH=self.ZH, sig_ZH=self.sig_ZH, + user=user + ) + outdict['class'] = self.__class__.__name__ + outdict['emlines'] = {} + for iline in self._emlines: + outdict['emlines'][iline.wrest.value] = iline.to_dict() + # Polish + outdict = ltu.jsonify(outdict) + # Return + return outdict + + def write_json(self, outfil=None): + """ Generate a JSON file from the system + + Returns + ------- + + """ + # Generate the dict + odict = self.to_dict() + # Write + if outfil is None: + outfil = self.name+'.json' + ltu.savejson(outfil, odict, overwrite=True, easy_to_read=True) + # Finish + print("Wrote {:s} system to {:s} file".format(self.name, outfil)) + + + def __repr__(self): + txt = '<{:s}: name={:s} type={:s}, {:s} {:s}, z={:g}'.format( + self.__class__.__name__, self.name, self.em_type, + self.coord.ra.to_string(unit=u.hour,sep=':',pad=True), + self.coord.dec.to_string(sep=':',pad=True,alwayssign=True), + self.zem) + # Finish + txt = txt + '>' + return (txt) + + +class GenericEmSystem(EmSystem): + """Class for Generic Emission Line System + """ + def __init__(self, radec, zem, **kwargs): + EmSystem.__init__(self, radec, zem, em_type='Generic', **kwargs) + self.name = 'Foo' + + def print_em_type(self): + """"Return a string representing the type of vehicle this is.""" + return 'Generic' diff --git a/linetools/isgm/setup_package.py b/linetools/isgm/setup_package.py index 983adfe9..d2925a3a 100644 --- a/linetools/isgm/setup_package.py +++ b/linetools/isgm/setup_package.py @@ -1,4 +1,4 @@ def get_package_data(): # Installs the testing data files. Unable to get package_data # to deal with a directory hierarchy of files, so just explicitly list. - return {'linetools.isgm.tests': ['files/*.fits', 'files/*.dat', 'files/*.json', 'files/*.all', 'files/*.fits.gz', 'files/*.VP']} + return {'linetools.isgm.tests': ['files/*.fits', 'files/*.dat', 'files/*.json', 'files/*.all', 'files/*.fits.gz', 'files/*.VP', 'files/*.out']} \ No newline at end of file diff --git a/linetools/isgm/tests/files/spec1d_J0018p2345_KASTb_coadd.mod.out b/linetools/isgm/tests/files/spec1d_J0018p2345_KASTb_coadd.mod.out new file mode 100644 index 00000000..c7d6fcb8 --- /dev/null +++ b/linetools/isgm/tests/files/spec1d_J0018p2345_KASTb_coadd.mod.out @@ -0,0 +1,140 @@ +# +# Generated by ALIS on 18/10/16 at 15:19:12 +# +# Running Time (hrs) = 0.001082 +# Initial Chi-Squared = 25517.886155 +# Bestfit Chi-Squared = 1077.242298 +# Degrees-of-Freedom = 231 +# Num. of Iterations = 15 +# Convergence Reason = The relative reduction in the sum of squares is less than atol + +run blind False +run nsubpix 20 +chisq ftol 1.0E-10 +chisq atol 0.001 +chisq miniter 10 +chisq maxiter 1000 +out fits False +out plots spec1d_J0018p2345_KASTb_coadd.pdf +plot dims 3x3 +plot labels True +plot fits False + +data read + spec1d_J0018p2345_KASTb_coadd.dat specid=0 fitrange=[4917.0,4957.0] loadrange=[4907.0,4967.0] resolution=vfwhm(282.946vtie) columns=[wave,flux,error] label=Hb4863 + spec1d_J0018p2345_KASTb_coadd.dat specid=1 fitrange=[4388.0,4428.0] loadrange=[4378.0,4438.0] resolution=vfwhm(282.946vtie) columns=[wave,flux,error] label=Hg4342 + spec1d_J0018p2345_KASTb_coadd.dat specid=2 fitrange=[4146.0,4186.0] loadrange=[4136.0,4196.0] resolution=vfwhm(282.946vtie) columns=[wave,flux,error] label=Hd4103 + spec1d_J0018p2345_KASTb_coadd.dat specid=3 fitrange=[4012.0,4052.0] loadrange=[4002.0,4062.0] resolution=vfwhm(282.946vtie) columns=[wave,flux,error] label=He3970 + spec1d_J0018p2345_KASTb_coadd.dat specid=4 fitrange=[3772.0,3802.0] loadrange=[3767.0,3807.0] resolution=vfwhm(282.946vtie) columns=[wave,flux,error] label=[OII]3727+3729 + spec1d_J0018p2345_KASTb_coadd.dat specid=5 fitrange=[5021.0,5051.0] loadrange=[5016.0,5056.0] resolution=vfwhm(282.946vtie) columns=[wave,flux,error] label=[OIII]4960 + spec1d_J0018p2345_KASTb_coadd.dat specid=6 fitrange=[5065.0,5105.0] loadrange=[5060.0,5110.0] resolution=vfwhm(282.946vtie) columns=[wave,flux,error] label=[OIII]5008 +data end + +model read +fix vfwhm value False +fix gaussian dispersion True + emission + legendre 136.86663193 -8.85011622 scale=1.0,1.0 specid=0 continuum=True + legendre 195.06532795 11.12511048 scale=1.0,1.0 specid=1 continuum=True + legendre 205.78720694 -7.67128769 scale=1.0,1.0 specid=2 continuum=True + legendre 238.51636060 4.99858582 scale=1.0,1.0 specid=3 continuum=True + legendre 254.20178241 -61.15994519 scale=1.0,1.0 specid=4 continuum=True + legendre 154.08236287 31.95230380 scale=1.0,1.0 specid=5 continuum=True + legendre 152.53327300 35.17673753 scale=1.0,1.0 specid=6 continuum=True + gaussian 7.85749594E+03 0.01540425ra 30.00000000dh specid=0 IntFlux=True wave=4862.683 + gaussian 3.59948935E+03 0.01540425ra 30.00000000dh specid=1 IntFlux=True wave=4341.684 + gaussian 1.65562745E+03 0.01540425ra 30.00000000dh specid=2 IntFlux=True wave=4102.892 + gaussian 1.56279300E+03 0.01540425ra 30.00000000dh specid=3 IntFlux=True wave=3971.195 + gaussian 3.46048372E+03 0.01540425ra 30.00000000do specid=4 IntFlux=True wave=3727.092 + gaussian 4.88103542E+03 0.01540425ra 30.00000000do specid=4 IntFlux=True wave=3729.874 + gaussian 4.79330012E+03orata 0.01540425ra 30.00000000dox specid=5 IntFlux=True wave=4960.295 + gaussian 1.43799004E+04oratb 0.01540425ra 30.00000000dox specid=6 IntFlux=True wave=5008.239 +model end + +link read +oratb(orata) = orata*3.0 +link end + +# +# Errors: +# +#emission +# legendre 4.83938676 15.35821602 scale=1.0,1.0 specid=0 continuum=True +# legendre 7.08376397 21.93070237 scale=1.0,1.0 specid=1 continuum=True +# legendre 7.53275066 22.28386888 scale=1.0,1.0 specid=2 continuum=True +# legendre 10.32401860 33.21308807 scale=1.0,1.0 specid=3 continuum=True +# legendre 14.91070787 41.98922735 scale=1.0,1.0 specid=4 continuum=True +# legendre 5.34940707 16.10995876 scale=1.0,1.0 specid=5 continuum=True +# legendre 4.65726318 13.63971131 scale=1.0,1.0 specid=6 continuum=True +# gaussian 1.07203070E+02 0.00000346ra 0.00000000dh specid=0 IntFlux=True wave=4862.683 +# gaussian 1.24076227E+02 0.00000346ra 0.00000000dh specid=1 IntFlux=True wave=4341.684 +# gaussian 1.34774396E+02 0.00000346ra 0.00000000dh specid=2 IntFlux=True wave=4102.892 +# gaussian 1.51170291E+02 0.00000346ra 0.00000000dh specid=3 IntFlux=True wave=3971.195 +# gaussian 2.13668992E+02 0.00000346ra 0.00000000do specid=4 IntFlux=True wave=3727.092 +# gaussian 2.14767652E+02 0.00000346ra 0.00000000do specid=4 IntFlux=True wave=3729.874 +# gaussian 9.65686646E+01orata 0.00000346ra 0.00000000dox specid=5 IntFlux=True wave=4960.295 +# gaussian 0.00000000E+00oratb 0.00000346ra 0.00000000dox specid=6 IntFlux=True wave=5008.239 + +# Convolution Models: +# vfwhm 282.946vtie +# +# Errors: +# +# vfwhm 2.315vtie + +# Shift Models: +# Ashift 0.0000000 +# +# Errors: +# +# Ashift 0.0000000 + + +################################################### +# # +# HERE IS A COPY OF THE INPUT MODEL # +# # +################################################### +# +# run blind False +# run nsubpix 20 +# chisq ftol 1.0E-10 +# chisq atol 0.001 +# chisq miniter 10 +# chisq maxiter 1000 +# out fits False +# out plots spec1d_J0018p2345_KASTb_coadd.pdf +# plot dims 3x3 +# plot labels True +# plot fits False +# data read +# spec1d_J0018p2345_KASTb_coadd.dat specid=0 fitrange=[4917.0,4957.0] loadrange=[4907.0,4967.0] resolution=vfwhm(150.0vtie) columns=[wave,flux,error] label=Hb4863 +# spec1d_J0018p2345_KASTb_coadd.dat specid=1 fitrange=[4388.0,4428.0] loadrange=[4378.0,4438.0] resolution=vfwhm(150.0vtie) columns=[wave,flux,error] label=Hg4342 +# spec1d_J0018p2345_KASTb_coadd.dat specid=2 fitrange=[4146.0,4186.0] loadrange=[4136.0,4196.0] resolution=vfwhm(150.0vtie) columns=[wave,flux,error] label=Hd4103 +# spec1d_J0018p2345_KASTb_coadd.dat specid=3 fitrange=[4012.0,4052.0] loadrange=[4002.0,4062.0] resolution=vfwhm(150.0vtie) columns=[wave,flux,error] label=He3970 +# spec1d_J0018p2345_KASTb_coadd.dat specid=4 fitrange=[3772.0,3802.0] loadrange=[3767.0,3807.0] resolution=vfwhm(150.0vtie) columns=[wave,flux,error] label=[OII]3727+3729 +# spec1d_J0018p2345_KASTb_coadd.dat specid=5 fitrange=[5021.0,5051.0] loadrange=[5016.0,5056.0] resolution=vfwhm(150.0vtie) columns=[wave,flux,error] label=[OIII]4960 +# spec1d_J0018p2345_KASTb_coadd.dat specid=6 fitrange=[5065.0,5105.0] loadrange=[5060.0,5110.0] resolution=vfwhm(150.0vtie) columns=[wave,flux,error] label=[OIII]5008 +# data end +# model read +# fix vfwhm value False +# fix gaussian dispersion True +# emission +# legendre 146.7779 0.01 scale=[1.0,1.0] specid=0 continuum=True +# legendre 222.0621 0.01 scale=[1.0,1.0] specid=1 continuum=True +# legendre 220.9873 0.01 scale=[1.0,1.0] specid=2 continuum=True +# legendre 278.4923 0.01 scale=[1.0,1.0] specid=3 continuum=True +# legendre 289.6434 0.01 scale=[1.0,1.0] specid=4 continuum=True +# legendre 169.9477 0.01 scale=[1.0,1.0] specid=5 continuum=True +# legendre 161.8335 0.01 scale=[1.0,1.0] specid=6 continuum=True +# gaussian 6547.75804488 0.01541ra 30.0dh wave=4862.721 specid=0 IntFlux=True +# gaussian 4038.49307181 0.01541ra 30.0dh wave=4341.684 specid=1 IntFlux=True +# gaussian 2312.54298859 0.01541ra 30.0dh wave=4102.891 specid=2 IntFlux=True +# gaussian 2603.38462088 0.01541ra 30.0dh wave=3971.195 specid=3 IntFlux=True +# gaussian 22165.9457491 0.01541ra 30.0do wave=3727.092 specid=4 IntFlux=True +# gaussian 7758.0810122 0.01541ra 30.0do wave=3729.874 specid=4 IntFlux=True +# gaussian 3654.59671255orata 0.01541ra 30.0dox wave=4960.295 specid=5 IntFlux=True +# gaussian 10963.7901376oratb 0.01541ra 30.0dox wave=5008.239 specid=6 IntFlux=True +# model end +# + diff --git a/linetools/isgm/tests/files/spec1d_J0018p2345_KASTr_coadd.mod.out b/linetools/isgm/tests/files/spec1d_J0018p2345_KASTr_coadd.mod.out new file mode 100644 index 00000000..ccd4c1c1 --- /dev/null +++ b/linetools/isgm/tests/files/spec1d_J0018p2345_KASTr_coadd.mod.out @@ -0,0 +1,106 @@ +# +# Generated by ALIS on 18/10/16 at 15:20:16 +# +# Running Time (hrs) = 0.000257 +# Initial Chi-Squared = 4369.046025 +# Bestfit Chi-Squared = 490.454547 +# Degrees-of-Freedom = 92 +# Num. of Iterations = 6 +# Convergence Reason = The relative reduction in the sum of squares is less than atol + +run blind False +run nsubpix 20 +chisq ftol 1.0E-10 +chisq atol 0.001 +chisq miniter 10 +chisq maxiter 1000 +out fits False +out plots spec1d_J0018p2345_KASTr_coadd.pdf +plot dims 3x3 +plot labels True +plot fits False + +data read + spec1d_J0018p2345_KASTr_coadd.dat specid=0 fitrange=[6630.0,6700.0] loadrange=[6625.0,6705.0] resolution=vfwhm(112.739vtie) columns=[wave,flux,error] label=Ha+[NII] + spec1d_J0018p2345_KASTr_coadd.dat specid=1 fitrange=[6804.0,6854.0] loadrange=[6794.0,6864.0] resolution=vfwhm(112.739vtie) columns=[wave,flux,error] label=[SII]6718+6732 +data end + +model read +fix vfwhm value False +fix gaussian dispersion True + emission + legendre 81.33345461 13.72812617 scale=1.0,1.0 specid=0 continuum=True + legendre 64.71638218 -0.61894027 scale=1.0,1.0 specid=1 continuum=True + gaussian 2.40482980E+04 0.01536343ra 30.00000000dh specid=0 IntFlux=True wave=6564.613 + gaussian 45.59333286nrata 0.01536343ra 30.00000000dn specid=0 IntFlux=True wave=6549.852 + gaussian 136.77999858nratb 0.01536343ra 30.00000000dn specid=0 IntFlux=True wave=6585.277 + gaussian 693.44286861 0.01536343ra 30.00000000ds specid=1 IntFlux=True wave=6718.294 + gaussian 497.13039835 0.01536343ra 30.00000000ds specid=1 IntFlux=True wave=6732.673 +model end + +link read +nratb(nrata) = nrata*3.0 +link end + +# +# Errors: +# +#emission +# legendre 2.86496029 6.92298470 scale=1.0,1.0 specid=0 continuum=True +# legendre 3.21791087 9.56936811 scale=1.0,1.0 specid=1 continuum=True +# gaussian 1.32380884E+02 0.00000117ra 0.00000000dh specid=0 IntFlux=True wave=6564.613 +# gaussian 15.78701551nrata 0.00000117ra 0.00000000dn specid=0 IntFlux=True wave=6549.852 +# gaussian 0.00000000nratb 0.00000117ra 0.00000000dn specid=0 IntFlux=True wave=6585.277 +# gaussian 49.00424871 0.00000117ra 0.00000000ds specid=1 IntFlux=True wave=6718.294 +# gaussian 72.80031637 0.00000117ra 0.00000000ds specid=1 IntFlux=True wave=6732.673 + +# Convolution Models: +# vfwhm 112.739vtie +# +# Errors: +# +# vfwhm 0.875vtie + +# Shift Models: +# Ashift 0.0000000 +# +# Errors: +# +# Ashift 0.0000000 + + +################################################### +# # +# HERE IS A COPY OF THE INPUT MODEL # +# # +################################################### +# +# run blind False +# run nsubpix 20 +# chisq ftol 1.0E-10 +# chisq atol 0.001 +# chisq miniter 10 +# chisq maxiter 1000 +# out fits False +# out plots spec1d_J0018p2345_KASTr_coadd.pdf +# plot dims 3x3 +# plot labels True +# plot fits False +# data read +# spec1d_J0018p2345_KASTr_coadd.dat specid=0 fitrange=[6630.0,6700.0] loadrange=[6625.0,6705.0] resolution=vfwhm(150.0vtie) columns=[wave,flux,error] label=Ha+[NII] +# spec1d_J0018p2345_KASTr_coadd.dat specid=1 fitrange=[6804.0,6854.0] loadrange=[6794.0,6864.0] resolution=vfwhm(150.0vtie) columns=[wave,flux,error] label=[SII]6718+6732 +# data end +# model read +# fix vfwhm value False +# fix gaussian dispersion True +# emission +# legendre 84.1884 0.01 scale=[1.0,1.0] specid=0 continuum=True +# legendre 79.6454 0.01 scale=[1.0,1.0] specid=1 continuum=True +# gaussian 19741.9266277 0.01542ra 30.0dh wave=6564.612 specid=0 IntFlux=True +# gaussian 174.005646513nrata 0.01542ra 30.0dn wave=6549.860 specid=0 IntFlux=True +# gaussian 522.016939538nratb 0.01542ra 30.0dn wave=6585.270 specid=0 IntFlux=True +# gaussian 863.173473057 0.01542ra 30.0ds wave=6718.290 specid=1 IntFlux=True +# gaussian 674.435966799 0.01542ra 30.0ds wave=6732.680 specid=1 IntFlux=True +# model end +# + diff --git a/linetools/isgm/tests/test_init_emsys.py b/linetools/isgm/tests/test_init_emsys.py new file mode 100644 index 00000000..15a8e855 --- /dev/null +++ b/linetools/isgm/tests/test_init_emsys.py @@ -0,0 +1,51 @@ +# Module to run tests on generating AbsSystem + +from __future__ import print_function, absolute_import, division, unicode_literals + +# TEST_UNICODE_LITERALS + +import numpy as np +import os + +import pytest +from astropy import units as u +from astropy.coordinates import SkyCoord + +from linetools.isgm.emsystem import EmSystem + +coord = "00:18:59.32 +23:45:40.32" +radec = SkyCoord(coord, unit=(u.hourangle, u.deg)) + +def data_path(filename): + data_dir = os.path.join(os.path.dirname(__file__), 'files') + return os.path.join(data_dir, filename) + +def test_init(): + emsys = EmSystem(radec, 0.0154) + assert np.isclose(emsys.zem, 0.0154) + + +def test_from_alis(): + # Tests from_dict too + alis_file = data_path('spec1d_J0018p2345_KASTb_coadd.mod.out') + emsys = EmSystem.from_alis(alis_file, 'J001859.32+234540.32') + # + assert np.isclose(emsys.zem, 0.01540425) + + +def test_add_lines_from_alis(): + # Tests from_dict too + alisb_file = data_path('spec1d_J0018p2345_KASTb_coadd.mod.out') + emsys = EmSystem.from_alis(alisb_file, 'J001859.32+234540.32') + # + alisr_file = data_path('spec1d_J0018p2345_KASTr_coadd.mod.out') + emsys.add_emlines_from_alis(alisr_file, chk_z=False) + assert len(emsys._emlines) == 13 + +def test_emsys_to_dict(): + emsys = EmSystem(radec, 0.0154) + todict = emsys.to_dict() + + assert len(todict) == 14 + np.isclose(todict['RA'],4.747166666666666) + np.isclose(todict['DEC'], 23.7612) \ No newline at end of file diff --git a/linetools/lists/sets/llist_v1.2.ascii b/linetools/lists/sets/llist_v1.2.ascii index 181b9ffb..b3d59575 100644 --- a/linetools/lists/sets/llist_v1.2.ascii +++ b/linetools/lists/sets/llist_v1.2.ascii @@ -1432,6 +1432,7 @@ | 6718.294 | [SII] 6716 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | | 6732.673 | [SII] 6731 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | | 7067.198 | HeI 7065 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | +| 7283.356 | HeI 7283 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | | 8752.876 | P12 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | | 8865.217 | Pa11 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | | 9017.385 | Pa10 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | diff --git a/linetools/setup_package.py b/linetools/setup_package.py new file mode 100644 index 00000000..016ce3e2 --- /dev/null +++ b/linetools/setup_package.py @@ -0,0 +1,4 @@ +def get_package_data(): + # Installs the testing data files. Unable to get package_data + # to deal with a directory hierarchy of files, so just explicitly list. + return {'linetools.tests': ['files/*']}