Eclipsing Binaries and Asteroseismology

Logo

Precise fundamental stellar parameters in the golden age of time-domain astronomy

  • Home
  • Registration
  • La Palma
  • Programme
  • Photos
  • Dynamic Asteroseismology

    Signal removal

    This notebook demonstrates the an example of quick and easy signal removal to search for periodicity underneath a stronger binary signal.

    For this notebook, you’ll need to download: tic049677785_lc.dat

    However, there are a number of other TESS lightcurves for you to try:
    tic0142874476_lc.dat
    tic0190693377_lc.dat
    tic0229585356_lc.dat
    tic0321947833_lc.dat
    tic0418292123_lc.dat
    tic048084398_lc.dat
    tic091369561_lc.dat
    tic097467902_lc.dat
    tic129764561_lc.dat
    tic152223725_lc.dat
    tic219238495_lc.dat
    tic231860752_lc.dat

    WARNING: Some lightcurves are in magnitude, others are in normalized flux – don’t blame me, blame past Cole

    %matplotlib notebook
    
    import numpy as np
    import matplotlib.pyplot as plt
    
    from pythia.utils.resampling import run_binning, run_mean_smooth
    from pythia.utils.conversions import time_to_ph
    from pythia.timeseries.periodograms import LS_periodogram
    from pythia.effective_models.binary import get_interp_model
    
    
    plt.rcParams.update({
        "text.usetex": True,
            "font.family": "sans-serif",
        "font.sans-serif": ["Helvetica"]})
    plt.rcParams.update({
        "pgf.rcfonts": False,
        "pgf.texsystem": "pdflatex",
        "pgf.preamble": "\n".join([
             r"\usepackage{amsmath}",
             r"\usepackage[utf8x]{inputenc}",
             r"\usepackage[T1]{fontenc}",
             r"\usepackage{cmbright}",
        ]),
    })
    
    plt.rcParams['xtick.labelsize']=18
    plt.rcParams['ytick.labelsize']=18
    
    times, signal = np.loadtxt('tic049677785_lc.dat').T
    
    fig, ax = plt.subplots(1,1,figsize=(6.6957, 6.6957), num=1)
    ax.plot(times, signal*1e3, 'kx')
    ax.set_xlabel(r'${\rm BJD-2457000}$',fontsize=18)
    ax.set_ylabel(r'${\rm Magnitude~[mmag]}$',fontsize=18)
    ax.set_ylim(ax.get_ylim()[::-1])
    fig.tight_layout()
    
    <IPython.core.display.Javascript object>
    

    nu, amp = LS_periodogram(times, signal-np.mean(signal))
    forb = nu[np.argmax(amp[nu<10.])]
    Aorb = amp[np.argmax(amp[nu<10.])]
    
    fig,ax = plt.subplots(1,1,figsize=(6.6957,6.6957),num=2)
    ax.plot(nu,amp*1e3,'k-')
    ax.plot(forb, Aorb*1e3, 'ro')
    ax.set_xlabel(r'$\nu\,\,{\rm [d^{-1}]}$',fontsize=14)
    ax.set_ylabel(r'${\rm Amplitude\,\,[mmag]}$',fontsize=14)
    ax.set_xlim(0.,45)
    fig.tight_layout()
    
    <IPython.core.display.Javascript object>
    

    %matplotlib inline
    import ipywidgets as wdg  
    from ipywidgets import interact
    
    def plot_folded_lc(period, t0, nbins):
        
    #     fig, ax = plt.subplots(1,1,figsize=(6.6957,6.6957), num=3)
        fig = plt.figure(constrained_layout=True, figsize=(6.6957*1.5,6.6957), num=4)
        axes = fig.subplot_mosaic(
        """
        AAA...
        AAACCC
        BBBCCC
        """)
        ph = time_to_ph(times, period, t0)
        bx, by, _ = run_binning(ph, signal, yerr=None, phStart=0., phStop=1., nbins=nbins)
        axes['A'].plot(ph, signal, 'k.')
        axes['A'].plot(ph+1, signal, 'k.')
        axes['A'].plot(bx, by, 'ro')
        axes['A'].plot(bx+1, by, 'ro')
        axes['A'].set_ylabel(r'${\rm Magnitude~[mmag]}$',fontsize=18)
        axes['A'].set_ylim(axes['A'].get_ylim()[::-1])
    
        
        ph_, x_, y_, ifunc = get_interp_model(times, signal, period, t0,
                                             phase_knots=[[0,1]], nbins=[nbins]
                                            )
        residuals = signal - ifunc(ph)
        axes['B'].plot(ph, residuals, 'k.')
        axes['B'].plot(ph+1, residuals, 'k.')
        axes['B'].set_ylabel(r'${\rm Residuals~[mmag]}$',fontsize=18)
        axes['B'].set_xlabel(r'${\rm Phase}$',fontsize=18)
        axes['B'].set_ylim(axes['B'].get_ylim()[::-1])
    
        idx = np.isfinite(residuals)
        nu_r, amp_r = LS_periodogram(times[idx], residuals[idx])
            
        axes['C'].plot(nu_r, amp_r*1e3)
        axes['C'].set_xlabel(r'${\rm Frequency~[d^{-1}]}$',fontsize=14)
        axes['C'].set_ylabel(r'${\rm Amplitude~[mmag]}$',fontsize=14)
        axes['C'].set_xlim(0.,45)
    
        fig.canvas.draw_idle()
        
    
    interact(plot_folded_lc, 
             period = wdg.FloatSlider(value=2./forb, description=r'$P_{\rm orb}$:',
                disabled=False, min=0, max=10., step=0.00001, readout_format='.6f'),
             t0 = wdg.FloatSlider( value=0., description='Reference Epoch:',
                disabled=False,min=-100, max=100., step=0.00001, readout_format='.6f') ,
             nbins = wdg.IntSlider(value=100, description='N bins: ', disabled=False, 
                                   min=50, max=1000, step=5, readout_format='d')
            )
    
    interactive(children=(FloatSlider(value=1.5932735895872152, description='$P_{\\rm orb}$:', max=10.0, readout_f…
    
    
    
    
    
    <function __main__.plot_folded_lc(period, t0, nbins)>
    

    Now let’s try it with a different rate of binning

    %matplotlib inline
    
    period_ = 1.594270
    t0_ = 0.7
    
    fig_ = plt.figure(constrained_layout=True, figsize=(6.6957*1.4,6.6957), num=5)
    axes_ = fig_.subplot_mosaic(
    """
    AAA...
    AAACCC
    BBBCCC
    """)
    phx = time_to_ph(times, period_, t0_)
    
    phase_knots = [ [0.,0.1], [0.11, 0.39], [0.4, 0.6], [0.61,0.89], [0.9,1.] ]
    nbins = [25,10,25,10,25]
    bx, by = [],[]
    
    for ii, knots in enumerate(phase_knots):
        bx_, by_, _ = run_binning(phx, signal, yerr=None, phStart=knots[0], phStop=knots[1], nbins=nbins[ii])
        bx.append(bx_)
        by.append(by_)
    
    bx,by = np.hstack(bx), np.hstack(by)
    
    
    axes_['A'].plot(phx, signal, 'k.')
    axes_['A'].plot(phx+1, signal, 'k.')
    axes_['A'].plot(bx, by, 'ro')
    axes_['A'].plot(bx+1, by, 'ro')
    axes_['A'].set_ylabel(r'${\rm Magnitude~[mmag]}$',fontsize=18)
    axes_['A'].set_ylim(axes_['A'].get_ylim()[::-1])
    
    
    ph_, x_, y_, ifunc = get_interp_model(times, signal, period_, t0_,
                                         phase_knots=phase_knots, nbins=nbins
                                        )
    residuals = signal - ifunc(phx)
    axes_['B'].plot(phx, residuals, 'k.')
    axes_['B'].plot(phx+1, residuals, 'k.')
    axes_['B'].set_ylabel(r'${\rm Residuals~[mmag]}$',fontsize=18)
    axes_['B'].set_xlabel(r'${\rm Phase}$',fontsize=18)
    axes_['B'].set_ylim(axes_['B'].get_ylim()[::-1])
    
    idx = np.isfinite(residuals)
    nu_r, amp_r = LS_periodogram(times[idx], residuals[idx])
    
    axes_['C'].plot(nu_r, amp_r*1e3)
    axes_['C'].set_xlabel(r'${\rm Frequency~[d^{-1}]}$',fontsize=14)
    axes_['C'].set_ylabel(r'${\rm Amplitude~[mmag]}$',fontsize=14)
    axes_['C'].set_xlim(0.,45)
    axes_['C'].tick_params(axis='y', which='both', left=True, right=True, labelleft=False, labelright=True)
    
    
    
    
    
    Adopted bin-width:  0.004
    Adopted bin-width:  0.028000000000000004
    Adopted bin-width:  0.007999999999999998
    Adopted bin-width:  0.028000000000000004
    Adopted bin-width:  0.003999999999999999
    Adopted bin-width:  0.004
    Adopted bin-width:  0.028000000000000004
    Adopted bin-width:  0.007999999999999998
    Adopted bin-width:  0.028000000000000004
    Adopted bin-width:  0.003999999999999999
    

    png

    The residuals

    Now that we’re confident that there’s some signal left over, lets see what we can do with it!

    # Thanks, Dan Hey!
    from echelle import interact_echelle
    
    interact_echelle(nu_r[nu<=45], amp_r[nu<=45], dnu_min=0.1, dnu_max=10.)
    
    /lhome/colej/software/miniconda3/envs/nominal3p9/lib/python3.9/site-packages/echelle-1.5.1-py3.9.egg/echelle/interact.py:201: UserWarning: You have multiple Jupyter servers open.             You will need to pass the current notebook to `notebook_url`.             i.e. interact_echelle(x,y,notebook_url='http://localhost:8888')