Code

Here you can find all of the code samples from A Student’s Guide to Python for Physical Modeling. The samples are available in three formats:

Some of the samples are drawn directly from the text. Others have been created to illustrate a point or technique. README.txt contains a description of all files in the folder.

These code samples are distributed under the BSD license. Be sure to read LICENSE.txt before using any of the code samples.

This Web page contains the HTML file. If you already know how to use the files, you can Jump to The Index.

Using the files

Please read LICENSE.txt before using any of the code samples. README.txt contains a brief description of all the available files.

All of the .py files and the plain text code has been tested with Python 3.4. Most programs also run correctly with Python 2.7. However, if you are using Python 2, you should execute the following two commands before running the programs to ensure they perform as intended:

from __future__ import division
input = raw_input

Two of the scripts, perrin.py and import_text.py, load data files. These data files are included in code_samples.zip and data_sets.zip. They should be copied into the current working directory before executing the scripts.

Code Archive

You can download a zipped folder full of code samples by following this link:

code_samples.zip

The link above will take you to a Google drive Web page where you can download a file called code_samples.zip. To save the file to your computer, click on the Download icon at the top of the page. This is a horizontal bar with a downward-pointing arrow above it.

If you click on code_samples.zip you will be able to see the folders it contains, and you can look inside these. However, you will not be able to download any files this way. You must download the entire archive by clicking on the Download icon.

Once you unzip the downloaded file, you should find a folder called code/. This folder contains 45 files: about 40 Python programs, a couple data files used by the Python files, and two files called LICENSE.txt and README.txt.

Text File

You can download a single text file that contains all of the code samples by following this link:

code_samples.txt

The link above will take you to a Google drive Web page where you can download a file called code_samples.txt. To save the file to your computer, click on the Download icon at the top of the page. This is a horizontal bar with a downward-pointing arrow above it.

code_samples.txt contains the text of about 40 Python programs, which can be copied and pasted into the text editor of your choice. (This is not a Python script or module. Trying to execute code_samples.txt with Python will result in an error.)

To use this file, first locate the file name you are looking for. Then, use your mouse to select all of the text between lines of asterisks. For example,

***************************************************************************
string_format.py
***************************************************************************

COPY ALL TEXT IN THIS REGION

***************************************************************************
string_percent.py
*************************************************************************** 
...

After you copy and paste this text into your editor, you should be able to run the program with a Python interpreter.

HTML File

On this page, you will find the text of about 40 Python programs, which can be copied and pasted into the text editor of your choice.

If you prefer to work offline, you can save the HTML file that contains the index and code samples to your own computer:

code_samples.html

Just click on the link of the file name you are looking for. Then, use your mouse to select all of the text between lines of asterisks. For example,

***************************************************************************
string_format.py
***************************************************************************

COPY ALL TEXT IN THIS REGION

***************************************************************************
string_percent.py
*************************************************************************** 
...

After you copy and paste this text into your editor, you should be able to run the program in a Python interpreter. Your browser’s “Back” button will return you to the top of this page after you follow a link.

Long lines of code may be wrapped in your browser, but they should paste correctly into your text editor.






Index

LICENSE.txt
README.txt
string_format.py
string_percent.py
for_loop.py
while_loop.py
vectorize.py
projectile.py
branching.py
nesting.py
import_text.py
save_load.py
perrin.py
print_write.py
simple_plot.py
graph_modifications.py
line3d.py
subplots.py
rotate.py
average.py
histogram.py
contour.py
matrix_inversion.py
quadrature.py
simple_oscillator.py
solve_ode.py
parametric_oscillator.py
quiver.py
gradient.py
streamlines.py
walker.py
html_movie.py
waves.py
sympy_examples.py
convolution.py
scope.py
name_collision.py
fancy_plot.py
legend.py
measurement.py
random_walk.py
surface.py
surprise.py









































Source

***************************************************************************
LICENSE.txt
***************************************************************************

Copyright (c) 2014, Jesse M. Kinder and Philip C. Nelson
All rights reserved.

Redistribution and use in source and binary forms, with or without modification,
are permitted provided that the following conditions are met:
    * Redistributions of source code must retain the above copyright
      notice, this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of the authors nor the names of other contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL JESSE M. KINDER OR PHILIP C. NELSON BE LIABLE FOR
ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

================================================================================
Additional license information for html_movie.py:
    html_movie.py was adapted from the scitools module created by Hans Petter
    Langtangen.  The scitools license is reproduced below.
================================================================================
Copyright (c) 2007-2009, Hans Petter Langtangen <hpl@simula.no> and
Simula Resarch Laboratory.

All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are
met: 

    * Redistributions of source code must retain the above copyright
      notice, this list of conditions and the following disclaimer.

    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in
      the documentation and/or other materials provided with the
      distribution. 

    * Neither the name of Simula Research Laboratory nor the names of
      its contributors may be used to endorse or promote products
      derived from this software without specific prior written 
      permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.






***************************************************************************
README.txt
***************************************************************************

Use of these code fragments is subject to the terms of LICENSE.txt.

CONTENTS:
------------------------------------------------------------------------- 
string_format.py
------------------------------------------------------------------------- 
Illustrate string formatting using the .format() method.

------------------------------------------------------------------------- 
string_percent.py
------------------------------------------------------------------------- 
Illustrate string formatting using the % method.

------------------------------------------------------------------------- 
for_loop.py
------------------------------------------------------------------------- 
Use a for loop to generate multiple solutions to the quadratic equation.

This script illustrates the funamental form of a for loop.  For alternate
solutions to the same problem, see while_loop.py and vectorize.py.

------------------------------------------------------------------------- 
while_loop.py
------------------------------------------------------------------------- 
Use a while loop to generate multiple solutions to the quadratic equation.

For alternate solutions to the same problem, see for_loop.py and
vectorize.py.

------------------------------------------------------------------------- 
vectorize.py
------------------------------------------------------------------------- 
Use vectorized operations to generate multiple solutions to the quadratic
equation.

For alternate solutions to the same problem, see for_loop.py and
while_loop.py.

------------------------------------------------------------------------- 
projectile.py
------------------------------------------------------------------------- 
Calculate how long an object is in the air when thrown from a spcified
height with a range of initial speeds assuming constant acceleration due
to gravity:
    0.5 * g * t**2 - v0 * t - y0 = 0

This script illustrates good coding practice in the solution of a simple
problem: parameters with descriptive names, comments, whitespace, and
blocking with '#%%' for debugging in Spyder.

-------------------------------------------------------------------------
branching.py
------------------------------------------------------------------------- 
This script illustrates branching with the use of multiple conditional
statements:
    if <condition1>:
        ...
    elif <condition2>:
        ...
    else:
        ...

------------------------------------------------------------------------- 
nesting.py
------------------------------------------------------------------------- 
Use nested for loops to fill a two-dimensional array of values.

This script illustrates "nesting"---one for loop inside of another.

------------------------------------------------------------------------- 
import_text.py
------------------------------------------------------------------------- 
Load data from a text file by reading the file line by line.

This script reads in data from a text file and stores it in a NumPy array.
It can be adapted to load data from files that are difficult or impossible
to load with NumPy's np.loadtxt function.

------------------------------------------------------------------------- 
save_load.py
------------------------------------------------------------------------- 
Save array data using NumPy's available methods, then load saved data.

This script demonstrates the various methods of saving and loading data
using NumPy arrays.

------------------------------------------------------------------------- 
perrin.py
------------------------------------------------------------------------- 
Generate figure displaying Perrin's experimental data on Brownian motion.
This script requires the data set 04brownian/g26perrindata.npy.

The script illustrates loading and plotting a data set and includes
LaTeX formatting of axis labels and grid lines.

You must copy g26perrindata.npy into the current folder to run the script.

------------------------------------------------------------------------- 
print_write.py
------------------------------------------------------------------------- 
Write same data to a file and print to display.

This script illustrates the similarites between writing to a text file and
printing to the screen.

------------------------------------------------------------------------- 
simple_plot.py
------------------------------------------------------------------------- 
Create and display a basic plot.

------------------------------------------------------------------------- 
graph_modifications.py
------------------------------------------------------------------------- 
This script creates a simple plot with two lines, then modifies several
features of the plot, including axis labels, labels and legend, line
style, tick labels, and title.

------------------------------------------------------------------------- 
line3d.py
------------------------------------------------------------------------- 
Create a three-dimensional parametric plot.

This script demonstrates how to create three-dimensional plots using the
Axes3D method from the mpl_toolkits.mplot3d module.  See also surface.py.

------------------------------------------------------------------------- 
subplots.py
------------------------------------------------------------------------- 
Create four plots in the same figure.

This script demonstrates PyPlot's subplot method, which can be used to
display several plots side-by-side in the same figure.

------------------------------------------------------------------------- 
rotate.py
------------------------------------------------------------------------- 
Define function to rotate a vector in two dimensions.

------------------------------------------------------------------------- 
average.py
-------------------------------------------------------------------------
Compute and return the cummulative average of an array.

This script illustrates the principles of functional programming.

------------------------------------------------------------------------- 
histogram.py
-------------------------------------------------------------------------
Create histograms of random numbers.

This script illustrates how to use NumPy and PyPlot to create histograms
and bar plots.

------------------------------------------------------------------------- 
contour.py
-------------------------------------------------------------------------
Create a labeled contour plot.

This script illustrates how to generate a grid of coordinates for contour
and surface plots.  It also demonstrates some options of plt.contour and
shows how to label contour lines.

------------------------------------------------------------------------- 
matrix_inversion.py
------------------------------------------------------------------------- 
Invert a simple matrix to solve a system of linear equations.

This script illustrates the use of a special method from the SciPy linear
algebra library, scipy.linalg.

------------------------------------------------------------------------- 
quadrature.py
------------------------------------------------------------------------- 
Integrate two functions using quad.

This script demonstrates numerical integration using the quad method of
scipy.integrate.  The first function is a built-in NumPy funciton whose
integral can be computed with pencil and paper for comparison.  The second
is a user-defined function.

------------------------------------------------------------------------- 
simple_oscillator.py
------------------------------------------------------------------------- 
Define function to use in solution of differential equation for a simple
harmonic oscillator.

This script illustrates how to write a function that generates the array
required to integrate a second-order ordinary differential equation.  It
is imported and used in solve_ode.py.

------------------------------------------------------------------------- 
solve_ode.py
------------------------------------------------------------------------- 
Solution of ODE for harmonic oscillator.

This script imports the function F(y,t) in simple_oscillator.py then uses
the odeint method of scipy.integrate to solve the ordinary differential
equation defined by F(y,t).

------------------------------------------------------------------------- 
parametric_oscillator.py
------------------------------------------------------------------------- 
Define a parametric function that accepts 4 parameters then integrate it
using odeint.

This script illustrates two methods for using scipy.integrate's odeint
methods to integrate a function that accepts more than two parameters.

------------------------------------------------------------------------- 
quiver.py
------------------------------------------------------------------------- 
Create a quiver plot.

This script illustrates the use of PyPlot's quiver method.

------------------------------------------------------------------------- 
gradient.py
------------------------------------------------------------------------- 
Calculate and display the gradient of a two-dimensional Gaussian.

This script illustrates the use of NumPy's gradient function and
demonstrates how to display a vector field.  It displays the gradient as a
quiver plot superimposed on a filled contour plot of the Gaussian.

------------------------------------------------------------------------- 
streamlines.py
------------------------------------------------------------------------- 
Create streamlines from a vector field.

This script demonstrates the use of PyPlot's streamplot method for
visualizing solutions to a differential equation defined by a vector
field.

------------------------------------------------------------------------- 
walker.py
------------------------------------------------------------------------- 
Make a movie out of the steps of a two-dimensional random walk.

This script demonstrates the use of the FuncAnimation method of
Matplotlib's animation module to create a movie.  If ffmpeg or mencoder is
installed on this computer, the script will save the movie to an mp4 file.

------------------------------------------------------------------------- 
html_movie.py
------------------------------------------------------------------------- 
Module to generate an HTML document from a collection of images.  When
viewed in a Web browser, the document will display a movie whose frames
are the individual images.

This module is adapted from the scitools library developed by Hans Petter
Langtangen.

------------------------------------------------------------------------- 
waves.py
------------------------------------------------------------------------- 
Create an HTML animation of a moving Gaussian waves.

This script illustrates a method for combining a series of plots into an
animation using HTML and Javascript.  It uses the html_movie.py module,
which is  adapted from the scitools library developed by Hans Petter
Langtangen.

------------------------------------------------------------------------- 
sympy_examples.py
------------------------------------------------------------------------- 
Demonstrate some useful methods available in the SymPy module.

Python may complain about undefined variables if you attempt to run this
script.  init_session() defines several variables, but Python may not be
aware of this.  It is better to run the commands one at a time from the
command line.

------------------------------------------------------------------------- 
convolution.py
------------------------------------------------------------------------- 
This script creates an eLoG (elongated Laplacian of Gaussian) filter that
emphasizes long, vertical lines in a figure.  The effect of the filter is
demonstrated on a plus sign.

------------------------------------------------------------------------- 
scope.py
-------------------------------------------------------------------------
Demonstrate Python's rules of scope.

------------------------------------------------------------------------- 
name_collision.py
-------------------------------------------------------------------------
Illustrate how Python's rules of scope prevent name collisions.

------------------------------------------------------------------------- 
fancy_plot.py
-------------------------------------------------------------------------
Add a title and axis labels to a simple plot.

------------------------------------------------------------------------- 
legend.py
-------------------------------------------------------------------------
Create a plot with a legend to distinguish multiple curves.

------------------------------------------------------------------------- 
measurement.py
------------------------------------------------------------------------- 
Functions to calculate distance between points using different metrics.

This script illustrates the fundamental form of user-defined functions as
well as keyword arguments and default values.

------------------------------------------------------------------------- 
random_walk.py
------------------------------------------------------------------------- 
Monte Carlo simulation of a two-dimensional random walk.

This script illustrates the use of a random number generator to create a
time series for a random walk.

------------------------------------------------------------------------- 
surface.py
------------------------------------------------------------------------- 
Create a three-dimensional surface plot.

This script demonstrates how to create three-dimensional plots using the
Axes3D method from the mpl_toolkits.mplot3d module.  See also line3d.py.

------------------------------------------------------------------------- 
surprise.py
------------------------------------------------------------------------- 
This script will create a familar but interesting image.
It may take about a minute to run.






***************************************************************************
string_format.py
***************************************************************************

# string_format.py  (Python 3.4)
# -------------------------------------------------------------------------
# Illustrate formatting strings using the .format() method.
# ------------------------------------------------------------------------- 
import numpy as np

print("The value of pi is approximately " + str(np.pi))
print("The value of {} is approximately {:.5f}".format('pi', np.pi))

s = "{1:d} plus {0:d} is {2:d}"
print(s.format(2,4,2+4))

print("Every {2} has its {3}.".format('dog','day','rose','thorn'))
print("The third element of the list is {0[2]:g}.".format(np.arange(10)))






***************************************************************************
string_percent.py
***************************************************************************

# string_percent.py  (Python 3.4)
# -------------------------------------------------------------------------
# Illustrate string formatting using the % method.
# ------------------------------------------------------------------------- 
import numpy as np

print("The value of pi is approximately " + str(np.pi))
print("The value of %s is approximately %.5f" % ('pi', np.pi))

s = "%d plus %d is %d"
print(s % (2, 4, 2 + 4))






***************************************************************************
for_loop.py
***************************************************************************

# for_loop.py
# -------------------------------------------------------------------------
# Use a for loop to generate multiple solutions to the quadratic equation.
# ------------------------------------------------------------------------- 
import numpy as np

b, c = 2, -1
for a in np.arange(-1, 2, 0.3):
    x = (-b + np.sqrt(b**2 - 4*a*c)) / (2*a)
    print("a= {:.4f}, x= {:.4f}".format(a,x))






***************************************************************************
while_loop.py
***************************************************************************

# while_loop.py
# -------------------------------------------------------------------------
# Use a while loop to generate multiple solutions to the quadratic equation.
# ------------------------------------------------------------------------- 
import numpy as np

a, b, c = 2, 2, -1
while (b**2 - 4*a*c >= 0):
    x = (-b + np.sqrt(b**2 - 4*a*c)) / (2*a)
    print("a = {:.4f}, x = {:.4f}".format(a,x))
    a = a - 0.3
print("done!")






***************************************************************************
vectorize.py
***************************************************************************

# vectorized.py
# -------------------------------------------------------------------------
# Use vectorized operations to generate multiple solutions to the
# quadratic equation.
# ------------------------------------------------------------------------- 
import numpy as np

b, c = 2, -1
a = np.arange(-1, 2, 0.3)
x = (-b + np.sqrt(b**2 - 4*a*c)) / (2*a)






***************************************************************************
projectile.py
***************************************************************************

# projectile.py (Python 3.4)
# Jesse M. Kinder -- 2014 Nov 11
"""
Calculate how long an object is in the air when thrown from a spcified height
with a range of initial speeds assuming constant acceleration due to gravity:
    0.5 * g * t**2 - v0 * t - y0 = 0
"""
import numpy as np

#%% Initialization of variables.
initial_speed = 0.0         # v0 = initial vertical speed of ball in [m/s]
impact_time = 0.0           # t = time of impact in [s] (computed in loop)

#%% Initialization of parameters.
g = 9.8066                  # gravitational acceleration in [m/s^2]
initial_height = 2.0        # y0 = height ball is thrown from in [m]
speed_increment = 5.0       # how much to increase speed in [m/s] for each iteration
cutoff_time = 10.0          # stop computing after impact time exceeds cutoff

#%% Calculate and display impact time.  Increment initial speed each step.
#   Repeat until impact time exceeds cutoff.
while impact_time < cutoff_time:
    # Use quadratic equation to solve kinematic equation for impact time:
    impact_time = (np.sqrt(initial_speed**2 + 2 * g * initial_height) + initial_speed) / g
    print("speed= {} m/s; time= {:.1f} s".format(initial_speed, impact_time))
    initial_speed += speed_increment
print("Calculation complete.")






***************************************************************************
branching.py
***************************************************************************

# branching.py (Python 3.4)
"""
This script illustrates branching.
"""
import numpy as np

for trial in range(5):
    userInput = input('Pick a number: ')
    userNumber = float(userInput)
    if userNumber < 0:
        print('Square root is not real.')
    else:
        print('Square root of {} is {:.4f}'.format(userNumber, np.sqrt(userNumber)))
    userAgain = input('Another [y/n]? ')
    if userAgain != 'y':
        break

if trial == 4:
    print('Sorry, only 5 per customer.')
elif userAgain == 'n':
    print('Bye!')
else:
    print('Sorry, I did not understand that.')






***************************************************************************
nesting.py
***************************************************************************

# nesting.py
# -------------------------------------------------------------------------
# Use nested for loops to fill a two-dimensional array of values.
# ------------------------------------------------------------------------- 
import numpy as np

# Set dimensions of array.
rows = 3
columns = 4

# Create empty array then fill with values.
A = np.zeros((rows, columns))
for i in range(rows):
    for j in range(columns):
        A[i, j] = i**2 + j**3






***************************************************************************
import_text.py
***************************************************************************

# import_text.py
# -------------------------------------------------------------------------
# Load data from a text file by reading the file line by line.
# ------------------------------------------------------------------------- 
import numpy as np

f = open("HIVseries.csv")
temp_data = []

for line in f:
    print(line)
    x,y = line.split(',')
    temp_data += [ (float(x), float(y)) ]

f.close()

data_set = np.array(temp_data)






***************************************************************************
save_load.py
***************************************************************************

# numpy_saves.py
# -------------------------------------------------------------------------
# Save array data using NumPy's available methods, then load saved data.
# ------------------------------------------------------------------------- 
import numpy as np

# Generate array data.
x = np.linspace(0, 1, 1001)
y = 3*np.sin(x)**3 - np.sin(x)

# Save data to text, .npy, and .npz files.
np.save('x_values', x)
np.save('y_values', y)
np.savetxt('x_values.dat', x)
np.savetxt('y_values.dat', y)
np.savez('xy_values', x_vals=x, y_vals=y) 

# Load saved data.
x2 = np.load('x_values.npy')
y2 = np.loadtxt('y_values.dat')
w = np.load('xy_values.npz')

# Check equality of original and saved data.
print((x2 == x).all())
print((y2 == y).all())
print((w['x_vals'] == x2).all())
print((w['y_vals'] == y2).all())






***************************************************************************
perrin.py
***************************************************************************

# perrin.py
# -------------------------------------------------------------------------
# Generate figure displaying Perrin's experimental data on Brownian motion.
# This script requires the data set 04brownian/g26perrindata.npy.
# ------------------------------------------------------------------------- 
import numpy as np
import matplotlib.pyplot as plt

data = np.load('g26perrindata.npy')

plt.figure(figsize=(10,10))
plt.plot(data[:,0], data[:,1], 'b.', ms=16)
plt.axis([-20,20,-20,20])
plt.grid(b=True,which='major',ls='-',lw=1.5,c='g')
plt.xlabel('$\\Delta x$ [$\\mu$m]', fontsize=24)
plt.ylabel('$\\Delta y$ [$\\mu$m]', fontsize=24)
ax=plt.gca()
ax.set_xticklabels(ax.get_xticks(), weight='bold', fontsize=16)
ax.set_yticklabels(ax.get_yticks(), weight='bold', fontsize=16)
plt.show()






***************************************************************************
print_write.py
***************************************************************************

# print_write.py    (Python 3.4)
# -------------------------------------------------------------------------
# Write same data to a file and print to display.
# ------------------------------------------------------------------------- 
f = open('power.txt','w')

print(" N \t\t2**N\t\t3**N")        # print labels for columns
f.write(" N \t\t2**N\t\t3**N\n")    # write labels to file
print("---\t\t----\t\t----")        # print separator
f.write("---\t\t----\t\t----\n")    # write separator to file

#%% loop over integers from 0 to 10 and print/write results
for N in range(11):
    print("{:d}\t\t{:d}\t\t{:d}".format(N, pow(2,N), pow(3,N)))
    f.write("{:d}\t\t{:d}\t\t{:d}\n".format(N, pow(2,N), pow(3,N)))
f.close()






***************************************************************************
simple_plot.py
***************************************************************************

# simple_plot.py
# -------------------------------------------------------------------------
# Create and display a basic plot.
# ------------------------------------------------------------------------- 
import numpy as np
import matplotlib.pyplot as plt

num_points = 26
x_min, x_max = 0, 4

x_values = np.linspace(x_min, x_max, num_points)
y_values = x_values**2

plt.figure()
plt.plot(x_values, y_values)
plt.show()






***************************************************************************
graph_modifications.py
***************************************************************************

# graph_modifications.py
# -------------------------------------------------------------------------
# This script creates a simple plot with two lines, then modifies many
# features of the plot, including axis labels, labels and legend, line
# style, tick labels, and title.
# ------------------------------------------------------------------------- 
import numpy as np
import matplotlib.pyplot as plt

# Generate data for plot.
num_points = 26
x_min, x_max = 0, 4
x_values = np.linspace(x_min, x_max, num_points)
y_values = x_values**2

# Create empty figure.
plt.figure()

# Plot data.
plt.plot(x_values, y_values, label="Population 1")      # label when plot is created
plt.plot(x_values, x_values**3, label="Population 2")   # label when plot is created
plt.legend()

# Gain control of current Axes object.
ax = plt.gca()

# Give plot a title.
ax.set_title("My First Plot", family='monospace', size=24, weight='bold')

# Label the axes.
ax.set_xlabel("Time [days]")
ax.set_ylabel("Population")

# Change tick labels and font
ax.set_xticklabels(ax.get_xticks(), family='monospace', fontsize=10)
ax.set_yticklabels(ax.get_yticks(), family='monospace', fontsize=10)

# Change the legend
lines = ax.get_lines()                      # returns a list of line objects
lines[0].set_label("Infected Population")   # change labels using line objects
lines[1].set_label("Cured Population")      # change labels using line objects
ax.legend()                                 # display legend in plot

# Make the first line a thick, red, dashed line.
plt.setp(lines[0], linestyle='--', linewidth=3, color='r')

# Change the legend again.
ax.legend(("Healthy", "Recovered"))         # change labels using Axes object

plt.show()






***************************************************************************
line3d.py
***************************************************************************

# line3d.py
# -------------------------------------------------------------------------
# Create a three-dimensional parametric plot.
# ------------------------------------------------------------------------- 
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D     # import 3D plotting tool

fig = plt.figure()                  # create a new figure
ax = Axes3D(fig)                    # create 3D plotting object attached to figure
t = np.linspace(0, 5*np.pi, 501)    # define parameter for parametric plot

ax.plot(t * np.cos(t), t * np.sin(t), t)        # generate 3D plot
plt.show()






***************************************************************************
subplots.py
***************************************************************************

# subplots.py
# -------------------------------------------------------------------------
# Create four plots in the same figure.
# ------------------------------------------------------------------------- 
import numpy as np
import matplotlib.pyplot as plt
from numpy.random import random

t = np.linspace(0, 1, 101)

plt.figure()
plt.subplot(2, 2, 1); plt.hist(random(20))
plt.subplot(2, 2, 2); plt.plot(t, t**2, t, t**3 - t)
plt.subplot(2, 2, 3); plt.plot(random(20), random(20), 'r*')
plt.subplot(2, 2, 4); plt.plot(t*np.cos(10*t), t*np.sin(10*t))
plt.show()






***************************************************************************
rotate.py
***************************************************************************

# rotate.py
# -------------------------------------------------------------------------
# Define function to rotate a vector in two dimensions.
# ------------------------------------------------------------------------- 
import numpy as np

def rotate_vector(vector, angle):
    """
    Rotate two-dimensional vector through given angle.
        vector = (x,y)
        angle = rotation angle in radians (counterclockwise)
    Returns the image of vector under rotation as a NumPy array.
    """
    rotation_matrix = np.array([[ np.cos(angle), -np.sin(angle) ],
                                [ np.sin(angle),  np.cos(angle) ]])
    return np.dot(rotation_matrix, vector)






***************************************************************************
average.py
***************************************************************************

# average.py
# -------------------------------------------------------------------------
# Compute and return the cummulative average of an array.
# ------------------------------------------------------------------------- 
import numpy as np

def running_average(x):
    """
    Return cummulative average of an array.
    """
    y = np.zeros(len(x))                    # empty array to store result
    current_sum = 0.0                       # running sum of elements of x
    for i in range(len(x)):
        current_sum += x[i]                 # increment sum
        y[i] = current_sum / (i + 1.0)      # update running average
    return y






***************************************************************************
histogram.py
***************************************************************************

# histogram.py
# -------------------------------------------------------------------------
# Create histograms of random numbers.
# ------------------------------------------------------------------------- 
import numpy as np
import matplotlib.pyplot as plt
from numpy.random import random as rng

#%% Generate data and view PyPlot's default histogram.
data = rng(1000)

plt.figure()
plt.hist(data)

#%% Get binned data from NumPy and make a colorful histogram where the
#   width of each bin is proportional to the number of elements in it.
counts, bin_edges = np.histogram(data)
bin_size = bin_edges[1] - bin_edges[0]
new_widths =  bin_size * counts / counts.max()

plt.figure()
plt.bar(bin_edges[:-1], counts, width=new_widths, color=['r','g','b'])

#%% Provide logarithmically spaced bin edges rather than using defaults.
log2bins = np.logspace(-8, 0, num=9, base=2)
log2bins[0] = 0.0           # set first bin edge to zero instead of 1/256

plt.figure()
plt.hist(data, bins=log2bins)

plt.show()






***************************************************************************
contour.py
***************************************************************************

# contour.py
# -------------------------------------------------------------------------
# Create a labeled contour plot.
# ------------------------------------------------------------------------- 
import numpy as np
import matplotlib.pyplot as plt

#%% Create a grid of x and y coordinates
x_vals = np.linspace(-3, 3, 21)
y_vals = np.linspace(0, 10, 11)
X, Y = np.meshgrid(x_vals, y_vals)

#%% Generate function values. 
Z = np.cos(X) * np.sin(Y)

#%% Plot and label contours.
plt.figure()
cs = plt.contour(X, Y, Z, 10, linewidths=3, colors='k')
plt.clabel(cs,fontsize=10)
plt.show()






***************************************************************************
matrix_inversion.py
***************************************************************************

# matrix_inversion.py
# -------------------------------------------------------------------------
# Invert a simple matrix to solve a system of linear equations.
# ------------------------------------------------------------------------- 
import numpy as np
from scipy.linalg import inv

#%% Set up and solve C.x = a
a = np.array([-1, 5])
C = np.array([[1, 3], [3, 4]])
x = np.dot(inv(C), a)

#%% Check solution.
error = np.dot(C,x) - a






***************************************************************************
quadrature.py
***************************************************************************

# quadrature.py
# -------------------------------------------------------------------------
# Integrate two functions using quad.
# ------------------------------------------------------------------------- 
import numpy as np
from scipy.integrate import quad

#%% Integrate a built-in function.
upper_limit = np.linspace(0,3*np.pi,16)
cos_integral = np.zeros(upper_limit.size)
for i in range(upper_limit.size):
    cos_integral[i], error = quad(np.cos, 0, upper_limit[i])

#%% Now integrate a user-defined function.
def integrand(x): return np.exp(-x**2/2)

upper_limit = np.linspace(0, 5, 51)
gauss_integral = np.zeros(upper_limit.size)
for i in range(upper_limit.size):
    gauss_integral[i], error = quad(integrand, 0, upper_limit[i])






***************************************************************************
simple_oscillator.py
***************************************************************************

# simple_oscillator.py
# -------------------------------------------------------------------------
# Define function to use in solution of differential equation for a simple
# harmonic oscillator.
# ------------------------------------------------------------------------- 
def F(y, t):
    """
    Return derivatives for 2nd order ODE y'' = -y.
    """
    dy = [0, 0]         # preallocate list to store derivatives
    dy[0] = y[1]        # first derivative of y(t)
    dy[1] = -y[0]       # second derivative of y(t)
    return dy






***************************************************************************
solve_ode.py
***************************************************************************

# solveODE.py
"""
Solution of ODE for harmonic oscillator.
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# Import function to integrate:
from simple_oscillator import F

# array of time values to study
t_min = 0; t_max = 10; dt = 0.1
t = np.arange(t_min, t_max+dt, dt)

# two sets of initial conditions
initial_conditions = [ (1.0, 0.0), (0.0, 1.0) ]

plt.figure()    # Create figure; then, add plots.
for y0 in initial_conditions:
    y = odeint(F, y0, t)
    plt.plot(t, y[:, 0], linewidth=2)

skip = 5
t_test = t[::skip]                          # compare at a subset of points
plt.plot(t_test, np.cos(t_test), 'bo')      # exact solution for y0 = (1,0)
plt.plot(t_test, np.sin(t_test), 'go')      # exact solution for y0 = (0,1)

plt.show()






***************************************************************************
parametric_oscillator.py
***************************************************************************

# parametric_oscillator.py
# -------------------------------------------------------------------------
# Define a parametric function that accepts 4 parameters then integrate it
# using odeint.
# ------------------------------------------------------------------------- 
import numpy as np
from scipy.integrate import odeint

def F(y, t, spring_constant=1.0, mass=1.0):
    """
    Return derivatives for harmonic oscillator:
        y'' = -(k/m) * y
    y = displacement in [m]
    k = spring_constant in [N/m]
    m = mass in [kg]
    """
    dy = [0, 0]         # array to store derivatives
    dy[0] = y[1]        
    dy[1] = -(spring_constant/mass) * y[0]  
    return dy

# -------------------------------------------------------------------------
# Integrate parametric function using two different methods.
# ------------------------------------------------------------------------- 
y0 = (1.0, 0.0)                     # initial conditions
t = np.linspace(0, 10, 101)         # times at which y(t) will be evaluated

# Method 1 -- dummy function
def G(y, t): return F(y, t, 2.0, 0.5)
yA = odeint(G, y0, t)

# Method 2 -- keywords
yB = odeint(F, y0, t, args=(2.0, 0.5))






***************************************************************************
quiver.py
***************************************************************************

# quiver.py
# -------------------------------------------------------------------------
# Create a quiver plot.
# ------------------------------------------------------------------------- 
import numpy as np
import matplotlib.pyplot as plt

coords = np.linspace(-1, 1, 11)
X, Y = np.meshgrid(coords, coords)
Vx, Vy = Y, -X

plt.figure()
plt.quiver(X, Y, Vx, Vy)
plt.show()






***************************************************************************
gradient.py
***************************************************************************

# gradient.py
# -------------------------------------------------------------------------
# Calculate and display the gradient of a two-dimensional Gaussian.
# ------------------------------------------------------------------------- 
import numpy as np
import matplotlib.pyplot as plt

coords = np.linspace(-2, 2, 101)
X, Y = np.meshgrid(coords[::5], coords[::5])    # coarse grid for vector field
R = np.sqrt(X**2 + Y**2)
Z = np.exp(-R**2)
x, y = np.meshgrid(coords, coords)              # fine grid for contour plot
r = np.sqrt(x**2 + y**2)
z = np.exp(-r**2)

ds = coords[5] - coords[0]                      # coarse grid spacing
dX, dY = np.gradient(Z, ds)                     # calculate gradient

plt.figure()
plt.contourf(x, y, z, 25)                   
plt.set_cmap('coolwarm')
plt.quiver(X, Y, dX.transpose(), dY.transpose(), scale=25, color='k')
plt.axis('tight')
plt.show()






***************************************************************************
streamlines.py
***************************************************************************

# streamlines.py
# -------------------------------------------------------------------------
# Create streamlines from a vector field.
# ------------------------------------------------------------------------- 
import numpy as np
import matplotlib.pyplot as plt

# generate grid of points
lower = -2; upper = 2; step = 0.2
coords = np.arange(lower, upper+step, step)
X, Y = np.meshgrid(coords, coords)

# define vector field
Vx, Vy = Y, -X

# display streamlines defined by vector field
plt.figure()
plt.streamplot(X, Y, Vx, Vy, linewidth=2)
plt.show()






***************************************************************************
walker.py
***************************************************************************

# walker.py (Python 3.4)
# Jesse M. Kinder -- 2014 Nov 16
""" Make a movie out of the steps of a two-dimensional random walk. """
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation
from numpy.random import random as rand

# Set number of steps for each random walk.
num_steps = 100

# Create an empty figure of the desired size.
bound = 20
fig = plt.figure()      # must have figure object for movie
ax = plt.axes(xlim=(-bound, bound), ylim=(-bound, bound))

# Create a line and a point with no data.  They will be updated during each
# frame of the animation.
(my_line,) = ax.plot([], [], lw=2)              # line to show path
(my_point,) = ax.plot([], [], 'ro', ms=9)       # dot to show current position

# Generate the random walk data.
x_steps = 2*(rand(num_steps) < 0.5) - 1     # generate random steps +/- 1
y_steps = 2*(rand(num_steps) < 0.5) - 1
x_coordinate = x_steps.cumsum()             # sum steps to get position
y_coordinate = y_steps.cumsum()

# This function will generate each frame of the animation.
# It adds all of the data through frame n to a line
# and moves a point to the nth position of the walk.
def get_step(n, x, y, this_line, this_point):
    this_line.set_data(x[:n+1], y[:n+1])
    this_point.set_data(x[n], y[n])

# Call the animator and create the movie.
my_movie = animation.FuncAnimation(fig, get_step, frames=num_steps, \
                    fargs=(x_coordinate,y_coordinate,my_line,my_point) )

# Save the movie in the current directory.
# *** THIS WILL CAUSE AN ERROR UNLESS FFMPEG OR MENCODER IS INSTALLED. ***
# my_movie.save('random_walk.mp4', fps=30)






***************************************************************************
html_movie.py
***************************************************************************

"""
This module provides a single function called movie that will create a Web page
to display a sequence of images as a movie.

The input files can be given in a list or tuple of file names or specified by a
regular expression that will generate the filenames.  I.e.,
input_files="frame_*.jpg" will use all .jpg files in the current directory with
a name of the form "frame_001.jpg", "frame_002.jpg", etc.
"""

# -------------------------------------------------------------------------
# This module is adapted from the scitools library developed by Hans Petter
# Langtangen.  The source from which this module was derived may be found in the
# file <scitools/lib/scitools/easyviz/movie.py> of the current distribution of the
# scitools module (as of May 30, 2015), available at
# 
#   https://code.google.com/p/scitools/
# 
# The MovieEncoder class and the movie function of <movie.py> from the scitools
# library have been merged into a single function called movie.  The html_movie
# function has been renamed get_html.
# 
# The movie function in the current module takes a collection of input files and
# generates an HTML file that displays the images as a movie.  The output is
# nearly identical to that of scitools.std.movie called with the default 'html'
# encoder.  This module is self-contained and does not require any of the other
# functions or libraries of the scitools module.  The print statements have also
# been modified to work with both Python 2 and Python 3.
# 
# scitools:
# ---------
#   Copyright (c) 2007-2009, Hans Petter Langtangen <hpl@simula.no> and
#   Simula Resarch Laboratory.
# 
#   All rights reserved.
# 
# Both the current module and scitools are distributed under the BSD license.
# See LICENSE.txt for details.
# ------------------------------------------------------------------------- 

# Modules that provide tools to search for file names using regular expressions.
import os
import glob
import re

# Specify types of image files that can be used to create a movie.
_legal_file_types = 'png gif jpg jpeg'.split()



# ------------------------------------------------------------------------- 
def movie(input_files, output_file="movie.html", fps=25):
    """
    Take a list or tuple of image file names or a regular expression that will
    generate those file names and then create an HTML file that uses JavaScript
    to display those images as a movie.

    Images must be PNG, JPEG, or GIF files, with extension .png, .jpg, .jpeg, or
    .gif, respectively.

        input_files:    a list or tuple of file names, or regular expression
        output_file:    name of the output HTML movie file
        fps:            frame rate, in frames per second
    
    The sequence of images will be determined by their order in input_files if
    input_files is a list or tuple.  If input_files is a regular expression like
    "frame_*.jpg", the sequence of images will be the alphanumeric order of all
    filenames that match the search pattern.  Because "frame_10.jpg" comes
    before "frame_2.jpg" in alphanumeric order, it is important to include
    leading zeros in file names: "frame_002.jpg" comes before "frame_010.jpg".
    """

    # Print blank lines before function output statements.
    print('\n\n')

    # Determine the file type of the input files.
    if isinstance(input_files, (tuple,list)):
        file_ = input_files[0]
    elif isinstance(input_files, str):
        file_ = input_files
    else:
        raise ValueError("The input files must be given as either a "\
                         "list or tuple of strings or a string, not '%s'" % \
                         type(input_files))

    # Check that the input files do exist.
    if isinstance(input_files, str):
        # Are input_files on ffmpeg/mpeg2enc format or Unix wildcard format?
        ffmpeg_format = r'(.+)%\d+d(\..+)'
        m = re.search(ffmpeg_format, input_files, re.DOTALL)
        if m:
            wildcard_format = m.group(1) + '*' + m.group(2)
        else:
            wildcard_format = input_files
        all_input_files = glob.glob(wildcard_format)
        if not all_input_files:
            print('No files of the form %s exist.' % input_files)
        else:
            print('Found %d files of the format %s.' % \
            (len(all_input_files), input_files))
    else:
        # User provided a list or tuple of specific filenames.
        all_input_files = input_files
    error_encountered = False
    for f in all_input_files:
        if not os.path.isfile(f):
            print('Input file %s does not exist.' % f)
            error_encountered = True
    if error_encountered:
        raise IOError('Some input files were not found.')

    fname, ext = os.path.splitext(file_)
    if not ext:
        raise ValueError("Unable to determine file type from file name.")
    file_type = ext[1:] # remove the . (dot)
    if not file_type in _legal_file_types:
        raise TypeError("File type must be %s, not '%s'" % \
                        (_legal_file_types, file_type))

    # Create an html file that can play image files
    files = input_files
    if isinstance(files, str):
        files = glob.glob(files)
        files.sort()
    print('\nMaking HTML code for displaying ' + ', '.join(files))

    # Turn frame rate into a pause, in milliseconds, between images.
    interval_ms = 1000.0/fps

    outf = output_file
    if outf is None:
        outf = 'movie.html'
    else:
        # Ensure .html extension
        outf = os.path.splitext(outf)[0] + '.html'
    casename = os.path.splitext(outf)[0]

    # Get the necessary HTML text to create the movie.
    header, jscode, form, footer, files = \
        get_html(files, interval_ms, casename=casename)

    # Write the HTML text to the output file.
    f = open(outf, 'w')
    f.write(header + jscode + form + footer)
    f.close()
    print("\n\nmovie in output file " + outf)

    return None

# ------------------------------------------------------------------------- 
def get_html(plotfiles, interval_ms=300, width=800, height=600,
               casename=None):
    """
    Takes a list plotfiles, such as::

        'frame00.png', 'frame01.png', ...

    and creates javascript code for animating the frames as a movie in HTML.

    The `plotfiles` argument can be of three types:

      * A Python list of the names of the image files, sorted in correct
        order. The names can be filenames of files reachable by the
        HTML code, or the names can be URLs.
      * A filename generator using Unix wildcard notation, e.g.,
        ``frame*.png`` (the files most be accessible for the HTML code).
      * A filename generator using printf notation for frame numbering
        and limits for the numbers. An example is ``frame%0d.png:0->92``,
        which means ``frame00.png``, ``frame01.png``, ..., ``frame92.png``.
        This specification of `plotfiles` also allows URLs, e.g.,
        ``http://mysite.net/files/frames/frame_%04d.png:0->320``.

    If `casename` is None, a casename based on the full relative path of the
    first plotfile is used as tag in the variables in the javascript code
    such that the code for several movies can appear in the same file
    (i.e., the various code blocks employ different variables because
    the variable names differ).

    The returned result is text strings that incorporate javascript to
    loop through the plots one after another.  The html text also features
    buttons for controlling the movie.
    The parameter `iterval_ms` is the time interval between loading
    successive images and is in milliseconds.

    The `width` and `height` parameters do not seem to have any effect
    for reasons not understood.

    The following strings are returned: header, javascript code, form
    with movie and buttons, footer, and plotfiles::

       header, jscode, form, footer, plotfiles = html_movie('frames*.png')
       # Insert javascript code in some HTML file
       htmlfile.write(jscode + form)
       # Or write a new standalone file that act as movie player
       filename = plotfiles[0][:-4] + '.html'
       htmlfile = open(filename, 'w')
       htmlfile.write(header + jscode + form + footer)
       htmlfile.close

    This function is based on code written by R. J. LeVeque, based on
    a template from Alan McIntyre.
    """
    # Alternative method:
    # http://stackoverflow.com/questions/9486961/animated-image-with-javascript

    # Start with expanding plotfiles if it is a filename generator
    if not isinstance(plotfiles, (tuple,list)):
        if not isinstance(plotfiles, (str,unicode)):
            raise TypeError('plotfiles must be list or filename generator,\
                            not %s' % type(plotfiles))

        filename_generator = plotfiles
        if '*' in filename_generator:
            # frame_*.png
            if filename_generator.startswith('http'):
                raise ValueError(   'Filename generator %s cannot contain *;\
                                    must be like\
                                    http://some.net/files/frame_%%04d.png:0->120'\
                                    % filename_generator)

            plotfiles = glob.glob(filename_generator)
            if not plotfiles:
                raise ValueError('No plotfiles on the form' %
                                 filename_generator)
            plotfiles.sort()
        elif '->' in filename_generator:
            # frame_%04d.png:0->120
            # http://some.net/files/frame_%04d.png:0->120
            p = filename_generator.split(':')
            filename = ':'.join(p[:-1])
            if not re.search(r'%0?\d+', filename):
                raise ValueError(   'Filename generator %s has wrong syntax;\
                                    missing printf specification as in\
                                    frame_%%04d.png:0->120' % filename_generator)
            if not re.search(r'\d+->\d+', p[-1]):
                raise ValueError(   'Filename generator %s has wrong syntax;\
                                    must be like frame_%%04d.png:0->120'\
                                    % filename_generator)
            p = p[-1].split('->')
            lo, hi = int(p[0]), int(p[1])
            plotfiles = [filename % i for i in range(lo,hi+1,1)]

    # Check that the plot files really exist, if they are local on the computer
    if not plotfiles[0].startswith('http'):
        missing_files = [fname for fname in plotfiles
                         if not os.path.isfile(fname)]
        if missing_files:
            raise ValueError('Missing plot files: %s' %
                             str(missing_files)[1:-1])

    if casename is None:
        # Use plotfiles[0] as the casename, but remove illegal
        # characters in variable names since the casename will be
        # used as part of javascript variable names.
        casename = os.path.splitext(plotfiles[0])[0]
        # Use _ for invalid characters
        casename = re.sub('[^0-9a-zA-Z_]', '_', casename)
        # Remove leading illegal characters until we find a letter or underscore
        casename = re.sub('^[^a-zA-Z_]+', '', casename)

    filestem, ext = os.path.splitext(plotfiles[0])
    if ext == '.png' or ext == '.jpg' or ext == '.jpeg' or ext == '.gif':
        pass
    else:
        raise ValueError(   'Plotfiles (%s, ...) must be PNG, JPEG, or GIF\
                            files with extension .png, .jpg/.jpeg, or .gif'\
                            % plotfiles[0])

    header = """\
<html>
<head>
</head>
<body>
"""
    no_images = len(plotfiles)
    jscode = """
<script language="Javascript">
<!---
var num_images_%(casename)s = %(no_images)d;
var img_width_%(casename)s = %(width)d;
var img_height_%(casename)s = %(height)d;
var interval_%(casename)s = %(interval_ms)d;
var images_%(casename)s = new Array();

function preload_images_%(casename)s()
{
   t = document.getElementById("progress");
""" % vars()

    i = 0
    for fname in plotfiles:
        jscode += """
   t.innerHTML = "Preloading image ";
   images_%(casename)s[%(i)s] = new Image(img_width_%(casename)s, img_height_%(casename)s);
   images_%(casename)s[%(i)s].src = "%(fname)s";
        """ % vars()
        i = i+1
    jscode += """
   t.innerHTML = "";
}

function tick_%(casename)s()
{
   if (frame_%(casename)s > num_images_%(casename)s - 1)
       frame_%(casename)s = 0;

   document.name_%(casename)s.src = images_%(casename)s[frame_%(casename)s].src;
   frame_%(casename)s += 1;
   tt = setTimeout("tick_%(casename)s()", interval_%(casename)s);
}

function startup_%(casename)s()
{
   preload_images_%(casename)s();
   frame_%(casename)s = 0;
   setTimeout("tick_%(casename)s()", interval_%(casename)s);
}

function stopit_%(casename)s()
{ clearTimeout(tt); }

function restart_%(casename)s()
{ tt = setTimeout("tick_%(casename)s()", interval_%(casename)s); }

function slower_%(casename)s()
{ interval_%(casename)s = interval_%(casename)s/0.707; }

function faster_%(casename)s()
{ interval_%(casename)s = interval_%(casename)s*0.707; }

// --->
</script>
""" % vars()
    plotfile0 = plotfiles[0]
    form = """
<form>
&nbsp;
<input type="button" value="Start movie" onClick="startup_%(casename)s()">
<input type="button" value="Pause movie" onClick="stopit_%(casename)s()">
<input type="button" value="Resume movie" onClick="restart_%(casename)s()">
&nbsp;
<input type="button" value="Slower" onClick="slower_%(casename)s()">
<input type="button" value="Faster" onClick="faster_%(casename)s()">
</form>

<p><div ID="progress"></div></p>
<img src="%(plotfile0)s" name="name_%(casename)s" border=2/>
""" % vars()
    footer = '\n</body>\n</html>\n'
    return header, jscode, form, footer, plotfiles






***************************************************************************
waves.py
***************************************************************************

# waves.py (Python 3.4)
# Jesse M. Kinder --- 2015 Feb 26
""" Generate frames for an animation of a moving Gaussian waves. """
import numpy as np
import matplotlib.pyplot as plt
# get html_movie module from physicalmodelingwithpython.blogspot.com/
from html_movie import movie 

# Generate waves for each frame.
# Return a gaussian with specified center and spread using array s.
def gaussian(s, center=0.0, spread=1.0):
    return np.exp(-2 * (s - center)**2 / spread**2)

# All lengths are in [m], all times are in [s], and all speeds are in [m/s].
# Define the range of values to display.
xmin, xmax  = -4.0, 4.0
ymin, ymax  = -3.0, 3.0
# Define array of positions.
dx = 0.01
x = np.arange(xmin, xmax + dx, dx)

# Define the duration and number of frames for the simulation.
tmin, tmax  = 0.0, 4.0
num_frames = 100
t = np.linspace(tmin, tmax, num_frames)

# Define the initial position and speed of gaussian waves.
r_speed = 2.0       # speed of right-moving wave
r_0 = -4.0          # initial position of right-moving wave
l_speed = -2.0      # speed of left-moving wave
l_0 = 4.0           # initial position of left-moving wave

# Generate a figure and get access to its axes object.
plt.close('all')
fig = plt.figure(figsize=(10,10))
ax = plt.axes(xlim=(xmin, xmax), ylim=(ymin, ymax))

# Create three empty line objects and grab control.
# The loop below will update the lines in each frame.
ax.plot([], [], 'b--', lw=1)        # line for right-moving wave
ax.plot([], [], 'r--', lw=1)        # line for left-moving wave
ax.plot([], [], 'g-', lw=3)         # line for sum of waves
lines = ax.get_lines()              # get list of 3 line objects in plot

# It is essential that the frames be named in alphabetical order.
# {:03d} will display integers with three digits and insert zeros if needed:
# '000_movie.jpg', '001_movie.jpg', etc.
file_name = "{:03d}_movie.jpg"

# Generate frames and save each figure as a separate .jpg file.
for i in range(num_frames):
    r_now = r_0 + r_speed * t[i]        # update center of right-moving wave
    l_now = l_0 + l_speed * t[i]        # update center of left-moving wave
    yR =  gaussian(x, r_now)            # get current data for right-moving wave
    yL = -gaussian(x, l_now)            # get current data for left-moving wave
    lines[0].set_data(x,yR)             # update right-moving wave
    lines[1].set_data(x,yL)             # update left-moving wave
    lines[2].set_data(x,yR + yL)        # update sum of waves
    plt.savefig(file_name.format(i))    # save current plot

# Use html movie encoder adapted from scitools to create an HTML document that
# will display the frames as a movie.  Open movie.html in Web browser to view.
movie(input_files='*.jpg', output_file='movie.html')






***************************************************************************
sympy_examples.py
***************************************************************************

# sympy_examples.py
# -------------------------------------------------------------------------
# Demonstrate some useful methods available in the SymPy module.
# ------------------------------------------------------------------------- 
from sympy import *
init_session()

expand( (x + y)**5 )
factor( x**6 - 1 )
pi.n(100)
plot(besselj(0, x), besselj(1, x), (x, 0, 10))
diff( x*sin(y), x, y )
integrate(cos(x)**2, x)
integrate(exp(-x**2), (x, -oo, oo))
Sum(k**3, (k, 0, m)).doit().factor()
dsolve(f(x).diff(x) + f(x) - x, f(x)).simplify()






***************************************************************************
convolution.py
***************************************************************************

# convolution.py
# -------------------------------------------------------------------------
# This script creates an eLoG (elongated Laplacian of Gaussian) filter that
# emphasizes long, vertical lines in a figure.  The effect of the filter is
# demonstrated on a plus sign.
# ------------------------------------------------------------------------- 
import numpy as np
import matplotlib.pyplot as plt
import scipy.ndimage as sim

#%% Create a grid of points for the Gaussian filter.
v = np.arange(-25, 26)
X, Y = np.meshgrid(v, v)

#%% Create elongated Gaussian filter, apply Laplacian.
gauss_filter = np.exp(-0.5*(X**2/2 + Y**2/45))
laplace_filter = np.array( [[0, -1, 0], [-1, 4, -1], [0, -1, 0]])
combined_filter  = sim.convolve(gauss_filter, laplace_filter)

#%% Create a plus sign '+' to demonstrate effect of filter.
plus = np.zeros((51, 51))
plus[23:28, 25] = 1.0
plus[25, 23:28] = 1.0

plt.figure()
plt.imshow(plus)
plt.gray()

#%% Apply filter to '+' and display resulting image.
#   Use vmin/vmax to emphasize features within a restricted range of intensity.
cplus = sim.convolve(plus, combined_filter)

plt.figure()
plt.imshow(cplus, vmin=0, vmax=0.5*cplus.max())
plt.gray()

plt.show()






***************************************************************************
scope.py
***************************************************************************

# scope.py 
# -------------------------------------------------------------------------
# Demonstrate Python's rules of scope.
# ------------------------------------------------------------------------- 

def illustrate_scope():
    s_enclosing = 'E'
    def display():
        s_local = 'L'
        print( "Local --> {}".format(s_local) )
        print( "Enclosing --> {}".format(s_enclosing) )
        print( "Global --> {}".format(s_global) )
        print( "Built-in --> {}".format(abs) )
    display()

s_global = 'G'

illustrate_scope()






***************************************************************************
name_collision.py
***************************************************************************

# name_collision.py
# -------------------------------------------------------------------------
# Demonstrate how Python's rules of scope prevent name collisions.
# ------------------------------------------------------------------------- 

def name_collisions():
    x, y = 'E', 'E'
    def display():
        x = 'L'
        print( "Inside display() ..." )
        print( "x= {}\ny= {}\nz= {}".format(x,y,z) )
    display()
    print( "Inside name_collision() ..." )
    print( "x= {}\ny= {}\nz= {}".format(x,y,z) )

x, y, z = 'G', 'G', 'G'
name_collisions()
print( "Outside function ..." )
print( "x= {}\ny= {}\nz= {}".format(x,y,z) )






***************************************************************************
fancy_plot.py
***************************************************************************

# fancy_plot.py
# -------------------------------------------------------------------------
# Add a title and axis labels to a simple plot.
# ------------------------------------------------------------------------- 
import numpy as np
import matplotlib.pyplot as plt

x_min, x_max = -4, 4
num_points = 51
x_list = np.linspace(x_min, x_max, num_points)
y_list = x_list**2

plt.figure()
plt.plot(x_list, y_list, 'r', linewidth=3)

ax = plt.gca()
ax.set_title("A Second Order Polynomial", fontsize=16)
ax.set_xlabel("$x$", fontsize=24)
ax.set_ylabel("$y = x^2$", fontsize=24)

plt.show()






***************************************************************************
legend.py
***************************************************************************

# legend.py
# -------------------------------------------------------------------------
# Create a plot with a legend to distinguish multiple curves.
# ------------------------------------------------------------------------- 
import numpy as np
import matplotlib.pyplot as plt

num_curves = 3
x = np.linspace(0, 1, 51)
y = np.zeros((x.size, num_curves))
for n in range(num_curves):
    y[:, n] = np.sin((n+1) * x * 2 * np.pi)

plt.figure()
plt.plot(x, y, linewidth=2)

ax = plt.gca()
ax.legend(  ("sin(2$\\pi$x)", "sin(4$\\pi$x)", "sin(6$\\pi$x)") )

plt.show()






***************************************************************************
measurement.py
***************************************************************************

# measurement.py
"""
Functions to calculate distance between points using different metrics.
"""

import numpy as np

def crow(pointA, pointB):
    """
    Distance between points A and B "as the crow flies."
        pointA = (x1, y1)
        pointB = (x2, y2)
    returns sqrt( (x2-x1)**2 + (y2-y1)**2 )
    """
    interval = np.sqrt( (pointA[0] - pointB[0])**2 + \
                        (pointA[1] - pointB[1])**2 )
    return interval


def taxicab(pointA, pointB):
    """
    Distance between points A and B "as the cab drives."
        pointA = (x1, y1)
        pointB = (x2, y2)
    returns |x2-x1| + |y2-y1|
    """
    interval =  abs(pointB[0] - pointA[0]) + \
                abs(pointB[1] - pointA[1])
    return interval

def distance(pointA, pointB=(0,0), method='taxi'):
    """
    Return distance between points A and B. If method is 'taxi', use taxicab
    metric. Otherwise, use Euclidean distance.
        pointA = (x1, y1)
        pointB = (x2, y2)
    """
    if method == 'taxi':
        return taxicab(pointA, pointB)
    else:
        return crow(pointA, pointB)






***************************************************************************
random_walk.py
***************************************************************************

# random_walk.py
# -------------------------------------------------------------------------
# Monte Carlo simulation of a two-dimensional random walk.
# ------------------------------------------------------------------------- 
import numpy as np
import matplotlib.pyplot as plt
from numpy.random import random as rng

num_steps = 500
x_step = rng(num_steps) > 0.5
y_step = rng(num_steps) > 0.5
x_step = 2*x_step - 1
y_step = 2*y_step - 1
x_position = np.cumsum(x_step)
y_position = np.cumsum(y_step)

plt.figure()
plt.plot(x_position, y_position)
plt.axis('equal')

# A more succinct alternative:
x = ( 2*(rng(num_steps) > 0.5) - 1 ).cumsum()
y = ( 2*(rng(num_steps) > 0.5) - 1 ).cumsum()

plt.figure()
plt.plot(x,y)
plt.axis('equal')

plt.show()






***************************************************************************
surface.py
***************************************************************************

# surface.py
# -------------------------------------------------------------------------
# Create a three-dimensional surface plot.
# ------------------------------------------------------------------------- 
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# define grid of points
points = np.linspace(-1, 1, 101)
X, Y = np.meshgrid(points, points)
Z = X**2 + Y**2

# create and display surface plot
ax = Axes3D(plt.figure())
ax.plot_surface(X, Y, Z)

# Adjust rstride and cstride to use more or fewer points.
# The following command will use all points in X, Y, and Z:
# ax.plot_surface(X, Y, Z, rstride=1, cstride=1)

plt.show()






***************************************************************************
surprise.py
***************************************************************************

# surprise.py
# -------------------------------------------------------------------------
# This script will create a familar but interesting image.
# It may take about a minute to run.
# ------------------------------------------------------------------------- 
import numpy as np
import matplotlib.pyplot as plt

# This determines the number of colors used in the plot.
# The larger the value, the longer the script will take.
max_iterations = 32

# These parameters define the boundaries of the plot
x_min, x_max = -2.5, 1.5
y_min, y_max = -1.5, 1.5

# This parameter controls the grid spacing.  A smaller value gives better
# resolution, but the script will take longer to run.
ds = 0.002

X = np.arange(x_min, x_max + ds, ds)
Y = np.arange(y_min, y_max + ds, ds)
data = np.zeros( (X.size, Y.size), dtype='uint')

for i in range(X.size):
    for j in range(Y.size):
        x0, y0 = X[i], Y[j]
        x, y = x0, y0
        count = 0
        while count < max_iterations:
            # Update x and y simultaneously.
            x, y = (x0 + x*x - y*y, y0 + 2*x*y)
            # Exit loop if (x,y) is too far from the origin.
            if (x*x + y*y) > 4.0: break
            count += 1
        data[i, j] = max_iterations - count

plt.imshow(data.transpose(), interpolation='nearest', cmap='jet')
plt.axis('off')
plt.show()






No comments:

Post a Comment

To avoid duplication and inappropriate content, your comment will be reviewed before being published. Thank you for your input!


-- Jesse & Phil