import numpy as np
def diff_flux(energy, gamma):
return energy ** (- gamma)
def int_flux(energy_min, energy_max, gamma):
return (energy_max ** (1 - gamma) - energy_min ** (1 - gamma)) / (1 - gamma)
def energy_lafferty_power_law(energy_min, energy_max, spectral_index):
term0 = 1 - spectral_index
term1 = energy_max - energy_min
term2 = 1. / term0
flux_lw = term2 / term1 * (energy_max ** term0 - energy_min ** term0)
energy_lw = np.exp(- np.log(flux_lw) / spectral_index)
return energy_lw, flux_lw
for spectral_index in [1.1, 2, 3, 4]:
energy_min, energy_max = 1, 2
#energy_min, energy_max = 10, 11
flux_min, flux_max = diff_flux(energy_min, spectral_index), diff_flux(energy_max, spectral_index)
energy_lw, flux_lw = energy_lafferty_power_law(energy_min, energy_max, spectral_index)
print('spectral_index = {0}'.format(spectral_index))
print('energy:', energy_min, energy_max, energy_lw)
print('flux:', flux_min, flux_max, flux_lw)
print(flux_lw * (energy_max - energy_min), int_flux(energy_min, energy_max, spectral_index))
x_lw = np.log(energy_lw / energy_min) / np.log(energy_max / energy_min)
print('x:', x_lw)
spectral_index = 1.1 ('energy:', 1, 2, 1.4398202389003543) ('flux:', 1.0, 0.4665164957684037, 0.6696700846319253) (0.6696700846319253, 0.6696700846319253) ('x:', 0.52588870289285417) spectral_index = 2 ('energy:', 1, 2, 1.4142135623730949) ('flux:', 1.0, 0.25, 0.5) (0.5, 0.5) ('x:', 0.49999999999999994) spectral_index = 3 ('energy:', 1, 2, 1.3867225487012693) ('flux:', 1.0, 0.125, 0.375) (0.375, 0.375) ('x:', 0.47167916642628116) spectral_index = 4 ('energy:', 1, 2, 1.3607498666342404) ('flux:', 1.0, 0.0625, 0.29166666666666663) (0.29166666666666663, 0.2916666666666667) ('x:', 0.44440189466588803)
2 ** (1. / 3)
1.2599210498948732