Spikes in the light curve

[mvelavok]

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

[amiszuda]

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'})

[mvelavok]

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

[kecnry]

@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.