Constraint problem on spectroscopic parameters

[mpogacea]

Hello developers,

I have a question concerning spectroscopic binaries (SBs) alone. I would like to create a bundle of spectroscopic parameters only. In particular, I am trying to introduce the semi-amplitude parameters K1 and K2. I am doing as follows:

import phoebe
from phoebe import u 

b = phoebe.default_binary()

K1_def = phoebe.parameters.FloatParameter(qualifier='K1', default_unit= u.km / u.s, value=1, description='semi-amplitude of the primary star')
K1_param, created = b.get_or_create('K1', K1_def, context='component', component='primary')

### Parameters constraining K1
a1_sini = b.get_parameter(qualifier='asini', component='primary', context = 'component')
orbital_period = b.get_parameter(qualifier='period', component='binary', context = 'component')
q = b.get_parameter(qualifier='q', component='binary', context = 'component')
ecc = b.get_parameter(qualifier='ecc', component='binary', context = 'component')

### Defining the constraint for K1
lhs = b.get_parameter(qualifier='K1', component='primary')
rhs = ( 2*np.pi * a1_sini ) / ( orbital_period * (1 + q) * ( 1 - ecc**2 )**0.5 ) 

b.add_constraint(lhs, rhs)

print(b['K1@primary@custom@constraint'])

This is what I get in return, accompanied by an error message:

Constrains (qualifier): K1
Expression in solar units (value): (6.28318530717958623e+00 * {asini@primary@component}) / (({period@binary@component} * (1.00000000000000000e+00 + {q@binary@component})) * ((({ecc@binary@component} ** 2.000000) - 1.00000000000000000e+00) ** 5.00000000000000000e-01))
Current Result (result): nan km / s

/home/matei/venv_phoebe/lib/python3.12/site-packages/astropy/units/quantity.py:658: RuntimeWarning: invalid value encountered in sqrt
  result = super().__array_ufunc__(function, method, *arrays, **kwargs)
<string>:1: RuntimeWarning: invalid value encountered in scalar power
<string>:1: RuntimeWarning: invalid value encountered in scalar power
<string>:1: RuntimeWarning: invalid value encountered in scalar power
<string>:1: RuntimeWarning: invalid value encountered in scalar power

From my understanding the error message arises from the fact that Phoebe is trying to compute the square root of the negative number (1 - ecc^2). And so, my question is, why does the constraint defined as (1 - ecc^2) in my code appear as (ecc^2 - 1) when I ask Phoebe to explicitly write the K1 constraint ?

Thank you in advance for your help !

[kecnry]

Thanks for reaching out! There are some limitations in the constraint math (caused by python's limitations for __sub__ vs __sub__ on builtin objects) that I suspect may be the problem here - if so we can try to revisit some workarounds or at least mention this better in the documentation for advanced use-cases like what you're trying here. Can you try defining as -1*ecc**2 + 1 and see if that gives the results you expect?