Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Fit blackbody spectrum to data in python

  1. Apr 7, 2016 #1
    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!
     

    Attached Files:

  2. jcsd
  3. Apr 7, 2016 #2
    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:
    Code (Python):

    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!
     
    Last edited by a moderator: Apr 7, 2016
  4. Apr 7, 2016 #3

    jedishrfu

    Staff: Mentor

    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.
     
  5. Apr 7, 2016 #4
    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.
     
  6. Apr 7, 2016 #5

    jedishrfu

    Staff: Mentor

    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.
     
  7. Apr 7, 2016 #6

    mfb

    User Avatar
    2016 Award

    Staff: Mentor

    Did you check the argument of the exponential? What is the result of h*np.power(10,x)/(k*T)?
     
  8. Apr 7, 2016 #7
    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.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted



Similar Discussions: Fit blackbody spectrum to data in python
  1. Data fitting (Replies: 1)

Loading...