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)

2 ** (1. / 3)