I need to compare a model to data to find the best fitting input. My data is surface reflectance of a body (= scattered light + thermal emission) at different wavelengths, and I simulate the reflectance by modeling the thermal emission part with Planck's law. I need to find the temperature at which the modelled emission (and thus reflectance) fits the observations.
To do it, I define Planck's law as:
h = 6.626070e-34 # the Planck constant
c = 2.997924e+8 # the speed of light in vacuum
k = 1.380649e-23 # the Boltzmann constant
def planck(T): # Planck's law for black body radiation
intensity = (2*h*c**2) / ( ((wvleng*1e-9)**5) * (np.exp(h*c/((wvleng*1e-9)*k*T)) - 1.0) )
thermal = 0.95 * (intensity*np.pi*D**2)/(SSI*1e9)
return thermal
which returns the thermal emission of the body at different wavelengths. Then, I defined the function that will qualify the fit (sum of squared difference between model and data at each wavelength):
def residuals(T): #
S = scat_light + 0.95*planck(T) # simulated reflectance
return np.sum((reflectance - S)**2)
where T is the black body temperature. Finally, I use optimize.minimize
to find which T leads to the best fit between S (simulated reflectance) and data reflectance:
x0 = 345.0 # initial guess, temperature is the only variable
res = scipy.optimize.minimize(residuals, x0, method='SLSQP') # optimization
fitted_temperature = res.x
Here's the problem: it gives me a temperature way below (~220K) than what it should be (around 340K), not fitting the data at all. Does anyone have an idea about why it doesn't converge properly? Thanks in advance.
According to the
scipy.optimize.minimize
documentation x0 should be, "Initial guess. Array of real elements of size (n,), where ānā is the number of independent variables." So you're guessing that a temperature of T=345 will minimize the sum of squared differences. This is then returning a single best T for all T's in your data (assuming you have more than one sample) to minimize this function. Since apparently you know the answer should be around 340K I would imagine that there is an error in your formulas leading into or out of theplanck()
function somewhere.It also wouldn't hurt to look at the other attributes of the
res
output to see if there can be additional clues garnered there.