<jupyter_start><jupyter_text>Fitting a sine curve to a pceb light curve with astropy fitting. https://docs.astropy.org/en/stable/modeling/index.html<jupyter_code>#usual bunch of imports

import numpy as np
import matplotlib.pyplot as plt
import astropy.io.fits as fits
import os
import glob
from astropy.table import Table
from astropy.io import ascii
import astropy.units as u
import lightkurve as lk
from matplotlib.gridspec import GridSpec



#matplotlib set up
%matplotlib inline
from matplotlib import rcParams
rcParams["figure.figsize"] = (14, 5)
rcParams["font.size"] = 20

from astropy.modeling import models, fitting

#grab a light curve and do all the normalisation etc. Using the PCEB CC Cet that I'm quite familar with.

search_result = lk.search_lightcurve('CC Cet')
search_result
lc = search_result[0].download(download_dir='data/') #there is a 20s light curve, but I'll use the 2min one as it's the data you're familar with
# lcn = lc.flatten(10001)
lcn = lc.normalize()
lcn.flux = lcn.flux.value -1
lcn.flux = lcn.flux*100
lcn.flux_err = lcn.flux_err.value*100
lcn=lcn.remove_nans()
lcn.scatter()
pg = lcn.to_periodogram(oversample_factor=10, minimum_period=4*u.min, maximum_period = 10*u.d)
plt.plot(pg.period, pg.power)
period = pg.period_at_max_power
print(pg.max_power)
print(pg.period_at_max_power)
plt.xscale('log')<jupyter_output>3.4612823715569814
0.2868948066810429 d<jupyter_text>Very easy period with no complicated extra harmonics. First we fit a sine curve to the whole light curve to get a precise peirod. Then we fold the light curve onto the period and fit for the amplitude and phase.<jupyter_code># -3.695494014784872
-2.6953779471917727  % 1

fitter = fitting.LevMarLSQFitter()

sin_mod =  models.Sine1D(amplitude=pg.max_power.value, frequency= pg.frequency_at_max_power.value, phase=0.5)#, bounds=dict(phase=[0,5]))
sin_fit = fitter(sin_mod, lcn.time.value, lcn.flux.value, weights = 1/lcn.flux_err.value, maxiter=100000)
print(sin_fit)
sin_fit_e = np.sqrt(np.diag(fitter.fit_info['param_cov']))
plt.scatter(lcn.time.value, lcn.flux,s=5,alpha=0.5) 
plt.plot(lcn.time.value, sin_fit(lcn.time.value), c='C1', lw=2)
print('Amplitude =', sin_fit.amplitude.value, sin_fit_e[0], '%')
print('Phase =', sin_fit.phase.value, sin_fit_e[2])
#now fold and refit


lcf = lcn.fold(1/sin_fit.frequency.value, normalize_phase=True, epoch_time=eph1-2457000)



sin_mod2 =  models.Sine1D(amplitude=sin_fit.amplitude.value, frequency= 1,phase=-0.5)# sin_fit.phase.value % 1), fixed=dict(frequency=True)) #note frequency for the fit is now fixed
sin_fit2 = fitter(sin_mod2, lcf.time.value, lcf.flux.value, weights = 1/lcf.flux_err.value, maxiter=1000)
print(sin_fit2)
sin_fit_e2 = np.sqrt(np.diag(fitter.fit_info['param_cov']))
plt.scatter(lcf.time.value, lcf.flux,s=5,alpha=0.5) 
plt.plot(lcf.time.value, sin_fit2(lcf.time.value), c='C1', lw=2)
print(pg.period_at_max_power.to(u.h))
print('Period =', (1/sin_fit.frequency.value)*24, (sin_fit_e[1]/(sin_fit.frequency.value**2))*24, 'h')
print('Amplitude =', sin_fit2.amplitude.value, sin_fit_e2[0], '%')
print('Phase =', sin_fit2.phase.value, sin_fit_e2[1])
period = 1/sin_fit.frequency<jupyter_output><empty_output><jupyter_text>Referee wants a zoomed in version<jupyter_code>plt.scatter(lcn.time.value, lcn.flux,s=5,alpha=0.5) 
plt.plot(lcn.time.value, sin_fit(lcn.time.value), c='C1', lw=2)
plt.xlim(2145, 2146)
t0 = lcn.time[0].value + (1/sin_fit.frequency)*((sin_fit2.phase)+0.75) #0.5 = inferior conjunction
print(t0)
tmid = np.mean(lcn.time.value)
print(tmid)
tmid2 = lcn.time[0].value + (0.5 * (lcn.time[-1].value - lcn.time[0].value))
print(tmid2)
# eph = t0 + tmid2- lcn.time[0].value
# eph = tmid2 + (1/sin_fit.frequency)*(sin_fit2.phase)
nhalf = (0.5 * (lcn.time[-1].value - lcn.time[0].value))*sin_fit.frequency.value #how many periods
# print(int(nhalf)/sin_fit.frequency.value)

eph = t0 + int(nhalf)/sin_fit.frequency.value
print(eph)
# eph = 2157.16489

plt.xlim(2156, 2159)
plt.axvline(eph)
#plt.axvline(tmid)
# print(sin_fit.phase, sin_fit_e[2])
plt.scatter(lcn.time.value, lcn.flux,s=5,alpha=0.5) 
timefit = np.linspace(lcn.time[0].value, lcn.time[-1].value, 1000)

sin_fit3 =  models.Sine1D(amplitude=sin_fit2.amplitude.value, frequency= sin_fit.frequency.value, phase= sin_fit2.phase.value)
plt.plot(timefit, sin_fit(timefit), c='C1', lw=2)

print(eph+2457000)
print(1/sin_fit.frequency.value*sin_fit_e2[1])

print(lcn.time[0], t0)
#old fashoned way
timefit = np.linspace(tmid2 - 0.5*period, tmid2 + 0.5*period, 1000000)
plt.plot(timefit,  sin_fit(timefit))
print(timefit[np.argmin(sin_fit(timefit))])
print(timefit[np.argmin(sin_fit(timefit))] + 2457000)
eph1 = timefit[np.argmin(sin_fit(timefit))] + 2457000
lcf = lcn.fold(period, epoch_time= 2459157.146- 2457000)
lcf.bin(10*u.min).scatter()
print(period*u.d.to(u.hr))
print(pg.period_at_max_power)
plt.scatter(lcn.time.value, lcn.flux,s=5,alpha=0.5) 
plt.plot(lcn.time.value, sin_fit(lcn.time.value), c='C1', lw=2)
plt.xlim(2144, 2146)
plt.axvline(t0)
fake_peaks = np.load('fake_peaks5.npy')
fa99 = np.percentile(fake_peaks, 99.9)
print(fa99)
lcf = lcn.fold(1/sin_fit.frequency.value, epoch_time=t0).bin(0.1*u.hr)
lcf.errorbar()
plt.figure(figsize=(8, 14))
gs = GridSpec(3, 1)
plt.subplot(gs[1, :])
pg20h = pg.period.to(u.hour)
plt.plot(pg20h, pg.power)
plt.xscale('log')
plt.axhline(fa99, c='C1', ls='--', lw=2)

plt.xlim(0.05, 14*24)#pg20h[0].value)
plt.ylim(0.0001)
plt.ylabel('Amplitude (\%)', size=20)
plt.xlabel('Period (h)', size=20)
#axins = inset_axes(ax, width='40%', height='30%', loc=2)
#axins.plot(pg1m[mask], pg1.power[mask])
plt.annotate('P\,=\,{0}\,$\pm$\,{1}\,h \n A\,=\,{2}\,$\pm$\,{3}\,\%'.format(6.88233,0.00045, 3.461,  0.039), (0.1, 3.), xycoords='data')
#6.882326693607637 0.0004561525821496103
#-3.4613593519684556 0.039138214912684814

# plt.annotate('0.1\% False Alarm Probability = 0.243\%', (0.95, fa99+0.1), xycoords=('axes fraction', 'data'), ha='right', va='bottom')

plt.subplot(gs[0, :])


plt.scatter(lcn.time.value, lcn.flux,s=5,alpha=0.5) 
plt.plot(lcn.time.value, sin_fit(lcn.time.value), c='C1', lw=2)




plt.ylabel('Normalised Flux (\%)', size=20)
plt.xlabel('Time (BJD - 2457000)', size=20)
plt.ylim(-39.9, 39.9)

gs1 = gs[2].subgridspec(3, 1, hspace=0)
plt.subplot(gs1[:2, :])
lcf = lcn.fold(1/sin_fit.frequency.value, epoch_time=eph).bin(0.1*u.hr)


# plt.scatter(lcf.time.value*sin_fit.frequency.value, lcf.flux,s=5,alpha=0.5) 
# plt.plot(lcf.time.value*sin_fit.frequency.value, sin_fit2(lcf.time.value), c='C1', lw=2)

# lcf_mod =  models.Sine1D(amplitude=fold_fit.amplitude, frequency= 1, phase=0, fixed={'amplitude': True, 'frequency':True})
# lcf_fit = fitter(lcf_mod, lcf.time, lcf.flux, weights = 1/lcf.flux_err, maxiter=100000)
#two cycles
plt.errorbar(lcf.time.value*sin_fit.frequency.value, lcf.flux, yerr=lcf.flux_err, ls='none', marker='o', c='C0')
plt.plot(lcf.time.value*sin_fit.frequency.value, sin_fit2(lcf.time.value), c='C1', lw=2)
# plt.errorbar(lcf.time+1, lcf.flux, yerr=lcf.flux_err, ls='none', marker='o', c='C0')
# plt.plot(fold_time+1, lcf_fit(fold_time), c='C1', lw=2, zorder=10)

plt.ylabel('Normalised Flux (\%)', size=20)
plt.xticks(visible=False)

xlims = [-0.49, 1.49] #storing here to update both plots at once
plt.xlim(xlims[0], xlims[1])

plt.subplot(gs1[2:, :])
mod = sin_fit2(lcf.time.value)
plt.errorbar(lcf.time.value*sin_fit.frequency.value, lcf.flux-mod, yerr=lcf.flux_err, ls='none', marker='o', c='C0')
plt.errorbar(lcf.time.value*sin_fit.frequency.value+1, lcf.flux-mod, yerr=lcf.flux_err, ls='none', marker='o', c='C0')
plt.xlim(xlims[0], xlims[1])
plt.ylabel('O-C (\%)', size=20)

plt.xlabel('Phase', size=20)
plt.tight_layout()
#plt.subplots_adjust(hspace=0.18)
# ax.subplots_adjust(hspace=0.7)
#ax1.set_top=0.99
#ax.set_bottom=0.01
# plt.savefig('plots/cc_cet_tess_20s.pdf', dpi=150)
time = np.linspace(0, 1, 100)
sine = models.Sine1D(1, 1, 0)
plt.plot(time, sine(time))
time = np.linspace(0, 1, 100)
sine = models.Sine1D(1, 1, 0.75)
plt.plot(time, sine(time))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.18362004487658937, 'fraction_numerical': 0.0731737721266517, 'mean_word_length': 5.964409722222222, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 12, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '229198', 'license': 'MIT', 'n_tokens_mistral': 3918, 'n_tokens_neox': 3540, 'n_words': 573, 'path': 'cc_cet/.ipynb_checkpoints/cc_cet_20s_tess-checkpoint.ipynb', 'repo_name': 'davidjwilson/pceb', 'size': '8022'}