Phoenix atmosphere out of bounds error

[erlwez]

I am trying to fit a binary consisting of a MS star and a pulsar. When using phoenix atmosphere for the companion I get a warning during run_compute()

Wed, 23 Oct 2024 09:05 UNIVERSE WARNING Atmosphere parameters out of bounds: ldatm=phoenix, teff=[5319.67964699], logg=[3.69438033], abun=[0.]

and the error

ValueError: Could not compute ldint with ldatm='phoenix'. Try changing atm to a table that covers a sufficient range of values. If necessary, set atm to 'blackbody' and/or ld_mode to 'manual' (in which case coefficients will need to be explicitly provided via ld_coeffs). Enable 'warning' logger to see out-of-bound arrays.

According to the linked paper on the Phoebe website these values should all be within range. I have looked through the datasets manually to make sure neither teff or logg goes out of range.

If I try setting the teff@secondary (MS star) to 10000 K and above there is no issue running the program. I can also NOT set asini@primary (pulsar) before run_compute() and the program runs fine with teff@secondary = 5630 K like I want it to be. In that case the asini@primary is calculated to be 6.15 It-s, instead of 0.47 like it should be.
If I try setting the atmosphere to either blackbody or ck2004 the program runs as expected. Could this be a bug or am I missing something?

[Raindogjones]

Can you post a minimal example that leads to the error? It's not obvious to me why phoenix should give and error but ck2004 not, so a code snippet that leads to the error might help us track it down.

As for asini, this is a constrained parameter so, without flipping the constraint with another parameter, you shouldn't be able to set it.

[erlwez]

Added code below. I am flipping for asini and have had the same error when setting asini@primary and if i set the calculated sma@binary instead.

import phoebe
import numpy as np
import matplotlib.pyplot as plt
import math
import astropy.units as u
from datetime import datetime
import os

logger = phoebe.logger('warning')

b = phoebe.default_binary()
period_ = 0.1725
times = np.linspace(0,2*period_,100)  

# datasets
b.add_dataset('lc', times=times, dataset = 'lc01', ld_mode = 'manual', passband='SDSS:g')
b.add_dataset('lc', times=times, dataset = 'lc02', ld_mode = 'manual', passband='SDSS:r')
b.add_dataset('lc', times=times, dataset = 'lc03', ld_mode = 'manual', passband='SDSS:i')
b.add_dataset('mesh', compute_times=period_/4, columns=['teffs', 'intensities@lc01'], ld_mode='manual', 
              dataset = 'mesh01', overwrite = True)
b.add_dataset('orb', compute_times=times, ld_mode='manual', dataset='orb01')
b.add_dataset('rv', compute_times=times, ld_mode='manual', dataset='rv01', passband='SDSS:g')

# flips
b.flip_constraint('asini@primary@constraint', 'sma')
b.flip_constraint('ebv', solve_for='Av')

# system
b.set_value('q@binary', value = 0.33 / 2.27)
b.set_value('incl@binary', value = 63.9)
b.set_value('period@binary', value = period_)
b.set_value('eclipse_method', value = 'only_horizon')
b.set_value('distance', value = 2.9, unit = 'kiloparsec')
b.set_value_all('pblum_mode', value = 'absolute')
b.set_value('ebv', value = 0.38)

# primary
b.set_value('requiv@primary@component', value = 10, unit = 'km')
#b.set_value('teff@primary', value = 3.831e6)
b.set_value('teff@primary', value = 3.0e6)
b.set_value('atm@primary', value = 'blackbody')
b.set_value_all('ld_mode@primary', value = 'manual')
b.set_value('ld_mode_bol@primary', value = 'manual')
#b.set_value('asini@primary@component', value = 0.46814, unit = 'lightsecond')

# secondary
b.set_value('teff@secondary', value = 5630)
b.set_value('requiv@secondary', value = 0.9963*b.get_value('requiv_max@secondary@component'))
b.set_value('gravb_bol@secondary', value = 0.10)
b.set_value('atm@secondary@phoebe01', value = 'phoenix')
b.set_value('irrad_frac_refl_bol@secondary', value = 1)
b.set_value_all('ld_mode@secondary', value = 'lookup')
b.set_value('ld_mode_bol@secondary@component', value = 'lookup')
b.set_value_all('ld_func@secondary',  'quadratic')
b.set_value_all('ld_func_bol@secondary',  'quadratic')

b.run_compute( 'phoebe01', irrad_method = 'wilson', ntriangles = 700, model = 'refl_true', overwrite = True)

[kecnry]

Distortion and/or irradiation may be causing some elements to go out of bounds. Can you enable the warning logger to see which values are out of bounds?

[Raindogjones]

Nope, it looks like it is an interpolation bug. I get the same error running pb.interpolate_ldcoeffs() for most temperatures around 5320 K and phoenix. Paging @aprsa for assistance.

[Raindogjones]

@erlwez The good thing is that you should be safe setting the atmosphere to 'ck2004' instead and going ahead with your fitting (unless you need to sample much lower temperatures for the secondary).

[erlwez]

Thanks for looking into it so quickly. I am trying to reproduce results from a different program so ideally I should be using phoenix, but can make do with ck2004 if it does not work.