Skip to content

Commit

Permalink
Fixes blackbody, but breaks all other passbands...
Browse files Browse the repository at this point in the history 
Not quite sure what causes the breaks yet, but it has to do with Passband.load() method.
  • Loading branch information
@aprsa
aprsa committed Mar 4, 2025
1 parent f1318f0 commit 7ffe0f6
Show file tree
Hide file tree
Showing 3 changed files with 78 additions and 110 deletions.
114 changes: 37 additions & 77 deletions phoebe/atmospheres/passbands.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -458,8 +458,8 @@ def save(self, archive, overwrite=True, update_timestamp=True, export_to_pre25=F
bb_func = Table({'teff': bb_func_energy[0], 'logi_e': bb_func_energy[1], 'logi_p': bb_func_photon[1]}, meta={'extname': 'BB_FUNC'})
data.append(fits.table_to_hdu(bb_func))
else:
data.append(fits.ImageHDU(self.ndp[atm.name].table['inorm@energy']['grid'], name=f'{atm.prefix.upper()}NEGRID'))
data.append(fits.ImageHDU(self.ndp[atm.name].table['inorm@photon']['grid'], name=f'{atm.prefix.upper()}NPGRID'))
data.append(fits.ImageHDU(self.ndp[atm.name].table['inorm@energy']['grid'], name=f'{atm.prefix}negrid'))
data.append(fits.ImageHDU(self.ndp[atm.name].table['inorm@photon']['grid'], name=f'{atm.prefix}npgrid'))

content.append(f'{atm.name}:Inorm')

Expand Down Expand Up @@ -553,33 +553,39 @@ def load(cls, archive, load_content=True):
self.ptf_photon = lambda wl: interpolate.splev(wl, self.ptf_photon_func)

if load_content:
stored_atms = set([content.split(':')[0] for content in self.content])

# TODO: replace with < parse('2.5.0') when 2.5.0 is released
if parse(self.phoebe_version) != parse('2.4.17.dev+feature-blending'):
if 'blackbody:Inorm' in self.content:
# 2.4.17+ passbands include bb_teffs; older versions do not.
if 'bb_teffs' not in hdul:
bb_teffs = Table({'teff': hdul['bb_func'].data['teff']}, meta={'extname': 'bb_teffs'})
hdul.append(fits.table_to_hdu(bb_teffs))

if 'blackbody' in stored_atms:
# blackbody atmospheres are reworked, so we need to
# recompute the intensities:
self.compute_intensities(
atm=models.BlackbodyModelAtmosphere(),
include_mus=False,
include_ld=False,
include_extinction='blackbody:ext' in self.content,
verbose=False
)
ndp = self.ndp['blackbody'] # initialized by compute_intensities()

hdul['bb_teffs'].data.columns.change_name('teff', 'teffs')
# store axes:
hdul['bb_teffs'] = fits.table_to_hdu(Table({'teffs': ndp.axes[0]}, meta={'extname': 'bb_teffs'}))
if 'blackbody:ext' in self.content:
hdul['bb_ebvs'].data.columns.change_name('ebv', 'ebvs')
hdul['bb_rvs'].data.columns.change_name('rv', 'rvs')
hdul['BBEGRID'].header['EXTNAME'] = 'bbxegrid'
hdul['BBPGRID'].header['EXTNAME'] = 'bbxpgrid'
hdul['bb_ebvs'] = fits.table_to_hdu(Table({'ebvs': ndp.axes[1]}, meta={'extname': 'bb_ebvs'}))
hdul['bb_rvs'] = fits.table_to_hdu(Table({'rvs': ndp.axes[2]}, meta={'extname': 'bb_rvs'}))

# store tables:
hdul.append(fits.ImageHDU(self.ndp['blackbody'].table['inorm@energy']['grid'], name='bbnegrid'))
hdul.append(fits.ImageHDU(self.ndp['blackbody'].table['inorm@photon']['grid'], name='bbnpgrid'))
if 'blackbody:ext' in self.content:
hdul.append(fits.ImageHDU(self.ndp['blackbody'].table['ext@energy']['grid'], name='bbxegrid'))
hdul.append(fits.ImageHDU(self.ndp['blackbody'].table['ext@photon']['grid'], name='bbxpgrid'))

stored_atms = set([content.split(':')[0] for content in self.content])
# clean up:
del self.ndp['blackbody']

# rename columns in old tables:
if 'ck2004' in stored_atms:
hdul['ck_teffs'].data.columns.change_name('teff', 'teffs')
hdul['ck_loggs'].data.columns.change_name('logg', 'loggs')
Expand Down Expand Up @@ -607,52 +613,52 @@ def load(cls, archive, load_content=True):
self.extern_wd_idx = header['wd_idx']

# Model atmospheres in the passband file:
atms = set([entry.split(':')[0] for entry in self.content])
stored_atms = set([entry.split(':')[0] for entry in self.content])

# We have to iterate over available atms rather than stored atms because
# the stored atms may not be available in the current version of PHOEBE.
available_atms = models.ModelAtmosphere.__subclasses__()
for atm in available_atms:
if atm.name not in atms:
if atm.name not in stored_atms:
continue

if atm.external:
continue

basic_axes = tuple([np.array(list(hdul[f'{atm.prefix}_{name}'].data[name])) for name in atm.basic_axis_names])
basic_axes = tuple(np.asarray(hdul[f'{atm.prefix}_{name}'].data[name]) for name in atm.basic_axis_names)
self.ndp[atm.name] = ndpolator.Ndpolator(basic_axes=basic_axes)

if f'{atm.name}:Inorm' in self.content:
self.ndp[atm.name].register('inorm@photon', None, hdul[f'{atm.prefix}npgrid'].data)
self.ndp[atm.name].register('inorm@energy', None, hdul[f'{atm.prefix}negrid'].data)
self.ndp[atm.name].register(name='inorm@photon', associated_axes=None, grid=hdul[f'{atm.prefix}npgrid'].data)
self.ndp[atm.name].register(name='inorm@energy', associated_axes=None, grid=hdul[f'{atm.prefix}negrid'].data)

if f'{atm.name}:Imu' in self.content:
atm_photon_grid = hdul[f'{atm.prefix}FPGRID'].data
atm_energy_grid = hdul[f'{atm.prefix}FEGRID'].data
atm_photon_grid = hdul[f'{atm.prefix}fpgrid'].data
atm_energy_grid = hdul[f'{atm.prefix}fegrid'].data

# normal passband intensities:
self.ndp[atm.name].register('inorm@photon', None, atm_photon_grid[..., -1, :])
self.ndp[atm.name].register('inorm@energy', None, atm_energy_grid[..., -1, :])
self.ndp[atm.name].register(name='inorm@photon', associated_axes=None, grid=atm_photon_grid[..., -1, :])
self.ndp[atm.name].register(name='inorm@energy', associated_axes=None, grid=atm_energy_grid[..., -1, :])

# specific passband intensities:
self.ndp[atm.name].register('imu@photon', (atm.mus,), atm_photon_grid)
self.ndp[atm.name].register('imu@energy', (atm.mus,), atm_energy_grid)
self.ndp[atm.name].register(name='imu@photon', associated_axes=(atm.mus,), grid=atm_photon_grid)
self.ndp[atm.name].register(name='imu@energy', associated_axes=(atm.mus,), grid=atm_energy_grid)

if f'{atm.name}:ld' in self.content:
self.ndp[atm.name].register('ld@photon', None, hdul[f'{atm.prefix}legrid'].data)
self.ndp[atm.name].register('ld@energy', None, hdul[f'{atm.prefix}lpgrid'].data)
self.ndp[atm.name].register(name='ld@photon', associated_axes=None, grid=hdul[f'{atm.prefix}legrid'].data)
self.ndp[atm.name].register(name='ld@energy', associated_axes=None, grid=hdul[f'{atm.prefix}lpgrid'].data)

if f'{atm.name}:ldint' in self.content:
self.ndp[atm.name].register('ldint@photon', None, hdul[f'{atm.prefix}iegrid'].data)
self.ndp[atm.name].register('ldint@energy', None, hdul[f'{atm.prefix}ipgrid'].data)
self.ndp[atm.name].register(name='ldint@photon', associated_axes=None, grid=hdul[f'{atm.prefix}iegrid'].data)
self.ndp[atm.name].register(name='ldint@energy', associated_axes=None, grid=hdul[f'{atm.prefix}ipgrid'].data)

if f'{atm.name}:ext' in self.content:
# associated axes:
ebvs = np.array(list(hdul[f'{atm.prefix}_ebvs'].data['ebvs']))
rvs = np.array(list(hdul[f'{atm.prefix}_rvs'].data['rvs']))

self.ndp[atm.name].register('ext@photon', (ebvs, rvs), hdul[f'{atm.prefix}XEGRID'].data)
self.ndp[atm.name].register('ext@energy', (ebvs, rvs), hdul[f'{atm.prefix}XPGRID'].data)
self.ndp[atm.name].register(name='ext@photon', associated_axes=(ebvs, rvs), grid=hdul[f'{atm.prefix}xegrid'].data)
self.ndp[atm.name].register(name='ext@energy', associated_axes=(ebvs, rvs), grid=hdul[f'{atm.prefix}xpgrid'].data)

return self

Expand Down Expand Up @@ -2519,52 +2525,6 @@ def get_passband(passband, content=None, reload=False, update_if_necessary=False

return _pbtable[passband]['pb']

def Inorm_bol_bb(Teff=5772., logg=4.43, abun=0.0, atm='blackbody', intens_weighting='photon'):
r"""
Computes normal bolometric intensity using the Stefan-Boltzmann law,
Inorm_bol_bb = 1/\pi \sigma T^4. If photon-weighted intensity is
requested, Inorm_bol_bb is multiplied by a conversion factor that
comes from integrating lambda/hc P(lambda) over all lambda.
Input parameters mimick the <phoebe.atmospheres.passbands.Passband.Inorm>
method for calling convenience.
Arguments
------------
* `Teff` (float/array, optional, default=5772): value or array of effective
temperatures.
* `logg` (float/array, optional, default=4.43): IGNORED, for class
compatibility only.
* `abun` (float/array, optional, default=0.0): IGNORED, for class
compatibility only.
* `atm` (string, optional, default='blackbody'): atmosphere model, must be
`'blackbody'`, otherwise exception is raised.
* `intens_weighting`
Returns
---------
* (float/array) float or array (depending on input types) of normal
bolometric blackbody intensities.
Raises
--------
* ValueError: if `atm` is anything other than `'blackbody'`.
"""

if atm != 'blackbody':
raise ValueError('atmosphere must be set to blackbody for Inorm_bol_bb.')

if intens_weighting == 'photon':
factor = 2.6814126821264836e22/Teff
else:
factor = 1.0

# convert scalars to vectors if necessary:
if not hasattr(Teff, '__iter__'):
Teff = np.array((Teff,))

return factor * sigma_sb.value * Teff**4 / np.pi


if __name__ == '__main__':
# This will generate bolometric and Johnson V passband files. Note that
Expand Down
73 changes: 40 additions & 33 deletions phoebe/backend/universe.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@
from math import sqrt, sin, cos, acos, atan2, trunc, pi
import sys, os
import copy
from astropy.constants import sigma_sb

from phoebe.atmospheres import models, passbands
from phoebe.distortions import roche, rotstar
Expand Down Expand Up @@ -335,11 +336,9 @@ def handle_reflection(self, **kwargs):
logger.debug("reflection: computing bolometric intensities")
fluxes_intrins_per_body = []
for starref, body in self.items():
if body.mesh is None: continue
abs_normal_intensities = passbands.Inorm_bol_bb(Teff=body.mesh.teffs.for_computations,
atm='blackbody',
intens_weighting='energy')

if body.mesh is None:
continue
abs_normal_intensities = sigma_sb * body.mesh.teffs.for_computations**4 / np.pi # bolometric intensities
fluxes_intrins_per_body.append(abs_normal_intensities * np.pi)

fluxes_intrins_flat = meshes.pack_column_flat(fluxes_intrins_per_body)
Expand All @@ -356,48 +355,47 @@ def handle_reflection(self, **kwargs):
normals_per_body = list(meshes.get_column('vnormals').values())
areas_per_body = list(meshes.get_column('areas').values())
irrad_frac_refl_per_body = list(meshes.get_column('irrad_frac_refl', computed_type='for_computations').values())
teffs_intrins_per_body = list(meshes.get_column('teffs', computed_type='for_computations').values())

ld_func_and_coeffs = [tuple([_bytes(body.ld_func['bol'])] + [np.asarray(body.ld_coeffs['bol'])]) for body in self.bodies]
logger.debug("irradiation ld_func_and_coeffs: {}".format(ld_func_and_coeffs))
fluxes_intrins_and_refl_per_body = libphoebe.mesh_radiosity_problem_nbody_convex(vertices_per_body,
triangles_per_body,
normals_per_body,
areas_per_body,
irrad_frac_refl_per_body,
fluxes_intrins_per_body,
ld_func_and_coeffs,
_bytes(self.irrad_method.title()),
support=_bytes('vertices')
)
fluxes_intrins_and_refl_per_body = libphoebe.mesh_radiosity_problem_nbody_convex(
vertices_per_body,
triangles_per_body,
normals_per_body,
areas_per_body,
irrad_frac_refl_per_body,
fluxes_intrins_per_body,
ld_func_and_coeffs,
_bytes(self.irrad_method.title()),
support=_bytes('vertices')
)

fluxes_intrins_and_refl_flat = meshes.pack_column_flat(fluxes_intrins_and_refl_per_body)

else:
logger.debug("handling reflection (general case), method='{}'".format(self.irrad_method))

vertices_flat = meshes.get_column_flat('vertices') # np.ndarray
triangles_flat = meshes.get_column_flat('triangles') # np.ndarray
normals_flat = meshes.get_column_flat('vnormals') # np.ndarray
areas_flat = meshes.get_column_flat('areas') # np.ndarray
vertices_flat = meshes.get_column_flat('vertices') # np.ndarray
triangles_flat = meshes.get_column_flat('triangles') # np.ndarray
normals_flat = meshes.get_column_flat('vnormals') # np.ndarray
areas_flat = meshes.get_column_flat('areas') # np.ndarray
irrad_frac_refl_flat = meshes.get_column_flat('irrad_frac_refl', computed_type='for_computations') # np.ndarray

ld_func_and_coeffs = [tuple([_bytes(body.ld_func['bol'])] + [np.asarray(body.ld_coeffs['bol'])]) for body in self.mesh_bodies] # list
ld_inds_flat = meshes.pack_column_flat({body.comp_no: np.full(fluxes.shape, body.comp_no-1) for body, fluxes in zip(self.mesh_bodies, fluxes_intrins_per_body)}) # np.ndarray

fluxes_intrins_and_refl_flat = libphoebe.mesh_radiosity_problem(vertices_flat,
triangles_flat,
normals_flat,
areas_flat,
irrad_frac_refl_flat,
fluxes_intrins_flat,
ld_func_and_coeffs,
ld_inds_flat,
_bytes(self.irrad_method.title()),
support=_bytes('vertices')
)


fluxes_intrins_and_refl_flat = libphoebe.mesh_radiosity_problem(
vertices_flat,
triangles_flat,
normals_flat,
areas_flat,
irrad_frac_refl_flat,
fluxes_intrins_flat,
ld_func_and_coeffs,
ld_inds_flat,
_bytes(self.irrad_method.title()),
support=_bytes('vertices')
)

teffs_intrins_flat = meshes.get_column_flat('teffs', computed_type='for_computations')

Expand Down Expand Up @@ -1897,6 +1895,15 @@ def _populate_lc(self, dataset, ignore_effects=False, **kwargs):
blending_method=blending_method
)

# print(f'universe: {pb.pbname=}')
# print(f'universe: {atm_model=}')
# print(f'universe: {ldatm_model=}')
# print(f'universe: {extinct=}')
# print(f'universe: {ignore_effects=}')
# print(f'universe: {query_table=}')
# print(f'universe: {query_pts.shape=}')
# print(f'universe: {abs_intensities=}')

# Beaming/boosting
if boosting_method == 'none' or ignore_effects:
boost_factors = 1.0
Expand Down
1 change: 1 addition & 0 deletions phoebe/frontend/bundle.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10950,6 +10950,7 @@ def compute_pblums(self, compute=None, model=None, pblum=True, pblum_abs=False,
if system is not None:
system.reset(force_recompute_instantaneous=True)

# print(f'{pblums_abs=} {pblums_scale=} {pblums_rel=} {pbfluxes=}')
return system, pblums_abs, pblums_scale, pblums_rel, pbfluxes

# users will see the twig dictionaries with the exposed values based on
Expand Down

0 comments on commit 7ffe0f6

Please sign in to comment.