I'm new to using mathematica, but I want to do a plot of this equation below:

$$ F(\nu) = \frac{4 \pi h \nu^{4}}{c^{2}} \int_{r_{i}}^{r_{o}} \frac{r d r}{\exp \left[\left(\frac{3 \pi}{2}\right)^{1 / 4} \frac{h \nu}{k T_{*}}\left(\frac{r}{R_{*}}\right)^{3 / 4}\right]-1} $$ where $R_* = 2.5(696.34*10^6)$, $r_i = 6*R_*$. $r_0 = 7*10^4*R_*$, $T=4000$, $h$ is de Planck constant, $c$ is speed of light and $k$ is the Bolztmann constant. The plot should be from $\nu_1 = 10^{10}Hz$ to $\nu_2 = 10^{16}Hz$. All values here are with the right units when using SI, so don't bother with it.

My main problem is that I do not know how can I set a function like this inside Mathematica because it is a function of $\nu$ and it has an integral in it. I know how to do integrals, but in this case, I have to compute the integral to each one of the $\nu$ values of my domain. I have no clue where to begin. Can someone help me?

Oh, it is supposed to look like this

EDIT: I followed the tip from the first comment but now I have another problem: The mathematica does not compute my integral. As it is indicate in the comment, I did:

fun[a_?NumericQ] :=  NIntegrate[((9.3*(10^(-50))) (a^4)*x)/(e^(6.89*(10^(-24))*a*x^(3/4)) - 1), {x,6*2.5*(696.34*(10^6)), 7*(10^4)*2.5*(696.34*(10^6))}]; Plot[fun[a], {a, 10^10, 10^16}]

but I get a black graph. It got me back a message saying:

The integrand has evaluated to non-numerical values for all sampling points in the region with boundaries

"Further output of NIntegrate::inumr will be suppressed during this calculation"

  • 1
    $\begingroup$ fun[a_?NumericQ] := NIntegrate[a*x, {x, 0, 1}] ; Plot[fun[a], {a, 1, 5}] $\endgroup$
    – I.M.
    May 7, 2021 at 2:42
  • $\begingroup$ Thanks! It was a valuable tip! $\endgroup$
    – kplt
    May 7, 2021 at 3:18
  • $\begingroup$ Post the formulas by Mathematica code. $\endgroup$
    – cvgmt
    May 7, 2021 at 6:36

1 Answer 1


R = 25 (69634*^3);
ri = 6 R;
r0 = 7*^4 R;
T = 4000;

h = QuantityMagnitude[Quantity[1, "PlanckConstant"] // UnitConvert];

c = QuantityMagnitude[Quantity[1, "SpeedOfLight"] // UnitConvert];

k = QuantityMagnitude[Quantity[1, 
       "BoltzmannConstant"] // UnitConvert];

F[v_?NumericQ] := 4 Pi h v^4/c^2 * NIntegrate[
  r/(Exp[(3 Pi/2)^(1/4) h v/(k T) (r/R)^(3/4)] - 1),
   {r, ri, r0}]

LogLogPlot[F[v], {v, 10^10, 10^16}]

enter image description here


Plot[Log10[F[10^v]], {v, 10, 16},
 Frame -> True,
 FrameLabel -> {"log v  (Hz)", "log F(v)"}]

enter image description here


Your Answer

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

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