Chapter 3Theory support

3.2.1 gnuplot
3.2.2 IDL
3.2.3 DX
3.2.4 PGPLOT

3.1 Computer algebra

Starlink provides access to computer algebra by supporting the Maple package. You might have access to Maple on your own Starlink node, but if not, you may use it on the machine star.rl.ac.uk, if you have an account there. If you are a Starlink user, you should apply for an account by mailing star@star.rl.ac.uk.

Maple allows you to enter mathematical expressions using a fairly natural syntax, substitute into them, simplify them, differentiate and (with limits) integrate them, and finally go on to graph them. As an added bonus, you can also produce output in C, Fortran and LaTeX.

For example, consider the following example.

star:nxg> maple
|\^/|     Maple V Release 4 (Rutherford Appleton Laboratory)
._|\|   |/|_. Copyright (c) 1981-1996 by Waterloo Maple Inc. All rights
\  MAPLE  /  reserved. Maple and Maple V are registered trademarks of
<____ ____>  Waterloo Maple Inc.
|       Type ? for help.
> gi := amp * exp(-(xparam^2/sa^2)/2);
2
xparam
gi := amp exp(- 1/2 -------)
2
sa

We start up Maple, and enter an expression for a gaussian. Maple makes an attempt to display the result intelligibly. If we had ended the expression with a colon rather than a semicolon, Maple would have suppressed the display. Note that undefined variables represent themselves, and that Maple knows that exp is the exponential function, so that it knows, for example, how to differentiate it.

Now define the variable xparam, and redisplay gi.

> xparam:= cos(theta)*(xc-x0);
xparam := cos(theta) (xc - x0)

> gi;
2          2
cos(theta)  (xc - x0)
amp exp(- 1/2 ----------------------)
2
sa

Then differentiate the gaussian, and assign the result to gid. The result is something you’re happy not to have had to work out yourself.

> gid := diff (gi, theta);
2          2
2                      cos(theta)  (xc - x0)
amp cos(theta) (xc - x0)  sin(theta) exp(- 1/2 ----------------------)
2
sa
gid := ----------------------------------------------------------------------
2
sa

If that the purpose of this was to do a calculation somewhere, you might want to code this expression in Fortran. Doing this by hand would be error-prone, but Maple can produce output in Fortran as well as this ‘prettyprinted’ style.

> fortran (gid,optimized);
t1 = cos(theta)
t4 = (xc-x0)**2
t6 = sa**2
t7 = 1/t6
t10 = t1**2
t15 = amp*t1*t4*t7*sin(theta)*exp(-t10*t4*t7/2)

The optimized argument tells Maple to try to produce Fortran code without repeated subexpressions. You can save this to a file with the expression fortran (gid, filename=‘gaussian.f‘, optimized); You can produce output in C as well, though because the identifierC is potentially such a common one, you must explicity load the C library first.

> C([gf=gid],optimized);
t1 = cos(theta);
t4 = pow(xc-x0,2.0);
t6 = sa*sa;
t7 = 1/t6;
t10 = t1*t1;
gf = amp*t1*t4*t7*sin(theta)*exp(-t10*t4*t7/2);

There are two things to note here. The first is that we have renamed the expression gid on the fly. The second is that the expression for t4 is not the most efficient — it is very bad to use the pow() function for raising expressions to small integer powers: much better would be t4a=xc-x0; t4=t4a*t4a;, as has happened automatically for t10.

You can also produce results in LaTeX

> latex(gid);
{\it amp}\,\cos(\theta)\left ({\it xc}-{\it x0}\right )^{2}\sin(\theta
){e^{-1/2\,{\frac {\left (\cos(\theta)\right )^{2}\left ({\it xc}-{
\it x0}\right )^{2}}{{{\it sa}}^{2}}}}}{{\it sa}}^{-2}

Maple has done the correct thing with the cosine and sine functions, and with the $\theta$ variable, and it has got all the braces matching correctly, but it has expressed the exponential as a simple e-to-the-power which will look rather ugly (as well, the exponential should be written with \mathrm{e}).

Leave Maple by giving the command quit.

SUN/107 provides an introduction to Maple, and SGP/47 is a comparison of Maple and Mathematica. Also, the Maple manual and tutorial are very clear. There is help within Maple (type ?intro), and this gives enough information to get you going. Maple’s web pages are at http://www.maplesoft.com/, but they don’t currently (December 1998) give a lot of tutorial help. See also the example Maple program in A.7.

There’s also a GUI for maple, which you can invoke with xmaple.

As a final point, don’t fall into the common trap of thinking that because you’ve produced your result using computer algebra, it must be right. This is as false of computer algebra as it is of numerical programming — be inventive in thinking of cross-checks.

3.2 Data visualisation

Starlink supports two data visualisation packages, IDL and DX, but for relatively simple graphing of simple results, gnuplot will probably produce acceptable results in rather less time.

The package SM, which is a descendent of Mongo, has numerous adherents, but I don’t propose to discuss it, partly because I’ve never used it, partly because it’s not part of Starlink’s base set and so is not available at all sites, but mostly because if gnuplot runs out of steam, you might as well go straight to IDL, which is easily available through Starlink.

SG/8, “An Introduction to Visualisation Software for Astronomy”, is an overview of visualisation systems, which mentions both IDL and DX.

3.2.1 gnuplot

Gnuplot is valuable because it’s so simple — easy things can be done easily. If you want to graph an equation or plot some data, and produce postscript output, then you’ll probably do it faster, from a standing start, than someone using one of the beefier packages. Its weaknesses are that it can’t easily do very complicated or fancy visualisations (that is, it doesn’t try to take over the world), and it deals naturally only with ASCII data in columns. It is scriptable, but I wouldn’t fancy programming anything complicated with it.

Start it up with the simple command gnuplot. Once it’s started, you can give it commands as simple as

gnuplot> plot sin(x)

to produce a plot of the sine function with default ranges, or you can give it a range and specify a line type as follows

gnuplot> plot [x=0:3.14] sin(x)*cos(x**2)**2 with impulses

Gnuplot can also graph data files. The example file gausssine.dat consists of two columns of $64×64$ values. It can be plotted in gnuplot with the commands:

gnuplot> set output ’gausssine-gnuplot.eps’
gnuplot> set terminal postscript eps
Terminal type set to ’postscript’
Options are ’eps monochrome dashed "Helvetica" 14’
gnuplot> splot ’gausssine.dat’ using 1 with lines
gnuplot> set terminal x11    # set the terminal type back to the (default) X
gnuplot> set output          # close the output file

This produces the EPS file shown in 3.1. The file has a blank line after each block of 64 numbers. Gnuplot interprets this as a signal that this is the end of a ‘row’ of a matrix (a standard gnuplot gotcha is that is must be a completely blank line, with no whitespace). You can plot the surface represented by the second column, with the clause using 2 to splot. Gnuplot can read more generally formatted data, but that’s already stepping towards advanced usage.