<jupyter_start><jupyter_text>C1 = BJD_UTC-2450000C2 = fluxC3 = errorC4 = XC5 = YC6 = fwhmC7 = fwhm-xC8 = fwhm-yC9 = backgroundC10 = airmass, irrelevant hereC11 = exposure time (subarray)> To model this light curve, I use a linear function of X, Y, fwhm-x, and fwhm-y, plus a transit model.<jupyter_code>bjd, flux, err, x, y, fwhm, fwhmx, fwhmy, bg, airmass, exptime = np.loadtxt('phot0002.txt', unpack=True)
bjd += 2450000
plt.errorbar(bjd, flux/np.median(flux), err, fmt='.', color='k', ms=1, ecolor='silver')

transit_model_g = transit_model(bjd, g)
oot = transit_model_g == 1

plt.plot(bjd, transit_model_g)
# plt.plot(bjd[oot], transit_model_g[oot], '.')
from toolkit import transit_duration

g.inc
X_all = np.vstack([x, y, fwhmx, fwhmy]).T

X = X_all[oot, :]

omega = np.diag(err[oot]**2)
omega_inv = np.linalg.inv(omega)

V = np.linalg.inv(X.T @ omega_inv @ X)
beta = V @ X.T @ omega_inv @ flux[oot]

regressed_lc = flux - (X_all @ beta) + 1
plt.plot(bjd, transit_model_g)

from scipy.optimize import fmin_powell

def minimize(p):
    return abs(np.sum((regressed_lc[oot] - transit_model_g[oot])**2 / 
                      (p[0] * err[oot])**2)/len(regressed_lc[oot]) - 1)

err_scale = fmin_powell(minimize, [1])

plt.errorbar(bjd, regressed_lc, err_scale*err, fmt='.')

np.savetxt('lightcurve.txt', np.vstack([bjd, regressed_lc, err_scale*err]).T)
plt.errorbar(bjd, flux, err, fmt='.', label='raw flux')
#plt.errorbar(bjd, regressed_lc, err_scale*err, fmt='.', label='Detrended flux')
plt.plot(bjd, (X_all @ beta), '.', label='detrending vector')
plt.legend()
plt.savefig('detrending.png')<jupyter_output><empty_output><jupyter_text>Compare with the lightcurve from Michael: ![](./gillon.png)<jupyter_code>def quadratic_to_nonlinear(u1, u2):
    a1 = a3 = 0
    a2 = u1 + 2*u2
    a4 = -u2
    return (a1, a2, a3, a4)

quadratic_to_nonlinear(*g.u)
import celerite
from celerite import terms
from scipy.optimize import minimize
from celerite.modeling import Model
from copy import deepcopy

original_params = g
times = bjd 
fluxes = regressed_lc
errors = err

class MeanModel3Param(Model):
    parameter_names = ['amp', 'depth', 't0']

    def get_value(self, t):
        params = deepcopy(trappist1('b'))
        params.rp = self.depth**0.5
        params.t0 = self.t0 + original_params.t0
        return self.amp * transit_model(t, params)

initp_dict = dict(amp=1, depth=original_params.rp**2,
                  t0=0)#t0=original_params.t0)

parameter_bounds = dict(amp=[0.9*np.min(fluxes), 1.3*np.max(fluxes)],
                        depth=[0.9 * original_params.rp**2,
                               1.1 * original_params.rp**2],
                        t0=[-0.05, 0.05])

mean_model = MeanModel3Param(bounds=parameter_bounds, **initp_dict)

bounds = dict(log_a=(-30, 30))#, log_c=(np.log(4), np.log(8)))

log_c_median = 1.98108915

kernel = terms.RealTerm(log_a=-2, log_c=log_c_median,
                        bounds=bounds)
kernel.freeze_parameter('log_c')

gp = celerite.GP(kernel, mean=mean_model, fit_mean=True)
gp.compute(times - original_params.t0, errors)

# Define a cost function
def neg_log_like(params, y, gp):
    gp.set_parameter_vector(params)
    return -gp.log_likelihood(y)

def grad_neg_log_like(params, y, gp):
    gp.set_parameter_vector(params)
    return -gp.grad_log_likelihood(y)[1]

# Fit for the maximum likelihood parameters
initial_params = gp.get_parameter_vector()
bounds = gp.get_parameter_bounds()
soln = minimize(neg_log_like, initial_params, #jac=grad_neg_log_like,
                method="L-BFGS-B", bounds=bounds, args=(fluxes, gp))
gp.set_parameter_vector(soln.x)

mu, var = gp.predict(fluxes, times - original_params.t0, return_var=True)
std = np.sqrt(var)

tmid = int(times.mean())
fig, ax = plt.subplots(2, 1, figsize=(6, 8), sharex=True)
ax[0].errorbar(times - tmid, fluxes, errors, fmt='.', color='k', ecolor='silver')
ax[0].fill_between(times - tmid, mu-std, mu+std, color='r', zorder=10, alpha=0.3)
ax[0].plot(times - tmid, mu, color='r', zorder=10)
ax[0].plot(bjd - tmid, transit_model_g)

ax[1].errorbar(times - tmid, fluxes - transit_model_g, errors, fmt='.', color='k', ecolor='silver')
ax[1].fill_between(times - tmid, mu-std-transit_model_g, mu+std-transit_model_g, color='r', zorder=10, alpha=0.3)
ax[1].plot(times - tmid, mu - transit_model_g, color='r', zorder=10)
ax[1].grid()
ax[1].set_xlabel('BJD - {0}'.format(tmid))
for axis in ax: 
    for j in ['right', 'top']:
        axis.spines[j].set_visible(False)

fig.tight_layout()
fig.savefig('gp.png', dpi=200, bbox_inches='tight')
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16100217864923747, 'fraction_numerical': 0.029411764705882353, 'mean_word_length': 4.31981460023175, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '495596', 'license': 'MIT', 'n_tokens_mistral': 1918, 'n_tokens_neox': 1818, 'n_words': 359, 'path': 'lightcurve.ipynb', 'repo_name': 'bmorris3/trappist1g_spots', 'size': '4590'}