From: W. Trevor King Date: Fri, 22 Feb 2013 16:19:51 +0000 (-0500) Subject: Update scatter package for assignment 5 X-Git-Url: http://git.tremily.us/?a=commitdiff_plain;h=f5c110992e6c2a7d13e1a972a032c99ebe5a7c70;p=assignment-template.git Update scatter package for assignment 5 --- diff --git a/Makefile b/Makefile index 0cc5151..ec3bc3e 100644 --- a/Makefile +++ b/Makefile @@ -21,12 +21,14 @@ # $ 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) @@ -65,8 +67,8 @@ TAR = tar 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 @@ -138,34 +140,61 @@ $(CXX_PROGRAMS): % : $$($$(*)_OBJECTS) # 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 diff --git a/README b/README index e63b768..dd3fc98 100644 --- a/README +++ b/README @@ -1,10 +1,10 @@ -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 ======== diff --git a/angle-vs-y0-zoom.gp b/angle-vs-y0-zoom.gp new file mode 100644 index 0000000..218f6e4 --- /dev/null +++ b/angle-vs-y0-zoom.gp @@ -0,0 +1,8 @@ +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 diff --git a/angle-vs-y0.gp b/angle-vs-y0.gp new file mode 100644 index 0000000..2de36c5 --- /dev/null +++ b/angle-vs-y0.gp @@ -0,0 +1,8 @@ +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 diff --git a/energy.gp b/energy.gp deleted file mode 100644 index 8137aa6..0000000 --- a/energy.gp +++ /dev/null @@ -1,11 +0,0 @@ -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 diff --git a/error.gp b/error.gp deleted file mode 100644 index 7975aee..0000000 --- a/error.gp +++ /dev/null @@ -1,11 +0,0 @@ -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' diff --git a/hw4.txt b/hw4.txt deleted file mode 100644 index e477e45..0000000 --- a/hw4.txt +++ /dev/null @@ -1,76 +0,0 @@ -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. diff --git a/hw5.txt b/hw5.txt new file mode 100644 index 0000000..cbe3742 --- /dev/null +++ b/hw5.txt @@ -0,0 +1,68 @@ +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. diff --git a/phase.gp b/phase.gp deleted file mode 100644 index c6861bc..0000000 --- a/phase.gp +++ /dev/null @@ -1,11 +0,0 @@ -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 diff --git a/plot_image.py b/plot_image.py new file mode 100755 index 0000000..ef8f057 --- /dev/null +++ b/plot_image.py @@ -0,0 +1,314 @@ +#!/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 + + """ + 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() diff --git a/one_gaussian_bump.c b/scatter.c similarity index 50% rename from one_gaussian_bump.c rename to scatter.c index 081ab20..ed164ae 100644 --- a/one_gaussian_bump.c +++ b/scatter.c @@ -19,31 +19,33 @@ along with this program. If not, see . #include /* for exit(), EXIT_*, malloc(), free() */ #include /* for stderr, *printf() */ -#include /* for memcpy() */ +#include /* for memcpy(), memset() */ #include /* for getopt_long() */ #include /* for assert() */ -#include /* for exp(), sqrt() */ +#include /* 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 { @@ -60,27 +62,50 @@ typedef double (*energy_fn_t)(state_t *state, void *system); /* 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) @@ -91,42 +116,57 @@ double gaussian_energy(gaussian_state_t *state, double x, double y) * 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; } @@ -234,11 +274,13 @@ void rk4_step(state_t *state, system_t *system, double dt) * 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]); @@ -249,22 +291,47 @@ void print_state(state_t *state, system_t *system) 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"); @@ -275,7 +342,7 @@ void parse_args(int argc, char **argv, double *dt, double *t_max, 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': @@ -285,54 +352,142 @@ void parse_args(int argc, char **argv, double *dt, double *t_max, 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; } diff --git a/time-vs-y0-zoom.gp b/time-vs-y0-zoom.gp new file mode 100644 index 0000000..52ff47e --- /dev/null +++ b/time-vs-y0-zoom.gp @@ -0,0 +1,8 @@ +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 diff --git a/time-vs-y0.gp b/time-vs-y0.gp new file mode 100644 index 0000000..cba9802 --- /dev/null +++ b/time-vs-y0.gp @@ -0,0 +1,8 @@ +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 diff --git a/trajectories.gp b/trajectories.gp new file mode 100644 index 0000000..afe2e65 --- /dev/null +++ b/trajectories.gp @@ -0,0 +1,25 @@ +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