import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from urllib.request import urlopen
from io import StringIO
from scipy.io.wavfile import write

# Fetch the data
url = "https://ngdc.noaa.gov/stp/space-weather/solar-data/solar-indices/sunspot-numbers/american/lists/list_aavso-arssn_daily.txt"
raw_data = urlopen(url).read().decode('utf-8')

# Skip the header rows
data_lines = raw_data.split('\n')[3:]

# Join the lines back into a single string and convert to a StringIO object
data = StringIO('\n'.join(data_lines))

# Define the column names
columns = ['Year', 'Month', 'Day', 'SSN']

# Read the data into a pandas DataFrame, specifying the delimiter as any number of spaces
df = pd.read_csv(data, delim_whitespace=True, names=columns)

# Create a 'Date' column from the 'Year', 'Month', and 'Day' columns
df['Date'] = pd.to_datetime(df[['Year', 'Month', 'Day']])

# Set the 'Date' column as the index
df.set_index('Date', inplace=True)

# Remove some entries where SSN does not exist. Fortunately, these all appear in the end
# The data is very clean like that.
df = df.dropna(subset=['SSN'])

# Generate a noise matrix
noise_matrix = np.random.normal(size=(len(df), 100))

# Multiply the SSN series with the matrix (broadcasting)
output_matrix = df['SSN'].values[:, np.newaxis] * noise_matrix

# Flatten the matrix
output_flattened = output_matrix.flatten()

# Convert it into a wav file with 36000 sampling rate.
loudness = 0.01
write("sunspot_sound.wav", 36000, output_flattened * loudness)

ssn = df['SSN']
n = len(ssn)

# Compute the one-dimensional n-point discrete Fourier Transform
yf = np.fft.fft(ssn)
# Compute the power spectral density
power_spectrum = np.abs(yf[:n//2]) ** 2

# Compute the frequencies associated with the elements of an FFT output
xf = np.fft.fftfreq(n, d=1) # Assume the time step d=1

# Only take the first half of the frequencies (the non-negative part)
xf = xf[:n//2]

# Create a figure with 2 subplots
fig, axes = plt.subplot_mosaic("A", figsize=(10, 6))

# Plot SSN over time in the first subplot
ax = axes["A"]
ax.scatter(df.index, df['SSN'], s=0.5)
ax.set_title('Sunspot Number Over Time')
ax.set_xlabel('Date')
ax.set_ylabel('SSN')

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
from matplotlib.animation import FuncAnimation

# Assuming 'df' is your DataFrame with sunspot data

ssn = df['SSN']

# Create a figure with a single subplot
fig, ax = plt.subplots(figsize=(10, 6))

# Convert dates to numeric format for the plot
dates = mdates.date2num(df.index.to_pydatetime())

# Create a scatter plot object
scat = ax.scatter([], [], s=0.5)

# Add a vertical line object
vline = ax.axvline(x=0, color='r')

# Set labels
ax.set_title('Sunspot Number Over Time')
ax.set_xlabel('Date')
ax.set_ylabel('SSN')

# Set the plot limits
ax.set_xlim(dates.min(), dates.max())
ax.set_ylim(ssn.min(), ssn.max())

# Format the x-axis to display dates
ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d'))

# Animation update function
def update(frame):
    # Convert frame number to date
    date = df.index.min() + pd.Timedelta(days=18*frame)
    date_num = mdates.date2num(date.to_pydatetime())

    # Update the scatter plot
    scat.set_offsets(np.c_[dates[dates <= date_num], ssn[df.index <= date]])

    # Update the vertical line
    vline.set_xdata(date_num)

# Create the animation
ani = FuncAnimation(fig, update, frames=range(len(df)//18), interval=50)

# Save the animation
ani.save('sunspot_animation.webm', fps=20, writer='ffmpeg', codec='vp9')

plt.close(fig)  # Close the figure

!ffmpeg -i sunspot_animation.webm -i sunspot_sound.wav -c:v copy -c:a libvorbis -af "volume=0.4" output.webm