Merge pull request #141 from gianninapr/specinputfix

Fixed input to  production rate functions so as to follow sbpy rules

Author: mommermi

docs/sbpy/spectroscopy.rst

@@ -14,10 +14,12 @@ JPLSpec Constants and Line Intensity Integral Conversion
 `sbpy.spectroscopy` has a function called ``molecular_data`` which takes care of
 querying the JPL Molecular Spectral Catalog through the use of
 `astroquery.jplspec` and calculates all the necessary constants needed for both
-production rate and Einstein coefficient calculations. The function
-``intensity_conversion`` takes care of converting the intensity line integral at
-300 K found in JPL Spec catalog and convert it closer to the temperature given
-by the user.
+production rate and Einstein coefficient calculations. ``molecular_data``
+returns an `sbpy.data.phys` instance with quantities in the order of: Transition
+Frequency, Temperature, Integrated line intensity at 300 K, and Partition
+function at 300 K. The function ``intensity_conversion`` takes care of
+converting the intensity line integral at 300 K found in JPL Spec catalog and
+convert it closer to the temperature given by the user.
 
 .. code-block:: python
 
@@ -82,18 +84,21 @@ section.
 .. code-block:: python
 
   >>> from sbpy.spectroscopy import prodrate_np
+  >>> from astropy.time import Time
+  >>> from sbpy.data import Ephem
   >>> temp_estimate = 33. * u.K
-  >>> target = '900918'
+  >>> target = '103P'
   >>> vgas = 0.8 * u.km / u.s
   >>> aper = 30 * u.m
   >>> b = 1.13
   >>> mol_tag = 27001
   >>> transition_freq = 265.886434 * u.MHz
   >>> spectra = 1.22 * u.K * u.km / u.s
-  >>> time = '2010-11-3 00:48:06'
+  >>> time = Time('2010-11-3 00:48:06', format='iso')
+  >>> ephemobj = Ephem(target, epochs=time.jd, id_type='id')
   >>> q = prodrate_np(spectra, temp_estimate, transition_freq,
-                            mol_tag, time, target, vgas, aper,
-                            b=b, id_type='id')
+                            mol_tag, ephemobj, vgas, aper,
+                            b=b)
 
   >>> q
   <Quantity 1.0432591198553935e+25 1 / s>
@@ -113,6 +118,8 @@ model does account for the effects of photolysis.
 .. code-block:: python
 
   >>> from sbpy.activity.gas import Haser
+  >>> from astropy.time import Time
+  >>> from sbpy.data import Ephem
   >>> coma = Haser(Q, v, parent)
   >>> Q = spec.production_rate(coma, molecule='H2O')
 
@@ -126,14 +133,15 @@ model does account for the effects of photolysis.
   >>> vgas = 0.5 * u.km / u.s
 
   >>> time = '2017-12-22 05:24:20'
+  >>> ephemobj = Ephem(target, epochs=time.jd)
   >>> spectra = 0.26 * u.K * u.km / u.s
 
   >>> parent = photo_timescale('CO') * vgas
 
   >>> coma = Haser(Q_estimate, vgas, parent)
 
   >>> Q = spec.production_rate(coma, spectra, temp_estimate,
-                               transition_freq, mol_tag, time, target,
+                               transition_freq, mol_tag, ephemobj,
                                aper=aper, b=b)
 
   >>> print(Q)
@@ -169,7 +177,7 @@ The names of the built-in sources are stored as an internal array.  They can be
 
   >>> from sbpy.spectroscopy import Sun
   >>> Sun.show_builtin()
-      name                                description                           
+      name                                description
   ------------ -----------------------------------------------------------------
      E490_2014                E490-00a (2014) reference solar spectrum (Table 3)
    E490_2014LR E490-00a (2014) low resolution reference solar spectrum (Table 4)

sbpy/spectroscopy/core.py

@@ -222,24 +222,17 @@ def einstein_coeff(temp_estimate, transition_freq, mol_tag):
     return au / u.s
 
 
-def photod_rate(time, time_scale, target, id_type, observatory, time_format,
-                mol_tag):
+def photod_rate(ephemobj, mol_tag):
     # imported here to avoid circular dependency with activity.gas
     from ..activity.gas import photo_timescale
+    from ..data import Ephem
 
-    epoch = Time(time, scale=time_scale, format=time_format)
-    obj = Horizons(id=target, epochs=epoch.jd, location=observatory,
-                   id_type=id_type)
+    if not isinstance(ephemobj, Ephem):
+        raise ValueError('ephemobj must be a `sbpy.data.ephem` instance.')
 
-    try:
-        orb = obj.ephemerides()
-
-    except ValueError:
-
-        raise
-
-    delta = (orb['delta'].data * u.au).to('m')
-    r = (orb['r'].data)
+    orb = ephemobj
+    delta = (orb['delta']).to('m')
+    r = (orb['r'])
     cat = JPLSpec.get_species_table()
     mol = cat[cat['TAG'] == mol_tag]
     name = mol['NAME'].data[0]
@@ -258,8 +251,7 @@ def photod_rate(time, time_scale, target, id_type, observatory, time_format,
 
 
 def total_number(integrated_line, temp_estimate, transition_freq, mol_tag,
-                 aper, b, time, target, time_scale, id_type, observatory,
-                 time_format):
+                 aper, b, ephemobj):
     """
     Basic equation relating number of molecules with observed integrated flux.
     This is given by equation 10 in
@@ -285,8 +277,8 @@ def total_number(integrated_line, temp_estimate, transition_freq, mol_tag,
                                                   transition_freq,
                                                   mol_tag))).decompose()
 
-    photod = photod_rate(time, time_scale, target, id_type, observatory,
-                         time_format, mol_tag)
+    photod = photod_rate(ephemobj, mol_tag)
+
     beta = photod["beta"]
 
     sigma = (1./2. * beta * b * con.c / (transition_freq * aper)).value
@@ -456,10 +448,8 @@ def fit(self, spec):
         """
 
     def prodrate_np(self, spectra, temp_estimate, transition_freq,
-                    mol_tag, time, target, vgas=1 * u.km/u.s,
-                    aper=25 * u.m, observatory='500', b=1.2,
-                    time_format='iso', time_scale='utc',
-                    id_type='designation'):
+                    mol_tag, ephemobj, vgas=1 * u.km/u.s,
+                    aper=25 * u.m, b=1.2):
         """
         | Returns production rate based on Drahus 2012 model referenced. Includes
         | no photodissociation
@@ -485,14 +475,8 @@ def prodrate_np(self, spectra, temp_estimate, transition_freq,
         temp_estimate : `~astropy.units.Quantity`
             Estimated temperature in Kelvins
 
-        time : str
-            Time of observation of any format supported by `~astropy.time`
-
-        target : str
-            | Target designation, if there is more than one aparition you
-            | will be prompted to pick a more specific identifier from a
-            | displayed table and change the parameter id_type to 'id'.
-            | Look at `~astroquery.jplhorizons` for more information.
+        ephemobj: `~sbpy.data.Ephem` object
+            An `~sbpy.data.Ephem` object that holds ephemerides information
 
         mol_tag : int or str
             Molecule identifier. Make sure it is an exclusive identifier.
@@ -503,27 +487,11 @@ def prodrate_np(self, spectra, temp_estimate, transition_freq,
         aper : `~astropy.units.Quantity`
             Telescope aperture in meters. Default is 25 m
 
-        observatory : str
-            | Observatory identifier as per `~astroquery.jplhorizons`
-            | Default is geocentric ('500')
-
         b : int
             | Dimensionless factor intrinsic to every antenna. Typical
             | value, and the default for this model, is 1.22. See
             | references for more information on this parameter.
 
-        time_format : str
-            | Time format, see `~astropy.time` for more information
-            | Default is 'iso' which corresponds to 'YYYY-MM-DD HH:MM:SS'
-
-        time_scale : str
-            | Time scale, see `~astropy.time` for mor information.
-            | Default is 'utc'
-
-        id_type : str
-            | ID type for target. See `~astroquery.jplhorizons` for more.
-            | Default is 'designation'
-
         Returns
         -------
         q : `~astropy.units.Quantity`
@@ -533,11 +501,15 @@ def prodrate_np(self, spectra, temp_estimate, transition_freq,
         --------
         >>> import astropy.units as u  # doctest: +SKIP
 
+        >>> from astropy.time import Time # doctest: +SKIP
+
+        >>> from sbpy.data import Ephem # doctest: +SKIP
+
         >>> from sbpy.spectroscopy import prodrate_np  # doctest: +SKIP
 
         >>> temp_estimate = 33. * u.K  # doctest: +SKIP
 
-        >>> target = '900918'  # doctest: +SKIP
+        >>> target = '103P'  # doctest: +SKIP
 
         >>> vgas = 0.8 * u.km / u.s  # doctest: +SKIP
 
@@ -551,11 +523,12 @@ def prodrate_np(self, spectra, temp_estimate, transition_freq,
 
         >>> spectra = 1.22 * u.K * u.km / u.s  # doctest: +SKIP
 
-        >>> time = '2010-11-3 00:48:06'  # doctest: +SKIP
+        >>> time = Time('2010-11-3 00:48:06', format='iso')  # doctest: +SKIP
+
+        >>> ephemobj = Ephem(target, epochs=time.jd, id_type='id') # doctest: +SKIP
 
         >>> q = prodrate_np(spectra, temp_estimate, transition_freq, # doctest: +SKIP
-                            mol_tag, time, target, vgas, aper,
-                            b=b, id_type='id')
+                            mol_tag, ephemobj, vgas, aper, b=b)
 
         >>> q  # doctest: +SKIP
         <Quantity 1.0432591198553935e+25 1 / s>
@@ -577,17 +550,6 @@ def prodrate_np(self, spectra, temp_estimate, transition_freq,
         assert isinstance(vgas, u.Quantity)
         assert isinstance(aper, u.Quantity)
         assert isinstance(spectra, u.Quantity)  # K * km / s
-        assert isinstance(time, str), "Input time as string, i.e. '2018-05-14'"
-        assert isinstance(time_scale, str), "Input time scale as string, i.e.\
-                                             'utc' see astropy.time.Time for \
-                                              more info"
-        assert isinstance(target, str), "Object name should be a string and\
-                                         should be identifiable through \
-                                         astroquery.jplhorizons"
-        assert isinstance(observatory, str), "Observatory should be a string\
-                                              and identified through \
-                                              astroquery.jplhorizons,\
-                                              i.e. '500' (geocentric)"
 
         if type(transition_freq) == list:
 
@@ -631,8 +593,8 @@ def prodrate_np(self, spectra, temp_estimate, transition_freq,
             gu = sum(gu_list) / float(len(gu_list))
             energy_J = sum(energy_J_list) / float(len(energy_J_list))
 
-            photod = photod_rate(time, time_scale, target, id_type, observatory,
-                                 time_format, mol_tag)
+            photod = photod_rate(ephemobj, mol_tag)
+
             delta = photod["delta"]
 
             calc = ((16*np.pi*k*t_freq.decompose() *
@@ -662,8 +624,8 @@ def prodrate_np(self, spectra, temp_estimate, transition_freq,
 
             au = einstein_coeff(temp_estimate, transition_freq, mol_tag)
 
-            photod = photod_rate(time, time_scale, target, id_type, observatory,
-                                 time_format, mol_tag)
+            photod = photod_rate(ephemobj, mol_tag)
+
             delta = photod["delta"]
 
             calc = ((16*np.pi*k*t_freq.decompose() *
@@ -678,10 +640,8 @@ def prodrate_np(self, spectra, temp_estimate, transition_freq,
         return q
 
     def production_rate(self, coma, integrated_line, temp_estimate,
-                        transition_freq, mol_tag, time, target,
-                        aper=25 * u.m, observatory='500', b=1.2,
-                        time_format='iso', time_scale='utc',
-                        id_type='designation'):
+                        transition_freq, mol_tag, ephemobj,
+                        aper=25 * u.m, b=1.2):
         """
         Calculate production rate for `GasComa`
 
@@ -702,39 +662,17 @@ def production_rate(self, coma, integrated_line, temp_estimate,
         mol_tag : int or str
             Molecule identifier. Make sure it is an exclusive identifier.
 
-        time : str
-            Time of observation of any format supported by `~astropy.time`
-
-        target : str
-            | Target designation, if there is more than one aparition you
-            | will be prompted to pick a more specific identifier from a
-            | displayed table and change the parameter id_type to 'id'.
-            | Look at `~astroquery.jplhorizons` for more information.
+        ephemobj: `~sbpy.data.Ephem` object
+            An `~sbpy.data.Ephem` object that holds ephemerides information
 
         aper : `~astropy.units.Quantity`
             Telescope aperture in meters. Default is 25 m
 
-        observatory : str
-            | Observatory identifier as per `~astroquery.jplhorizons`
-            | Default is geocentric ('500')
-
         b : int
             | Dimensionless factor intrinsic to every antenna. Typical
             | value, and the default for this model, is 1.22. See
             | references for more information on this parameter.
 
-        time_format : str
-            | Time format, see `~astropy.time` for more information
-            | Default is 'iso' which corresponds to 'YYYY-MM-DD HH:MM:SS'
-
-        time_scale : str
-            | Time scale, see `~astropy.time` for mor information.
-            | Default is 'utc'
-
-        id_type : str
-            | ID type for target. See `~astroquery.jplhorizons` for more.
-            | Default is 'designation'
-
         Returns
         -------
         Q : `~astropy.units.Quantity`
@@ -743,6 +681,9 @@ def production_rate(self, coma, integrated_line, temp_estimate,
         Examples
         --------
         >>> from sbpy.activity.gas import Haser # doctest: +SKIP
+        >>> from sbpy.data import Ephem # doctest: +SKIP
+        >>> from astropy.time import Time # doctest: +SKIP
+
         >>> coma = Haser(Q, v, parent) # doctest: +SKIP
         >>> Q = spec.production_rate(coma, molecule='H2O') # doctest: +SKIP
 
@@ -755,15 +696,16 @@ def production_rate(self, coma, integrated_line, temp_estimate,
         >>> b = 0.74 # doctest: +SKIP
         >>> vgas = 0.5 * u.km / u.s # doctest: +SKIP
 
-        >>> time = '2017-12-22 05:24:20' # doctest: +SKIP
+        >>> time = Time('2017-12-22 05:24:20', format = 'iso') # doctest: +SKIP
+        >>> ephemobj = Ephem(target, epochs=time.jd) # doctest: +SKIP
         >>> spectra = 0.26 * u.K * u.km / u.s # doctest: +SKIP
 
         >>> parent = photo_timescale('CO') * vgas # doctest: +SKIP
 
         >>> coma = Haser(Q_estimate, vgas, parent) # doctest: +SKIP
 
         >>> Q = spec.production_rate(coma, spectra, temp_estimate, # doctest: +SKIP
-                                     transition_freq, mol_tag, time, target,
+                                     transition_freq, mol_tag, ephemobj,
                                      aper=aper, b=b) # doctest: +SKIP
 
         >>> print(Q) # doctest: +SKIP
@@ -785,8 +727,7 @@ def production_rate(self, coma, integrated_line, temp_estimate,
         # integrated_line = self.integrated_flux(transition_freq) - not yet implemented
 
         molecules = total_number(integrated_line, temp_estimate,
-                                 transition_freq, mol_tag, aper, b, time, target,
-                                 time_scale, id_type, observatory, time_format)
+                                 transition_freq, mol_tag, aper, b, ephemobj)
 
         model_molecules = coma.total_number(aper)