73 lines
2.2 KiB
Python
73 lines
2.2 KiB
Python
import numpy as np
|
|
from matplotlib import pyplot as plt
|
|
from mpl_toolkits.basemap import Basemap
|
|
|
|
|
|
def data_read(file):
|
|
with open(file, 'r') as f:
|
|
data = f.readlines()
|
|
|
|
row_header = int(data[0].split(',')[1].strip())
|
|
row_variable = int(data[0].split(',')[5].strip())
|
|
row_column = int(data[0].split(',')[6].strip())
|
|
|
|
var = data[row_header - row_variable:row_header]
|
|
data = [tuple(s.replace('\n', '').strip() for s in row.split(',')) for row in data[row_header:-1]]
|
|
|
|
dtype = []
|
|
for row in var:
|
|
row = [s.replace("'", '').strip() for s in row.split(',')]
|
|
num = int(row[2])
|
|
type_s = 'float' if num > 0 else 'string'
|
|
for k in range(num):
|
|
name_s = row[0] + ('_{:d}'.format(k + 1) if num > 1 else '')
|
|
dtype.append((name_s, type_s))
|
|
|
|
assert len(dtype) == row_column
|
|
|
|
data = np.array(data, dtype=dtype)
|
|
|
|
return dtype, data
|
|
|
|
|
|
def orbit():
|
|
_, data = data_read('assets/orbit.csv')
|
|
lat, lon = data['Latitude'], data['Longitude'] - 180
|
|
|
|
m = Basemap()
|
|
m.drawcoastlines()
|
|
m.fillcontinents(color='white', lake_color='lightskyblue')
|
|
m.drawmapboundary(fill_color='skyblue')
|
|
m.scatter(lon, lat, s=1)
|
|
plt.savefig('docs/orbit.png', bbox_inches='tight')
|
|
|
|
|
|
def GCR(M=1):
|
|
_, data = data_read('assets/CREME86-M{:d}.csv'.format(M))
|
|
E = data['Energy']
|
|
proton_i, proton_d = data['IFlux_1'], data['DFlux_1']
|
|
alpha_i, alpha_d = data['IFlux_2'], data['DFlux_2']
|
|
|
|
fig, ax1 = plt.subplots(1, 1)
|
|
ax1.plot(E, proton_i, label=r'$p$')
|
|
ax1.plot(E, alpha_i, label=r'$\alpha$')
|
|
ax1.set_yscale('log')
|
|
ax1.set_ylabel(r'$Integrated\ Flux\ (\mathrm{m^{-2}sr^{-1}s^{-1}})$')
|
|
ax1.legend(loc="upper left")
|
|
|
|
ax2 = ax1.twinx()
|
|
ax2.plot(E, proton_d, linestyle=':', label=r'$p$')
|
|
ax2.plot(E, alpha_d, linestyle=':', label=r'$\alpha$')
|
|
ax2.set_yscale('log')
|
|
ax2.set_ylabel(r'$Differential\ Flux\ (\mathrm{m^{-2}sr^{-1}s^{-1}MeV^{-1}})$')
|
|
ax2.legend(loc="upper right")
|
|
|
|
ax1.set_xlabel(r'$Energy\ (\mathrm{MeV})$')
|
|
ax1.set_xlim(1, 1e5)
|
|
plt.title('Average spectra - M = {:d}'.format(M))
|
|
plt.xscale('log')
|
|
plt.savefig('docs/spectra-M{:d}'.format(M))
|
|
|
|
GCR()
|
|
GCR(3)
|