From the readme:
Sunspot umbral position and area information were digitized from the Mount Wilson daily white-light solar images some years ago (Howard, Gilman, and Gilman, 1984). These photographic images exist in a series that extends from 1917 through the present time. The digitized data extend from 1917 through 1985. Details about the observations, the measurement procedure, and the analysis techniques used earlier may be obtained from the earlier reference (Howard, Gilman, and Gilman, 1984).
These data were first published in Howard, Gilman, and Gilman, 1984. The data are accessible online here.
%matplotlib inline
import numpy as np
import matplotlib
matplotlib.rcParams['font.size'] = 14
import matplotlib.pyplot as plt
from glob import glob
paths = glob('/Users/bmmorris/data/Mt_Wilson_Tilt/*/gspot??.dat')
from astropy.time import Time
import astropy.units as u
import string
from astropy.table import Table
def split_interval(string, n, cast_to_type=float):
return [cast_to_type(string[i:i+n]) for i in range(0, len(string), n)]
base_time = Time('1915-01-01')
all_years_array = []
header = ("jd n_spots_leading n_spots_following n_spots_day_1 n_spots_day_2 "
"rotation_rate latitude_drift area_weighted_latitude_day_1 area_weighted_longitude_day_1 "
"area_weighted_longitude_day_2 area_day_1 area_day_2 tilt_day_1 delta_polarity_separation "
"area_weighted_longitude_day_1_leading area_weighted_longitude_day_1_following "
"area_weighted_latitude_day_1_leading area_weighted_latitude_day_1_following "
"area_leading area_following area_weighted_longitude_day_2_leading "
"area_weighted_longitude_day_2_following delta_tilt").split()
for path in paths:
f = open(path).read().splitlines()
n_rows = len(f) // 3
n_columns = 23
yearly_array = np.zeros((n_rows, n_columns))
for i in range(n_rows):
# First five ints specify time, afterwards specify sunspot data
int_list = split_interval(f[0+i*3][:18], 2, int)
month, day, year_minus_1900, hour, minute = int_list[:5]
year = year_minus_1900 + 1900
jd = Time("{year:d}-{month:02d}-{day:02d} {hour:02d}:{minute:02d}"
.format(**locals())).jd
row = [jd] + int_list[5:] + split_interval(f[1+i*3], 7) + split_interval(f[2+i*3][1:], 7)
yearly_array[i, :] = row
all_years_array.append(yearly_array)
table = Table(np.vstack(all_years_array), names=header)
WARNING: ErfaWarning: ERFA function "dtf2d" yielded 1 of "dubious year (Note 6)" [astropy._erfa.core]
table
jd | n_spots_leading | n_spots_following | n_spots_day_1 | n_spots_day_2 | rotation_rate | latitude_drift | area_weighted_latitude_day_1 | area_weighted_longitude_day_1 | area_weighted_longitude_day_2 | area_day_1 | area_day_2 | tilt_day_1 | delta_polarity_separation | area_weighted_longitude_day_1_leading | area_weighted_longitude_day_1_following | area_weighted_latitude_day_1_leading | area_weighted_latitude_day_1_following | area_leading | area_following | area_weighted_longitude_day_2_leading | area_weighted_longitude_day_2_following | delta_tilt |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 |
2421232.78125 | 1.0 | 1.0 | 2.0 | 3.0 | 13.151 | 0.417 | 13.142 | -11.249 | 0.314 | 19.615 | 13.818 | -25.35 | -0.372 | -10.382 | -11.392 | 13.542 | 13.076 | 2.778 | 16.837 | 0.919 | 0.269 | -2.74 |
2421232.78125 | 1.0 | 1.0 | 2.0 | 2.0 | 13.568 | 0.464 | 13.207 | -4.112 | 7.847 | 9.468 | 20.938 | 13.468 | 0.11 | -3.352 | -4.209 | 13.029 | 13.229 | 1.068 | 8.4 | 7.863 | 6.944 | -36.03 |
2421232.78125 | 1.0 | 3.0 | 4.0 | 12.0 | 12.176 | 1.184 | 13.705 | 18.61 | 29.247 | 46.216 | 58.587 | -2.183 | 0.829 | 19.616 | 17.678 | 13.742 | 13.67 | 22.228 | 23.988 | 30.838 | 28.035 | -0.8 |
2421233.83056 | 2.0 | 1.0 | 3.0 | 2.0 | 13.864 | 0.778 | 13.538 | 0.314 | 13.619 | 13.818 | 9.298 | -28.093 | 0.147 | 0.919 | 0.269 | 13.852 | 13.515 | 0.949 | 12.869 | 14.261 | 13.404 | 12.33 |
2421233.83056 | 1.0 | 1.0 | 2.0 | 3.0 | 14.073 | 0.093 | 13.648 | 7.847 | 21.368 | 20.938 | 10.639 | -22.565 | -0.434 | 7.863 | 6.944 | 13.654 | 13.283 | 20.572 | 0.366 | 21.397 | 21.17 | -43.1 |
2421233.83056 | 2.0 | 10.0 | 12.0 | 5.0 | 14.615 | 0.466 | 14.829 | 29.247 | 43.327 | 58.587 | 33.791 | -2.981 | -0.788 | 30.838 | 28.035 | 14.909 | 14.768 | 25.322 | 33.264 | 44.099 | 42.146 | -9.04 |
2421234.79792 | 0.0 | 0.0 | 1.0 | 4.0 | 13.76 | 0.132 | -22.199 | -28.179 | -15.905 | 3.542 | 2.573 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | -16.143 | 0.0 |
2421234.79792 | 1.0 | 1.0 | 2.0 | 3.0 | 14.281 | -0.386 | 14.341 | 13.619 | 26.394 | 9.298 | 20.175 | -15.763 | -0.251 | 14.261 | 13.404 | 14.516 | 14.282 | 2.331 | 6.967 | 26.754 | 26.136 | 4.41 |
2421234.79792 | 1.0 | 2.0 | 3.0 | 3.0 | 14.197 | -0.751 | 13.744 | 21.368 | 34.063 | 10.639 | 14.355 | -65.662 | 0.493 | 21.397 | 21.17 | 13.805 | 13.319 | 9.298 | 1.34 | 34.45 | 33.594 | 30.02 |
2421234.79792 | 2.0 | 3.0 | 5.0 | 7.0 | 14.066 | -0.506 | 15.31 | 43.327 | 55.895 | 33.791 | 79.416 | -12.021 | 1.774 | 44.099 | 42.146 | 15.469 | 15.068 | 20.437 | 13.355 | 58.337 | 55.392 | -27.64 |
... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
2446388.25486 | 1.0 | 1.0 | 2.0 | 3.0 | 14.494 | -0.733 | -8.952 | 21.784 | 35.237 | 1.748 | 0.15 | -27.515 | 0.925 | 21.816 | 21.414 | -8.969 | -8.762 | 1.608 | 0.139 | 35.482 | 35.114 | -47.08 |
2446388.25486 | 0.0 | 0.0 | 1.0 | 1.0 | 14.435 | 0.252 | 1.451 | -0.993 | 12.401 | 0.92 | 1.266 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
2446389.25 | 1.0 | 3.0 | 4.0 | 3.0 | 15.528 | -0.143 | -9.635 | 42.245 | 56.022 | 6.3 | 3.703 | -4.859 | 0.125 | 42.495 | 40.972 | -9.656 | -9.528 | 5.266 | 1.035 | 56.066 | 54.43 | 13.51 |
2446404.21042 | 3.0 | 1.0 | 4.0 | 6.0 | 15.263 | -0.414 | 19.574 | 26.951 | 41.31 | 2.359 | 4.32 | 30.273 | 0.329 | 29.087 | 25.792 | 18.399 | 20.211 | 0.829 | 1.53 | 42.915 | 39.0 | -10.74 |
2446413.32014 | 0.0 | 0.0 | 1.0 | 18.0 | 11.877 | -1.278 | 4.156 | -9.648 | 2.268 | 4.277 | 16.219 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | -0.601 | 0.0 |
2446413.32014 | 2.0 | 1.0 | 3.0 | 9.0 | 14.186 | -0.071 | -9.377 | -34.565 | -20.123 | 1.252 | 8.443 | -16.907 | 1.291 | -32.389 | -35.735 | -10.03 | -9.026 | 0.438 | 0.814 | -17.019 | -21.647 | 1.29 |
2446414.41389 | 6.0 | 12.0 | 18.0 | 17.0 | 13.829 | -0.557 | 2.757 | 2.268 | 15.65 | 16.219 | 9.338 | -25.858 | -0.389 | 4.293 | -0.601 | 3.738 | 1.369 | 9.507 | 6.712 | 17.701 | 12.998 | 4.56 |
2446414.41389 | 4.0 | 5.0 | 9.0 | 15.0 | 13.839 | -0.02 | -9.299 | -20.123 | -6.729 | 8.443 | 23.571 | -15.618 | -0.885 | -17.019 | -21.647 | -10.155 | -8.879 | 2.781 | 5.663 | -4.84 | -8.71 | 23.51 |
2446415.45556 | 11.0 | 6.0 | 17.0 | 11.0 | 14.713 | -0.722 | 2.177 | 15.65 | 26.664 | 9.338 | 24.766 | -21.296 | 0.958 | 17.701 | 12.998 | 2.976 | 1.144 | 5.267 | 4.072 | 29.399 | 23.858 | -1.34 |
2446417.19236 | 2.0 | 1.0 | 3.0 | 2.0 | 14.81 | 0.264 | 0.965 | 35.635 | 49.569 | 2.055 | 0.1 | 45.77 | -1.193 | 36.524 | 35.439 | 0.052 | 1.166 | 0.371 | 1.684 | 49.669 | 49.469 | -102.12 |
fig, ax = plt.subplots(1, 2, figsize=(12, 5), sharey=True)
ax[0].plot(table['area_weighted_longitude_day_1'],
table['area_weighted_latitude_day_1'], 'k.', alpha=0.2)
ax[1].hist(table['area_weighted_latitude_day_1'], 100,
histtype='stepfilled', color='k', orientation='horizontal')
ax[0].set(xlabel='Longitude', ylabel='Latitude')
ax[1].set(xlabel='Frequency')
fig.subplots_adjust(wspace=0.05)
plt.show()
plt.hist(table['area_day_1'], 50, log=True, histtype='stepfilled', color='k')
plt.xlabel('Spot group area [$\mu$Hem]')
plt.ylabel('Frequency')
plt.show()
img, xedges, yedges = np.histogram2d(table['area_day_1'],
np.abs(table['area_weighted_latitude_day_1']), 50)
plt.imshow(np.log(img), cmap=plt.cm.viridis, interpolation='nearest', origin='lower')
plt.xlabel('| Latitude |')
plt.ylabel('Area [$\mu$Hem]')
plt.show()
/Users/bmmorris/anaconda/lib/python3.5/site-packages/ipykernel/__main__.py:3: RuntimeWarning: divide by zero encountered in log app.launch_new_instance()
unique_dates = np.sort(np.unique(table['jd']))
fig = plt.figure(figsize=(8,8))
ax = fig.add_subplot(111)#, projection='hammer')
n_days = 50
cmap = lambda x: plt.cm.viridis(x/n_days)
for i, day in enumerate(unique_dates[:n_days]):
day = table[table['jd'] == day]
ax.plot(day['area_weighted_longitude_day_1'],
day['area_weighted_latitude_day_1'], 'o', color=cmap(i))
ax.set_xlim([-60, 60])
ax.set_ylim([-60, 60])
ax.grid()
ax.set(xlabel='longitude', ylabel='latitude')
[<matplotlib.text.Text at 0x11741ecf8>, <matplotlib.text.Text at 0x115393be0>]
area_per_spot = table['area_day_1'].data/table['n_spots_day_1'].data
not_inf = ~np.isinf(area_per_spot)
plt.hist(area_per_spot[not_inf], 100, log=True, histtype='stepfilled', color='k')
plt.xlabel('Mean area per spot group [$\mu$Hem]')
plt.ylabel('Frequency')
plt.show()
/Users/bmmorris/anaconda/lib/python3.5/site-packages/ipykernel/__main__.py:1: RuntimeWarning: divide by zero encountered in true_divide if __name__ == '__main__':
alledges = np.logspace(np.log10(1), np.log10(80-10), 100)
bin_centers = 0.5*(alledges[:-1] + alledges[1:])
n_spots_per_area_bin = np.zeros_like(bin_centers)
n_spots_per_area_bin_err = np.zeros_like(bin_centers)
for i in range(len(bin_centers)):
within_bin = ((alledges[i] < area_per_spot[not_inf]) &
(alledges[i+1] > area_per_spot[not_inf]))
n_spots_per_area_bin[i] = np.count_nonzero(within_bin)
n_spots_per_area_bin_err[i] = np.sqrt(n_spots_per_area_bin[i])
fig, ax = plt.subplots()
#ax.set_yscale("log")#, nonposy='clip')
ax.set_xscale("log")
ax.errorbar(bin_centers, n_spots_per_area_bin, n_spots_per_area_bin_err, fmt='.')
<Container object of 3 artists>
from scipy.optimize import leastsq
def fit_model(p, x):
a, x0, sigma = p
return a * np.exp(-0.5 * (x - x0)**2 / sigma**2)
log_bins = np.log(bin_centers)
def chi2(p, x):
return (fit_model(p, x) - n_spots_per_area_bin) / n_spots_per_area_bin_err
initp = [600, 1.0, 1.0]
bestp, success = leastsq(chi2, initp, args=(log_bins,))
fig, ax = plt.subplots()
ax.errorbar(log_bins, n_spots_per_area_bin, n_spots_per_area_bin_err, fmt='k.')
ax.plot(log_bins, fit_model(bestp, log_bins), color='r', lw=2)
ax.set(xlabel="log(Mean area per spot group [$\mu$Hem])", ylabel="log(frequency)")
[<matplotlib.text.Text at 0x113b89c50>, <matplotlib.text.Text at 0x116cdea58>]
n_samples = 100000
def draw_samples_from_area_dist(n_samples, lower_limit=1, upper_limit=10**4.5):
samples = []
while len(samples) < n_samples:
sample = np.random.lognormal(mean=bestp[1], sigma=bestp[2])
if sample > lower_limit and sample < upper_limit:
samples.append(sample)
return np.array(samples)
fig, ax = plt.subplots()
ax.hist(np.log(draw_samples_from_area_dist(n_samples)), 100)
ax.set(xlabel="log( Spot area [$\mu$Hem] )", ylabel="freq")
plt.show()
Area in a hemisphere:
$$A_{hemisphere} = 2\pi R_*^2$$Area in a spot in units of spot radii:
$$A_{spot} = \pi R_{spot}^2$$Area in a spot in units of hemispheres:
$$ A_{sp,hem} = \frac{\pi R_{spot}^2}{2\pi R_*^2} = \frac{1}{2} \left(\frac{R_{spot}}{R_*}\right)^2$$Radius of a spot in units of stellar radii given the spot area in hemispheres:
$$ R_{spot} / R_* = \sqrt{2 A_{sp,hem}} $$rspot_rstar = np.sqrt(1e-6 * 2 * np.exp(log_bins))
rp_rs = 0.058330305324663184
fig, ax = plt.subplots()
ax.errorbar(rspot_rstar / rp_rs, n_spots_per_area_bin, n_spots_per_area_bin_err, fmt='k.')
ax.plot(rspot_rstar / rp_rs, fit_model(bestp, log_bins), color='r', lw=2)
ax.set(xlabel="$R_{sp}/R_{planet}$", ylabel="frequency")
[<matplotlib.text.Text at 0x117003668>, <matplotlib.text.Text at 0x117c5cfd0>]
fig, ax = plt.subplots()
ax.errorbar(rspot_rstar, n_spots_per_area_bin, n_spots_per_area_bin_err, fmt='k.')
ax.plot(rspot_rstar, fit_model(bestp, log_bins), color='r', lw=2)
ax.set(xlabel="$R_{sp}/R_{*}$", ylabel="frequency")
[<matplotlib.text.Text at 0x117419eb8>, <matplotlib.text.Text at 0x11700f828>]
from astropy.constants import R_sun
import astropy.units as u
fig, ax = plt.subplots()
r_spot_km = (rspot_rstar * R_sun).to(1000 * u.km).value
ax.errorbar(r_spot_km, n_spots_per_area_bin, n_spots_per_area_bin_err, fmt='k.')
ax.plot(r_spot_km, fit_model(bestp, log_bins), color='r', lw=2)
ax.set(xlabel="$R_{sp}$ [1000 * km]", ylabel="frequency")
# for l in ax.get_xticklabels():
# l.set_rotation(45)
# l.set_ha('right')
[<matplotlib.text.Text at 0x11638ab00>, <matplotlib.text.Text at 0x1166a7978>]
latitude_grid = np.linspace(-45, 45, 80)
time_grid = Time(np.linspace(Time('1917-01-01').jd, Time('1985-01-01').jd, 200), format='jd')
butterfly_grid = np.zeros((len(latitude_grid)-1, len(time_grid)))
for i in range(len(time_grid)-1):
within_time_bin = ((table['jd'] > time_grid[i].jd) &
(table['jd'] < time_grid[i+1].jd))
spots_in_time_bin = table[within_time_bin]
freq, edges = np.histogram(spots_in_time_bin['area_weighted_latitude_day_1'], bins=latitude_grid)
butterfly_grid[:, i] = freq
fontsize = 14
fig, ax = plt.subplots(figsize=(14, 4))
cax = ax.imshow(np.log10(1 + butterfly_grid), cmap=plt.cm.viridis,
interpolation='nearest', origin='lower')
cbar = fig.colorbar(cax, pad=0.01)#, ticks=[-1, 0, 1])
cbar.ax.set_ylabel('$\log_{10} \, N_{spots}$', fontsize=fontsize)
xtickskip = 15
ax.set_xticks(list(range(len(time_grid)))[::xtickskip])
ax.set_xticklabels([i.strftime("%Y")
for i in time_grid[::xtickskip].datetime])
desired_yticks = np.arange(-45, 60, 15)
ax.set_yticks(np.interp(desired_yticks, latitude_grid, np.arange(len(latitude_grid))))
ax.set_yticklabels(desired_yticks)
ax.grid(color='silver', ls=':')
ax.set_xlabel('Year', fontsize=fontsize)
ax.set_ylabel('Latitude [degrees]', fontsize=fontsize)
fig.savefig("butterfly.png", bbox_inches='tight', dpi=200)
plt.show()
WARNING: ErfaWarning: ERFA function "dtf2d" yielded 1 of "dubious year (Note 6)" [astropy._erfa.core] WARNING: ErfaWarning: ERFA function "d2dtf" yielded 9 of "dubious year (Note 5)" [astropy._erfa.core]