English: The daily sunspot number from 1945-01-02 to 2017-06-30, and its power spectral analysis. There are two prominent peaks corresponding to its 27-day cycle and 11-year cycle.

See "The 27-day signal in sunspot number series and the solar dynamo" (JL Le Mouël, MG Shnirman, EM Blanter - Solar Physics, 2007)

data from NOAA plotted in Python ```python import pandas as pd import matplotlib.pyplot as plt from urllib.request import urlopen from io import StringIO

1. 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') print("data retrieved")

1. Skip the header rows

data_lines = raw_data.split('\n')[3:]

1. Join the lines back into a single string and convert to a StringIO object

data = StringIO('\n'.join(data_lines))

1. Define the column names

columns = ['Year', 'Month', 'Day', 'SSN']

1. 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)

1. Create a 'Date' column from the 'Year', 'Month', and 'Day' columns

df['Date'] = pd.to_datetime(df'Year', 'Month', 'Day')

1. Set the 'Date' column as the index

df.set_index('Date', inplace=True)

1. Remove some entries where SSN does not exist. Fortunately, these all appear in the end
2. The data is very clean like that.

df = df.dropna(subset=['SSN'])

import numpy as np

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

1. Compute the one-dimensional n-point discrete Fourier Transform

yf = np.fft.fft(ssn)

1. Compute the power spectral density

power_spectrum = np.abs(yf[:n//2]) ** 2

1. Compute the frequencies associated with the elements of an FFT output

xf = np.fft.fftfreq(n, d=1) # Assume the time step d=1

1. Only take the first half of the frequencies (the non-negative part)

xf = xf[:n//2]

1. Create a figure with 2 subplots

fig, axes = plt.subplot_mosaic("A;B;B", figsize=(10, 16))

1. 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')

1. Plot the power spectrum in log-log scale in the second subplot

ax = axes["B"] ax.set_xscale('log') ax.set_yscale('log')

ax.scatter(xf, power_spectrum, s=1) num_segments = 100

1. Plot the power spectrum in log-log scale in the second subplot

num_segments = 100 min_log_freq = np.log(np.min(xf[xf > 0])) # Exclude zero frequency max_log_freq = np.log(np.max(xf)) log_freq_range = max_log_freq - min_log_freq

for i in range(num_segments):

   start = i * len(xf) // num_segments
   end = (i + 1) * len(xf) // num_segments
   avg_log_freq = np.mean(np.log(xf[start:end][xf[start:end] > 0]))  # Exclude zero frequency
   alpha = 1 - ((avg_log_freq - min_log_freq) / log_freq_range)  # Calculate alpha relative to log-frequency range
   alpha = alpha**2
   ax.loglog(xf[start:end], power_spectrum[start:end], color='b', alpha=alpha)

solar_cycle = (10.7 * 365) ax.vlines(1/solar_cycle, ymin=10**2, ymax=10**12, linestyle='--', linewidth=0.5) ax.text(1/solar_cycle, 10**12, 'Solar cycle = 10.7 years', rotation=0, verticalalignment='center', fontsize=12)

lunar_cycle = 27.3 ax.vlines(1/lunar_cycle, ymin=10**2, ymax=10**12, linestyle='--', linewidth=0.5) ax.text(1/lunar_cycle, 10**12, 'Monthly cycle = 27.3 days', rotation=0, verticalalignment='center', fontsize=12)

1. Filter out zero frequency

nonzero_indices = xf > 10**(-4) xf_nonzero = xf[nonzero_indices] power_spectrum_nonzero = power_spectrum[nonzero_indices]

1. Define the weights

weights = 1 / xf_nonzero

1. Calculate the logarithms of xf and power_spectrum

log_xf = np.log(xf_nonzero) log_power_spectrum = np.log(power_spectrum_nonzero)

1. Create the matrix A and vector b

A = np.vstack([log_xf, np.ones(len(log_xf))]).T b = log_power_spectrum

1. Solve the weighted least squares problem

x, residuals, rank, s = np.linalg.lstsq(A * np.sqrt(weights[:, np.newaxis]), b * np.sqrt(weights), rcond=None)

1. Extract the parameters from the solution

b_weighted = x[0] log_a_weighted = x[1] a_weighted = np.exp(log_a_weighted)

xf_fit = np.linspace(np.min(xf_nonzero), np.max(xf_nonzero), 1000) power_spectrum_fit_weighted = a_weighted * xf_fit ** b_weighted ax.loglog(xf_fit, power_spectrum_fit_weighted, 'r-', label=f'S(f) = {a_weighted:.2e} * f^{{{b_weighted:.2f}}}',

         color="black", linewidth=1, alpha=0.5, linestyle="--")

ax.legend()

ax.set_title('Power Spectrum of SSN') ax.set_xlabel('Frequency [day^(-1)]') ax.set_ylabel('Power Spectrum') ax.grid(True)

1. Adjust the spacing between subplots

plt.tight_layout()

plt.show()

```