solar_insolation_v1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# GEOL 2001: Solar Radiation Lab\n",
    "\n",
    "#### Written by Abby Eckland, October 2020\n",
    "#### Questions? Contact: abigail.eckland@colorado.edu"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Learning Goals:\n",
    "- Become familiar with terms: zenith angle, solar flux/insolation, radiation, energy, and power.\n",
    "- Compare and contrasts the time of year, time of day, and the latitude for impacts on the amount of solar insolation a location will receive (hint: will need this later for discussing glacier growth and melt). \n",
    "- ***New/Revised:*** Devise Python code to calculate solar flux on variations of Earth’s tilts, latitude and time of day for the following days (a) winter solstice, (b) summer solstice, and (c) Fall/Spring equinox.\n",
    "- ***New/Revised:*** Evaluate and apply lessons learned to calculate solar flux using Python Code in Jupyter Notebooks."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Skills:\n",
    "- Dimensional analysis\n",
    "- Working with data\n",
    "- Unit conversions\n",
    "- ***New/Revised:*** Figure making in Python\n",
    "- ***New:*** Building pandas dataframes to analyze data\n",
    "- ***New:*** Visualizing data in Jupyer Notebooks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook, ***you are not expected to write your own code***, but instead, run each line of code line by line, and see what it does. The goal is to increase your comfort level with code and introduce you to some basic coding skills.\n",
    "\n",
    "### <font color = orange >***To run each line of code, line by line, click \"shift-enter\"***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 1: Modeling the Solar Flux by Latitude\n",
    "    \n",
    "We want to model the solar flux at various latitudes on earth. To do this, we will use Python code and the following constants and information.\n",
    "\n",
    "Constants:\n",
    "- Solar constant = 1370 W/m$^2$\n",
    "- Axial Tilt = 23.45${^\\circ}$\n",
    "\n",
    "When the sun's radiation hits a surface at an angle, such as the surface of the earth, the angle between incoming radiation that is perpendicular to the surface is called the **zenith angle, Z**. \n",
    "\n",
    "![Zenith Angle Concept](data/zenith_angle.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can use the following equations to calculate the midday (noon) solar flux at any zenith angle **(Eq. 1)**:\n",
    "\n",
    "$$q(Z) = S*cos(Z)$$ \n",
    "\n",
    "where q is the **insolation**, Z is the **zenith angle**, and S is the **solar constant**. \n",
    "\n",
    "Let's break down this equation. The solar constant is how much incoming radiation is received by earth from the sun, which is always constant. However, the zenith angle controls how much of this incoming radiation is received at different latitudes, which is not always constant. Think of it this way: at higher latitudes (far up north or down south), these latitudes receive less incoming solar radiation and hence, have a smaller insolation. This is why annual temperatures at the poles are so cold!\n",
    "\n",
    "Recall high school trigonometry. For the cosine function, the angle $\\theta_s$ is calculated by dividing the adjacent length by the hypotenuse (remember SOH, CAH, TOA?). The adjacent length is controlled by latitude, and the hypotenuse by the incoming solar flux. This is why we multiply the incoming solar flux by cos(Z)!\n",
    "\n",
    "The angle of the earth's tilt relative to the ecliptic will vary based on time of year. In summer in the northern hemisphere, the earth is tilted toward the sun by 23.45${^\\circ}$, which is why temperatures in the summer are warm. In winter in the northern hemisphere, the earth is tilted away from the sun by 23.45${^\\circ}$, so it is cold. At the equinoxes, the earth is neither tilted toward nor away from the sun, so the tilt does not impact our insolation calculations at these times. See the figures below to visualize axial tilt as it changes in relation to the sun depending on the time of year.\n",
    "\n",
    "![Earth's Annual Rotation](data/earth_annual_rotation.png)\n",
    "![Earth's tilt](data/axial_tilt.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### **Step 1**\n",
    "Write code to model the insolation, or solar flux, at latitudes every 10${^\\circ}$ between -90${^\\circ}$ and 90${^\\circ}$N. To start, we will consider the solar flux at local noon - when the sun is highest in the sky. \n",
    "\n",
    "Create a **Pandas dataframe**, a data structure similar to a table or chart in excel, to make calculations and view the results. Pandas dataframes have columns and rows. We will create a column called \"latitude\" of latitude values between -90${^\\circ}$ and 90${^\\circ}$N to get started."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Blue text following a hashtag like this are comments to help you understand the code, but are not executed as code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import packages to run python code and assign shortcuts\n",
    "import pandas as pd\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create pandas dataframe called Q_summer\n",
    "Q_summer = pd.DataFrame()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Numpy** is a package in Python that is helpful for creating data lists and arrays, among many other capabilities. The function \"arange\" allows us to make a list of values that are equally spaced out by a specified amount. We want to make a column with values from -90 to 90, increasing by 10, so we will use the arange function. See how this is done in the cell below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create column called \"latitude\" and fill with values from -90 to 90 deg N\n",
    "Q_summer['Latitude (N)'] = np.arange(-90, 100, 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Preview the entire dataframe by typing Q_summer\n",
    "Q_summer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### **Step 2**\n",
    "\n",
    "Use Python to calculate insolation at noon during both equinoxes (spring and fall) and both solstices (summer and winter). \n",
    "\n",
    "You will need to break down calculations and equations into smaller steps and go column by column to get the insolation for each of these times of the year at all latitudes. You will need to repeat your calculations for summer solstice, winter solstice, and the equinoxes. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make a new column called Axial Tilt (deg) and assign 23.45 as all values \n",
    "Q_summer['Axial Tilt (deg)'] = 23.45"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Preview the first 5 lines of the dataframe by using the head() function\n",
    "Q_summer.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the zenith angle by subtracting the Axial Tilt from Latitude, both of which are in degrees\n",
    "Q_summer['Zenith Angle (deg)'] = Q_summer['Latitude (N)'] - Q_summer['Axial Tilt (deg)']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the zenith angle in radians. We will use a numpy function to help us with this\n",
    "Q_summer['Zenith Angle (rad)'] = np.radians(Q_summer['Zenith Angle (deg)'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Assign the solar constant as a new column\n",
    "Q_summer['Solar Constant (W/m^2)'] = 1370"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Finally, calculate the insolation during the summer at all latitudes using equation 1\n",
    "Q_summer['q Summer (W/m^2)'] = np.cos(Q_summer['Zenith Angle (rad)']) * Q_summer['Solar Constant (W/m^2)']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This code will get rid of all negative values, since negative insolation is not possible\n",
    "Q_summer['q Summer (W/m^2)'] = [0 if i < 0 else i for i in Q_summer['q Summer (W/m^2)']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# View the results in our pandas dataframe!\n",
    "Q_summer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's consolidate the above lines of code into one line to calculate the insolation during the winter solstice and the equinoxes. ***Remember that you will need to account for how the zenith angle will change in the different seasons!*** For example, in the winter, the zenith angle will be the latitude + axial tilt (in the summer, this was latitude - axial tilt). This is because in the winter, the earth is tilted further away from the sun in the northern hemisphere than in the summer. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Q_winter = pd.DataFrame()\n",
    "Q_winter['Latitude (N)'] = np.arange(-90, 100, 10)\n",
    "Q_winter['Axial Tilt (deg)'] = 23.45\n",
    "Q_winter['Zenith Angle (deg)'] = Q_winter['Latitude (N)'] + Q_winter['Axial Tilt (deg)']\n",
    "Q_winter['Zenith Angle (rad)'] = np.radians(Q_winter['Zenith Angle (deg)'])\n",
    "Q_winter['Solar Constant (W/m^2)'] = 1370\n",
    "Q_winter['q Winter (W/m^2)'] = np.cos(Q_winter['Zenith Angle (rad)']) * Q_winter['Solar Constant (W/m^2)']\n",
    "Q_winter['q Winter (W/m^2)'] = [0 if i < 0 else i for i in Q_winter['q Winter (W/m^2)']]\n",
    "Q_winter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Q_fall = pd.DataFrame()\n",
    "Q_fall['Latitude (N)'] = np.arange(-90, 100, 10)\n",
    "Q_fall['Axial Tilt (deg)'] = 23.45\n",
    "Q_fall['Zenith Angle (deg)'] = Q_fall['Latitude (N)']\n",
    "Q_fall['Zenith Angle (rad)'] = np.radians(Q_fall['Zenith Angle (deg)'])\n",
    "Q_fall['Solar Constant (W/m^2)'] = 1370\n",
    "Q_fall['q Fall (W/m^2)'] = np.cos(Q_fall['Zenith Angle (rad)']) * Q_fall['Solar Constant (W/m^2)']\n",
    "Q_fall['q Fall (W/m^2)'] = [0 if i < 0 else i for i in Q_fall['q Fall (W/m^2)']]\n",
    "Q_fall"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Q_spring = pd.DataFrame()\n",
    "Q_spring['Latitude (N)'] = np.arange(-90, 100, 10)\n",
    "Q_spring['Axial Tilt (deg)'] = 23.45\n",
    "Q_spring['Zenith Angle (deg)'] = Q_spring['Latitude (N)']\n",
    "Q_spring['Zenith Angle (rad)'] = np.radians(Q_spring['Zenith Angle (deg)'])\n",
    "Q_spring['Solar Constant (W/m^2)'] = 1370\n",
    "Q_spring['q Spring (W/m^2)'] = np.cos(Q_spring['Zenith Angle (rad)']) * Q_spring['Solar Constant (W/m^2)']\n",
    "Q_spring['q Spring (W/m^2)'] = [0 if i < 0 else i for i in Q_spring['q Spring (W/m^2)']]\n",
    "Q_spring"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### STOP:\n",
    "\n",
    "**Q1: What do you notice about the insolation for the fall and spring equinoxes? Write your observation in one sentence in the cell below.**\n",
    "\n",
    "#### Your answer goes here:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step 3\n",
    "\n",
    "Plot all your solar insolation data for the solstices and equinoxes as a function of latitute. All four should be on ***one*** plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "fig, ax = plt.subplots(figsize = (10,8))\n",
    "plt.plot(Q_summer['Latitude (N)'], Q_summer['q Summer (W/m^2)'], label = 'Summer Solstice')\n",
    "plt.plot(Q_winter['Latitude (N)'], Q_winter['q Winter (W/m^2)'], label = 'Winter Solstice')\n",
    "plt.plot(Q_fall['Latitude (N)'], Q_fall['q Fall (W/m^2)'], label = 'Fall Equinox')\n",
    "plt.plot(Q_spring['Latitude (N)'], Q_spring['q Spring (W/m^2)'], label = 'Spring Equinox')\n",
    "plt.legend()\n",
    "plt.title('Solar Insolation at Solstices and Equinoxes by Latitude', size = 18)\n",
    "plt.xlabel('Latitude ($^\\circ N$)', size = 14)\n",
    "plt.ylabel('Insolation (W/m$^2$)', size = 14)\n",
    "plt.xlim(-90, 90)\n",
    "plt.ylim(0, 1400)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### STOP: \n",
    "\n",
    "**Q2: Was the observation that you made above in Q1 confirmed by the plot? Why or why not?**\n",
    "\n",
    "#### Your answer goes here:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 2: Modeling the Solar Flux Throughout the Day"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recall **Eq. 1** for calculating insolation:\n",
    "\n",
    "$$q(Z) = S * cos(Z)$$\n",
    "\n",
    "At **midday**, you can simply use the latitude and the axial tilt to calculate the zenith angle. However, to get insolation at **any hour of the day**, you must use the following equation **(Eq. 2)**: \n",
    "\n",
    "$$cos(Z) = sin(\\phi)sin(\\delta) + cos(\\phi)cos(\\delta)cos(h)$$ \n",
    "\n",
    "where $\\phi$ is **axial tilt**, $\\delta$ is the **latitude in degrees**, and h is the **hour angle**.\n",
    "\n",
    "The **hour angle** is a measure of the angle the sun is hitting the surface of the earth relative to that at noon. At night, the insolation is zero because the sun is below the horizon. At noon, the hour angle is zero since the sun is directly above in the sky. At any time other than noon, the hour angle is measured as +/- degrees from noon. There are 360$^\\circ$ in a circle and 24 hours in a day, so to calculate the hour angle, you can use the following formula **(Eq. 3)**:\n",
    "\n",
    "$$h = 15^\\circ / \\mbox{ hour} * (\\mbox{hours before or after noon})$$\n",
    "\n",
    "If before noon, the hour angle should be negative. If after noon, the hour angle should be positive. \n",
    "\n",
    "Let's see an example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 10:30 am local time in hour angle is -22.5 degrees\n",
    "\n",
    "hours = 12-10.5 \n",
    "print(hours, '= hours between noon and 10:30 am')\n",
    "h = 15*(-hours) # negative hours since 10:30 is in the AM\n",
    "print(h, '= hour angle in degrees')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### **Step 1**\n",
    "\n",
    "Using the equations, skills, and knowledge gained from Part 1, in addition to the new equations (Eqs. 2 and 3), model the solar insolation at one latitude over a 24-hour period on each solstice and equinox. For latitude, please use the latitude of Boulder (40$^\\circ$N).\n",
    "\n",
    "It will be helpful to follow the steps and set-up for Part 1 again, only this time your x variable is not latitude (think about what it would be instead). \n",
    "\n",
    "Create a new Pandas dataframe with columns 'Hour' and 'Hour Angle'. Use the numpy \"arange\" function to fill the columns with values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Q_summer_40 = pd.DataFrame()\n",
    "Q_summer_40['Hour'] = np.arange(0, 25, 1)\n",
    "Q_summer_40['Hour Angle'] = np.arange(-180, 195, 15)\n",
    "Q_summer_40"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### **Step 2**\n",
    "\n",
    "Let's calculate cos(Z) using Eq. 2, creating new columns along the way. This time, instead of making new columns with values of latitude, tilt, etc., we will assign a new variable with each of these constants. This will help to declutter our new dataframe and avoid making columns we don't need (that just contain constants)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "latitude = 40\n",
    "tilt_rads = np.radians(23.45)\n",
    "latitude_rads = np.radians(latitude)\n",
    "Q_summer_40['Hour Angle (rad)'] = np.radians(Q_summer_40['Hour Angle'])\n",
    "Q_summer_40['Cos(Z)'] = (np.sin(tilt_rads) * np.sin(latitude_rads)) + (np.cos(tilt_rads) * np.cos(latitude_rads) * np.cos(Q_summer_40['Hour Angle (rad)']))\n",
    "Q_summer_40.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### **Step 3**\n",
    "\n",
    "Now we can calculate solar insolation using Eq. 1, the same method as in Part 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "S = 1370\n",
    "Q_summer_40['q Summer (W/m^2)'] = S * Q_summer_40['Cos(Z)']\n",
    "Q_summer_40['q Summer (W/m^2)'] = [0 if i < 0 else i for i in Q_summer_40['q Summer (W/m^2)']]\n",
    "Q_summer_40"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### **Step 4**\n",
    "\n",
    "Repeat the steps above to calculate the insolation over a day on the remaining solstice and equinoxes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's calculate the insolation over the day during the winter solstice, then the equinoxes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Q_winter_40 = pd.DataFrame()\n",
    "Q_winter_40['Hour'] = np.arange(0, 25, 1)\n",
    "Q_winter_40['Hour Angle'] = np.arange(-180, 195, 15)\n",
    "Q_winter_40['Hour Angle (rad)'] = np.radians(Q_winter_40['Hour Angle'])\n",
    "Q_winter_40['Cos(Z)'] = (np.sin(-tilt_rads) * np.sin(latitude_rads)) + (np.cos(-tilt_rads) * np.cos(latitude_rads) * np.cos(Q_winter_40['Hour Angle (rad)']))\n",
    "Q_winter_40['q Winter (W/m^2)'] = S * Q_winter_40['Cos(Z)']\n",
    "Q_winter_40['q Winter (W/m^2)'] = [0 if i < 0 else i for i in Q_winter_40['q Winter (W/m^2)']]\n",
    "Q_winter_40"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## STOP:\n",
    "\n",
    "**Q3: What variable needed to change between the summer and winter solstice insolation calculations? In other words, what controls the insolation at a specific latitude, i.e., between summer and winter months? Hint: see line 5 in the code cell above.**\n",
    "\n",
    "#### Your answer goes here:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's calculate the solar insolation over a day during the equinoxes at Boulder, CO."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Q_fall_40 = pd.DataFrame()\n",
    "Q_fall_40['Hour'] = np.arange(0, 25, 1)\n",
    "Q_fall_40['Hour Angle'] = np.arange(-180, 195, 15)\n",
    "Q_fall_40['Hour Angle (rad)'] = np.radians(Q_fall_40['Hour Angle'])\n",
    "Q_fall_40['Cos(Z)'] = (np.sin(0) * np.sin(latitude_rads)) + (np.cos(0) * np.cos(latitude_rads) * np.cos(Q_fall_40['Hour Angle (rad)']))\n",
    "Q_fall_40['q Fall (W/m^2)'] = S * Q_fall_40['Cos(Z)']\n",
    "Q_fall_40['q Fall (W/m^2)'] = [0 if i < 0 else i for i in Q_fall_40['q Fall (W/m^2)']]\n",
    "Q_fall_40"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## STOP: \n",
    "\n",
    "**Q4: Why was the tilt *zero* for the fall equinox insolation calculation for Boulder? Do we need to repeat the steps above for the spring equinox?**\n",
    "\n",
    "#### Your answer goes here:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Turns out, we do not need to complete the same steps as above for the Spring equinox, since solar insolation will again be controlled by an axial tilt of zero (same as for the fall equinox)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### **Step 5**\n",
    "\n",
    "Finally, let's plot the insolation at Boulder, CO (40$^\\circ$ N) during the solstices and equinoxes. Please plot each dataset on ***one*** plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize = (10,8))\n",
    "plt.plot(Q_summer_40['Hour'], Q_summer_40['q Summer (W/m^2)'], label = 'Summer Solstice')\n",
    "plt.plot(Q_winter_40['Hour'], Q_winter_40['q Winter (W/m^2)'], label = 'Winter Solstice')\n",
    "plt.plot(Q_fall_40['Hour'], Q_fall_40['q Fall (W/m^2)'], label = 'Fall & Spring Equinox')\n",
    "plt.legend()\n",
    "plt.title('Solar Insolation at Solstices and Equinoxes in Boulder, CO (40$^\\circ$N)', size = 18)\n",
    "plt.xlabel('Hour of Day', size = 14)\n",
    "plt.ylabel('Insolation (W/m$^2$)', size = 14)\n",
    "plt.xlim(0, 24)\n",
    "plt.ylim(0, 1400)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### **Step 6**\n",
    "\n",
    "One unit used when determining if a building could benefit from solar panels is called **peak sun hours**. \n",
    "\n",
    "One peak sun hour = 1000 W/m$^2$ of sunlight. When calculating the total amount of peak sun hours received at any location, you don’t just consider hours with 1000 W/m$^2$ of **solar radiation**. Instead, you need to sum the **total amount of solar insolation** received by the location, which is latitude dependent. You then express that in terms of the equivalent number of hours with 1000 W/m$^2$.  \n",
    "\n",
    "Calculate the sun peak hour values for Boulder, CO."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Peak_Sun_Hours = pd.DataFrame()\n",
    "Peak_Sun_Hours['Season'] = ['Summer', 'Fall', 'Winter', 'Spring']\n",
    "Peak_Sun_Hours['Cumulative Q (W/m^2)'] = (Q_summer_40['q Summer (W/m^2)'].sum(), Q_fall_40['q Fall (W/m^2)'].sum(), Q_winter_40['q Winter (W/m^2)'].sum(), Q_fall_40['q Fall (W/m^2)'].sum())\n",
    "Peak_Sun_Hours['Cumulative Q (kWh)'] = Peak_Sun_Hours['Cumulative Q (W/m^2)']/1000\n",
    "Peak_Sun_Hours"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The peak sun hours for Boulder are ~12, 8, and 4 kWh for summer, fall/spring, and winter, respectively."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## STOP: \n",
    "\n",
    "**Q5: List and briefly describe three factors that might impact how much actual solar powered electricity someone gets. Think about the incoming radiation, the panels themselves, and any other things that might impact net electricity production. Be creative!**\n",
    "\n",
    "E.g., Snow Cover - in winter these may be even more covered. (Especially for Boulder, CO).\n",
    "\n",
    "#### Your answer goes here:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 3: Modeling the Solar Flux for Different Tilts\n",
    "\n",
    "For this part, you will want to think about how the tilt of the earth impacts the incoming radiation. The axial tilt of the earth varies from 22$^\\circ$ - 24.5$^\\circ$ on a timescale of 41,000 years. This is one component of Milankovitch cycles, which are orbital patterns that influence how much solar flux Earth receives. The other two orbital patterns are variation in the earth orbit's eccentricity and precession of the tilt.\n",
    "\n",
    "We want to explore what the insolation would be like if the earth had a larger tilt (>23.45$^\\circ$) or a smaller tilt (<23.45$^\\circ$). These are not purely hypothetical scenarios! The earth had tilts like this in the past, and it will continue to change in the future.\n",
    "\n",
    "#### **Step 1**\n",
    "\n",
    "Write functions that allow you to change the tilt, but keep other variables the same (i.e., solar constant, latitude, etc.). \n",
    "\n",
    "Note that for winter in the northern hemisphere, tilt will be big(ger). In summer, tilt will be small(er). During the equinoxes, tilt will be 0. This is why we will write three different functions for each time of year."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def varying_tilt_winter(tilt):\n",
    "    latitude = np.arange(-90, 100, 10)\n",
    "    zenith = latitude + tilt\n",
    "    zenith_rad = np.radians(zenith)\n",
    "    S = 1370\n",
    "    return S * np.cos(zenith_rad)\n",
    "\n",
    "def varying_tilt_summer(tilt):\n",
    "    latitude = np.arange(-90, 100, 10)\n",
    "    zenith = latitude - tilt\n",
    "    zenith_rad = np.radians(zenith)\n",
    "    S = 1370\n",
    "    return S * np.cos(zenith_rad)\n",
    "\n",
    "def varying_tilt_fall_spring(tilt):\n",
    "    latitude = np.arange(-90, 100, 10)\n",
    "    zenith = latitude\n",
    "    zenith_rad = np.radians(zenith)\n",
    "    S = 1370\n",
    "    return S * np.cos(zenith_rad)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### **Step 2** \n",
    "\n",
    "Now that the functions are written, we can use them to calculate insolation with different tilts for each time of year! Calculate the solar insolation for each solstice and equinox with tilts of 22 (small), 23.45 (now), and 24.5 (big). Creating a variable for each calculation will allow us to easily plot the results later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# assign a variable, tilt_now_winter, equal to the function, varying_tilt_winter, and assign a tilt, 23.45 (now)\n",
    "tilt_now_winter = varying_tilt_winter(tilt = 23.45)\n",
    "# let's check to see that the function worked\n",
    "tilt_now_winter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a new variable for small tilt during the winter\n",
    "tilt_small_winter = varying_tilt_winter(tilt = 22)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a new variable for big tilt during the winter\n",
    "tilt_big_winter = varying_tilt_winter(tilt = 24.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### **Step 3**\n",
    "\n",
    "Let's plot varying tilt during the winter at all latitudes!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "latitude = np.arange(-90, 100, 10)\n",
    "\n",
    "fig, ax = plt.subplots(figsize = (10,8))\n",
    "plt.plot(latitude, tilt_now_winter, label = 'Tilt = 23.45', color = 'green')\n",
    "plt.plot(latitude, tilt_small_winter, label = 'Tilt = 22', color = 'red')\n",
    "plt.plot(latitude, tilt_big_winter, label = 'Tilt = 24.5', color = 'purple')\n",
    "plt.legend()\n",
    "plt.title('Winter Solstice Solar Insolation with Varying Axial Tilt', size = 18)\n",
    "plt.xlabel('Latitude ($^\\circ$N)')\n",
    "plt.ylabel('Insolation (W/m$^2$)')\n",
    "plt.xlim(-90, 90)\n",
    "plt.ylim(0, 1400)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## STOP:\n",
    "\n",
    "**Q6: Calculate the solar insolation during the summer solstice and equinoxes using the functions defined above.**\n",
    "\n",
    "Hint: you have to create new variables, e.g., \"tilt_now_summer\" and set this equal to the name of the function, e.g., \"varying_tilt_summer(tilt = 23.45)\". Remember you have to create a new variable for each different tilt. Make sure to enter the tilt value as an argument in the function (arguments are specified in the parenthesis, and tilt is the current argument we care about).\n",
    "\n",
    "To create the plots, copy and paste the plotting code from the above cell. All you have to change is the y variable. For example, \"plt.plot(latitude, **tilt_now_summer**, label = 'Tilt = 23.45', color = 'green')\". Also, don't forget to change the title of your figure to specify the time of year (i.e., \"**Summer Solstice** Solar Insolation with Varying Axial Tilt\"). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# assign a variable, tilt_now_summer, equal to the function, varying_tilt_summer, and assign a tilt, 23.45 (now)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a new variable for small tilt during the summer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a new variable for big tilt during the summer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot varying summer tilt here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# assign a variable, tilt_now_fall_spring, equal to the function, varying_tilt_fall_spring, and assign a tilt, 23.45 (now)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a new variable for small tilt during the fall/spring"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a new variable for big tilt during the fall/spring"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot varying fall/spring tilt here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## STOP:\n",
    "\n",
    "**Q7: In one to two sentences, describe the changes in solar flux between the latitude range -30$^\\circ$ to 30$^\\circ$N and the range 60$^\\circ$ to 90$^\\circ$(N or S) latitude among the three tilts.**\n",
    "\n",
    "**Your answer goes here:**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q8: What part of the earth system might these affect most strongly and at what latitude?**\n",
    "\n",
    "**Your answer goes here:**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q9: Can you think of any global climate events that the earth's tilt might be related to? Explain your answers.**\n",
    "\n",
    "**Your answer goes here:**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Congratulations! You have finished the lab."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}