Fit blackbody spectrum to data in python

[Silviu]

Hi! I have to fit a blackbody spectrum to some data points. The y-axis is in mJy and the x-axis is in log_10(freq). My code looks like this:

from __future__ import division

import matplotlib.pyplot as plt

import numpy as np

from scipy.optimize import curve_fit

h = 6.63*10**(-34)

c = 3.0*10**8

k = 1.38*10**(-23)

pi = 3.14159265

scale_y = 1.5*10**13

def func(x,T):

f=scale_y*8*pi*h/(c**3)*(np.power(10,x))**3/(np.exp(h*np.power(10,x)/(k*T))-1)

return f

popt,pcov=curve_fit(func,frec, Flux)

print (popt)

freq and Flux are the lists with the x and y data and I need to find T. Just by eye T should be around 10000K. But when I run the program I get this error

RuntimeWarning: overflow encountered in exp
f=scale_y*8*pi*h/(c**3)*(np.power(10,x))**3/(np.exp(h*np.power(10,x)/(k*T))-1)

and I am not sure why, because for T=10000k the value of exp should be just fine. Any idea? Thank you!

 

[Silviu]

Hi! I have to fit a blackbody spectrum to some data points. The y-axis is in mJy and the x-axis is in log_10(freq). My code looks like this:

from __future__ import division

import matplotlib.pyplot as plt

import numpy as np

from scipy.optimize import curve_fit

h = 6.63*10**(-34)

c = 3.0*10**8

k = 1.38*10**(-23)

pi = 3.14159265

scale_y = 1.5*10**13

def func(x,T):

    f=scale_y*8*pi*h/(c**3)*(np.power(10,x))**3/(np.exp(h*np.power(10,x)/(k*T))-1)

    return f

popt,pcov=curve_fit(func,frec, Flux)

print (popt)

freq and Flux are the lists with the x and y data and I need to find T. Just by eye T should be around 10000K. But when I run the program I get this error

RuntimeWarning: overflow encountered in exp
f=scale_y*8*pi*h/(c**3)*(np.power(10,x))**3/(np.exp(h*np.power(10,x)/(k*T))-1)

and I am not sure why, because for T=10000k the value of exp should be just fine. Any idea? Thank you!

 

[jedishrfu]

Try breaking up the equation into smaller factors and print them out to see which ones may be the culprit.

Also you need to indent the equation and the subsequent return statement to be a part of the function you've defined. I did that too in your listing.

Lastly I added code tags i.e. (code=python) and (/code) around your code for syntax coloring. The parens are actually square brackets in the tags.

 

[Silviu]

Try breaking up the equation into smaller factors and print them out to see which ones may be the culprit.

Also you need to indent the equation and the subsequent return statement to be a part of the function you've defined. I did that too in your listing.

Lastly I added code tags i.e. (code=python) and (/code) around your code for syntax coloring. The parens are actually square brackets in the tags.

Thank you so much. I tried to break the equation in parts and the error seems to be in the np.exp() part, which I think it means that T goes to 0, but I don't know why.

 

[jedishrfu]

So now you should print out the factors in the np.exp() call to see if there is anything out of place.

Another thing to consider is that overflows and underflows can occur from intermediate results

As an example, consider a * b * (1/c) where a, b and c are big numbers

a * b--> causes an overflow exception

but if you restructure it to: ( a * (1/c) ) * b then the factor ( a * (1/c) ) creates a not so big number so that overflow might not occur when you finally multiply by b.

 

[mfb]

Did you check the argument of the exponential? What is the result of h*np.power(10,x)/(k*T)?

 

[rumborak]

Another option, since your function mostly involves multiplications and exponentiation, would be to operate in the log domain, I.e.take the logarithm of your input values, thus reducing all those operations to mostly additions or multiplications. At the end you exponentiate again to return to the original domain.
This trick is used a lot in machine learning where probabilities can be exceedingly low.