Create a B1259-63 (aka PSR J1302-6350) light curve using Bayesian Blocks

Special data check

clear = True
#update_cache(); 

class B1259Periastron[source]

B1259Periastron(**kwargs) :: WtLike

weekly = self=B1259Periastron(clear= clear)
SourceData: photons and exposure for PSR B1259-63: Saving to cache with key "PSR_B1259-63_data"
loading all filesSourceData: Assembling photon data and exposure for source PSR B1259-63
	from folder "/home/burnett/wtlike_data/data_files", with 670 files,
	Weights from file PSR_B1259-63_weights.pkl
..............................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................
SourceData: Source PSR B1259-63 with:
	 data:       881,665 photons from 2008-08-04 to 2021-06-07
	 exposure: 3,370,383 intervals,  average flux 3477 cm^2 for 100.6 Ms
	 rates:  source 8.18e-09/s, background 2.51e-06/s, S/N ratio 3.26e-03
CellData: Bin photon data into 669 1-week bins from 54683.0 to 59366.0
LightCurve: select 662 cells for fitting with e>1 & n>2
def B1259( clear=False):
    r"""
    ## Fit to all data

    {date}
    
    Create a `WtLike` object with all the data
    
    {outp}
 
    ## The full weekly-interval light curve, showing the BB partitions
    {out2}
    {fig2}
    Table of fits (note that the "flux" column is relative to the 12-year count flux measurement
    6.7e-9 cm**-2 s**-1.
    {bbf}
    
    ## Expand about each periastron
 
    #### Periastron dates

    Assuming {period}-day orbital period, the MJD and UTC values are:
    
    {utc}
 
    Expand the above, with 1-day bins, following those dates:

    {fig3}
    
    ## Recent detail 
    {out4}{fig4}
    
    """
    global weekly #  make availlable for follow-up cells

 
    pd.set_option('display.precision', 4)#, 'display.colheader_justify','left')

    
    with capture_print('Output from analysis: create cells, fit each, run BB partition, fit partitions.') as outp:
        if weekly is None or clear:
            weekly = B1259Periastron(time_bins=(0,0,7), clear=clear)#.bb_view()
        else:
            print('(already done)')
        bb = weekly.bb_view(0.5)
    
    bbf = monospace(str(bb.fluxes), 'BB fit table', open=False)

    # fig 2: full light curbe
    plt.rc('font', size=16)
    
    fig2 = figure( bb.plot(fignum=2, figsize=(15,4), title='Full weekly light curve'),
                  width=600)
    
    # fig 3 -- stack light curves
    utc = monospace(str(weekly.date_info())) 
    with capture_print('Output from periastron analyses') as out2:
        fig3 = figure( weekly.stacked_plots( ylim=(4,300), fignum=3), width=600)
    
    # fit 4 -- last two weeks
    with capture_print(f'Output, with light curve table, from selecting the last two weeks') as out4:
        recent_wk = weekly.view(-14,0,1)
        print(recent_wk.fluxes)
        fig4=figure(
            recent_wk.plot(show_flux=True, ylim=(-.1,2), title='Last 2 weeks',
                       xlabel=f'MJD, until {UTC(weekly.stop)}',fignum=4),
            caption=None, width=600)
  
    return locals()

nbdoc(B1259, clear=False)

Fit to all data

2021-06-07 07:08

Create a WtLike object with all the data

Output from analysis: create cells, fit each, run BB partition, fit partitions.
(already done)
LightCurve: select 662 cells for fitting with e>1 & n>2
Bayesian Blocks: using penalty 0.5
Partitioned 662 cells into 31 blocks, using LikelihoodFitness
LightCurve: Loaded 31 / 31 cells for fitting

The full weekly-interval light curve, showing the BB partitions

Output from periastron analyses
CellData: Bin photon data into 125 1-day bins from 55544.7 to 55669.7
LightCurve: select 119 cells for fitting with e>1 & n>2
CellData: Bin photon data into 125 1-day bins from 56781.4 to 56906.4
LightCurve: select 125 cells for fitting with e>1 & n>2
CellData: Bin photon data into 125 1-day bins from 58018.1 to 58143.1
LightCurve: select 125 cells for fitting with e>1 & n>2
CellData: Bin photon data into 125 1-day bins from 59254.9 to 59379.9
LightCurve: select 118 cells for fitting with e>1 & n>2
LightCurve: select 119 cells for fitting with e>1 & n>2
Bayesian Blocks: using penalty 0.5
Partitioned 119 cells into 12 blocks, using LikelihoodFitness
LightCurve: Loaded 12 / 12 cells for fitting
LightCurve: select 125 cells for fitting with e>1 & n>2
Bayesian Blocks: using penalty 0.5
Partitioned 125 cells into 20 blocks, using LikelihoodFitness
LightCurve: Loaded 20 / 20 cells for fitting
LightCurve: select 125 cells for fitting with e>1 & n>2
Bayesian Blocks: using penalty 0.5
Partitioned 125 cells into 16 blocks, using LikelihoodFitness
LightCurve: Loaded 16 / 16 cells for fitting
LightCurve: select 118 cells for fitting with e>1 & n>2
Bayesian Blocks: using penalty 0.5
Partitioned 118 cells into 15 blocks, using LikelihoodFitness
LightCurve: Loaded 15 / 15 cells for fitting

Figure 1 at images/B1259_fig_01.png
Table of fits (note that the "flux" column is relative to the 12-year count flux measurement 6.7e-9 cm**-2 s**-1.
BB fit table
          t      tw       n     ts      flux           errors     limit
0 55103.0 840.0 145438 0.0 0.0000 (0, 0.039) 0.2348
1 55544.0 42.0 8515 37.3 12.1020 (-2.155, 2.17) 15.6960
2 55568.5 7.0 2302 1.6 1.9816 (-1.715, 2.362) 7.2135
3 55575.5 7.0 1923 64.3 40.5665 (-5.752, 5.841) 50.3100
4 55582.5 7.0 1746 253.9 97.9705 (-7.73, 7.855) 111.0829
5 55593.0 14.0 2654 228.5 70.7162 (-5.669, 5.743) 80.2762
6 55607.0 14.0 2879 89.3 39.0787 (-4.729, 4.789) 47.0466
7 55621.0 14.0 3895 11.8 8.9265 (-2.849, 3.011) 14.1380
8 56083.0 910.0 168666 0.0 0.0000 (0, 0.113) 0.4192
9 56660.5 245.0 40193 0.0 0.0000 (0, 0.367) 1.0271
10 56797.0 28.0 7546 18.4 9.1142 (-2.223, 2.277) 12.9447
11 56814.5 7.0 2803 95.7 39.4052 (-4.654, 4.713) 47.2486
12 56828.5 21.0 7207 623.0 71.0512 (-3.433, 3.46) 76.7843
13 56842.5 7.0 1899 2.6 7.2258 (-4.534, 4.712) 15.3769
14 56849.5 7.0 1622 121.6 64.9796 (-7.021, 7.138) 76.9030
15 56867.0 28.0 3914 10.3 9.4505 (-3.08, 3.193) 14.8815
16 57420.0 1078.0 201548 0.0 0.0000 (0, 0.055) 0.2680
17 57987.0 56.0 10849 6.1 4.0526 (-1.674, 1.725) 6.9756
18 58029.0 28.0 5319 63.1 20.6815 (-2.901, 2.927) 25.5374
19 58050.0 14.0 2516 0.1 0.8438 (-0.844, 3.572) 7.8938
20 58060.5 7.0 2238 164.4 62.2969 (-5.884, 5.967) 72.2418
21 58067.5 7.0 2181 53.8 34.2157 (-5.223, 5.298) 43.0462
22 58074.5 7.0 1844 373.8 117.5356 (-7.934, 8.059) 130.9832
23 58085.0 14.0 2509 169.6 57.2343 (-5.268, 5.339) 66.1248
24 58116.5 49.0 11403 11.0 5.1192 (-1.627, 1.689) 7.9951
25 58680.0 1078.0 206040 0.0 0.0000 (0, 0.045) 0.2675
26 59264.5 91.0 21018 31.3 6.4596 (-1.25, 1.255) 8.5331
27 59313.5 7.0 1159 33.1 33.4177 (-6.656, 6.788) 44.7878
28 59331.0 28.0 4825 324.2 60.6652 (-3.991, 4.03) 67.3529
29 59352.0 14.0 3080 406.4 93.8585 (-5.829, 5.899) 103.6710
30 59362.5 7.0 1028 34.9 37.3580 (-7.322, 7.476) 49.8937
## Expand about each periastron #### Periastron dates Assuming {period}-day orbital period, the MJD and UTC values are:
         MJD               UTC
0 55544.694 2010-12-14 16:39
1 56781.418 2014-05-04 10:02
2 58018.143 2017-09-22 03:25
3 59254.867 2021-02-09 20:48
Expand the above, with 1-day bins, following those dates:
Figure 2 at images/B1259_fig_02.png

Recent detail

Output, with light curve table, from selecting the last two weeks
CellData: Bin photon data into 14 1-day bins from 59358.0 to 59372.0
LightCurve: select 14 cells for fitting with e>1 & n>2
t tw n ts flux errors limit
0 59358.5 1.0 168 3.5 35.4451 (-19.778, 21.782) 74.8803
1 59359.5 1.0 179 9.9 45.9444 (-17.017, 19.28) 81.4914
2 59360.5 1.0 46 7.2 111.0989 (-47.705, 55.802) 216.9795
3 59361.5 1.0 138 5.4 37.4053 (-17.361, 19.39) 72.7737
4 59362.5 1.0 160 9.6 50.5254 (-18.247, 20.17) 86.9032
5 59363.5 1.0 174 1.4 19.0946 (-16.269, 18.613) 54.7079
6 59364.5 1.0 165 2.3 21.0946 (-14.533, 16.71) 52.7782
7 59365.5 1.0 166 4.2 35.4112 (-18.263, 20.138) 71.8151
8 59366.5 1.0 198 13.3 62.4958 (-19.596, 21.385) 100.6115
9 59367.5 1.0 119 0.0 0.0000 (0, 16.027) 40.9647
10 59368.5 1.0 106 0.0 0.0000 (0, 14.509) 42.0175
11 59369.5 1.0 106 1.9 25.2428 (-19.161, 22.18) 67.7728
12 59370.5 1.0 115 0.0 0.0000 (0, 8.212) 32.6119
13 59371.5 1.0 149 0.4 12.4996 (-12.5, 21.473) 54.6739
Figure 3 at images/B1259_fig_03.png

mdf = d.query('59319>t>59093').copy()
mdf.loc[:,'ts']= mdf.fit.apply(lambda x: x.ts)
mdf.loc[:,'flux']=mdf.fit.apply(lambda x: round(x.flux,1) )
mdf.loc[:,'errors']= mdf.fit.apply(lambda x: (np.array(x.errors)-x.flux).round(1))
mdf.loc[:,'limit']= mdf.fit.apply(lambda x: x.limit.round(1))
#mdf['t ts flux errors limit'.split()].describe(percentiles=[])

### Use the last one, at 59318.5, +63, to make tentative normalization


f = mdf.iloc[-1].fit
print(f'Compare point at {mdf.iloc[-1].t}: TS: me {f.ts:.1f}, Tyrel {tdf.iloc[-1].TS:.1f}')
a,b = f.flux, tdf.iloc[-1].flux; a,b, 
print(f'flux ratio check: {a/b:.3e}')

fig, (ax1,ax2,ax3) = plt.subplots(1,3,figsize=(15,5))
ax1.plot(mdf.ts, tdf.TS, 'o');
ax1.set(xlabel='wtlike TS', ylabel='gtlike TS')
ax1.grid(alpha=0.5)
mtscut= (mdf.ts>4).values
ax2.plot(mdf.flux, tdf.flux, 'o');
ax2.plot(mdf.flux[mtscut], tdf[mtscut].flux, 'or', label='wtlike TS>4');
ax2.plot([0,70], [0,70], '--', color='grey');
ax2.set(xlabel='wtlike flux', ylabel=f'gtlike flux x {flux_factor:.2e}', xlim=(0,70), ylim=(0,70))
ax2.legend();ax2.grid(alpha=0.5);

ttscut=(tdf.TS>4).values
ax3.plot(mdf.flux, tdf.flux, 'o');
ax3.plot(mdf.flux[ttscut], tdf[ttscut].flux, 'or', label='gtlike TS>4');
ax3.plot([0,70], [0,70], '--', color='grey');
ax3.legend()
ax3.set(xlabel='wtlike flux', ylabel=f'gtlike flux x {flux_factor:.2e}', xlim=(0,70), ylim=(0,70))
ax3.legend();ax3.grid(alpha=0.5);
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-4-e4555ecb1020> in <module>
     24 
     25 
---> 26 lc=daily.lc;
     27 mdf = lc.dataframe.query('59319>t>59093').copy()
     28 mdf.loc[:,'ts']= mdf.fit.apply(lambda x: x.ts)

NameError: name 'daily' is not defined