# $ make VERSION=2 dist
# will generate phys305-hw4-2.tar.gz
COURSE = phys305
-PACKAGE = hw4
+PACKAGE = hw5
VERSION = 1
RELEASE = $(COURSE)-$(PACKAGE)-$(VERSION)
RUN_PROGRAM =
-SCRIPTS =
-C_PROGRAMS = one_gaussian_bump error_analysis
+SCRIPTS = plot_image.py angle-vs-y0.gp time-vs-y0.gp \
+ angle-vs-y0-zoom.gp time-vs-y0-zoom.gp trajectories.gp
+
+C_PROGRAMS = scatter
CXX_PROGRAMS =
PROGRAMS = $(C_PROGRAMS) $(CXX_PROGRAMS)
all: $(PROGRAMS)
# Dependency declarations
-one_gaussian_bump_OBJECTS = one_gaussian_bump.o safe-strto.o
-one_gaussian_bump.o safe-strto.o: config.h safe-strto.h
+scatter_OBJECTS = scatter.o safe-strto.o
+scatter.o safe-strto.o: config.h safe-strto.h
# target: help - display callable targets
# Use `grep` to search this file for target comments
# gnuplot plot.gp
# where plot.gp was a gnuplot script for plotting data generated by
# RUN_PROGRAM.
-run: error.png phase.png energy.png
+run: angle-vs-y0.png time-vs-y0.png trajectories.png \
+ angle-vs-y0-E.png angle-vs-y0-E.png
.SECONDEXPANSION:
-error-%.data: $(PROGRAMS)
- echo -e "# dt\terror" > "$@"
- for dt in 0.00004 0.0002 0.001 0.005 0.025; do \
- echo -en "$$dt\t" >> "$@"; \
- ./one_gaussian_bump --time-step "$$dt" --$(@:error-%.data=%) | \
- ./error_analysis >> "$@"; \
+final-vs-y0.data: scatter
+ echo -e "# y0\ttime\tangle\tenergy" > "$@"
+ for y0 in $$(seq 0 0.0008 0.8); do \
+ echo -en "$$y0\t" >> "$@"; \
+ ./scatter --y0 "$$y0" --final >> "$@"; \
done
-error.png: error.gp error-euler.data error-midpoint.data error-rk4.data
- gnuplot "$<" > "$@"
+.SECONDEXPANSION:
+final-vs-y0-zoom.data: scatter
+ echo -e "# y0\ttime\tangle\tenergy" > "$@"
+ for y0 in $$(seq 0.2772 0.0000008 0.2780); do \
+ echo -en "$$y0\t" >> "$@"; \
+ ./scatter --y0 "$$y0" --final >> "$@"; \
+ done
+
+.SECONDEXPANSION:
+trajectories.data: scatter
+ echo -e "# t\tx\ty\tv_x\tv_y\tenergy" > "$@"
+ for y0 in 0.024 0.0312 0.277672 0.277776 0.2784 0.3024; do \
+ ./scatter --y0 "$$y0" >> "$@"; \
+ echo '' >> "$@"; \
+ echo '' >> "$@"; \
+ done
-phase.data:
- > "$@"
- for E in 0.1 0.5 1 1.5 2 2.5 3 3.5 4 4.5 4.9 5 5.1 5.5; do \
- ./one_gaussian_bump --vx0-energy "$$E" >> "$@"; \
- echo >> "$@"; \
- echo >> "$@"; \
+# for y0 in $$(seq 0.0 0.00078125 0.8); do \
+# for E in $$(seq 0.1 0.00107421875 1.2); do
+.SECONDEXPANSION:
+final-vs-y0-E.data: scatter
+ echo -e "# y0\tinitial energy\ttime\tangle\tenergy" > "$@"
+ for y0 in $$(seq 0.0 0.1 0.8); do \
+ for E in $$(seq 0.1 0.1 1.2); do \
+ echo -en "$$y0\t$$E\t" >> "$@"; \
+ ./scatter --y0 "$$y0" --vx0-energy "$$E" --final >> "$@"; \
+ done \
done
-phase.png: phase.gp phase.data
+%-vs-y0.png: %-vs-y0.gp final-vs-y0.data
+ gnuplot "$<" > "$@"
+
+%-vs-y0-zoom.png: %-vs-y0-zoom.gp final-vs-y0-zoom.data
gnuplot "$<" > "$@"
-energy.png: energy.gp phase.data
+trajectories.png: trajectories.gp trajectories.data
gnuplot "$<" > "$@"
+time-vs-y0-E.png: plot_image.py final-vs-y0-E.data
+ cut -d ' ' -f 1,2,3 final-vs-y0-E.data | python2.7 ./plot_image.py --xyz -t 'Time vs impact parameter and energy'
+
+angle-vs-y0-E.png: plot_image.py final-vs-y0-E.data
+ cut -d ' ' -f 1,2,4 final-vs-y0-E.data | python2.7 ./plot_image.py --xyz -t 'Angle vs impact parameter and energy'
+
# Pattern rule for compiling object files from C++ source
# There is an implicit rule for this in GNU make
# http://www.gnu.org/software/make/manual/html_node/Catalogue-of-Rules.html
-Assignment 4
+Assignment 5
============
-This package solves `assignment 4 for physics 305`__ (see ``hw4.txt``
+This package solves `assignment 5 for physics 305`__ (see ``hw5.txt``
for a plain-text version), plotting the Mandelbrot set.
-__ http://physics.drexel.edu/~valliere/PHYS305/phys305.w13/assignments/hw4/hw4.pdf
+__ http://physics.drexel.edu/~valliere/PHYS305/phys305.w13/assignments/hw5/hw5.pdf
Building
========
--- /dev/null
+set term png
+set termoption enhanced
+set title "Scattering angle vs. impact parameter"
+set xlabel "y_0"
+set ylabel "θ"
+unset key
+set style data points
+plot 'final-vs-y0-zoom.data' using 1:3
--- /dev/null
+set term png
+set termoption enhanced
+set title "Scattering angle vs. impact parameter"
+set xlabel "y_0"
+set ylabel "θ"
+unset key
+set style data points
+plot 'final-vs-y0.data' using 1:3
+++ /dev/null
-set term png
-set termoption enhanced
-set title 'Particle energy'
-set xlabel 'x'
-set ylabel 'V(x) + K(x)'
-unset key
-k = 10
-b = 1
-V(r) = 0.5 * k * b**2 * exp(-r**2/b**2)
-set style data lines
-plot V(x) title 'V(r)', for [i=0:13] 'phase.data' index i using 2:6
+++ /dev/null
-set term png
-set title 'Error scaling'
-set xlabel 'dt'
-set ylabel 'relative error'
-set key left top
-set logscale x
-set logscale y
-set style data linespoints
-plot 'error-euler.data' title 'Euler', \
- 'error-midpoint.data' title 'Midpoint', \
- 'error-rk4.data' title '4th order Runge-Kutta'
+++ /dev/null
-PHYS 305 - Assignment #4
-Due: Friday, February 8th
-
-Make sure your name is listed as a comment at the beginning of all your work.
-
-Purpose: Illustrate the accuracy of the Differential Equations
-solvers. Study a Phase Space Portrait. Write an analysis tool.
-
-The one Gaussian potential hill solution
-========================================
-
-Understand and check your version of the code we did in class. Call
-this code one_gaussian_bump.c. Make sure your version defaults to the
-following:
-
-* Height depth = 5, mass m = 1.0
-* E = 0.5
-* x = 3.5, v = −sqrt(2E/m)
-* tmin = 0.0, tmax = 7.0, dt = 0.001
-
-Make sure that the execution line modifiers, i.e., -d 0.001, allow you
-to reset the initial values of x and v, as well as dt (there could be
-more options).
-
-Accuracy - Order of the solver
-==============================
-
-Hint: Do not modify the code one_gaussian_bump.c in this section!
-
-You want to check the order of the differential equation solver. For
-this you need to find α in
-
- error = dt^α
-
-based on the energy which ought to be a constant. Define
-
- error = max((Ecalculated - Einitial)/Einitial)
-
-where the max means to find the maximum of the error over the entire
-time domain over which you integrate.
-
-* Write a small analysis program, error_analysis.c, which reads in
- stdin from the stdout of one_gaussian_bump.c via a pipe
-
- ./one_gaussian_bump | ./error_analysis
-
- and computes error.
-
-* Write a shell script ( bash or tcsh ) to run the codes in the pipe
-
- ./one_gaussian_bump -d0.001 | ./error_analysis
-
- over a range in dt, [ 0.00004, 0.0002, 0.001, 0.005, 0.025 ]
-
-* Use the Euler, Mid-Point and RK4 (will be explained later in class)
- methods in solving for the solution of the system of diff. eq.
-
-* Plot error for the three methods of solution versus dt in the
- appropriate way to extract the power α.
-
-Phase Space Portrait
-====================
-
-Hint: Do not modify the code one_gaussian_bump.c in this section!
-
-* Use the script written above to run ./one_gausian_bump over a range
- of initial energies, [ 0.1, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3,5, 4.0,
- 4.5, 4.9, 5.0, 5.1, 5.5 ]. Use the default dt = 0.001 and the RK4
- method of solution.
-* Record the output data from all the runs in a single (big) file.
-* Plot all the trajectories, v versus x, at once in a single plot
- (with dots).
-* Plot E versus x corresponding to the different initial energies at
- once in a single plot. Superpose a plot of the potential function
- onto the latter (with dots).
-* Comment on what you see.
--- /dev/null
+PHYS 305 - Assignment #5
+Due: Monday, February 18, 2013
+
+Make sure your name is listed as a comment at the beginning of all
+your work. Please make sure that every graph and chart you turn in is
+well-labeled with axes and a title.
+
+Purpose: Develop a physical intuition about Chaotic scattering.
+Complete all exercises on chaotic scattering in the course web page
+(reproduced below):
+
+Chaotic Scattering
+==================
+
+Use the program we wrote in class to solve this scattering
+problem. Use a scaled mass m = 1.0. Use a time step parameter dt =
+0.01. Check that the code compute each trajectory until the distance
+of the particle from the origin exceeds a distance of 4.1 in scaled
+units. Choose v₀ = 0.15, k = 1.0, b = 0.25, and compute the
+trajectory (y(t) versus x(t)) for y₀ = 0.0 and y₀ = 0.2. Use x₀ =
+-3.9. Plot the two trajectories. Plot E(t) as a check. Use the code
+potential.c as a guide to compute the latter. Check the of the
+scattering angle θ and the scattering time T.
+
+* Compute the trajectories for y₀ = 0, …, 0.8 in 1000 equal
+ increments. Plot separately the scattering angle and the scattering
+ time for these trajectories versus y₀. You should find that, for
+ some ranges in y₀, both θ and T vary smoothly with y₀. However, you
+ should also find certain irregular regions where both quantities
+ change greatly from one trajectory to the next. This sensitive
+ dependence of the outcome on small changes in the initial conditions
+ is the hallmark of chaos. The behavior of the system in the
+ irregular regions is known as chaotic scattering.
+
+* Choose one of the irregular regions in the previous question and
+ “zoom in” on it by successively narrowing the range in y₀ and
+ covering the new range in the same number (1000) of steps. Do you
+ see regular regions within the chaotic band? Continue to zoom in on
+ the new chaotic regions until you have decreased the range of y₀ by
+ a factor of 1000 below that used initially. You should find
+ self-similar structure–the same pattern of regularity and
+ irregularity keeps recurring on smaller and smaller scales.
+
+* Identify some of the impact parameter y₀ values that yield the
+ longest scattering times and plot the corresponding trajectories.
+ You may want to draw small circles centered on the three repellor
+ positions to guide the eyes. See for instance: draw circles.c.
+ Better yet is to draw the interception of a plane at the appropriate
+ energy E with the potential surface. See for instance: potential
+ contour.c
+
+* Can you identify any critical trajectories that seem to separate
+ regular from irregular motion?
+
+Energy Dependence
+=================
+
+The scattering functions, scattering angle and escape time, in reality
+depend on both the impact parameter and the energy of the incoming
+particle.
+
+Form two color images (1024x1024) of the scattering angle and escape
+time, having the incoming kinetic energy (vertical axis) and the
+impact parameter (horizontal axis) as labels. Use an impact parameter
+range b = 0.0, …, 0.8 and an energy range corresponding to 0.1 to 1.2
+the hill height. Write a small feeder code setup scatter.c to
+generate initial conditions evenly distributed within these
+ranges. These are to be fed in the differential equation solver code.
+++ /dev/null
-set term png
-set termoption enhanced
-set title 'Phase space'
-set xlabel 'x'
-set ylabel 'v_x'
-unset key
-k = 10
-b = 1
-V(r) = 0.5 * k * b**2 * exp(-r**2/b**2)
-set style data lines
-plot V(x) title 'V(r)', for [i=0:13] 'phase.data' index i using 2:4
--- /dev/null
+#!/usr/bin/env python
+
+"""Generate an image from ASCII data.
+
+Step 1: make this file executable::
+
+ chmod +x plot_image.py
+
+Step 2: pipe data into python script::
+
+ ./gen_data | ./plot_image.py -s nx,ny -c nc -t 'image title'
+
+Data should be one ASCII float per line, in the following order::
+
+ z[x1,y1]
+ z[x2,y1]
+ ...
+ z[xN,y1]
+ z[x1,y2]
+ ...
+ z[xN,y2]
+ ...
+ z[xN,yN]
+
+where `x` increases from `x1` to `xN` and `y` increases from `y1`
+through `yN`.
+
+You can use the `--xyz` option to read an alternative data format,
+which is::
+
+ x1 y1 z[x1,y1]
+ x2 y2 z[x1,y1]
+ ...
+
+If you use this method, the ordering of lines in the data file is
+irrelevant, and `Nx` and `Ny` are extracted from the data file itself.
+However, you still need to use a rectangular grid (i.e. for every
+`xi`, you need to have entries for every `yi`).
+
+For usage info, run::
+
+ ./plot_image.py --help
+
+When run with interactive output (i.e. no `-o ...` option), the
+interactive figure is displayed using `pyplot.show`, which means that
+you'll have to kill `plot_image.py` using ^C or similar [1]_.
+
+For other ideas, see the Matplotlib website [2]_.
+
+.. [1] http://matplotlib.sourceforge.net/faq/howto_faq.html#use-show
+.. [2] http://matplotlib.sourceforge.net/
+"""
+
+import optparse
+import sys
+
+import matplotlib
+import matplotlib.cm
+import matplotlib.image
+import numpy
+
+# Depending on your Matplotlib configuration, you may need to adjust
+# your backend. Do this before importing pylab, matplotlib.pyplot, or
+# matplotlib.backends.
+#matplotlib.use('Agg') # select backend that doesn't require X Windows
+#matplotlib.use('GTKAgg') # select backend that supports pyplot.show()
+
+# Alternatively, you can configure the backend using
+# ~/.matplotlib/matplotlibrc:
+# $ cat ~/.matplotlib/matplotlibrc:
+# backend : Agg
+# #backend : GTKAgg
+
+import matplotlib.pyplot
+
+
+_DOC = __doc__
+
+
+def read_data_1d(stream, nx, ny):
+ """Read in data, one entry per line.
+
+ >>> from StringIO import StringIO
+ >>> s = StringIO('\\n'.join(map(str, range(10)))+'\\n')
+ >>> X,Y,Z = read_data_1d(s, 5, 2)
+ >>> X
+ array([[0, 1, 2, 3, 4, 5],
+ [0, 1, 2, 3, 4, 5],
+ [0, 1, 2, 3, 4, 5]])
+ >>> Y
+ array([[0, 0, 0, 0, 0, 0],
+ [1, 1, 1, 1, 1, 1],
+ [2, 2, 2, 2, 2, 2]])
+ >>> Z
+ array([[ 0., 1., 2., 3., 4.],
+ [ 5., 6., 7., 8., 9.]])
+ """
+ X,Y = numpy.meshgrid(range(nx+1), range(ny+1))
+ Z = numpy.loadtxt(stream)
+ assert Z.size == nx*ny, 'Z.size = %d != %d = %dx%d' % (
+ Z.size, nx*ny, nx, ny)
+ Z = Z.reshape([x-1 for x in X.shape])
+ return (X,Y,Z)
+
+def read_data_3d(stream):
+ """Read in data, one `(x, y, z)` tuple per line.
+
+ >>> from StringIO import StringIO
+ >>> lines = []
+ >>> for x in range(5):
+ ... for y in range(2):
+ ... lines.append('\t'.join(map(str, [x, y, x+y*5])))
+ >>> s = StringIO('\\n'.join(lines)+'\\n')
+ >>> X,Y,Z = read_data_3d(s)
+ >>> X
+ array([[ 0., 1., 2., 3., 4., 5.],
+ [ 0., 1., 2., 3., 4., 5.],
+ [ 0., 1., 2., 3., 4., 5.]])
+ >>> Y
+ array([[ 0., 0., 0., 0., 0., 0.],
+ [ 1., 1., 1., 1., 1., 1.],
+ [ 2., 2., 2., 2., 2., 2.]])
+ >>> Z
+ array([[ 0., 1., 2., 3., 4.],
+ [ 5., 6., 7., 8., 9.]])
+ """
+ XYZ = numpy.loadtxt(stream)
+ assert len(XYZ.shape) == 2 and XYZ.shape[1] == 3, XYZ.shape
+ Xs = numpy.array(sorted(set(XYZ[:,0])))
+ Ys = numpy.array(sorted(set(XYZ[:,1])))
+ Z = numpy.ndarray((len(Ys), len(Xs)), dtype=float)
+ xyz = {} # dict of z values keyed by (x,y) tuples
+ for i in range(XYZ.shape[0]):
+ xyz[(XYZ[i,0], XYZ[i,1])] = XYZ[i,2]
+ for i,x in enumerate(Xs):
+ for j,y in enumerate(Ys):
+ Z[j,i] = xyz[x,y]
+ # add dummy row/column for pcolor
+ dx = Xs[-1] - Xs[-2]
+ dy = Ys[-1] - Ys[-2]
+ Xs = numpy.append(Xs, Xs[-1] + dx)
+ Ys = numpy.append(Ys, Ys[-1] + dy)
+ X,Y = numpy.meshgrid(Xs, Ys)
+ return (X,Y,Z)
+
+def plot(X, Y, Z, full=False, title=None, contours=None, interpolation=None,
+ cmap=None):
+ """Plot Z over the mesh X, Y.
+
+ >>> X, Y = numpy.meshgrid(range(6), range(2))
+ >>> Z = X[:-1,:-1]**2 + Y[:-1,:-1]
+ >>> plot(X, Y, Z) # doctest: +ELLIPSIS
+ <matplotlib.figure.Figure object at 0x...>
+ """
+ X_min = X[0,0]
+ X_max = X[-1,-1]
+ Y_min = Y[0,0]
+ Y_max = Y[-1,-1]
+
+ fig = matplotlib.pyplot.figure()
+ if full:
+ axes = fig.add_axes([0, 0, 1, 1])
+ else:
+ axes = fig.add_subplot(1, 1, 1)
+ if title:
+ axes.set_title(title)
+ axes.set_axis_off()
+
+ if contours:
+ cset = axes.contour(X[:-1,:-1], Y[:-1,:-1], Z, contours, cmap=cmap)
+ # [:-1,:-1] to strip dummy last row & column from X&Y.
+ axes.clabel(cset, inline=1, fmt='%1.1f', fontsize=10)
+ else:
+ # pcolor() is much slower than imshow.
+ #plot = axes.pcolor(X, Y, Z, cmap=cmap, edgecolors='none')
+ #axes.autoscale_view(tight=True)
+ plot = axes.imshow(Z, aspect='auto', interpolation=interpolation,
+ origin='lower', cmap=cmap,
+ extent=(X_min, X_max, Y_min, Y_max))
+ if not full:
+ fig.colorbar(plot)
+ return fig
+
+
+def get_possible_interpolations():
+ try: # Matplotlib v1.0.1
+ return sorted(matplotlib.image.AxesImage._interpd.keys())
+ except AttributeError:
+ try: # Matplotlib v0.91.2
+ return sorted(matplotlib.image.AxesImage(None)._interpd.keys())
+ except AttributeError:
+ # give up ;)
+ pass
+ return ['nearest']
+
+def test():
+ import doctest
+ results = doctest.testmod()
+ return results.failed
+
+
+def main(argv=None):
+ """Read in data and plot it.
+
+ >>> from tempfile import NamedTemporaryFile
+ >>> i = NamedTemporaryFile(prefix='tmp-input', suffix='.dat')
+ >>> i.write('\\n'.join([str(x) for x in range(10)])+'\\n')
+ >>> i.flush()
+ >>> o = NamedTemporaryFile(prefix='tmp-output', suffix='.png')
+ >>> main(['-i', i.name, '-s', '5,2', '-o', o.name, '-m', 'binary'])
+ Plot_image
+ Title: Some like it hot
+ Image size: 5 2
+ False color
+ X range: 0 4 (6 steps)
+ Y range: 0 1 (3 steps)
+ Z range: 0.0 9.0
+ >>> img = o.read()
+ >>> img.startswith('\\x89PNG')
+ True
+ >>> i.close()
+ >>> o.close()
+ """
+ if argv == None:
+ argv = sys.argv[1:]
+
+ usage = '%prog [options]'
+ epilog = _DOC
+ p = optparse.OptionParser(usage=usage, epilog=epilog)
+ p.format_epilog = lambda formatter: epilog+'\n'
+
+ p.add_option('-i', '--input', dest='input', metavar='FILE',
+ help='If set, read data from FILE rather than stdin.')
+ p.add_option('-o', '--output', dest='output', metavar='FILE',
+ help=('If set, save the figure to FILE rather than '
+ 'displaying it immediately'))
+ p.add_option('-s', '--size', dest='size', default='%d,%d' % (16, 16),
+ help='Data size (columns,rows; default: %default)')
+ p.add_option('-3', '--xyz', dest='xyz', default=False, action='store_true',
+ help=('If set, read (x,y,z) tuples from the input data rather'
+ 'then reading `z` and calculating `x` and `y` from '
+ '`--size`.'))
+ p.add_option('-c', '--contours', dest='contours', type='int',
+ help=('Number of contour lines (if not set, draw false color '
+ 'instead of contour lines; default: %default)'))
+ p.add_option('-f', '--full-figure', dest='full', action='store_true',
+ help=('Set axes to fill the figure (i.e. no title or color '
+ 'bar'))
+ p.add_option('-t', '--title', dest='title', default='Some like it hot',
+ help='Title (%default)')
+ p.add_option('--test', dest='test', action='store_true',
+ help='Run internal tests and exit.')
+ interpolations = get_possible_interpolations()
+ p.add_option('--interpolation', dest='interpolation', default='nearest',
+ help=('Interpolation scheme (for false color images) from %s '
+ '(%%default)') % ', '.join(interpolations))
+ maps=[m for m in matplotlib.cm.datad if not m.endswith("_r")]
+ maps.sort()
+ p.add_option('-m', '--color-map', dest='cmap', default='jet',
+ help='Select color map from %s (%%default)' % ', '.join(maps))
+
+ options,args = p.parse_args(argv)
+
+ if options.test:
+ sys.exit(test())
+
+ nx,ny = [int(x) for x in options.size.split(',')]
+ try:
+ cmap = getattr(matplotlib.cm, options.cmap)
+ except AttributeError:
+ raise Exception('no color map named %s in %s'
+ % (options.cmap, ', '.join(maps)))
+
+ print('Plot_image')
+ print('Title: {0}'.format(options.title))
+ if not options.xyz:
+ print('Image size: {0} {1}'.format(nx, ny))
+ if options.contours:
+ print('# countour lines: {0}'.format(options.contours))
+ else:
+ print('False color')
+ if options.input:
+ fin = open(options.input, 'r')
+ else:
+ fin = sys.stdin
+
+ if options.xyz:
+ X,Y,Z = read_data_3d(fin)
+ else:
+ X,Y,Z = read_data_1d(fin, nx, ny)
+
+ if options.input:
+ fin.close()
+
+ Z_min = numpy.min(Z.flat)
+ Z_max = numpy.max(Z.flat)
+ print('X range: {0} {1} ({2} steps)'.format(
+ X[0,0], X[0,-2], X.shape[1]))
+ print('Y range: {0} {1} ({2} steps)'.format(
+ Y[0,0], Y[-2,0], Y.shape[0]))
+ print('Z range: {0} {1}'.format(Z_min, Z_max))
+
+ fig = plot(X, Y, Z, full=options.full, title=options.title,
+ contours=options.contours, interpolation=options.interpolation,
+ cmap=cmap)
+
+ if options.output:
+ fig.savefig(options.output)
+ else:
+ matplotlib.pyplot.show()
+
+
+if __name__ == '__main__':
+ main()
#include <stdlib.h> /* for exit(), EXIT_*, malloc(), free() */
#include <stdio.h> /* for stderr, *printf() */
-#include <string.h> /* for memcpy() */
+#include <string.h> /* for memcpy(), memset() */
#include <getopt.h> /* for getopt_long() */
#include <assert.h> /* for assert() */
-#include <math.h> /* for exp(), sqrt() */
+#include <math.h> /* for exp(), sqrt(), atan2 */
#include "config.h" /* enabling safe_strtoc() */
#include "safe-strto.h" /* for safe_strtol() and safe_strtoc() */
static struct option long_options[] = {
- {"max-time", required_argument, 0, 't' },
{"time-step", required_argument, 0, 'd' },
+ {"max-time", required_argument, 0, 't' },
+ {"max-distance", required_argument, 0, 'r' },
{"x0", required_argument, 0, 'x' },
{"y0", required_argument, 0, 'y' },
{"vx0", required_argument, 0, 'X' },
- {"vx0-energy", required_argument, 0, 'E' },
+ {"vx0-energy", required_argument, 0, 'e' },
{"vy0", required_argument, 0, 'Y' },
{"mass", required_argument, 0, 'm' },
{"bump-magnitude", required_argument, 0, 'k' },
{"bump-width", required_argument, 0, 'b' },
- {"euler", no_argument, 0, 'e' },
+ {"euler", no_argument, 0, 'E' },
{"midpoint", no_argument, 0, 'M' },
- {"rk4", no_argument, 0, 'r' },
+ {"rk4", no_argument, 0, 'R' },
+ {"final", no_argument, 0, 'f' },
{}
};
-static char *short_options = "t:d:x:y:X:E:Y:m:k:b:eMr";
+static char *short_options = "d:t:r:x:y:X:e:Y:m:k:b:EMRf";
/* dynamic system parameters */
typedef struct state_struct {
/* static potential parameters */
typedef struct gaussian_state_struct {
- double k; /* magnitude */
- double b; /* width */
+ double k; /* magnitude */
+ double b; /* width */
+ double *x; /* position, size of state_x.n/2 */
} gaussian_state_t;
+typedef struct gaussians_state_struct {
+ size_t n; /* number of gaussians */
+ gaussian_state_t *gaussian; /* array of gaussians */
+} gaussians_state_t;
+
/* static system parameters */
typedef struct system_struct {
double m; /* particle mass */
dx_dt_fn_t dx_dt_fn; /* &single_gaussian_dx_dt, etc. */
energy_fn_t energy_fn; /* &single_gaussian_energy, etc. */
- void *dx_dt_state; /* gaussian_state_t, etc. */
+ void *dx_dt_state; /* &gaussians_state_t, etc. */
} system_t;
+/* simulation parameters */
/* function type for integration */
typedef void (*step_fn_t)(state_t *state, system_t *system, double dt);
+typedef struct simulation_struct {
+ double dt; /* timestep */
+ double t_max; /* maxium simulation time */
+ double r_max; /* maxium particle distance from the origin */
+ double _r_max_squared; /* to avoid an expensive sqrt() */
+ step_fn_t step_fn; /* &euler_step, etc. */
+} simulation_t;
+
/* Potential: V(r) = 1/2 kb^2 exp(-r^2/b^2) */
-double gaussian_energy(gaussian_state_t *state, double x, double y)
+double gaussian_energy(state_t *state, gaussian_state_t *gaussian)
{
- double k = state->k;
- double b = state->b;
- return 0.5 * k*b*b * exp(-(x*x+y*y)/b*b);
+ double k = gaussian->k;
+ double b = gaussian->b;
+ double r2 = 0; /* accumulate distance squared */
+ double x;
+ size_t i;
+
+ for (i = 0; i < state->n / 2; i++) {
+ x = state->x[i] - gaussian->x[i]; /* distance from gaussian */
+ r2 += x * x;
+ }
+ return 0.5 * k*b*b * exp(-r2/(b*b));
}
/* Potential: V(r) = 1/2 kb^2 exp(-r^2/b^2)
* Acceleration: a(r) = F(r)/m
* = k/m * exp(-r^2/b^2) * (x ihat + y jhat)
*/
-void gaussian_acceleration(gaussian_state_t *state, double m,
- double x, double y,
- double *dx_dt, double *dy_dt)
+void gaussian_acceleration(state_t *state, double m,
+ gaussian_state_t *gaussian, double *dx_dt)
{
- double k = state->k;
- double b = state->b;
- double mag = k/m * exp(-(x*x+y*y)/b*b);
- *dx_dt = mag*x;
- *dy_dt = mag*y;
+ double k = gaussian->k;
+ double b = gaussian->b;
+ double r2 = 0; /* accumulate distance squared */
+ double *x = (double *) malloc(sizeof(gaussian->x));
+ double mag;
+ size_t i, n = state->n / 2;
+
+ assert(x); /* check for successful malloc() */
+ for (i = 0; i < n; i++) {
+ x[i] = state->x[i] - gaussian->x[i]; /* distance from gaussian */
+ r2 += x[i] * x[i];
+ }
+ mag = k/m * exp(-r2/(b*b));
+ for (i = 0; i < n; i++)
+ dx_dt[i + n] += mag * x[i];
+ free(x);
return;
}
-void single_gaussian_dx_dt(state_t *state, void *system, double *dx_dt)
+void gaussians_dx_dt(state_t *state, void *system, double *dx_dt)
{
system_t *s = (system_t *)system;
+ gaussians_state_t *g = (gaussians_state_t *)s->dx_dt_state;
+ size_t i;
- assert(state->n == 4);
dx_dt[0] = state->x[2]; /* dx/dt = v_x */
dx_dt[1] = state->x[3]; /* dy/dt = v_y */
/* dv_x/dt = a_x
* dv_y/dt = a_y
*/
- gaussian_acceleration((gaussian_state_t *)s->dx_dt_state, s->m,
- state->x[0], state->x[1],
- &dx_dt[2], &dx_dt[3]);
+ dx_dt[2] = dx_dt[3] = 0;
+ for (i = 0; i < g->n; i++) /* add accel. from each gaussian */
+ gaussian_acceleration(state, s->m, &g->gaussian[i], dx_dt);
}
-double single_gaussian_energy(state_t *state, void *system)
+double gaussians_energy(state_t *state, void *system)
{
system_t *s = (system_t *)system;
+ gaussians_state_t *g = (gaussians_state_t *)s->dx_dt_state;
+ size_t i;
+
double kinetic, potential;
assert(state->n == 4);
kinetic = 0.5 * s->m * (state->x[2]*state->x[2] + state->x[3]*state->x[3]);
- potential = gaussian_energy((gaussian_state_t *)s->dx_dt_state,
- state->x[0], state->x[1]);
+ potential = 0;
+ for (i = 0; i < g->n; i++) /* add energy from each gaussian */
+ potential += gaussian_energy(state, &g->gaussian[i]);
return kinetic + potential;
}
* If system->energy_fn is defined, also print "\tenergy" before the
* newline.
*/
-void print_state(state_t *state, system_t *system)
+void print_state(state_t *state, system_t *system, int print_final)
{
size_t i;
double energy;
+ if (print_final)
+ return; /* don't print the trajectory */
printf("%0.10g", state->t);
for (i = 0; i < state->n; i++)
printf("\t%0.10g", state->x[i]);
printf("\n");
}
-void parse_args(int argc, char **argv, double *dt, double *t_max,
- state_t *state, system_t *system, step_fn_t *step_fn)
+/* prints "T\ttheta\n"
+ * If system->energy_fn is defined, also print "\tenergy" before the
+ * newline.
+ */
+void print_final_state(state_t *state, system_t *system, int print_final)
+{
+ size_t i;
+ double energy;
+
+ if (!print_final)
+ return; /* don't print the final state */
+ assert(state->n == 4); /* feel safe about extracting v_x and v_y */
+ printf("%0.10g\t%.10g", state->t, atan2(state->x[3], state->x[2]));
+ if (system->energy_fn) {
+ energy = (*system->energy_fn)(state, (void *)system);
+ printf("\t%0.10g", energy);
+ }
+ printf("\n");
+}
+
+void parse_args(int argc, char **argv, state_t *state, system_t *system,
+ simulation_t *simulation, int *print_final)
{
int ch;
- gaussian_state_t *gstate = (gaussian_state_t *)system->dx_dt_state;
+ gaussians_state_t *gstate = (gaussians_state_t *)system->dx_dt_state;
+ double dtemp;
+ size_t i;
while (1) {
ch = getopt_long(argc, argv, short_options, long_options, NULL);
if (ch == -1)
break;
switch(ch) {
+ case 'd':
+ simulation->dt = safe_strtod(optarg, "time-step");
+ break;
case 't':
- *t_max = safe_strtod(optarg, "max-time");
+ simulation->t_max = safe_strtod(optarg, "max-time");
break;
- case 'd':
- *dt = safe_strtod(optarg, "time-step");
+ case 'r':
+ simulation->r_max = safe_strtod(optarg, "max-distance");
break;
case 'x':
state->x[0] = safe_strtod(optarg, "x0");
case 'X':
state->x[2] = safe_strtod(optarg, "vx0");
break;
- case 'E':
+ case 'e':
state->x[2] = -sqrt(2*safe_strtod(optarg, "vx0-energy")/system->m);
break;
case 'Y':
system->m = safe_strtod(optarg, "mass");
break;
case 'k':
- gstate->k = safe_strtod(optarg, "bump-magnitude");
+ dtemp = safe_strtod(optarg, "bump-magnitude");
+ for (i = 0; i < gstate->n; i++)
+ gstate->gaussian[i].k = dtemp;
break;
case 'b':
- gstate->b = safe_strtod(optarg, "bump-width");
+ dtemp = safe_strtod(optarg, "bump-width");
+ for (i = 0; i < gstate->n; i++)
+ gstate->gaussian[i].b = dtemp;
break;
- case 'e':
- *step_fn = &euler_step;
+ case 'E':
+ simulation->step_fn = &euler_step;
break;
case 'M':
- *step_fn = &midpoint_step;
+ simulation->step_fn = &midpoint_step;
break;
- case 'r':
- *step_fn = &rk4_step;
+ case 'R':
+ simulation->step_fn = &rk4_step;
+ break;
+ case 'f':
+ *print_final = 1;
break;
case '?':
exit(EXIT_FAILURE);
default:
- printf("?? getopt returned character code 0%o ??\n", ch);
+ fprintf(stderr, "?? getopt returned character code 0%o ??\n", ch);
}
}
return;
}
+void setup(state_t *state, system_t *system, simulation_t *simulation)
+{
+ gaussians_state_t *gstate;
+ size_t i;
+
+ state->n = 4;
+ state->x = (double *) malloc(sizeof(double) * state->n);
+ assert(state->x); /* check for successful malloc() */
+ memset((void *)state->x, 0, sizeof(state->x));
+ state->x[0] = -3.9;
+ state->x[2] = 0.15;
+
+ system->m = 1;
+ system->dx_dt_fn = &gaussians_dx_dt;
+ system->energy_fn = &gaussians_energy;
+ gstate = (gaussians_state_t *) malloc(sizeof(gaussians_state_t));
+ assert(gstate); /* check for successful malloc() */
+ system->dx_dt_state = (void *)gstate;
+ gstate->n = 3;
+ gstate->gaussian = (gaussian_state_t *) malloc(sizeof(gaussian_state_t) *
+ gstate->n);
+ assert(gstate->gaussian); /* check for successful malloc() */
+ for (i = 0; i < gstate->n; i++) {
+ gstate->gaussian[i].k = 1;
+ gstate->gaussian[i].b = 0.25;
+ gstate->gaussian[i].x = (double *) malloc(sizeof(double) * state->n / 2);
+ assert(gstate->gaussian[i].x); /* check for successful malloc() */
+ memset((void *)gstate->gaussian[i].x, 0, sizeof(gstate->gaussian[i].x));
+ switch (i) {
+ case 0:
+ gstate->gaussian[i].x[1] = 0.5;
+ break;
+ case 1:
+ gstate->gaussian[i].x[1] = -0.5;
+ break;
+ case 2:
+ gstate->gaussian[i].x[0] = 0.5*sqrt(3);
+ break;
+ default:
+ fprintf(stderr, "no position for gaussian %d?\n", (int) i);
+ }
+ }
+
+ simulation->dt = 0.001;
+ simulation->t_max = 0;
+ simulation->r_max = 4.1;
+ simulation->step_fn = &rk4_step;
+}
+
+/* free any memory allocated in setup() */
+void teardown(state_t *state, system_t *system, simulation_t *simulation)
+{
+ gaussians_state_t *gstate;
+ size_t i;
+
+ if (state->x) {
+ free(state->x);
+ state->x = NULL;
+ }
+ if (system->dx_dt_state) {
+ gstate = (gaussians_state_t *)system->dx_dt_state;
+ if (gstate->gaussian) {
+ for (i = 0; i < gstate->n; i++) {
+ if (gstate->gaussian[i].x) {
+ free(gstate->gaussian[i].x);
+ gstate->gaussian[i].x = NULL;
+ }
+ }
+ free(gstate->gaussian);
+ gstate->gaussian = NULL;
+ }
+ free(gstate);
+ system->dx_dt_state = NULL;
+ }
+}
+
+int continue_simulation(state_t *state, system_t *system,
+ simulation_t *simulation)
+{
+ assert(state->n == 4);
+ return ((simulation->t_max <= 0 ||
+ state->t < simulation->t_max) &&
+ (simulation->_r_max_squared <= 0 ||
+ state->x[0]*state->x[0] + state->x[1]*state->x[1] <
+ simulation->_r_max_squared));
+}
+
int main(int argc, char **argv) {
- double t_max = 7;
- double dt = 0.001;
state_t state = {};
- double x[4] = {3.5, 0, 0, 0};
system_t system = {};
- gaussian_state_t gaussian = {};
- step_fn_t step_fn = &euler_step;
-
- state.x = (double *)x;
- state.n = sizeof(x) / sizeof(double);
- system.m = 1;
- system.dx_dt_fn = &single_gaussian_dx_dt;
- system.energy_fn = &single_gaussian_energy;
- system.dx_dt_state = (void *)&gaussian;
- gaussian.k = 10;
- gaussian.b = 1;
- state.x[2] = -sqrt(2*0.5/system.m); /* v = -sqrt(2E/m), with E = 0.5 */
-
- parse_args(argc, argv, &dt, &t_max, &state, &system, &step_fn);
-
- print_state(&state, &system);
- while (state.t < t_max) {
- (*step_fn)(&state, &system, dt);
- print_state(&state, &system);
+ simulation_t simulation = {};
+ int print_final = 0;
+
+ setup(&state, &system, &simulation);
+
+ parse_args(argc, argv, &state, &system, &simulation, &print_final);
+ simulation._r_max_squared = simulation.r_max * simulation.r_max;
+
+ print_state(&state, &system, print_final);
+ while (continue_simulation(&state, &system, &simulation)) {
+ (*simulation.step_fn)(&state, &system, simulation.dt);
+ print_state(&state, &system, print_final);
}
+ print_final_state(&state, &system, print_final);
+
+ teardown(&state, &system, &simulation);
+
return 0;
}
--- /dev/null
+set term png
+set termoption enhanced
+set title "Scattering time vs. impact parameter"
+set xlabel "y_0"
+set ylabel "T"
+unset key
+set style data points
+plot 'final-vs-y0-zoom.data' using 1:2
--- /dev/null
+set term png
+set termoption enhanced
+set title "Scattering time vs. impact parameter"
+set xlabel "y_0"
+set ylabel "T"
+unset key
+set style data points
+plot 'final-vs-y0.data' using 1:2
--- /dev/null
+set term png
+set title "Long trajectories"
+set xlabel "x"
+set ylabel "y"
+unset key
+
+r(x, y, x0, y0) = sqrt((x-x0)**2 + (y-y0)**2)
+k = 1
+b = 0.25
+Vgauss(x, y, x0, y0) = 0.5*k*b**2 * exp(-r(x,y,x0,y0)**2 / b**2)
+V(x, y) = Vgauss(x, y, 0, 0.5) + Vgauss(x, y, 0, -0.5) \
+ + Vgauss(x, y, 0.5*sqrt(3), 0)
+
+set xrange [-1.5:1.5]
+set yrange [-1.5:1.5]
+
+# set things up so we can plot a grayscale potential under our trajectories
+set view map
+set style data dots
+set samples 200; set isosamples 200
+set pm3d at b
+set palette gray
+set cbrange [0:0.5*k*b**2]
+
+splot V(x, y) with pm3d nosurface, for [i=0:5] 'trajectories.data' index i using 2:3:1
projects / assignment-template.git / commitdiff