#!/usr/bin/env python
# coding: utf-8

# # Introduction to programming for Geoscientists (through Python)

# # Lecture 3 solutions
# ## Gerard J. Gorman (g.gorman@imperial.ac.uk) http://www.imperial.ac.uk/people/g.gorman

# ## <span style="color:blue">Exercise 3.1: Make a table of function values</span>
# 
# Write a program that prints a table with $t$ values in the first column and the corresponding $y(t) = v_0 t − 0.5gt^2$ values in the second column.
# Use $n$ uniformly spaced $t$ values throughout the interval [0, $2v_0/g$]. Set $v0 = 1$, $g = 9.81$, and $n = 11$.

# In[2]:


# Make a table of function values

v0 = 1.0
g = 9.81
n = 11

# The step size between each uniformly spaced value. 
# The end-points of the interval are included.
step = (2*v0/g)/(n-1) 

for i in range(0, n):
    t = step*i
    y = v0*t - 0.5*g*t**2
    print("%f, %f" % (t, y))


# Modify the program so that the $t$ and $y$ values are stored in two lists *t* and *y*. Thereafter, transverse the lists with a *for* loop and write out a nicely formatted table of *t* and *y* values using a *zip* construction.

# In[7]:


# Use list comprehension to build the lists
t = [step*i for i in range(n)]
y = [v0*ti - 0.5*g*ti**2 for ti in t]

for (T,Y) in zip(t,y):
    print("%f, %f" % (T, Y))


# ## <span style="color:blue">Exercise 3.2: Implement a Gaussian function</span>
# 
# Make a Python function *gauss*( *x*, *m*=0, *s*=1) for computing the Gaussian function 
# $$f(x)=\frac{1}{\sqrt{2\pi}s}\exp\left(-\frac{1}{2} \left(\frac{x-m}{s}\right)^2\right)$$
# Call the function and print out the result for x equal to −5, −4.9, −4.8, ..., 4.8, 4.9, 5, using default values for *m* and *s*.

# In[7]:


# Implement a Gaussian function

from math import exp, pi, sqrt

def gaussian(x, m = 0, s = 1):
    # Note: 'm' is the mean, and 's' is the standard deviation.
    coefficient = 1/(sqrt(2*pi)*s)
    result = coefficient*exp(-0.5*(float(x - m)/s)**2) # Be careful here. x, m, and s are supplied as integers - let's cast the numerator to a float to ensure we don't encounter integer division problems.
    return result

x = -5.0
while x <= 5.0:
    print("For x = %f, gaussian(x) = %f" % (x, gaussian(x)))
    x = x + 0.1 # Increment x by a step size of 0.1.


# ## <span style="color:blue">Exercise 3.3: Express a step function as a Python function</span>
# 
# The following "step" function is known as the Heaviside function and
# is widely used in mathematics:
# $$H(x)=\begin{cases}0, & \text{if $x<0$}.\\\\
# 1, & \text{if $x\ge 0$}.\end{cases}$$
# Write a Python function H(x) that computes H(x).

# In[8]:


# Express a step function as a Python function

def step(x):
    if(x < 0):
        H = 0
    else: # Otherwise x must be greater than or equal to zero
        H = 1
    return H # Return the result

x = 0.5
print(step(x)) # The value that is returned by the function will be printed out here. For the case of x = 0.5, this should print '1'.