ch2-ex1-SolarFluxIntensity.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Flux, Intensity, Wavelength, and Numerical integration\n",
    "In this exercise, we will examine the spectral distribution of solar radiation. We will use numerical methods to compute definite integrals and to convert between spectral intensity, intensity, spectral flux density, and flux density."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import modules that we will need\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.integrate import quad"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerical integration \n",
    "### Background\n",
    "Numerical integration is the process of computing the value of a definite integral\n",
    "within some acceptable error.\n",
    "\\begin{equation}\n",
    "\\int_{a}^{b}f(x)\\, dx\n",
    "\\end{equation}\n",
    "\n",
    "Algorithms for numerical integration divide the interval $[a,b]$ into small slices then add up the area the slices (Riemann sum), as depicted below. \n",
    "<center><img src=\"img/RiemannSum.png\"></center>\n",
    "\n",
    "Unlike analytical integration, numerical integration does not require finding the antiderivative (or indefinite integral) function $ F(x) = \\int f(x)\\,dx$. \n",
    "This means we can compute the integral of complex functions whose antiderivatives are unknown or hard to determine. The disadvantage of numerical integration compared to analytical integration is that the results always have some error, meaning that they differ from the true \n",
    "analytic solution. Some numerical integration algorithms, however, provide an estimate of the error with the integral value.  \n",
    "\n",
    "### Example\n",
    "As an example, let's integrate the function $f(x) = x^3$ over the interval $[1,2]$:\n",
    "\\begin{equation}\n",
    "\\int_{1}^{2}f(x)\\, dx\n",
    "\\end{equation}\n",
    "For simple integrands like this, we could compute the solution, which is 3.75 here. We will compare this to the numerical result below. \n",
    "\n",
    "The first step in numerical integration is to define a Python function for the integrand. The function should have one input argument and one output argument."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(x):\n",
    "    '''Integrand function'''\n",
    "    return x**3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `quad` function (short for 'quadrature') from the SciPy package can evaluate one-dimensional integrands of the form $f(x)$. \n",
    "`quad` was already imported at the top of this program. It's syntax is\n",
    "\n",
    "`value, error = quad(f,a,b)`\n",
    "\n",
    "where `f` is the function to be integrated, `a` is the lower limit of integration, and `b` is the upper limit. `quad` returns two numbers: the value of the integral and an estimate of its absolute error. The true value of the integral should be in the interval `[value–error, value+error]`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Evaluate the integral on the interval [1,2]\n",
    "value, error = quad(f, 1, 2)\n",
    "\n",
    "# Print the result\n",
    "print(f'The numerical value is {value:.16f} ± {error:.16f}')\n",
    "print(f'The true solution is 3.75. The numerical solution has an error of {value-3.75:g}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On my computer, the numerical result from `quad` differs from the true answer in the 16th decimal place. This is within the error estimate from `quad`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Spectral Intensity\n",
    "Above Earth's atmosphere, the spectral intensity of sunlight is \n",
    "\\begin{equation}\n",
    "I_\\lambda = \\frac{a_1}{\\lambda^5 (e^{a_2/\\lambda} - 1)}\n",
    "\\end{equation}\n",
    "with constants $a_1 = 1.191 \\times 10^8$ ${\\rm W\\, m^{-2}\\, µm^4\\, sr^{-1}}$ and $a_2 = 2.493$ µm. When using this equation, the wavelength ($\\lambda$) should be specified in units of µm. The spectral intensity ($I_\\lambda$) will then have values of ${\\rm W\\, m^{-2}\\, µm^{-1}\\, sr^{-1}}$.\n",
    "\n",
    "#### DO...\n",
    "1. Write a Python function called `spectral_intensity` for the spectral intensity of sunlight. The function should have one input argument (wavelength) and one output argument (spectral intensity from the equation above).\n",
    "2. Plot the intensity as a function of wavelength over the range 0.1 to 4 µm. If you have done things correctly, your figure should look like this\n",
    "\n",
    "<center><img src=\"img/SolarIntensity.png\" width=\"300\"></center>\n",
    "\n",
    "*Turn in:* Your Python function and plotting commands."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Broadband Intensity\n",
    "The broadband intensity is \n",
    "\\begin{equation}\n",
    "I(\\lambda_1,\\lambda_2)=\\int_{\\lambda_1}^{\\lambda_2}I_\\lambda d\\lambda\n",
    "\\end{equation}\n",
    "\n",
    "#### DO...\n",
    "1. Use the `quad` to integrate your function `spectral_intensity` over the wavelength range $\\lambda_1=$ 0.1 µm to $\\lambda_2=$ 10 µm. In other words, compute the value of $I{\\rm(0.1\\,µm, 10\\,µm)}$. This is the broadband intensity of sunlight. \n",
    "2. Compare your value from 1 to the value of solar intensity that we calculated in class ($I_s=2.0\\times 10^7\\, {\\rm W\\, m^{-2}\\, sr^{-1}}$).\n",
    "3. Change your `quad` command to integrate over the range 0.1 to 20 µm. Does your answer change? What does this mean about the intensity of sunlight in the wavelength range 10-20 µm?\n",
    "\n",
    "*Turn in:* Value from 1, Answers to 2, 3. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Spectral Flux Density\n",
    "Since the sun has small apparent size ($\\omega_s = 6.8 \\times 10^-5$ sr), as seen from Earth, we can compute the radiation flux density along the solar beam using the simple equation (integration over solid angle not required)\n",
    "\\begin{equation}\n",
    "F_s=I_s\\omega_s.\n",
    "\\end{equation} \n",
    "The flux density and intensity in this expression are broadband quantities, but the spectral flux density and spectral have an identical relationship\n",
    "\\begin{equation}\n",
    "F_{\\lambda,s}=I_{\\lambda,s}\\omega_s.\n",
    "\\end{equation}\n",
    "\n",
    "#### DO...\n",
    "1. Write a Python function called `spectral_flux` for the spectral flux density of sunlight. The function should have one input argument (wavelength) and one output argument (spectral flux density). Your function can use `spectral_intensity` inside to avoid repeating any code.\n",
    "2. Plot flux density as a function of wavelength over the range 0.1 to 4 µm. If you have done things correctly, the peak should be about 1750 W $\\rm m^{-2} µm^{-1}$\n",
    "\n",
    "*Turn in*: your Python function and your plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Flux Density\n",
    "\n",
    "The broadband flux density is\n",
    "\\begin{equation}\n",
    "F(\\lambda_1,\\lambda_2) =\\int_{\\lambda_1}^{\\lambda_2}F_\\lambda d\\lambda\n",
    "\\end{equation}\n",
    "\n",
    "#### DO...\n",
    "1. Use numerical integration to compute the broadband flux density over the wavelength range 0.1 µm to 10 µm.\n",
    "2. Compare your value to the measured solar constant (1360 W m<sup>-2</sup>). \n",
    "3. Repeat the integration using several values for the upper limit (e.g. 0.5, 1, 2, 4, 10, and 20 µm). \n",
    "4. What range of wavelengths contribute at least 90% of the total solar flux? What can you conclude about the importance of infrared wavelengths ($\\lambda$ > 0.76 µm) to the solar flux?\n",
    "\n",
    "*Turn in*: your answers to 1, 2, 4\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Extra: Double Integrals \n",
    "\n",
    "This is section is optional and extra.\n",
    "\n",
    "Double (and triple) integrals can also be evaluated numerically. In this example, you will compute the downward radiation flux density at the surface from the intensity under cloudy, overcast skies. \n",
    "\n",
    "Let $I(\\theta,\\phi)$ be the radiation intensity from direction $(\\theta,\\phi)$, where $\\theta$ is the zenith angle and $\\phi$ is the azimuth angle. The surface flux is related to the intensity by\n",
    "\\begin{equation}\n",
    "F^\\downarrow = \\int_0^{2\\pi} \\int_0^{\\pi/2} I(\\theta,\\phi)\\sin{\\theta}\\cos{\\theta}\\, d\\theta\\, d\\phi \n",
    "\\end{equation}\n",
    "\n",
    "If radiation were completely isotropic in overcast conditions, then the intensity would be a constant $I(\\theta,\\phi) = I_0$ and the flux would be $\\pi I_0$, as we derived in class. Measurements under overcast skies, however, show that the intensity is not truly isotropic. Rather, the intensity decreases modestly from the zenith to the horizon, following\n",
    "\\begin{equation}\n",
    "I(\\theta,\\phi) = I_0 \\left( \\frac{1 + b \\cos{\\theta}}{1+b} \\right)\n",
    "\\end{equation}\n",
    "where $I_0$ is the zenith intensity and $b \\approx 1.25$ is an empirical parameter that depends on surface albedo. The figures below show the overcast intensity in two equivalent graphs. The left panel shows the dependence on zenith angle $\\theta$ while the right panel shows the intensity across the entire hemisphere of the sky.\n",
    "![Alt text](img/OvercastIntensity.png)\n",
    "\n",
    "#### DO...\n",
    "0. Import `dblquad` from `scipy.integrate`.\n",
    "1. Define a function `I(theta,phi)` for the intensity. Set $I_0=1$ W/m2/sr for the zenith intensity.\n",
    "2. Define a second function `f(theta,phi)` for the integrand $I(\\theta,\\phi)\\sin{\\theta}\\cos{\\theta}$.\n",
    "3. Use the function `dblquad` to evaluate the integral and find the surface flux density. The command `dblquad(f,phi_0,phi_1,theta_0,theta_1)` will evaluate the inner integral of $\\theta$ over the interval `[theta_0,theta_1]` followed by the outer integral of $\\phi$ over the interval `[phi_0,phi_1]`. \n",
    "4. Compare your value from 3 to what you would get in isotropic conditions. In the isotropic case with $I_0=1$ W/m2/sr, the flux would be $F^\\downarrow=\\pi$ W/m2. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, the double integral could be simplified to a single integral by recognizing that the $\\phi$ integral is separable and $\\int_0^{2\\pi}d\\phi=2\\pi$, so \n",
    "\\begin{equation}\n",
    "F^\\downarrow = 2\\pi \\int_0^{\\pi/2} I(\\theta,\\phi)\\sin{\\theta}\\cos{\\theta}\\, d\\theta \n",
    "\\end{equation}\n",
    "This simplified integral could be evaluated with `quad` instead of `dblquad`, yielding the same value.\n",
    "\n",
    "`dblquad` can also evaluate more complex integrals, such as the following. Suppose we want to compute the downward flux from the large circular patch shown below. \n",
    "<center><img src=\"img/image.png\" width=\"200\"></center>\n",
    "\n",
    "The perimeter of the circle is $\\theta(\\phi)=\\sin(\\phi)$ for $\\phi\\in[0,\\pi]$, so the limits of integration change to\n",
    "\\begin{equation}\n",
    "F_{\\circ}^\\downarrow = \\int_0^{\\pi} \\int_0^{\\sin{\\phi}} I(\\theta,\\phi)\\sin{\\theta}\\cos{\\theta}\\, d\\theta\\, d\\phi \n",
    "\\end{equation}\n",
    "\n",
    "To accomplish this with `dblquad`, the upper and lower limits of the inner integral can be functions. In this case the function `np.sin` can be used in place of `theta_1`. \n",
    "\n",
    "#### DO...\n",
    "1. Create a new `dblquad` and modify the integration limits to integrate over the circle shown. Hint: You will need to change the upper limits `phi_1` and `theta_1`, but not the lower limits. Your answer should be around 0.56 W/m2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code here"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "main",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}