{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The two-layer grey radiation model (aka the leaky greenhouse)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### I think this notebook needs updating to account for the reversal of the vertical pressure axis that happened awhile back in `climlab`. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Surface temperature is $T_s$\n",
    "- Atm. temperatures are $T_1, T_2$ where $T_1$ is closest to the surface.\n",
    "- absorptivity of atm layers is $\\epsilon_1, \\epsilon_2$\n",
    "- Surface emission is $\\sigma T_s^4$\n",
    "- Atm emission is $\\epsilon_1 \\sigma T_1^4, \\epsilon_2 \\sigma T_2^4$ (up and down)\n",
    "- Abs = emissivity for atm layers\n",
    "- Transmissivity for atm layers is $\\tau_1, \\tau_2$ where $\\tau_i = (1-\\epsilon_i)$\n",
    "\n",
    "### Emission\n",
    "$$ E_s = \\sigma T_s^4 $$\n",
    "\n",
    "$$ E_1 = \\epsilon_1 \\sigma T_1^4 $$\n",
    "\n",
    "$$ E_2 = \\epsilon_2 \\sigma T_2^4 $$\n",
    "\n",
    "### Incident radiation\n",
    "\n",
    "$$ F_s = \\tau_1 E_2 + E_1 $$\n",
    "\n",
    "$$ F_1 = E_s + E_2 $$\n",
    "\n",
    "$$ F_2 = \\tau_1 E_s + E_1 $$\n",
    "\n",
    "### Net radiation\n",
    "(absorptivity) * incident - emission\n",
    "\n",
    "$$ R_s = F_s - E_s $$\n",
    "\n",
    "$$ R_1 = \\epsilon_1 F_1 - 2 E_1 $$\n",
    "\n",
    "$$ R_2 = \\epsilon_2 F_2 - 2 E_2 $$\n",
    "\n",
    "### OLR\n",
    "\n",
    "$$ OLR = \\tau_2 F_2 + E_2 $$\n",
    "\n",
    "$$ = \\tau_1 \\tau_2 E_s + \\tau_2 E_1 + E_2 $$\n",
    "$$ = \\tau_1 \\tau_2 \\sigma T_s^4 + \\tau_2 \\epsilon_1 \\sigma T_1^4 + \\epsilon_2 \\sigma T_2^4 $$\n",
    "\n",
    "### Net radiation in terms of emissions\n",
    "\n",
    "$$ R_s = \\tau_1 E_2 + E_1 - E_s $$\n",
    "\n",
    "$$ R_1 = \\epsilon_1 (E_s + E_2) - 2 E_1 $$\n",
    "\n",
    "$$ R_2 = \\epsilon_2 (\\tau_1 E_s + E_1) - 2 E_2 $$\n",
    "\n",
    "### Net radiation in terms of temperatures\n",
    "\n",
    "$$ R_s = \\tau_1 \\epsilon_2 \\sigma T_2^4 + \\epsilon_1 \\sigma T_1^4 - \\sigma T_s^4 $$\n",
    "\n",
    "$$ R_1 = \\epsilon_1 (\\sigma T_s^4 + \\epsilon_2 \\sigma T_2^4) - 2 \\epsilon_1 \\sigma T_1^4 $$\n",
    "\n",
    "$$ R_2 = \\epsilon_2 (\\tau_1 \\sigma T_s^4 + \\epsilon_1 \\sigma T_1^4) - 2 \\epsilon_2 \\sigma T_2^4 $$\n",
    "\n",
    "### Net radiation in terms of temperatures and absorptivities\n",
    "\n",
    "$$ R_s = (1-\\epsilon_1) \\epsilon_2 \\sigma T_2^4 + \\epsilon_1 \\sigma T_1^4 - \\sigma T_s^4 $$\n",
    "\n",
    "$$ R_1 = \\epsilon_1 (\\sigma T_s^4 + \\epsilon_2 \\sigma T_2^4) - 2 \\epsilon_1 \\sigma T_1^4 $$\n",
    "\n",
    "$$ R_2 = \\epsilon_2 ((1-\\epsilon_1) \\sigma T_s^4 + \\epsilon_1 \\sigma T_1^4) - 2 \\epsilon_2 \\sigma T_2^4 $$\n",
    "\n",
    "## Solve for radiative equilibrium\n",
    "Need to add the solar energy source. We assume atm is transparent, solar is all absorbed at the surface.\n",
    "\n",
    "$$ R_1 = R_2 = 0$$\n",
    "\n",
    "$$ R_s = - (1-\\alpha) Q $$\n",
    "\n",
    "Introduce useful notation shorthand:\n",
    "\n",
    "$$ (1-\\alpha) Q = \\sigma T_e^4 $$\n",
    "\n",
    "This gives a 3x3 system which is **linear in $T^4$** (divide through by $\\sigma$)\n",
    "\n",
    "$$ - T_s^4 + \\epsilon_1 T_1^4 + (1-\\epsilon_1) \\epsilon_2  T_2^4 + T_e^4 = 0 $$\n",
    "\n",
    "$$ \\epsilon_1 T_s^4 - 2 \\epsilon_1 T_1^4 + \\epsilon_1 \\epsilon_2  T_2^4 = 0$$\n",
    "\n",
    "$$ \\epsilon_2 (1-\\epsilon_1) T_s^4 + \\epsilon_1 \\epsilon_2 T_1^4 - 2 \\epsilon_2 T_2^4 = 0$$\n",
    "\n",
    "Here we use the `sympy` module to solve the algebraic system symbolically."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAApcAAAA0BAMAAAAtVQ7LAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEImZRO/dIma7q80y\ndlRAyO8oAAAACXBIWXMAAA7EAAAOxAGVKw4bAAALg0lEQVR4Ae1ab4hcVxU/M/tnZmdnJ+OWWkgI\nGRJESglZjX9aUTr4oZ/ExKLVWLBb0QQt2sFibD5lSoRUimTFD9oqZv1glS1xV5CK3WIGAlIaMBui\nWJClW0orim3TxFqzVcdzz7nn/nnvvjvv7Y4bP/R+eO/c8++e89v37rx79gDkHHM59bZKbTjxTGxV\nuN46jde96XWfDCme15pDz2R6ZqDL6f8zMIcTT2nfkMA81hYER+/qCJl5v8hgVrqZGvkF1cHLDXSm\n4xmoF1eY2sVgPhhXS0pPJhlQvSZ/lVOtlFAYlTceW7y2cLZVf4HBPCaCjd4fuO2HAD/dqDUk49mg\no+OHFvqPPXEFfsdgVuYK+KnPQC39NBy9l11UrmW7KrdgbBZGuuUSg/mrbNVcktI87GvDjly6IaVk\nPCGdHLzdAHcB3AotBnMyDU62k/E5KC2lxJV1ZtUi22EN4PgMVNsnGczqTMpLktFYSXJo/gpza00Y\nOwjlsE7Q0Gcm4vGF+WfzULkK8M1qh8E8MdjSpvWHOYC9aYPbe8TbtpQWCecZgFNrMPWOFQZzvCOC\nzLtd1VPRYOJjXvsXlA56sgKTRDwFLF3Veheq/wG4uBMYzCOuMEybtCqPIJg3ppUOdIm3az4tEk4P\nYLEJo439+692kXdU+Nl3s6qvosGcfF2BCb/0hflnPT+e/IaeZh1AvY7tF/fffisKKrOeNDgxaZUn\nEczJdkoJX2A1YmCiGDcXHBXaCy4RzZfS3ksdpqZXP2j5ZlWPCxpM1BtHT/utfmHKjaewsTHAV0SN\n4028TDEOxAATts1PCUxaJxWYI64F22kU42CqzQXHt9bX8Pp9RepxR6fSInL0CHxNmHZVn+uAeQCt\nHrH6RSkvnqLGVp8fpOrpnyNromv5NmyTHwkFzNKKArO6ZC00lQtMtbmY8VVDQeVzcEMP6uj162ul\nD0P5tvewTFb1uQ6YX0a9l62fopQXT1Fjq4+/BGaMOzSFDefxVaT89ny3yWqSVhkUmI1ZYyxELjC9\n3/pPApSePIzjvl7tM6sPQf1D6PUn336uCe+GJ9BteXn5ybuXl+eRdLhwy/LyF5eXn0KufkMuELmh\nixfPhjyQ0aIGSU3G2gA7VFaHn+awzz2Om6rKr9yp4dPipvU3ArN+Wdl5Iwhm8uivNxc2RDBlbJtX\nlPoTrSviCFzoqLvZXHyufTJp242CmYyA3JqLF4/hhom5MFtxeedluQJTBoWt3hzKb7xTf5Ml+sks\ntQqB+VpT/PJd/0rxxHnNt7UUS4F5RRFfgONtdTdg+lwDZmMGvgFwB6lmXJIR+GpePL4oOYsUQ/TO\nyxbua05hE5gtFI63K39nHQ1mY2Fh8bMrMDrLXOcaejJTR393c3F/gMZaACsE5p34e4heNQayufhc\nA+Y5gD/Gf4BSETgRI+nF44uSs0gxxNt5J7rWksPGJ5PzA/mDSFoowOc95w+QHP2N+9Puk7rHsGFq\nCW7oEJh/gpL6NPoKy2RVnytgVj6/W22z77V+UlQqAl/Di8cXJWeRYsgEfelpg/KMteSwEUzOD8ZX\nWCZp4fuPYI60rIWmQk+mPvqLbvkH/afQ3fQq+3Rfz0vPdfk1n/rIw6guEcmqPlfAnOz3+7MAT4v/\nwD0Rga+RiMcXJmZSnEmw1fTms2+dwd/Svaskc19aDhvBBMoP5PGRtKB8er0Lkz0ydC8hMPXR31XD\nVX+BBQo1UsfJxizx8SLnW7MqCYQrYLJ29MgRjEBW4buJx2cnZlKcSbDtdEdzO0+OWB5TCkwaphjh\npxU4Cbpg6qKUHP3FGd9rZ9Szh8P45ik/mUSPtkZ7RJQ6IsO74cJNDjdc6IhG4FrjZ4vE47P1TMp0\nujgT1CHmgnq31DB/cp4638EfhY8xz0srVOhwwdRFKTn6i1e+22+H7/kCaBzUjPOv/qWZkOE0zA2X\n4KIR+K5tPD6fZ2UuGz6kizMhFeaZD70J9wFQspe0Uf03r35Hk+4tUIKTQ/mFFirqopQ++ruWSG/r\ngUbqmC+pfGB9hTn39Pu+iGZhLiT/JKQbjcD3bePx+TzTZbqP6+JMSIV5WL5iIlkcfuB+3FHVGOv3\nLzPlXVPgo1T9LuE41Qao66KUPvoT314m1vSeCZWu5W6YCv7bIh6Bv5aNx+fz7BkuGzZ1cSakwrwf\ngd4zYdP/tsBC0kFyehivUpQCPvonAiidWW0mWEOfxiPwl4vH0+MynRRnfFN3tn31ojvdHN34p7Iv\nv8lexmb5TtfxeWeyVeTwInAPi7BFuSx2EKd3tRks75x249JWIeisM7QIvMMibFEuNfw92/Fjnc6p\nNSev60IOLQLvsLhlqby0Au/UdUj17wn8lvGW3obnlf/VuOytRBMVgamLslgqY7FI0o5CZbqYh83m\nyLn8dgXg/T2ORu0zqpB3/YaKgOuim4zB23w36Su3+Yj6bB35NOnzPvNybtvhK1IEpi66Gf/e5rsZ\nR0Vs9/WU9n1kwvvM9QSTIjB1UQpqg5ehbb4F1h/9Nyk/31Y3LkoJmFvxC7hriZY3F10W03XRAhGU\n3jI+mKAyXXl1jWcFPCX8FJlO8MZJx0kuSpn/dG3Ft9n4vV6wOgKpixaJ4HDT9aTLdO+rzzG3iCfX\nTzFat3LsmrFm8mRaTpSq7H+0l6EQEaHFOWX155Cp1EVDsjSP+sCg0UlL4NEXPIgDGj6LXfk8PYvm\ngo1carhVI21WEMztUNb/I9H29hYRodKXlGJLXZLjRJIRm2PJYF8bv6dCOs/2eyF2Fk+7CoqjuWAj\nlxqbBxM/+D8VXB4gJHpFdKcOIVVakalztxVQh5lJ1qgPDHYGFBqdWsh/QJNZ2lVQHs1FNXLhCIAp\nhbygzzQTvwwXe2m24oREBswHf40a1ZBhRgU0pIq8MeoDg25APN4b6QXYmSztKiiP5UKNXGiVBtMU\n8oI+08yzzUwwQyIDZkuBGRz3BOuiQVVk6j6wkLj68O9D7ExexBXEcqFGLvSaBjO1VEYbk6N3f9P2\nOjlsRaZEAma1kwlmwoWdmkgsiyjVB5ZvGA+GSNihK79by5Fn5kKNXKg4GExuY1IHdmpjMp03dpXR\nK2CVLFtRIrLHbQFzJxQGkxaxnuxKB1qWjlImzJgrSrO+lHKUmQs3cqH+YDCpjUkd2KmNSXfeeCvh\nR5xR8gT4r0wtssdtAbNVHExaxHqyS6k+sFzDhBlxRWlSA1XCZWYu3MiF2oPBpO4r9SVPbUy680Za\nt2g97LQ0SsiQKo+SiUgft233VqNXHExaRHtyFwl0SlJc6YsJM+KK0qTOCtuhRp4yc+FGLtQZDOY6\neUIwuY3JdN6YUKdmdOOWOTulRfa4rZ/MaSgOJkViPZlVLhlqEGFyibiiNDWYrj+Tpo2AcylxIxfq\n7ponA30jOnG5QnMFZoso6bwxai9CpWOUDJcIK5I2JOlE+OvCwhtnfOWBM17EeBL9BvWBySx6N2FG\nXHGajdmko8xcdCMX6g8G807qyUIwx1qqTct03shieJSa6hgl4dLdEclxW8BE+SFPN8eEFzGexOIc\n9YHJLHo3YUZccZopMKO5qEYuHPpQHnkyk21MexLhHt2953EwSp7QEZnjNr8aSu2ap5tjwosYT9pC\n94HlsAcbZsTVlOlG81xGc9H/MNc11LH0l4C4SrQxpbpjzvb7/wCrJGbq7ohOCN+AeXP/Z8LLeedF\njCdtpfvAcvkwYcZcmW40z2UsF2rkQm2qvUkl07P2J/ia8zCdN8Jw7kbJ4TFpj9sGzJROPob1lE8/\noCVhxlylXnPrx5olc+FCO9SvWuUgJQf2jM4bthGltAd73L4pLSzEsZ4KmbnKEmbMVeOga+HR1iyR\nS0k2rufxpyUyzIE9o/OGTI1S2lGx43ba3nI278mEGXFlG6jsykJlmh39hFaprjdF+e37xhBo3G0g\n3D6/MRdvWwkC57uK+i8FwG/fuFDP+wAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$\\left \\{ T_{1}^{4} : \\frac{T_{e}^{4} \\left(- e_{1} e_{2} + e_{2} + 2\\right)}{e_{1} e_{2} - 2 e_{1} - 2 e_{2} + 4}, \\quad T_{2}^{4} : - \\frac{T_{e}^{4}}{e_{2} - 2}, \\quad T_{s}^{4} : \\frac{T_{e}^{4} \\left(- e_{1} e_{2} + 4\\right)}{e_{1} e_{2} - 2 e_{1} - 2 e_{2} + 4}\\right \\}$$"
      ],
      "text/plain": [
       "⎧        4                             4                4                ⎫\n",
       "⎪  4   Tₑ ⋅(-e₁⋅e₂ + e₂ + 2)     4  -Tₑ        4      Tₑ ⋅(-e₁⋅e₂ + 4)   ⎪\n",
       "⎨T₁ : ───────────────────────, T₂ : ──────, T_s : ───────────────────────⎬\n",
       "⎪     e₁⋅e₂ - 2⋅e₁ - 2⋅e₂ + 4       e₂ - 2        e₁⋅e₂ - 2⋅e₁ - 2⋅e₂ + 4⎪\n",
       "⎩                                                                        ⎭"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import sympy\n",
    "sympy.init_printing()\n",
    "T_s, T_1, T_2, T_e, e_1, e_2 = sympy.symbols('T_s, T_1, T_2, T_e, e_1, e_2', positive=True )\n",
    "system = [-T_s**4 + e_1*T_1**4 + e_2*(1-e_1)*T_2**4 + T_e**4,\n",
    "              e_1*T_s**4 - 2*e_1*T_1**4 + e_1*e_2*T_2**4,\n",
    "              e_2*(1-e_1)*T_s**4 + e_1*e_2*T_1**4 - 2*e_2*T_2**4]\n",
    "out1 = sympy.solve( system, [T_s**4, T_1**4, T_2**4])\n",
    "out1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAw8AAAA0BAMAAAA+pxkWAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEImZRO/dIma7q80y\ndlRAyO8oAAAACXBIWXMAAA7EAAAOxAGVKw4bAAALwklEQVR4Ae1ab6hcRxU/u2/f7vu7b03BlgTJ\nJWmVojGvVK0tkSwt+KmQ16DVGDDrn6aYWrpYjAaFPIl/AiXkFT9ZpV0RW3klvJcPFsxTsxgo2nzI\nliAG5DWPSiyKaZI2bZrX1PXMzDkzd2bn3r139+7zS+fDnZlzzpw/89s7s3fmAKQrW9sDLemc0dLJ\nvLqs5VM18gONuN2jV7VUMayV8EC9Kq5VFADrppPaKs0llVxLucF6dWDAoRyosoHCg3VudqvLrW4S\n/w/+YL26OOCQRq5XyMLRILGpoWZi0QSCj9/9iwRS3UWy9cq1t80l9NCfmI4ZtP8rilm6HiPksDJd\nL3MN2Fp1DPTUjfLqAv/UetLKg2rc6KMuzsUMLq0q5tiVGCGHtcHpd+v6X2uijlVgeKabhiR8v1e5\nXy1kAURuMYkLXhmzZP41DgjY3pTDp1JYetVrMJoYC8RwDcZuRI9Nzony6nwWQJSD5I44khqI0hOx\nQOyYlQM3NpzxMd1bY3g+ViwQ41cyAiLKq0yAyFd9gSWiaSDy47FAHJyW6tIAcTzCgdyWs3XFWrd8\nV0iGgTB8wWQqQDHFshhS6zSjvMoEiPGKY83T1UHbcWogfhIPBCGQBohFjxeC9Nl6KYAJ5BYehm9B\n/u47SI6nXPI3/5RiYirAjoAE+6qivMoGiO6uyaDhNAracTIQuVbWQBSm/V6Vvgw3NSfuqQF8eyW3\nDT4Cz5McTbnk5+tjNGMGiH1+femoUV5BJkDc190ZGfSpZwGcOBmIPGQNxOSK7VXuhb1YHmmOfXH5\nBwDlGsCvj5yrwMNwpo6CH11aemhp6QS2JL9Yn7hmUQGGIoBFsRTF9UoPzQSIV7S6jsYGEf3e36ug\n4Z92nPmlpRd2Ly01cNR/MgditNLhjCRMCXMSiFXR+DocrIpa7waSX6yWrlpUgLOq3+czyqts3ogn\nu3sngxZAOHHSG5ELMgeiGOHUVCAY4o14QzQALhFitAgpPpRpa+alqTwN31HyfT2jvMoGiFp331TQ\nAohACnOcBER5fn7hS60YNbRLUxUjKFnid/dShNBwANCSQOwEmEShR0mQplzxoUjeMBCnAP4WoTEp\nOc6r/oHAt2Gi0d0VFTQC4cTJewQy5uK0pAPiRwjEuQh1k4twU10C8XfI4d/X/DQJ0pQrPmy2qaWv\nbhL7e18lzqu+gRh5AH9Ws939U0EjEE6cBoip7IAovIhAPBfl1Nlz6G+5hq585jDKHGI5/u1L/ogD\nzzjeuOCQfkqsV/fsOTHbj3K4H0/jxurdVaigEQiw49RA5J9ZjfMk/EYc3DXffup5WuGV6dJbTy1c\nnz8ZqN4tRQSiptoRTwGELIWg0FQtBkL27oX7PVRF8j+z8MqvWVA92i1SbhaBKEaPdzgCCFl0nBoI\n5kTUYSA2ATwI8OmwZD6A4RoMzSpaC4HocgFTnqHhp1//F6Imys2qks+JP75O/0DC1JBAZzMLrzq1\nMsWj3SJNAgJxgKW71hdIwsSZS/A2iUFhIBpQehPgh2FrY/ibmYaRqqRNbLrtYzDSCvPddulTq8Tf\n0267TOwPt1Pfz2bglccRJnm0W6Q/bdrZhCMs3a1+/LFjSiR9nATEmQD/G8zCyH8BXg5b+zPA0RWY\npN82DFdgsBcwYduqPVivPNpdEm7WWVwLdUZmU2grP1pFIHBXugKALVOaAAsVKBAht/WuFOul0dJP\na7BeebQ7pHVv1rtsi/1EZ8aOz8j2XkXB/UCXfdTCbSNc/BcwYYnM24P1yqM9TMJo+rgWSj4X5XeE\nbP6aGiEPxcUJIlrfrUhi2wiXV8OdtWkLr+gwdwBe0UWACIW1h0iC3Me1kBiesCzUUfC2qpLG/QDE\nCWKoiG2DyvdFfSv31q4WXpnDXGk3Q6+EdqcYUqmJrHwVHwMvY58D2PBLMoP7AYiDq1AJ3Wb/QZCP\nh3hr1BRe0WEuW/R5xeegLJOwljHbsoY03kJOkmshW0FPvQst+OAdNFLuBzYQZr0srwqpRfEYZAKi\n0G8X4ZU5zJU8r1c8Ksa7yyxjamcPFAxDutTA7rigTcVo7ZelvHoRQf9kU9iSnxHuGyHXy3XLKHTL\nQyhSmJaSa/mgXepSJWQ0O6+U9i2vhJRLUm7LMm4aP55B+n0h3uCaQ7g0wdAXpAG18tpvhFgvc8dz\nDYCV7TgVeeyucaH94NGw2ey8ktqHqn8JaZekDZX1uE0P4f95+HiIN7jm1qbQ/Yg0MCrM0h7B/yCe\nwdkfO3ZYMET0o+HfpSBmXzYu2jqVV+owt3+vvNpHn/6AsEnapcF5cYQpv3DhScEbdCm8Jy2cr2J1\n+8l3xRc6vRHyOyL/8/aJFgwLLh52BGnOv+SQXh5FPN0JFfLqkCL17ZVX+8TCDalealcGxUKBi7U4\nAq2J1qDL6GVpQRxxcLGXJkGdaoJ4EUQW2kssNcD6NY9uc5irmT165dP+GtAZkVYONyAnOu8muxYy\n43ptUYrfxmmjgE8QDWV0RewRMI6gRV0LGWHRKt35s6ZNMb3YrONTQi4wwrplDnM1Ka1XZNin/Zud\nx0lPA+4RACeTXQuhZExccdMBlJtMh35UCdP6BFF0VMkdWxZvRB7f3+eYFluvh/zVCAGEdGs1ggfw\nDeTkWh72ns7D3JRekWGv9sNHOlxav/yycAO30CTXQuh1TFwx04GL/ZwMtxMISfY+xHFIzctxifiB\n+HmXRn1f1rG+PJnchVIjESMjyEm9IsMptV+aTXjM6YuLPfZNx0VmUm5yGiByq92uhVg5fhEtNLlj\n1/gp1pF1rIH47u9s4SS9pF75DHfXj/OT7FooTr1vOhgIzk1OAwR+ccZfC+m4TlYigfBlHWsggh6A\nSOqVz7B2OLKB83skkhlmxKn3TQcDwbnJqYBYqPC1kEm01bm3Ya+w/VgFoliYdWzGi2EMxEi9FyCM\nV+yCrZ2puBzL7yTTT9DCQ61tJKbDiVEfxeqcDgaCc5NTAXG+WiSnohKOdWiFN1QuskzL1VTV2BFE\nJOp+CHoBwnjFduw0YKb2lO48eZW3xagsY0u9nhdDlS2ejo7UbJ2bnAqIMy26FopMONb2iw2Vi+wc\nqkv+vqhE3aAnILRXbNxJA2YygPxeM90krYn3+FooKsvYaKG4ZCK2ocoWT4c5zac3QucmpwJiao6u\nhbwJx5yNLC3fSbnI9J0ePqLGrONwQrJJ1C03ewJCe8XBh7U7hlkkeX2dr4VkarUIJ0a9mRf8P0vJ\n2dZ0dKRm69zkVECMz3xCRTDVEHW5BrAqGs4ZNVImp4nFByZCjMpZ/EZvYNskJNMesQ56AkJ7xQYc\n7UzuLd35t0NVpUBGKsKJUW/mRRuVDT0dZqbUG5HTucmpgBi6tqgMTAWiFkBQOpp1Ro2sf0Cprljo\nuV3KmHWsxuuEZALi3/Pzbx2zhZP0tFcs7GhnsjDM7eT19u+RrA4nRr2ZF9uAmQ4nNdvkJqcCYuQd\nAmI4iEg4Jvv42T5Z34kvhj5CNI6dwqxjNV4nJPO/JoBdRi5xS3vFIxztTBaGuZ28PvoEyepwYtSb\nebEMhKbDTc3Wucl0yER4WMM7O6W3ccURZTIq4VixYf+mzc+CTssloqpk1rEarxOSDRDXLdlkHe0V\nizvaidxjuvPBr9F4HU6MesUSK4VVQtPhpmbjSjcnZeVVF8JCP3VrfGenvUK0qIRjYp9st9+mXGR3\naVJZxxGJure3f9NptCtFe8WStnai9pjuPMXJE1FZxpZ6PS/siqxD03GIGRepwbnJZwJJUJcvLBRZ\ni1siUzTynjNqLeUCoRkAqRN1Q2PDTdsr5mjtTOitdmZGhxOjXs9Lp0UzUwwEy+DpoigTTvaSJHY+\nPmyRyjPU9ZxRa8HOQ3VmpU/U5ZFObXtFTKPdkU7ZzdesARxOnHo9L9ZI2TEzdbPNzPGifL5lM/y9\ne8PkbgnHUtZzqM460ifq8kintrxiXlbaCwFrFLUOJ0a9mZfwSNX2nOYrxv4HSHhk1Vp1OlW8Txnk\nDJR36+lf3xikofd1x8/A6VnB/x9BpXmbbZ/C4QAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$\\left \\{ T_{1} : T_{e} \\sqrt[4]{\\frac{- e_{1} e_{2} + e_{2} + 2}{e_{1} e_{2} - 2 e_{1} - 2 e_{2} + 4}}, \\quad T_{2} : T_{e} \\sqrt[4]{- \\frac{1}{e_{2} - 2}}, \\quad T_{s} : T_{e} \\sqrt[4]{\\frac{- e_{1} e_{2} + 4}{e_{1} e_{2} - 2 e_{1} - 2 e_{2} + 4}}\\right \\}$$"
      ],
      "text/plain": [
       "⎧           _________________________             ________              ______\n",
       "⎪          ╱     -e₁⋅e₂ + e₂ + 2                 ╱  -1                 ╱      \n",
       "⎨T₁: Tₑ⋅4 ╱  ─────────────────────── , T₂: Tₑ⋅4 ╱  ────── , T_s: Tₑ⋅4 ╱  ─────\n",
       "⎪       ╲╱   e₁⋅e₂ - 2⋅e₁ - 2⋅e₂ + 4          ╲╱   e₂ - 2           ╲╱   e₁⋅e₂\n",
       "⎩                                                                             \n",
       "\n",
       "___________________⎫\n",
       "  -e₁⋅e₂ + 4       ⎪\n",
       "────────────────── ⎬\n",
       " - 2⋅e₁ - 2⋅e₂ + 4 ⎪\n",
       "                   ⎭"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "quarter = sympy.Rational(1,4)\n",
    "out2 = {}\n",
    "for var4, formula in out1.items():\n",
    "    var = (var4)**quarter\n",
    "    out2[var] = sympy.simplify(formula**quarter)\n",
    "out2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhUAAAA/BAMAAAClYic0AAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAIqt2Zs0QmTK73URU\n74mR/c/RAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAH2ElEQVR4Ae2bTYgcRRSA3+zs/Ozu7Ox6EPWg\nO0kUT2FHEuMlsAMbEfSwS8TgwbBtcAUFs2sUSeIhQ1CCB00IGDwEMlFQUTRziyczwZOgZi8GjJAd\nyckfspuY/7+1frtfVdd0V/W0l0wK7Hrv1auvXr3pru2uigCkZNbQa++W/EvB3Ke2BHJPShMVOe3s\nVSn1al28IWe+/4CU4uvySlxZjocYPZbiwCvGbvHGWC6J+FRFcI6NxgOlx6AUUq8bqRMF0AY8cIs7\nD110iOIhB18n14zDzZk+uHidMwuXHNg/OPg6uZaqTu72znbgRY8R5f1hhd9p5ZXAqVBP0Mmmix14\nhg+fu2mDFD4NB18n1/9tIYoEP7DrEx6lzMU1+6Dz05G+P3uRzVGNs6bG0mmT1c1mBAtEvgJLTSYn\nyEWpHRFI/vFxL6I5uukXQ/O2HVcMVkeTCSwRfR4MLzAlQS7KPIuSpdcTnm6x1v82eWZTyIURLAYb\nbkAffyoS5CLnmUL2bV3kwvgZkEYujGARce5KN7nwp02E7evexCqRu8iFcSFKIxcYHI4YBvmtl+C+\nOI4mP3QUHkEqFZPnYmheQzE1hVxgsCFiONFmIyXIxUYU8cOt/E6kUjF5LrItDcXUFHKBwYaI4VM+\nboJcHCI9Z7+i5UVYO7neQyplJs/FAEWFSgq5oOCIiAvifkyQiwYK1/AZkzwX5jeiFHKBwYaInxIT\ncswFuSeKFZSLO0gWYrJc0JttbxhGLF3mQgeHIy7Nw31sZLdc9J8msdVQyGeJjlQqJsoFBcN6jcTV\n7nIRAocj3gfwBBvKLRfffQ7QN4pCfhLyqfxNpWBYhcCB2F0uQuBQxJlvz+9osOGccpGvkZDxwwfZ\ntzYFUTNpx+LmmmaKVxkYGibHwtrLj5rsdrYwOBRxjux58ZFxLjK3D45fXz3TxsOM/bZ65eBj4iHL\nAslFkp0chcLxiomBk+zkKJRuwSIXg/SNvNwG8nJeqHEmv54HOAPwNlc+On+2DpO42VJWKJ3A/VVL\nGnJLFSxyMbxMBugDGJuH/iYaCyqQuQrwnDSRJU5/t5JNUbVGoa6aiYDtNlzUYTRKd+DxFoPPNUj1\nHsBUC7IeEWUp1qD/FsAGoW+/Omp+qqW/udYo1EkzUfCguXOUVaN0CZ6aZ2NNtUlVByCbD0PMIC5F\ncrdcAUC3SsxODu7sy2EKGEyzvr+1YKAYTLbguQU28DdieLI2+EW8pQ83fAsVSm1FtVUwxQyGqA2X\nzsOkCC6TJ4A8quwKbG2A4jQbOf8rD4AsIaJk6kSI2cmRvnpNKeVd71BzGAwPUnvUhgttN5cgvO7B\nWz0yxkSLD0TXhqJ48+AWtoQIMVclQs6TDU41WYjgD3gM96EmUZ6n9RapOdWIIvshkxv4dfLKUPhd\nYOjaAKWGhLI62L88VSGGHDWORJ/JLVMftVDKUdg/iqwBuMQ+mKZJm/ya7DBABzCCUtEE7gCU8xDg\n3XUYbgsae/i0XPhLSP7ZBeJ2XLg6VpTyNYw1UTcfDK/Q5QpvuCCvODGg+J6ByRX8YRXyL3gcxB4+\nNRf09QLy59aR+6VA75p/uKfjlVHI8a0X9AvA0NpD7OVW0GYvcfC5jahHcnCWnh7SB4UW9qSpuaBL\nCMx628iViYeon3NhXQE+Qx0RmI074KFGa5FRCs0PUIfk4IkawWTEseoRjygiF2K5H6A3w+r1NXLN\n0M+SBvkvvsxNqz6MAuV5ag2DYaytffOpvZFmBA8cfpW6dAUm/Yt8n2cPfZDLH69srvq54LuAf85c\nXgMg1tbL2k4ODcBcBsWNJlo5BZ7hKnu/UMAjJHV7zSTNagQXx+nnFBgitgeT7jn6swPgf4qiPiOs\n+RrkaT2j7eSwNnTxTyYBTiKzFIfaQ3Upi1qAc8uddnICfwE3gU8C+b3UYglGEQOcqDDGCAlGltKC\nlPz6MND1Apbq6k6O386F4GQSoK21UfX9n370NLMAl8nvukpr0lQJb2t2qn4Z/mC0A0soZ7Llgtwe\nF/whMrsvVn1FCNvWbaDSqVrkB1SffzIJ+RCCdF9cCf0TIwEukQW8QUfoXATcCN402dQ72oFRxARw\npMUoOf7A6URNn6tE7uSQtxNxMgn9Ws8YNX8R4nZyBDxdMIqYRDjDM2qXC9J1MmJWZO2RuYjwMjad\ngbidnITwaLAKFXs5drkgr+g70VQ6n0wiJztx3CvUVc+U4Co4GuqUi+xN/FRHnEyq07LQJpqDqlda\ncAUcA3XKRfFfvJMTcTKpzstC21+dVb3SgivgGKhTLuA63skxnKXKk0l1XhbayAFtJyctuAKOgbrl\nYmsB/fGKOJm0mL3qklt4VzWkBVfAMVC3XOy5H0VMv07U4p9MqmYLrXBhWvVKC66AY6BuuZh6GkV8\nNnSWuk+eTCIvO7H/kpaLtOAKOAbqlouxL9DUIk4mkZedmLkxrzqmBVfAMVC3XIzQzQxZIk4mpYt9\nvdJSfVODY3AM1C0XfBtCjTodje2cpINSKQ5gt1yUG+pA6Wl/pYdSSQ5gt1wMtdWB0tO+Tw+lkhzA\nbrlQh7nbtHu5CH7RY3Um9+E/EEFrb0mLo2y+9/43RLKNfpv/9LLurRtBnW32jtCXqmpDD2pvnBaT\nttvYuqszdKQpppcXG5939WwjJ/caPUzlhZ2rSqUX6zP1YNYvB2JPSvwJ+Q/EcC3Z3XqUHQAAAABJ\nRU5ErkJggg==\n",
      "text/latex": [
       "$$\\left \\{ T_{1} : T_{e} \\sqrt[4]{\\frac{- e - 1}{e - 2}}, \\quad T_{2} : T_{e} \\sqrt[4]{- \\frac{1}{e - 2}}, \\quad T_{s} : T_{e} \\sqrt[4]{\\frac{- e - 2}{e - 2}}\\right \\}$$"
      ],
      "text/plain": [
       "⎧           ________             _______              ________⎫\n",
       "⎪          ╱ -e - 1             ╱  -1                ╱ -e - 2 ⎪\n",
       "⎨T₁: Tₑ⋅4 ╱  ────── , T₂: Tₑ⋅4 ╱  ───── , T_s: Tₑ⋅4 ╱  ────── ⎬\n",
       "⎪       ╲╱   e - 2           ╲╱   e - 2           ╲╱   e - 2  ⎪\n",
       "⎩                                                             ⎭"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#  The special case of equal absorptivities\n",
    "e = sympy.symbols('e')\n",
    "out3 = {}\n",
    "for var4, formula in out1.items():\n",
    "    var = (var4)**quarter\n",
    "    simple_formula = sympy.cancel(formula.subs([(e_2, e),(e_1, e)]))\n",
    "    out3[var] = sympy.simplify( simple_formula**quarter )\n",
    "out3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The solution is\n",
    "\n",
    "\\begin{align} \n",
    "T_s^4 &= T_e^4 \\frac{4 - \\epsilon_1 \\epsilon_2}{4 + \\epsilon_1 \\epsilon_2 - 2 \\epsilon_1 - 2 \\epsilon_2}  \\\\\n",
    "T_1^4 &= T_e^4 \\frac{2 -\\epsilon_1 \\epsilon_2 + \\epsilon_2}{4 + \\epsilon_1 \\epsilon_2 - 2 \\epsilon_1 - 2 \\epsilon_2} \\\\\n",
    "T_2^4 &= T_e^4 \\frac{ 1}{2 - \\epsilon_2}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the special case $\\epsilon_1 = \\epsilon_2$ this reduces to\n",
    "\n",
    "\\begin{align} \n",
    "T_s^4 &= T_e^4 \\frac{2+\\epsilon}{2-\\epsilon}  \\\\\n",
    "T_1^4 &= T_e^4 \\frac{1+\\epsilon}{2-\\epsilon} \\\\\n",
    "T_2^4 &= T_e^4 \\frac{ 1}{2 - \\epsilon}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAJ8AAAAOBAMAAADDD9M1AAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAIpm7MhCriUTv3c12\nVGZoascqAAAACXBIWXMAAA7EAAAOxAGVKw4bAAACl0lEQVQ4EVWUS2gTURiFv0kySSavjgURdGGI\njwoirdSKCELowq2zUfABDQXBZ1sqtPigCS51UQsWqi6MG7tw0WyqG6EPfCwUOrhwaWZlKRS11mhr\nqeN/70RM7uKcM3P+Ofefe+8MRq4zj/l+sIRWELDZtpiHbcPH4NzbfXUw24Zskq8jB6B116G6/dLp\ndhsDThFbIwRX0UoCNSdd8wlGlqpHlqgbwGnPGMPy/QHMAsVSYD/2f9EYcB4e8Al2oBUEvAQ3SW3S\nUk57GBUNXIRekkNzEHWw7mmbF1dsGgMeQn/+GpxBKwnUfBiqTnyCnkqsA7OsgQ0oEpUSrDLRH9qm\noK4bAmYcCRztYgStxNT8E0Zduah68VUnFoC5Cgt5HZisSaC2g8CmAKacqL/oiS1KofDJ7xI4C+ak\nvMCfkTo8kg5L0S/zeSkhUQvsvuXdsgSNAXFZ1OKGilJKjfgvQ9SJAYwbeyDmT9ehKO17YcdYU0Wj\ndbuLnaXmgESB9IWpCakRpUeiYEiHEgjXS3y8s16H5Gy66MrNo6rqqQKxIVxuDsjBZ1Lr0qIoPXKY\n9VfGGs8UWAgAll8V81LRbstByapSa1wgtNIUkBJnTPpxUUoNxbIpVdewCW9aNpEAlNfvXFYLCfPy\nBSg7WlEHuTHgA8bWb7KsHqJkajR3wlenpSZP9MidPg3Ke47sTLtDOssWbVsVQjWzISCSJWXLBGFb\nK3lG8xLsx+ogsSIdclADS3ZkRb4ZnkE3zGk7JEtfUR3+C9g+PHiJuw7LaNVTCTjhmvfJeFRn49Nk\nAuCdc9ajlcwYxuTw3rK2jazM0xDAjO//JtIrPwetZKM1m7fe5OF27ggcb5Ofg4ZYTpTR1uuQ9H2/\nHNitQ3K0/gf8BR6JDv4GdzlSAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$282.203889523582$$"
      ],
      "text/plain": [
       "282.203889523582"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "out2[T_s].subs([(T_e, 255), (e_1, 0.4), (e_2, 0.4)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "246.627893584128\n",
      "282.203889523582\n",
      "226.730624779963\n"
     ]
    }
   ],
   "source": [
    "for var, formula in out2.items():\n",
    "    print(formula.subs([(T_e, 255), (e_1, 0.4), (e_2, 0.4)]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Coding up the analytical solutions for radiative equilibrium\n",
    "#  These use the analytical results returned by sympy and wrap them in callable functions\n",
    "\n",
    "def Ts(Te, e1, e2):\n",
    "    #return Te*((4-e1*e2)/(4+e1*e2-2*(e1+e2)))**0.25\n",
    "    return out2[T_s].subs([(T_e, Te), (e_1, e1), (e_2, e2)])\n",
    "def T1(Te, e1, e2):\n",
    "    #return Te*((2+e2-e1*e2)/(4+e1*e2-2*(e1+e2)))**0.25\n",
    "    return out2[T_1].subs([(T_e, Te), (e_1, e1), (e_2, e2)])\n",
    "def T2(Te, e1, e2):\n",
    "    #return Te*(1/(2-e2))**0.25\n",
    "    return out2[T_2].subs([(T_e, Te), (e_1, e1), (e_2, e2)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "from climlab import constants as const\n",
    "import climlab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "mycolumn = climlab.GreyRadiationModel( num_lev=2 )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "climlab Process of type <class 'climlab.model.column.GreyRadiationModel'>. \n",
      "State variables and domain shapes: \n",
      "  Tatm: (2,) \n",
      "  Ts: (1,) \n",
      "The subprocess tree: \n",
      "top: <class 'climlab.model.column.GreyRadiationModel'>\n",
      "   LW: <class 'climlab.radiation.greygas.GreyGas'>\n",
      "   SW: <class 'climlab.radiation.greygas.GreyGasSW'>\n",
      "   insolation: <class 'climlab.radiation.insolation.FixedInsolation'>\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(mycolumn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Integrating for 3652 steps, 3652.422 days, or 10.0 years.\n",
      "Total elapsed time is 9.99884460229 years.\n"
     ]
    }
   ],
   "source": [
    "mycolumn.integrate_years(10.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 287.84605967]\n",
      "[ 229.43643402  252.95019305]\n"
     ]
    }
   ],
   "source": [
    "print(mycolumn.Ts)\n",
    "print(mycolumn.Tatm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.477374247427 0.477374247427\n"
     ]
    }
   ],
   "source": [
    "(e1, e2)= mycolumn.subprocess['LW'].absorptivity\n",
    "print(e1, e2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "254.865280431\n"
     ]
    }
   ],
   "source": [
    "ASR = (1-mycolumn.param['albedo_sfc'])*mycolumn.param['Q']\n",
    "Te = (ASR/const.sigma)**0.25\n",
    "print(Te)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Check numerical versus analytical results\n",
    "\n",
    "Use a tolerance value to test if the results are the same."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "tol = 0.01\n",
    "\n",
    "def test_2level(col):\n",
    "    (e1, e2)= col.subprocess['LW'].absorptivity\n",
    "    ASR = (1-col.param['albedo_sfc'])*col.param['Q']\n",
    "    Te = (ASR/const.sigma)**0.25\n",
    "    print('Surface:')\n",
    "    num = col.Ts\n",
    "    anal = Ts(Te,e1,e2)\n",
    "    print('  Numerical: %.2f   Analytical: %.2f    Same:' %(num, anal) , abs(num - anal)<tol)\n",
    "    print('Level 1')\n",
    "    num = col.Tatm[0]\n",
    "    anal = T1(Te,e1,e2)\n",
    "    print('  Numerical: %.2f   Analytical: %.2f    Same:' %(num, anal) , abs(num - anal)<tol)\n",
    "    print('Level 2')\n",
    "    num = col.Tatm[1]\n",
    "    anal = T2(Te,e1,e2)\n",
    "    print('  Numerical: %.2f   Analytical: %.2f    Same:' %(num, anal) , abs(num - anal)<tol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Surface:\n",
      "  Numerical: 287.85   Analytical: 287.85    Same: [ True]\n",
      "Level 1\n",
      "  Numerical: 229.44   Analytical: 252.95    Same: False\n",
      "Level 2\n",
      "  Numerical: 252.95   Analytical: 229.44    Same: False\n"
     ]
    }
   ],
   "source": [
    "test_2level(mycolumn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Integrating for 3652 steps, 3652.422 days, or 10.0 years.\n",
      "Total elapsed time is 19.9976892046 years.\n"
     ]
    }
   ],
   "source": [
    "e1 = 0.3\n",
    "e2 = 0.4\n",
    "mycolumn.subprocess['LW'].absorptivity = np.array([e1,e2])\n",
    "mycolumn.integrate_years(10.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Surface:\n",
      "  Numerical: 278.53   Analytical: 278.53    Same: [ True]\n",
      "Level 1\n",
      "  Numerical: 223.20   Analytical: 243.87    Same: False\n",
      "Level 2\n",
      "  Numerical: 241.15   Analytical: 226.61    Same: False\n"
     ]
    }
   ],
   "source": [
    "test_2level(mycolumn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Integrating for 3652 steps, 3652.422 days, or 10.0 years.\n",
      "Total elapsed time is 29.9965338069 years.\n",
      "Surface:\n",
      "  Numerical: 297.52   Analytical: 297.52    Same: [ True]\n",
      "Level 1\n",
      "  Numerical: 234.30   Analytical: 263.52    Same: False\n",
      "Level 2\n",
      "  Numerical: 263.52   Analytical: 234.30    Same: False\n"
     ]
    }
   ],
   "source": [
    "e1 = 0.6\n",
    "e2 = 0.6\n",
    "mycolumn.subprocess['LW'].absorptivity = np.array([e1,e2])\n",
    "mycolumn.integrate_years(10.)\n",
    "test_2level(mycolumn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.652920962199 0.652920962199\n",
      "Integrating for 3652 steps, 3652.422 days, or 10.0 years.\n",
      "Total elapsed time is 9.99884460229 years.\n",
      "Surface:\n",
      "  Numerical: 301.92   Analytical: 301.92    Same: [ True]\n",
      "Level 1\n",
      "  Numerical: 236.57   Analytical: 268.24    Same: False\n",
      "Level 2\n",
      "  Numerical: 268.24   Analytical: 236.57    Same: False\n"
     ]
    }
   ],
   "source": [
    "col1 = climlab.GreyRadiationModel(num_lev=2, abs_coeff=1.9E-4)\n",
    "(e1, e2) = col1.subprocess['LW'].absorptivity\n",
    "print(e1, e2)\n",
    "col1.integrate_years(10.)\n",
    "test_2level(col1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Integrating for 3652 steps, 3652.422 days, or 10.0 years.\n",
      "Total elapsed time is 9.99884460229 years.\n",
      "Surface:\n",
      "  Numerical: 307.85   Analytical: 307.85    Same: [ True]\n",
      "Level 1\n",
      "  Numerical: 242.68   Analytical: 272.49    Same: False\n",
      "Level 2\n",
      "  Numerical: 276.56   Analytical: 236.57    Same: False\n"
     ]
    }
   ],
   "source": [
    "col1 = climlab.GreyRadiationModel(num_lev=2, abs_coeff=1.9E-4)\n",
    "lw = col1.subprocess['LW']\n",
    "(e1, e2) = lw.absorptivity\n",
    "e1 *= 1.2\n",
    "lw.absorptivity = np.array([e1, e2])\n",
    "col1.integrate_years(10.)\n",
    "test_2level(col1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Conclusion: The Two-level model works"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Three-layer model\n",
    "\n",
    "Extend the analysis to three layers.\n",
    "Start numbering the layers from 0 to be consistent with array indexing.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "T_s, T_0, T_1, T_2, T_e, epsilon_0, epsilon_1, epsilon_2, sigma = \\\n",
    "    sympy.symbols('T_s, T_0, T_1, T_2, T_e, epsilon_0, epsilon_1, epsilon_2, sigma', positive=True )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the transmissivities $\\tau_i$ for layers $i=0, 1, \\dots, N-1$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARkAAAAUBAMAAABL11ROAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAIma7zZnddlTvRIkQ\nMqvFy5UvAAAACXBIWXMAAA7EAAAOxAGVKw4bAAACbklEQVRIDc1WPYgTQRT+EpPcutmcgmB7e4VY\na22xwnKdR0CwUUIUiT9YBMHqChULtVJQUOTgUgnXSBSDhSAjCCIeXgp7w6GIjWdxByLi+WZnxp13\nO7lsGskrdua9931vPt7O7CxQCDEZNitl7JsMLcD0XlJydKSa22IkhAO8Re6P9iQjiIBifRj2pkoE\nX9bFMAiPe23lf2j85omhHme8ACqyP07TaoCrwpnPBE1tTI2rRjHeAfszRU3gv6spdXDFLE5j8GrN\n8vKoYYxcvXn/PUzX4IxyHYfTHJY7wMKCCeRRIxm1SDN4bVOGj/4joHr8lJNRjHA+RRd+0qbuegMd\nyaFGMmqNpibkUTNzDXiKz06G1wRpDeZaZGf705svw0oHyRZ8HscX4ng+YWV3sc0AFSGrxvHc/Tge\nyLljF7+RS7SOYf3kM5zGJ3l0MozaD6nG2O4jK1jt46H2c/RGMrQaIuXpzb0VgYuY6ag1OIPUWG9q\nD316SPWZ/GokYzw1m0nxO0KtwdX4TXsXlwakRuCcQiJHbyRjPDXq3ZsluJpihJ5emwZ/EQH15oGO\n2GqCjRRmzSRjiJrVyMKl02UBgapsqTSuptzFYxVPnk8ut2nf/NIRo6axNB+iJXSUD8TIqil/21hD\n5RJHas//ehA4YFJGTcJAqb/9ZqjoM0V4oyahevIUuE2dKcqZ2gp23Y2mqN/1+yrJGXQtbLs1p+r/\nvjev7XJV2+Fz7672A6a4y1GW9/HWDaFczqBbE4csHE17PcEDynvrCiaxwok/bUcycAUVbmlry0FI\n/ihy/m2FLv5OseJOSWdul2xtIaTHJNgs8BegurIq71jkKAAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$\\left ( - \\epsilon_{0} + 1, \\quad - \\epsilon_{1} + 1, \\quad - \\epsilon_{2} + 1\\right )$$"
      ],
      "text/plain": [
       "(-ε₀ + 1, -ε₁ + 1, -ε₂ + 1)"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tau_0 = (1-epsilon_0)\n",
    "tau_1 = (1-epsilon_1)\n",
    "tau_2 = (1-epsilon_2)\n",
    "tau_0, tau_1, tau_2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that if the atmosphere has $N$ layers then $\\epsilon_N = 0$\n",
    "\n",
    "Define the emissions for each layer:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAS0AAAAcBAMAAADGnPRhAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAMkS7zRCZdiKJ71Rm\nq90icBAQAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAELElEQVRYCc2WX4hUdRTHv3d2786dPzu7ZPRQ\n0d4eU0LBFEHEfaqHEK/QKvbi9G/FsJwXiyxknqSndiBKBCNJfApyCRIxsqEIMwIXCnywbf1DPkmt\npYv9kemcc+/5/e7O/nZ+87gH5v6+vzOfe+bc7+/sZYF8DOQ3Lu0FzE3rjPKJPsivE08RL6D3V+6q\n8q3Lkg+YO8OPPX15AVNqd999LSWLUiWwPhb3ePryAqavg333tZQszHCZQdvLI76+vID2FXzRb18u\n8jiX2au1gJavLy+gtUphv325yBpZFVGBqZOznS0//11tLteXF9B+on+2zP01u731vbevXmSxDpTo\ncwz4DNi/C8v15QW0r1IL5asY+q7h7asXGc0D5WlgBtEC8M2ZyQP79QcWr15A8WEyv47qmsnJhbbm\n3GtP8hYwliBoo/ovcJBq0sE6wguYex4FLk2jmMh4mKxL9CRfAn4AAmCYpixGdesVVwk/YO4aBeYS\nFIBv79M59Iqe5CHgZbmZhkJj5+lOp64bs+aA3bPtNO0maVI1lAxf7XRoULrDkuFHv6RfZuTaBOsl\nQUORRXBzzY13dWNXCxQ+AYpHW2Szk+RJzYJJPEifx148v6mZJe2SI1fHwIcbE0OOjWKbgDQUWTw/\njUuqc6sFpgi9iK8AN8mTmgWTD/0KhJ+iOq9Ju1oyontKzeFxQ07FuCkgDUUWdK7k4pJgQFzC3L7j\nuI6RhCfAQfKkIhjnAkziPDA0jgJnu0JIcWl44Vi71gzuGZKeKPXLHvV/wOGuArxlQFzCyYkkvM2v\nFzfJgxi8RhcwKX2NtFG8SvuuYFJcQnnDBGpx9CeUJL+kL3vU0W3g6a4CtBVAXALNDz1Zrb0MKYNY\noZ8Ek9LX1CgG66wXB5PiEkZIAfQPj5JjTfzEKXvU9J4OTnFqcTCQusQPEtxFjd/oLlIGMe1LHpnO\ncU+Chx2jwaS4hMEZ/rVaw5A0H/SqAAbk+I++TXIz3oup1DjJXDAgLgGrE+y4h3IjJcM7OYrlVu4g\n7YvIhOdrYDq6To/lJNklFE4hpD9IJeWlyqP7/vY7l2k841VUcuebH9C19htdbAggLlGRyxfIuZE4\nI3/nPkyUnupsa2hfRMo5hnvfYshJskvAk+80UK1T7ykJet+XY8pLDGyeUAmctVIV+UUuSaSTJrLS\nTFOLrqlfkiK/snCS7FIaT+BxlXiGXGrpLpijmdcwWU3IfJFLEm/w+yuNkor8muvrnMk7SXZJInj2\nyw2ZBP8zEpzQ3VnQYWYRZsboXlbrUmnfK/rNLhX51fZ15A9T1Ekal8qdzryWqM6QuqK7a3hBJapG\n5UTOJZttW2lV5YTVRrWNsiLvks3WRkk/p/sfL8Qq3WvOJTeg2ej1+y6/9evcmnfJpi+yrNRtYmWo\naJ308fnK6MZ2MZRYvdLU/5AdX4C+IcH2AAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$\\left ( T_{s}^{4} \\sigma, \\quad T_{0}^{4} \\epsilon_{0} \\sigma, \\quad T_{1}^{4} \\epsilon_{1} \\sigma, \\quad T_{2}^{4} \\epsilon_{2} \\sigma\\right )$$"
      ],
      "text/plain": [
       "⎛   4      4         4         4     ⎞\n",
       "⎝T_s ⋅σ, T₀ ⋅ε₀⋅σ, T₁ ⋅ε₁⋅σ, T₂ ⋅ε₂⋅σ⎠"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E_s = sigma*T_s**4\n",
    "E_0 = epsilon_0*sigma*T_0**4\n",
    "E_1 = epsilon_1*sigma*T_1**4\n",
    "E_2 = epsilon_2*sigma*T_2**4\n",
    "E_s, E_0, E_1, E_2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the longwave fluxes incident on each layer, $F_{i}$\n",
    "\n",
    "Note that if the atmosphere has $N$ layers then $F_{N}$ is the OLR (emission to space)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAACaoAAAAcBAMAAAB4hgdcAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAMkS7zRCZdiKJ71Rm\nq90icBAQAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAQ20lEQVR4Ae1cXWxlVRVet9Pb6c+dTqPEB/+m\nPvggBiVBDAkSimmMiRmoCT+CMVNRhmBAGxMgqNGbkBgfCNNEnZFE40QkPpDARIMSiHCjRNBgIJGI\nEXAKChIxCMgQ/0hda+291v65e5+zd29b5qHnofecfb79rfV9a519zj2dDoC/ne4fNOyX4hoohk7t\nGhp5IwZKlZXiajS0OtAKsNE2G1cqojRuKR/jSp0uxVUFLwRvifDC2AzbkgRKSUtxNYI2ETtzooys\nFFfGZlH3L1XBtwZcqqwUV5VlqwOtgEorS/lKVWw2H8UtdboUV6qlCrcVwt/4BEpVleKqFI0KfpMS\nXFq4qpXilLhgp/v9k2FVK1VWiisQrpBWB1oBlmqzcZphy05p3Baa8HSp06W4kH1zjrZEeE1qW5JA\nKWkprkbQKNgVmtxxj+7XFq5qpbia3HYfOBlWtVJlpbhNdaDUos3GlYoojVvKx7hSp0txVcELwVsi\nvDA2w7YkgVLSUlyNoFGwp9HkcV1NOneVrWoZ3OSgIJdePwN6W3pVKyJ9yHLelOH2h7MJECijzJ/P\n+xlcVbJDpBkHHK4VYKEbw42WPIXOxB2pKhmnnSl2L8QVaZFOKAJLiw0FpoFNED5aEiMkAFllGdIh\nB0JcTbWrRJcQs5ipAaZ4JUz+55zj/3zqgtWpbvOq1oK7cUhvauAP0eChW59aP+exf8NqelUrIt3T\nZ9JJ99QZBfEP4wToXIsynd6Cq0lWOaHFAWgHWK5WokbcxpInyua4G6tKi9NWSa5yRVrAdkIR2LaY\nxuWdTRS+sSQ2I4EhZc2kzoEkrqraNaLLiFnMs9gWJ2BqFabXYGLwUMuq1oL7qBOc2bsHHwsvjM7d\nAnAHwNW9fnpVayWduQ2gu8Ck433+aPiRSoDgLcqUsQVXk6xyQosD0A6wXK1EjbiNJU+UzXHH+4Rp\n2lJVaXFa6dK4Vi3gxWwFey2mcXlnE4VvLIlRE0gqayZ1DiRx430HSO9t0Pl2YhVzHV7Qy7AHb7fL\n0Hv3Sm5V+5nJrxnXW06rALDTu78/jqva1IqF2dFHYfI1gF9cAslVLUs6Y2k+dQ09Xx5hzl/VJeCh\nm5UhcHMc8JN1pM0OIK4NIF6MhKty2rOuNb/KqhQ57cxLVi6rRerot2IWLLZGVfO0NxteIzybhGQc\nJVF08WQTSCuT0WZV+Ftpe+klcW0xi5xPiW4jBufQIYDpY/B2gEeOwe5PHjz42sCrmdu1HjbjZvtu\nQrhnpwN8bQmfq47ak2a0M4DefwGu/enBL1xtT/gfWVKxFnbTqnYZz3nen+nvpxPwEM3KELhJDnjJ\nKmmLA/iyr9kibbPRcFVOe9a15ldZlSKn1bx072a1SB39VsyCUy22VcKzSWjGG2idrPNpZXZ0tDZq\nj+kWgSrR7cTq0OwA9i3BHMDxJRjjb6NB2eTAdloz7mZBD30Giwp82J43ox2APbguzePjIi55Q1uW\nNCzMOM6HybWh6XYgnYCHblaGwE1ygI3nZCl4kQO4ajRb5Fa1kXA1TnvO4W5zfrVVKXJazUv3blaL\n1NGsaqYVs+BEi22Z8GwSmnF96+SdTyuTVW2UNiqISatai/MJ0QXE6tDE6fAbrtQd/POXrx8DuPSp\nAR90P7e+jt8NadNVgV6B4ZbE4bdZ3brfe1L3veks6KA9o6TTazTSO/cJfJCTeRI8SxoWZmIZGXbT\nD9lEBR9rqCABgZpPp2woCQQog+fAEC6brE5n4zlZCqqkvgPOdZMX//QAKuziO9fXl/msekGvR3Fj\nK0txQpNNXtnD5Dmy/GiIW1sVNaXJaYybxFnzslrcLNcJWfB2Cs8moRmH7qt6z3ntR1vSvPNpZTo6\nQhsVxCxwPiG6gFgd2n0UPkudSa+27Db2A7vzzsvvO6tv9tXDJtyH7Dz6OHXeO9Ac+RsoPGBPKSm+\n05NN50nwLKmWgLX0FpBh10BoAFQFD2kodlQScGjc85QNJYGnlaEJl01Wp7tkKbiSeg6IcDrvNgdQ\nYZ3n3vvnbxiEekGvR+1WilOabPLKHiYvgfizIW5tVdSUJqcxaAon5mW1uFmuE7Lg7RSeTUIzDt1X\n9c55veikpHnn08p0dIQ2KohZ4HxCdAGxOjR2As6gtqRXW3Y7hI9ru4+sQvd26L1kx9TDCAen4OOV\n4F4QBlwiiO6GG2RAp7Og+2JSfKdnN553+INLBaRaAtYys4YMs8qDX2dpH7MzWzoBOcufThkn0Vnw\nlCFAGSJcY7JEYjaZ7pKlcRmlt5p2UzdlwHw6AAmj8sBlx+ARC1IvPCLCTX3milac0mTLp+xh8kGC\nDfnVVkVNiZwuqIial9XiLHetGIO1atspPE6C+sps4kfovox6Fde+lZLGznM7MGlamY56pLVtFMfk\nVg1jFjivZXKiY+ImMZ2XYT/FpFdbdjt+1Xfh1/BzmFgAXPNw+9bi4t8WFxkW497yJ3A4fMbrnv8i\nbn+Z2/PaLYPe6syj8XQWtC8ePS4FBJo31d+z0EY6tbh4/g8XF5mfhXdeQtLpefx3I5TAi+cBqaDs\ncPPzdwnwKe+Hc4CS6Fyz5inzGUJcY7JMEibgkg3Scg6o615quOsAUh56xn4/Ged7EeG+A49BG05o\nIF0+n91LPsyuMb+qquSdLqmImpfWEljuOiEC26pts/AoCe4rNNnzw3PfG/Uqrn0rJY2d53YIG8aR\n+no9Umq3mjaKY9JKEjVpm/Np0TFxk5jJV+A5isrfzs3D1a0XLcEzsHdp7wB2r+Ep2vTOwDhzMyMc\n4HOX4rAwsk2feRHMzoMslTrdCfJJ6Y2WWXp5Xr/zrwJSvbG4wpBw2SQ7c5xOQLD86ZRREkAPf6oM\nAcrAOL6TtibLJCaITHfJ0riM2jdINOTHpGO7kUXmtiflgf8BfNWeVS+YyFhJuOdhXz+PM5OVJls+\nZQ+Tt7HNR0N+tVVRU6oqwmmoeVktznLXihHYVW07hUdJzNJFYDbxI3RfRk3r8JcS7Vspaew8twOT\nppXp6AhtFMeklcToUPYC57VMTnRM3CQGn1r3U1D6dm4frvANW/dl/Pceh+ZgHEd5Uw8JZ29m/CYO\nVzXFvWDB+LEXYQfm4HY7otNZ0PuiUX5/YpZemjc7j0ttO6maxMLH1pDUf0iV7EysdALmnP3plFES\n3NqaBB4rA+HMnbQ1WXd96HSXLHLqqPcGyY9JELMxwNz2pDyTLwN8xJ5WLzwryYBn4NB8HseT3els\n+ZQ9TN7GNiwULZNfbVUCp02vtVfEJKPmZbWo5eYVL7diBHZV207hURJ8ERhR4kfovoxyxc2XEnFJ\nSxo7T+1gONPKdHSENopi8koSxXSLQCTaOa9lcqIjYu7tiFjBuIDwqkbvRezDFR7jfWJ2cGAJ3mqX\nWQ1ivnDPrBEdz8NVTXHeC8/xR4EeE/5uwrrpLOiBaJTfn5ill+bhinKigFRLwFp6CzjPf6Eo2RGf\nWz7CX1eYc/YnvxliZZwE7akyn4Gd4jtpa7K5KnGyFFZa071BCmISxGwMMLc9KQ/+g+nObfa0euFZ\nyQbA/Ut5HE92p7PlU3bntI2rH0351VZFTPF6rb0iJhMtWFaLWu51QgR2VdtO4XES7j9gEj9C92XU\ntA5ehvj33HTxYN9qSYecx3YwRqWV6egIbRTF5JUkiukWgVg0X36MFnlOdESMqLwYjPo7ojkX5dqH\nq1OX4OMnYHZl17HJZzgC/pAgjLNlR9wSfQNV3GFB428gb4PuviX4qx3R6Szoqmh0F31RNfcRmofL\n60oBqZaAhU+sIsfUsmXGD8nOjAwn0H3VYXmPHDDKOAkyWJXhGWUgnLmTtibrrg+d7pKloELKDtCA\nxDywwAf6gwD2tiflgbPh5nkLUC88K9kAqkAWZybr6Wz5lN0lX5NfTVUoJTHF67X2ihgtWrCsFsfu\nWjECu6ptp/A4CboIzCZ+OPdpXEZN6/Cqpi5JSYec1wsyoQw5dXSENopi8kpidCh7gfMqz4mOiPEq\nlNVFiRWMvw/4Ei4HH1jfvyIPV2OPP4hr3fRK98ob6ELnzXpocLbsiOP3aorjx3k74bSvrOCz2o/C\n6XDNnfsHAOeFo9+84NXHaeR+iobzALDC7aSiZeKPrz6J96k5nDy2hj/sJtmZQ+kBL4EXVRxBAmWU\nBK1qmgQCQgfMnbQtWXd9yHQvWYpqSa0DR76MQybm7LN0WjcG2NuelAcuvv7bAhAvfCvZAGqFPI4D\n6uls+YTdS74mv6qqoKLQaXP7bq/I1PXHxDzkyGoRdr8VI7Cr2nYKj5Ogi8Bs1g/PfRoPW4dXNb54\nyC8p6ZDzujIklCGnjI7SRlFMXkkoX8fuXYOxaFNsAg+Ljoi9ZxhJ2znUO2p+kUZM+/ThCh8L9s7T\nkGw2iD10welZTbfZvu7SDj76vWIHgumTa6lRt/RCb9kC6CNLKloM9mb+iP+oQrNLJDDTNzODn06Z\n2zOAgIEfJ814U7Lu+pAqmTkmWdwPSCfm32xO08+73a7sdegBmg+GyhN64VkZ/e1ciAsCVjpdk19N\nVVBfYIpnYXNFLvf+l8AGLSG7acWoxVzI0C5Tta0SHifh+irwI9063mUYuBQ7r+2QUubWHek4eRgC\nnWfOhLOjNgpjeq0azGpxPixT+uLWpAJiIDD+bcH0vFXhPVy531yYc5+wEPPhOafrBp5xhWDYrPsd\naDB9ylyaAMGo9/XxPfAuFy5L2u07UO6v2zW7IJRJYMqfLvtOmdsz5wIGfpw0403Juusj1HpEwgWk\nu87GX73Ktio77tO77cXlCb1wVo6tjs05AnyW6ftHQcCG8oWzTPI1+Wn72eBNVUFIYIpnYXNFzrhr\nydOW7ZqQ3XRCBHYht1N4nITrq8CPdOvkVrXIedcOKWVxe2ywjaKYrlWDmC3Oh2UyoiPiJjHjK8Cv\npKgjvIer68y/MvH6xN/1+usef/ws/wD/fmnm0WDAHlyYGgR3P+h87N4zPUgRaXeBZ+zqexNxN8hO\nT5kELtFjb8cpmznqDQ/tag82JgtpEptsTNo5jr/VtFt3Rfbcp3fbay6Ps/KUe3/iX+uOi/f8gDRQ\n4XRVfjVViVLE54c1GUqbKWffsT4nu/RZpEX+U6wQ3Fi1rRMeJhH1lUjLtI6/qh0VLH5GzmfaIUeq\ny0hmng0UtVEUM9OqdhEIRTc6XyPmwBJ0xIdZ93A1ddUVnjnxruu0r//jce/kjd4+7h4+nLykIiV2\njlt6p9fXX/KIikitldF/LBdmp5wmgYEeezuqbPLzr6944/Gu3kkbk82QRHUX7rvBedmTQf/T3faa\ny+OsxD8U9QmifT8gnapwuiq/mqpEKXqrWsZMO2GmvyeoV5EWWftCcCaQqdrWCQ+TiPpKbMm0jlvV\nwuQj5zPtkCHdYBtFMTOtaheBUHSYvGiW5SwibhLzW5z7hJ2fe7hSetnJ3TUnB4Jo+Oz1kyez94Mi\n0ocs501J7nAwkwCDcspCBsjcSfEvxQYRMnUoyUbnnoZPRyPRYea2F6Hw78Qan9AcPA44SvLEms1v\nlKoUVmR2bmLOKSsshHTCSSF8pCTSX0oASpzP/Q/fG22jkpgbcr6E2Ih5GltBvxJmHq78bqH9zJIa\nw+qOM0tvHcmo6FJlmTvpaOEffnC+mSBz2xuaVGpla8Ah5uaB0vyaWcKzpRXpPfzFcOJ2Hm2F8Kr8\nM19KqjiGwG9UGw0lUj3QW8UpM8vV83Ym7Diw48COAyenA+Z77Y9PzuR2stpxYMeBHQeqHWj6nUA1\n2c6EHQd2HNhx4GRw4P/J+Io/pbSwhQAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$\\left ( T_{0}^{4} \\epsilon_{0} \\sigma + T_{1}^{4} \\epsilon_{1} \\sigma \\left(- \\epsilon_{0} + 1\\right) + T_{2}^{4} \\epsilon_{2} \\sigma \\left(- \\epsilon_{0} + 1\\right) \\left(- \\epsilon_{1} + 1\\right), \\quad T_{1}^{4} \\epsilon_{1} \\sigma + T_{2}^{4} \\epsilon_{2} \\sigma \\left(- \\epsilon_{1} + 1\\right) + T_{s}^{4} \\sigma, \\quad T_{0}^{4} \\epsilon_{0} \\sigma + T_{2}^{4} \\epsilon_{2} \\sigma + T_{s}^{4} \\sigma \\left(- \\epsilon_{0} + 1\\right), \\quad T_{0}^{4} \\epsilon_{0} \\sigma \\left(- \\epsilon_{1} + 1\\right) + T_{1}^{4} \\epsilon_{1} \\sigma + T_{s}^{4} \\sigma \\left(- \\epsilon_{0} + 1\\right) \\left(- \\epsilon_{1} + 1\\right), \\quad T_{0}^{4} \\epsilon_{0} \\sigma \\left(- \\epsilon_{1} + 1\\right) \\left(- \\epsilon_{2} + 1\\right) + T_{1}^{4} \\epsilon_{1} \\sigma \\left(- \\epsilon_{2} + 1\\right) + T_{2}^{4} \\epsilon_{2} \\sigma + T_{s}^{4} \\sigma \\left(- \\epsilon_{0} + 1\\right) \\left(- \\epsilon_{1} + 1\\right) \\left(- \\epsilon_{2} + 1\\right)\\right )$$"
      ],
      "text/plain": [
       "⎛  4          4                    4                             4          4 \n",
       "⎝T₀ ⋅ε₀⋅σ + T₁ ⋅ε₁⋅σ⋅(-ε₀ + 1) + T₂ ⋅ε₂⋅σ⋅(-ε₀ + 1)⋅(-ε₁ + 1), T₁ ⋅ε₁⋅σ + T₂ ⋅\n",
       "\n",
       "                    4      4          4           4                4          \n",
       "ε₂⋅σ⋅(-ε₁ + 1) + T_s ⋅σ, T₀ ⋅ε₀⋅σ + T₂ ⋅ε₂⋅σ + T_s ⋅σ⋅(-ε₀ + 1), T₀ ⋅ε₀⋅σ⋅(-ε₁\n",
       "\n",
       "          4           4                          4                            \n",
       " + 1) + T₁ ⋅ε₁⋅σ + T_s ⋅σ⋅(-ε₀ + 1)⋅(-ε₁ + 1), T₀ ⋅ε₀⋅σ⋅(-ε₁ + 1)⋅(-ε₂ + 1) + \n",
       "\n",
       "  4                    4           4                                ⎞\n",
       "T₁ ⋅ε₁⋅σ⋅(-ε₂ + 1) + T₂ ⋅ε₂⋅σ + T_s ⋅σ⋅(-ε₀ + 1)⋅(-ε₁ + 1)⋅(-ε₂ + 1)⎠"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F_s = E_0 + tau_0*E_1 + tau_0*tau_1*E_2\n",
    "F_0 = E_s + E_1 + tau_1*E_2\n",
    "F_1 = tau_0*E_s + E_0 + E_2\n",
    "F_2 = tau_1*tau_0*E_s + tau_1*E_0 + E_1\n",
    "F_3 = tau_2*tau_1*tau_0*E_s + tau_2*tau_1*E_0 + tau_2*E_1 + E_2\n",
    "\n",
    "F_s, F_0, F_1, F_2, F_3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now define the net absorbed longwave radiation (flux divergence) in each layer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAACFMAAAAcBAMAAABVUX/PAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAMkS7zRCZdiKJ71Rm\nq90icBAQAAAACXBIWXMAAA7EAAAOxAGVKw4bAAARx0lEQVR4Ae1ca4xkRRU+PTs9z57eUYk/fG37\nwx9iQAxqTMQ4mvGRGKRNBAVDduSxBCM6MRFFjXZiYvxhmEnURRKNG5CYSIITDWggwkQJItGwiShG\nwB1REBURlCXiazx1qs6pU3Wr7q3b2yv82Pvjdt2q8/i+r86trnvnAaCP0/RFTbvUriZEZWhPpef/\n3dGOVTtr4VJKs9ROApc0xsI8JpJSt0nbleiQtilFEnuPpWocJLouxVJqF4U/5svFo2UhSu3Kojmr\n24atzCdv3I5VO2uPtpRmqZ2P3NwaD/OYSErdJm3XLEPOohRJ5D+eqlGQ+LIUS6ldHH+86+eK27mF\nS0WpnQQuaHS//kwvFe1YtbL2EpfSLLUrUNabeMzrvrOp1QpJZ8DhSt0mbcf53edsdF1zWYqEQnQG\nEsmrKl3H3CjFUmp3DIBIQStjx++fPlq4VJTatQE4u/+ZXirasWpjrSQupVlq10Zi8JhPKfdrheR5\nErfUbdJ2AsA2pg5HHfnLUiQUwRNVquZjtx0pxVJq1za/sicFrYzTcot2bixbKjJ2c9sqQa7ZG2VG\nXpheKoqC3ulifiETW3dnAWRYaV/VbmXtJYYMTRXZNjN2k2I4v13JmOvIIEmbv0W6S93Gsysvi68K\noqZGKRKK44lmKqEcYQpXKZbQrqQ6oPXNQgrS6RKY++cbjvz9gbM257v1S0WD3WdTnCt9v4p6Nq59\nYPcN9zwNm+mloijo0oiCzvn9UZREX8YAGlgp189d+KpirZQbwCUADTTFvN5ucgwfkoyJhiaam5iE\nG8DM+oSI+uhpQcrLoi/fhD6mtOYOnL5cjFi8sIFEG+qmHKGOm2arLWw7aVdWHa1vFlLQnOaOwvwm\nLOzAzPadDUtFg93bq4yinpsx37ujvqsBbgD4YG+UXioagy5eB9BdoaDTI/qoOaUANLDy0bqH4euD\nYmvvZySGBppiXW83PRLDTKOU4eUZf9MdEM1NTNJ/Yzghoj56WpDysphd87EqrfNh/m/FiLU3Em2o\nhHKEOm6arbaw7aTd9KhqGPaMdbOQguY0vwZL+KW3Br2XreeWiu/ZhPV2PQyWPpx795dHcKmYX3dG\nrvcwzD0F8MP3QHKpyAZddGHOv8zshK6imHek0wOkATjrelYq5NIQFg41WjMw5YgSQz1NAHart8sy\ndO5piVOYNxS+uKmJDjIT4xFr74PNRHky6oli0Lr6aFEWc49rgFH7FwDfKp4a7YtEU6p6kxYIjVPZ\n9DfYNVUHjHezkILmtLAFLwK4ewtm33fgwFPbnqxquVmrt+uPlEfQdO4AnxniV9YhN2Z7O9vQ+xfA\nR2868OEPugH9kQ3K2sKsWSrOI59HtKdupwE4i3pWKszCDiw90WgtwLwnStxAU2qgwS7LULImJE7O\nb3/b44tbmugwMzGCOHC+vpmoWwIaiGLUuvpoUxaPBgjDC9zPHnlOfQXmiCZV9dHbIEQvN3+NqtTa\nNVfHeDcLKYinfUNYRsGGMEUPI56sarlbrd7uSuUQNoM7Fd7sBl0p4PKMN/sANza4jlSObFC5N4j9\nNPrD3E7F3XWkAbjBelYq5PRRXCoarQWY90SJO/U0fa3U2uUZSlazVEQSJ+d35jSPL26FRNMTI4i1\nc+doM1FeAmqJmqB19dGmLC7UCKP2WUNcKpqQiLbeGYkmVfUWbRCil8txTGVSUB3j3SykIJ7uIno3\n0PlH/9kCOPeBbbrofmB3Fx8NzCG3mnmtgEfS7nIasqfu1+5XV+JOdXzAjUgvfonh0XvjfbjlYD9O\nng0q80fsZ9YwQvBUyixMaI8/AEBDdPKsOK8fC1p9s4Xx1imtGJhQYYkVTQF3zrd3d9coA7vhNm/H\ndJAcsV2eobgHDEVijdnmnD1ksmQPRTSJBB0lpaCE3iETUBEQETxRPxnj2LmabFMWnzCY8sdfcW1V\nSISMR5wlqlQVolxAbRAiOMlRikXbOahBdQggw1yij3ezkIJ4usjEMq8L3DF1jWu85IJbXzeybSm5\nOrsznJ/5OHmgLnxxUB3f7oYk6MaaGIsfJ88GDdn3VjDEnm2JA8KCuiRVAECsFSvOK2Nh45WbjVox\nMKHiJDavg9wh4DoPv+J3n7ed7EavjXJ2eYbiHjAU3oqhyzllFr38kSOaQCxs8OcCKyaiJyp1oNx8\nNYxhx9PTpixeiWtB/ph6OkAsZBRi1lbGmKhSVWZ7HISIgHNo9SRfAou242FdHaI8MZfotFS0vllI\nQTydboKZ1wXu2MCNxexVm9C9HnqPuz4pucgOTsKNANv9mSPg3WTCXXEFd4g71fGtcVB8T+IO8jv4\n2mFB0JD94g5G6EscFNK0EZ090gB4VLEXMjIWNv6Cl5EGsVYOmJaAJDavg9wh4M7bgrtdn/CJ7OYv\nvBhfxTi7PENxT0usMLtYnScYTfIzQZTkZCToxCmFDS7Wp5lgnqiIoNz8UhHZmXlvsJPpydYag7Kv\nsKgs9i0nGbrO/uEAsSETSI6jHFPGmKhXVYhmEXZWXEKOphCqHFo9ky8nurZj0XR1CCDKmsipjSUR\nGVdvFlIQT2eacfO6wB1HLv0q/AR+YL4h7DfPl1ZX/7S6Smax3fN/QwustcPdSPedj+Hx++Wlp67e\n7m0umlkA7U51vC/uPTJ0qcH4zY+WVnzydND51dV3fmN1leLTQtl5HEMsDPBHsQbAY28Cw8KgwyMN\ngA2NhWcvpE03H8wKua5hn7emLIFWHpiWwGrnaQq4iwBoxfZuGD+0+wrcQ9sSsssw1O5piRVml3Pu\nb0zPf9YTJTkTiIUNwDTNiCcgIrBbMBmhHc17g90M12SbstgYeIbc8rNPj8MeiSETSK60lTEm6lUV\nojmEnct2TG4VzS4VVLi6N5r+nOi6TFg0Ux08hQIolzNTSsmbhRTE08OGwsIOnuw24Nqzh/Ag7B3u\n3YZZ02sOWWjIzi6Qxg5whyB2/OCC9guvORv6A+D1R9x9HeugN+AFreTWb9T5R0HQcKGUpQJD0cHo\n7FUagDPlR1VzKWRkLGjQA2iDVg6YloAktg+2tNERcP8G+LRLIHzIzsph7B6BfSNgOzO/fEgQ6hD3\ntMSE2Xq6WPS1w7Eqn55oONmCBD04pQLSXzORzHyavZaqAyZghmUyEvOetTOOanqytcag1I1I35nW\nv3qeXTN9HklFchzlmDKGG1hy85Ugsy0FFCEE2uCYXBxNIdS9rcrExAOZE10dAogsEjm1sZpBNJf5\nkVIiBfF0pglmHhvdNgDfWnSfwB+hbizDNPbSId7Gzi2Q9HYDlwqx+7Mzxo+9aLZ/Ga53PeJOyU+N\neumJj1Zy8usP8AuvOWjIfmoHg+o9FaOzudIAHA792Cd5ZUw3FtcA3yw0aOWAaQlIYqJpNzoMbg4f\nAt7qEggfJYexexA2BmKXZyjuSYnV6wOOldpVCFdPNJps9jaWnJLZYBftKoiA2Wv5OtBuUoqKqBHL\nzHveDoPhIdOTrTUGZW9EKouNAfmmTzfB3Mi+qbMVWJEc3Tgmj2GX3T75SpDZziHMLRWEUOWwWErL\nhDiJaLo6BBCZMAOlijZWM4jm1ZuFFMQT1fHdW2auaBuA1/i13t/eP4QX4L6BDvE2do41+eFSIXZn\nOGP8MEril+GjrkfcqY5vj3rpiY++PN0M4B/1NgcN2fdWMOiebRcZPxid7UkDEGtiRVeSV8Z04/kA\nH7cPiTVaOWBaAgJjH2xRMA8Of+Otc51LIHyUHOQHtw3FLs9Q3JMSq9cHHAvnOH94otFks7dx5ZRK\narqDiIDZl/o60G5SioqovfMWj+btLFaZnmytMSh7U1BZYCFmj84azI7s2ydbgRXJ0ZVj8hh2WcD+\nrpHZziHMLRWEUOVwb8IKy4R4iWi6OgQQmTADpYo2VjOI5tWbhRTE089NsDfivLptwMlDeNdR6K/v\n2Zp7kPJob2PnWKPd0DyAiN1BtgaYug66+4bwB9cTJr806t1jnlPMl6f1wzVrvSBoyH5mE2PMr+HJ\nHYzOXqYBsC2xtxeWTPdJGdKNuT9ebZ44G7RywLQEJDHRNILhz4ecdPB6uHJgOvAQPkoOsjMqsl2e\nobh/xsxQLLFiyLHo/VIz0Wiy2VsjFjb4Guc0HDAEaF+q6oAJGD+eDEXUiEXznrUzjiW1JjrYHwua\nsqC3PPtXjH/luPLqg7+xiF0FViRHF44pY46oqgSZ7dzdkFsqqHBVDlKPHuxLysTyYdF0dQggMmEG\ndqlofbOQgnjCn5fOv3r3zHXeBkzd+2P8xllY715yhak6OtzsWjvHGu2Iktidytb4ecqn1nEV+qbr\n4eK47NtnbgO8Kez94llP3mt6bjPZ0A8A15zmoMx+5tdP3o+r/DI6T+3gyR2Mzl6mAThTx2r+Y1uS\n9zGDpXJM7+7u7jRqxcCUBObn0o4mLRUC7pyPfZmzsJuWg+xMAbBdniG7JyXWDDlW75DJ3Eg0mmz2\nNr6cUtjgG4pDjijtS9GGRWACxs9NhiYq8563u+qT6DpWWdxlkvYfMufKcdbu7n95aqgCK5KjS0DU\n3o+GaFAJTDSHsLpUqMJVOVqWCanCounqEOWJMjNQObWxmkE0r94spCCeaMkwAf02AL8T9g5MFx/s\nba8Xd7jffkm6q/5Iuk0DNyn8kj1wn9txZkEvyBYEemvOwHxkgzJ7a3slfTxiL+RM96W5ClIJALGj\nxgXqn0osjsKxylWNVgoYSyASu68LCifgbHDlRh28I4M7dPIsw8BdGAa8QTME+9uajUR9iWvsBClI\nadnIL2vgbzP2cdE3B4tgr/AcgsrMe2w3M3ieBMBGm7K4kBy/r93T7bTkfqkgLzsdQlTdNSHRCKHX\nMZANbOFGObTUtWUSqRJVhweUyhkZ+1s5mB8qJVIQTwsDJ5zaBtgnTS/oe30TW5mlorcWWPX9T0AC\n93lXQhD0qqeHl8NLfaBs0O7IG+X+XExkDlIJAB0A4PQb/VZiPhyqXtVopYCxBCJxTQ0oN5Nufs0l\nndqcWnZN/AjWDbwWhoG7MAx4BwxhmqahkWhusjF3kNIBucZBpX2pbbMIbgTK5h2tA/B7Xn+2BMBG\nm7J4GzluavdkOyN5SJSng4mavbjbPIdEI4Rex0A293eOYQ5EJ/OqWgQ69I5UiarDAwq9rqJIkbFP\nFOhOpUQK4omflvTyfzm9v04Kip1qqbhZ27xOX+CvWS8eDjrcxbtTneom6Lzjltcok6Kg3RXy2DNS\njtgM0MlQBsCLd/0t+R4xzjT8gg01WrEEIrGugTQ4zicTedIt3x1yJz5Vj3zbtNJBChi6v+NtJJqd\n7BAHA+GXJGpfyiJEDnwpRKN553H32TkS/spYeVngWz88uutRwOqlIAklDw15jImquyYiGiKExUNh\nIHvlCrcy5JeK9Aw7h0iVqDoiQJzE5YyMeQbZyn2aUiIFzalzyPX6RQjmL73YdaY+/FLxub/eqww+\nq9rYPHhQ1bgfiiR0A7xa4y5nd/dxbw1FQR3t6F97hOgkZhrA4mjJV9O2GGcahVo5CURitVRkwLl8\nXg78IxGFYWIM4WcUdVvFTjdzkx1aM5sXcLfal2bqwFp6otG8cyD3+X3QlQYtyqJ32IToRfGqlx5J\nKHloyWNCVFVCSDQs3LkP/ceXlw8Z3688IksFC8sD4WekSlQdmTuw9c1CCtLpPpc+swiF4MxVeoHE\nX+fertpWenqjSpfp4NW6MlgU9E7n9oWKe7UjA6C/PLNcNc71FGtlA7DEuW1AJU1Wjokx/G0lZ7oj\nN9lp66Ut11+z1wo8s0QDK4DfwvuDnvKy6BdOaykSC0OIZiuhHGFAzF7U7iW8faxKSXW0/4d5pCCd\nZLsarooeUNTKLJCRVctLXq1buk3OvPfTj7QJVqiVCykS139JeADHQ46AYW/TJ6trtZzs7ooLVr8v\n9RlLif70xwPv1Kr1k0LrUiQ2nBDNfHUXJs2YlZbJ+KpkEie7SUE6La4lDU50Tk6BZ5/E4fZ4ckxv\nmlyoyUSKN+WTiQrwrCM6KWKVOKSgk/E7ldETHRNW4Fkncd27qGPhvjg4Fu/j4DszPA5BMeSzjujx\noYlRScHjJeNxQ30i8AkFTijwTCnwPzkrnpMe7HSBAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$\\left ( T_{0}^{4} \\epsilon_{0} \\sigma + T_{1}^{4} \\epsilon_{1} \\sigma \\left(- \\epsilon_{0} + 1\\right) + T_{2}^{4} \\epsilon_{2} \\sigma \\left(- \\epsilon_{0} + 1\\right) \\left(- \\epsilon_{1} + 1\\right) - T_{s}^{4} \\sigma, \\quad - 2 T_{0}^{4} \\epsilon_{0} \\sigma + \\epsilon_{0} \\left(T_{1}^{4} \\epsilon_{1} \\sigma + T_{2}^{4} \\epsilon_{2} \\sigma \\left(- \\epsilon_{1} + 1\\right) + T_{s}^{4} \\sigma\\right), \\quad - 2 T_{1}^{4} \\epsilon_{1} \\sigma + \\epsilon_{1} \\left(T_{0}^{4} \\epsilon_{0} \\sigma + T_{2}^{4} \\epsilon_{2} \\sigma + T_{s}^{4} \\sigma \\left(- \\epsilon_{0} + 1\\right)\\right), \\quad - 2 T_{2}^{4} \\epsilon_{2} \\sigma + \\epsilon_{2} \\left(T_{0}^{4} \\epsilon_{0} \\sigma \\left(- \\epsilon_{1} + 1\\right) + T_{1}^{4} \\epsilon_{1} \\sigma + T_{s}^{4} \\sigma \\left(- \\epsilon_{0} + 1\\right) \\left(- \\epsilon_{1} + 1\\right)\\right)\\right )$$"
      ],
      "text/plain": [
       "⎛  4          4                    4                               4          \n",
       "⎝T₀ ⋅ε₀⋅σ + T₁ ⋅ε₁⋅σ⋅(-ε₀ + 1) + T₂ ⋅ε₂⋅σ⋅(-ε₀ + 1)⋅(-ε₁ + 1) - T_s ⋅σ, - 2⋅T₀\n",
       "\n",
       "4           ⎛  4          4                     4  ⎞        4           ⎛  4  \n",
       " ⋅ε₀⋅σ + ε₀⋅⎝T₁ ⋅ε₁⋅σ + T₂ ⋅ε₂⋅σ⋅(-ε₁ + 1) + T_s ⋅σ⎠, - 2⋅T₁ ⋅ε₁⋅σ + ε₁⋅⎝T₀ ⋅ε\n",
       "\n",
       "        4           4            ⎞        4           ⎛  4                    \n",
       "₀⋅σ + T₂ ⋅ε₂⋅σ + T_s ⋅σ⋅(-ε₀ + 1)⎠, - 2⋅T₂ ⋅ε₂⋅σ + ε₂⋅⎝T₀ ⋅ε₀⋅σ⋅(-ε₁ + 1) + T₁\n",
       "\n",
       "4           4                      ⎞⎞\n",
       " ⋅ε₁⋅σ + T_s ⋅σ⋅(-ε₀ + 1)⋅(-ε₁ + 1)⎠⎠"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "R_s = F_s - E_s\n",
    "R_0 = epsilon_0*F_0 - 2*E_0\n",
    "R_1 = epsilon_1*F_1 - 2*E_1\n",
    "R_2 = epsilon_2*F_2 - 2*E_2\n",
    "\n",
    "R_s, R_0, R_1, R_2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solve for radiative equilibrium\n",
    "\n",
    "Use `sympy.solve` to automatically solve the algebraic system.\n",
    "\n",
    "We will solve for the **radiative equilibrium temperatures** in two steps:\n",
    "\n",
    "- First solve for $T_i^4$, which is a purely linear problem.\n",
    "- Then take the fourth roots to solve for the temperatures."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABdUAAAA0BAMAAACeIczVAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEImZRO/dIma7q80y\ndlRAyO8oAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAVYElEQVR4Ae1cbYic13U+M9rdWe2uVmNDErBj\ntEiUxi1CW6stcShkKDg/SqkVSFNUk2bTpnabtPVSg2L/0rYuyCUYb+iPpmmJt39So6CuDCWlVqkW\nDCWtQ7W2Cc0f11tCEmqa2IrsJFKcbu89955z7rlfc2c/ZH29sPOee+85z/mcd9535u4BaDr2NnHd\nfEzL15LLu5uka8rV3Qz76/3dRL9usWcvXkum72qSri1Xtx722+eHyHaO3Kr1XIhuv5ZqfXeTdE25\nmstFbe6xAa2OfXSRyMJ53wFX648W1m+E6VNbcOJFV+u9pS3IapHJYSnQ7LkRJSm3tv25nXN1+7Yo\nhJa8TV6ma/VTc0o4GJw8fnrzC1/+PvyHq/XeMq/1jv7VGg+GEY/c+zfDWKCBp3fw7o1hOGWYqsUz\n8zDVXm29t76wevn0+bmZ/3a1/tgwq+rraPPf1Xkqq1GSKpyjLWHACHyHXB3NguHcrXk78QmH1btc\nxDwI8FGA98Ocq/VpKYc7oPtmUSxa6KzAkUE0Fw9beB6B2R/FctG4AlO1eGIZOmcjrPKwOwfjC7Bn\nqdtxtf6PZdaGFWfznQ2ceZYoSXmmLcxiwAh8Z1zdghV1kda89a44nKnyTecK9C4B/Nnkoqv1x0Xx\n3wL8uowK1Ow6Lkz1YfxYgQW+6xZqPB4GXt2AH5Rg/HwFpmrx15cBDg/BluUpgJPzMDk45Wp9cl6W\n8hR5oFeV610XK83RNIqS1CTTwoQB8+DQ7Coj532m5foqVQVxF8/NefvgGmLsP1uCmlmCyf8DePEu\ncLX+kDCay/2qE5e5hPIOmWvg1NvJop/wCa/xUFzuG3R+XILx8xWYmsW9z5paf9cQbFn+KsBTG7Dv\ntnVX6xOLspSnyAO9qlzvHNOLzaM4Sc2CwxhtwL7mK6DT7Cqj5n2m5fpqa6235+3+JdR8YIUMiM8z\nAPaiP/jm0Q++36z1FoThfL+91qcvDq/1Go/EZeg9TAWmZnF32tT69EDcq1Nrpgr6MDZ79OilJcN5\nos5tVsWDkNXXurf5H8KlEegoSSNIDmG1Afuar4B2Vxk07zMt11dba709b+Zz2B7lWjeL5kJpj5N9\n87LP8eOEeXm4D7efXnKjlLDz4tDERej8zDOOlwkc+oQbusgjMO+eg4xOB0uvIQzN+XNosVo6ZWt9\nj3ZPMSQDc9EzRw9v/14OVtk5Juwie8BxsrPKdThqp7Z2hEnaGkJBygSMKqDFVeUd+cyBYMIqo9VC\nOik0LMOEsrQ9b77Iq7Xu3g+TT/+90bF3KVQ09n0Y+4ifSAlcYIfg/jm4b+CZmcAxeQVlHob5498x\nV9NEp4elU6iK5tw5tFitdNZtrU+eVZPVgX2MMcefX9kwr5+3pD/YOSbsAnnAxiO3ch0+6xG2cAqT\ntAXxoogJGD6ZWIYGV7V35DMHggkLR6ssw4Rd5csAyzCBy/5lhLy11Lq5L+VjIqDNZXgFTpqJfYfn\nTDzsygvmzxLde3/OUPYghwA+BT1z4w8vvWSuhZaY4bqShMc8h/7CXFXsITDdj5RUOU776mFEAy95\ni8U+WumCrfXZBRoPP9vHGD7+iClxDr1MPNDhCV0H+LagjEqFSRpVtsZvAoZPJsITu6oSqr3zWSvk\nnXJqZbiEJDM+NFwrOpxkzwh5a6n1VV9vFn58AHDng/b4ZzMyn7mr97wP/hRec8TzXzKTduan4cuG\ngu65c1954Ny5FUub24OpSweXJudmV5CY+cCCnYafPXfud8+dew7pmKe7OGXfECGM+e6nbzXEqsQo\nVkUawiVvsbev8xXryafXjI7/xVqfeQPtaHpR312ZTxoCs16iakukHgThiVwHuNCkOMsUJinLsMVJ\nvK1S4JGrOqGBd5K1XN5lVZeQz4yEJhdOCrXxaYS8ZWs92kDk7ktdrGyt82Fv3o/f1odPwv4+Enhh\nsjMPwYVFx0ZvXnjZvFE+dxtMDOAiEsE1lC9uMc/E4swPNcxnAJ7esBoyqtguDxNooCWyWOzzK525\n4bW+TCju7O+Q3YDuqszIeomq0d3EAx2e0PUhtR5lRRtjfwNpPepAGsU9nynwyFWdUO2dT34h71Qa\nVoZLSDLjQ1MIpzez05A38ihb66+bUpLD35e6CXUP803oLV4yF7S3YXwDDOEK0BKfhJMDO5abj9l5\n+Mz+efMUvAa/BZYIKpESnvBMDHpvIgrfw2z2Ta0XVDlOA+xVBRpoyVsc2OdXZk+fXv3NdRhbIM70\nHO988o/1jjH4YGfnLJF6oMMTug7woVSrzOisyDxSKknRWjysAkXMNmD4A4vMR67qhGrvfDUX8k61\nbmSkhKRyfGgK4fT2NOSNLc/VeueIqnV1X7p3iUXB/Di7b9HcfJhr78QSGMLVOhJAASWHngf4z+kV\nc+lahN8DSwSVSAlPeYAKjGCOm+9+rM6cKjtrD4IJNLgFttgMyT6/Yk7jy/Vn03jnk7pD/rzgsHNI\nJB7o8ISuQ/XZNMqK6ENKJSlai4Z1IM2MKcYfWGQ+clUnVHvns1bIO+XUyAQlxJnxoSmEU+wZljfm\nPLCCpD+56WgD0V78Us1LdO0V2R8nDh76EtzXh/dehIl1S/TxxgIJ+EPP5B3q/fbBDyyMfRg6F/pm\nyRJBJXqvMjwWGA+KyxMw+5OSKq+RYQINfoksNu9lss+vmNN+U+t75mQcU37nE08/bVD4OMSUOIde\nJh7o8ISuA9wjKAkVZSVaV0mK1qJhHUgzY8BAg0eu6oRq73zWCnmnnBqZoIQ4Mz40XCs6nGLnsLwx\nZ67W/S4vz3P3+R+fMR8yh1/Bcfghf35z8wcwduZJ86YcNx//Z550F1sk+D3hHZre3NxcgF/5xrq5\nDHwMLBFUovcqwwMUWIrL2NF7ByVV3mCGCTT4JbLYPO7O+yk+dZ++sgTTazyOCdr55Oe7f7353Lr5\npv8V82KO8O6DnbNE6oGNE6tXroN93i8dOisRl05StBgNq0ARLwYsAo9c1QnNJr+Qd86pKqE4NFIr\nKpxs6LC8MSP9iKSu636XlzAZ6s7+HW78kJrGgbnZ2j9w09+m1ceJIIdobG7v6BZ8geZ8wmlo7+k9\nD28ziWHc2wolWBXLEzG7QFRyLgidSBh5gnY+8YQlOs+aXVv2mFjEE7+w6rwHrF65Hv4ozUhE5LJC\na/7MSYrm9bABSAu4EYNHrkqy2Gf2Lshahi1YDUqIhcPQMDSHU5tYyZswhtd1v3eTdnkJk6FOf2PJ\njdmUYBkfonFMtT42N7bmGDpREeBjOy6x/fAex8qvwvPL8Kt5GKl1UcXyRIgGmqFzSegwMciZNu/6\nnU+yYKmpM0+4iTgHrDrrgahXrmf3ftWyoo2RJEXzbjgCUE6eKyByVZJFPot3QfJTNghW3fcwVq0I\nh6EhaOBwahMzedMMdhTW+kG3e5d2eSlm/qJpb1S7lutl+/06Ht/y5xe+9z99T0anffPm+3U8Zo9F\nSzxknpl/+d7neDYihquCsoaCfbk9vd053Lz7J37nk7ZCvoL9S71AqvMeFNRn9/TWsqJ1AicpmnfD\nEYBy8gKuXeVkcbiz3qVsSgmXUFaYofPhbNyL7Wv9wpzR7Pdu+l1eyhR42zxM4tFb1gt21L3H/0r6\nyMNn3OrHNzdTNjdz6JB7F/R+8cr6MJ7xzc03CjwNqioaCvbl3sh+8+6v+U1e2pz9a+aBHI/H1AKr\nzntQUK+LyAPWsqJ0giQpmnfDEYBy8lwBoF2FJKF57xI2pYRLKCs8JJyQy5uCx4F9iDXHUwPzvY/f\nu2l+5ff5wyX38kXw9+vwaDB7o5GnMg591W3e7ftNXppj70Znxc30lvTKFkaTixmhalYifklStGCH\nowBlxEHAd8DVnIJtzOXylsJNH8O5B80rbQwFt8tL897xyot64qYZrbnNu7TzSfvdOfNKel3QLNsc\nVbMSYVeTNApQhIvDKnhO4Jqbc7vBuz90ho0vBAZOrASDm5pUP5F/6qqHYseysmNAVz0EO6Nw1X5w\n/tTAgamfvd91dmc0XPco6vf3zgNX3Z8dy8qOAV31EOyMwinzfH2n+bdCPJ7a8MStUxCBEX5/D6R2\njtyxrOwY0M75dnWRvrUO7/Zfo9j/KnM7icWE/Zs32SGuE4Wbd2f0p9yuBeUN0ipnmxXZ1m3neU9r\nJTsFIIFFqgKway6+E8AuGv+6DvALay4E9sYUN6O74a1XjIC9zaW98O9ESGxW/Lbu7alXzx3bg7ou\npffY3wj2/Abajjem8iPodenPLhiNt7n8y90uKKhDYlZkW3edubaqnjtqjDfq2pE169mn0T17Y+p2\nEuPw1ouLAN7mvnO1jo8Lsq1761l5p587tm75zkiO/QRxXh3Yk9276XYS29HN+j1M0oAGN+9Sre/+\n9zAHztroy+F31L7ex6kRspIHOvxfDnoEILHluqb2upt23COAezdncDO69emm/X79QVdVPq9+8y7V\nuvm/7V0+Jj6hFLgdtbzhfoSsZIH2DP7NwY8ApOy5fge+4deBeXLBbUan0YjnctdQAWrhsf9AlD+a\npO1/JhWOIvBdL/08icwuEhWcudaDuTJZaZJa9QAN/04Ol7d15xbjOa8jB7T3i7fF3C3jikMs3sKD\n/zXGEiHRJl3MXxkYKLXhPkfUHOwkDi1pos3ekCODIZwtPOYzZTkP0yYNv5+XNrMlYDBdDja8VDcn\nPVqtl5uk1j1Aw+dy+h/PTRbmSEcOaGb17YJUdbrskIi18EAxN23SxfyVgTm1Sa3LTmJxopWaqjUo\n9SA1Htm8//VCrdekuXkO7DseWTwceHYJeme91F2RNA5Hq3Xz41yhrWvOg++SQjS8s07D4CzbuoPJ\nEul1ZIG+A2dKYrX5skMiVePRLoqMUDXpSv6GA3Nq01rnncRiRitlvoguNyj1IDUedgn7Uea01qSl\n1h/9p0h2OHDXNKRe8FJLkTQOaUN6bi2dM99kr66l02Ym5wHnCw2fzMnlt3XnOEVHFugP4JcKUtXp\nskMiVuPRLoqMUDXpSv6GA3Nq01rnncRiRitV6RrKEDUedgn7UbJIQNSkpdbnirVeBB673O8uBYpi\nkndQxwv58fliW9ecB5yvxHBG/3j5HwKYh4mcDlp84skBkaOcyw4JSo1nuIs16UphDAfm1Ka1LrYH\nFLeNZEL1qBRO1TW0gYfxLAS7dGra3MOMKs21Prloan1U4NUfvU+8aKAYn4lI6OGgrateUjGyS5Qv\nNFzzDh2VozRU1DEwABEtDrXwEB6qUS6OLF0pjAZgSm1jrd83cI0Ysfuo3Rpi20zqPRoucve7BqW2\nx59rMfqCmw9fQx4LnLQ8xH6UVgN2+QtFDe2lcX+KkuaSuQtMrdsl3sNCsWLg1PTu5rORovqQ8bUN\nLERNUsv+s9vacEZoIip5aJLHJGGYsXGotbbukEsJe6+VeKcZT8JMJcm5wQrJSmNTUO1Xkj/kQekE\nWMJK6JTatlq3bSOxb58lcGuI7TeZ26PhuoYiL7YYxbaLpNSfAx6L51semkVyCftRWg25vTlOGo3Q\n0lwyc6bW2VDUmQCnpr/2obcGyNv2wviRDSxtvr62HpT9F7eV4QzQRpTz0CaPRmKYydohDiEvex8p\n8U4znoSZSpJy4ypEi4ch034l+ZOwxsBBWAmdUttW67Z/5IRtxGgJbLtie1T6PRq86c5g73ENSpEX\ne1L61rNh+9CQB4Fdy8OgmSX2o7Qa/N6cjAY0IpAOWqDOrplaZ0PNe0kaqDJwYvrkCrz6YQpOw5nx\nAxtCM6mtq/c/XCL/fadH6dKJhjfoDlmCPEQxDrkqNIfZEtbawKGgSy33qcWUsPfBlkur46jrXIs8\ngWWRi6jBVpM5QqO9NIYskM7lz4c1B+zDKsCc2rZat/0jD9hGjNxI0vWbPDmw9gbHy651J/JKN9OA\nwZAhj+tqGTVt7GA/SqMhtzfHS2Ota2m6PN4OptbZUNTtrwsMnPZzHF+EGfMZ03wwfmQDAeybd21d\nfVJo2p69B9LpURkecjbQ5Tw0CFsWDrNFstbWHXIpYe+1Eu+049GW+csv58ZViBJXIdPSSf4krDFw\nEFaPzqltq3XbP/KCbcTIjSQvIZDfo8Emz/quocgr3Ux53RKKB/HiloeuH6XRIHtzGIGksda1NJXM\na6dPv3WGDUVJHysGNpOR6QfM1DOsZTjB+JENJElNUiUptMIesNvKcGJrPJfz0A7gwoxIptbrDjle\n9l4r8U4HeBxmX5KcG1chSlyFTPuV5s8Y6o4Y2BQYfmDQumkLY0hMbVutT6+YWreNGC2BZfacOfMe\nDaTty/PUoNTySjdTXk94EC9peWguLstWWvbmMAJpQCO0NJWM4T0uhqKkj5WhHbAhjH3hYd78kG1Z\nETIFNAcissGzmJ+ysa2rtGxiWfaA3VaGM18bUc5Dm7yEGZFMCdUdcilh75USctrxaMt8SRp+lxus\npqy0C5mWTvMX1zoDm5/G1xUucGr9RpgDK3o9GmHzSduIkRtJYo/KeI8Gdw0179qPmQd66mYaomke\n14/ykGcQl0wrDyP9XmwUmZXGWtfSQclcFkNROgbuS0NFjz72LExGMQr1JjQHIrLBM0qTVE6KX2H/\npdOjMjzRVJ8o56Eux6scZspW3SGXEvaeYSxBTjsebZnUusuNq5BQnKSDBrhcX2n+OKwxcBBWj86p\n9f9uO342VJvSpm2kucN6U5qOYo/KxyNG7hrqeLmbacgW8dh+lJPznoFdsv0ojXS6N4elsdaxrSVL\nS8ncvfmMdLy02DEwQGw6/Ps9vPcrNLdIW8NnF8yy8oDYz/u2rul1XTwoGE4QbedKHpoBXJg5W1WH\nfErY+1AJOR3gcZi5JH1uXIXkpF3ItF9p/tJap6RLWAmdUnthDmd012FiUmf33YqdwhSbc3mPhvCm\n96sEKjzco49dIp7K/0iREUGHP44niTNPAlw2nWRbzozPHqRSnJRkSYQSwxPe6sS2nZEwk7ViW6xZ\neNn7mIW+PgsrJHZRsp9Ic8mIX0n+mEcucIRTNv3IGvLMuGcB4s+d04Z85T0awgvUdjGBZB7p0ddZ\njLgqe3NmjzlekU5aoHIXQEiAy6ZHFlSHGRsS/qL/NcMTlOrEtp2RMHtrA9tizcJL3sccQX9Psew9\nERdnP5q3QwqZSCf5Y54k6WXTO5e9rlfXM0r1VNKQr7JHg3iB2y5qLDsinnzLQ+Qv783h/SkVaeZJ\ndVdMT5lLM4xfsaHsf0WopLAwv21nOMxkbcU25mXvU6uYp2IZZT+VJiOgIs08iXjZ9BP0+8nklX4i\ndmviVgRunAjMPsAVfsfKjePWLU9uRSCJwAtLdur/AU2Q+ZgFSVKXAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$\\left \\{ T_{0}^{4} : \\frac{T_{e}^{4} \\left(- 2 \\epsilon_{0} \\epsilon_{1} \\epsilon_{2} + 2 \\epsilon_{0} \\epsilon_{1} + 2 \\epsilon_{0} \\epsilon_{2} + 3 \\epsilon_{1} \\epsilon_{2} - 2 \\epsilon_{1} - 2 \\epsilon_{2} - 4\\right)}{\\epsilon_{0} \\epsilon_{1} \\epsilon_{2} - 2 \\epsilon_{0} \\epsilon_{1} - 2 \\epsilon_{0} \\epsilon_{2} + 4 \\epsilon_{0} - 2 \\epsilon_{1} \\epsilon_{2} + 4 \\epsilon_{1} + 4 \\epsilon_{2} - 8}, \\quad T_{1}^{4} : \\frac{T_{e}^{4} \\left(- \\epsilon_{1} \\epsilon_{2} + \\epsilon_{2} + 2\\right)}{\\epsilon_{1} \\epsilon_{2} - 2 \\epsilon_{1} - 2 \\epsilon_{2} + 4}, \\quad T_{2}^{4} : - \\frac{T_{e}^{4}}{\\epsilon_{2} - 2}, \\quad T_{s}^{4} : \\frac{2 T_{e}^{4} \\left(- \\epsilon_{0} \\epsilon_{1} \\epsilon_{2} + \\epsilon_{0} \\epsilon_{1} + \\epsilon_{0} \\epsilon_{2} + \\epsilon_{1} \\epsilon_{2} - 4\\right)}{\\epsilon_{0} \\epsilon_{1} \\epsilon_{2} - 2 \\epsilon_{0} \\epsilon_{1} - 2 \\epsilon_{0} \\epsilon_{2} + 4 \\epsilon_{0} - 2 \\epsilon_{1} \\epsilon_{2} + 4 \\epsilon_{1} + 4 \\epsilon_{2} - 8}\\right \\}$$"
      ],
      "text/plain": [
       "⎧       4                                                                     \n",
       "⎪  4  Tₑ ⋅(-2⋅ε₀⋅ε₁⋅ε₂ + 2⋅ε₀⋅ε₁ + 2⋅ε₀⋅ε₂ + 3⋅ε₁⋅ε₂ - 2⋅ε₁ - 2⋅ε₂ - 4)    4  \n",
       "⎨T₀ : ─────────────────────────────────────────────────────────────────, T₁ : \n",
       "⎪      ε₀⋅ε₁⋅ε₂ - 2⋅ε₀⋅ε₁ - 2⋅ε₀⋅ε₂ + 4⋅ε₀ - 2⋅ε₁⋅ε₂ + 4⋅ε₁ + 4⋅ε₂ - 8        \n",
       "⎩                                                                             \n",
       "\n",
       "   4                             4                       4                    \n",
       " Tₑ ⋅(-ε₁⋅ε₂ + ε₂ + 2)     4  -Tₑ        4           2⋅Tₑ ⋅(-ε₀⋅ε₁⋅ε₂ + ε₀⋅ε₁ \n",
       "───────────────────────, T₂ : ──────, T_s : ──────────────────────────────────\n",
       "ε₁⋅ε₂ - 2⋅ε₁ - 2⋅ε₂ + 4       ε₂ - 2        ε₀⋅ε₁⋅ε₂ - 2⋅ε₀⋅ε₁ - 2⋅ε₀⋅ε₂ + 4⋅ε\n",
       "                                                                              \n",
       "\n",
       "                             ⎫\n",
       "+ ε₀⋅ε₂ + ε₁⋅ε₂ - 4)         ⎪\n",
       "─────────────────────────────⎬\n",
       "₀ - 2⋅ε₁⋅ε₂ + 4⋅ε₁ + 4⋅ε₂ - 8⎪\n",
       "                             ⎭"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "out1 = sympy.solve([R_s + sigma*T_e**4, R_0, R_1, R_2],\n",
    "            [T_s**4, T_0**4, T_1**4, T_2**4])\n",
    "out1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABoIAAAA0BAMAAABcCzflAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEImZRO/dIma7q80y\ndlRAyO8oAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAVz0lEQVR4Ae1db4xdx1U/7+3ue5v9+2qJtEqI\nvLJb2gDFm6QFEkXkqZX5gCrFQWmRqUoW2gTaKsqKSJQKJK9aoJGqyBuBhFoQ2QpoK1eW7Q+p2hjw\nikgIpR+8TaiaL8ZLqhgRkSaOkzh20i5nzsw58+fOzJ371rte9HylfffMzPn7mzPv3fve3LMAzY59\nG0NzNANGuLcSoVfESiOiXTRnAyoPPNnK8IvC2D6mIPLC5kIh3/Cy7UCEOts4Gzsw/G2Iftd8qZHu\ncinnsPLtRIQ+v32TsRPD36rou3/Fmkc/tshk3XlmrY5j2Md3IkIvbd+k7MTwtyz6w48Z1Yfnim2M\nrIas3du+UukLeaT98O1/K3SCKNHX3XPzekJcutOmshamij+OxZJPVBHyxxu0so420AN3NmHeHO8V\nDH9zjgwg/VRTme4DPRLpXiqXrF5R3wDt10rlWyuwr1/DXKLvYZh5s0ZNxlTWQme5RnHdcBWhOonk\neNbRpFRkYCHSh10v6PmPDw7YewXDH9CDwcU+3Vj00AqJTJwvl7yxwvo1gI9WOsMOcxkx0YOxA+EY\ntw1PTh9fIpxdhzdYLHHOmMpZgB9sdgVVEUp4SN3x66sCKHJKw7HW8bBHtVt/f2wLVlCz8GN+Rfvi\nODEr5wW3g3NeWJinDwpZSoy8TZyzUYDjSp6vdH8M4NhqpTfoMDGMLcCEthmMq6bhyeljpPb3W29F\nNLhdGVM5C90vb3YFVRFy3Qrp+OQWQBEqyrVn5uKjZ7dgBTULP+5XpDeOEzNyXnA7OOeFhfmPvytk\nKTGqL4V2r5QKALy7wnqqV76CJs/Xr6CcPotU7VVcxlTOQntysyuoilAFMqcjPrmmN+eoo6OWbPfj\nLGe3YAU1Cz/uV6Q3jhMz2rzgHu+cFxbWueYrqHuBpJusoBNizyEe6sGuI0u6o0qofhtD5zy0fv6b\nmlcIalqeqD7FY5G6fg4iNkmPvLimpFMTrgVv6EubXkFRhNBGHpgMFJ5/AzUmEytlK1ZQKvw6x6vw\nSI8S5eQQnLxRyQvuFTYyy8I8aufC9Wp8sfkKal0mDU1W0HHXpqFHX4XRewGm3z8HRMD3cED1tG+/\nxbBwDAB3z8H+PsAzz+CIIqZEo/CE+qwaQeoPP9lLmTIW8WRMWQsyZCxYxTzSWtv0CpJ4WKc+Kzwy\nwCgo9v4lZ7oLha9moNZkQmorVlAi/IQH0l2aQJw7Aidp4LwQLVE4RUYR1bm/CbZjBY3OS8yW6KzA\noXWAL8CLQMRTX8cx1fNe+JbhkpSAz0D3pwDjczMrREzdsWBY5G0GQn1WDSOFANybMsXqwJhyLMiQ\nsWAV80gbNruCogiheoVHGhiFSXtxgtNP4EJHr8Dx4YSOLVhBqfATHkh3YQJJ7gicpIHzgrXE4RQZ\nRVTnfm5bVtA02jbHjQ+o45+wdRt+l3Dr++BTMIu3REjAOexUxO/D6UUkf+HkyftPnnyS5EbmYeLC\nnqVOH/CWCAmYWaB+lyfUZ9S0T5789sdPnlwhfvhoL2bKOgVgTLEFd8hYMIpb31aRfHYVFf/vpleQ\ngxA5ymZj3tqgCZPFqYtKxPYStKRlcy//mRDfghUUhp+wLN0uPAUJJLkjcOI7j80L1au0xOEUGUWY\nuWf7eIewui0r6LqexC7ENH4uHXxHr/U2jK0rQq8gRXwKDvU1m7ypPgsw9tg7YPcq/A4RnN/IxjwV\nfVYNv9d8DuDx9YQp8cuYcizwkLHg+GdGWnObXkExhJT2hLcmaIVJp999zfjhQmG6NnN6LCE82ApS\nER6i7ZsRtanwI6xeF8JTkkCSOwInaeG8YC1xOEVGETapjB+7YFtWUMcLWzd+BN3FC3hDcxE6S4CE\nXkFEwMtmxXFKzMzD52ZxxeFn0x+AIpz8Zp6qPlHDSG30cAUlTGmnULEx5VjgIWMBm+yfGZk5cuTY\nb68x20DnGEJKUcJbEzRBATPnjUkXioG88IUW/Ka0BltBf45z+g+iIyBS4QdslSbCU5JAkjsCJ2ni\nvGAtcThFhohw7l88cuT1oxW/ajqafZOg3l6ermrEbTDTi3iJNnUeOmtA12rnkEtftD1o+DklngL4\n4eQKrqAePAiKcPLb8ET0Ia8+GKmDAA8pm3qx+qYMK7Apx4IZYgvYZMUshJ+Qy5ZuSiUQIjUJb03Q\nBIXCTx8uFE2dCPjx82dqJejj5kAraPTfetA9zirccy58ly9GIzwlCSS5I3CSMs4L1hKHU2SIiMw9\nplXDo9kKUu89z1Ut/NGevV+H/T342YswtqaIHt0HEdGeN/wmJbq/t+eOhdF7oIWfQZ8ARTj5bXgi\n+kQNI/UIzPwkZcpYFFOOBTPEFnogim1Ys8uWbkolECI1eWAICtjLBl0ouG+w8/g9+CXpUkJ2oBX0\nrk4PRtZjGnPhx/jdPoSnJIEkdwROUsJ5wVricIoMEZG5v+R6VEQ3WkHqvQe+UdV7amPjDRg9+ihe\nxs72FaE/GIj4IrOblJjEy+cF+I3n1vA+CK/5kXDy2/BE9IkaRmr0ttuTpoxJMeVYMENsAUAUs5/Q\nfvzykjQaEimESE0eGIJinN9v+JaQHG3oRMD+kd8FmFgMOk3zjvueXIqP5HrXcAVdH2PIhh8TcPsK\nE0hyR+AkJZIXnIYqsypwigwR1bm/ecP8Tuk6lqfdFXTo4JGNr37rVVeg+/pXj106cmpO96n3Hlhw\nxwOavkqhPnUVp47RudFVIjgldANf6bs4as0scKdZQdykL1WoYdUwUsJDH3eqZXnsmKGshXAoIxSy\nYnvzCEW8dYP+EHzEmHV7I55Ilz9Duttzs7WEK2jAe5OYcrzW2vOeX4RfA9o27nPUJoi4nSAKE8jm\nDueZ80s72LSJw8kyzeY+4TG4K2gPAO4X+1WXtT0HYwswsqT71HtPd9kdD+hn1e9BdLxgzt/78f/g\nqlPHO/VJXqfn8fcgOmYOcGfII/qsmtYiM/O5aopH5GwtSJchrOJwJNbePEJQ9dYJeupffsxfmjm9\nMU+kz58h3e25OQ24gj4v/I2ImHKlYKzXWgDaNu5z1CZInXWZ8DxONncYTnDzgrUk4GSZZnOfct1d\nQfgDJ35F8Wcu6wS+787DeJ/66L1nfM0dD+j2rbfonocfOqqJ+zY2Ah5p7t2r11b3ly+ndIq+jJoC\nUxkLGcXiqCU2j1De27GNjVestSLKmyEj4bn5r3t+cxUeLdJVYYopR6bWvl9pLwHtb/c46hOkYiHo\nkAnP4wScO8Lm6WEtcThFptncexachllBp+fww3kJxn8K8H1nFP4d4PA6TOtUV+89MLLqjg8XvRMR\nCmZITUjoJn6TcOdgExVRbhThbRDtbw84hjFB9JcPh/sIPN5wngdAyh6ruLMAt5+ZDnzvGfSK2qr8\nf0ztRIT8GSJwAzd3XViEhcFQjyg3ivAb2VNqR77PMZQJ8gRB8oAGBu955PiMofCtxj226PEp18RO\npnciQs4MdVY0eK6b2BN/vq4I56ry0T4KLpDwQ+rixOFQnUOXIC8vYdTtiyp2uufR+6sR9I/rLvP4\ng27g6/NCDSWBd4V6f/XOQYhmyGzq/pnjelbITTtBqefrLEeKcqffKP8B2pheUwK4vx3o3tmVHroE\nmVCL5z19jcHhdbO/WjfVq7o1MsefqvO7uTWcZ4UQ7a+24V9thJR9Z1M3Oabc1Ed3Fc/tPr4MdDjB\nGfnuTfgzXoca9InncAxpgvxjH27Er1XowHse+2W67nIqKPyz6jmhu4f1VSEkv+loEK42Qsp+x2zq\n5mkhN6kxqT4t1PN1zvZz5io4O8Exdwtvnb9AjdvUq8MxpAky8wm4/haDDl7S6u2xpo0ne0U9c1n1\nHlcvG0NyqFj9gy76z3l9MYRmtwwgzzQ1lH1nUzf12XuTl1ewY5I6zUvGt1dcPqJtcHborkX1ExNe\nys2rV8thEyRjYsuAuRqKNV5TDwJ8cFVhQZe0enssNemFrqh3nVkDeNf92DFKqNnhIaP0bcE5L+qr\njZC+55FN3co3crP1/jP4hvgXB7D9Yc/hJg1S3j6z7socmh+fU221v13fO8NQJ8gH+/gk2m8pROie\nR2+PpSa9qCvq1onWCsD6XT28osbmEB/6ot9fQVcbIX3PI5u61eyQmzf2bsB97yN4jQe/NPCUkfIP\nTC27CibOq91dQPvb6ffC4U6Q9tsKnM8SQtch2vgZNLaGLf4u7nEEa+LoI2pcoXmdwm54jt3H/VgV\nQnwfdJUQqhT5UjOEFRbUC5ivy8jNI88tYQ8tpsfUWMkRV/6V/yIT/EXf1Ks/Urpofzs+54hjw5sg\niMPdawqNs318ufnUW0fVfdCsamCJATzaf7Px5BqM9RUNh+b4SxhqDsNLhy74JVJCiFfQ1ULoAUpn\n9knPEO9Cpm/H9ETCvcTSVd84LzB37Tmq/OmNVSWolSPxxAnVPKV25A97giAOjyswQO3qMYfdHss9\nMLuKT/wAqLKMT0vncBD/HQvTv4pTHNuJ0MxixCe7C1kG38ZHr/B4C68rVqSzjogpn1mcWPPkXp73\nmtsbfmD66jcRYDx2W0x4Y6vj2nXr6j4IJl+JPl/nMIZkuvS7cJaVVsdr7vhRYIGu2OPSmaGbnvmA\nkpmLCfLmXjvWFKGM25khdeuB1wXWrFDOLmTp+zvcPY3HqczzdcKLhLYbU95ZDTZDTi66gkg3DV/E\nc8EaprIEeUpUBkSBhUwW4I9/gT5pmgTxnm7Qg7yxVVjxnujomR421S3TN5zuWhLX3b5+DVdZafXO\nclxNiQW8/EhIo87kEBYUWsfI1yJ2ZXOvHWuIUMbtzJDx9iZrVqjYLuQbznxfje9bTT5fJ+JIGLsx\n5eOP/IfLGaEbhi8assEarrIE+bTo9IkSC5ksAEgpBp0gCBz9yAO7V3zD0ZYqs7sQHUl0Zkq/i8TX\nMlXr7dN0qWLwWQsvsZGKdL3imSWqBDDOKsrOpQhl3I4MBd4ulfnCXLhtq8N05mzsLmVYaodKwxdF\nkWBljImiBElWjM9ZCHBlg+osQ0nFJkGarSBcbdnn61wXiB5byNTINtz429+x1Yqk7pA4ksXgsxZ4\nBVWl6xW3l6G1kHAr3V2KUMbtyFC9t2mP8AJ9pej5uojdnNboWGn4IlxitChBkhXjcxYyuMpQUrEk\nSJPPINyFm32+ToBhIlP6nVn0Hnlp+YTEkSwGn7XAK6gqXa949FKvveR7U9IqRCjjdmSo3tuca5hD\nj+bGzVjEboFUwFIYvkiVGD2V+bcGgszcd0WnT+QsiHQmQZKKJUEaraBjPb6llLLe+UreeIWZrDJv\nI8U98nE1EuKXJvFORnispKLSFngFVaULFB97832+naKWRUi8FcJX4Lrtj3gRqSHf24C5ton71u4M\nmNIuBYxNm7nw00YlmYQI7NYmCFWMT0m7SHs8Pq6efzyUUcwJ0mgFne13TGz7+7pyfE0lbyn9Hqky\nLyjVV5mnYvBSUlwENXG3rmM/dRybyq1K8XaRpqrvJMQAyVC1CHl7g3700DaKXwWhUmDI7UC9iUgC\n8bwNeOub06+FN6/xqus0V/Xqshxe+Ll/Q2C1mGA3lyBUMZ5z0urWVDJBPFxlwkiIh0SxTIZo5wRp\ntIJOr5nHp6rVvw+txyp5c+n3WJV5caWzQjXZpWq9LQjOcVAxeGWBeERQE7qOPRWXj1cbF2n7VEJF\nsbXJ2l/89df7TJefBSHlbaTKvNXkuG07NaWH7AMLnrchc2176ifh83XpuapVVsPghS8TWoDDJhNE\nVYyXnAx9dJD2E8TDVSaMxHmIFdvJEO2cII1W0Ozy81pBtfp3tJK3lH7vhFXmbbXvatX6apV5KgbP\nJcXRA1d6RNexp5Jw8WrjIm2eSnDKlMtQpQj5+AqcvUfQKiYEoXyJcwGGK9lFIjIPLITeFrsijJfC\n5+vScyUyAxJe+PTTfCEOYYJ4j2LgQxQy+Y7vFpkZVTFechJ9j8AZJoiV1v9rQBzFX22kjr1VbJ4e\nsYolQRqtoMkD9Csj7mjnyvH5St4Azzq8SogzxpkhtUeey4VTcXZbENy8E7SoGLzwOLJIGgukOFpt\nXKRlN47cWdgha9NoH1uEKftsoemsPwlChcBE8OCI7AMLHgz1PoQcT4z0/S7lWyVgQtLna95yw89W\nkRfVVyRBqGK85KTo1kQ6QTxcZcJIygyJYjsZrF0SpNEKGrmobjbwUGW9T6vK8RdU6xz+ERFW8pbS\n78SrhCIZU19lXheDRwvhkxdKn1tcniyExdtFurqCnKHQ9d2o+5v41/AQhAqBieAhEUkgZjoHrYl/\n158EQaTnKmBs3HTDz1aRZ80c7OYShCrGS06ybn1mC4S0nyAerjJhJGaGRDEm2nlfryRIoxU0/qZZ\nQZMrpnL8k0qtWkFEVCp545aIHwqvIqoZgxtVwqr1osbEgWJYDB4thE9e4ABbIMVkoVK8naWrK0gU\nI/Eg/jkHvsXAXzvtQlIQKgSmioeNSALxYCh0xGE7/GWnocj0XAWMjZtu+Nkq8qyZp+905N8QGJ6y\nBDlI/8SA1LBufWYLkQTxcJUJIzE7pBXjV75rvl5JkJZ+Q9q94o/HW9031McIHtXq3/tVufm2GdZM\n+KBXrsq8YeIa8E7RcVFj48Bi8GhB87B2dRYLBFC82jh+ZJK0LodP0qHiquujJxr+9qW9EoQID/XW\nkgemuoJsRFKF3vNW22nyeuiTAXfapYCxcdMN35nQehzwM6jybwiM9bIEuWRz0vPawjmzoPOWHwTB\nd3NZFJkE0YqtEGuXBDG1WMaO80j2vLFuhrlyfL6St5R+3x2pMm800R55VONUrf8i+yAhqmLwwsOj\n6iwWdCpGq41TKXlxlKRDxZEC9E/fau75XHP1NCMk9oiQiIyCwG1XrR0a5/cjz1uXt4yeDe/nMi6V\naUxzOeE7E6pWUMYoJ1PwbwiMlaIEURXjRY3rnoVTrSC/ur+Hq0wYScuQUcxPjziqOUGm1OMjuLk2\nuMpzOF1SP8yle+jrEyIVQOpIV/K2vBSHZg9e5akKq0biYFbh4Q57FsV2n/9LdlRT7Kjz7mNYrM1Q\npmHbRUgqkmS0i9tVOxJIBYYqb64nOrcZl3K66sbc8O1kMe5po1coQayaqqOCdAZXdrSaICJUUawf\nUYUpfTFXGQ46fs5pV6t/pyt5W96ZA44Kj5SnKqya1qLHob4uejHokSYrdvb5v1MGDSFPJVQUW5uh\nTMO2i1CkenpFG7tdGQAbSMXbKnOup70QGb1iAfu63fDtZDHuaaNXKEGsGt8t1WKkc7iyo14deyVs\nhVTLO3Yv6+bZNa870fiQ21+p/n3flleZjzx5YTyS4vKxff6GJ/JUAseTcZ1Zys4uQmIvrV3crmrP\nBFJlzvWMzkVG0y5FmMu73PBlsgpwkFLyGUBEX8Z3zsmqx6I4g6s4WhFPC3Xv72nu8cuGqAhf67iG\nwDUEkgi8cIKHblhh6tr5GgLXEChEYOo7ivH/ALncszBaCOKuAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$\\left \\{ T_{0} : T_{e} \\sqrt[4]{\\frac{- 2 \\epsilon_{0} \\epsilon_{1} \\epsilon_{2} + 2 \\epsilon_{0} \\epsilon_{1} + 2 \\epsilon_{0} \\epsilon_{2} + 3 \\epsilon_{1} \\epsilon_{2} - 2 \\epsilon_{1} - 2 \\epsilon_{2} - 4}{\\epsilon_{0} \\epsilon_{1} \\epsilon_{2} - 2 \\epsilon_{0} \\epsilon_{1} - 2 \\epsilon_{0} \\epsilon_{2} + 4 \\epsilon_{0} - 2 \\epsilon_{1} \\epsilon_{2} + 4 \\epsilon_{1} + 4 \\epsilon_{2} - 8}}, \\quad T_{1} : T_{e} \\sqrt[4]{\\frac{- \\epsilon_{1} \\epsilon_{2} + \\epsilon_{2} + 2}{\\epsilon_{1} \\epsilon_{2} - 2 \\epsilon_{1} - 2 \\epsilon_{2} + 4}}, \\quad T_{2} : T_{e} \\sqrt[4]{- \\frac{1}{\\epsilon_{2} - 2}}, \\quad T_{s} : \\sqrt[4]{2} T_{e} \\sqrt[4]{\\frac{- \\epsilon_{0} \\epsilon_{1} \\epsilon_{2} + \\epsilon_{0} \\epsilon_{1} + \\epsilon_{0} \\epsilon_{2} + \\epsilon_{1} \\epsilon_{2} - 4}{\\epsilon_{0} \\epsilon_{1} \\epsilon_{2} - 2 \\epsilon_{0} \\epsilon_{1} - 2 \\epsilon_{0} \\epsilon_{2} + 4 \\epsilon_{0} - 2 \\epsilon_{1} \\epsilon_{2} + 4 \\epsilon_{1} + 4 \\epsilon_{2} - 8}}\\right \\}$$"
      ],
      "text/plain": [
       "⎧           _________________________________________________________________ \n",
       "⎪          ╱   -2⋅ε₀⋅ε₁⋅ε₂ + 2⋅ε₀⋅ε₁ + 2⋅ε₀⋅ε₂ + 3⋅ε₁⋅ε₂ - 2⋅ε₁ - 2⋅ε₂ - 4    \n",
       "⎨T₀: Tₑ⋅4 ╱  ─────────────────────────────────────────────────────────────── ,\n",
       "⎪       ╲╱   ε₀⋅ε₁⋅ε₂ - 2⋅ε₀⋅ε₁ - 2⋅ε₀⋅ε₂ + 4⋅ε₀ - 2⋅ε₁⋅ε₂ + 4⋅ε₁ + 4⋅ε₂ - 8  \n",
       "⎩                                                                             \n",
       "\n",
       "            _________________________             ________                    \n",
       "           ╱     -ε₁⋅ε₂ + ε₂ + 2                 ╱  -1           4 ___       ╱\n",
       " T₁: Tₑ⋅4 ╱  ─────────────────────── , T₂: Tₑ⋅4 ╱  ────── , T_s: ╲╱ 2 ⋅Tₑ⋅4 ╱ \n",
       "        ╲╱   ε₁⋅ε₂ - 2⋅ε₁ - 2⋅ε₂ + 4          ╲╱   ε₂ - 2                 ╲╱  \n",
       "                                                                              \n",
       "\n",
       "_________________________________________________________________⎫\n",
       "              -ε₀⋅ε₁⋅ε₂ + ε₀⋅ε₁ + ε₀⋅ε₂ + ε₁⋅ε₂ - 4              ⎪\n",
       " ─────────────────────────────────────────────────────────────── ⎬\n",
       " ε₀⋅ε₁⋅ε₂ - 2⋅ε₀⋅ε₁ - 2⋅ε₀⋅ε₂ + 4⋅ε₀ - 2⋅ε₁⋅ε₂ + 4⋅ε₁ + 4⋅ε₂ - 8 ⎪\n",
       "                                                                 ⎭"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "quarter = sympy.Rational(1,4)\n",
    "out2 = {}\n",
    "for var4, formula in out1.items():\n",
    "    var = (var4)**quarter\n",
    "    out2[var] = sympy.simplify(formula**quarter)\n",
    "out2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now wrap these analytical radiative equilibrium solutions in callable functions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def Ts(Te, e0, e1, e2):\n",
    "    return out2[T_s].subs([(T_e, Te), (epsilon_0, e0), (epsilon_1, e1), (epsilon_2, e2)])\n",
    "def T0(Te, e0, e1, e2):\n",
    "    return out2[T_0].subs([(T_e, Te), (epsilon_0, e0), (epsilon_1, e1), (epsilon_2, e2)])\n",
    "def T1(Te, e0, e1, e2):\n",
    "    return out2[T_1].subs([(T_e, Te), (epsilon_0, e0), (epsilon_1, e1), (epsilon_2, e2)])\n",
    "def T2(Te, e0, e1, e2):\n",
    "    return out2[T_2].subs([(T_e, Te), (epsilon_0, e0), (epsilon_1, e1), (epsilon_2, e2)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compare numerical and analytical solutions for radiative equilibrium\n",
    "\n",
    "Define a function that takes a `climlab.GreyRadiationModel` object (which should be first integrated out to equilibrium), and compares the numerical solution to our analytical solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "tol = 0.01\n",
    "\n",
    "def test_3level(col):\n",
    "    (e0, e1, e2)= col.subprocess['LW'].absorptivity\n",
    "    ASR = (1-col.param['albedo_sfc'])*col.param['Q']\n",
    "    Te = (ASR/const.sigma)**0.25\n",
    "    print('Surface:')\n",
    "    num = col.Ts\n",
    "    anal = Ts(Te,e0,e1,e2)\n",
    "    print('  Numerical: %.2f   Analytical: %.2f    Same:' %(num, anal) , abs(num - anal)<tol)\n",
    "    print('Level 0')\n",
    "    num = col.Tatm[0]\n",
    "    anal = T0(Te,e0,e1,e2)\n",
    "    print('  Numerical: %.2f   Analytical: %.2f    Same:' %(num, anal) , abs(num - anal)<tol)\n",
    "    print('Level 1')\n",
    "    num = col.Tatm[1]\n",
    "    anal = T1(Te,e0,e1,e2)\n",
    "    print('  Numerical: %.2f   Analytical: %.2f    Same:' %(num, anal) , abs(num - anal)<tol)\n",
    "    print('Level 2')\n",
    "    num = col.Tatm[2]\n",
    "    anal = T2(Te,e0,e1,e2)\n",
    "    print('  Numerical: %.2f   Analytical: %.2f    Same:' %(num, anal) , abs(num - anal)<tol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Integrating for 3652 steps, 3652.422 days, or 10.0 years.\n",
      "Total elapsed time is 9.99884460229 years.\n",
      "Surface:\n",
      "  Numerical: 287.85   Analytical: 287.85    Same: [ True]\n",
      "Level 0\n",
      "  Numerical: 224.73   Analytical: 256.29    Same: False\n",
      "Level 1\n",
      "  Numerical: 242.05   Analytical: 242.05    Same: True\n",
      "Level 2\n",
      "  Numerical: 256.29   Analytical: 224.73    Same: False\n"
     ]
    }
   ],
   "source": [
    "col = climlab.GreyRadiationModel( num_lev=3 )\n",
    "col.integrate_years(10.)\n",
    "test_3level(col)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Integrating for 3652 steps, 3652.422 days, or 10.0 years.\n",
      "Total elapsed time is 19.9976892046 years.\n",
      "Surface:\n",
      "  Numerical: 291.71   Analytical: 291.71    Same: [ True]\n",
      "Level 0\n",
      "  Numerical: 223.20   Analytical: 262.65    Same: False\n",
      "Level 1\n",
      "  Numerical: 247.60   Analytical: 242.93    Same: False\n",
      "Level 2\n",
      "  Numerical: 264.53   Analytical: 220.04    Same: False\n"
     ]
    }
   ],
   "source": [
    "e0 = 0.3\n",
    "e1 = 0.6\n",
    "e2 = 0.2\n",
    "col.subprocess['LW'].absorptivity = np.array([e0,e1,e2])\n",
    "col.integrate_years(10.)\n",
    "test_3level(col)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Conclusion: The three-level model works"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAwQAAAA0BAMAAADGR+LhAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEImZRO/dIma7q80y\ndlRAyO8oAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAMPklEQVR4Ae1bXWhcxxU+d7XalaXVai1oE6wG\nLXJSmraulKTpTwhkafBLKdgJbYob0m5/bFqnIYJA29CCltKfQDCSoS9NS7ylNAkqQvJDQmslaGmg\ntOmDNzElfnEsUpxCaIgVJ7GluFHP/P/cmdl7d++qL5mHO2fOOfOdM3vunbl35ixAujK93b9yMZ0r\nQjuXyKO+gifyIK4kRpCurqdT3wntQj+N9BVcc3x8RmsEyeJCUPx/ET7cT6t9BS/+Sriev2dWkJ3q\ncruTxs7L3+inyb6Cw/xx7vt8NfEgBlqJVTsqPvT533bUSaJwexKlbnV6By8FppjikQp1rLiZ3L8M\np8aoCdO15Jb9mnW36AIbnVuYmOsBT9wfoLAQUJ5rUuHwRkDHEk1Y7XDTPW3xZ3u4AoMHw/0TSaMV\nl1r0++UsQuAGdxm0eGro/wyFYOAq7TfmHIMFyZuvutkervJDV+AhGKzDMHNAF6any1V3n/NZhMAH\n7japceXQi4+GQpC/QvtMNrWuHcjrO8hNsfTDYPMQjGxkE4JczUCXjUxC4AOXVnyEHHpuJBSC4iUK\nkCYEJz0mxxcbTCIJ0hR+RB9/iknplYcA6UKKKVADMMkRz92eSQh84JoLcnCSIEIxdPhlMATRFkVK\nE4IVzbZG5r8MMLqvCkAI+IeQCD/21wBKsqsKwQHs0XMZ8SBkEwIPuMYmg3vpJWQYoxRDj9oZhyA/\no9nWyLl1gJ/C6wCEeP4JIeF+FN/HCNxWF1wVgqOC1Ut9p6dzJiHwgSubZHBD1XITzFGKEOQg4xCM\nrivbhJo4QsqzsHzzx+A7MFahBLxGlXKrq8/cu7raBBi+NNWAcp1yP7G6enh19RSlBzwBpcLEl1c8\nmpmEwAcOED1Dhn5/iwyuUANc2NQo1dDhP1mHYFfFPd5DuyvRVRhcBySAhwA1+a0weHw30nXRVT4F\nZwSnp/q4p3cmIfCBK5tkcJMt+DqYo+RDj6pZh6CgTBsUruuly1BoAF3g2VOACtyPMXK3x0NQnoEf\nGijdNeqebpmEwAeubJLBnZ6F74I5Sj708uLi8tfaStum0i3H5Al4wYbgbZxYShtQaAOdYewQjDRR\nLR6C5wFe9uAlZeNNWiLgrtJzCELgyiAZ3OkKPADmKHkIUG9wQSnHqHQh+DmG4GwMgzH2V+Ajl2Gw\nDUhUYhNR/i6I4iEofmtKLdEe3A7sobvwRazhUeo1BEFwZZMMDp+C+8AcpQrBWGYhyP8VQ/CkMm1Q\n+aVjuBaM1QCJ+FoAXzzbjodgBI8w6gZK6saXvolr/ay7223fONVwSxJyQ+A6BA4O14K3wRylDEHu\nxFbAD/0pmDu0uP3YH9/SsYvvPLa8ubhWZbxrCxiCui63aPpGRHn2RESZ8YnI6h9rmvaZ2HAzamAI\nfMtTDM5iGEg9gtM3IoohRylDYJm1mnoIpgDuAficrpGrwmAdBhqM18YQBA9szpDvAlou8Boi7Q4t\nHxTcawTRoTbtM2XDzVHAEDzcAcUnNpB6BB+dwe8CWuQo9aEzkfOqhwA/LfBt5me63jB+ac3AUI3y\nSlM3fBKG2rrconM338Q4Dz24ZImwWfzMVqhzvANOMbp9rmC4+Zepu1uAE19XxUBiCAYrFfjevRUK\nkXqUPASnq/he0YCh9wFe1EfzN4D5dRhl4LiyVyDLAxvdkpu27BMl201cjrs8U7GRMgV3j8fJ3aTc\n+RqODW+6DQCkVGkBLFcgzxnR9Ge7nncVZhrKtE97Wm6OX5oNLk8BaxYS0bRYPYAH7NqipynjCGPj\nvC/LUU7h8qCXdAc2es8uac1+ockwdDeR0/WZCj7WdYaI1+zBJXQH4s0GKuQuMy2c99luJ47rXsbi\nu9msgddXJbUzBLXPt1g/tMJsEje10vWZCl1oBFD24AK5Uz1Mfv0bakwN532228ma5EqWB15+Qurr\nRWuHamJf22KlVombrBRbWOdqeOmuKCTZX7F6BZeQnYg/1GDid1wJ532226k6aafKzxHuSSXaEYra\nl6/azCRxk5WRNtbkTEVsWgpBwlohyQ6KJcGlrE9E+T748E0cG+ddttupbKnZsrxFuCvkEs/Ly4hz\nkaAbhdq3QqCWhzebqGwc2Iz5PYmDk+8guyiWAg+A+s0lklDrpQcAbm1Rksy7bLdT+UXn3fFzbYBr\nDyM3P6NEO0JR+2YI6PIQ7TuH98svDqITd3btCFvo9r2iAWQHroF2IG+tAQx8lSrReZfudqo+ZGqM\nTkZNgPU7KjjvYnNHC52azRDQ5Wmisge3XgfwLRo+1bVDFGmg9ncNIDtwDTRM5mgSyf1UaRcOCJ8C\n3O3En52/EZ3A33146REiJ7+G78CGyLMo0XsWCrEv9vf4SwtxExbPNvBKf7DjSCQrkyumHkXa9fhu\nwu0Z3IRO0TpAfm84X8PLjWvvLZG1YIw0gH4X5H6zfaoNg5QBc9Xud8QIYqLC0/u4LrMvQsBe3amb\nQHIEcM+DbCrWCZWoFHBDSSsMqbRM78KewTXgdOQJqk42KHhRu52CA2MtPAEAIPleL0hmn4jyrAPY\nnIiowlU8f8CCD433wMYB9G8nb8nmdgduoyRtswd/ckbqq91Oydq1TtYCGLnoP7CRuoQIJOsGRFAi\nTuQMJN6Qm49K+DjgWgCwFjiwUcrSparO4/T343tMKcGLt/y65QAmrIAI80xoJ32nlKHI3U7WJNdo\n6VwFK7JsPKm4XgrjNV1zSwMinOQWsNN1jo6uzcc9514kmtMt74GNDsTtRm2dyelHjtVsbjpwvBly\neFzjLAERwPdol3gInEiUWb6SbN4NJOu6RG8IizT5tSFayWrcXykk0OR2hxKoaioJwQG/bL+iddNJ\nh0ie5IweopppQoC6wQMbYRq/pnzJui6RCEE4+VWAW/VkM9GBjcuuheRoJgQnH3jLLUd/ZDlEMgQ/\n+jPtkyYEiBc8sBFOBJJ1XSIRgnDyqwC3avxtj1ksV9Nl16Vn8hKCwxrmr7XMrqLlEMkQVLsIwXJF\nHNjI9FUjhVfYxWl9A6SK4lIqJhIhCCe/WiiiiVtItwua136XLMWOTQ1cgPpG9WAFhIoNa4tECIZm\nuwjB+VqB4++vsTxWM4VXGT9QpRmuNNVVcSnFRXz/GVk8BB2SXy0U0Rx92/4sCLgkOiWtFbjMWSYD\nV65LoPxbNJuZpjVLJiO4KJ7ofB10EYLT7QkGK/NYzRReZfwozXClqa6KySgm0vafeQg6JL/aMLxd\n+q99YON3yQPhZytwAkpyls28XdUTv+uEimIyiotiic5Q7SYEYwv8wEbmsZJcXp62pW8XY7KuVEFH\nROYv8YmLxAevSvHtkPzKxhO/btoHNppLlt14504cCS5ylrW8XWNUt+CCLNKa8SWep/tSeC6KJTqX\nW92EYOTgp5nTMo/VTOGVIzqDR4Ii1VUyGcFFIgTIZE9B1CH51YKRzacHapKmhN8lUy9JS4KLnGUz\nb1dCjOJHpVCRTEZwkbhRkcnXgnHoJgQDl1cYrsxjvUTarzGmupYxWVeqKDahhCgWgk7JryaMat3x\nY0VTyu+SpZigKcERlO7ik1FprguIf0FxVqgIHq+5SPuVeAheX1x8Z4kopXopHbrCQzDS5HmspwhG\nLAT44f2yVCEaqgiRNg6+FuCDs6D0ElPzj1qqfpcsxQRNCY6gNGeZjEpznUPgzsrorFAxYYVI+5XE\nGxE+OFQ3ojcNTDbNnu5W8V1yE2CReaxmCi8TAk3WlSqcySop0sYhQxBMfjVgtMbct7UGIf0uWYoJ\nmhIcQWnOspm3KxB+MLX3CWKXpTULLq2FKJ7ojOJNqsLOjWCQ395G93hje53zRB6rmcLLhSxZV6gY\nKFLkCEE4+dWA0RpjKr2AcQMuad2SkRJc5iyb2ckcZW17+12SzSw3+hW6ELmeghu3nyKKJZbFS88u\nVEcfRc9PuLBA/tlDS2wiEt2ViuCoulwXtHwKBCNV7fbc61IqbDDA5S6+ct1Gkyq2ANvSJTURMS12\naAYlNh85ehqsj2otlcd6QeMapFIx2LSh9p+viQtTcHJ1l7LXJZeyn2eAy1185brdU6rYAmxLl+xs\n30m+BJ5vO3rFWF/QOSKP1ZnCyxSFit6N0a7957hWAk6+6lAKuOTQ9rMMcLGLH3BdqDgQvS4VD1eY\n+tAWJxy9P2D18xe4cFKg72kK6oN6J3+B0p+Itf8BpXtNMcGumNsAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$\\left \\{ T_{0} : T_{e} \\sqrt[4]{\\frac{- \\epsilon_{0} \\epsilon_{1} + \\epsilon_{1} + 2}{\\epsilon_{0} \\epsilon_{1} - 2 \\epsilon_{0} - 2 \\epsilon_{1} + 4}}, \\quad T_{1} : T_{e} \\sqrt[4]{- \\frac{1}{\\epsilon_{1} - 2}}, \\quad T_{s} : T_{e} \\sqrt[4]{\\frac{- \\epsilon_{0} \\epsilon_{1} + 4}{\\epsilon_{0} \\epsilon_{1} - 2 \\epsilon_{0} - 2 \\epsilon_{1} + 4}}\\right \\}$$"
      ],
      "text/plain": [
       "⎧           _________________________             ________              ______\n",
       "⎪          ╱     -ε₀⋅ε₁ + ε₁ + 2                 ╱  -1                 ╱      \n",
       "⎨T₀: Tₑ⋅4 ╱  ─────────────────────── , T₁: Tₑ⋅4 ╱  ────── , T_s: Tₑ⋅4 ╱  ─────\n",
       "⎪       ╲╱   ε₀⋅ε₁ - 2⋅ε₀ - 2⋅ε₁ + 4          ╲╱   ε₁ - 2           ╲╱   ε₀⋅ε₁\n",
       "⎩                                                                             \n",
       "\n",
       "___________________⎫\n",
       "  -ε₀⋅ε₁ + 4       ⎪\n",
       "────────────────── ⎬\n",
       " - 2⋅ε₀ - 2⋅ε₁ + 4 ⎪\n",
       "                   ⎭"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# The 3-layer solution reduces to two layer solution if we set e_2 = 0\n",
    "out3 = {}\n",
    "for var, formula in out2.items():\n",
    "    if var is not T_2:\n",
    "        out3[var] = sympy.simplify(formula.subs(epsilon_2,0))\n",
    "out3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generalizing to N layers\n",
    "\n",
    "\\begin{align}\n",
    "E_i &= \\epsilon_i \\sigma T_i^4  \\\\\n",
    "F_i &= ... + \\tau_{i+2} \\tau_{i+1} E_{i+3} + \\tau_{i+1} E_{i+2} + E_{i+1} + E_{i-1} + \\tau_{i-1} E_{i-2} + \\tau_{i-1} \\tau_{i-2} E_{i-3} + ... \\\\\n",
    "F_i &= \\sum_{n=2}^{N-i} \\bigg( \\prod_{j=1}^{n-1} \\tau_{i+j} \\bigg) E_{i+n} + E_{i+1} + E_{i-1} + \\sum_{n=2}^{i-1} \\bigg( \\prod_{j=1}^{n-1} \\tau_{i-j} \\bigg) E_{i-n}  \\\\\n",
    " &= \\sum_{n=1}^{N-i} \\bigg( \\prod_{j=0}^{n-1} \\tau_{i+j} \\bigg) E_{i+n} / \\tau_{i}  + \\sum_{n=1}^{i-1} \\bigg( \\prod_{j=0}^{n-1} \\tau_{i-j} \\bigg) E_{i-n} / \\tau_i \\\\\n",
    " &= \\frac{1}{\\tau_i} \\left\\{ \\sum_{n=1}^{N-i} \\bigg( \\prod_{j=0}^{n-1} \\tau_{i+j} \\bigg) E_{i+n} + \\sum_{n=1}^{i-1} \\bigg( \\prod_{j=0}^{n-1} \\tau_{i-j} \\bigg) E_{i-n}  \\right\\} \n",
    "\\end{align}\n",
    "\n",
    "Now substitute $n \\rightarrow -n$ in the second summation.\n",
    "\n",
    "\\begin{align}\n",
    "F_i &= \\frac{1}{\\tau_i} \\left\\{ \\sum_{n=1}^{N-i} \\bigg( \\prod_{j=0}^{n-1} \\tau_{i+j} \\bigg) E_{i+n} + \\sum_{-n=1}^{i-1} \\bigg( \\prod_{j=0}^{-n-1} \\tau_{i-j} \\bigg) E_{i+n}  \\right\\}  \\\\\n",
    " &= \\frac{1}{\\tau_i} \\left\\{ \\sum_{n=1}^{N-i} \\bigg( \\prod_{j=0}^{n-1} \\tau_{i+j} \\bigg) E_{i+n} + \\sum_{n=1-i}^{-1} \\bigg( \\prod_{j=0}^{-n-1} \\tau_{i-j} \\bigg) E_{i+n}  \\right\\} \n",
    "\\end{align}\n",
    "\n",
    "And substitude $j \\rightarrow -j$ in the second product.\n",
    "\n",
    "\\begin{align}\n",
    "F_i &= \\frac{1}{\\tau_i} \\left\\{ \\sum_{n=1}^{N-i} \\bigg( \\prod_{j=0}^{n-1} \\tau_{i+j} \\bigg) E_{i+n} + \\sum_{n=1-i}^{-1} \\bigg( \\prod_{-j=0}^{-n-1} \\tau_{i+j} \\bigg) E_{i+n}  \\right\\} \\\\\n",
    "  &= \\frac{1}{\\tau_i} \\left\\{ \\sum_{n=1}^{N-i} \\bigg( \\prod_{j=0}^{n-1} \\tau_{i+j} \\bigg) E_{i+n} + \\sum_{n=1-i}^{-1} \\bigg( \\prod_{j=0}^{1+n} \\tau_{i+j} \\bigg) E_{i+n}  \\right\\} \\\\\n",
    "  &= \\frac{1}{\\tau_i} \\left\\{ \\sum_{n=1}^{N-i} \\bigg( \\prod_{j=0}^{|n|-1} \\tau_{i+j} \\bigg) E_{i+n} + \\sum_{n=1-i}^{-1} \\bigg( \\prod_{j=0}^{1- |n|} \\tau_{i+j} \\bigg) E_{i+n}  \\right\\} \\\\\n",
    "  &= \\frac{1}{\\tau_i} \\left\\{ \\sum_{n=1}^{N-i} \\bigg( \\prod_{j=0}^{|n|-1} \\tau_{i+j} \\bigg) E_{i+n} + \\sum_{n=1-i}^{-1} \\bigg( \\prod_{j=0}^{sign(n) (|n|-1)} \\tau_{i+j} \\bigg) E_{i+n}  \\right\\} \n",
    "\\end{align}\n",
    "\n",
    "Now combine both sums together\n",
    "\n",
    "\\begin{align}\n",
    "F_i &= \\frac{1}{\\tau_i} \\left\\{ \\sum_{n=1}^{N-i} \\bigg( \\prod_{j=0}^{|n|-1} \\tau_{i+j} \\bigg) E_{i+n} + E_i + \\sum_{n=1-i}^{-1} \\bigg( \\prod_{j=0}^{sign(n) (|n|-1)} \\tau_{i+j} \\bigg) E_{i+n}  \\right\\} - \\frac{E_i}{\\tau_i} \\\\\n",
    "   &= \\frac{1}{\\tau_i} \\left\\{ \\sum_{n=1-i}^{N-i} \\bigg( \\prod_{j=0}^{sign(n) \\big(|n|-1\\big)} \\tau_{i+j} \\bigg) E_{i+n}  - E_i \\right\\} \n",
    "\\end{align}\n",
    "\n",
    "with the convention that $\\prod_{j=0}^{-1} = 1$.\n",
    "\n",
    "Alternatively if we set $\\prod_{j=0}^{-1} = 0$ then\n",
    "\n",
    "$$ F_i = \\frac{1}{\\tau_i} \\left\\{ \\sum_{n=1-i}^{N-i} \\bigg( \\prod_{j=0}^{sign(n) \\big(|n|-1\\big)} \\tau_{i+j} \\bigg) E_{i+n} \\right\\} $$\n",
    "\n",
    "Let's use this notation since it simplifies our expressions.\n",
    "\n",
    "### Net radiation\n",
    "(absorptivity) * incident - emission\n",
    "\n",
    "\\begin{align}\n",
    "R_i &= \\epsilon_i F_i - 2 E_i  \\\\\n",
    "  &= \\frac{\\epsilon_i}{\\tau_i}  \\left\\{ \\sum_{n=1-i}^{N-i} \\bigg( \\prod_{j=0}^{sign(n) \\big(|n|-1\\big)} \\tau_{i+j} \\bigg) E_{i+n} \\right\\}  - 2 E_i\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Alternative...\n",
    "\n",
    "We will define the **transmissivity between layer i and layer i+n **\n",
    "\n",
    "$$ T_{in} = \\left\\{ \\begin{array}{cc} \n",
    " \\prod_{j=1}^{n-1} \\tau_{i+j} &  n > 1 \\\\\n",
    " 1 &  n = 1 \\\\\n",
    " 0  &  n = 0 \\\\\n",
    " 1 & n = -1 \\\\\n",
    " \\prod_{j=1}^{-n-1} \\tau_{i-j} &  n < -1\n",
    " \\end{array} \\right\\}\n",
    " $$\n",
    "\n",
    "Then the incident flux follows directly\n",
    "\n",
    "$$ F_i = \\sum_{n=1-i}^{N-i} T_{in} E_{i+n} $$\n",
    "\n",
    "and the net radiation is\n",
    "\n",
    "\\begin{align}\n",
    "R_i &= \\epsilon_i F_i - 2 E_i  \\\\\n",
    " &= \\epsilon_i \\sum_{n=1-i}^{N-i} T_{in} E_{i+n} - 2 E_i \n",
    "\\end{align}\n",
    "\n",
    "Now make the substitution $i+n \\rightarrow m$\n",
    "\n",
    "\\begin{align}\n",
    "F_i &= \\sum_{m=1}^{N} T_{im} E_{m}  \\\\\n",
    "T_{im} &= \\left\\{ \\begin{array}{cc} \n",
    " \\prod_{j=1}^{m-i-1} \\tau_{i+j} &  m > 1+i \\\\\n",
    " 1 &  m = i+1 \\\\\n",
    " 0  &  m = i \\\\\n",
    " 1 & m = i-1 \\\\\n",
    " \\prod_{j=1}^{i-m-1} \\tau_{i-j} &  m < i-1\n",
    " \\end{array} \\right\\}\n",
    "\\end{align}\n",
    "\n",
    "Or using the Einstein summation notation, since $m$ is a repeated index, we can just write\n",
    "\n",
    "$$ F_i = T_{im} E_{m} $$\n",
    "\n",
    "and the net radiation is\n",
    "\n",
    "$$ R_i = \\epsilon_i F_i - 2 E_i $$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reformulating in terms of flux between layers\n",
    "\n",
    "Let the upwelling flux be a vector ${\\bf{U}} = [U_0, U_1, ..., U_{N-1}, U_N]$.\n",
    "\n",
    "If there are $N$ levels then $\\bf{U}$ has $N+1$ elements. We will number the layers starting from 0 following `numpy` index conventions.\n",
    "\n",
    "- $U_0$ is the upwelling flux from surface to layer 0.\n",
    "- $U_1$ is the upwelling flux layer 0 to layer 1, etc.\n",
    "- $U_N$ is the upwelling flux from layer N-1 (the top level) to space.\n",
    "\n",
    "Same for the downwelling flux ${\\bf{D}} = [D_0, D_1, ..., D_N]$. So $D_N$ is the flux down from space and $D_0$ is the backradiation to the surface.\n",
    "\n",
    "The absorptivity vector is ${\\bf{\\epsilon}} = [\\epsilon_0, \\epsilon_1, ..., \\epsilon_{N-1}]$ ($N$ elements)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "epsilon, epsilon_i, N = sympy.symbols('epsilon, epsilon_i, N', nonnegative=True )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Will do the 3 layer version first"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAEYAAABMCAMAAADnazHQAAAAPFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAo1xBWAAAAE3RSTlMA\nMquZdlQQQOkwRCLvu82J3WZsMMTa8AAAAAlwSFlzAAAOxAAADsQBlSsOGwAAAxRJREFUWAntWNuS\n2yAMlQGzW1+CXf7/XysJEPiiLMnOdKbT5SEmAp2DZOxjAaC3VR+6jAyRm7kMwBj91XiybMkbYIjW\nYRtP4wA+dMDs5GsiwQwXBDY43wHDM/cGZo7LakyM+POIAeHHt2ACugIsG/3aHSDAWzCc4zlagtkn\n2Oe3YAZcAPpHTrXz3mGO38hNyrNJnhOMFluk4Dpam2KenlKTPeeyGh+MOQB6u/B+4QxQDMcbnlOT\nYPY1Gg5xWia2uOD4NsCw2Wnc1nnOfGeYnJo8mi9zyhfsiDHSjfALoZNvbmeYnBoAIcaJJtoQMKaN\n1kRx2ge5TyXka1AlNUJM85fVe3T2vNtpIZHz5HiH0YRzbiQ1QkyT0lbCh4ZSse0wp6fnkfN7hSmp\nqcQE80gP/8BBbC5vp2Gpj/8pN2sOV4gJBSZ28Mm44LZcMbCZF8zjx6DC+ojxYXAarp+YkDi1ybid\nFkE7NGKePb5ZQl3LESb70KUSN8Z6p1oj909BlfFKXCx0dZhT3jetkfsKjEIcbLBtLAKnwejE4tp2\nNBhQiVtv6aswMqOr83dgeuXu+Wp65I6DfgrTJXdfw/zIXfkU+F/lDgZWFnytfUfuBmuSpvTLnRDj\nLi1yh+++BHMrdx+fv+RtUOROiGkkyx32MsxF7n5/Hj4FRO7Eg2Cy3AlMt9yJB6EUuRNj+nrqkDvx\nIJQid9XYK3fVg2AgyV019spd9SCU2nKKqyH39NfWrce4XBDY8BLMbLfIEn/BegKjEF8gyKDB6MQv\nwdxO1o3aanSP25EvYH7kLheJP9Xdcfs0++afq+6CXdN38bfkLmAhZEg3++VOiNGryN2GMFxN3Mpd\ne0ZR5E6I6Y4UuaPSZ6YK5iJ3zZ0iB5E7ISZrlTtUQSyR+uVOiAkmV3fU5RrrFblLxOyaqzvqB64v\nX5C7prircrenKvUFucvEvJzyM2E9P92dGZ1SXObjNRM3FvzkGcdxvavLVJg74nS60QKXvgajExfP\nw1WD0YkP7uWPBlPGO68JRjs+7AGR40NPB4DO0anKG42PD52DPzPNNPozuwviAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}T_{0}^{4} \\epsilon_{0} \\sigma\\\\T_{1}^{4} \\epsilon_{1} \\sigma\\\\T_{2}^{4} \\epsilon_{2} \\sigma\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡  4     ⎤\n",
       "⎢T₀ ⋅ε₀⋅σ⎥\n",
       "⎢        ⎥\n",
       "⎢  4     ⎥\n",
       "⎢T₁ ⋅ε₁⋅σ⎥\n",
       "⎢        ⎥\n",
       "⎢  4     ⎥\n",
       "⎣T₂ ⋅ε₂⋅σ⎦"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# vector of emissions\n",
    "E = sympy.Matrix([E_0, E_1, E_2])\n",
    "E"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlcAAABmCAMAAADRXESXAAAAP1BMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADFBd4eAAAAFHRS\nTlMAMquZdlQQQO0wRCLvu82J3WZ8bGJP3RYAAAAJcEhZcwAADsQAAA7EAZUrDhsAAA+ZSURBVHgB\n7V3porMoDKXr9033zvj+zzrsCUkAbdHbevFHiwrZOI2InKpU38ZH4Dq+aq/ZIzA2AofhOLbqr6+3\nGey2/fWBqAfgeO+4qkXp4eCk1GbY7fV2qDXo59X+2HFVg8HTgGk7GFxtanX7eRuBzaHjahwUnh1X\nhUCdh8t1ux0G/XEb7krdVcdVIVzoVMcVCgYr3jWWlLo8zOfuqZ7njisTihFbx1UpSPZW5jzsTJ3n\n6bhXHVelcKFzHVcoGLS4eZojz8He0eyPh53eBp23+laNQMdVIUTudmbrbgFPtuK53w8WAganOq4g\nFpmSG165k8/rsHXzMcfdxc7T2GtkpuUvPtxxVet8P7wi1TaP3enwuJ7P5HjfdRHouKohwQ+v0mrH\ni8laJnh9EyPQcSWGBR30wyt0RBd3N7N/6oOtNCyw13EFsZBLeHgVawz2pnBvJ7biwV6ACHRcQSzE\nkji8OrsnX7c+aBdjpg92XOUi44/Lwyt7Adxc+pqZXPQ6rnKR8cevMIa6w5qPqx62nx9uSqsi4Hee\n7rgq9fv9ehuG21Y/v9Hb7qQ2YbXoUa8suvdslY9dx1U+NuTM9Xo4diiRoOR2O65ykWHHN3rBTF+o\nxsIiH+i4kuPCjx701Lq/IvKT/QiJQMcVCUh296Fz1a7nq2x80hMdV2k88nvP53PfbwDz8UnPdFyl\n8eh7bSKwalyFWYE2oVKqtbxWdn2gnDXjqjWHtLW8D4RDM5NWjKvWHNLW8pr14ScKWjGuWnNIW8v7\nRDg0s2m9uGrNIS3JO7fiUmgm2Tq29eCqNYd0grxzuwH9diXAWg+uWnNIJ8i7vQ+Gg19zc7QLUb8/\nZ60HV29zSA9u2ULo0/HyntOX96W6jttdXO18T60I1nzb92pw9T6HNO3rCfIe09NVqkuDZh+WeR0v\n3wYh0d714Mq69w6HNO3r8ZzU8wur3FNdGFfqtoqHRavBlfvVYJJD5JAe79ttesMmkkpZX2uRsryU\nkrrPXQa5Xmek/mS6Yr5Su5y42PgbCuvClUhyOF1cBtjf9/bvYfQffomkUtbXeqmx+8uPtCNp62uK\n2ViZ642nSrh6rmLkvi5cSSSHs/vbDvXUoDrYwXiGVCrgSpLHWj8csR5A40pcL6rBdEG+OrxwWUWS\nP6S4LlzFuyrITWo77O53nVIsy8GOjhmp9H4z2+Nivyz0fO9Eeai3WGufDo+3R9xuenTG9VohGV2A\nq80qSNQGV3+GPyhu31wMwyHITXqIdD2aZelHu4bY0t8zpFKWQ9LhVYgLa+1xFc6Hb643nNHfTNfa\ncPXvmv4HMg6HIDcp5YdIm8HMBjx04sqRSllfi8Mr3jpzHWR6EaxKuDr16yCO1CeUw3AI5Salbu7C\ntrGXwEf8yz1OKuW4CvKwb+4PRnHrzLid6cVSmC7IV33cjgP1EeXAIYXcpM062SckR3fsYmazM6RS\n1tcqyNNtgJPKWt/liQGuF8WI6QJc7fxNK6r9hcX1jNsRhxRyk+mR03b/3Pvx1WCmHDKkUtLXSF7C\nSWWtcxcuphfBg+jabS+a/urmK8yY//u39eAK9QXKTXAUj7ngKCqRvkZndIorcVLLz3FEvXldr8ze\nJ6Z+xs4qceXGVzY3QZTNtLibv4JjSelYeMxX5KSWnzuLevO6+nPnpFM+a0fMEffdffciD77CSb0W\nEKmHZlP09nUyn4Wk1BoxR6RVJu1VOKkN1/WVITrJ6B+tvMrr4MQcUe+AGie1r0OmMVwprqibfX/h\nCHRcLRzwX6LO4Oq4f3E4+3aM2vEN3jaFC9iIC+xam9xaHvcjHllyhfNZK/ux9w9+NIE4/jdf7BdT\naG1ya3mJsWTnhJdqkHOtd3/yOvjZBOKblMRbm9xaXhkeZrnQQttP4qoFgXi2O7GdeNVoYTLu2ZK8\nGVx7SL8VbE+z8oK4okTPEoF4rH8NZ45SRqh7b1drk6fIm8O17FL8sfEeXW9BXBGiZ5OX2s7GCHWL\nFFqbPEXeHK6dF1uLuiCu7KjRL717nl55qS17Wlt+Mif9uIiILCPULixVrU2eIG8e18y6xkW25XBF\niJ6vvNSWgEIv/yw+mJMCyETIjFC30qa1yVPkzePabqlpjQVxZbs5EkdfeaktBcULa0qoCMQ0xozQ\np2UdE27q2yZPkDeTa/ul2NTL4cplj8BssHvxpbacwTmKOZodhnJ5IXeVcIUYoY50M8lk0WK+kJ2Q\nMSAEb5Ndx7i22JvtFsZVZDaYLouE5MDgVJvAyaTcT9fDNHKZleV66bGfKo/yXHv9SUXgfIVWll/h\nelE2ORLKZIsFfbI82rzmWtQbPRNUweLm4NpiU+AL40piImjai2N2bnZbz0Vh3E8ZVxkmDJc3LviI\nEepJD6Zd0eRIKMtYLHS2KI81r7gW9YJngirAVXDtTBY7ouZtiwvjCoieKJUEBqfOJh5XjPspszkr\njFCQZ2Imi8D5CjFCH8CFiCbjFBFMjgsImcVZfVEe7kjWvOJa1GuFjHct/OSw7lnKBld//v4zi2xB\naBheQWrSlTyDU5cCrhj304miF7FwuSOKuDyoQEVkcIXyVTA5SRGMdJqxWEgiQR7YpEusedm1hMgW\nBI1xbbF89d/fJZ87o7FFgJAJakwO/iDnfrrY0chlLhZMXoi8/qYiMK4QsQZwFU1OUgQlneYs5vqi\nPGSUQJUtu5YQ2YKgMa65hqHFjN/LXgfR2ALhCjrRH+TcTxcBGrnM4JbJQ/GjIjCuwuBWV9/GcXsw\nOU0RXkUklOUs5rgK8pBNhuRvHty9SXYd49pK7wcR0RPhyjM4dWTDQcb9dL1AI1dmhCJ5rrn5pCIw\nrhAjFOZ5gslpimCk04zFXF+Qp20BriunypZdE4lsY1xz83IQjtlKC+arhOgZIWQ8cwxOXQi4YtxP\n5z+NHLpwuQr+k8mDs1REhhHqf9fI5JianCxKOs1YTHCF5CVcV06VLbsmEtnGuLbC+XboW1sKEEoO\niwehBo1c9TkOl8dEgPRkipuNm8UUgf78COQkpYK+Ite14loy2PMKC6qiazdx9U9icZudBfNVajDv\ncn3+UH7MwNictYezXB4TAVYljFB0TXQ1xBShqoSygr4i11WVXZP0FlQF19w4Dlyer/RJuDrvHsPE\nN5IW6XbT5KWM0A1DuJQi3iGUVbiuqujai2TX52IrkX8OV6zjXvnxNFz8RvpxS//dUUoRr5gc2lS4\nrmoO14zOZbYfwtW0VFIIxQyLdZ22NH2ZY5P48AWT/aka11W1d81M7S60/RCuFvLuHTWHOFv7jpRP\natswA1bd6rjKhugg8gez1T//xKv/evKKZzVcxXnnV4TP3WYZ5mhLL2SLm7/gd8Feg4kLKOmIVXC1\nJG1ycv8twxydbFahgWzxSsiuCe21jKtlaZOFDhFPLcIcFTW/elC0WLWOcmt5ZW+B7AqlWr4q0SbL\n2uLZZnc1hsCTbLMxR5uZPM5i1SDKSWAK8pq5hnoDyK5QSq+DU2iTiSf5nZb3IEswR7UnDU2WLNaL\nYi7X7XYY9Mdt0Hf+7/JzJ8hr6JqKvgHLAEoprqbQJvNQSs40YFfqpzuO/p1OKc3DHNXGNzC5ZLFq\nHeUJ8hq4xnsDyK5QSnE1gTaZgAd26LPP8lMuaIdKRMTczFGiTo8Mps9bpTIqFr9Pdk3VTZD3tmtK\n9A3IrlDC4/YptEkEBFwkLleeyuOWsUxFoBVS+F2irZijTN10QihZDKM9iYQFbrF6P8qpyRPktXBN\n8A0W30ApwZXt25bM0bg+I6KmXkijZurHXpqDOUrVNTG5YLGaQE7NBCs1eby8Jq4JvQGLIKGEceXc\nSNb1A22SvnN0FA0TDeTSGOWJo/zHjzyZwhzlKkST007SIM5dBrm86BKVUbc4eRGrHsn7ASQNsl7v\ndxnMlhjF1GlDxF4jrZu4JvgGi5uhxHCVrOtnzNH77uqeBlAepY8xcTmzAD0SR6O42EXllcJoBXqN\nOcq5qbLJxGJVMxm4s2AzlYFiL1uskYSQEqIcTAZC2TiTtSGiPNq65hroLbgm+AZkVygxXInr+gPt\n7K6nkLaGk8x4lN4SEuEMq4SJAz/KuEKLAYEcoe896JoWE+jluK48xcJ1ULa4aDIQyjJRJkE2wZNC\nwFpXegP0ou7gurhvQB6DEsNVpE3iXBJpmBpXdkEu41HKNMwauzKIc46Mp1fq+hJzFOWSYDGsmWcm\ny+oqJoM8a7MsA/2mZa6rfpmq/+c8lCOCybCAcKTJ2pIoz0XSfrLWFddAb8k1wbfwG4ZfszDfHi7U\nMTUZJYHpaVRbTDIepTWFJRu2StxVY+J8a/tV+oWgXkL5KpiMya5BBeAgYzJVVzEZ5CGbqQwh9ro2\nsjgOh3COYFxXzlZ1Opk6MrzyljGHy66lRLbgHdcF+Sr0BmQpKNF8FS/U5p7UcQW0CjwaMI+ix9Iw\nM5mXigtO2O+SJ4ilAr0UTcZ9Dir8lShnMlVXNRld2YLdVAbClWgxDIdwjvAmu6B/7ItdBd8iTAAw\n6byojlO8UMfUZIIHnagvQHo9BOdR+hCTCGdGilScb+y+iAh9EH4haBTMmaO6JvQ5qPDHciZTdVWT\nQUc0m8qoWRyjnOQIbzIQysaajHotmiT1Udk10IuE8KGj4BvcBUKJ5itEm9TdFEYBwBzVq3GN2pE0\nzDK70gjiNVgnIU8QSwamSpDJ0OdgcTiWMZmq4wYZM8NbV3UpyLOH3QeVUbM4vogVchOocMc+9sWu\ngm9AdkUl/N7whDapPTWpyW2R6ekfBYyjYSp0FQiS7DcVh0/SThrNHNVCUJ9HFeFYxmSqrmpykIds\nJjLKFisU5TRHUK4rZ6s6nUQdkpeQXZnDZddkIhvRpSTfYJYdSjRfoWgJuUSd9PP3U4HUQc2oPDmQ\nxFERyKJkwlgahQp9jrGGRMUiU1cxWZLHZETpFYtRbopN3uK6LvhiVxV9A7IrlAq4Ep5Sbm6Hw+Hq\nb5EhElCi3EhBBlTW0wKCOCoC1Q/0SnsIXRNjFRFXZUIZU1c2WedELo/JiAZVLBZzhJkW/4IXu+q/\nl3BuAtkVSmzcHiMipib3UAHq1EuElZc2mCguXSfDmaPJddBreoFQVjR5mryqxfh+MMbmDUJZhexa\ndG0ikS36BmRXKOVxJeaS6ProQsuVZCQqjDkq5pLRlsaKDU2uWlzNTdGqcYUK2bWha8DHBrIrlPK4\nmphLsm7PsfLVKYs/mKB7Wi4Jrfh3M5PpOmRm8eJk12auwTpkILtCSZhv52H+2CPfxxz9PournQ8p\nEEq6kRu327UYi/0nRNXWsRW+jzn6fRbX+gLIrqH0sGiy70vdm21t9N5aRPr5eSLwtGjaq/8BKzXQ\nJDvqpwwAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}T_{s}^{4} \\sigma\\\\T_{0}^{4} \\epsilon_{0} \\sigma + T_{s}^{4} \\sigma \\left(- \\epsilon_{0} + 1\\right)\\\\T_{1}^{4} \\epsilon_{1} \\sigma + \\left(- \\epsilon_{1} + 1\\right) \\left(T_{0}^{4} \\epsilon_{0} \\sigma + T_{s}^{4} \\sigma \\left(- \\epsilon_{0} + 1\\right)\\right)\\\\T_{2}^{4} \\epsilon_{2} \\sigma + \\left(- \\epsilon_{2} + 1\\right) \\left(T_{1}^{4} \\epsilon_{1} \\sigma + \\left(- \\epsilon_{1} + 1\\right) \\left(T_{0}^{4} \\epsilon_{0} \\sigma + T_{s}^{4} \\sigma \\left(- \\epsilon_{0} + 1\\right)\\right)\\right)\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡                                    4                                    ⎤\n",
       "⎢                                 T_s ⋅σ                                  ⎥\n",
       "⎢                                                                         ⎥\n",
       "⎢                         4           4                                   ⎥\n",
       "⎢                       T₀ ⋅ε₀⋅σ + T_s ⋅σ⋅(-ε₀ + 1)                       ⎥\n",
       "⎢                                                                         ⎥\n",
       "⎢             4                  ⎛  4           4            ⎞            ⎥\n",
       "⎢           T₁ ⋅ε₁⋅σ + (-ε₁ + 1)⋅⎝T₀ ⋅ε₀⋅σ + T_s ⋅σ⋅(-ε₀ + 1)⎠            ⎥\n",
       "⎢                                                                         ⎥\n",
       "⎢  4                  ⎛  4                  ⎛  4           4            ⎞⎞⎥\n",
       "⎣T₂ ⋅ε₂⋅σ + (-ε₂ + 1)⋅⎝T₁ ⋅ε₁⋅σ + (-ε₁ + 1)⋅⎝T₀ ⋅ε₀⋅σ + T_s ⋅σ⋅(-ε₀ + 1)⎠⎠⎦"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# upwelling flux\n",
    "fromsurface = E_s\n",
    "U = sympy.Matrix([fromsurface, tau_0*fromsurface + E_0, tau_1*(tau_0*fromsurface + E_0) + E_1, \n",
    "                  tau_2*(tau_1*(tau_0*fromsurface + E_0) + E_1) + E_2])\n",
    "U"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAawAAABlCAMAAAAmhidWAAAAP1BMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADFBd4eAAAAFHRS\nTlMAMquZdlQQQO0wRCLvu82J3WZ8bGJP3RYAAAAJcEhZcwAADsQAAA7EAZUrDhsAAAsgSURBVHgB\n7V3rorMoDLTa+u32Zrvr+z/rEpAkkEhrUats/XFqUYYkIwFxjq2qxHZLHNvmobktnhsvM2qH3m61\nAnPqW6V0I0WHs2LI3BbPjaeYnCg6+mMPx1FVHfrmaLaTP0Cfbbdhsg7aNT+3xXPjUWzf2jv7HvQE\nhuoeyDqM1Dy22WRdniPYU4ufl6jGXev0M1gcNDMNb35nOx6+Jyfr0l9vdd335s+97wyLp2yyLtrl\nH4Tj7S91yFYDKWJui/PwlnD2wS7JgKzOEFRV1wf8bQylXZVN1j2MMCBP3U5XZ2975zUvts/PbXEe\n3hLOHhvyOiDLJshLbw8/z5VJPBPJOuGAOLTwZE1Ro6m9CKKtm9qn4o6jdxY52+IqbC8LbxFnLyb1\n+Y2TdbD58dnbmcaxbU1scsl6TO5YYfDAzKMnq716q83nFazMtzgkKw9vGWcfNGoFZNlYDBfyuTo1\nZushHb69xZG+2Iz6dnU4MYZgZFV3mqsfLINuYpRjcdheFt5CzjY06nOyXFDdkDUE2I0Mpod1dR2y\n1jZXO/UP8lwcaZ5wB0T3IfH84RiCk9VQY0/Wy1SLpcmqxSFZzogP8RZy9kieCrKGIcuZ/bz1NaSb\n85WuaXfk8GjOp8ftEqa5ONK3kGBXVcXDQzEEJ+tJU4yGdlWLpcm6xQpZn+It5OzZjwJVJcgahiwM\nHuxc3DBm4tYd7YSxau2QAZWDLY70Q7nRVvEIJYbgZJ0oq954clBaQZOrg6N1xGKFLC0CFIKuuTUw\nO5V4r5z1lpCvSuM4QFfeWXYXLMjCuRfDrPumg7uzp2HqZCdM7sJmpHd32B5X+2FPgfpDh2zvD9zu\nZlyQeHCyDsHJOtDFccc2DBq7FQEg2HwTh6Z2DAuLR9tL4nUmldRAv8R74SxaYu1739lLj1lNkIUJ\nG3uRifmtbSEgD6hmQ9Pb/HakK91aIK4UkT3daRLPlcPfVM9iZD1o+EKL+aXrmzCAzsgRi5X2EM93\nIjDL48GM79CbPxLvhbNoCcAN2zvOUo4QaRATNvUiY5iLTGuXpSABXtwC1Z0i5lqPGx/JDALPG28+\nYwjes850cVDPQouDS3doAgBtpTGLZXuIh50IrPN4cL3Cta7gpZ1FSwDNb+84m+hZmLCpF1XVEBl7\nRVUw73e3X4dhacE3LSM9MuYKPEKQwaP7rIpNMIgstNjzYsHouCNrzGLZHuJhJwJAwjPNmOSi4KWd\nNRh+FLIG2j8psryzLur29DgN3oYBgHqROe1sWWndvc0D1hFu0L0snxbF/4kbd8sM/ih+Sjw8JIPH\nyGrc7AZOrnGC4S02hSwaQxNUOGKxbA/xfCeC5nwI7P5IBNLOmprMPItjSviajC2jCYZ3lk0MArK6\n273v7zUgUC8CjHN9fB6Hsiscbs1DlU6O63HjLG8BDG4CD48I+5v6akxytwAwORm24e6DWWwO8Gi4\nJqhwxOKovRDPIOLkBfEqR4rESztLlgwewEccL81ZdkcZkMVgWC/CUtfbaHKCB3AnbrxKrsCoeAIC\nwU1Xpn12vVEhJwtL1UI8KuJFR2DPdqKwaHwFMOmsAZGWvONscgXDmeZ6lu1FZCsfx6iU7bXhPbKZ\n7MdTEHYyn11SsYCgQ8FCrjb3ktEwldVCAk20Z+4mpPlnk4rP1MMJ6JWzmiWJxtHZO6XKsZ6lXvWw\noOLus7iNyf1bTB8/expe+IjEp3SOp/JyouUafu47+8qldrifTqebHAAsXNJZE7oJlqCzfCl9jCx2\nT8Xc6prO3r6zohe76edxk/DCUBwUzxWyLs2jt4PwCzu1w1oncgui2tmmLOnsNEvQ2Sfd/Yv7LDRj\n2lWP1cTO/E+6hybsomXY3JRLN6ypfUt2Iq2CYctNhPSDU0pJw/BgKXe0Z5l8PbUXTbEm/1xMFB5q\n2qXra41/JjvReLV5j8DaBG7jZOEpW905ydF/q6Z+bFeYWL9LFt7ZfuTNCVc4ffU8PI+yoc9whvBV\nsuiWc574zI03j1XzoXyTrLkFlHPjzRflmZC+SdY0AeVrhyfhLTBxe21h5hkrkpUnoJR+ZuGFI7cE\nn1ASiU8n1Jx66opk5QkopWNZeEvoMaWJM5esSJa9FR8e7c0hIc3BUxaSXgQ2WnMdFZ++gMk6vB5Z\neQJKcDKMVxbeqwVyGdOwcTiOj54C8amsOV/JimRZo3MEmWG83DLMZ3j8YcubsQwbh0pIFhefvon2\n2WnrkeXsQzUKfM2UkALEZ3gL6THBoCW3lclCNQr4JCSkg8TPPKd+S0JqIJJ4pE6K8UbkEqQMRUsw\n9qme5eUSePJCOyuThWoU5o7XWpE6SQoo7ekyXuahgRR4ejxSJwm8ESGSr0mWkJmycUqDyqMZqjjj\n3spkaQJKr8ekR6lCQDkmiVQFnh6P1EkC74UekyyBSI81TmQxPeOM1EgoIOtP/0ceWKaEhhhKNF4/\nSSGSAkprjby4aciilEd6zEHipwgyNU2AaUJaQlGQja9P1j+GLKamJuMW2cMhhicar59EshQBpbVG\nxgvxKOWRHhPq2KVdiTeSBoUlLAiycSJrTNfEqs+yu24aZEMM5XnSTw5lioDS+irjhXiU8kI9plUn\nSbyRCYawhIVYNk5klTnBQAEl9iITDqHHfFdCarSmXtbHBZmE59VJQuCpqJaAF6pJl5LnK0WWJt7x\n9eb8XLFnhQJKFg7UT/oyKaC0LkfxCvGGlAcnIp5fVBJ4Y3kLa3pLKNRR45Wmx6Szl9lbkazQARkO\n3tvCc4dvcbyik4QgU1MnDXVeLTdJ6xKNf7AeEpn+5tdtkaXIy5gfCUmkOUuktpQ6yfc5hh7sSp1U\nonHUYwYQC3zZEFl56iQZ/qQ6CXV5WkynWSJkVhrkLGUbIivLn0TKU3FnfPiY5F1t/NPC75GVTnkT\n/UmlPB3q91hfj4tSOi3RKABRUTLlRefu9+vXetZ+Q/Y9y4Gs9jjyXxHfMyuj5eKEnhSLi/nvn/XW\nBqndxfZKFnqWlgaLFnqWRtYkoedi3Xsh4J2TlSX0XCiky8HunKwsoedyUV0Ieedk5Qg9F4rogrD7\nJitL6LlgVBeC3jlZNiqfCT0XCuiSsPsmy0WGVDjmuxeOvvsmzyWDOzd2AWShagZi44WjqNdE2VMs\n9Jw7kivgFUAWqmZ4uLxeE2VPQujJz97JfgFkoXAUO5GJvVB6CqHnTgjiZgJZf/76mxftbd8PWdiJ\nwAGv10TZ04hwdFfO/vvX3hdyccji2kF88yaQAWu7Uui5K5qcsftPgzhkYScCz0iv6d5DJ4WeP7K+\nEAEUekLb+ICE9JqD7EkIPb9gam6TO+9ZsdCTtIOo1xxkT0LomRu5L9TfOVlRxIR20Mhzx1/mGFXe\n/teiyJLaQfOq98TLHLdPT2hhSWRpnago2VNBZBXVicIuNXwriKyiOlHpZKn+FVVYUM8qihfVmR9Z\nali2Wfgja5u8qFb9yFLDss3CH1nb5EW1qlCy2qZ7duz19arvuysslKy7eaVT++q/vH9kbSIC7jen\nOvp53E1YlW1EmT3rZpW6+HArO0obASiTrN6SdVbeb7eRsH9mRpFkte6NWeeefifss+hsrFaRZF16\n+9s4B/exsYhnmFMoWfYHgH5kZVwYq1X9pcHVQj1DQ26CcfpNMGaI5eIQN/t+had/G+Hi7a3UQJFj\nVuVuip28faVArtFMmWRV8BLc9lra4mChZLXd8ViXxtX4T9+u0a1/bUyLQKE9a1oQ9nL2j6y9MGXs\n/JH1I2tHEdiRqb+etTuyetjsM6Admf6/MvVhKbJv8jzCJn78/H8VjY07+7QUHav/AG1jh1YO9V/1\nAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}T_{0}^{4} \\epsilon_{0} \\sigma + \\left(- \\epsilon_{0} + 1\\right) \\left(T_{1}^{4} \\epsilon_{1} \\sigma + T_{2}^{4} \\epsilon_{2} \\sigma \\left(- \\epsilon_{1} + 1\\right)\\right)\\\\T_{1}^{4} \\epsilon_{1} \\sigma + T_{2}^{4} \\epsilon_{2} \\sigma \\left(- \\epsilon_{1} + 1\\right)\\\\T_{2}^{4} \\epsilon_{2} \\sigma\\\\0\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡  4                  ⎛  4          4               ⎞⎤\n",
       "⎢T₀ ⋅ε₀⋅σ + (-ε₀ + 1)⋅⎝T₁ ⋅ε₁⋅σ + T₂ ⋅ε₂⋅σ⋅(-ε₁ + 1)⎠⎥\n",
       "⎢                                                    ⎥\n",
       "⎢             4          4                           ⎥\n",
       "⎢           T₁ ⋅ε₁⋅σ + T₂ ⋅ε₂⋅σ⋅(-ε₁ + 1)            ⎥\n",
       "⎢                                                    ⎥\n",
       "⎢                        4                           ⎥\n",
       "⎢                      T₂ ⋅ε₂⋅σ                      ⎥\n",
       "⎢                                                    ⎥\n",
       "⎣                         0                          ⎦"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# downwelling flux...\n",
    "fromspace = 0\n",
    "D = sympy.Matrix([ tau_0*(tau_1*(tau_2*fromspace + E_2) + E_1) + E_0, \n",
    "                 tau_1*(tau_2*fromspace + E_2) + E_1, tau_2*fromspace + E_2, fromspace])\n",
    "D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlcAAABmCAMAAADRXESXAAAAP1BMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADFBd4eAAAAFHRS\nTlMAMquZdlQQQO0wRM3dIu+7iWZ8bO+SI+0AAAAJcEhZcwAADsQAAA7EAZUrDhsAABYhSURBVHgB\n7V3rYvKqEo1a3fu0tuo5vv+zHm5zH0hISGv3/vxRMYE1a2ACJGGVadrvcxsMPRpvPb3Du1t2NMHR\neC7p2sG32ollxw/P9Dkuy92T6/y89GSfzTsab9ZgNcPBb/DRBEfjVf1xT7yvColrDqdpOjxPb+Fz\ndrG3HLzcx8bVaLwNvn24F8xogqPxeh2+P3pLhPyPGEzHZ4yrw4riC4q8XcbGVTfe55p68Rx7fMqj\nJ3+I6CYoUc2vLrw9nL26l4/haQ88xsXV5/Prdjw+n+HPx/MeTB3O2+JKA3bjffqDla2F+SNHEVif\n+XrZTFDZ3YS3i7NvJ0Vx6c+BcXWPsTR9XePfU+wo7tO2uNKA3XgfIhgir+7P+StfsZcPXvSeq3sz\nQY4Z0pvwdnH2M3Q7qz4D4ypN8j6fqcof4W4pjB2dcXWWo4sC7MZ7dF9sisDleDrCSH7n3L7yXHQr\nwUnZ24K3k7PXlTOJcXF1SAwez1Tlb5fpEhpiU1wpwH68a3d3pdo5XKlvEFeXL7puD/ngZoIqrjbh\n7eTsaeVUYmBcpXov13fors6n8HmmAZFapJ2SzZpvJhCwG+8zjchti+qsJBBPYlxNH/TA6pFjbCvB\nUEe8E5y24O3l7Bu7nFRlNX+Oi6tsJk+vwGSZ306X+/EoO9TL6Ss95hBDlaxmBxDxvNJgFL6rc07L\nBopYAhRXJ6J64pMt4TESHO3wHN5ezr5Dfw1VtPB7cFyV6VUx/rg9j3FYfP/K1/rb/S3N7cOt4vX0\nfr7ePuVIZZt1EoCA55fWHt9kJONpywZPWQIUVw8KphsbHFyCYGI6QCGfsrXnOwx499PtlO8jNN6c\ns8Skz9m1z6AGx1WZXiH3lPjMc67pEYLqnOamlzTxjbbFx6lmD7BSWkCFH1f/Qa9lQwUtAYqrMw2r\nH8mJXMwjCCYOp2MpVKFs7U0tvHu4DI8pUg3ejLPEpNPZzycN/1R0PrUhri4fV/x8lCerePtEXVOo\niufpHp/cXiPF1K/mgYT1sfeP+Ll+pS/WaqGsfTJnSntMwhOPXCP6pGUTa6lGgOLqQFfBlYZERhD7\nJnQ4XEYlrgzlmj3Co74J8eLc/PCMfbzBm3GWmHQ6C5dILNbz2RBXnhmYbFDXFHJ93S6XEB2X9Fg/\n9VTPNES9UQ+QsZzLFwC5sUppniWmS1Wbw4YN5bAE3Lji/RUQ5D1CcTjgQlxVKFt75flfCHTsm7AC\n04WZOxCD13aWMel09gf6K2KIKZxsUNcUzuVHWuVKiw9EPvOLow922ScIW80IiCbqpVmWlKwMDYYN\nK2cJUFy901XA4ooRhBhChwNwObbCYeqbOF5AjB24xWs7S0x6nc3dIyu1MDm2v4LJAeuaAo/SDPmp\nzxUfax3Ko2xiapsVACkPPBSzpXmemK5MZQ0bVs4SoLhi8/YjzdsZQYorirtyLD/Hs5StPcSLVyZ0\nFoQXJhPx0YTFazsbihA78Ncad5xlcxUot+h7bFzdymwoRzk8q31PAXTJB79ivdzClPoz9WmCpPUU\nAEO2O664qJQWULGA7g5zBsuGCloCVNWncisbcrOHOowgtVwxEXLCsQpla4/hlb4poBAeOGXw2s4y\nJiFZPta442x5VAdlFn8PjKv77eP5/DjGuKGuKRF5P7493sr8Kt1fXMLCnLudkCtPGWCYqb5PsOap\nUlr7zAYuccqwobOKwHQ6fgWP8vMKuDUJueEiFgQxhmKG5HBIQFxVKCt7Eq/0TZEd4sHbGoPXdpYx\niXD5o4z7zv7883agG79Z18QOi0kXO45J7SmeCInb7Rwn/z2f9qsNj02DgHig7c6SIYY4Re8YO9+w\nF3LZPug99Jnv5dabwcRk21mKcCrWME7OfsRuYsVnYH/FrOf5lX70EZ8J5+dXLCtPXuItdOVzCCtw\n/BqtFAgvvv2BsOT32DQIiPfObEwk614MnduvQRr2PPqHj/P5fPOvrrazoeINk4ZxdLbzBS9Vxj5x\nxR5Vkal4Ad7L42J+dFH6HG+7e6+dWyNOO9nIdTIH00jBCRtXn6frs5s0VIfTN+VXX5BBfTed7WNC\nzj74s0RlsPlzp7jyOoMmj7mT19BXnTr7q2ngUjfVaun1lOJsewSVoetnq29ygXZxNlb7qs9OcbWl\na3L9eDweb/1vFPZYmpvo0RUNbPt6BChV/272TW6xHZyNj7fXffaKq3Vsfk+pc3Pu9nv8aDLd0AX+\niatmzdZPnvt7zzrYi55ZOxkO7uwdV/RkekzdvTreGC8dlNGOOyZGHto5rkbrKl8db2TTCKzRjgvw\nHX7sG1ejdZWvjrdDA2XI0Y7vRhSA942rLl0lUGp8/yjeDvdbDVflqZbjP8lLsmS/hsbVJl0lIwXJ\nl8LbcHME/sC3FLnCUfbd4/h38mIUZ5JD42qTrtIh+lJ4u+g+HafjoR7Hv5NXha5zeGhcpYf+ZaVb\nUKZ2C0k1v1fCm3v/prmH3/LFblXkmkrKvB2O9/OSpqY2L8etRYdGxtUmXaXD9qXw5tYLOPxVAzIx\nIhe55oIib4/j/byEqWQd111ZXo5biw4Njatkcb2QVBPO76ZeA49WjmiW9d+mAbH9uMg1lxd5Oxxf\nwUuYStYbvOretc+MjKtsCYQE6VenTtPj+hp4K3SfahwMvlH7MZFr9tk2NiooREVOStC7gpc11eDl\ntciSY8PjigkJwtP8TmGqw7iJR2pLrdN0kPKhJh5pqzReZf04Sm6JCVk2DUjtxxbL5/wmry9QNYLe\nOV5cb1eYWVMNXuRNX2p4XOHCf8YDRGik/jK6SpZbJlt4pK0ag0faKoNX0buAZ8SEsTcNSO1n1mqZ\nvK5AderlRTVOvKypBi8q1pcaHlcoJGVXipWCGl2lVo+iF4jHegTAo7V0K/Cobxqu+6ypXKn9mMi1\nplAlx7Eu+vWocrF1zZTLi1ldkYxx9dfzrxUlK0VgOsSvlMXCVAcT8HiPsFj32cCjvmkv3WdrfsXi\nKnO0nQg4Lnzo1KNKvV1BsqZ2iKv/hrha+78dhMflB05f+JVipKBWV+lhxWOIR30T12mW8WQFXrw/\nB9Fl4Zc4jNN9tuLKyGdMYzPHqW6sn5Xx2dQ4YVha7H7C8GLlupKjx0GYDokrxUhBra6yRhrwwnma\nk5BOsxxbgRfjfl/dp21A6hfm5+3Mcaob62dl3m5qnDAsLRZXhhcr15UcHVegqxTKVCsFNbrKGmnA\nC+cprkinCcfW4AXEIqMlPNBWGTyruUqMqSQwYY6YPojiygh6TF7mOClyraC3zcvV2xlTLK4ML+ZO\nV3JoXDFdpVSmGimo0VW6pBleOM9aDnWacGwN3t66T90xVESu2XPZ2MJxpsgNT6+UoLc2cJkapwqW\npsK9gC++pQJrUkPjihFwrxQ+52J5FyYhhnh27xg/30zba93RVgHC3PsSh4luQICK/0OA0jlVzzuj\nyG3z8mq8Ycry0jyX/t4rrlxl6jb1l9NyjtpyqeP76z6nuvATdZ/Itp53aity2++dvRpvmLK8kGBn\nYq+48pWpG4SpYhwsTm7SVjl9U1NbpRSEsp77mFiVmEQTv+YUuU1efXq7Ll6CpPmxW1x5V4qx3nXg\nn6j7XFABc4rcgev62iG6gCtl2S2u+q4UIlRL9fUINRQ63uybKBtL/dB631lF7g/xYjXjJPeLK8fY\nn0P/mhr4E1f/mqb+VkdjXF3CriNjP8NFlMMBB/rr75w6mvFovEYF9P7bHg/qM4CMfD+YbMBzbM/g\nqmPDAVex8AvBfxGUZ0czHo0n2cpf67ZJlRh7jIPDRZTDAWUdbPvl7pw6mvFovLbLq7ZJlZB7xFVL\nRCmtL/w1AHC3myZ/59QBjEXdtPB2cG31NqlIekBc9Ygo0W4jofE277sabA18yCNFpXkBv6bcvbOr\nqo4evD1cq66aVzzrPwfEVY+Isk6Ezmi8zfuuBujdxJv5HaOm3L2zK3mfUj14e7i2eptU9GNAXHWI\nKNFsK6HwVshbzZvV9ks0j4yCqIo3886pivIKQa6014G3j2vwr/e9qll0bHtc9Ygol1BSeNv3XQ1G\n2y/9PVaynWMOXDzFxZt5OZCi3L+za3j5yW/ue/D2cW3tv23HqhwQVwlrN/Xo9n1XvWUp6H8tIdo5\nZcK44qLSvB1Dh450kb0OvBUrW5a4xnbUqDFuH98eVxlfrPIHNerSfVIdigIvLBjOT27tTqdKppmR\ndM1Vp6EWD6hoCNZfTUxUmnVAXVXgMpb9VRNPbhW7k2uwowZUR/f3oLgSq/xBjeqINrXcs8ZX4Bl5\nK0m0fDwdFJVF4A4/5KMheFyxReBs51RB2VQBqt58xk5c+Xi6+JxraBc9c0xRVwyubX5UPiiu3FX+\nVrRpZJXkrUw18UiiVcHTQVERrVh+REJD8LhiCwxJweHrSMEEqt4qjJ3GdqvAFJ9xDe2SZ44piitw\nDRQlrFhfclBckYiSdg1FuSctyVssH0U86prYvqso0TJ4vvJyr81Eww0B/bttpMy7CFDQ4npgw7im\nYaWNU3mDmuIzrqHdBOLXDr9kQNcI1wO33ZWOcfXX3//pKuNkhukQV4/CPqkhO1wGRlbpQKVDgEdd\nUzgMctRYWfmCquDpzsbdJonwiB/R0RBe5YfcrL8CyqKLMIrcCmOnEwE84hRSpnjbNaG2A6Alrm3u\nr/7394j3zmwuACEUKwEv5nLQyirBV/WNePEuGtSjDC9Egr9taMbRNVcZLAw/xkJD8LhiGhiKK6Ic\n/7F7uc8Ayqh6q9aAsYd4jFS4mtKtIttotu0a2uUgxhQ9QpnANaxzXrAnPWYcZHMBFldU6eWglVVW\nqCIedU0hJ+EViVYNT9dcZXJLeIx0YaQheFzB5DZkpZ1TgbLsIooJVL3VGNv+CvBEDdnibdfQLkdZ\n4tqL3A8yESVrIivaNHJP7i9LM7zSNcWThFeVj2YMXXNW0JXyER4jnRFsO7OLmok36TkPUJZdRDGR\nD8atYms1oClPgBf4/IAude02qaX2huw3IUSUOJWKFtbJRyUeU48SHry8MDLN7JduJOjd0euSMPwo\ng4aoiDfLdc0oqy5C60MrjFUcMzyxU+y36VJf4Hk7tUVK2Us/HHYPqoK1n05n40i0RGkdFLPvcSw/\nA0EGxCNuM2+mrolKhGuDT7r4CUg37LV3im2/x/HsNkyha2u3SQV3xu+PY5so2Nog0YKuCRmHiXxj\n29CUzSgvHRCG5/EzEJRfiDfZmJhzuILcaVb11rD3A7rUPI8jl/tTY+btzK4TV1skWl7X1C/RmprK\nuD5+Urxpd071uojOzVlZdYaYD7fErU1Xm671qe3QtdXbpCLz8XHlbVGL5noTs13TQsCBi99UO5qd\nU2e7poWUIdtP6FJXb5MKpEePg32XPrKoJlZ0TT7WDot1syG8xtHupn8WgCiY+AFdany0u/EzvL/a\nyOf3Ff/n7Zw6onP/E1ebI/kft3Pqhm1SsTLn4uob5ZDIaWPiD+Xp56tgJq6+Uw7ZH06u0vilKbuM\n8f9R9teAX+I7q4DWT1MqsGrH1ffKIf1Kqh91lcYvTdllHFbVwmvquq89Z0bjNW2TOJpSoUA7rlpy\nyKY1OjnsRizqcuTHVRq/EuVljKcBlEXFNPB2aA0SR1NKxVWPHFJ4Uv8x4t4C0KUidEpK49em7DEO\ni12+bsfj8xn+fDzDHf03ilh3aQ0SR1NK9lc9ckho7JnvAarJ8BaoqCY+uLGspXhNyi3GXZuhcodr\n6Y4q2KU1SLpBKRlXHXLIipP6nebcizkHRkFUFaH5ffRmyspc6MFxNaJDzj8kMWYYT6MpL8fb7Jq7\nuSqJoykl4qpHDulXsFruEd7j6zlRpRw7LBspnsBl/VwROiWl8XbKxtwIyg3G02jKHXgjXHNag8TR\nlOLz9g45JIsDkVSthOsuRKb2DwURMmMrcUVoURonrC2iWG1uCOUG42l7LUvKy/GGuOa0Bi3WohSP\nq9zeYr0+CEKtwnSRvJINuDKY6oJQ0+XxuGKKUL6i0aVsTbiUZSMFY7Vh0OKhSxqD1X2Fsb8ZqjWx\niHLk4VfBZg3rotagRbOUMnEl1usbeSXKrrQ+stSxquHK6msUhCIcNlFYFyKer8UTdPWzleVZ8ZTK\nuZTfYcUdCs98ytrcHGW2DSKS1hizjMM9IQtfqGWgTEKxZZQDDxdPl55zjeyiZ4tagxbDU8rElbte\nH+RkKLsy+sjCRNVwRS1i4MiPdlyxxV1MaexRBhMkPKtQVoynGcqExzhrDBZXPuMwrT2z8jkJlEko\ntpByKO7hmdIzrpFdRs245vhG4mhKmbhCeSUDR4VpnPilhbZGH+nLK0ufoTdDRbUmwGVby2WTIT+J\naVwFJ5gIcVr2ojGUfXMzlAkvcfYxWN2D0FMyblKmhYELKQdsr9VM6RnXyG7LNcc3EhtSysQVDtR8\njDKKUKOPTFRMZwNjUTkLXwYOTsTv1hXCWokpjZEyG6PABMVBhbI2N0OZ8BhnjeHUfcjNGNN0iI09\nRsNqVajZpjGnpleFmXG47ZoUqIF31hbNSqA1oKONwzF0wzqucKDGIS9a4LOB+FLT6iMLD0Wj0vNq\nOHAifSuIcIw8YcIa6q+QMh+jyETpr2qUtblZytD/MdIaY44xTYf42FMok1BsKWU1vSrEbOm2a2SX\neda8yqE1qJeilI4rHKhxyItWqBGz7MrqIwsXVcOVmaKG4440PWHzdsJAyrwvodMlDmqUFeNplnJf\nXLmMaTrEx55CmYRiSylXplfpZfahvKuIVdx2jezy5tDVw69y8C1HZCxFKR1XKIeMHmP4kYKzKEIX\nyisdjVY0b+DiQfi0PGHqF1IaI2UeV2QC4qBCWZubpQx4wDd8awzWX7mMUXQqxh6jYa2JWI05xAtc\nSMRqHG67lmMiamf5x9qi0QN8o7tASom4EnLIgE7reFDBWV4FLJNXotyfM41pDcfPa08qitDyqERS\nZm2OJuBYhbI2B507pyQoAx7LoDDajCdGWY49WsNqVajZpjLH8ISI1Tjcds0XqClb7uaq9CiRpVr7\n0V9V8IaACO/f3/MTXlavlNQ0Zt4ceHAagsDFNqN0abAMTpvj/SDLxpPG3AxlD89goAH+iNth7I49\nfGxEIJ6omwtj3e18qe5K03bNtduwhb7RU3ZKif6Kkw9p22/Oyq60vLL9ptOF0xCMlFCEenc3bly1\nhWfGXJty6MQtnsFAzjOM3bEnPvE/p5fJCCMTdXNhhhMW4FSv+7Zrrt2GLfSNxNGUasSVw6JfdqXU\ndrKCOuGkpApGdw5p42qF8KxJuQ9vjrE/9mwQis2IWJuu9WlYJ/SNxNGUEusZeAvND3kyd+3XyJVk\nslas0jhc5mY/7hqvxvGBlGcZu2NPg9vcqRkR60DXSEJO4mhK1ePKHaPm3HLO77DytVgxSmN3jHIo\nzR0aRlmvQ7aM3bFnjl/j/JyIdZhrbFU4iaMpVY+rzjGq4etep7AjBgN9YxSU+sZvw3jL/234Rt5N\nUySOphSsb3/GT2uq2ET+qZO/T2n8+xjPti0NrZC6pmhK+6W+xU/8F02/6/P7lMa/j/FcRJA4GlKP\nFE3mEdUc0J/zf2pgUQ38H2szOb82cVvTAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}- T_{0}^{4} \\epsilon_{0} \\sigma + T_{s}^{4} \\sigma - \\left(- \\epsilon_{0} + 1\\right) \\left(T_{1}^{4} \\epsilon_{1} \\sigma + T_{2}^{4} \\epsilon_{2} \\sigma \\left(- \\epsilon_{1} + 1\\right)\\right)\\\\T_{0}^{4} \\epsilon_{0} \\sigma - T_{1}^{4} \\epsilon_{1} \\sigma - T_{2}^{4} \\epsilon_{2} \\sigma \\left(- \\epsilon_{1} + 1\\right) + T_{s}^{4} \\sigma \\left(- \\epsilon_{0} + 1\\right)\\\\T_{1}^{4} \\epsilon_{1} \\sigma - T_{2}^{4} \\epsilon_{2} \\sigma + \\left(- \\epsilon_{1} + 1\\right) \\left(T_{0}^{4} \\epsilon_{0} \\sigma + T_{s}^{4} \\sigma \\left(- \\epsilon_{0} + 1\\right)\\right)\\\\T_{2}^{4} \\epsilon_{2} \\sigma + \\left(- \\epsilon_{2} + 1\\right) \\left(T_{1}^{4} \\epsilon_{1} \\sigma + \\left(- \\epsilon_{1} + 1\\right) \\left(T_{0}^{4} \\epsilon_{0} \\sigma + T_{s}^{4} \\sigma \\left(- \\epsilon_{0} + 1\\right)\\right)\\right)\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡         4           4               ⎛  4          4               ⎞     ⎤\n",
       "⎢     - T₀ ⋅ε₀⋅σ + T_s ⋅σ - (-ε₀ + 1)⋅⎝T₁ ⋅ε₁⋅σ + T₂ ⋅ε₂⋅σ⋅(-ε₁ + 1)⎠     ⎥\n",
       "⎢                                                                         ⎥\n",
       "⎢         4          4          4                     4                   ⎥\n",
       "⎢       T₀ ⋅ε₀⋅σ - T₁ ⋅ε₁⋅σ - T₂ ⋅ε₂⋅σ⋅(-ε₁ + 1) + T_s ⋅σ⋅(-ε₀ + 1)       ⎥\n",
       "⎢                                                                         ⎥\n",
       "⎢        4          4                  ⎛  4           4            ⎞      ⎥\n",
       "⎢      T₁ ⋅ε₁⋅σ - T₂ ⋅ε₂⋅σ + (-ε₁ + 1)⋅⎝T₀ ⋅ε₀⋅σ + T_s ⋅σ⋅(-ε₀ + 1)⎠      ⎥\n",
       "⎢                                                                         ⎥\n",
       "⎢  4                  ⎛  4                  ⎛  4           4            ⎞⎞⎥\n",
       "⎣T₂ ⋅ε₂⋅σ + (-ε₂ + 1)⋅⎝T₁ ⋅ε₁⋅σ + (-ε₁ + 1)⋅⎝T₀ ⋅ε₀⋅σ + T_s ⋅σ⋅(-ε₀ + 1)⎠⎠⎦"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Net flux, positive up\n",
    "F = U - D\n",
    "F"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABAQAAABNCAMAAAAGoUZqAAAAPFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAo1xBWAAAAE3RSTlMA\nMquZdlQQQOkwRM3dIruJ72ZsPtw6XAAAAAlwSFlzAAAOxAAADsQBlSsOGwAAGo9JREFUeAHtXeti\nwrwNDVC6rdBSxvu/63yJLFmSZTuxof1WfoAT7KOjSxQnBGtZnvr6Giztp+MNVhfhRiuOyLtaP4jW\n4UPTZDjB4YAa68K+t8L+7t2HR3gduwduGXB+vG8ZVhzz0/GKxPd+MVrxvXzW8T+I1kE9OocTHA7Y\n44mPnUftNR77y3J4nN7c69wjfGvf99vYJPDT8bbaqTputOJVgW0dfhKti3a6GU5wOGCboaHX7Q6t\nTZ93f+QfHz4JHDYBbBn09j42CbwU73OfA9B+909sN7YMxYfRWvp5GbQKig1jy8me1KlyP8ECb9jd\nBThB2auW6oBc2+d9RhL4PJ2+vnxi+Xx8fx2Pj4d7uzxubsfhvC8JcMCX4n2q0802u7Nex1oW6FB8\nIK2lxquDFtMYNgeyzcl+hrPNfoJANH5yvL6InqHs2ylnuGFrRhJ499cpt4e7K3PzB/7yffXvJ3/a\nvC37kgAHfCnepXbkerXt1/k75vH3i91PWNJQfACtpZXXfn8MYKuTvYWDYz/B3DEcry+iZyj76U7j\nO18zksAxBPb397KEuxafj+COu0sKbsbWmQTO+ZyOAb4U796dgpky78fTES6ObrmewqvtivfTWnp4\n5X3baQmN4o5+tjmBpWjE73CPazdBZhyG1xfRc5S97r4onZEE4pn/+DgcAr37I7jj7X15d5G+Kwkw\nwNfiXbsnAix+3WHwBkng3aVM49WheD8tFueeRplXpkMHrYJy/WwzAgFVJXsIlt1PMDcOw+uM6DnK\nnnZflc5IApdgf58Ego/Ws52bCJxP7vUI1wWFmBC7c5czwJfifYaLHEHY2pEr43um+F0u6m/agNau\n+AZaeZwHkWVemQ7ttEAR9rmBbUagTPYesupugi5o6RyN4fVF9CRl3+zzB7O4tjkjCUQ5VzjJxYkB\nCI/3a9yM4HY85vOY99N3+NEym2VnLlgxMsCX4RVvyEjNQHepDB5sp0xtGJF/6ornZttAK4/zILLM\nS+qw3vNZqYI/FtWduT4uCZa03m3EE7nLotptRwRmeO7u93pfpxbRk5T9gAONG7d5e1oS+Ai/Bjge\n6y2BldH963H0Vwcf3/HM93Z7CzcP3Q8H19PH+fr1mU+ylZjLADnecgDnj8G7nb5O0ccc7ytPYcni\noBkySV9JZfBguwPv1Fs0VMV7aaHBEb+Hl+yr0iq4E2XGVs2ICltJQDXiF06SdYLST9yUkaKUl+Et\nPAKfHjH7f+GflgSu4RaKs+N6SyB3/2e8T7DcXQY4h57v4UaOp5O9pAtUQMA7nI7rLH0M3s2lpGM4\nPAXeVX+4SjJBfaQyGL/n+tWFZsleWmhwpGXOBDgvqYNGaxG8iDjSrBhRYysJqEa8QPxVInBDxKgK\ng9+fHzGf/ne4Xa8RSeD9ck2vS7xqWnB2m26A06x+fJxu/lmnq+cfpjNx+kamNreLf12/wwd61B2U\nafqTzvpuX8RzKWU9mDbgYQ5PeP5WzuHhJycCb53KcOUlE++fkjIYv4eU/jiiHx9eqDjs6aeFBg8Y\nPbxKfTVa0ly6VhUj5mxLBDQjLle80ECCGDDJwxsiJuGRgEl4z48YSD8YFL2tEUlAkUkuf+ACKsvq\n31/v726W/R6eVQxzgEeYXb/x86HM++kSFHO4exQh4jkmkAT68TCHI56PwphpBd4av1x5yQR7SGXU\n+MUBeQssSfd20iIGJyg9vGRfjdYieBFxpGkbUWUrCahGJDMBIEgDZk/EAB4NmIT3/Ij5GTMB4tS1\nefZpOF7dpwuoLKvHRwfWk6z/ofMzPrp8wfwdoaTLEyAe8G4ykcatSWADHuZwiuek+JmHxCvMZAUT\nYh2pDMbvB09/ZFxsEsXTd7204qyG/7Lcw0v01Wgp5kqUs4ZtRJWtIEB+YiFGxCRACMIZwnEQfpKm\njEyFvIRHA4bgPTtiopUys3ZuTJkJHMIh+RZu8cEFVJ7VVxfFH3Ov6fGBw/oAHSohXEBvCaBP0eXr\nvvg4Qhce5vBlQTx3yeJ/I5J4hXtaOBLZgTpSGUwC9RuDYElA85+9tNDgFKWHl+ir0VJ4UXnYto2o\nshUESBIgRjymG4OEIHGJ8JM0ZaQp5CU8GjAvjBhyCY2G7WrNSAKf19vtBnfovtYr+Dyrf4Sj/T3u\n/PZH2Ze7z/YZZgsZf+GCBQBdN/Tpikf2bcFbc7hDQbwlPn4q+a37M7bZSGQHfaQymARO648k0Fd+\nEsVv6T+fQk2bFjE4EdDDS/RVaUlzEXGkuYGtIECSADEi/npOCBKXoIdhpzBl5CnkETwMmBdGTHwi\nghi1uzkjCVzif5Td/Pb25dqXoz/IWVb/OL7d39Z7AuHu5rv7O/Mt/hhHlWAuIICuF7jPNSMe2bcF\nbz3re/kJDx71FHhk6ukHpFcaSditXzJlltPx21kn/tQId1QTTt7IFD99LPBf+U5acT7Gbyf38Mr7\nlmi5pwRUd+ZKOTMXLoJEeODAnIC7NaobcT0/ZgRpwKCHwU8Fyrm8HI8EDOI9O2J+5hOD6DDS0s9B\n2Y0C0js1cxek3bEB7qO7tX3kexMPzvqk/4c7QX+sv3iQ3b7pLwmtl2RiCO96muzr6+xvq+ovm5Zq\n8B5eRl+Tlk62YkSNrUEgM6J2z1G6xNFSdyJdQ547y6V7UWnA0yPm4s+xu14zZgIqIf0c5H9FiM8J\nqIPcRa91pGnuO9sPUZp4kMORy+FyPp+/9ANO9sZxviWZGMJrfyDKoA/u39l6YnLdbFqqwXt4GX1N\nWpkCuNHP1iCQGZFcGiRxWsAofkr9fcOQp9n66RET72VklHs3npYEyCMBlOPtdFufyKN7G9vSp5+n\n6yNcfTQiZN2UHB4fZM564caXlaD6mFT/SoxSXXZxcg0lTVru3NVj8B5eFVpUA9IeyDYne1DOBjJg\nlj4/EeKuqQTM8vSIudNnaHJ+rVvPSwLqOaiVpt5Pnmv1fm17rRyuIgxcIsI+FJj0q5sFnIozgWUg\nraWHV4UWUwI2B7JlZMPj6SAmfr44YEa6Jinrzb7z9bwk0HkOquu1J4dr6GYO1wYsExaLUuWwnff7\n/c16UnQYrb7lxWq0mBawOYwtX14snxg4ca8PmAkR45/B2/t6YhLYS/Vv/J8FuiwQHlnrGvH7Og+Z\nSP0lgd/n+D/GjRY4W/OlRowf3m37DTWi2OQkkB7bIiIHNUdDj8YbpKaEGU10NJ5k/FP3jNZ8NN6z\n7DY3CUysyzAaejTeNAeOJjoab5rie4DVWkSjNR+NZym8+9kACj41CUysyzAaejQetfHQ9miio/GG\nKjsKDJ6vzPBGaz4aLyPLN/ZWHsrwpiaBrroMGa3qxmjoLrxpN7SrarsFuQY8HELFGHjD1Gz7lWGi\nOLUWkaE5NVBz28Ibplv6CWRn5aFMq5FJgBdm2FsZxP2m8yOLmAy5Ixu9kNfLyDwTN0bbtANvoJrV\nIiZO2YniQi2iDs0VP8hdPXgDdUsRM6DyECg1MgnwwgxGgQwQb37+1CImM0pIlAwx2qYdeAPUbC5i\n4tSfJy6uBNqheckZ2f4evAG6iQorZN2ejNeGjZFJIDy/uC654EqN9FcGYfxHFjHh0H57I1X7cXcm\nKGyy/6AU62X4zqzvbptuxutXk3M39eR2migu/svndZbU/mLA1WfbzGuaJQdUHgKhA5MAK8ywoTII\nkFo/4zJOY4qY5NC7qNp/08sFxS3mUrczrSMgio7kfXcRDcI34/WryZOAqWe0DL5PFBeWr3uhJSt/\nlUQbYCv3mt8vI4avD4Wje1sjk0CQva7C+NFXl0GlPbKISS4gPm69jWr2f9UctrRluFQUHcn77iEa\n2WzF26CmlQSEnsxWE8XFlSxeZ0m3Vg5Ttr6Ze833T0kgWXL/OgLAY2ASiJCwCmPYaq7LAHSUzzFF\nTBTgtGJpH9XitdjuehmBh3Q/VHRdVei0aTvec4uYMIdssKqZc3C164WuvLMnOrdaclaFFVw7idmy\ne3N0EkirMHomvC4DLjqu13nQ2DcWMcFSH83QJlVcTZrjFZbFS+VUkElSR4YP5nWyLF7oL/vCvYsI\nx21aJhr7t+L1qom+TGqaRyXXE0eF1gartoojtYhMl6NG3Bb7LLnUdEPBySjSazJi9q8tCOJGJ4G0\nCiMIcJ+wMDouOt5YmsJjtBUxwdWk26EtqriatMArLJALSiITtIDlUv4Xd9mXLq2aMEGcQTT2bcTr\nVRN9mSiJm5ruGwxdricZ5psbrGomASKO1CKyXI4aCVtErlstWdMNBaNNpCxpyf2Vh0De6CSQCjOQ\n9Ab1OHCxKFHJQy9N4UjitC4h01MtQLsAXC+82qETIJ5MX1hColygxJFalzbqI1qq1YF4EATezqHM\nEjm5dJUF2V9cpSIO/Rs4l1TDIwVrudB1gJMl+6KzJK7NkktFNzwsvG4lWVK1/fUGgi3dm08C//r3\nf2Bz9ydcdNH01lxpREonl4qAnJ1qZamPxqoXThQA4sn0JxQdkeeArUSjNQ08am5hNm2JPm80UTWG\noEhZMnRJd9q0xbmekOTJoEZxpBYRuHxTdEpxyTOElFJ2xdYtX4t/RZKypCVhMkiFb2v/999+IeD9\ni5Os0tNFF01vzZVGhApaEZMsHppLSAhovNL2P07FZVBpSQon5iVFR6T70aZ9RKPKBh6xiay8UZif\nC18SEHN+XlpXeB1vi3OdOpMAEYcVBtCSqfpdcnmMgM11cJIZpCVLlwOWKaXXMAmAamNnAgOTAFx0\nZelt9QIuOl6q85AsuTa0IiZ5PKCD1yBphSZlKn2+AoMi3rqatMQr3ObBkTJcLZfyG2ayL9g0rNLY\nQTQa0cAj5u5VE31JQMwkwPWk41x7g1VbxaFnwJLbonOrJSu6qaaUsjAJgCXhxMUsuWFz8D0BKMyA\nedVxaq40wvirRUxcH3KUNZeQYNBuE6iGb+BfoIgHq0mLkhTKKtMeAkcSegHbERZ//ESX8kVxZd+t\nRKNwE++FRUxW08DHBqu2mhVrEYElt0XnVktCJIGq8CkPC/imLWLIDRwcuKk1MgmQwgx5ehOlJAp1\nHpgGahET14ceZaLURxs0oRqEhlJjvpXw4DFWgQezsTCOvKWRlF78nodPoV5G6Mz67iCqys7wXlnE\nhFjONzdYVSSBglnX39OJ5tui0/KMZcmKbupa/EyWWmGFPv/ArNm5OTIJENEx14b6YriX3ifAvf0t\neZQ5DHVnK7Y8D2mrSa9otQdcJRPuUkJLPE1m9HXDuogGMRaeWS3EVlP1pSFL6ElsEJq2ONW/jeLk\nGXNbdBriFtOSlceGNVMaspIlf/ATg9G5anrzt/rNSiM8MArb8ihzHfesJg1nfZRnLT8ue+M435JM\njPoVWb0MP9roqxUWsYh6OBPPrBZiq6n60uAu9AzkyJstTrOqpVomTtyd3xadhnaLacnKH4g0Uxqy\nkmr7Kw+B+SfNBPRSI32FL4Ci+FSSwK7VpJWzvrn8eFryXTBzO/qYiGWxNci0r5doGqg3KtVCTDU7\nF5Bv0NMU12fVJRfHb7u4e77k14Fkmx3RWbFkpXxDl+CkWryXm9jvacxKAlp628OTjpWnWvptd7t6\nMuWIA5eIMCOfy+0mygHYdqVayEA1K0dB4DVPnKxFNDo6K5YcuWBKipgBlYcgHmYlgc4zBdBp+Ow8\nKdQRzbO+Onz8YlGqGL6znyhHyLdr1UKGqfnq5cVkLaKuk29uNW2rZskJRUeWAZWHQJVpSQAE/H3+\nWeDVFkhT6FcTGSjfP/Q46vWXBEZZ8g/n51rgn1eLaOTFU/jvgPXE4Lx6CvOQpwXjaMqj8aYpPhH4\nOTYYWYtoNONNeEMqD4FbKzMBeJIOuo/7nIc8gKNaqyL+mWAAOkA80wTimUUg8bTP59h0qDoq5dFe\nG41nWQCjAFuuv50E5tVTmIdsGaHxO7VWxTKa8mg8U7mhxSpMSYUvn2PTgvBtu1XKo702Gs9UFaMA\nW26AnQSsegqmtOqXA5CH3bxOBR2AtVqr4leW/UDV2opV/EKbPpnygMCFMAufFt541TAKsMWSQE89\nhUyT0kaxesjuuiTuqZxN11Iq01TQIX4balW4PxZ+fx2Pj4d7uzzcrdi9lDvwpqjWUqxiiuBgVLX+\nx26bOugnU35iFExRDaMAW/lMoKeegno0sZ3l6iF765I4QTMKOgT+60qeo43RgTdFNf+ATO01RbAX\nqtf/WH5+GHR4rWbc8H0P3gxvYBRgK08Cu+szMDsUq4f01yUR/6ioPWvOqPjNHEMr6BAGrX/S2W2M\nXNzSjjdHtYZiFf2CmY4VmzITbChPw+Rp/6ZQHJ/tYhB9lPsDl0UdM4GFN8UbGAXYypLA/voMmbVh\nAS9ZPWRDXRLmOvf4t1tmp/MlMNK/+rMSIKFWxbLfGLm4DrxJqtWLVfQLznX0/jBsykywvCYM9lDe\nwDhPAswEJt4cb2AUYIveGNxfn4EdlaXqIeeTez1Od9bd2uSuS3+otAax7zgGBmwq6OBGxD+bL/uN\nkYtrx5ulWvWvpxsE5zp6e6ckIG3KTLClPA2X92TKGwKXJYEQkut6p3aBnkmqYRRgiyaBQDBfPLGz\n0kVEYO969ZBOZO59ckWTiyuXAMnd4UdhwOKixlmtil3G4JS9RFjp0rfXC2X3b98JZT8U1arFKjbY\nVOqoCHaT9u+gcHjLTIA2uB2P+UmBWSUCcHkvoAyB634yZpRVxjLqClEwIAy4dbQIxyjAlkgCaTFG\nb3Ve6SJ6wr/r9Rnwe2zp1UN6kbl+hTXpsAQIEkgtjkGSAKzb5vqSWhUuRMndtL2UHbaKx005SzW5\nuEayTGxsECxMWrVpZoLuAOPynk8ZokBGGvfjal1OuRAF4oiapBpGAbZEEoDFGGmA4OLGh7A8vTt5\nhQtnP7j6MquHJOS0pH4BmZuysDptwqPFCYAixyABS5YoILUq3B2TMwxOn0lEWru+kbJD0PDE6Jpq\nSW6ixO95ui/whJxUs54OD1A1wclHKFiYVBPs1nPEn3M1E7jkCHauBBiXV6OshAGHUG3VQRncIfy4\nWknKU00ghtdUA7noDGXSIcMAowBbIgmk+gwEHEp84JL/olZFV/UQhAZkXPtfIOvVGNbVYrhYwEOm\nQZaOQQJWr1WhFpcAEbh2fSNlx0QzrhhdUQ3ldqoGCxW7YdxsAalWJAN95LsXTFq1qWYCZ5fTzT+9\ngm4TVtHlVWyFeNMoJ3cIxkUTqSYQwyuqJbnBd7p1aH6DCMcowJZIAnjJBkk5VZxw4uC0ImpVBCrK\nG7loS8jkjALFQ/yd0LjyWwGZ59PVRFwi4CFT0oNj6AFLalXgJTxJuyACV6dppOyYJBMQVmJ0RTWU\nS0BaVMPTLRlImzXByUc4SMit2hRNgPEFJU0cbiXAuLwKZcQbQplEQaKc3CH8uIrklPUokDVLKqol\nuaiZORPAJACTLowHngTSJRtNonhlvPpIqbBAuJCmVj2EnlEA2asUMlMJmZuyMFkCPM370kY4WSKL\n3eIy9XgJT9PuKgLXrm+lzC4GVyvJ0bZqKJeYuUk1zPx0JGnbgrOyBzCKu4Wee1SbqvG1CLdJq0SJ\nXF6FshYGHKKdMo0CoJzcUWIsXZNMAEb0n3K4rVqSS0FaVMMowBZPAuR6BZKyezgvPOLgxa37ZK0K\nygXbavUQPOtTZIftK/6UkLl+hdsmgilycQKyv065bzAJkBuDCIGX8DTtrt/j2vWtlAu3BEKhwcP3\nWm7QsbJVQ7m9qsWRdBRr24JDZ/6XN2HSmk3V+GoPMC6vShmDGJTlENUwSJRpFADl5I5SFMioS3jA\nyH/K4bZqSS4FaVENowBbPAlAfQaHjfaThTVESQ7KJbX16iHprO/6IfJa8WcpIHP95NrbQSriIftE\nh2OQgCVrUWKtilSfJEu7qWaEP2xDvYJGyqzeCdASo23V8jXzAaRFNbwbDKPYpy04dE71GdahQm7N\npmp8kTAAtwmrRIFcXpUy4KGqHIImATUMgHIWBUAZ3VFgLJMA4CEl3xLDbdVQLoFpUQ2jAFtZEiD1\nGRw0sZ8orCFKchAq2CxVD1nP+r5jQoYl9QvIXD8y10R5FI+whw4Mw65VsRBj5Gk3Uo77QmWFNsoE\nzxHC2j9itK0akQt6uc8W1eiP9WQoNm3Bvp8ITCZXLZLhxsn6H24n9VAKA9gprOLFCz0rdT38gLVY\ndRgd3jZTzqNgjVx0R4Exc00pCtxTAm9vN5wO1sqxoFzUjIlyv3Yfvx+X9fGLS3xUizyygfGQJQGC\n55rSfvq+fFTbFj+juKky+UVeweCuqz42rLAXGEkOfToLE2T6Wku78bzw8GtXF15lcc49H4v6b/UA\nZT4vqssty0LV8AmxAuGqTaWPynIXFOwCGoqrU8mKh6oBJuSZtnLSpBABgZxsyloUuGm8P7g2hoEZ\nBbY3VLktqmEUYOs1SUCcURZlSX10jmuJagwyILP+ivclRhqRCjr4PfK2rJp2syvEhEQagjL5zqxY\nY6umyi3LQtXqxSpswZqPynLdXIfoK22qHZ9ugDxoCcqQMNhMWY0CvYQB5WzIM6OgclrUwsAQlbyB\nUYCt3iRAnv+kmva1ZbBtWFI/LcCuy+4pTpCvRksuDFdoNe36Hz8311OyK9aYqvXJRdXirSfdWLDX\nFNznIxTswaVNC8d7j9sCa5OyE9IRsBXKahQsfe4AQ8dPOwrsYg19cpNqGAXYyv5FmDNUkvKgJf+V\ns/6GJfXN1ST6mOaBJGtV6Ol+x9r1lYo1pmp9FR1QtZZiFabgPh+hYB9Wik2V+Oos3xTi1aQ8Ngy0\nk6+/T3LbuOpnJQoqi4p0yU3ewCjAlpkEOpIoyx/mZt8ZpQw1fvGlVZasVdGXdsuU4ZtqxZr8vzQw\nrP8TlxdrK1bxRJt2naQNzZ9H+dlRMLdkCSleUrwc6EuihpPEV31nFDF8/o40d0JRXWkXh5Va1Yo1\npYHb9/snXV74EjadF1+jtBSUd5z0VU4vjQIaDzEJPPwrPRGkMv5/2vnPq1VRmVo+wbm/0Ka/kHLF\nkXjtBK1rOPTdHwHf38LL+LWrgv2P+3pkrYqfYZyNl60Dyf9Cm/5CyrbDMAqgdY/Hvj3s79s/C/xZ\n4J9vgf8B2rWZDp5L7TAAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}- 2 T_{0}^{4} \\epsilon_{0} \\sigma + T_{1}^{4} \\epsilon_{1} \\sigma + T_{2}^{4} \\epsilon_{2} \\sigma \\left(- \\epsilon_{1} + 1\\right) - T_{s}^{4} \\sigma \\left(- \\epsilon_{0} + 1\\right) + T_{s}^{4} \\sigma - \\left(- \\epsilon_{0} + 1\\right) \\left(T_{1}^{4} \\epsilon_{1} \\sigma + T_{2}^{4} \\epsilon_{2} \\sigma \\left(- \\epsilon_{1} + 1\\right)\\right)\\\\T_{0}^{4} \\epsilon_{0} \\sigma - 2 T_{1}^{4} \\epsilon_{1} \\sigma - T_{2}^{4} \\epsilon_{2} \\sigma \\left(- \\epsilon_{1} + 1\\right) + T_{2}^{4} \\epsilon_{2} \\sigma + T_{s}^{4} \\sigma \\left(- \\epsilon_{0} + 1\\right) - \\left(- \\epsilon_{1} + 1\\right) \\left(T_{0}^{4} \\epsilon_{0} \\sigma + T_{s}^{4} \\sigma \\left(- \\epsilon_{0} + 1\\right)\\right)\\\\T_{1}^{4} \\epsilon_{1} \\sigma - 2 T_{2}^{4} \\epsilon_{2} \\sigma + \\left(- \\epsilon_{1} + 1\\right) \\left(T_{0}^{4} \\epsilon_{0} \\sigma + T_{s}^{4} \\sigma \\left(- \\epsilon_{0} + 1\\right)\\right) - \\left(- \\epsilon_{2} + 1\\right) \\left(T_{1}^{4} \\epsilon_{1} \\sigma + \\left(- \\epsilon_{1} + 1\\right) \\left(T_{0}^{4} \\epsilon_{0} \\sigma + T_{s}^{4} \\sigma \\left(- \\epsilon_{0} + 1\\right)\\right)\\right)\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡            4          4          4                     4                  4 \n",
       "⎢      - 2⋅T₀ ⋅ε₀⋅σ + T₁ ⋅ε₁⋅σ + T₂ ⋅ε₂⋅σ⋅(-ε₁ + 1) - T_s ⋅σ⋅(-ε₀ + 1) + T_s ⋅\n",
       "⎢                                                                             \n",
       "⎢         4            4          4                    4           4          \n",
       "⎢       T₀ ⋅ε₀⋅σ - 2⋅T₁ ⋅ε₁⋅σ - T₂ ⋅ε₂⋅σ⋅(-ε₁ + 1) + T₂ ⋅ε₂⋅σ + T_s ⋅σ⋅(-ε₀ + \n",
       "⎢                                                                             \n",
       "⎢  4            4                  ⎛  4           4            ⎞             ⎛\n",
       "⎣T₁ ⋅ε₁⋅σ - 2⋅T₂ ⋅ε₂⋅σ + (-ε₁ + 1)⋅⎝T₀ ⋅ε₀⋅σ + T_s ⋅σ⋅(-ε₀ + 1)⎠ - (-ε₂ + 1)⋅⎝\n",
       "\n",
       "              ⎛  4          4               ⎞      ⎤\n",
       "σ - (-ε₀ + 1)⋅⎝T₁ ⋅ε₁⋅σ + T₂ ⋅ε₂⋅σ⋅(-ε₁ + 1)⎠      ⎥\n",
       "                                                   ⎥\n",
       "               ⎛  4           4            ⎞       ⎥\n",
       "1) - (-ε₁ + 1)⋅⎝T₀ ⋅ε₀⋅σ + T_s ⋅σ⋅(-ε₀ + 1)⎠       ⎥\n",
       "                                                   ⎥\n",
       "  4                  ⎛  4           4            ⎞⎞⎥\n",
       "T₁ ⋅ε₁⋅σ + (-ε₁ + 1)⋅⎝T₀ ⋅ε₀⋅σ + T_s ⋅σ⋅(-ε₀ + 1)⎠⎠⎦"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# The absorption is then simply the flux convergence in each layer\n",
    "\n",
    "# define a vector of absorbed radiation -- same size as emissions\n",
    "A = E.copy()\n",
    "\n",
    "# Get the convergence\n",
    "for n in range(3):\n",
    "    A[n] = -(F[n+1]-F[n])\n",
    "\n",
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAABoAAABLCAMAAABZRmeuAAAAPFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAo1xBWAAAAE3RSTlMA\nMquZdlQQQOkwRIlmzd0i77ts7uXj/QAAAAlwSFlzAAAOxAAADsQBlSsOGwAAAMtJREFUOBHtVcsS\nwiAMDOVhBVqK+f9/lRClZETPztBcSndZCLAQUFhjgRY7IwAKtSlhGwOR/hckSp3o2YqCcjrE8O4n\nqa0M6vbESkHFg8CwDahck7ToKidUWCmPnGtPOdTU3aP5UCUMhCn+QK9KrBpRPwYETsMO0oCcaa44\nSp6XrEdLhp026nhtYp9hgUM5i/H20lQtpKrB1Jiemsait/XOJy8t+li7u3xZVDxff3JTprFoqynS\nouKNuixaLfq1ZDsq0Mb407u1ZBsDT+aoGdQmGjgBAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}0\\\\0\\\\0\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡0⎤\n",
       "⎢ ⎥\n",
       "⎢0⎥\n",
       "⎢ ⎥\n",
       "⎣0⎦"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# this should reduce to zero if I did it right\n",
    "sympy.simplify(A - sympy.Matrix([R_0, R_1, R_2]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## So that works. I can formulate the tests against numerical code this way\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [default]",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}