# Planck's Law 3D Plot with units outputting blank graph

I'm trying to plot Planck's Law for black body radiation on a wavelength-temperature graph where the output (z-axis) is the intensity. The function has multiple variables so I tried Plot3D.

https://en.wikipedia.org/wiki/Planck's_law#Different_forms

ClearAll["Global*"]
Subscript[B, 1][λ_, T_] :=
2 Quantity[1, "PlanckConstant"] Quantity[1, "SpeedOfLight"]^2 / Quantity[λ, "Meters"]^5 *
1/(Exp[Quantity[1, "PlanckConstant"] Quantity[1, "SpeedOfLight"] / (Quantity[λ, "Meters"] Quantity[1, "BoltzmannConstant"] Quantity[T, "Kelvins"])] - 1)

Plot3D[Subscript[B, 1][Quantity[λ, "Meters"], Quantity[T, "Kelvins"]],
{λ, 0, 3*10^(-6)}, {T, 0, 1000},
AxesLabel -> Automatic]


This is the first time I'm trying to use units with the Quantity function, but I can't seem to get the output to show anything...

Any help is appreciated as always!

The plot can be obtained with the following changes to the plot statement:

Plot3D[QuantityMagnitude@Subscript[B, 1][λ, T],
{λ, 0, 3*10^(-6)}, {T, 0, 1000},
AxesLabel -> Automatic]


However, it plots faster in the following version:

Substituting SI values for the quantities, and using simple variable names.

First obtaining the physical constant values.

entities = EntityValue[EntityList[
EntityClass["PhysicalConstant", "SIExact"]],
{"Symbol", "Name", "Value"}];

{h, c, k} = N[QuantityMagnitude@*Last /@ Lookup[Association[Thread[
entities[[All, 1]] -> entities]], {"h", "c", "k"}]]


{6.62607015*^-34, 2.99792458*^8, 1.380649*^-23}

Then using the OP's function with substitutions.

b[λ_, t_] := 2 h c^2/λ^5 * 1/(Exp[h c/(λ k t)] - 1)

Plot3D[b[λ, t], {λ, 0, 3*10^(-6)}, {t, 0, 1000},
AxesLabel -> Automatic]


You could use the built-in PlanckRadiationLaw:

Plot3D[PlanckRadiationLaw[Quantity[T, "Kelvins"],
Quantity[\[Lambda], "Meters"]], {\[Lambda], 10^(-9), 3*10^(-6)}, {T, 10^(-3), 1000},
AxesLabel -> Automatic, PlotRange -> All]
`