import numpy as np
import matplotlib.pyplot as plt

# 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])

# Guess parameters
guess_norm = 1.3e22
guess_temp = 5770
y_guess = func(x, guess_norm, guess_temp)

# Calculate and print the chi^2 value
chisq = np.sum((y - y_guess)**2)
print chisq

# Plot the data and the guess
plt.plot(x, y, 'ob')
plt.plot(x, y_guess, '-g')
plt.xlabel("Lambda (Angstrom)")
plt.ylabel("Flux (Arbitrary)")

# Plot the residual
for i in xrange(len(x)):
	plt.plot([x[i],x[i]], [y[i],y_guess[i]], '-r')

plt.show()