Coverage for src/galsbi/ucat_sps/util_sps.py: 100%
164 statements
« prev ^ index » next coverage.py v7.17.2, created at 2026-11-09 18:05 +0000
1# Copyright (C) 2025 LMU Munich
2# Author: Luca Tortorelli
3# created: Sep 2026
5"""General numerical and model-evaluation utilities for GalSBI-SPS."""
7import h5py
8import numpy as np
9from promage import ProMage
10from promage.exceptions import ManifestError
11from scipy import interpolate
12from scipy.integrate import trapezoid
14from galsbi.ucat.spectrum_util import extinction_coefficient, spline_ext_coeff
16GALAXY_TYPE_SEED_IDS = {"blue": 1, "red": 2}
17PROPERTY_SEED_IDS = {
18 "stellar_mass_function": 1,
19 "star_formation_history": 2,
20 "gas_metallicity": 3,
21 "gas_ionization": 4,
22 "dust_attenuation": 5,
23 "dust_emission": 6,
24 "size": 7,
25 "ellipticity": 8,
26 "sersic_index": 9,
27 "velocity_dispersion": 10,
28 "agn": 11,
29 "position": 12,
30}
31MASS_COMPONENT_SEED_IDS = {None: 0, "lowmass": 1, "highmass": 2}
34def galaxy_property_seed(par, pixel, galaxy_type, property_name, mass_component=None):
35 """Return a reproducible and independent seed for one sampled property.
37 The seed depends on the catalogue seed, HEALPix pixel, galaxy population,
38 physical property and, for the stellar mass function, Schechter component.
39 It does not depend on the position of a pixel in the input pixel list, so
40 reordering that list does not alter the realization assigned to each pixel.
42 :param par: (obj) GalSBI parameter object containing the catalogue seed.
43 :param pixel: (int) HEALPix pixel being sampled.
44 :param galaxy_type: (str) ``"red"`` or ``"blue"`` galaxy population.
45 :param property_name: (str) name of the sampled physical property.
46 :param mass_component: (str or None) ``"lowmass"`` or ``"highmass"`` for
47 the corresponding Schechter component; otherwise ``None``.
48 :return seed: (int) deterministic 32-bit seed for the requested random stream.
49 """
50 seed_sequence = np.random.SeedSequence(
51 [
52 int(par.seed),
53 int(pixel),
54 GALAXY_TYPE_SEED_IDS[galaxy_type],
55 PROPERTY_SEED_IDS[property_name],
56 MASS_COMPONENT_SEED_IDS[mass_component],
57 ]
58 )
59 return int(seed_sequence.generate_state(1, dtype=np.uint32)[0])
62def stellar_mass_function_seed(par, pixel, galaxy_type, mass_component):
63 """Return a distinct seed for each stellar-mass-function component.
65 The red and blue populations each contain low- and high-mass Schechter
66 components. These four components must use different random streams so
67 their Poisson counts and sampled masses are statistically independent.
69 :param par: (obj) GalSBI parameter object containing the catalogue seed.
70 :param pixel: (int) HEALPix pixel being sampled.
71 :param galaxy_type: (str) ``"red"`` or ``"blue"`` galaxy population.
72 :param mass_component: (str) ``"lowmass"`` or ``"highmass"`` Schechter
73 component.
74 :return seed: (int) deterministic seed for the requested population and
75 Schechter component.
76 """
77 return galaxy_property_seed(
78 par,
79 pixel,
80 galaxy_type,
81 "stellar_mass_function",
82 mass_component=mass_component,
83 )
86def get_age_of_the_universe_from_redshift(redshifts, cosmology):
87 """Return the age of the Universe at each input redshift.
89 The calculation uses the flat matter-plus-cosmological-constant solution
90 and the parameters stored in the configured PyCosmo object.
92 :param redshifts: (float or array_like[n_gal,]) galaxy redshifts.
93 :param cosmology: (obj) configured PyCosmo cosmology object.
94 :return age: (float or array_like[n_gal,]) age of the Universe in years.
95 """
97 redshifts = np.asarray(redshifts, dtype=np.float64)
98 H0 = cosmology.params.H0
99 omega_m = cosmology.params.omega_m
100 omega_l = cosmology.params.omega_l
101 age = np.sqrt(omega_l * (1 + redshifts) ** (-3)) + np.sqrt(
102 omega_m + omega_l * (1 + redshifts) ** (-3)
103 )
104 age /= np.sqrt(omega_m)
105 age = np.log(age)
106 age *= 2 / (3 * H0 * np.sqrt(omega_l))
107 age *= 3.0857e19
108 age /= 3.1558150e7
109 return age
112def get_lookbacktime_from_redshift(redshifts, cosmology):
113 """Return the observer lookback time at each input redshift.
115 :param redshifts: (float or array_like[n_gal,]) galaxy redshifts.
116 :param cosmology: (obj) configured PyCosmo cosmology object.
117 :return lookback_time: (float or array_like[n_gal,]) lookback time in years.
118 """
120 age_at_redshift_zero = get_age_of_the_universe_from_redshift(0, cosmology)
121 return age_at_redshift_zero - get_age_of_the_universe_from_redshift(
122 redshifts, cosmology
123 )
126def _normalise_filter_names(filter_names):
127 """Return configured filter names as an ordered list.
129 Some GalSBI configurations store ProSpect filter names as a list, while
130 older configurations use a dictionary whose values contain the full
131 ProSpect names. This helper provides one representation without changing
132 the order supplied by the configuration.
134 :param filter_names: (sequence[str] or dict) configured ProSpect/ProMage
135 filter names.
136 :return normalised_filter_names: (list[str]) ordered full filter names.
137 """
139 if isinstance(filter_names, dict):
140 return list(filter_names.values())
141 return list(filter_names)
144def _local_filter_names(par, filter_names):
145 """Translate full ProSpect filter names to GalSBI catalogue band names.
147 The explicit ``filters_prospect_to_galsbi`` mapping is used whenever it is
148 available. The prefix before the first underscore is retained only as a
149 compatibility fallback for older configurations.
151 :param par: (obj) GalSBI-SPS parameter object containing the filter maps.
152 :param filter_names: (sequence[str]) full ProSpect/ProMage filter names.
153 :return local_filter_names: (list[str]) corresponding GalSBI band names.
154 """
156 prospect_to_galsbi = getattr(par, "filters_prospect_to_galsbi", {})
157 filter_full_names = getattr(par, "filters_full_names", {})
158 local_filter_names = []
159 for filter_name in filter_names:
160 candidate = prospect_to_galsbi.get(
161 filter_name, str(filter_name).split("_", maxsplit=1)[0]
162 )
163 if candidate in filter_full_names:
164 local_filter_names.append(candidate)
165 elif candidate.lower() in filter_full_names:
166 local_filter_names.append(candidate.lower())
167 else:
168 local_filter_names.append(candidate)
169 return local_filter_names
172def _h5_filter_name(par, h5_file, local_filter_name, prospect_filter_name):
173 """Find the HDF5 group containing one configured filter transmission.
175 GalSBI filter files may identify a band by the value in
176 ``par.filters_full_names``, by its full ProSpect name, or by the short
177 GalSBI name. The candidates are checked in that order.
179 :param par: (obj) GalSBI-SPS parameter object containing filter names.
180 :param h5_file: (h5py.File) open HDF5 file containing filter curves.
181 :param local_filter_name: (str) short GalSBI band name.
182 :param prospect_filter_name: (str) full ProSpect/ProMage filter name.
183 :return h5_filter_name: (str) matching HDF5 group name.
184 :raises KeyError: if none of the configured names is present in the file.
185 """
187 filter_full_names = getattr(par, "filters_full_names", {})
188 if (
189 isinstance(filter_full_names, dict)
190 and local_filter_name in filter_full_names
191 and filter_full_names[local_filter_name] in h5_file
192 ):
193 return filter_full_names[local_filter_name]
194 if prospect_filter_name in h5_file:
195 return prospect_filter_name
196 if local_filter_name in h5_file:
197 return local_filter_name
198 raise KeyError(
199 f"Could not find filter {prospect_filter_name!r} in {par.filters_file_name}"
200 )
203def _filter_effective_wavelengths(par, filter_names):
204 """Compute transmission-weighted effective wavelengths for selected bands.
206 Each value is evaluated as the integral of wavelength times transmission
207 divided by the integral of transmission, using the curves in the filter
208 HDF5 file configured by ``par.filters_file_name``.
210 :param par: (obj) GalSBI-SPS parameter object containing filter paths and
211 name mappings.
212 :param filter_names: (sequence[str]) full ProSpect/ProMage filter names.
213 :return effective_wavelengths: (array_like[n_bands,]) effective wavelengths
214 in Angstrom.
215 """
217 local_filter_names = _local_filter_names(par, filter_names)
218 effective_wavelengths = np.empty(len(filter_names), dtype=float)
219 with h5py.File(par.filters_file_name, "r") as h5_file:
220 for index, (local_filter_name, prospect_filter_name) in enumerate(
221 zip(local_filter_names, filter_names, strict=True)
222 ):
223 h5_filter_name = _h5_filter_name(
224 par, h5_file, local_filter_name, prospect_filter_name
225 )
226 wavelength = h5_file[h5_filter_name]["lam"][:]
227 transmission = h5_file[h5_filter_name]["amp"][:]
228 effective_wavelengths[index] = trapezoid(
229 wavelength * transmission, wavelength
230 ) / trapezoid(transmission, wavelength)
231 return effective_wavelengths
234def apply_extinction_to_magnitudes(
235 magnitudes, excess_b_v, par, filter_names=None, extinction_spline=None
236):
237 """Apply Milky Way extinction to observed-frame band magnitudes.
239 :param magnitudes: (array_like[n_gal, n_bands]) observed-frame magnitudes.
240 :param excess_b_v: (array_like[n_gal,]) Milky Way E(B-V) values.
241 :param par: (obj) GalSBI-SPS parameter object containing filter data.
242 :param filter_names: (sequence[str], optional) ProSpect/ProMage filter names;
243 defaults to ``par.filters_full_names_prospect``.
244 :param extinction_spline: (obj, optional) extinction-law interpolation
245 spline. It is constructed when omitted.
246 :return extincted_magnitudes: (array_like[n_gal, n_bands]) observed-frame
247 magnitudes including Milky Way extinction.
248 """
250 if extinction_spline is None:
251 extinction_spline = spline_ext_coeff()
252 filter_names = _normalise_filter_names(
253 par.filters_full_names_prospect if filter_names is None else filter_names
254 )
255 effective_wavelengths = _filter_effective_wavelengths(par, filter_names)
256 extinction = extinction_coefficient(
257 effective_wavelengths / 10000, np.asarray(excess_b_v), extinction_spline
258 )
259 return np.asarray(magnitudes) + extinction
262def reorder_prospect_magnitudes_to_filters(par, magnitudes, filter_names=None):
263 """Reorder ProSpect or ProMage magnitudes to match ``par.filters``.
265 :param par: (obj) GalSBI-SPS parameter object containing the explicit
266 ProSpect-to-GalSBI filter mapping.
267 :param magnitudes: (array_like[n_gal, n_bands]) magnitudes in the supplied
268 ProSpect filter order.
269 :param filter_names: (sequence[str], optional) order of the magnitude
270 columns; defaults to ``par.filters_full_names_prospect``.
271 :return reordered_magnitudes: (array_like[n_gal, n_bands]) magnitudes in
272 ``par.filters`` order.
273 """
275 filter_names = _normalise_filter_names(
276 par.filters_full_names_prospect if filter_names is None else filter_names
277 )
278 column_index = {
279 par.filters_prospect_to_galsbi[filter_name]: index
280 for index, filter_name in enumerate(filter_names)
281 }
282 reorder = [column_index[filter_name] for filter_name in par.filters]
283 return np.asarray(magnitudes)[:, reorder]
286def validate_recent_burst_properties(has_burst, fburst, n_galaxies):
287 """Prepare recent-burst arrays for magnitude and mass calculations.
289 Omitting both inputs represents a purely smooth population, preserving the
290 magnitude-function interface used before the burst model was introduced.
291 The population sampler is responsible for enforcing the physical relation
292 between the burst indicator and burst formed-mass fraction.
294 :param has_burst: (array_like[n_gal,] or None) recent-burst indicators.
295 :param fburst: (array_like[n_gal,] or None) recent-burst formed-mass
296 fractions.
297 :param n_galaxies: (int) expected number of galaxies.
298 :return has_burst: (array_like[n_gal,]) Boolean recent-burst indicators.
299 :return fburst: (array_like[n_gal,]) recent-burst formed-mass fractions.
300 """
302 if has_burst is None and fburst is None:
303 return (
304 np.zeros(n_galaxies, dtype=bool),
305 np.zeros(n_galaxies, dtype=np.float64),
306 )
307 if has_burst is None or fburst is None:
308 raise ValueError("has_burst and fburst must be supplied together.")
309 has_burst = np.asarray(has_burst, dtype=bool)
310 fburst = np.asarray(fburst, dtype=np.float64)
311 if has_burst.shape != (n_galaxies,) or fburst.shape != (n_galaxies,):
312 raise ValueError("has_burst and fburst must contain one value per galaxy.")
313 return has_burst, fburst
316def build_promage_physical_properties(
317 z,
318 mSFR,
319 mpeak,
320 mperiod,
321 mskew,
322 Zgas_final,
323 logU,
324 dust_attenuation_params,
325 dust_emission_params,
326):
327 """Build the smooth-SFH physical properties required by ProMage.
329 :param z: (array_like[n_gal,]) galaxy redshifts.
330 :param mSFR: (array_like[n_gal,]) smooth SFH amplitudes in Msun/yr.
331 :param mpeak: (array_like[n_gal,]) smooth SFH peak times in Gyr.
332 :param mperiod: (array_like[n_gal,]) smooth SFH durations in Gyr.
333 :param mskew: (array_like[n_gal,]) dimensionless smooth SFH skew values.
334 :param Zgas_final: (array_like[n_gal,]) gas metallicities at observation.
335 :param logU: (array_like[n_gal,]) base-10 gas ionisation parameters.
336 :param dust_attenuation_params: (array_like[n_gal, 3]) ProSpect
337 ``logtaubirth``, ``logtauscreen``, and ``Eb`` parameters.
338 :param dust_emission_params: (array_like[n_gal, 2]) ProSpect
339 ``alphabirth`` and ``alphascreen`` parameters.
340 :return physical_properties: (dict[str, array_like[n_gal,]]) properties in
341 the order used to train the smooth ProMage model.
342 """
344 return {
345 "z": np.asarray(z),
346 "logmSFR": np.log10(mSFR),
347 "mpeak": np.asarray(mpeak),
348 "logmperiod": np.log10(mperiod),
349 "mskew": np.asarray(mskew),
350 "logZfinal": np.log10(Zgas_final),
351 "logtaubirth": np.asarray(dust_attenuation_params)[:, 0],
352 "logtauscreen": np.asarray(dust_attenuation_params)[:, 1],
353 "alphabirth": np.asarray(dust_emission_params)[:, 0],
354 "alphascreen": np.asarray(dust_emission_params)[:, 1],
355 "logU": np.asarray(logU),
356 }
359def load_promage_models(par, include_recent_bursts=False):
360 """Load the smooth ProMage model and the recent-burst model when needed.
362 :param par: (obj) GalSBI-SPS parameter object.
363 :param include_recent_bursts: (bool) load the model trained with
364 ``massfunc_snorm_burst_trunc`` in addition to the smooth model.
365 :return smooth_model: (ProMage) magnitude model for smooth SFHs.
366 :return burst_model: (ProMage or None) magnitude model for recent-burst
367 SFHs.
368 """
370 from galsbi.ucat_sps.config.common_sps import (
371 PROMAGE_BURST_PROPERTIES,
372 PROMAGE_BURST_SFH_MODEL,
373 PROMAGE_SMOOTH_PROPERTIES,
374 PROMAGE_SMOOTH_SFH_MODEL,
375 PROSPECT_SNORM_TRUNCATION,
376 )
378 if not np.isclose(
379 getattr(par, "mtrunc", PROSPECT_SNORM_TRUNCATION),
380 PROSPECT_SNORM_TRUNCATION,
381 ):
382 raise ValueError("The ProMage models require mtrunc=2.")
383 smooth_sfh_model = getattr(par, "promage_sfh_model", PROMAGE_SMOOTH_SFH_MODEL)
384 if smooth_sfh_model != PROMAGE_SMOOTH_SFH_MODEL:
385 raise ValueError(
386 "promage_sfh_model must be 'massfunc_snorm_trunc' for GalSBI-SPS."
387 )
388 requested_models = [smooth_sfh_model]
389 if include_recent_bursts:
390 burst_sfh_model = getattr(
391 par, "promage_burst_sfh_model", PROMAGE_BURST_SFH_MODEL
392 )
393 if burst_sfh_model != PROMAGE_BURST_SFH_MODEL:
394 raise ValueError(
395 "promage_burst_sfh_model must be "
396 "'massfunc_snorm_burst_trunc' for recent-burst galaxies."
397 )
398 requested_models.append(burst_sfh_model)
400 loaded_models = []
401 for sfh_model in requested_models:
402 try:
403 model = ProMage(
404 resource_dir=par.ProMage_model_filepath,
405 sfh_model=sfh_model,
406 device="cpu",
407 on_out_of_range="raise",
408 lazy=True,
409 )
410 except ManifestError as error:
411 if "legacy ProMage manifest" not in str(error):
412 raise
413 raise RuntimeError(
414 f"GalSBI-SPS requested the {sfh_model!r} ProMage model, but "
415 f"{par.ProMage_model_filepath!r} contains old ProMage resources "
416 "that do not identify their training SFH. Refresh the cached "
417 "ProMage_res data."
418 ) from error
419 expected_properties = (
420 PROMAGE_BURST_PROPERTIES
421 if sfh_model == PROMAGE_BURST_SFH_MODEL
422 else PROMAGE_SMOOTH_PROPERTIES
423 )
424 if model.sfh_model != sfh_model:
425 raise ValueError(
426 f"The loaded ProMage model describes {model.sfh_model}, not "
427 f"{sfh_model}."
428 )
429 if tuple(model.required_properties) != expected_properties:
430 raise ValueError(
431 f"The {sfh_model} ProMage model expects "
432 f"{tuple(model.required_properties)}, not {expected_properties}."
433 )
434 loaded_models.append(model)
435 return loaded_models[0], loaded_models[1] if include_recent_bursts else None
438def predict_promage_magnitudes(
439 physical_properties,
440 has_burst,
441 fburst,
442 bands,
443 frame,
444 role,
445 smooth_model,
446 burst_model=None,
447):
448 """Predict ProMage magnitudes and restore the original galaxy order.
450 :param physical_properties: (dict[str, array_like[n_gal,]]) common smooth
451 and burst ProMage inputs.
452 :param has_burst: (array_like[n_gal,]) recent-burst indicators.
453 :param fburst: (array_like[n_gal,]) recent-burst formed-mass fractions.
454 :param bands: (sequence[str]) requested ProSpect filter names.
455 :param frame: (str) ``observed`` or ``absolute`` magnitude frame.
456 :param role: (str) ``selection`` or ``magnitude`` model purpose.
457 :param smooth_model: (ProMage) model trained with smooth SFHs.
458 :param burst_model: (ProMage, optional) model trained with recent bursts.
459 :return magnitudes: (array_like[n_gal, n_bands]) predicted magnitudes in
460 input galaxy order.
461 """
463 n_galaxies = np.asarray(physical_properties["z"]).size
464 has_burst, fburst = validate_recent_burst_properties(has_burst, fburst, n_galaxies)
465 bands = list(bands)
466 magnitudes = np.empty((n_galaxies, len(bands)), dtype=np.float64)
467 for model, selected in ((smooth_model, ~has_burst), (burst_model, has_burst)):
468 if not np.any(selected):
469 continue
470 if model is None:
471 raise ValueError("A burst ProMage model is required for recent bursts.")
472 selected_properties = {
473 name: np.asarray(values)[selected]
474 for name, values in physical_properties.items()
475 }
476 if model is burst_model:
477 selected_properties["logfburst"] = np.log10(fburst[selected])
478 prediction = np.asarray(
479 model.predict(
480 properties=selected_properties,
481 bands=bands,
482 frame=frame,
483 role=role,
484 return_format="array",
485 )
486 )
487 expected_shape = (int(np.count_nonzero(selected)), len(bands))
488 if prediction.shape != expected_shape:
489 raise ValueError(
490 f"The {model.sfh_model} ProMage model returned shape "
491 f"{prediction.shape}; expected {expected_shape}."
492 )
493 magnitudes[selected] = prediction
494 return magnitudes
497def apply_galaxy_observed_magnitude_cut(par, galaxy_properties):
498 """Select galaxies brighter than the ProMage observed-magnitude limit.
500 Smooth and recent-burst galaxies are evaluated with their corresponding
501 ProMage selection models. Their predictions are restored to catalogue
502 order before Milky Way extinction and the common magnitude cut are
503 applied. The configured limit must be supported by every SFH family used
504 by the selected population model.
506 :param par: (obj) object containing the GalSBI-SPS parameters.
507 :param galaxy_properties: (obj) galaxy catalogue containing the ProMage
508 physical properties, ``has_burst``, ``fburst``, and ``excess_b_v``.
509 :return selected_objects: (array_like[n_gal,]) Boolean mask selecting
510 objects that pass the magnitude cut.
511 """
513 from galsbi.ucat_sps.config.common_sps import CONDITIONAL_SNORM_BURST_TRUNC
515 physical_properties = build_promage_physical_properties(
516 galaxy_properties.z,
517 galaxy_properties.mSFR,
518 galaxy_properties.mpeak,
519 galaxy_properties.mperiod,
520 galaxy_properties.mskew,
521 galaxy_properties.gas_metallicity,
522 galaxy_properties.gas_ionization,
523 galaxy_properties.dust_attenuation,
524 galaxy_properties.dust_emission,
525 )
526 has_burst, fburst = validate_recent_burst_properties(
527 getattr(galaxy_properties, "has_burst", None),
528 getattr(galaxy_properties, "fburst", None),
529 len(galaxy_properties.z),
530 )
531 include_recent_bursts = (
532 np.any(has_burst)
533 or getattr(par, "sfh_population_model", None) == CONDITIONAL_SNORM_BURST_TRUNC
534 )
535 smooth_model, burst_model = load_promage_models(
536 par, include_recent_bursts=include_recent_bursts
537 )
538 selection_filter = par.filters_full_names_prospect_selection
539 selection_limits = {
540 model.sfh_model: model.selection_threshold(selection_filter, frame="observed")
541 for model in (smooth_model, burst_model)
542 if model is not None
543 }
544 faintest_supported_limit = min(selection_limits.values())
545 if par.promage_selection_mag_threshold > faintest_supported_limit:
546 raise ValueError(
547 "promage_selection_mag_threshold cannot exceed the ProMage "
548 f"selection limit ({faintest_supported_limit}) for {selection_filter}; "
549 f"the SFH-family limits are {selection_limits}."
550 )
552 magnitude_selection = predict_promage_magnitudes(
553 physical_properties,
554 has_burst,
555 fburst,
556 bands=[selection_filter],
557 frame="observed",
558 role="selection",
559 smooth_model=smooth_model,
560 burst_model=burst_model,
561 )[:, 0]
562 magnitude_selection = apply_extinction_to_magnitudes(
563 np.asarray(magnitude_selection)[:, np.newaxis],
564 galaxy_properties.excess_b_v,
565 par,
566 filter_names=[selection_filter],
567 )[:, 0]
568 return magnitude_selection < par.promage_selection_mag_threshold
571def compute_ab_magnitudes(par, observed_wavelengths_A, observed_fluxes_erg_s_cm2_A):
572 """Integrate observed galaxy SEDs to obtain AB magnitudes.
574 :param par: (obj) GalSBI-SPS parameter object containing filter curves.
575 :param observed_wavelengths_A: (array_like[n_gal, n_lambda]) observed-frame
576 wavelengths in Angstrom.
577 :param observed_fluxes_erg_s_cm2_A: (array_like[n_gal, n_lambda]) observed
578 flux densities in erg/s/cm2/Angstrom.
579 :return observed_magnitudes: (array_like[n_gal, n_bands]) observed-frame
580 AB magnitudes in ``par.filters`` order.
581 """
583 from galsbi.ucat_sps.config.common_sps import LIGHTSPEED_ANGSTROM_S
585 observed_magnitudes = np.zeros((len(observed_fluxes_erg_s_cm2_A), len(par.filters)))
586 with h5py.File(par.filters_file_name, "r") as h5_file:
587 for band_index, filter_name in enumerate(par.filters):
588 filter_wavelength = h5_file[par.filters_full_names[filter_name]]["lam"][:]
589 filter_transmission = h5_file[par.filters_full_names[filter_name]]["amp"][:]
590 interpolate_transmission = interpolate.interp1d(
591 filter_wavelength,
592 filter_transmission,
593 bounds_error=False,
594 fill_value=0,
595 )
596 for galaxy_index in range(len(observed_fluxes_erg_s_cm2_A)):
597 transmission = interpolate_transmission(
598 observed_wavelengths_A[galaxy_index]
599 )
600 numerator = trapezoid(
601 observed_fluxes_erg_s_cm2_A[galaxy_index]
602 * transmission
603 * observed_wavelengths_A[galaxy_index]
604 / LIGHTSPEED_ANGSTROM_S,
605 x=observed_wavelengths_A[galaxy_index],
606 )
607 denominator = trapezoid(
608 transmission / observed_wavelengths_A[galaxy_index],
609 x=observed_wavelengths_A[galaxy_index],
610 )
611 observed_magnitudes[galaxy_index, band_index] = (
612 -2.5 * np.log10(numerator / denominator) - 48.6
613 )
614 return observed_magnitudes