Given a certain numerical value for the energy density flux $F$ of a star, e.g. the solar constant $s_\odot = 1362~\frac{W}{m^2}$, and a certain wavelength range corresponding to the irradiated energy $(\lambda_{min},\lambda_{max})$, how can one find the corresponding black body temperature?

For this task, I aimed to find a solution in Mathematica which would allow me to fit the planck law function

\begin{equation} B_\lambda(T) = \frac{2\,h\,c^2}{\lambda^5}\,\frac{1}{\exp(\frac{h\,c}{\lambda\,k\,T})-1} \quad, \end{equation}

to the numerical output of the energy density $F$. Assuming the star to have an isotopically emitting surface (i.e. approximation as black body), the energy density flux for a given wavelength window can be computed via 1

\begin{equation} F = \pi \, \int^{\lambda_{max}}_{\lambda_{min}} B_\lambda(\color{red}{T})\,d\lambda \quad. \end{equation}

In Mathematica I have tried the following -- so far -- unfruitful approach (assuming the irradiation power of $F_{emphir} = 1340~\frac{W}{m2}$ was measured in the range of $(0.1~nm, 2300~nm)$:

  • First defining any variables

    refF = 1340.(*W/m^2*); c = 299792458(*m/s*); 
    h = 6.62607015*10^-34(*Js*); k = 1.380649*10^-23(*J/K*);
  • Defining the planck law function $B(\lambda,\,T)$and the energy density flux $F(T)$

    B[\[Lambda]_, T_] := (2*h*c^2)/\[Lambda]^5*1/(
    Exp[(h*c)/(\[Lambda]*k*T)] - 1);
    F[T_?NumericQ] := \[Pi]*NIntegrate[
    B[\[Lambda], T], {\[Lambda], 0.1*10^-9, 
  • Lastly, trying to find the corresponding effective temperature $T$ to the integrated area under the black body function, given a certain range of wavelengths from which the energy is emitted.

    FindFit[F, {a*
    F[T_], {0.99 < a < 1.01}}, a, T]

enter image description here

Unfortunately, the FindFit[data, expr, pars, vars] function does not work with a single numerical value as data and also needs parameters and values. In this particular case, the parameter I aim to tweak for the fit to be successful is also the variable. However, I did not find any other (for me obvious) helper operators in Mathematica which would solve a problem like this.

(I used a(n unnecessary) pre-factor $a$ to have both, variables and parameters in the FindFit operator.)

Hence, my question: How would one go about this in Mathematica?

Update #1

I now also tried to see if I can simply solve the problem through the operator NSolve but the system could not be solved given the methods available in this operator.

enter image description here

P.S. I also tried to directly integrate the planck law, given the effective temperature of the Sun for a certain wavelength range to see what wavelength range does contribute the most energy of the total solar energy flux, i.e. the solar constant. This at first did not work due to a wrong value of the speed of light or due to the extensive usage of subscript to label variables, which has now been amended. The cross-check calculation below yields $\approx 1313~\frac{W}{m^2}$, which -- given the large range of wavelength covered should be responsible for the emission of the bulk energy of the Sun -- already corresponds to the solar constant very well.

enter image description here

Footnotes & References:

1 Rybicki, G. B. & Lightman, A. P.: Radiative Processes in Astrophysics. Wiley-VCH Verlag, 2004, p. 7ff.

  • 4
    $\begingroup$ Among other issues, the speed of light is wrong by several orders of magnitude, the correct value is c=299 792 458 m/s. Also, using variables with subscripts is ill-advised (they tend to break things) which might be the source of the issues you see. $\endgroup$
    – Hans Olo
    Sep 6, 2022 at 9:59
  • $\begingroup$ I made an attempt at an answer, please check if this is what you want. $\endgroup$
    – Hans Olo
    Sep 6, 2022 at 14:23
  • $\begingroup$ Btw, FindFit is only suitable to fit a model to data, not solve equations. Using the help, eg ?FindFit you can get more info ;-) $\endgroup$
    – Hans Olo
    Sep 6, 2022 at 14:29
  • 1
    $\begingroup$ Thank you very much for your input, calling my attention to some overlooked mishaps and providing a valuable solution to my question! :) Oh, and thanks for the advice on the 'Help' Shortcut -- definitely helpful, as I hope to use Mathematica more frequently in the future. $\endgroup$ Sep 7, 2022 at 0:28

1 Answer 1


Now that the bugs seem to have been sorted out, I think the following should work. The point is to find the effective temperature, such that the flux at 1 AU matches the expected value. For this FindRoot[] is ideal, given our functions are numeric. Specifically, let's define the functions (and some new constants):

refF = 1340;
c = 299792458;
h = 6.62607015*10^-34; 
k = 1.380649*10^-23;
Rsun = 696340000;
AU = 149597870700;

where Rsun and AU are the radius of the Sun and 1 astronomical Unit (AU). Also,

B[\[Lambda]_,T_] := (2*h*c^2)/\[Lambda]^5*1/(Exp[(h*c)/(\[Lambda]*k*T)] - 1);
F[T_?NumericQ] := \[Pi]*NIntegrate[B[\[Lambda], T], {\[Lambda], 0.1*10^-9, 2300.0*10^-9}];

Then it is enough to say:

FindRoot[(Rsun/AU)^2 F[T] == refF, {T, 5800}]
(* {T -> 5808.6} *)

Which is reasonably close to the temperature of the Sun $(T_{\odot}\simeq 5778 K)$.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.