phoebe-project / phoebe2

PHOEBE - Eclipsing Binary Star Modeling Software
http://phoebe-project.org
GNU General Public License v3.0
76 stars 28 forks source link

When flipping ecc and per0 to esinw and ecosw, per0 is returned in the wrong quadrant. #905

Open kmhambleton opened 2 weeks ago

kmhambleton commented 2 weeks ago

For example, when omega is set to 170 degrees, after fitting to esinw and ecosw, it is returned as 10 degrees. (180 degree offset)

Example script that shows where the conversion issues are:

import phoebe
import numpy as np
import matplotlib.pyplot as plt

b = phoebe.default_binary()

### When flipping ecc and per0, per0 is returned in the wrong quadrent. 

b['ecc'] = 0.3
b['per0'] = 170

print(b['per0'],b['esinw'],b['ecosw'],b['ecc'])

b.flip_constraint('esinw', solve_for='per0')
b.flip_constraint('ecosw', solve_for='ecc')

print(b['per0'],b['esinw'],b['ecosw'],b['ecc'])

### Looking at the ecosw and esinw distributions

b.add_distribution({'esinw@binary': phoebe.uniform(-0.99,0.99), 'ecosw@binary': phoebe.uniform(-0.99,0.0)},overwrite_all=True, 
                   distribution='mydist')

### Both quadrents are covered here but I'm not convinced that the values after 0 are expected (90--180 and 270--360 deg)

dist = b.get_parameter('per0', component='binary', context='component').get_distribution('mydist')
plt.clf()
figure = plt.figure(figsize=(4,4))
_ = dist.plot(plot_uncertainties=False)
_ =plt.tight_layout()

dist = b.get_parameter('per0', component='binary', context='component').get_distribution('mydist')
plt.clf()
figure = plt.figure(figsize=(4,4))
plt.xlim(0,60)
_ = dist.plot(plot_uncertainties=False)
_ =plt.tight_layout()

### Now looking when ecosw is positive

b.add_distribution({'esinw@binary': phoebe.uniform(-0.99,0.99), 'ecosw@binary': phoebe.uniform(0.0,0.99)},overwrite_all=True, 
                   distribution='mydist2')

### Only one quadrent is visible:

dist = b.get_parameter('per0', component='binary', context='component').get_distribution('mydist2')
plt.clf()
figure = plt.figure(figsize=(4,4))
_ = dist.plot(plot_uncertainties=False)
_ =plt.tight_layout()

### For all ranges of e and omega:

b.add_distribution({'esinw@binary': phoebe.uniform(-0.99,0.99), 'ecosw@binary': phoebe.uniform(-0.99,0.99)},overwrite_all=True, 
                   distribution='mydist3')

### The third quadrent is still missing and the first quadrant has double the density.

dist = b.get_parameter('per0', component='binary', context='component').get_distribution('mydist3')
plt.clf()
figure = plt.figure(figsize=(4,4))
_ = dist.plot(plot_uncertainties=False)
_ =plt.tight_layout()

### There are two distinct ranges of ecosw and esinw that give very similar resulting distributions of omega, as shown below:

b.add_distribution({'esinw@binary': phoebe.uniform(-0.99,0.00), 'ecosw@binary': phoebe.uniform(-0.,0.99)},overwrite_all=True, 
                   distribution='mydist4')

dist = b.get_parameter('per0', component='binary', context='component').get_distribution('mydist4')
plt.clf()
figure = plt.figure(figsize=(4,4))
_ = dist.plot(plot_uncertainties=False)
_ =plt.tight_layout()

b.add_distribution({'esinw@binary': phoebe.uniform(-0.,0.99), 'ecosw@binary': phoebe.uniform(-0.,0.99)},overwrite_all=True, 
                   distribution='mydist5')

dist = b.get_parameter('per0', component='binary', context='component').get_distribution('mydist5')
plt.clf()
figure = plt.figure(figsize=(4,4))
_ = dist.plot(plot_uncertainties=False)
_ =plt.tight_layout()

### Eccentricity distribution 

### The eccentricity distribution also does not work if the negative values of ecosw and esinw are more negative than the positive values (or equal):

b.add_distribution({'esinw@binary': phoebe.uniform(-.99,0.99), 'ecosw@binary': phoebe.uniform(-.99,0.99)},overwrite_all=True, 
                   distribution='mydist6')

dist = b.get_parameter('ecc', component='binary', context='component').get_distribution('mydist6')
plt.clf()
figure = plt.figure(figsize=(4,4))
_ = dist.plot(plot_uncertainties=False)
_ =plt.tight_layout()

### It works for some ranges

b.add_distribution({'esinw@binary': phoebe.uniform(-0.6,0.99), 'ecosw@binary': phoebe.uniform(-0.6,0.99)},overwrite_all=True, 
                   distribution='mydist7')

### There are places where the eccentricity distribution does work. It seems to be related to when ecc_low that results from the conversion of esinw and ecosw to ecc is higher than the ecc_high. Due to the negative sign on ecosw and esinw being squared away in the conversion to ecc, sometimes the resulting ecc_low is higher than ecc_high. (I also tried adding a "wrap_at" function to ecc, but it is not considered (I didn't quite expect it to work, but I figured I'd try ;) ). 

### There are some occasions when I would expect the conversion to work but it doesn't, but I expect this has to do with the missing quadrent.

dist = b.get_parameter('ecc', component='binary', context='component').get_distribution('mydist7')
plt.clf()
figure = plt.figure(figsize=(4,4))
_ = dist.plot(plot_uncertainties=False)
_ =plt.tight_layout()