phoebe-project / phoebe2

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

Spikes in the light curve #579

Closed mvelavok closed 1 year ago

mvelavok commented 2 years ago

Hello PHOEBE team!

I found strange behavior of PHOEBE 2.3 during modelling of the TESS light curve. LC shows 0.15% spikes, which originates from changes in spot shape. I attach script that will generate plots, where such spikes are clearly visible at times 0.147,0.176,0,189 d,

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

b = phoebe.default_binary()
qq=0.2106537891#1./4.75
b.set_value('q', qq)

b.add_constraint('semidetached', 'secondary')

period=7.054905
b.set_value('period@binary', period)
b.set_value('teff', component='primary', value=5592.)
b.set_value('teff', component='secondary', value=4769.)
b.flip_constraint('asini@binary', solve_for='sma@binary')

b.add_spot(component='primary', feature='spot01')
b.add_constraint(b.get_parameter('requiv@primary'),24./25.*b.get_parameter("requiv_max@secondary@component"))

#solution with spot on primary and default gde and albedo
b["t0_supconj@binary@orbit@component"]=0.051233149902048816#,0.022060250657176605
b["gravb_bol@primary@star@component"]=0.32
b["gravb_bol@secondary@star@component"]=0.32
b["incl@binary@orbit@component"]=39.80638565338866#49.79894994015926
b["irrad_frac_refl_bol@primary@star@component"]=0.5#0.09293901423605994
b["irrad_frac_refl_bol@secondary@star@component"]=0.5#0.9262508771090747
b["asini@binary@orbit@component"]=12.201941180329381
#spot
b['colat@spot01']=93.5191972757996
b['long@spot01']=7.53262022863316
b['radius@spot01']=25.759455456205373
b['relteff@spot01']=0.9209290488374553

data=np.array([
[0.06408226820591345, 0.9602671797842749, 0.9594889996490671],
[0.07217726820055859, 0.9602935274506741, 0.961169541515931],
[0.08508226820655374, 0.9603522735186585, 0.9601078045256515],
[0.09217726820463312, 0.9604130847103299, 0.9612580725791999],
[0.10608226819991806, 0.9605309462466235, 0.9593122715214302],
[0.1131772682052734, 0.9605714350758937, 0.9616122783836647],
[0.12708226820055835, 0.9607094935255753, 0.9593122715214302],
[0.1341772682059137, 0.9607607863257719, 0.9613466117968567],
[0.14708226820463288, 0.9624753486875852, 0.9604615864776328],# x
[0.15517726820655398, 0.9610153621768093, 0.9627643488891963],
[0.16808226820527317, 0.9611696674479447, 0.9591355759453193],
[0.1761772681999183, 0.9627179390744364, 0.9621438318271684],# x
[0.18908226820591345, 0.9629045969450848, 0.9607270084891201],# x
[0.1971772682005586, 0.9614860681118392, 0.9632078203620027],
[0.21008226820655374, 0.9616477664000752, 0.9610810186062989],
[0.21717726820463312, 0.9617170359138295, 0.9640065838447285],
[0.23108226819991806, 0.9619185236930053, 0.9608154987917139],
[0.2381772682052734, 0.9620368162918097, 0.963474001290798],
[0.25208226820055835, 0.9622600967878215, 0.9616122783836647],
[0.2591772682059137, 0.962434085232735, 0.9640065838447285]
])

timesli=data[:,0]
models=data[:,1]
sapi=data[:,2]
esap=sapi*0.0+1e-3

phaseli=((timesli-b.get_value('t0_supconj@binary@orbit@component'))%period)/period

b.add_dataset('lc', overwrite=True,passband= 'TESS:T',
              compute_phases=phaseli,
              times=timesli, 
              fluxes=sapi, 
              sigmas=esap, 
              dataset='lc01')
b.add_dataset('mesh', times=timesli, columns=['teffs', 'intensities@lc01'])
b.set_value_all('pblum_mode', 'dataset-scaled')
b.run_compute(irrad_method='horvat')
fluxes_horvat = b.get_value('fluxes', context='model')
fig1=plt.figure()#figsize=[8,12])

b['lc01@model'].plot( fig=fig1,save="time_lc.png")
plt.clf()
for i in range(20):
  print(timesli[i],fluxes_horvat[i],sapi[i])
  b["mesh@model"].plot(time=timesli[i],fc='intensities@lc01',ec="None", fig=fig1,save="time_%.3f.png"%(timesli[i]))
  plt.clf()

I think that it is due to the finite mesh resolution and very dense phase coverage. Is it possible to make spot to always use same set of triangles in the mesh?

regards,
Mikhail time_lc times_li0 155 times_li0 147

amiszuda commented 2 years ago

The reason behind sudden spikes in the flux is, as you mentioned, the finite mesh; these are simply the numerical effects you see. Try to increase the number of triangles, it will smooth your flux.

As for your second question, no I don't think you can always use the same set of triangles in the mesh. As the binary orbits around the barycentre, the visibility not only of the spot but also of the Roche distortion changes in time. Therefore, PHOEBE will need to recalculate the mesh for each time point you provide. To visualise that, I changed your script a bit. I hope it helps.

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

b = phoebe.default_binary()
qq=0.2106537891#1./4.75
b.set_value('q', qq)

b.add_constraint('semidetached', 'secondary')

period=1#7.054905
b.set_value('period@binary', period)
b.set_value('teff', component='primary', value=5592.)
b.set_value('teff', component='secondary', value=4769.)
b.flip_constraint('asini@binary', solve_for='sma@binary')

b.add_spot(component='primary', feature='spot01')
b.add_constraint(b.get_parameter('requiv@primary'),24./25.*b.get_parameter("requiv_max@secondary@component"))

#solution with spot on primary and default gde and albedo
b["t0_supconj@binary@orbit@component"]=0.051233149902048816#,0.022060250657176605
b["gravb_bol@primary@star@component"]=0.32
b["gravb_bol@secondary@star@component"]=0.32
b["incl@binary@orbit@component"]=39.80638565338866#49.79894994015926
b["irrad_frac_refl_bol@primary@star@component"]=0.5#0.09293901423605994
b["irrad_frac_refl_bol@secondary@star@component"]=0.5#0.9262508771090747
b["asini@binary@orbit@component"]=12.201941180329381
#spot
b['colat@spot01']=93.5191972757996
b['long@spot01']=7.53262022863316
b['radius@spot01']=25.759455456205373
b['relteff@spot01']=0.9209290488374553

b['ntriangles@primary']   = 1000 
b['ntriangles@secondary'] = 1000 

# data=np.array([
# [0.06408226820591345, 0.9602671797842749, 0.9594889996490671],
# [0.07217726820055859, 0.9602935274506741, 0.961169541515931],
# [0.08508226820655374, 0.9603522735186585, 0.9601078045256515],
# [0.09217726820463312, 0.9604130847103299, 0.9612580725791999],
# [0.10608226819991806, 0.9605309462466235, 0.9593122715214302],
# [0.1131772682052734, 0.9605714350758937, 0.9616122783836647],
# [0.12708226820055835, 0.9607094935255753, 0.9593122715214302],
# [0.1341772682059137, 0.9607607863257719, 0.9613466117968567],
# [0.14708226820463288, 0.9624753486875852, 0.9604615864776328],# x
# [0.15517726820655398, 0.9610153621768093, 0.9627643488891963],
# [0.16808226820527317, 0.9611696674479447, 0.9591355759453193],
# [0.1761772681999183, 0.9627179390744364, 0.9621438318271684],# x
# [0.18908226820591345, 0.9629045969450848, 0.9607270084891201],# x
# [0.1971772682005586, 0.9614860681118392, 0.9632078203620027],
# [0.21008226820655374, 0.9616477664000752, 0.9610810186062989],
# [0.21717726820463312, 0.9617170359138295, 0.9640065838447285],
# [0.23108226819991806, 0.9619185236930053, 0.9608154987917139],
# [0.2381772682052734, 0.9620368162918097, 0.963474001290798],
# [0.25208226820055835, 0.9622600967878215, 0.9616122783836647],
# [0.2591772682059137, 0.962434085232735, 0.9640065838447285]
# ])
# timesli=data[:,0]
# models=data[:,1]
# sapi=data[:,2]
# esap=sapi*0.0+1e-3

timesli = np.linspace(0,1,100)
sapi = np.linspace(0.96,0.97,100)
esap=sapi*0.0+1e-3

phaseli=((timesli-b.get_value('t0_supconj@binary@orbit@component'))%period)/period

b.add_dataset('lc', overwrite=True,passband= 'TESS:T',
              compute_phases=phaseli,
              times=timesli, 
              fluxes=sapi, 
              sigmas=esap, 
              dataset='lc01')

b.set_value('atm', component='primary', value='blackbody')
b.set_value('atm', component='secondary', value='blackbody')
b.set_value('ld_mode', dataset='lc01', component='primary', value='manual')
b.set_value('ld_mode', dataset='lc01', component='secondary', value='manual')

b.add_dataset('mesh', times=timesli, columns=['teffs', 'intensities@lc01'])
b.set_value_all('pblum_mode', 'dataset-scaled')
b.run_compute()
fluxes_horvat = b.get_value('fluxes', context='model')
fig1=plt.figure()#figsize=[8,12])

b['lc01@model'].plot( fig=fig1,save="time_lc.png")
plt.clf()
b["mesh@model"].plot(x='us', y='vs', fc='intensities', 
                      animate=True, save='single_spots_1.gif', save_kwargs={'writer': 'imagemagick'})

single_spots_1

mvelavok commented 2 years ago

thank you! Can you point me how to change the resolution of the mesh?

kecnry commented 2 years ago

@amiszuda's suggestion of changing the mesh resolution (through the ntriangles parameter on the star with the spot) should help... but I also don't think the mesh needs to be regenerated or the location of the spot recomputed in this case since you have a circular orbit and synchronous rotation, so there's a chance that logic could be optimized and remove the numerical noise without having to increase the resolution (which will also increase the time cost).

UPDATE: the mesh is currently recomputed if there are any spots, but that logic can probably be extended to avoid remeshing in the synchronous case.