Notebook

Introduction to Programming in Python: Functions¶

This practical is the final one in our exploration of basic concepts in programming. It builds on the previous 3. It is based on an almagamation of course and examples, including the following (which you are welcome to explore on your own!):

Introduction to programming for Geoscientists (with Python) by Gerard Gorman and Christian Jacobs: http://ggorman.github.io/Introduction-to-programming-for-geoscientists/lecture_series/
Introduction to scientific programming in Python by the UCL graduate school: http://www.cs.ucl.ac.uk/scipython/index.html
Programming with Python by Software Carpentry: http://swcarpentry.github.io/python-novice-inflammation/
CS For All: Introduction to Computer Science and Python Programming by HarveyMuddX at edX: https://www.edx.org/course/cs-all-introduction-computer-science-harveymuddx-cs005x-0

Recommended Reading: Think Python, Sections 3.1, 3.5-3.9, 6.1 and SciPy Lecture Notes, Sections 1.2.4.1-1.2.4.3

MAKE SURE TO EXECUTE THIS CELL BEFORE YOU START YOUR WORK!

In [ ]:

import numpy
import matplotlib.pyplot as pyplot
%matplotlib inline

EXERCISE 1a

The function we defined in the main notebook (F_to_C) is copied below. First, execute the function cell so that the function has been defined in this notebook.

In [ ]:

def F_to_C(Temp_in_F):
    Temp_in_C = (5.0/9.0)*(Temp_in_F-32)
    return Temp_in_C

Now, convince yourself that you understand how to call a function by using the function we have already defined (F_to_C) to calculate the temperature in degrees Celsius when the temperature in Fahrenheit is:

Try using different variable names each time.

*You do not need to redefine the function. Just use it with different input values for F.*

Check your answer using Google.

In [ ]:

EXERCISE 1b

Let's take a closer look at the syntax of function definitions by - as usual - finding out what kind of errors we get when we make a mistake.

Try out each of the examples below without changing anything so you see the error message, then fix the error. Also add a comment explaining what the error was.

i. Original:

In [ ]:

# Define the function
def F_to_C(Temp_in_F)
    Temp_in_C = (5.0/9.0)*(Temp_in_F-32)
    return Temp_in_C

# Test the function
print( F_to_C(100) )

i. Fixed:

In [ ]:

ii. Original:

In [ ]:

# Define the function
def F_to_C(Temp_in_F):
    Temp_in_C = (5.0/9.0)*(F_Temp-32)
    return Temp_in_C

# Test the function
print( F_to_C(100) )

ii. Fixed:

In [ ]:

iii. Original:

In [ ]:

# Define the function
def F_to_C(Temp_in_F):
    Temp_in_C = (5.0/9.0)*(Temp_in_F-32)
    return Temp_in_C

# Test the function
print( F_to_C_conversion(100) )

iii. Fixed:

In [ ]:

iv. Original:

In [ ]:

# Define the function
def F_to_C(Temp_in_F):
    Temp_in_C = (5.0/9.0)*(Temp_in_F-32)
    return

# Test the function
print( F_to_C(100) )

iv. Fixed:

In [ ]:

v. Original:

In [ ]:

# Define the function
def F_to_C(Temp_in_F):
    Temp_in_C = (5.0/9.0)*(Temp_in_F-32)
return Temp_in_C

# Test the function
print( F_to_C(100) )

v. Fixed:

In [ ]:

EXERCISE 2a

Recall that once we know the solar flux $S$ at a given distance $r$, we can also calculate the effective temperature of a planet $T_e$ in Kelvin (in the absence of an atmosphere): $$T_e(r) = \left (\dfrac{S}{4\times\sigma}\times(1-A) \right )^{0.25}$$

where:

$\sigma$ is the Stefan-Boltzmann constant, equal to $5.67\times 10^{-8}$ W/m$^2$/K$^4$
$A$ is the albedo (the amount of solar energy that is reflected by the planet's surface), and is unitless. For now, you can assume $A$ = 0.3.

Adapt the function solar_flux (copied below), so that *after* it calculates $S$, it then calculates $T_e$. (Hint: leave the existing lines intact but add an extra line after the line that starts S=...)
Give the function a new name (suggested name: effective_temperature), and make sure that it takes $r$ as input and returns $T_e$ as output.
Test your function by calling (i.e. using) it for $r$ values between 1 and 20. Remember, after you *define* the function you need to *use* the function to get output!
Plot the result with $T_e$ on the y-axis and $r$ on the x-axis.

Hint: We programmed this equation in the Week 2 practical. If you want to save yourselves a headache now, copy the function from that practical. (Side note: when you're programming, you should always look for shortcuts like this that reduce the risk of errors.)

Extra: You can actually call one function within another function - so you could leave the function solar_flux intact, and use it inside your new effective_temperature function.

In [ ]:

# Function to compute the solar flux at a given distance
def solar_flux(r):

    # Define constant variables
    r0 = 1.00  # Distance from sun to Earth in AU
    S0 = 1366  # Solar flux at Earth
    
    # Add Albedo and sigma
    A = 0.3           # albedo, unitless
    sigma = 5.67e-8   # Stefan-Boltzmann constant, in W/m2/K4

    # Use the equation to calculate the solar flux
    S = S0*(r0/r)**2
    
    # return solar flux S
    return S

In [ ]:

# Call the function

In [ ]:

# Plot the result

EXERCISE 2b

1. Copy your function from Exercise 2a and adapt it so that it returns both $S$ and $T_e$ as outputs. (Hint: You only need to change ONE line!)

Call (use) your function with $r$ values between 1 and 20, storing the $S$ and $T_e$ output in two new variables .

Plot the output with $S$ on the x-axis and $T_e$ on the y-axis.

In [ ]:

# Copy and adapt the function

In [ ]:

# Call the function

In [ ]:

# Plot the result

EXERCISE 2c

1. Copy your function from Exercise 2b and adapt it so that the albedo (A in the equation) becomes an additional input argument (instead of being set as a constant in the function).

Call (use) your function with $A$ = 0.3 and $r$ values between 1 and 20, storing the $S$ and $T_e$ output and plotting $T_e$ as a function of $r$ (this should look just like the plot in 2a, with $T_e$ on the y-axis and $r$ on the x-axis).

Repeat step 2, but this time with A=0.9 (an icy planet!). Does the temperature look different? (Hint: it should -- check the y-axis values!)

In [ ]:

# Copy and adapt the function

In [ ]:

# Call the function with A = 0.3 and plot the result

In [ ]:

# Call the function with A = 0.9 and plot the result

EXERCISE 2d

1. Copy your function from Exercise 2c and adapt it so that the albedo (`A`) has a default value of 0.3.

Call (use) the function for $r$ values between 1 and 20, and plot temperature ($T_e$) as a function of distance ($r$). Do this three different times:

A. Specify in your call that A = 0.3
B. Do not specify a value for A
C. Specify that A = 0.6

The temperatures that you find for A and B should be the same. The temperature that you find for C should be colder!

In [ ]:

# Copy and adapt the function

In [ ]:

# A. Call the function specifying A = 0.3 and plot the result

In [ ]:

# B. Call the function without specifying A and plot the result

In [ ]:

# C.Call the function specifying A = 0.6 and plot the result

EXERCISE 2e

Copy the code you used to make the final plot from Exercise 2d here and modify it so that it:
- Uses any color **except** red or blue (bonus points for creativity!).
- Shows points and the line through them
- Zooms in on the region between 5 and 15 AU.

In [ ]:

EXERCISE 3: BRINGING IT ALL TOGETHER

This final exercise is designed to give you a chance to practice everything you have learned so far today (and in the last few weeks). Make sure to spend time thinking about what your answers mean.

This week, we will continue developing an understanding of Daisyworld (look back to the previous practicals and/or the textbook if you don't remember!). Two weeks ago, we made a plot that showed how the temperature of the Daisyworld planet changed as the number of daisies increased. Last week, we made a plot that showed how daisy growth changes as temperature increases.

This time we will start putting these two ideas together.

EXERCISE 3a

We will start by writing a function that calculates the temperature of the planet when the fraction covered by daisies is provided as an input.

Luckily, we have already written most of the code! Back in Week 3 (Exercise 5), we determined that if we know what fraction of the planet is covered by flowers (frac_flower) we can calculate the temperature (Te):

# Define the albedos of flowers and the albedo of soil (unitless) albedo_flower = 0.75 albedo_soil = 0.4 # The fraction covered by soil is the rest of the planet, or 1.0 - the fraction covered by flowers frac_soil = 1.0 - frac_flower # The albedo is the sum of the two albedos after each has been multiplied by its fractional coverage. albedo = frac_soil * albedo_soil + frac_flower * albedo_flower # Define the solar flux on Daisyworld, in W/m2 S = 3700.0 # Define the Stefan-Boltzmann constant, in W/m2/K4 sigma = 5.67e-8 # Define the temperature as a function of the albedo Te = ((S*(1-albedo))/(4*sigma))**0.25

For this exercise:

Copy the code above and modify it so that it becomes a function that takes the fraction of the planet covered by daisies as an input and returns the temperature as an output. (Hint: what does a function need? Think about function names, indentations, colons, arguments and return values...)
To test that it is working, call (use) the function to find the temperature when:

a. 50% of the planet is covered by daisies (i.e., the fraction covered is 0.5).
b. 80% of the planet is covered by daisies (i.e., the fraction covered is 0.8).

Hint: If in step 2 above a and b give you the same value, or you get more than one temperature value as output, then you will know something has gone wrong.

*Do not delete the %reset in the cell below* (this makes sure your program isn't just "remembering" something you did before).

In [ ]:

%reset -s -f
import numpy
import matplotlib.pyplot as pyplot
%matplotlib inline

# Define your function here

In [ ]:

# Call your function here

EXERCISE 3b

Now we will write a function that calculates the growth rate of daisies when the temperature of the planet is provided as an input.

Recall from Week 4 that the growth rate of the *population* of flowers (in units of fraction of the planet covered by flowers / year) is defined as: $$1-0.005 \times (295.5-T_e)^2$$ and can never be below zero (this would represent daisy death, which we will deal with separately).

Good news! We already wrote most of this code in Week 4 (Exercise 5). If we know the temperature of the planet in Kelvin (Te), we can calculate the growth rate (growth_rate):

# Test condition for temperature and calculate growth rate if Te <= 281: growth_rate=0 elif Te < 310: growth_rate=1-0.005*(295.5-Te)**2 else: growth_rate=0

(Some of you might notice this is simpler than the version we used before - that is because we were working with a whole array of temperature values, but now we will just need to use one at a time.

For this exercise:

Copy the code above and modify it so that it becomes a function that takes the temperature of the planet as an input and returns the growth rate as an output.
To test that it is working, call (use) the function to find the growth rate when:

a. the temperature is 289 K
b. the temperature is 269 K

At which temperature are the daisies growing faster, and why?

Hint: If the two different temperature inputs in 3 give you the same value, or you get more than one value as output, then you will know something has gone wrong.

In [ ]:

# Define your function here

In [ ]:

# Call your function here

In [ ]:

# Answer the question here

EXERCISE 3c

You might remember from last week that flowers die at a constant rate of 0.3 (fraction/year).

Using your two functions, you are now going to experiment with different values for the fraction of the planet covered by daisies. For each experiment you should:

First calculate the temperature of the planet.
Use that temperature to calculate the flower growth rate.

Note that the way we have set up these functions, you can only use ONE input value -- they won't work with arrays. So here you are just going to pick some random values for the fraction of the planet covered by flowers and use the functions you've already written to test them, one at a time.

Here's an example, using what we've already found in 3a and 3b. First, we wanted to understand what happens when we start with 0.5 of the planet covered by flowers. In 3a, we used the function defined in 3a with 0.5 as input to calculate the effective temperature. We found that the effective temperature in that scenario was about 289 K. Then in 3b we used the function defined in 3b with 289 K as input to calculate the growth rate. We found the growth rate was about 0.79. In other words, by using our two functions one after another, we found that with 0.5 fractional coverage, the temperature is 289 K and flowers are growing at a rate of 0.79.

We did the same thing starting with a starting fraction of 0.8. We called the first function with 0.8 as input (3a), found an effective temperature of 269 K, and used that as input to the second function (3b) to find a growth rate of 0. At 0.8 fractional coverage, the temperature is 269 K and the growth rate is 0 (i.e. flowers are not growing).

Your job is to repeat this process choosing different numbers as your starting point. We already know that at 0.5 coverage flowers are growing rapidly, and at 0.8 coverage they are not growing at all. Test out some different numbers (your choice!) to answer the following questions:

For what range of initial daisy coverage values are the flowers increasing in area (i.e. growth > death)? Why?
In that range (when flowers are increasing), is the temperature of the planet likely to get warmer or colder?
For what range of daisy coverage values are the flowers decreasing in area? Why? What does that mean for the temperature?

Use comments to show your thinking as you experiment. You can more cells if you like.

Hint: DO NOT re-define your functions. All you need to do now is call them with different input values.

In [ ]:

# Call your functions here.

In [ ]:

# Extra cell to run more experiments if desired

In [ ]:

# Extra cell to run more experiments if desired

In [ ]:

# Extra cell to run more experiments if desired

In [ ]:

# Extra cell to run more experiments if desired

# To add even more, use the + symbol above (underneath "File")

In [ ]:

# Answer the questions here

# Q1
#---
# 
#
#
#

# Q2
#---
# 
#
#
#

# Q3
#---
# 
#
#
#