.. _chap-advanced-plotting: =================== Advanced plotting =================== [status: mostly-written] .. rubric:: Motivation A quick plot is OK when you are exploring a data set or a function. But when you present your results to others you need to prepare the plots much more carefully so that they give the information to someone who does not know all the background you do. .. rubric:: Prerequisites * The 10-hour "serious programming" course. * The "Data files and first plots" mini-course in :numref:`chap-data-files-and-first-plots`. * The "Intermediate plotting" mini-course in :numref:`chap-intermediate-plotting`. .. rubric:: Plan We will go through some topics in plotting that might come up when you need to prepare a *polished* plot for publication, or an animation for the accompanying materials in a publication. This chapter will have a collection of examples that you might adapt to use in your papers or supporting web materials. They will all be based on the physical system of the *harmonic oscillator*, in this case a mass bouncing on a spring. (The harmonic oscillator equations are discussed and solved in more detail in :numref:`chap-differential-equations`.) After showing simple oscillator line plots we will work through examples of adding touches that convey information about the axes and lines, as well as showing example captions for plots. Then we will add another dimension: line plots that depend on a parameter can be viewed in interesting ways that give insight. It's hard to go beyond a surface plot, but there are some techniques to visualize more complex systems. Techniques that do this are volume rendering, isosurfaces and animations, and we will give examples of these. Our setup ========= The same simple harmonic oscillator equations can describe very many physical systems. Ours will be a mass bouncing sideways on a spring as shown in :numref:`fig-plotting-advanced-SHO`. .. _fig-plotting-advanced-SHO: .. figure:: ../differential-equations/Simple-harmonic-oscillator.* :width: 40% A mass bouncing on a spring. The zero position on the x axis is when the spring is at rest and does not push or pull the mass. When we teach this in introductory physics classes we assume that there is no friction as the mass bounces to the right and left. We call that the *simple harmonic oscillator* (SHO). If there is friction it is a *damped harmonic oscillator*. The spring has a *spring constant* :math:`k` which represents how stiff it is. The force that returns the spring toward the rest position after you pluck it a distance :math:`x` is :math:`F = -k \times x`. The mass is :math:`m`. We are interested in the *position* of this mass as a function of *time*. The equation for this is: .. math:: :label: SHO-solution x(t) = A \cos(\omega t) This equation has two variables x and t, and two constants, A and :math:`\omega`. :math:`A` is the *amplitude* of the oscillation: the maximum displacement :math:`x` will ever reach, both to the right and to the left. The angular frequency of the cosine wave is given by :math:`\omega = \sqrt{\frac{k}{m}}`. A few simple things we can say about this formula are: * The formula comes from solving Newton's law: :math:`F = m a` where the force is :math:`-k x`, so we have :math:`a + \frac{k}{m} x = 0`. Since the accelartion :math:`a` is the second derivative of the position :math:`\frac{d^2x}{dt^2}` or :math:`\ddot{x}`, so we have .. math:: \frac{d^2x}{dt^2} + \frac{k}{m} x = 0 which is a differential equation which is solved by equation :eq:`SHO-solution`. I don't discuss how to solve differential equations here (see :numref:`chap-differential-equations`) since we are only interested in the plotting of solutions, so you can ignore that: we will focus on plotting the solution. I only mention the equation it comes from so you can start getting excited about the amazing world of mathematical physics that lies ahead of you. .. note:: For teachers: this mention of of differential equations should be very rapid, or you could skip it altogether. Mentioning it gives a vision of depth, but you need to slip it in casually and without spending too much time on it, or it can be daunting for some students. You could say "you will see this toward the end of high school or in college; for now just look at it and let's move on to the solution." * If you pluck the spring to the right and then let it go, the amplitude :math:`A` is how far you plucked it initially. * The initial velocity is zero: you just let it go, without a push in either direction. * The frequency :math:`f = \omega / (2 \pi)` is greater with a stiffer spring, and smaller with a larger mass. The period is the inverse of the frequency: :math:`T = 1/f = 2 \pi / \omega`. A first line plot ================= Let us now jump in to generating data for this system and plot it. A very simple program will generate the trajectory. Let us pick amplitude, spring constant and mass: :math:`A = 3, k = 2, m = 1`, and we have :math:`x(t) = 3 \cos(2 t)`. This program will generate the data: .. literalinclude:: generate-SHO.py :language: python :caption: generate-SHO.py - generate data for the simple harmonic oscillator solution. Run this program with .. code-block:: console $ chmod +x generate-SHO.py $ ./generate-SHO.py > SHO.dat You can then plot it with something like: .. code-block:: console gnuplot> plot 'SHO.dat' using 1:2 with linespoints What you see is a plot of :math:`x(t)` versus :math:`t`. This plot gives you some insight: it shows you that if you pluck the spring to the right by 3 units it will bounce back and forth regularly, with an oscillation period of about 4 and a half units of time (more precisely :math:`T = 2 \pi / \sqrt{2} \approx 4.44`. This plot is OK for you to get quick insight into the output of your ``generate-SHO.py`` program. But you would certainly not want to use it to present your results to someone else because of several problems: * There is no title in the plot. * The x (time) and y (x(t)) axes are not labeled. * The line has a very generic legend. It probably says something like ``'SHO.dat' using 1:2``, which is not informative. In some situations this default plot might have other problems, like: * The ``+`` signs for points might not be what you want (maybe a filled circle, or a hollow box are better). * You might want a bit of padding on the y axis since the plot touches the top. * You might want to change the *aspect ratio* of the plot. The default is a rectangular overall shape, but sometimes you might want a different shape of rectangle, and sometimes you might want a square. Examples of a square aspect ratio are in :numref:`sec-the-drunkards-walk`. * There may be too many numbers on the x or y axes, so that they crowd in to each other. .. _sec-polishing-line-plots: Polishing line plots ==================== Let us address the first two issues. We can choose labels for the x and y axes like this: .. code-block:: console gnuplot> set xlabel 'time t (seconds)' gnuplot> set ylabel 'spring displacement x(t) (cm)' gnuplot> plot 'SHO.dat' using 1:2 with linespoints With these instructions I have also communicated the units of measure I am using for these plots. Without units the plots are junk. As for the legend for the line, we use the *title* directive in gnuplot's plot command. Something like: .. code-block:: text gnuplot> set xlabel 'time t (seconds)' gnuplot> set ylabel 'spring displacement x(t) (cm)' gnuplot> plot 'SHO.dat' using 1:2 with linespoints title 'simple harmonic oscillator' We can also give the plot an overall title. This is not too interesting when we have a single line drawing, but if we have several lines then an overall title is useful. .. code-block:: text gnuplot> set title 'Position versus time for harmonic oscillators' gnuplot> set xlabel 'time t (seconds)' gnuplot> set ylabel 'spring displacement x(t) (cm)' gnuplot> plot 'SHO.dat' using 1:2 with linespoints title 'simple harmonic oscillator' Gnuplot has very many options you can manipulate, and there is a vast collection of demos at http://gnuplot.sourceforge.net/demo/ -- it is quite rewarding to go through all the demos and adapt them to the things you might need to plot. Increasing the challenge: parametrized solutions ================================================ The simple harmonic oscillator is not a terribly rich system, so we will add a bit more to it by looking at the *damped harmonic oscillator*. The equation for this is: .. math:: \frac{d^2x}{dt^2} + c \frac{dx}{dt} + \frac{k}{m} x = 0 where :math:`c` is a constant which expresses how much the system is damped by friction. The solution is: .. math:: :label: damped-oscillator-solution x(t) = A \times e^{-ct/(2 m)} \times \cos(\omega t) where the frequency :math:`\omega` is: .. math:: :label: damped-oscillator-omega \omega = \sqrt{k/m - c^2/(4 m^2)} Let's do a quick function plot to try to underestand what a negative exponential times a cosine function looks like. In gnuplot you can type these instructions: .. literalinclude:: sample-damped-cos.gp :language: gnuplot :caption: Gnuplot instructions for a sample damped cosine function. The result should look like: .. _fig-sample-damped-cos: .. figure:: sample-damped-cos.* :width: 50% Now that we have an idea of the shape of a damped cosine curve, let us look briefly at equations :eq:`damped-oscillator-solution` and :eq:`damped-oscillator-omega` and imagine in our minds how they might behave. Overall the shape of the damped oscillator solution should look like :numref:`fig-sample-damped-cos` The exponential factor :math:`e^{-c/(2 m)}` is what makes the oscillations get smaller, and if the friction :math:`c` is big then it will damp much faster, while if the mass is big it will damp less (more inertia). The cosine factor :math:`\cos(\omega t)` means that you have oscillations (at least until it has damped a lot), and equation :eq:`damped-oscillator-omega` tells us that the friction has made the frequency of the oscillations somewhat smaller than the natural frequency :math:`\sqrt{k/m}` of the undamped oscillator. All of this depends on the friction constant :math:`c`. But how do we visualize this curve for different values of :math:`c`? Let us start by writing a program which prints out several columns for :math:`x(t, c)` with a few different values of :math:`c`. .. literalinclude:: damped-oscillator-columns.py :language: python :caption: damped-oscillator-columns.py - generate several columns of data with different damping constants. You can run this with .. code-block:: bash chmod +x damped-oscillator-columns.py ./damped-oscillator-columns.py > damped-oscillator-columns.dat and then plot it with: .. literalinclude:: damped-oscillator-columns.gp :language: gnuplot :caption: Gnuplot instructions to plot the output columns of damped-oscillator-columns.py The result should look like: .. _fig-damped-oscillator-columns: .. figure:: damped-oscillator-columns.* Damped harmonic oscillator with separate line plots for different values of the damping constant :math:`c`. :numref:`fig-damped-oscillator-columns` gives us some information if we pick through it. The legend tells us which colored line has :math:`c=0.0`, and that line should look like there is no damping, which is correct. Then we see the line with :math:`c=0.5` and we notice that (a) the peaks of the oscillation get smaller (the damping), (b) if you follow it further out you see that the peaks don't line up with the peaks of the :math:`c=0` curve (the frequency is a bit lower than the undamped curve, which validates equation :eq:`damped-oscillator-omega`). Then we see that the lines with more damping (c = 1.0, 1.5, 2.0) damp more quickly. But we had to really pick through this plot and strain our eyes and patience to see those effects. Let us see if we can make the information about the :math:`c` parameter come out more clearly. Adding a dimension: surface plots ================================= .. literalinclude:: damped-oscillator-surface.py :language: python :caption: damped-oscillator-surface.py - generate damped harmonic oscillator data for different values of the damping constant :math:`c`. Sections of data are separated by a blank line, which is what gnuplot needs to make a surface plot. We can run this program with .. code-block:: bash chmod +x damped-oscillator-surface.py ./damped-oscillator-surface.py > damped-oscillator-surface.dat .. literalinclude:: damped-oscillator-surface.gp :language: gnuplot :caption: damped-oscillator-surface.gp - make a surface plot of the damped oscillator displacement as a function of both time and the damping constant :math:`c`. The surface plot looks like: .. _fig-damped-oscillator-surface: .. figure:: damped-oscillator-surface.* :width: 80% The main difference between :numref:`fig-damped-oscillator-surface` and the figure with several lines: with the surface you can see the behavior of x versus t for very many values of c without having to scrunch your eyes. With a single glance you can see the effect of the damping constant: a downslope in the oscillatory behavior. Adding a dimension: color ========================= Spetrograms ----------- If you want to represent the information about a voice, or music, or an animal call, or a radio signal. Think carefully of what kind of information is in this type of signal. You will need to represent the *intensity* of the signal (for example how loud the music is or how strong a signal is) for each frequency, and all of this will depend on time. One tool to visualize such signals is the *spectrogram*. Let's start with picture: .. _fig-spectrogram-sample: .. figure:: Dolphin1.* :width: 60% Spetrogram of a dolphin's vocalizations. You can see chirps, clicks and harmonizing. (From Wikipedia) In :numref:`fig-spectrogram-sample` we see that the x axis is time, the y axis is frequency, and the loudness of the sound (in decibels) is represented by the color: black and dark red are more quiet, bright yellow and white are louder. What you see visually is that the various features of the sound are seen in features of the spetrogram: the chirps are the inverted V shapes, the clicks are vertical lines, and the harmonizing is seen in the horizontal striations. Once you get used to reading spectrograms you can quickly get a feeling for not just the volume of sound as a function of time, but also for *in which frequencies* those louder sounds occur. Thus the use of color helps us get a quick-look view of the entire picture, and it helps us find distinguishing features in the data. .. _sec-spectrograms-for-standard-acoustic-files: Spectrograms for standard acoustic files ---------------------------------------- Let us download a brief violin sample, specificially a single F note. :: wget --continue http://freesound.org/data/previews/153/153595_2626346-lq.mp3 -O violin-F.mp3 This gives us the file :file:`violin-F.mp3`. It's hard to load this file in as data, but we can do a couple of conversions with the programs `ffmpeg` (which manipulates audio and video formats) and `sox`, which manipulates some audio formats and can turn audio files into human-readable columns of data. :: ffmpeg -i violin-F.mp3 violin-F.aif sox violin-F.aif violin-F.dat We can look at :file:`violin-F.dat` and see that it has three columns of data. The first is time, the second is the left stereo channel, and the third is the right stereo channel. For simplicity we will analyze the left stereo channel. We use this .. _fig-spectrogram-violin-F: .. figure:: violin-F_spec.* :width: 90% Spetrogram of a violin playing an F (Fa) note. Note the fundamental frequency at the bottom, and the many harmonics above that. .. _fig-spectrogram-tuningfork440: .. figure:: tuningfork440_spec.* :width: 90% Spetrogram of the note of a tuning fork. This produces a 440 Hz A (La) sound. Note that there are almost no harmonics: the tuning fork is designed to put out a sound close to a pure sin wave. Let us record our own signal so that we can experiment. Type :: rec myvoice.dat then speak in to it for no more than two seconds and hit control-C to abort. You now have a file with time and one or two amplitude channels. You can look at it and plot it. You can also record yourself playing a musical instrument, or clapping. To play the sample back you can type: :: play myvoice.dat and you should hear back your brief audio sample. But remember: we are scientists, not just users of gadgets, so when we think of an audio file we think of it as a *data*. You can look at the file with :: less myvoice.dat and you will see three columns of numbers. The first is time, the others are the amplitude of the sound in the left and right stereo channels. .. literalinclude:: music2spectrogram.py :language: python :caption: music2spectrogram.py - Takes an audio file and makes a spectrogram of it. You can also download it at this link: :download:`music2spectrogram.py` We can run this program with .. code-block:: bash chmod +x music2spectrogram.py ./music2spectrogam.py violin-F.dat ## or to save it to a file: ./music2spectrogram.py violin-F.dat violin-F_spectrogram.png Ideas for further exploration ----------------------------- Earthquakes : color is magnitude Loops in gnuplot ================ Look at examples from here: http://gnuplot.sourceforge.net/demo/bivariat.html in particular the fourier series. Adding a dimension: animation ============================= .. _sec-animation-in-gnuplot: Animation in gnuplot -------------------- [FIXME: incomplete] http://personalpages.to.infn.it/~mignone/Numerical_Algorithms/gnuplot.pdf also: http://www.gnuplotting.org/tag/do/ .. _sec-animation-in-matplotlib: Animation in matplotlib ----------------------- simple-animated-plot.py .. literalinclude:: simple-animated-plot.py :language: python :caption: simple-animated-plot.py -- ... simple_anim.py .. literalinclude:: simple_anim.py :language: python :caption: simple_anim.py -- ... using matplotlib's ``FuncAnimation``: .. code-block:: python import matplotlib.animation as animation import matplotlib.pyplot as plt import numpy as np fig = plt.figure() def updatefig(i): fig.clear() p = plt.plot(np.random.random(100)) plt.draw() anim = animation.FuncAnimation(fig, updatefig, 100) anim.save("/tmp/test.mp4", fps=24) And from the matplotlib docs: https://matplotlib.org/gallery/animation/simple_anim.html Further reading =============== .. seealso:: clean plots: http://triclinic.org/2015/04/publication-quality-plots-with-gnuplot/ http://gnuplot.sourceforge.net/demo/ https://www.jstor.org/stable/pdf/2683253.pdf more about sox and audio formats: http://sox.sourceforge.net/sox.html https://newt.phys.unsw.edu.au/jw/musical-sounds-musical-instruments.html object oriented approach in matplotlib: https://python4astronomers.github.io/plotting/advanced.html more on matplotlib animation https://nickcharlton.net/posts/drawing-animating-shapes-matplotlib.html cool matplotlib resources: https://realpython.com/python-matplotlib-guide/ https://python4astronomers.github.io/plotting/advanced.html