pythonscipyscipy-optimizeastronomy

Why does my Python code using scipy.curve_fit() for Planck's Radiation Law produce 'popt=1' and 'pcov=inf' errors?

Covariance not estimated in SciPy's Curvefit

Here's my dataset:

frequency (Hz)  brightness (ergs/s/cm^2/sr/Hz)  brightness (J/s/m^2/sr/Hz)
float64 float64 float64
34473577711.372055  7.029471536390586e-16   7.029471536390586e-19
42896956937.69582   1.0253178228238486e-15  1.0253178228238486e-18
51322332225.44733   1.3544045476166584e-15  1.3544045476166584e-18
60344529880.18272   1.6902073280174815e-15  1.6902073280174815e-18
68767909106.5062    2.0125779972022745e-15  2.0125779972022743e-18
77780126454.10146   2.3148004995630144e-15  2.3148004995630145e-18
... ... ...
489996752265.52826  3.201319839821188e-16   3.201319839821188e-19
506039097962.6759   2.5968748350997043e-16  2.596874835099704e-19
523273092332.3638   2.0595903864583913e-16  2.0595903864583912e-19
539918248580.7806   1.7237876060575648e-16  1.7237876060575649e-19
557158231134.7507   1.3879848256567381e-16  1.3879848256567383e-19
573803387383.1646   1.0521820452559118e-16  1.0521820452559118e-19
591049358121.42 9.178609330955852e-17   9.178609330955852e-20

I tried to use CurveFit to fit this to Planck's Radiation Law:

import numpy as np
from scipy.optimize import curve_fit

h=6.626*10e-34
c=3*10e8
k=1.38*10e-23
const1=2*h/(c**2)
const2=h/k

def planck(x,v):
    return const1*(v**3)*(1/((np.exp(const2*v/x))-1))

popt,pcov= curve_fit(planck, cmb['frequency (Hz)'],cmb['brightness (J/s/m^2/sr/Hz)'])
print(popt, pcov)

Warning:

/tmp/ipykernel_2500/4072287013.py:11: RuntimeWarning: divide by zero encountered in divide
  return const1*(v**3)*(1/((np.exp((const2)*v/x))-1))

I get popt=1 and pcov=nan. Now the exponential term in the function differs by several orders of magnitude. And some of the values don't permit to approximate the law mathematically. I tried using the logarithmic form of the law but that doesn't work either. How can I overcome this problem?

Solution

  • A lot of problems here, including that your variables were swapped, you're needlessly redefining physical constants, and your expression was highly numerically unstable. You need to use exp1m instead:

    import matplotlib.pyplot as plt
import numpy as np
from scipy.constants import h, c, k
from scipy.optimize import curve_fit

freq, brightness_erg, brightness_j = np.array((
    (34473577711.372055, 7.0294715363905860e-16, 7.0294715363905860e-19),
    (42896956937.695820, 1.0253178228238486e-15, 1.0253178228238486e-18),
    (51322332225.447330, 1.3544045476166584e-15, 1.3544045476166584e-18),
    (60344529880.182720, 1.6902073280174815e-15, 1.6902073280174815e-18),
    (68767909106.506200, 2.0125779972022745e-15, 2.0125779972022743e-18),
    (77780126454.101460, 2.3148004995630144e-15, 2.3148004995630145e-18),
    (489996752265.52826, 3.2013198398211880e-16, 3.2013198398211880e-19),
    (506039097962.67590, 2.5968748350997043e-16, 2.5968748350997040e-19),
    (523273092332.36380, 2.0595903864583913e-16, 2.0595903864583912e-19),
    (539918248580.78060, 1.7237876060575648e-16, 1.7237876060575649e-19),
    (557158231134.75070, 1.3879848256567381e-16, 1.3879848256567383e-19),
    (573803387383.16460, 1.0521820452559118e-16, 1.0521820452559118e-19),
    (591049358121.42000, 9.1786093309558520e-17, 9.1786093309558520e-20),
)).T


def planck(v: np.ndarray, T: float) -> np.ndarray:
    return 2*h/c/c * v**3 / np.expm1(h*v/k/T)

guess = 2.5,

(T,), _ = curve_fit(
    f=planck, xdata=freq, ydata=brightness_j, p0=guess, method='lm',
    # bounds=(0.1, np.inf),
)
print('T =', T)

fig, ax = plt.subplots()
v = np.linspace(freq.min(), freq.max(), 500)
ax.scatter(freq, brightness_j, label='data')
ax.plot(v, planck(v, *guess), label='guess')
ax.plot(v, planck(v, T), label='fit')
ax.legend()
plt.show()

    fit