henon.pdf builds via make

.gitignore

@@ -0,0 +1,5 @@
+*.pyc
+*.aux
+*.log
+/dots.pdf
+/henon.pdf

Makefile

@@ -0,0 +1,9 @@
+henon.pdf: henon.tex dots.pdf
+	pdflatex henon
+dots.pdf: henon.py plot.py
+	python plot.py --dots dots.pdf
+
+#---------------
+# Local Variables:
+# eval: (makefile-mode)
+# End:

factor.py

@@ -0,0 +1,31 @@
+'''A short script for factoring integers that illustrates:
+
+Function definition
+Recustion
+Looping
+Command line arguments
+Modules
+
+From a command line try:
+->python factor.py test
+->python factor.py 255
+->ipython
+In [1]: import factor
+In [2]: factor.factor(42)
+
+'''
+def factor(n):
+    if n == 1:
+        return []
+    for i in range(2,n/2+2):
+        if n%i == 0:
+            return [i] + factor(n/i)
+    return [n]
+if __name__ == "__main__":
+    import sys
+    if sys.argv[1] == 'test':
+        for n in (1,2,3,4,5,12,16,25,29,27,35):
+            print('prime factors of %d are %s'%(n,factor(n)))
+    else:
+        n = int(sys.argv[1])
+        print('prime factors of %d are %s'%(n,factor(n)))

henon.py

@@ -0,0 +1,46 @@
+'''henon.py Code for studying the Henon system.  See
+http://en.wikipedia.org/wiki/H%C3%A9non_map
+
+Tools and ideas:
+
+doctest http://docs.python.org/2/library/doctest.html
+
+'''
+import numpy as np
+def map(z,       # A sequence of 2 numbers
+        a=1.4,   # Parameter of the Henon system
+        b=0.3    # Parameter of the Henon system
+        ):
+    ''' Apply the Henon map to z and return the result as a numpy
+    array.  Example:
+    
+    >>> print(map((1,1)))
+    [ 0.6  0.3]
+    '''
+    x,y = z
+    return np.array([y+1-a*x**2, b*x])
+def iterate(z, n, a=1.4, b=0.3):
+    ''' Return a sequence of n x,y pairs.  Because this is implmented
+    as an iterator.  It returns the sequential values as they are
+    used, n may be larger than the memory of the computer.
+    
+    The following test passes on my system.  Since trajectories of the
+    Henon system are unstable, ie, chaotic, the test will not pass on
+    a system that does not implment floating point calculations
+    exactly like mine.
+    
+    >>> L = list(iterate((1,1),20))
+    >>> print('%7.3f %7.3f'%L[-1])
+     -0.553   0.302
+    '''
+    for i in xrange(n):
+        x,y = z
+        z = y+1-a*x**2, b*x
+        yield z
+
+def _test():
+    import doctest
+    doctest.testmod()
+
+if __name__ == "__main__":
+    _test()

henon.tex

@@ -0,0 +1,39 @@
+\documentclass[]{article}
+\usepackage{amsmath,amsfonts}
+\usepackage{amsmath,amsfonts}
+%\usepackage{showlabels}
+\usepackage[pdftex]{graphicx,color}
+\title{Using the Henon System to Illustrate Software Tools}
+
+\author{Andrew M.\ Fraser}
+\begin{document}
+\maketitle
+\begin{abstract}
+  This is a short abstract
+\end{abstract}
+
+\section{From Wikipedia}
+\label{sec:wikipedia}
+
+From http://en.wikipedia.org/wiki/Henon\_map I copied this equation
+\begin{align*}
+   x_{n+1} &= y_n+1-a x_n^2,\\
+  y_{n+1} &= b x_n.
+\end{align*}
+
+One can use gnu-make, \texttt{Makefile} and the code in \texttt{plot.py} and
+\texttt{henon.py} to make Fig.~\ref{fig:dots}
+\begin{figure}
+  \centering
+  \resizebox{0.75\textwidth}{!}{\includegraphics{dots.pdf}}
+  \caption{Five thousand iterations of the Henon system.}
+  \label{fig:dots}
+\end{figure}
+
+
+\end{document}
+
+%%%---------------
+%%% Local Variables:
+%%% eval: (TeX-PDF-mode)
+%%% End:

plot.py

@@ -0,0 +1,88 @@
+"""plot.py makes plots.  It imports stuff to do most of the
+calculations.
+
+"""
+import sys
+import matplotlib as mpl
+import numpy as np
+from numpy import random
+import matplotlib.pyplot as plt
+
+DEBUG = False
+plot_dict = {} # Keys are those keys in args.__dict__ that ask for
+               # plots.  Values are the functions that make those
+               # plots.
+def main(argv=None):
+    import argparse
+    global DEBUG
+
+    if argv is None:                    # Usual case
+        argv = sys.argv[1:]
+
+    parser = argparse.ArgumentParser(description='Make plots')
+    parser.add_argument('--N', type=int, default=5000,
+                       help='Number of Henon points')
+    parser.add_argument('--N_relax', type=int, default=100,
+                       help='Discard N_relax iterations first')
+    parser.add_argument('--debug', action='store_true')
+    # Plot requests
+    parser.add_argument('--dots', type=argparse.FileType('w'),
+                       help='Write figure to this file')
+    args = parser.parse_args(argv)
+
+    params = {'axes.labelsize': 18,     # Plotting parameters for latex
+              'text.fontsize': 15,
+              'legend.fontsize': 15,
+              'text.usetex': True,
+              'font.family':'serif',
+              'font.serif':'Computer Modern Roman',
+              'xtick.labelsize': 15,
+              'ytick.labelsize': 15}
+    mpl.rcParams.update(params)
+    if args.debug:
+        DEBUG = True
+        mpl.rcParams['text.usetex'] = False
+    else:
+        mpl.use('PDF')
+    import matplotlib.pyplot as plt  # must be after mpl.use
+
+    # Make requested plots
+    for key in args.__dict__:
+        if key not in plot_dict:
+            continue
+        if args.__dict__[key] == None:
+            continue
+        print('work on %s'%(key,))
+        fig = plot_dict[key](plt, args)
+        if not DEBUG:
+            fig.savefig(args.__dict__[key], format='pdf')
+    return 0
+
+def dots(plt, args):
+    import henon
+    fig = plt.figure(figsize=(6,6))
+    ax = fig.add_subplot(1,1,1)
+    xy = np.empty((2,args.N))
+    i = -args.N_relax
+    for _xy_ in henon.iterate((1,1),args.N_relax+args.N):
+        if i >= 0:
+            xy[:,i] = _xy_
+        i += 1
+    ax.plot(xy[0], xy[1],'b,')
+    ax.set_ylim(-.5,.5)
+    ax.set_yticks(np.arange(-.4, .41, .2),minor=False)
+    ax.set_xticks(np.arange(-1.0, 1.1, 1.0),minor=False)
+    return fig
+plot_dict['dots'] = dots
+
+if __name__ == "__main__":
+    rv = main()
+    if DEBUG:
+        import matplotlib.pyplot as plt
+        plt.show()
+    sys.exit(rv)
+
+#---------------
+# Local Variables:
+# eval: (python-mode)
+# End: