Introduction to the Maths Lab - Part 2

The METRIC Project

Introduction

In this session we will concentrate on the graphical capabilities of Mathematica; we will also show you how data can be imported into Mathematica from other applications and how it can be exported out of Mathematica.

1. Two Dimensional Plotting

A standard task in the physical laboratory is the plotting and analysis of two dimensional experimental data. In many cases this data will be held as a sequence of numbers in a file external to Mathematica (in fact, all computer-driven lab instruments can give you results directly in this form). The numbers can be stored in a variety of formats; for this exercise the file is simply a text file holding pairs of values of time and intensity separated by a space, each pair appearing on a new line. 

To see what is inside any file you simply type in two exclamation marks followed by the path to the file. The example from our computer was.

! ! "Hard disk:Desktop Folder:expt1.dat"

The path will be different for you, and the way the path is expressed will look different on a PC. So to find expt1.dat on your computer, first start Mathematica, open up a new notebook, type in two exclamation marks and then go to  the Input menu and select Get File Path ...  to locate the data file expt1.dat on your hard disk. When you click on Open the appropriate file path will appear on the input line. Press <shift-return> now and the contents of the file will appear.

To read the data into Mathematica we use the ReadList command. The command opens up the disk file and is told to read pairs of number into a list we have called data. Try this command, using your own file path, and compare it to the the contents you read above.

data = 
      ReadList["Hard disk:Desktop Folder:expt1.dat", {Number, Number}]

We can plot this using the following command.

ListPlot[data]

The data is from a kinetic experiment measuring current in μA as a function of time in s, so we should label the plot appropriately. We also choose to join up the data points, colour the resulting line, put the plot into a frame, add a title and change the font.

\(plot = 
      ListPlot[data, Frame -> True, 
        FrameLabel -> {"t/s", "I/μA"}, 
RotateLabel -> False, 
        PlotLabel -> "Kinetic Data", DefaultFont -> {"Helvetica", 9}, 
        
PlotStyle -> RGBColor[1, 0, 0]]\)

It is anticipated that the data should follow first order kinetics; we can check this by plotting the logarithm of the current against the time. To do this we first need to load up an add-on package which allows all sorts of graphical manipulations

Needs[ "Graphics`Graphics`"]

and then plot the intensity data on logarithmic scales.

\(logplot = 
      LogListPlot[data, Frame -> True, 

        FrameLabel -> {"t/s", "I/μA"}, RotateLabel -> False, 
        PlotLabel -> "Kinetic Data", 

        DefaultFont -> {"Helvetica", 9}, 
        PlotStyle -> {RGBColor[1, 0, 0]}]\)

Exercises and Explorations 1

1: Logging the data

You can estimate the gradient and intercept from this plot and use these values in the model of the current as a function of the time, t.

current[t_] := (Subscript[I, 0]) (E^(-λ t))

Notice that you have to define the function via the "colon-equals" sign and each variable must be declared (on the left hand side) by being followed by an underscore. However, when you want to evaluate the function, for example at t=10, you simply type Current[10]; there is no underscore. Below we show how the function has been used to plot the model for rather poor guesses at the kinetic parameters.

\(logmodel = 
      LogPlot[current[t]\  /. {I\_0 -> \ 1.5, λ\  -> \  .3}, {t, 0, 
          4}, Frame -> True, 
FrameLabel -> {"t/s", "I/μA"}, 
        RotateLabel -> False, 
DefaultFont -> {"Helvetica", 9}, 
        PlotStyle -> {RGBColor[0, 0, 1]}]\)

Here we have introduced the idea of replacement rules when we asked for Current[t] to be evaluated for {(Subscript[I, 0]) ->  1.5, λ  ->   .3}. Using replacement rules means that the values for (Subscript[I, 0]) and λ are only "remembered" during the calculation; we haven't asigned permanent values to them. Adjust the parameters until you find the best fit to the experimental data; checking your guesses by superimposing data and model.

Show[logplot, logmodel]

2: Saving a Mathematica graphic

With the best fit parameters plot the data and model on linear scales and call the result, fig1. You can transfer the resulting graphical output to a disk file using the following command, but with your own file path.

fig1  >>  Hard disk : Desktop Folder : figure . dat

Check that the file has been put in the place you expected to go. Open it up, by double-clicking on the icon, and you should see that instead of the plot the file contains the Mathematica code that produced it. Now close the file and return to your own notebook. To read the file back in and display it as a graphic is simple.

fig1again = <<  figure.dat

Show[fig1again]

3: Exporting a Mathematica plot

What you have done above is all fine and well if you want to do is save a figure in a compact form. However, Mathematica's great power lies, not in its plotting abilities, but in its mathematical prowess. You may therefore want to export the complex functions you have generated with Mathematica to more convenient plotting programs like Excel, SigmaPlot or Kaleidagraph. There are a number of ways of doing this, but the most straightforward generates a simple data file which such programs can read. 

The code to generate a file containing the time and the current generated by the theoretical model is shown below.

OutputForm[
        TableForm[
          Table[{FortranForm[t], FortranForm[current[t]]}  /. {
                (Subscript[I, 0]) ->  1.5, λ  ->   .3}, {t, 0, 4, 0.1}]]]  >> 
       "Hard disk:Desktop Folder:mytheory.dat"

If you do not use the function OutputForm what is exported to mytheory.dat? What does FortranForm do and under what circumstances would it matter if you left it out of the code? (Hint: try to evaluate FortranForm[2.50 10^6].)

There are many other forms of input and output and you may wish to look some of these up under Help. There are also very useful examples and ideas to be found at the Wolfram website; http://www.wolfram.co.uk/support/InputOutput/ExternalFiles/. However, for now we turn to the handling of three dimensional plots. Before doing this save your work from this part of the lab and quit the kernel, before opening a new Notebook to work on Section 2.

2. Surface Plots

A surface plot arises when we have a function z=f(x,y); that is, the height z of the surface is a function of position in the x-y plane.  A good chemical example of this arises when we consider the quantum mechanical wave function of a particle confined to a rectangular planar surface of side (Subscript[l, 1])in the x direction and (Subscript[l, 2])in the y direction. (This is a reasonable model of the wave function for a molecule adsorbed on to a crystal surface.) The wave function, ψ(x,y), for the particle is given by 

ψ(x,y)=(2/(Sqrt[(Subscript[l, 1]) (Subscript[l, 2])]))sin((((Subscript[n, 1]) πx)/(Subscript[l, 1])))sin((((Subscript[n, 2]) πy)/(Subscript[l, 2])))   (Subscript[n, 1])= 1, 2, ... and (Subscript[n, 2])= 1, 2, ...

where (Subscript[n, 1])and (Subscript[n, 2]) are the quantum numbers and the energy of each state is given by

E = \(TraditionalForm
      \`\(h\^2\/\(8  m\)\) \((n\_1\%2\/l\_1\%2 + n\_2\%2\/l\_2\%2)\)\)

Ignoring the scaling constants, you can plot the three dimensional surfaces using the following code. We start by looking at the wavefunction for a particle confined to a square surface.

ψ[x_, y_] := 
      Sin[(((Subscript[n, 1]) π x)/(Subscript[l, 1]))] Sin[(((Subscript[n, 2]) π y)/(Subscript[l, 2]))]

Plot3D[
      ψ[x, y]  /. {(Subscript[n, 1]) -> 1, (Subscript[n, 2]) -> 2, (Subscript[l, 1]) -> 1, (Subscript[l, 2]) -> 1}, {x, 0, 
        1}, {y, 0, 1}]

As with two dimensional plots we can change the appearance of the surface plot in a large number of ways. We increase the resolution of the plot, label the axes, change the viewpoint, the length of the edges of the box, within which we plot ψ, and colour code the height of the function. The latter makes high points red and low points orange. The view point was set using the 3D ViewPoint Selector ... under the Input pull down.

\(Plot3D[
      ψ[x, y]\  /. {n\_1 -> 1, n\_2 -> 2, l\_1 -> 1, l\_2 -> 1}, {x, 0, 
        1}, {y, 0, 1}, 
PlotPoints -> 30, ColorFunction -> Hue, 
      AxesLabel -> {"x", "y", "ψ"}, 

      ViewPoint -> {5.375, \ 2.227, \ 2.965}, BoxRatios -> {1, 1, 1}]
      \)

Exercises and Explorations 2

1: Plotting (ψ^2)

The probability of finding the particle at a position (x, y) is equal to (ψ^2). Plot (ψ^2) for the quantum numbers used above. (N.B. you should use the normalisation constant given above.)

2: Contours

Plot the function from exercise 1 as a contour plot, making sure the axes are labelled.

3: Degeneracy

Show that the wave functions for (Subscript[n, 1])= 1, (Subscript[n, 2])= 2 and for (Subscript[n, 1])= 2, (Subscript[n, 2])= 1 are degenerate when (Subscript[l, 1])= (Subscript[l, 2]). That is, you should show that the energy for these two states is the same, but the wave functions are different. What is the relationship between these degenerate wavefunctions?

4: Accidental degeneracy

Obvious symmetries do not always exist between degenerate states. Show that when (Subscript[l, 2]) = 2  (Subscript[l, 1]) there is degeneracy between states (Subscript[n, 1]) = 1, (Subscript[n, 2]) = 4 and  (Subscript[n, 1]) = 2, (Subscript[n, 2]) = 2, but that a simple symmetry operation does not connect the wavefunctions. If you have time, find another non-symmetrical pair of degenerate states for (Subscript[l, 2])= 2(Subscript[l, 1]). N.B. since (Subscript[l, 1])≠ (Subscript[l, 2])  now you will have to change the range over which you are plotting.

Once again, save your work from this part of the lab and quit the kernel, before opening a new Notebook to work on Section 3.

3. Three Dimensional  Surface Representations

The surface plots of section three are single valued functions of x and y, i.e. there is only one value of z for each position (x, y). Frequently we want to make three dimensional surface representations where the surface is a function of x, y and z as well. A good chemical example of this are the molecular orbitals of hydrogen.

In this case the spherical symmetry of the problem indicates that we will simplify the mathematical description if we use the spherical coordinates r, θ and ϕ. Here r is the radial distance from the origin, θ is the angle with respect to the z-axis and ϕ is the angle in the x-y plane with respect to the x-axis. (Hence x = r sinθ cosϕ, y = r sinθ sinϕ and z = r cosθ.) In these coordinates the hydrogenic wave function, Ψ, can be expressed as

Ψ(r,θ,ϕ) = R(r) Θ(θ) Φ(ϕ)

In other words the wave function is the product of separate functions of r, θ and ϕ. The angular parts, Θ(θ) Φ(ϕ), are functions called spherical harmonics and are what give us the hydrogenic molecular orbitals (henceforth MO). The MOs and their associated spherical harmonics are given below.

(OverBar[(UnderBar[({{MO, Spherical Harmonic}})])])
(UnderBar[({{s, Constant}, {(Subscript[p, z]), cosθ}, {(Subscript[p, y]), sinθ  sinϕ}, {(Subscript[p, x]), -}, {(Subscript[d, (z^2)]), 3 (cos^2) θ  -  1}, {(Subscript[d, yz]), cosθ sinθ  sinϕ}, {(Subscript[d, zx]), -}, {(Subscript[d, (x^2) - (y^2)]),  (sin^2) θ cos2ϕ}, {(Subscript[d, xy]), -}})])

We have created the three dimensional renditions of the spherical harmonics by loading  the following package

Needs["Graphics`ParametricPlot3D`"]

which contains code to plot surfaces in three dimensions using spherical coordinates. Below is the code to plot the (Subscript[p, z])orbital. The radial distance from the origin, r, is being plotted as a function of θ and ϕ. In this case according to the MO table r ∝ cosθ, however, because of the way in which SphericalPlot3D does its computations, it is necessary to make r ∝ |cosθ|. This is true for plotting all the spherical harmonics using this function. Opening up the Help file on ParametricPlot3D you may be able to work out why this is.

SphericalPlot3D[ Abs[Cos[θ]], {θ, 0, π}, 
       {ϕ,  0,   2 π}]

The code which follows labels things neatly, and increases the number of points being used to plot the surface to give it a smoother appearance.

\(SphericalPlot3D[\ Abs[Cos[θ]], {θ, 0, π}, 
      \ {ϕ, \ 0, \ \ 2\ π}, PlotRange -> All, 
	
      AxesLabel -> {"x", "y", "z"}, 
      PlotLabel\  -> \*""\!\(p\_z\)-orbital"", PlotPoints -> 25, 
	
      DefaultFont -> {"Helvetica", 8}, 
	
      AxesEdge -> {{\(-1\), \(-1\)}, {1, \(-1\)}, {\(-1\), \(-1\)}}]\)

Exercises and Explorations 3

1: Making molecular orbitals

Replicate the s, p and d orbitals shown above, making sure everything is correctly labelled. You will have to figure out how to produce (Subscript[p, x]), (Subscript[d, zx]), and (Subscript[d, xy]).

2: All the bells and whistles

Choose one of the MOs to be exported to a special Java program where you will be able to turn it around and zoom in and out on it. Render it in colour, using the following code (of course you should replace Abs[Cos[θ]] with the function you wish to plot) :

\(orbitalPlot1 = 
      SphericalPlot3D[\ Abs[Cos[θ]], {θ, 0, π}, 
        \ {ϕ, \ 0, \ \ 2\ π}, Boxed -> False, 
		Axes -> False, 
        PlotPoints -> 25, \ AmbientLight -> RGBColor[0.6, \ 0, \ 0.1]]\)

The reason for the "AmbientLight" here is prevent the MO having an unilluminated side when you load it into the Java program. (More sophisticated lighting and colouring effects are possible, if you want to venture into advanced Mathematica territory you might start by looking up AmbientLight in the help system.)

We suggest that you now save the work you have so far done for section 3 - deleting the graphical outputs first.

Now, locate the notebook "LiveGraphics3D.nb" on your hard disk (it's also on the Maths Lab web site) and open it. Follow the instructions there to set up the WriteLiveForm function. Then in your own notebook type:

WriteLiveForm["orbital1.html", orbitalPlot1]

This creates a file called "orbital1.html" containing a Java version of your orbitalPlot graphic (the filename is arbitrary, but should always end in ".html" or ".htm"). You should note that as written, the file will be placed in the "working directory". The command Directory[] will show what this isand you can use SetDirectory to change it.

 To view your orbitalPlot graphic start up a web browser. In Netscape you will now need to select Open Page from the File menu, whilst in Internet Explorer you will need to select Open File. Now locate the file you just created. All being well, your MO graphic should soon appear: just drag the mouse to rotate, and with the shift key pressed, drag down or up to zoom in or out. 

(Note for people with a personal computer: by default, these html files call up some special Java files from METRIC's web site. For instructions on how to use LiveGraphics3D without needing a network connection, see the "LiveGraphics3D" notebook.)

We'll be doing more with LiveGraphics3D in the next session. Note that, with a few exceptions, you can use it to look at any 3D Mathematica graphic.