<jupyter_start><jupyter_text>Reproduce [Chaplin 2010](https://ui.adsabs.harvard.edu/abs/2010ApJ...713L.169C/abstract) Figure 1<jupyter_code>%matplotlib inline
import matplotlib.pyplot as plt
import kplr
import numpy as np
import sys
sys.path.insert(0, '../')
from astropy.io import fits

# data = fits.getdata('ftp://ftp.pmodwrc.ch/pub/data/irradiance/virgo/'
#                     '1-minute_Data/VIRGO_1min_0083-7404.fits', cache=False)

data = fits.getdata('../data/VIRGO_1min_0083-7404.fits.gz', cache=False)
from shocksgo import interpolate_missing_data
import numpy as np
time = np.arange(len(data))
#flux = data[data != 99]

times, fluxes = interpolate_missing_data(time[data != -99], data[data != -99])

fluxes /= np.median(fluxes)
plt.plot(times, fluxes)
len(fluxes)//20
#from scipy.signal import periodogram
import os

from shocksgo import power_spectrum
freqs, powers = power_spectrum(fluxes[:len(fluxes)//10], d=60)
freqs *= 1e6
# periodogram_path = 'periodogram.npy'

# if os.path.exists(periodogram_path):
#     freqs, powers = np.load(periodogram_path)
# else: 
#     freqs, powers = periodogram(fluxes, fs=1/60)
#     freqs *= 1e6
#     np.save(periodogram_path, np.vstack([freqs, powers]))

plt.loglog(freqs, powers, marker=',', lw=0)
plt.xlim([0.1, 1e4])
plt.ylim([1e-10, 1e-2])
from scipy.stats import binned_statistic

cutoff_freq = 1e5
bs = binned_statistic(np.log(freqs[freqs != 0]), powers[freqs != 0], statistic=np.nanmedian, bins=1000)
bincenters = 0.5 * (bs.bin_edges[:-1] + bs.bin_edges[1:])
binned_power = bs.statistic[np.exp(bincenters) < cutoff_freq]
binned_freq = np.exp(bincenters)[np.exp(bincenters) < cutoff_freq]
plt.loglog(freqs, powers, ',', alpha=0.5)
plt.loglog(binned_freq, binned_power)
plt.xlabel('Freq [$\mu$Hz]')
plt.ylabel('Power')
plt.xlim([1e2, 6e3])
plt.ylim([1e-10, 1e-4])
from scipy.ndimage import gaussian_filter1d

plt.semilogy(freqs, powers, ',', alpha=0.5)
# plt.semilogy(freqs[np.argsort(freqs)], gaussian_filter1d(powers[np.argsort(freqs)], 100))
plt.semilogy(binned_freq, binned_power)

plt.xlabel('Freq [$\mu$Hz]')
plt.ylabel('Power')
plt.xlim([2000, 4000])
plt.ylim([1e-10, 1e-4])
from shocksgo import generate_solar_fluxes
from astropy.constants import M_sun, L_sun
import astropy.units as u

times, fluxes, kernel = generate_solar_fluxes(duration=10*u.min)
psd = kernel.get_psd(2*np.pi*freqs*1e-6) /2/np.pi
fig, ax = plt.subplots(1, 2, figsize=(8, 2.5))
# ax[0].loglog(binned_freq, binned_power)
ax[0].loglog(freqs, powers, marker=',', lw=0, alpha=0.3, rasterized=True, color='k')

ax[0].loglog(freqs, psd, color='r')
ax[0].set_xlim([1e-2, 1e4])
ax[0].set_ylim([1e-10, 1e-2])
ax[0].set_xlabel('Frequency [$\mu$Hz]')
ax[0].set_ylabel('Power [(flux)$^2$/Hz]')

# ax[1].semilogy(binned_freq, binned_power)
ax[1].semilogy(freqs, powers, marker=',', lw=0, alpha=0.3, rasterized=True, color='k')
ax[1].semilogy(freqs, psd, color='r')
ax[1].set_xlim([2000, 4000])
ax[1].set_ylim([1e-10, 1e-6])
ax[1].set_xlabel('Frequency [$\mu$Hz]')
ax[1].set_ylabel('Power [(flux)$^2$/Hz]')

for s in ['right', 'top']:
    for axis in ax: 
        axis.spines[s].set_visible(False)

fig.tight_layout()
fig.suptitle('Sun (SOHO/VIRGO SPM)', va='bottom')
fig.savefig('paper_plots/sun.pdf', bbox_inches='tight', dpi=300)
durations = []
for i in np.arange(2, 7):
    duration = %timeit -o -r 1 generate_solar_fluxes(duration=10**i * u.min)
    durations.append(duration)
number_points = 10**np.arange(2, 7)
duration_best = np.array([duration.best for duration in durations])

plt.figure(figsize=(3, 3))
plt.loglog(number_points, duration_best, color='k')
plt.xlabel('Sim. Duration [min]')
plt.ylabel('Runtime [s]')
plt.grid(ls=':')
for s in ['right', 'top']:
    plt.gca().spines[s].set_visible(False)
plt.savefig('paper_plots/runtime.pdf', bbox_inches='tight')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.18697698744769875, 'fraction_numerical': 0.04602510460251046, 'mean_word_length': 5.746031746031746, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 5, 'https://': 1, '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': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '447462', 'license': 'MIT', 'n_tokens_mistral': 1756, 'n_tokens_neox': 1592, 'n_words': 275, 'path': 'notebooks/validation_sun.ipynb', 'repo_name': 'bmorris3/shocksgo', 'size': '3824'}