What is the RGB curve for blackbodies? I created a program to convert the temperature (in Kelvin) of a blackbody to RGB color. However, it is slightly inaccurate, and the deviations increase for values greater than 10000K and less than 1500K. Is there an accurate mathematical function that generates RGB values for given temperature, ranging from ~500K to 20000+ K?
Ok, here's my take on calculating the color of a blackbody, or any spectrum in fact:
Disclaimer: I'm not a color theorist, and there may be more accurate methods. But the result, shown in the bottom, looks about right.

Spectrum
First note that since color is a function of the relative intensity in various wavelength bands, it doesn't matter whether we express our spectrum in flux, flux density, intensity, luminance, or something else that only differs by constant values.
Let's consider a spectrum $S(\lambda)$ as a function of wavelength $\lambda$. If you're interested in a Planck spectrum, first calculate this from the temperature $T$:
$$
S(\lambda) = B(\lambda,T) = \frac{2hc^2}{\lambda^5} \frac{1}{e^{hc/\lambda k_\mathrm{B}T} - 1},
$$
where $h$, $c$, and $k_\mathrm{B}$ are Planck's constant, the speed of light, and Boltzmann's constant, respectively.
CIE 1931 color space
The CIE 1931 color space was the first attempt to quantify the relation between the distribution of wavelengths and the color perceived by humans. It makes use of the so-called color matching functions $\bar{x}$, $\bar{y}$, and $\bar{z}$, which are functions of $\lambda$, seen here:

Credit: Wikipedia
You can find tabulated values of these functions on the internet.
The color (which is 3D because humans have three different cones) of any spectrum can be mapped to a 2D position $\{x,y\}$ in the CIE 1931 color space chromaticity diagram seen here:

Credit: Wikipedia
Tristimulus values
The pair $\{x,y\}$ are normalized versions of the tristimulus values $\{X,Y,Z\}$ defined by
$$
X = \int S \, \bar{x} \, d\lambda\\
Y = \int S \, \bar{y} \, d\lambda\\
Z = \int S \, \bar{z} \, d\lambda.
$$
(Actually they're defined using the spectral radiance), but this is equal to the intensity to within a constant factor.)
The lower-case $\{x,y\}$ are then given by:
$$
\begin{array}{rcl}
x &amp; = &amp; \frac{X}{X+Y+Z}\\
y &amp; = &amp; \frac{Y}{X+Y+Z}\\
z &amp; = &amp; \frac{Z}{X+Y+Z} = 1-x-y,
\end{array}
$$
and the last of these relations is the reason only two values are needed.
XYZ to RGB
To go from $\{x,y\}$ or $\{X,Y,Z\}$ to RGB is a linear transformation which can be conducted by matrix multiplication. There are many variations here; my understanding is that is has to do with whatever device you're using to display you colors, but there might be more to it.
The sRGB color space assumes that the $\{X,Y,Z\}$ set has been scaled to $\{0.9505, 1, 1.0890\}$ for the CIE illuminant D65 (&quot;white&quot;); that is, $Y\equiv 1$.
With the above values, this new set can be calculated as
$$
\begin{array}{rcl}
Y &amp; \rightarrow &amp; 1\\
X &amp; \rightarrow &amp; x/y\\
Z &amp; \rightarrow &amp; z/y,
\end{array}
$$
and the RGB values are then, using the conversion matrix from here which cites this,
$$
\left[
\begin{array}{l}
R\\
G\\
B
\end{array}
\right]
=
\left[
\begin{array}{rrr}
 1.656492 &amp; -0.354851 &amp; -0.255038 \\
-0.707196 &amp;  1.655397 &amp;  0.036152 \\
 0.051713 &amp; -0.121364 &amp;  1.011530
\end{array}
\right]
\left[
\begin{array}{l}
X\\
Y\\
Z
\end{array}
\right]
$$
Gamma correction
I'm not 100% sure, but I think this set of RGB values needs to be gamma-corrected, which is a non-linear correction that account for the fact the the device displaying the colors shows different relative luminances than the human eye perceives. At least my result, shown in the bottom, looks about right, whereas it doesn't if I don't perform this correction.
The gamma correction is given by the function
$$
\gamma(u) = \left\{\begin{array}{ll}
12.92 u                 &amp; \mathrm{if\,} u \leq 0.0031308 \\
1.055 u^{1/2.4} - 0.055 &amp; \mathrm{otherwise}
\end{array}\right.
$$
and the perceived RGB values should then be
$$
\begin{array}{rcl}
R &amp; \rightarrow &amp; R\times\gamma(R)\\
G &amp; \rightarrow &amp; G\times\gamma(G)\\
B &amp; \rightarrow &amp; B\times\gamma(B)
\end{array}
$$
These values are mostly in the range $[0,1]$, but if they fall outside, they're usually clipped. If you want your values to be in the range $[0,255]$ you just multiply them by 255.
Result: the color of the Universe through time
I implemented the above algorithm in Python, as a part of a larger code called timeline that calculates various properties of the Universe as a function of time after the Big Bang. You can find it here on GitHub.
With that code I calculated the color of the Universe as a function of time after the Big Bang, and hence as a function of temperature. Until star formation kicks in around 200 million years after the Big Bang, the color is given by the temperature of the cosmic microwave background.
When stars arrive, the average cosmic spectrum becomes more complicated and can't be calculated (yet) from first principles, but must be observed. I haven't found such observations at early times, but Baldry et al (2002) presents the spectrum at various times after an age of 11 billion years, at which time the Universe is a dull &quot;cosmic latte&quot;. So in this figure, I grayed out the color between $t=200\,\mathrm{Myr}$ and $t=11\,\mathrm{Gyr}$:

But please note the disclaimer in the top of this answer :)