flip constraint for overcontact binary

[minnieliao]

I want to simulate the overcontact binary system with parameters T_p = 5000K, T_s = 6000K, q = 0.35, inclination = 90 degrees, orbital period = 1 day, and fill-out factor = 0.2. However, I can't figure out how to fly the constraint to solve for the fill-out factor and run_compute.

here is my code

cb =phoebe.default_binary(contact_binary =True)

cb['teff@primary'] =5000
cb['teff@secondary'] =6000

cb['q@binary'] =0.35

cb['incl'] =90

cb['period@orbit'] =1.0

cb.flip_constraint('pot@contact_envelope@envelope@component', solve_for ='fillout_factor@component')
cb['fillout_factor@component'] = 0.2

cb.add_dataset('lc', compute_times =phoebe.linspace(0, 1.0, 101), dataset ='lc01', overwrite =True)

cb.add_dataset('mesh', compute_times=[0], dataset='mesh01', overwrite =True)

cb.run_compute()

I think there are some thing wrong for cb.flip_constraint('pot@contact_envelope@envelope@component', solve_for ='fillout_factor@component')

because it shows thee error like

ValueError: 0 results found for twig: 'pot@contact_envelope@envelope@component', {'context': 'constraint', 'check_visible': True, 'check_default': True, 'check_advanced': False, 'check_single': False}

Thanks!

[kecnry]

There are a few things going on here that make flipping this constraint complicated:

1. The error message isn't the most helpful here, but the reason it isn't returning any results is because you're passing context='component' (via the @component at the end) while flip constraint is expecting the constraint itself in context='constraint'. If you just leave of @component then that error message will be replaced by KeyError: 'fillout_factor@contact_envelope@component was not found in expression', which leads us to....

2. the constraint you're trying to flip is actually between the potential and equivalent radius, and is not able to be flipped directly for fillout-factor. For contact binaries, there are 5 intertwined parameters, only 2 of which can be free, which makes flipping them more complicated than usual. If you call print(b.filter(qualifier=['pot', 'fillout_factor', 'requiv', 'q'])) (directly after creating the bundle), you'll see those 5 parameters and 3 constraint expressions:

ParameterSet: 8 parameters
         requiv@primary@component: 1.5 solRad
C      requiv@secondary@component: 1.5000000000000122 solRad
               q@binary@component: 1.0
C  fillout_factor@contact_enve...: 0.6417897080770861
C  pot@contact_envelope@component: 3.4013774072298815
                   pot@constraint: requiv_to_pot_contact({requiv@primary@component}, {q@binary@component}, {sma@binary@component}, 1)
                requiv@constraint: pot_to_requiv_contact({pot@contact_envelope@component}, {q@binary@component}, {sma@binary@component}, 2)
        fillout_factor@constraint: pot_to_fillout_factor({q@binary@component}, {pot@contact_envelope@component})

this is showing that, by default, the equivalent radius of the primary star and the mass-ratio are the only parameters that are "free", with the equivalent radius of the secondary, the potential, and fillout-factor all being "constrained" (noted by the "C" in front) via the listed expressions in the constraint parameters. In order to make fillout-factor free instead, you'll ultimately need to flip the fillout_factor constraint. You could flip fillout_factor to solve for the mass-ratio, q, and then have fillout-factor and requiv as free parameters, but since you also want q to be free, that means you'll need to flip multiple constraints so that requiv@primary is constrained in its place (since flipping fillout_factor for pot will complain that pot is already constrained).

import phoebe

b = phoebe.default_binary(contact_binary=True)
b.flip_constraint('pot@contact_envelope', solve_for='requiv@primary')
# print(b.filter(qualifier=['pot', 'fillout_factor', 'requiv', 'q']))
b.flip_constraint('fillout_factor@contact_envelope', solve_for='pot@contact_envelope')
print(b.filter(qualifier=['pot', 'fillout_factor', 'requiv', 'q']))

which now gives you

ParameterSet: 8 parameters
C        requiv@primary@component: 1.5000000000000033 solRad
C      requiv@secondary@component: 1.5000000000000122 solRad
               q@binary@component: 1.0
   fillout_factor@contact_enve...: 0.6417897080770861
C  pot@contact_envelope@component: 3.4013774072298815
        requiv@primary@constraint: pot_to_requiv_contact({pot@contact_envelope@component}, {q@binary@component}, {sma@binary@component}, 1)
      requiv@secondary@constraint: pot_to_requiv_contact({pot@contact_envelope@component}, {q@binary@component}, {sma@binary@component}, 2)
                   pot@constraint: fillout_factor_to_pot({q@binary@component}, {fillout_factor@contact_envelope@component})

Hope this helps and makes a little more sense!