import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

# Function to calculate y, given x and 2 model parameters
def func(x, norm, temp):
  y = norm / x**5 / (np.exp(1.435e8 / (x * temp)) - 1)
  return y

# Data points, hard-wired into the code
x = np.array([1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000,
              11000, 12000, 13000, 14000, 15000], dtype=np.float)
y = np.array([0.00, 1.62, 13.43, 25.37, 28.97, 26.91, 22.81, 18.54, 14.82,
              11.79, 9.39, 7.52, 6.06, 4.92, 4.03])

# Run curve-fit to find best model parameters
popt, pcov = curve_fit(func, x, y, p0=[1e22, 4000.])
norm = popt[0]
temp = popt[1]
print norm, temp

# Plot the results
plt.plot(x, y, 'ob')
plt.plot(x, func(x, norm, temp), '-g')
plt.show()