Jump to content

File:Power spectrum of sunspot number, from 1945 to 2017.png

Page contents not supported in other languages.
This is a file from the Wikimedia Commons
From Wikipedia, the free encyclopedia

Original file (959 × 1,560 pixels, file size: 325 KB, MIME type: image/png)

Summary

Description
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()

```
Date
Source Own work
Author Cosmia Nebula

Licensing

I, the copyright holder of this work, hereby publish it under the following license:
w:en:Creative Commons
attribution share alike
This file is licensed under the Creative Commons Attribution-Share Alike 4.0 International license.
You are free:
  • to share – to copy, distribute and transmit the work
  • to remix – to adapt the work
Under the following conditions:
  • attribution – You must give appropriate credit, provide a link to the license, and indicate if changes were made. You may do so in any reasonable manner, but not in any way that suggests the licensor endorses you or your use.
  • share alike – If you remix, transform, or build upon the material, you must distribute your contributions under the same or compatible license as the original.

Captions

Add a one-line explanation of what this file represents

Items portrayed in this file

depicts

13 July 2023

image/png

File history

Click on a date/time to view the file as it appeared at that time.

Date/TimeThumbnailDimensionsUserComment
current04:37, 14 July 2023Thumbnail for version as of 04:37, 14 July 2023959 × 1,560 (325 KB)Cosmia NebulaUploaded while editing "Sunspot" on en.wikipedia.org

The following page uses this file:

Metadata