calculate the rgb of a very hot black-body

infinite_temperature.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "b5be4f28-d219-40c4-bc68-35fb3158af86",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from mpmath import mp, mpf, exp\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "ea9d967a-b90d-4d0b-928b-a501dae5364e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set the precision to 10 decimal places\n",
    "mp.dps = 10\n",
    "\n",
    "# High precision constants for Planck's law using exact values\n",
    "h = mpf(\"6.62607015e-34\")  # Planck constant (J·s)\n",
    "c = mpf(\"299792458\")  # Speed of light in vacuum (m/s)\n",
    "k = mpf(\"1.38064852e-23\")  # Boltzmann constant (J/K)\n",
    "\n",
    "# Load CIE XYZ data\n",
    "cie_xyz_data = pd.read_csv(\"lin2012xyz2e_fine_7sf.csv\")\n",
    "wavelengths = cie_xyz_data.iloc[:, 0].values * 1e-9  # Convert nm to meters\n",
    "X_values = cie_xyz_data.iloc[:, 1].values\n",
    "Y_values = cie_xyz_data.iloc[:, 2].values\n",
    "Z_values = cie_xyz_data.iloc[:, 3].values\n",
    "\n",
    "# Transformation matrices for RGB to XYZ and XYZ to RGB\n",
    "matrix_XYZ_to_RGB = np.array(\n",
    "    [\n",
    "        [3.240969941904523, -1.537383177570094, -0.498610760293003],\n",
    "        [-0.969243636280880, 1.875967501507721, 0.041555057407176],\n",
    "        [0.055630079696994, -0.203976958888977, 1.056971514242879],\n",
    "    ]\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "b37caf58-416f-401a-bf86-dee417b65905",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_rgb_square(rgb):\n",
    "    \"\"\"Plots a square filled with the specified RGB color and labels it with the RGB values.\"\"\"\n",
    "    normalized_rgb = [x / 255.0 for x in rgb]\n",
    "    square = np.ones((10, 10, 3)) * normalized_rgb\n",
    "    plt.imshow(square)\n",
    "    plt.axis(\"off\")\n",
    "    rgb_text = f\"rgb({rgb[0]}, {rgb[1]}, {rgb[2]})\"\n",
    "    plt.text(\n",
    "        5,\n",
    "        5,\n",
    "        rgb_text,\n",
    "        color=\"white\" if sum(rgb) < 400 else \"black\",\n",
    "        fontsize=12,\n",
    "        ha=\"center\",\n",
    "        va=\"center\",\n",
    "        weight=\"bold\",\n",
    "    )\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "def blackbody_spectrum(wavelengths, temperature):\n",
    "    \"\"\"Calculate the blackbody spectrum using mpmath for arbitrary precision.\"\"\"\n",
    "    return (2 * h * c**2 / wavelengths**5) / (\n",
    "        exp(h * c / (wavelengths * k * temperature)) - 1\n",
    "    )\n",
    "\n",
    "\n",
    "def temperature_to_xyz(temperature):\n",
    "    \"\"\"Convert a temperature in Kelvin to CIE XYZ values using mpmath.\"\"\"\n",
    "    spectrum = np.array([blackbody_spectrum(w, temperature) for w in wavelengths])\n",
    "\n",
    "    # Integrate the spectrum with the CIE XYZ color matching functions\n",
    "    X = np.trapz(spectrum * X_values, wavelengths)\n",
    "    Y = np.trapz(spectrum * Y_values, wavelengths)\n",
    "    Z = np.trapz(spectrum * Z_values, wavelengths)\n",
    "\n",
    "    # Normalize by Y to obtain chromaticity values\n",
    "    XYZ = np.array([X, Y, Z])\n",
    "    XYZ_normalized = XYZ / XYZ.sum()\n",
    "\n",
    "    return XYZ_normalized\n",
    "\n",
    "\n",
    "def xyz_to_rgb(XYZ):\n",
    "    \"\"\"Convert CIE XYZ values to linear RGB using the transformation matrix.\"\"\"\n",
    "    RGB = np.dot(matrix_XYZ_to_RGB, XYZ)\n",
    "\n",
    "    # Clip and normalize RGB values\n",
    "    RGB = np.clip(RGB, 0, None)\n",
    "    RGB /= RGB.max()\n",
    "\n",
    "    return RGB\n",
    "\n",
    "\n",
    "def gamma_correct(rgb):\n",
    "    \"\"\"Apply gamma correction based on the sRGB transfer function.\"\"\"\n",
    "    corrected_rgb = np.where(\n",
    "        rgb <= 0.0031308, 12.92 * rgb, 1.055 * np.power(rgb, 1 / 2.4) - 0.055\n",
    "    )\n",
    "    return np.clip(corrected_rgb, 0, 1)\n",
    "\n",
    "\n",
    "def temperature_to_rgb(temperature):\n",
    "    \"\"\"Convert a temperature in Kelvin to RGB color using mpmath and the sRGB transfer function.\"\"\"\n",
    "    XYZ = temperature_to_xyz(temperature)\n",
    "    linear_rgb = xyz_to_rgb(XYZ)\n",
    "    sRGB = gamma_correct(linear_rgb)\n",
    "    return np.round((sRGB * 255).astype(float))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "0e17f484-fe3b-47cb-a519-c865ead9803f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "rgb(139, 177, 255)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAGFCAYAAAASI+9IAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAUOklEQVR4nO3cedzVdZ338Tc7AoIoiygKuGCaOmpq7jCUlqVpptlUd2pTltpoRZZaMy7NbZnLTN0PtTJT27RFcytyScOVFNcxEJdAwJBF9n277j+u+Mg1XMAFuaA+n3+d5Xt+53e+5+K8fss5tErSEABI0vqNXgEANhyiAEARBQCKKABQRAGAIgoAFFEAoIgCAKVtSwdeeJPfuAG8mZ1+RKu1jrGnAEARBQCKKABQRAGAIgoAFFEAoIgCAEUUACiiAEARBQCKKABQRAGAIgoAFFEAoIgCAEUUACiiAEARBQCKKABQRAGAIgoAFFEAoIgCAEUUACiiAEARBQCKKABQRAGAIgoAFFEAoIgCAEUUACiiAEARBQCKKABQRAGAIgoAFFEAoIgCAEUUACiiAEARBQCKKABQRAGAIgoAFFEAoIgCAEUUACiiAEARBQCKKABQRAGAIgoAFFEAoIgCAEUUACiiAEARBQCKKABQRAGAIgoAFFEAoIgCAEUUACiiAEARBQCKKABQRAGAIgoAFFEAoIgCAEUUACiiAEARBQCKKABQRAGAIgoAFFEAoIgCAEUUACiiAEARBQCKKABQRAGAIgoAFFEAoIgCAEUUACiiAEARBQCKKABQRAGAIgoAFFEAoIgCAEUUACiiAEARBQCKKABQRAGAIgoAFFEAoIgCAEUUACiiAEARBQCKKABQRAGAIgoAFFEAoIgCAEUUACiiAEARBQCKKABQRAGAIgoAFFEAoIgCAEUUACiiAEARBQCKKABQRAGAIgoAFFEAoIgCAEUUACiiAEARBQCKKABQRAGAIgoAFFEAoIgCAEUUACiiAEARBQCKKABQRAGAIgoAFFEAoIgCAEUUACiiAEARBQCKKABQRAGAIgoAFFEAoIgCAEUUACiiAEARBQCKKABQRAGAIgoAFFEAoIgCAEUUACiiAEARBQCKKABQRAGAIgoAFFEAoIgCAEUUACiiAEARBQCKKABQRAGAIgoAFFEAoIgCAEUUACiiAEARBQCKKABQRAGAIgoAFFEAoIgCAEUUACiiAEARBQCKKABQRAGAIgoAFFEAoIgCAEUUACiiAEARBQCKKABQRAGA0vaNXgFeO7dfe07uuO7cJMlHT70qe73n+JY97rpzc8e156RLt14584qxad+hU5LkmcfvyMi7rsmEZx/Ky5OeS0NDQ5Lk8/95d7bdZXCTZUyfPC6//8kZmTTuicyZ8VIWLZiT9ht1Sc8tdsiu+x2dAw4/LW3btW/ymGcevyPDf3tRJjz7UBYvmp/uPbfOLvt+JEOOOSsdO3Vt8euePnls7rjuvDzz+O2ZN3tqOnftmYG7HZKDP3Z2Nu3dv8XLac66zMHlXx+cvz41fI3LW/G+TJ88Lt86ccBan//CmxpatJ6v1RwsWbww993y3YwddW8mTxiVubOmpHWbtuneq3/j+3rYaenYaeOV1mPNr+tdQ47Lx067uq6v/DfbnObmeXVaOgczp03MBZ/fLkuXLMqxp12TPYd8qkXLf6sSBZqYOXVC/nT9BUmSwUd9tYKQJKNH/i6PDf95i5Yz6+WJeeK+Xza5beG8WZnw7EOZ8OxDeeGZETnujOvrvgeHfT+//cHJ9SGbJNMmPZe7b7ggox/5fU759v1NPmxW56UXnsrlXx+U+XOm122zp/8tI++6OqNH3pqTzr8nvbfasUWvoTnrMgct0WGjtb+mGtuxS4vGvZZzsGDezPz+J2escvuksU9k0tgn8tjwn+ffLnyoRe/Va2ld5mCTHn2z13v/NQ8OuyzDfnJGdtn3qHTYqGVz/VYkCm8Cy5YuSVq1Sps2r/3bdc9Nl2TJ4gVp265D9nrvp5vc13fbd+V9n/hmttp+7/zumq9m0tgnVruczl175rBPX5x37HFouvfqlyWL5mf4jRfn7uu/nSR56sEbMn/ujHTq0j3z50zPLVcNTUNDQ9q275gTz7szffrtkpuuOC0j77o6L73wP7njl+fm8BMuWuv6//rSz9YHwWEnXJS93vvpPHznj3PrVV/JvNnT8ptLP5tTvn3fes/PuszBSf/3T6vcNmXimFz0hR3T0NCQLpv0zk57HZYk2bR3/2b3Aob99Kzc9ZtvJUn2PvgzLVrH13oOum3WNwd+6IvZZd+jsvEmm+fZJ/+Yn1/0sSxeOC9TJj6dh++8Mgd+6IurPG5d9laTlu8VNWdd52C/Q0/Kg8Muy+wZk/LwnT/OAYefut7P/WbnnMIG5LrvHp/Tj2iV049oladG3JjrL/t8zv1U75x5dIfMmjYxSTJ5/Khccc77c9ZHO+fsT/bIby49MROff7Qed/nXBze77Ibly3Lnr/4z53+2f874SIdcfOquefKB65uMWbJ4YR6+66okyTv2/GA6dene5P53/fP/yXs/+o3ssPshadeu4xpfS6++O2TQEV9O7612TPsOndK5a48MOfrMur9Vq1Zp27bx8NHYUfdlyaL5SZJtdx6cATvun46dujb5YBn5x6ua7EU056Xxf8n4MSOSJN179c+gI4emU5fuGXTk0HTv1T9JMm70/Zk8YfQal7Mm6zIHzbn35v+q17H/B7+Qtu06rHbs4kXzM+K2HyRJWrdukwMOP22ty3+t56Dzxpvla5c/k0FHDs2mvQekXYeNstNeh+Vd/3xcjZky8en1WvarZX3mYPN+O6dP/12TJA8Mu+x1X+cNiT2FDdRvLv1s5s2e1uS26ZPH5tIzD8iCuTOSJEsWzc+fb78iYx79w1qXd/u1Z2fWyy/W9Zde+J/87DvH5ONDr81uBx6bpPHDeeG8WUkaP5xfLQ0NDZk3a2ruufm/6rZ3H3Ji2nfsXK9jdY9bYf6c6Xn5pefTo892q32e8c/8uS736bdLk/v69NslM6aMaxw3ZsQ/dAhpfc2b/XIeufsnSZJ2HTpl3/eftMbxI/94dW3t7rLvR1p0LuC1noM2bdulTdt2q9y+8nu4Sc+tm33s764+Pddf9rm0adsuvfrumD2HHJ99Dz0prVs3v2163vF9Mm/2tGzUpXv6v2P/DDn6zGw9cO+1ruP6zsG2Ow/OpHFPZuqLYzJt0vPp0WfbtT7XW5EobKCWL1uaE75+c7bbdUhmTpuYLpv0yg2Xn1RBGLj7+3LsqVdl6ZKF+el3PpqZ0yascXlLFi/Myd+6N33675r7b/1/+cPPv5GGhobcetXQ7Lrf0Wndpk3GPzOixm8x4J9eldfxsws/tsq5hYOOHJrDjr/wlefaZve6/PxTf8q40Q9k8347575bvtvkcXNnTV1jFObOnFyXO3bepMl9HTt3q8tzVhr3enpg2GVZsnhBkmTPIcenc9fNVjt2+fLlufeW/67rgz78lRY9xxsxBxOeHZnH7vlFksbzHnsOOa7ZcSs2cpYtXZyJz43MxOdGZtzT9+cTQ3/R7Pg5M15qfNysqfnLn2/M6JG35vizbsqOe35gjeuzvnOw5bZ71OXxY0a8baPg8NEG6sAjvpyd9j487Tt2Tq++O6R9h0555vHb6/4PHv+ddN20TzbtPSDv+8Q317q8fd53YgbsdEA6duqaIceclW6bbZkkmfXyi5k8YVSSZPb0STW+c9eer/IresU9N16c3/7gC3W9V98d6vzF0sULc+kZ++ff/6VbRt51dZPHrTjc1DJNDzWtvNfRqlWrdV7nf9TSJYvy4N8PS7Rq3ToHHfGlNY4f9fAtmfa3Z5MkA3Y6MFttv9d6POtrPwdjR92XK845JMuWLk6btu3zidN/WX9bSdK+Y+cc8vHzctolj+Sb187OOT+dlsNOuKie//F7rs240Q/U+C232SPHnnZ1zvzhX3P+rxfk9EtHZ+BuhyRp3FC6+covruMatnwOVv6bnz1jUt6uRGED1Xfbd61y28qHk7r37FeXN/37cdI12WSl8a1atUq3HlvV9bmzpqznWq7dJ0+/Lhf8dlnOvmZyPvy5S9P67yfLHxx2WV786+M17uiTf5gPHved9Nxyh7Rp2z7dNuubfQ89OV033WKl19D8YYkVumzSuy4vmDezyX0L589qdtzr5dHhP68t33fufcQa93iSxnCuMOjIlu0lJK/vHDx5/2/yw7MPzoK5M9KhY5ec8I1bVtmK79KtZw4+9t/Td9s90rHTxuncdbMMOnJottv1PTVm3NOvROGd7/5Q9hxyXOP5ivYd06vvO/LRU39c90/727OrHFb93zbkv4M3A4ePNlArfxV0hS5de9YWzKxpE7LR33eFp08eu9blzZz6Ql1uaGjIrJUON3Xp1itJ0nXTPnXbvNlTk7w6x91bt26dLpv0yn4fODkjbvtBJo17Mkky9cUx2XKb3RrHtGmTwUednsFHnV6PmzJxTEb84fIkyRYDdkuXbmvee+k3cJ+6/NILTzW5b+XrWw989z/0etbHvSudTznoyKFrHDvxuUcydtS9SZKeWwzMTnsf3uLneb3m4N6b/7vxG2PLl6dr9z759H/8LluudBhwheXLlqV1mzZrXFarlc4pND++6RZ9q1Zr3pZd3zlo/JtvtHH3zdf4HG9l9hTeRLbf/ZC6POynZ2XOzMmZPnlcbvvFf6z1sX++40cZN/qBLJw/J3f9+vw66dxtsy3Te6udkjT9R/K3lbbiV1i8aH7mzZ6WebOnZfmypXX7gvmz6vYVbr/2nIy86yeZNun5LF2yKPNmv5wRt/2wDlUlabK1/PQjw/Lck3dn/pzpWbRgbp594o+55ltH1u7+kGPOarIu53+2f33jaoXeW++UrXdo/ECYMWVcht94cRbMnZnhN15cJxf777h/k5OLzS1nTdZlDlYY89jt9WG09cB3Z8CO+6/xOYbf9MpewoEf+tJqD/W8EXPQ0NCQW348NDdf+aU0LF+e3lu/M1+4cESzQUiSm350Wm7+0Zcy/pmHsmTRgsyfMz3Db7okzz35xxqzzU4H1eWLT905Dwy7PNMnj8vSJYszZeLT+dX3Tqj7N++3czptvOmrPgdJ8uLzj9blfjvsk7crewpvIod87OyMeujmLJg7I6MeviXnHde4NbPyMdxWaf4fdtu27XPpGat+GB12/EW1ZTZgx8ZzDgvnz87zfxm+yne1/3TDd5r9tek15x9Zl1d8t/z5p/60xl+m7jH4k+m73SuHyMY8dtsqJ5ZXGHTkV/JP+x+z2mWt7JhTrqgfLd161Vdy61WvHHrptPFmOfqUK1q0nNVZlzlYYV0OBc2cOiFP3v/rJEnnrj1We9J2TV7LOZgx5YXcc9MldX3y+L/k/M/0azJmm50H1W80Fi2cm0fuuqbJSfOV7f/Bf8tW2+9Z16dMfDq//f7JzY5t16FTjvr85S1az/WZg+efalznnlsMXOvhvbcyUXgT2bT3gJz8rXtzy5VfzthR96Zdh43yzr2PyK77H5Mrz2s8ltu5a49mH3vIv5yb2TMm5c+3X5E5M15Kzy0H5uBjz86u+x9dY9p12Ch7vueE3HfLd/P0I7/Pgrkzs1GXTdZrXfcY9Mm0bdshkyeMyrw509KwfFk6b9wjW2yzW3Y/6BPZfdDHm4wfsOMBmfjcyEx9cUwWzJuZjp26ZeuBe2e/D3xhrd82Wdnm/XbOaRePzO3XnZtnH7s98+ZMS+eNe2Tg7u9b5b83WLZsaX2bq7lzOK+Gl8b/pb4gsGnvAdl5nw+vcfx9t36v9kD2PfTktOuw0To/54Y0B/sdenI6btQ1Y0fdm1kvT6z3dstt98g+h5zY5O8vSY4+5Yo8/eiwTBr7RObOnJylSxenW4++2W6XIRl81FfTc4vtW/S86zIHSeP7tOKw5r4faD5Kbxet8r9Pz6/GP/LrQl49Yx69LVvvsE+dT5g7a2p+9b0TMnrk75I0/qN69yEt++Vrc2ZMHZ8LT35HlixekMM/fclavyXzZvb8U8Pz/a8PTvuOnfPFSx5Nzy0HvtGr9LozB41u+P4peXDYZdm4++b52uXPvmX/m4uWHCK0p/Amc/3ln8vMqePTqWuPtGndNnNmTU7D8uVJku13Ozh7rsN/I9Cc7j23zuCPfC13XHtO7r7hguzz/s81e9L7rWBFSI/4zHffth+G5qDxa9kP33llkuQDn7rgLRuElrKn8Cbzh599I6NG3pqZU17IogVz0rFzt2zeb5fsftDHs/d7/3Wt3/QA3r5asqcgCgBvEy2Jgq+kAlBEAYAiCgAUUQCgiAIARRQAKKIAQBEFAIooAFBEAYAiCgAUUQCgiAIARRQAKKIAQBEFAIooAFBEAYAiCgAUUQCgiAIARRQAKKIAQBEFAIooAFBEAYAiCgAUUQCgiAIARRQAKKIAQBEFAIooAFBEAYAiCgAUUQCgiAIARRQAKKIAQBEFAIooAFBEAYAiCgAUUQCgiAIARRQAKKIAQBEFAIooAFBEAYAiCgAUUQCgiAIARRQAKKIAQBEFAIooAFBEAYAiCgAUUQCgiAIARRQAKKIAQBEFAIooAFBEAYAiCgAUUQCgiAIARRQAKKIAQBEFAIooAFBEAYAiCgAUUQCgiAIARRQAKKIAQBEFAIooAFBEAYAiCgAUUQCgiAIARRQAKKIAQBEFAIooAFBEAYAiCgAUUQCgiAIARRQAKKIAQBEFAIooAFBEAYAiCgAUUQCgiAIARRQAKKIAQBEFAIooAFBEAYAiCgAUUQCgiAIARRQAKKIAQBEFAIooAFBEAYAiCgAUUQCgiAIARRQAKKIAQBEFAIooAFBEAYAiCgAUUQCgiAIARRQAKKIAQBEFAIooAFBEAYAiCgAUUQCgiAIARRQAKKIAQBEFAIooAFBEAYAiCgAUUQCgiAIARRQAKKIAQBEFAIooAFBEAYAiCgAUUQCgiAIARRQAKKIAQBEFAIooAFBEAYAiCgAUUQCgiAIApVWShjd6JQDYMNhTAKCIAgBFFAAoogBAEQUAiigAUEQBgCIKABRRAKD8fxDXQ6yQr1kuAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Calculate the RGB value for a very high temperature (e.g., 1e10 K)\n",
    "color_largest_temp = temperature_to_rgb(mpf(\"1e10\"))\n",
    "\n",
    "print(f\"rgb({', '.join(map(str, map(round, color_largest_temp.tolist())))})\")\n",
    "plot_rgb_square(color_largest_temp)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}