abscomponent and system

linetools/isgm/abscomponent.py

@@ -136,7 +136,7 @@ def from_dict(cls, idict, coord=None, **kwargs):
         # Return
         return slf
 
-    def __init__(self, radec, Zion, z, vlim, Ej=0./u.cm, A=None,
+    def __init__(self, radec, Zion, zcomp, vlim, Ej=0./u.cm, A=None,
                  Ntup=None, comment='', name=None, stars=None):
         """  Initiator
 
@@ -147,7 +147,7 @@ def __init__(self, radec, Zion, z, vlim, Ej=0./u.cm, A=None,
         Zion : tuple 
             Atomic number, ion -- (int,int)
             e.g. (8,1) for OI
-        z : float
+        zcomp : float
             Absorption component redshift
         vlim : Quantity array
             Velocity limits of the component w/r to `z`
@@ -175,7 +175,7 @@ def __init__(self, radec, Zion, z, vlim, Ej=0./u.cm, A=None,
         elif isinstance(radec, SkyCoord):
             self.coord = radec
         self.Zion = Zion
-        self.zcomp = z
+        self.zcomp = zcomp
         self.vlim = vlim
 
         # Optional
@@ -206,7 +206,8 @@ def __init__(self, radec, Zion, z, vlim, Ej=0./u.cm, A=None,
         # Other
         self._abslines = []
 
-    def add_absline(self, absline, tol=0.1*u.arcsec, chk_vel=True, chk_sep=True, **kwargs):
+    def add_absline(self, absline, tol=0.1*u.arcsec, chk_vel=True,
+                    chk_sep=True, vtoler=1., **kwargs):
         """Add an AbsLine object to the component if it satisfies
         all of the rules.
 
@@ -220,8 +221,12 @@ def add_absline(self, absline, tol=0.1*u.arcsec, chk_vel=True, chk_sep=True, **k
           Tolerance on matching coordinates.  Only used if chk_sep=True
         chk_vel : bool, optional
           Perform velocity test (can often be skipped)
+          Insist the bounds of the AbsLine are within 1km/s of the Component
+             (allows for round-off error)
         chk_sep : bool, optional
           Perform coordinate check (expensive)
+        vtoler : float
+          Tolerance for velocity in km/s
         """
         # Perform easy checks
         if chk_sep:
@@ -233,10 +238,11 @@ def add_absline(self, absline, tol=0.1*u.arcsec, chk_vel=True, chk_sep=True, **k
         testE = bool(self.Ej == absline.data['Ej'])
         # Now redshift/velocity
         if chk_vel:
+            dz_toler = (1+self.zcomp)*vtoler/3e5  # Avoid Quantity for speed
             zlim_line = (1+absline.attrib['z'])*absline.analy['vlim']/const.c.to('km/s')
             zlim_comp = (1+self.zcomp)*self.vlim/const.c.to('km/s')
-            testv = (zlim_line[0] >= zlim_comp[0]) & (
-                zlim_line[1] <= zlim_comp[1])
+            testv = (zlim_line[0] >= (zlim_comp[0]-dz_toler)) & (
+                zlim_line[1] <= (zlim_comp[1]+dz_toler))
         else:
             testv = True
         # Combine
@@ -248,8 +254,12 @@ def add_absline(self, absline, tol=0.1*u.arcsec, chk_vel=True, chk_sep=True, **k
         if test:
             self._abslines.append(absline)
         else:
-            warnings.warn('Input absline with wrest={:g} does not match component rules. Not appending'.format(absline.wrest))
-            pdb.set_trace()
+            warnings.warn("Failed add_absline test")
+            print('Input absline with wrest={:g} does not match component rules. Not appending'.format(absline.wrest))
+            if not testv:
+                print("Absline velocities lie beyond component\n Set chk_vel=False to skip this test.")
+            if not testc:
+                print("Absline coordinates do not match.  Best to set them")
 
     def build_table(self):
         """Generate an astropy QTable out of the component.
@@ -606,6 +616,26 @@ def to_dict(self):
         # Return
         return cdict
 
+    def copy(self):
+        """ Generate a copy of itself
+        Returns
+        -------
+        abscomp : AbsComponent
+
+        """
+        # Instantiate with required attributes
+        abscomp = AbsComponent(self.coord, self.Zion, self.zcomp, self.vlim)
+        # Add in the rest
+        attrs = vars(self).keys()
+        for attr in attrs:
+            if attr == '_abslines':
+                for iline in self._abslines:
+                    abscomp._abslines.append(iline.copy())
+            else:
+                setattr(abscomp, attr, getattr(self, attr))
+        # Return
+        return abscomp
+
     def __getitem__(self, attrib):
         """Passback attribute, if it exists
 

linetools/isgm/abssystem.py

@@ -230,11 +230,14 @@ def __init__(self, radec, zabs, vlim, zem=0., abs_type=None,
         # Refs (list of references)
         self.Refs = []
 
-    def add_component(self, abscomp, tol=0.2*u.arcsec, chk_sep=True, **kwargs):
+    def add_component(self, abscomp, tol=0.2*u.arcsec,
+                      chk_sep=True, chk_z=True, overlap_only=False,
+                      vtoler=1., debug=False, **kwargs):
         """Add an AbsComponent object if it satisfies all of the rules.
 
         For velocities, we demand that the new component has a velocity
         range that is fully encompassed by the AbsSystem.
+        We allow a small tolerance for round-off error
 
         Should check for duplicates..
 
@@ -246,25 +249,50 @@ def add_component(self, abscomp, tol=0.2*u.arcsec, chk_sep=True, **kwargs):
           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:
-            test = bool(self.coord.separation(abscomp.coord) < tol)
+            testcoord = bool(self.coord.separation(abscomp.coord) < tol)
         else:
-            test = True
+            testcoord = True
         # Now redshift/velocity
-        zlim_comp = (1+abscomp.zcomp)*abscomp.vlim/c_mks
-        zlim_sys = (1+self.zabs)*self.vlim/c_mks
-        test = test & (zlim_comp[0] >= zlim_sys[0]) & (zlim_comp[1] <= zlim_sys[1])
+        if chk_z:
+            dz_toler = (1+self.zabs)*vtoler/3e5  # Avoid Quantity for speed
+            zlim_comp = abscomp.zcomp + (1+abscomp.zcomp)*(abscomp.vlim/c_mks).decompose()
+            zlim_sys = self.zabs + (1+self.zabs)*(self.vlim/c_mks).decompose()
+            if overlap_only:
+                testz = True
+                if debug:
+                    pdb.set_trace()
+                if np.all(zlim_comp > np.max(zlim_sys+dz_toler)) or np.all(
+                                zlim_comp < np.min(zlim_sys-dz_toler)):
+                    testz = False
+            else:
+                testz = (zlim_comp[0] >= (zlim_sys[0]-dz_toler)) & (
+                    zlim_comp[1] <= (zlim_sys[1]+dz_toler))
 
         # Additional checks (specific to AbsSystem type)
-        test = test & self.chk_component(abscomp)
+        testcomp = self.chk_component(abscomp)
+        test = testcoord & testcomp & testz
 
         # Append?
         if test:
             self._components.append(abscomp)
         else:
-            warnings.warn('Input AbsComponent with does not match AbsSystem rules. Not appending')
+            warnings.warn('Input AbsComponent with Zion={} does not match AbsSystem rules. Not appending'.format(abscomp.Zion))
+        #
+        return test
 
     def chk_component(self, component):
         """Additional checks on the component
@@ -421,6 +449,16 @@ def update_component_colm(self, **kwargs):
         for comp in self._components:
             comp.synthesize_colm(**kwargs)
 
+    def stack_plot(self, **kwargs):
+        """Show a stack plot of the system, if spec are loaded
+        Assumes the data are normalized.
+
+        Parameters
+        ----------
+        """
+        from linetools.analysis import plots as ltap
+        ltap.stack_plot(self.list_of_abslines(), vlim=self.vlim, **kwargs)
+
     def to_dict(self):
         """ Write AbsSystem data to a dict that can be written with JSON
         """

linetools/isgm/tests/test_init_abscomp.py

@@ -12,9 +12,6 @@
 from linetools.isgm.abscomponent import AbsComponent
 from linetools.spectralline import AbsLine
 
-import pdb
-#pdb.set_trace()
-# Set of Input lines
 
 def test_init():
     # Simple properties
@@ -24,6 +21,7 @@ def test_init():
     np.testing.assert_allclose(abscomp.zcomp,1.0)
     print(abscomp)
 
+
 def test_init_failures():
     with pytest.raises(IOError):
         AbsComponent.from_abslines('blah')
@@ -32,6 +30,7 @@ def test_init_failures():
     with pytest.raises(IOError):
         AbsComponent.from_component('blah')
 
+
 def test_init_single_absline():
     # Single AbsLine
     lya = AbsLine(1215.670*u.AA)
@@ -43,6 +42,24 @@ def test_init_single_absline():
     np.testing.assert_allclose(abscomp.zcomp,2.92939)
     print(abscomp)
 
+
+def test_copy():
+    # Single AbsLine
+    lya = AbsLine(1215.670*u.AA)
+    lya.analy['vlim'] = [-300.,300.]*u.km/u.s
+    lya.attrib['z'] = 2.92939
+    abscomp = AbsComponent.from_abslines([lya])
+    # Copy
+    abscomp2 = abscomp.copy()
+    # Checks
+    attrs = vars(abscomp).keys()
+    attrs2 = vars(abscomp2).keys()
+    for attr in attrs:
+        assert attr in attrs2
+    np.testing.assert_allclose(abscomp._abslines[0].attrib['z'],
+                               abscomp2._abslines[0].attrib['z'])
+
+
 def test_init_multi_absline():
     # AbsLine(s)
     lya = AbsLine(1215.670*u.AA)

linetools/isgm/tests/test_init_abssys.py

@@ -21,6 +21,8 @@ def data_path(filename):
     data_dir = os.path.join(os.path.dirname(__file__), 'files')
     return os.path.join(data_dir, filename)
 
+
+
 def test_from_json():
     # Tests from_dict too
     HIsys = LymanAbsSystem.from_json(data_path('HILya_abssys.json'))
@@ -102,3 +104,31 @@ def test_multi_components():
     assert len(LLSsys._components) == 2
 
 
+def test_add_component():
+    radec = SkyCoord(ra=123.1143*u.deg, dec=-12.4321*u.deg)
+    # HI Lya, Lyb
+    lya = AbsLine(1215.670*u.AA)
+    lya.analy['vlim'] = [-300.,300.]*u.km/u.s
+    lya.attrib['z'] = 2.92939
+    lya.attrib['N'] = 1e17 /  u.cm**2
+    lyb = AbsLine(1025.7222*u.AA)
+    lyb.analy['vlim'] = [-300.,300.]*u.km/u.s
+    lyb.attrib['z'] = lya.attrib['z']
+    abscomp = AbsComponent.from_abslines([lya,lyb])
+    abscomp.coord = radec
+    # Instantiate
+    abssys = GenericAbsSystem.from_components([abscomp])
+    # New component
+    oi = AbsLine('OI 1302')
+    oi.analy['vlim'] = [-300.,300.]*u.km/u.s
+    oi.attrib['z'] = lya.attrib['z']
+    abscomp2 = AbsComponent.from_abslines([oi])
+    abscomp2.coord = radec
+    # Standard
+    assert abssys.add_component(abscomp2)
+    # Fail
+    abssys = GenericAbsSystem.from_components([abscomp])
+    abscomp2.vlim = [-400.,300.]*u.km/u.s
+    assert not abssys.add_component(abscomp2)
+    # Overlap
+    assert abssys.add_component(abscomp2, overlap_only=True)

linetools/spectralline.py

@@ -386,7 +386,16 @@ def to_dict(self):
         # Return
         return adict
 
-    # Output
+    def copy(self):
+        """ Generate a copy
+
+        Returns
+        -------
+        newline : SpectralLine
+          copy of the object
+        """
+        return copy.deepcopy(self)  # Cheat for now
+
     def __repr__(self):
         txt = '<{:s}:'.format(self.__class__.__name__)
         try:
@@ -397,8 +406,7 @@ def __repr__(self):
         txt = txt + '>'
         return (txt)
 
-# ###########################################
-# Class for Generic Absorption Line System
+
 class AbsLine(SpectralLine):
     """Class representing a spectral absorption line.
 
@@ -468,7 +476,6 @@ def fill_data(self, trans, linelist=None, closest=False, verbose=True):
         # Additional fundamental attributes for Absorption Line
         self.attrib.update(abs_attrib.copy())
 
-    # Voigt
     def generate_voigt(self, wave=None, **kwargs):
         """ Generate a Voigt profile model for the absorption line in
         a given spectrum.
@@ -501,7 +508,6 @@ def generate_voigt(self, wave=None, **kwargs):
         spec = lav.voigt_from_abslines(wave, self, **kwargs)
         return spec
 
-    # AODM
     def measure_aodm(self, nsig=3., normalize=True):
         """ AODM calculation
 
@@ -544,7 +550,6 @@ def measure_aodm(self, nsig=3., normalize=True):
         # Log
         laa.log_clm(self.attrib)
 
-    # Output
     def __repr__(self):
         txt = '<{:s}:'.format(self.__class__.__name__)
         # Name