Addition of a spectroscopic dataset ('sed').

phoebe/backend/backends.py

@@ -106,7 +106,7 @@ def _needs_mesh(b, dataset, kind, component, compute):
         return False
 
     # Note: 'spe' is included too (for future mesh-based computations)
-    if kind not in ['mesh', 'lc', 'rv', 'lp', 'spe']:
+    if kind not in ['mesh', 'lc', 'rv', 'lp', 'spe', 'sed']:
         return False
 
     # if kind == 'lc' and compute_kind=='phoebe' and b.get_value(qualifier='lc_method', compute=compute, dataset=dataset, context='compute')=='analytical':
@@ -218,7 +218,7 @@ def _extract_from_bundle(b, compute, dataset=None, times=None,
         dataset_compute_ps = b.filter(context='compute', dataset=dataset, compute=compute, **_skip_filter_checks)
         dataset_kind = dataset_ps.kind
         time_qualifier = _timequalifier_by_kind(dataset_kind)
-        if dataset_kind in ['lc', 'spe']:
+        if dataset_kind in ['lc', 'spe', 'sed']:
             # then the Parameters in the model only exist at the system-level
             # and are not tagged by component
             dataset_components = [None]
@@ -510,6 +510,7 @@ def run(self, b, compute, dataset=None, times=[], **kwargs):
 
         logger.debug("rank:{}/{} calling get_packet_and_syns".format(mpi.myrank, mpi.nprocs))
         packet, new_syns = self.get_packet_and_syns(b, compute, dataset, times, **kwargs)
+        print("new_syns = ", new_syns)  # dbg
 
         if mpi.enabled:
             # broadcast the packet to ALL workers
@@ -1193,6 +1194,11 @@ def _run_single_time(self, b, i, time, infolist, **kwargs):
                 else:
                     raise NotImplementedError("spe_method='{}' not supported".format(spe_method))
 
+            elif kind == 'sed' and dataset != previous:
+                wavelengths = b.get_value('wavelengths@'+dataset+'@dataset')
+                previous = dataset
+                k = 0
+
             # now check the kind to see what we need to fill
             if kind=='lp':
                 profile_func = b.get_value(qualifier='profile_func',
@@ -1309,6 +1315,18 @@ def _run_single_time(self, b, i, time, infolist, **kwargs):
                                  info,
                                  index=info['original_index']))
 
+            elif kind=='sed':
+
+                obs = spectroscopy.sed(b, system, wavelengths=wavelengths, info=info, k=k)
+
+                # Note: spectroscopy.sed_integrate() is used instead of system.observe()
+
+                packetlist.append(_make_packet('fluxes',
+                                 obs['flux']*u.W/u.m**3,
+                                 time,
+                                 info,
+                                 index=info['original_index']))
+
             elif kind=='orb':
                 # ts[i], xs[cind][i], ys[cind][i], zs[cind][i], vxs[cind][i], vys[cind][i], vzs[cind][i]
 

phoebe/backend/spectroscopy.py

@@ -18,6 +18,7 @@
 from phoebe.backend import pyterpolmini
 
 sg = None
+sg2 = None
 fluxes = None
 
 def planck(T, lambda_):
@@ -32,7 +33,7 @@ def planck(T, lambda_):
 
 def spe_simple(b, system, wavelengths=None, info={}, k=None):
     """
-    Compute monochromatic flux F_nu.
+    Compute relative monochromatic flux F_nu.
     A simple model of uniform disk(s).
 
     b           .. Bundle object
@@ -41,6 +42,8 @@ def spe_simple(b, system, wavelengths=None, info={}, k=None):
     info        .. dictionary w. 'original_index' to get wavelengths
     k           .. index to run computation
 
+    Note: All wavelengths are computed at once, but returned sequentially 0, 1, 2, ...
+
     Note: Applicable to detached non-eclipsing binaries.
 
     """
@@ -64,7 +67,7 @@ def spe_simple(b, system, wavelengths=None, info={}, k=None):
 
     for i, body in enumerate(system.bodies):
 
-        rv = (system.vzi[i]*u.solRad/u.day).to('km/s').value		# km/s
+        rv = -(system.vzi[i]*u.solRad/u.day).to('km/s').value		# km/s
         area = np.pi*(body.requiv*u.solRad.to('m'))**2			# m^2
         teff = body.teff						# K
         mass = body.masses[body.ind_self]				# M_S
@@ -87,7 +90,7 @@ def spe_simple(b, system, wavelengths=None, info={}, k=None):
 
         s = sg.get_synthetic_spectrum(props, angstroms, order=2, step=step, padding=20.0)
 
-        wave_ = pyterpolmini.doppler_shift(s.wave, -rv)
+        wave_ = pyterpolmini.doppler_shift(s.wave, rv)
         intens_ = pyterpolmini.rotational_broadening(wave_, s.intens, vrot)
         intens__ = pyterpolmini.interpolate_spectrum(wave_, intens_, angstroms)
 
@@ -105,13 +108,11 @@ def spe_simple(b, system, wavelengths=None, info={}, k=None):
 
 def spe_integrate(b, system, wavelengths=None, info={}, k=None):
     """
-    Compute monochromatic flux F_nu.
+    Compute relative monochromatic flux F_nu.
     A complex model w. integration over meshes.
 
     Note: See spe_simple().
 
-    Note: All wavelengths are computed at once, but returned sequentially 0, 1, 2, ...
-
     Note: Applicable to contact or eclipsing binaries.
 
     """
@@ -149,7 +150,7 @@ def spe_integrate(b, system, wavelengths=None, info={}, k=None):
     angstroms = wavelengths*1.0e10		# Ang
     fluxes = np.zeros(len(wavelengths))		# 1
 
-    for i in range(len(abs_intensities)):
+    for i in range(len(Lum)):
         if Lum[i] == 0.0:
             continue
 
@@ -169,6 +170,64 @@ def spe_integrate(b, system, wavelengths=None, info={}, k=None):
     return {'flux': fluxes[0]}
 
 
+def sed_integrate(b, system, wavelengths=None, bandwidths=None, info={}, k=None):
+    """
+    Compute absolute monochromatic flux F_nu.
+
+    Note: See spe_integrate().
+
+    """
+    global sg2
+    global fluxes
+
+    if sg2 is None:
+        sg2 = pyterpolmini.SyntheticGrid(flux_type='absolute', debug=False)
+
+    j = info['original_index']
+    if k > 0:
+        return {'flux': fluxes[j]}
+
+    meshes = system.meshes
+    components = info['component']
+    dataset = info['dataset']
+
+    visibilities = meshes.get_column_flat('visibilities', components)
+
+    if np.all(visibilities==0):
+        return {'flux': np.nan}
+
+    mus = meshes.get_column_flat('mus', components)
+    areas = meshes.get_column_flat('areas_si', components)
+    rvs = (meshes.get_column_flat("rvs:{}".format(dataset), components)*u.solRad/u.d).to(u.m/u.s).value
+    teffs = meshes.get_column_flat('teffs', components)
+    loggs = meshes.get_column_flat('loggs', components)
+    zs = 10.0**meshes.get_column_flat('abuns', components)
+
+    Lum = areas*mus*visibilities		# m^2
+
+    step = 0.1					# Ang
+    angstroms = wavelengths*1.0e10		# Ang
+    fluxes = np.zeros(len(wavelengths))		# 1
+
+    for i in range(len(Lum)):
+        if Lum[i] == 0.0:
+            continue
+
+        props = {'teff': teffs[i], 'logg': loggs[i], 'z': zs[i]}
+
+        s = sg2.get_synthetic_spectrum(props, angstroms, order=2, step=step, padding=20.0)
+
+        rv = rvs[i]*1.0e-3							# km/s
+        wave_ = pyterpolmini.doppler_shift(s.wave, rv)				# Ang
+        intens_ = pyterpolmini.interpolate_spectrum(wave_, s.intens, angstroms)	# erg s^-1 cm^-2 Ang^-1 
+        intens_ *= 1.0e7							# W m^-2 m^-1
+
+        fluxes += Lum[i]*intens_
+
+    return {'flux': fluxes[0]}
+
+
 spe = spe_integrate
+sed = sed_integrate
 
 

phoebe/backend/universe.py

@@ -641,6 +641,8 @@ def gaussian(sv):
 
         # Note: See spectroscopy.spe_integrate().
 
+        # Note: See spectroscopy.sed_integrate().
+
 
 class Body(object):
     """
@@ -1275,7 +1277,7 @@ def from_bundle(cls, b, component, compute=None,
             kwargs.setdefault('mesh_init_phi', b.get_value(qualifier='mesh_init_phi', compute=compute, component=component, unit=u.rad, mesh_init_phi=kwargs.get('mesh_init_phi', None), **_skip_filter_checks))
 
         # Note: 'spe' is included too (for future mesh-based computations)
-        datasets_intens = [ds for ds in b.filter(kind=['lc', 'rv', 'lp', 'spe'], context='dataset').datasets if ds != '_default']
+        datasets_intens = [ds for ds in b.filter(kind=['lc', 'rv', 'lp', 'spe', 'sed'], context='dataset').datasets if ds != '_default']
         datasets_lp = [ds for ds in b.filter(kind='lp', context='dataset', **_skip_filter_checks).datasets if ds != '_default']
         atm_override = kwargs.pop('atm', None)
         if isinstance(atm_override, dict):
@@ -1772,6 +1774,19 @@ def _populate_spe(self, dataset, **kwargs):
 
         return cols
 
+    def _populate_sed(self, dataset, **kwargs):
+        """
+        Populate columns necessary for a spectroscopic dataset
+
+        """
+        logger.debug("{}._populate_sed(dataset={})".format(self.component, dataset))
+
+        cols = self._populate_rv(dataset, **kwargs)
+
+        # Note: See observe().
+
+        return cols
+
     def _populate_lc(self, dataset, ignore_effects=False, **kwargs):
         """
         Populate columns necessary for an LC dataset

phoebe/frontend/bundle.py

@@ -12123,6 +12123,9 @@ def _scale_fluxes_cfit(fluxes, scale_factor):
                 for spe_param in ml_params.filter(qualifier=['wavelengths'], kind='spe', **_skip_filter_checks).to_list():
                     spe_param.set_value(self.get_value(spe_param.twig+'@dataset'), ignore_readonly=True)
 
+                for sed_param in ml_params.filter(qualifier=['wavelengths'], kind='sed', **_skip_filter_checks).to_list():
+                    sed_param.set_value(self.get_value(sed_param.twig+'@dataset'), ignore_readonly=True)
+
                 ml_addl_params += [StringParameter(qualifier='comments', value=kwargs.get('comments', computeparams.get_value(qualifier='comments', default='', **_skip_filter_checks)), description='User-provided comments for this model.  Feel free to place any notes here.')]
                 self._attach_params(ml_params+ml_addl_params, check_copy_for=False, **metawargs)
 

phoebe/parameters/dataset.py

@@ -805,6 +805,30 @@ def spe(syn=False, as_ps=True, **kwargs):
 
     return ParameterSet(params) if as_ps else params, constraints
 
+def sed(syn=False, as_ps=True, **kwargs):
+    """
+    Create a <phoebe.parameters.ParameterSet> for an SED dataset.
+
+    Note: See lc().
+
+    """
+
+    params, constraints = [], []
+
+    params += [FloatArrayParameter(qualifier='times', value=kwargs.get('times', []), required_shape=[None], readonly=syn, default_unit=u.d, description='Model (synthetic) times' if syn else 'Observed times')]
+    params += [FloatArrayParameter(qualifier='wavelengths', value=kwargs.get('wavelengths', []), required_shape=[None], readonly=syn, default_unit=u.m, description='Synthetic wavelengths' if syn else 'Observed wavelengths')]
+    params += [FloatArrayParameter(qualifier='fluxes', value=_empty_array(kwargs, 'fluxes'), required_shape=[None] if not syn else None, readonly=syn, default_unit=u.W/u.m**3, description='Synthetic fluxes' if syn else 'Observed fluxes')]
+
+    if not syn:
+        params += [FloatArrayParameter(qualifier='compute_times', value=kwargs.get('compute_times', []), required_shape=[None], default_unit=u.d, description='Times to use during run_compute. If empty, will use times parameter')]
+        params += [FloatArrayParameter(qualifier='sigmas', value=_empty_array(kwargs, 'sigmas'), required_shape=[None], default_unit=u.W/u.m**3, description='Observed uncertainty of flux')]
+
+    lc_params, lc_constraints = lc(syn=syn, as_ps=False, is_lc=False, **kwargs)
+    params += lc_params
+    constraints += lc_constraints
+
+    return ParameterSet(params) if as_ps else params, constraints
+
 
 # del _empty_array
 # del deepcopy

phoebe/parameters/figure/dataset.py

@@ -209,10 +209,28 @@ def spe(b, **kwargs):
     params += [ChoiceParameter(qualifier='x', value=kwargs.get('x', 'times'), choices=['times', 'wavelengths'], description='Array to plot along x-axis')]
     params += [ChoiceParameter(qualifier='y', value=kwargs.get('y', 'fluxes'), choices=['fluxes', 'residuals'], description='Array to plot along y-axis')]
 
-    params += _label_units_lims('times', visible_if='x:times', default_unit=u.d, is_default=True, **kwargs)
-    params += _label_units_lims('wavelengths', visible_if='x:wavelengths', default_unit=u.m, is_default=True, **kwargs)
-    params += _label_units_lims('fluxes', default_unit=u.dimensionless_unscaled, is_default=True, **kwargs)
+    params += _label_units_lims('x', visible_if='x:times', default_unit=u.d, is_default=True, **kwargs)
+    params += _label_units_lims('x', visible_if='x:wavelengths', default_unit=u.m, is_default=False, **kwargs)
+    params += _label_units_lims('y', visible_if='y:fluxes', default_unit=u.dimensionless_unscaled, is_default=True, **kwargs)
+
+    return ParameterSet(params)
+
+def sed(b, **kwargs):
+    params = []
+
+    params += [SelectParameter(qualifier='contexts', value=kwargs.get('contexts', '*'), choices=['dataset', 'model'], description='Contexts to include in the plot')]
+    params += [SelectParameter(qualifier='datasets', value=kwargs.get('datasets', '*'), choices=[''], description='Datasets to include in the plot')]
+    params += [SelectParameter(qualifier='models', value=kwargs.get('models', '*'), choices=[''], description='Models to include in the plot')]
+
+    params += [ChoiceParameter(qualifier='x', value=kwargs.get('x', 'times'), choices=['times', 'wavelengths'], description='Array to plot along x-axis')]
+    params += [ChoiceParameter(qualifier='y', value=kwargs.get('y', 'fluxes'), choices=['fluxes', 'residuals'], description='Array to plot along y-axis')]
+
+    params += _label_units_lims('x', visible_if='x:times', default_unit=u.d, is_default=True, **kwargs)
+    params += _label_units_lims('x', visible_if='x:wavelengths', default_unit=u.m, is_default=False, **kwargs)
+    params += _label_units_lims('y', visible_if='y:fluxes', default_unit=u.W/u.m**3, is_default=True, **kwargs)
 
     return ParameterSet(params)
 
 
+
+

phoebe/parameters/parameters.py

@@ -179,7 +179,7 @@
 _forbidden_labels += ['bol']
 
 # forbid all kinds
-_forbidden_labels += ['lc', 'rv', 'lp', 'sp', 'orb', 'mesh', 'spe']
+_forbidden_labels += ['lc', 'rv', 'lp', 'sp', 'orb', 'mesh', 'spe', 'sed']
 _forbidden_labels += ['star', 'orbit', 'envelope']
 _forbidden_labels += ['spot', 'pulsation']
 _forbidden_labels += ['phoebe', 'legacy', 'jktebop', 'photodynam', 'ellc']
@@ -3726,6 +3726,8 @@ def calculate_residuals(self, model=None, dataset=None, component=None,
             qualifier = 'rvs'
         elif dataset_kind == 'spe':
             qualifier = 'fluxes'
+        elif dataset_kind == 'sed':
+            qualifier = 'fluxes'
         else:
             # TODO: lp compared for a given time interpolating in wavelength?
             # NOTE: add to documentation if adding support for other datasets
@@ -4178,7 +4180,7 @@ def _handle_additional_calls(ps, kwargss):
                 return_ += this_return
             return _handle_additional_calls(ps, return_)
 
-        if len(ps.components) > 1 and ps.context in ['model', 'dataset'] and ps.kind not in ['lc', 'spe']:
+        if len(ps.components) > 1 and ps.context in ['model', 'dataset'] and ps.kind not in ['lc', 'spe', 'sed']:
             # lc, spe have per-component passband-dependent parameters in the dataset which are not plottable
             return_ = []
             for component in ps.filter(component=kwargs.get('component', filter_kwargs.get('component', None)), **_skip_filter_checks).exclude(qualifier=['*_phases', 'phases_*'], **_skip_filter_checks).components:
@@ -4768,6 +4770,12 @@ def _kwargs_fill_dimension(kwargs, direction, ps):
                         'z': 0}
             sigmas_avail = ['fluxes']
 
+        elif ps.kind == 'sed':
+            defaults = {'x': 'times',
+                        'y': 'fluxes',
+                        'z': 0}
+            sigmas_avail = ['fluxes']
+
         else:
             logger.debug("could not find plotting defaults for ps.meta: {}, ps.twigs: {}".format(ps.meta, ps.twigs))
             raise NotImplementedError("defaults for kind {} (dataset: {}) not yet implemented".format(ps.kind, ps.dataset))