<jupyter_start><jupyter_text>Random sampling- $0<t<100$, a hundred points in all- $\tilde y=\sin(2\pi t)$, along with a Gaussian white noise (amplitude = 0.1)<jupyter_code>from astropy.timeseries import LombScargle
import astropy.units as u
def amplitude(t, factor=2):
    if t % 4 < 1:
        return 1
    else:
        return factor * 1
best_freq = []
for factor in np.linspace(1, 10, 100):
    rand = np.random.RandomState()
    t = np.sort(100 * rand.rand(100)) #np.linspace(0, 5, 1000) # 
    y = [amplitude(item, factor=3) for item in t] * np.sin(2 * np.pi * t) + 0.1 * rand.rand(100)
    frequency, power = LombScargle(t, y).autopower()
    best_freq.append(frequency[np.argmax(power)])
plt.figure(figsize=(8, 6))
plt.scatter(np.linspace(1, 10, 100), best_freq)
plt.xlabel('ellipsoidal factor')
plt.ylabel('best frequency')
rand = np.random.RandomState()
t = np.sort(100 * rand.rand(100)) #np.linspace(0, 5, 1000) # 
y = [amplitude(item, factor=3) for item in t] * np.sin(2 * np.pi * t) #+ 0.1 * rand.rand(100)
plt.figure(figsize=(12, 8))
plt.scatter(t, y)
frequency, power = LombScargle(t, y).autopower()
plt.figure(figsize=(12,8))
plt.plot(frequency, power)
frequency[np.argmax(power)]
fre = np.linspace(0.5,2.5,10000)
power = LombScargle(t, y).power(fre)
plt.figure(figsize=(12,8))
plt.plot(fre, power)
fre, power = LombScargle(t, y).autopower(nyquist_factor=4)
plt.figure(figsize=(12,8))
plt.plot(fre, power)
fre, power = LombScargle(t, y).autopower(
    minimum_frequency=0.5, maximum_frequency=2.5, samples_per_peak=10)
plt.figure(figsize=(12, 8))
plt.plot(fre, power)<jupyter_output><empty_output><jupyter_text>- Consider uncertainty (error) in observation, $\sigma=0.2$<jupyter_code>dy = np.random.randn(100)/5
LS = LombScargle(t, y,fit_mean=True)
fre, power = LS.autopower()
best_frequency = fre[power.argmax()]
best_frequency<jupyter_output><empty_output><jupyter_text>- Fit: $y = a + b\sin(2\pi \hat ft)+c\cos(2\pi \hat ft)$- Draw the periodogram<jupyter_code>best_frequency = fre[np.argmax(power)]
t_fit = np.linspace(0, 1 / best_frequency, 1000)
y_fit = LS.model(t_fit, best_frequency)
plt.figure(figsize=(8, 4))
plt.plot(t_fit, y_fit)
a, b, c = LS.model_parameters(best_frequency)

plt.scatter((t) % (1 / best_frequency), y)
best_frequency
offset = LS.offset()
design_matrix = LS.design_matrix(best_frequency, t_fit)
#np.allclose(y_fit, offset + design_matrix.dot(theta))
best_frequency = fre[np.argmax(power)]
t_fit = np.linspace(0, 1 / best_frequency * 2, 1000)
y_fit = LS.model(t_fit, best_frequency)
plt.figure(figsize=(8, 4))
plt.plot(t_fit, y_fit)
plt.errorbar((t) % (1 / best_frequency * 2), y, yerr=dy, fmt='o')
raw_g = np.loadtxt('ex_g')
raw_r = np.loadtxt('ex_r')
raw_i = np.loadtxt('ex_i')
def period(raw, col):
    raw = raw[raw[:,0].argsort()]
    raw[:,0] -= raw[0,0]
    LS = LombScargle(raw[:,0], raw[:,1], dy = raw[:,2], fit_mean=True)
    fre, power = LS.autopower(minimum_frequency=1e-10, maximum_frequency=50)
    best_frequency = fre[power.argmax()]
    
    plt.figure(figsize=(16,6), dpi=150)
    plt.subplot(121)
    plt.title('Power-frequency under LombScargle')
    plt.plot(fre, power, color=col)
    plt.xlabel('Frequency (Hz)', fontsize = 15)
    plt.ylabel('Power', fontsize = 15)
    
    plt.subplot(122)
    plt.title('Lightcurves')
    plt.errorbar(raw[:,0]%(1/best_frequency*2), raw[:,1], yerr = raw[:,2], color=col, fmt = 'o')
    plt.xlabel('Days', fontsize = 15)
    plt.ylabel('Magnitude', fontsize = 15)
    plt.gca().invert_yaxis()
    
    t_fit = np.linspace(0,1/best_frequency*2,1000)
    mag_fit = LS.model(t_fit,best_frequency)
    plt.plot(t_fit, mag_fit, color = 'black')
    
    print('Period: {:10.6f} Days'.format(1/best_frequency))
    print('Model parameters are: {}'.format(LS.model_parameters(best_frequency)))
    return(1/best_frequency)<jupyter_output><empty_output><jupyter_text>g band- Model parameters $\theta_n$:$$   y(t ; f, \vec{\theta})=\theta_{0}+\sum_{n=1}^{\text { nterns }}\left[\theta_{2 n-1} \sin (2 \pi n f t)+\theta_{2 n} \cos (2 \pi n f t)\right]$$<jupyter_code>pg = period(raw_g,'green')<jupyter_output>Period:   0.155596 Days
Model parameters are: [0.0089739  0.00913693 0.07261083]<jupyter_text>r band<jupyter_code>pr = period(raw_r,'red')<jupyter_output>Period:   0.155590 Days
Model parameters are: [0.01712923 0.04817438 0.17392221]<jupyter_text>i band<jupyter_code>pi = period(raw_i,'pink')
P = (pg+pr+pi)/3
zero_point = min(raw_g[0,0], raw_r[0,0], raw_i[0,0])
plt.figure(figsize=(12,8), dpi=250)
plt.title('Lightcurves')
plt.errorbar((raw_g[:,0]-zero_point)%(P*2), raw_g[:,1], yerr = raw_g[:,2], color='green', fmt = 'o', label = 'g')
plt.errorbar((raw_r[:,0]-zero_point)%(P*2), raw_r[:,1], yerr = raw_r[:,2], color='red', fmt = 'o', label = 'r')
plt.errorbar((raw_i[:,0]-zero_point)%(P*2), raw_i[:,1], yerr = raw_i[:,2], color='pink', fmt = 'o', label = 'i')
plt.xlabel('Days', fontsize = 15)
plt.ylabel('Magnitude', fontsize = 15)
plt.gca().invert_yaxis()
plt.legend()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1906093906093906, 'fraction_numerical': 0.054945054945054944, 'mean_word_length': 5.22636815920398, 'pattern_counts': {'":': 0, '<': 26, '<?xml version=': 0, '>': 23, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '133452', 'license': 'MIT', 'n_tokens_mistral': 2244, 'n_tokens_neox': 2059, 'n_words': 379, 'path': '_site/Docs/LombScargle.ipynb', 'repo_name': 'AstroJacobLi/astrojacobli.github.io', 'size': '5005'}