diff --git a/phoebe/backend/backends.py b/phoebe/backend/backends.py index 9c5e2a878..cb263eefd 100644 --- a/phoebe/backend/backends.py +++ b/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] diff --git a/phoebe/backend/spectroscopy.py b/phoebe/backend/spectroscopy.py index 45c75a9f7..0f711cec3 100644 --- a/phoebe/backend/spectroscopy.py +++ b/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 diff --git a/phoebe/backend/universe.py b/phoebe/backend/universe.py index 6e32a511e..68ee20960 100644 --- a/phoebe/backend/universe.py +++ b/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 diff --git a/phoebe/frontend/bundle.py b/phoebe/frontend/bundle.py index 280360abd..4a89b7ed0 100644 --- a/phoebe/frontend/bundle.py +++ b/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) diff --git a/phoebe/parameters/dataset.py b/phoebe/parameters/dataset.py index 2e1c1c958..cb47fe223 100644 --- a/phoebe/parameters/dataset.py +++ b/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 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 diff --git a/phoebe/parameters/figure/dataset.py b/phoebe/parameters/figure/dataset.py index 481dc8c9c..53222d3de 100644 --- a/phoebe/parameters/figure/dataset.py +++ b/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) + + diff --git a/phoebe/parameters/parameters.py b/phoebe/parameters/parameters.py index 8526ff5f0..061bc3154 100644 --- a/phoebe/parameters/parameters.py +++ b/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))