Exercise 2: Hydrogen Ionization

Hi,

This exercise contains ready-to-run IDL, Fortran-90, Maple and Mathematica programs, and is a jump start to programming in these four environments.

As discussed in the lecture notes, there are several purposes with this exercise:

1. To provide a simple example of the ionization of atoms. Hands-on experience with the solution of physical problems is the fastest way to get a basic feeling for "how it works".
   
2. To illustrate that even a seemingly straightforward evaluation often contains "numerical traps", such as
   • the possibility that intermediate results in a calculation become too small (or too large) for the standard (IEEE) representation of floating point values in computers, or
   • loss of precision, because two numbers that are almost equal are subtracted from each other – the internal precision is limited (roughly 7 digits) and thus the difference between two numbers has correspondingly less precision.
     
3. To introduce you to the three (four) language environments that we will be using in the rest of the course:
   • IDL (Interactive Data Language), an efficient environment for working interactively with data, using line graphics, two-dimensional and three-dimensional visualization, statistical analysis, or just simple computations.
   • Fortran-90 (and its siblings Fortran95/2003 and High Performance Fortran); a simple and efficient language for computationally intensive applications.
   • Maple (and Mathematica), a system in which one can do both symbolic computations, arbitrary precision numerical computations, graphics and visualizations.

   Each of these three tools is useful to have in your mental "toolbox", and the overall purpose of all of these exercises is basically to help you learn how to use them efficiently.

We will be using CVS (Concurrent Versions System) to distribute the exercise files. Later, we will look at CVS in more detail, but for now it is only a "black box" method to get the files you need for the exercise.

To get the files, do the following

  cd ~/ComputerPhysics             (or whereever you have it)
  cvs update -d

NOTE: If you get a message such as "permission denied" you need to do "cvs login" and then repeat the "cvs update -d" command.

This updates the directory with the files needed for this exercise, which end up in the directory ComputerPhysics/2_H+He.

We will first use idlde, the 'Development Environment' for IDL (Interactive Data Language), to look at why there is a problem with computing one of the roots of a second degree polynomial with the standard high-school formula. Some of you all ready know this problem from Dat-f, it should therefore be relative easy to answer the questions and you should mainly spend the allocated time to learn about the new working environment.

• Start idlde from an terminal window command line, by doing

    cd ~/ComputerPhysics/2_H+He  # Try hitting the TAB key instead of typing the whole thing
    idlde &                      # The "&" starts the command "in the background"
  

  NOTE: If this results in an error message, just use "idl" instead of "idlde". No fancy windows, but you don't need them this time.

  Shortly after idlde asking where to locate a workspace. This is a directory where it stores information about the state of idlde. This information is used to restore the working environment you used last time running idlde. Choose to make a new directory (e.g. called workspace) inside the present exercise directory and use this one.

  A window with several panels, buttons, and menus appears. Take a few minutes, look around, and familiarize yourself with the various elements of the window, but don't be alarmed if you initially only have a vague idea about what many of them might be good for. Later on, you should spend time reading introductions and documentation, but for now we are just going to use it.

• At the bottom of the window you find the one-line command input pane, where one can type IDL commands. The outputs from the commands appear in the output log pane. There is also a variable trace pane--ignore it for now (you may temporarily close it from the "Windows" menu).
• We will look at the positive root of a second order polynomial with a=0.5, c=-0.5, for various values of b. Here is what you should type in the command line window pane (for clarity, the IDL> prompt is shown in a different color):

          IDL> a=0.5
          IDL> c=-0.5
          IDL> b=1.
          IDL> print, (-b+sqrt(b^2-4.*a*c))/(2.*a)
  

  • Try the up-arrow key on the key board, with the IDL command window panel (the lower most one) active, and note that this recalls old command lines. So, you seldom have to type an IDL command more than once. Also, if you make a typing mistake, you can fix it without having to retype the whole line
  • The usual mouse driven editing works, including copy / paste from the output log panel to the input line. Keyboard short-cuts may be found by looking in the menus. If you are a keyboard short-cut lover, you may find additional keyboard short-cuts for input line editing in the on-line documentation (menu "Help" -> Contents -> Using IDL -> Running IDL -> Input to IDL).
• It is possible, and often convenient, to type several commands on one line, separated by an ampersand (&). For example, we may set a value of b and print the root on a single line. Using keyboard short cuts, you can do this by hitting first the up-arrow, then Ctrl-a, then typing "b=1e1 & ". The resulting line is

          IDL> b=1e1 & print, (-b+sqrt(b^2-4.*a*c))/(2.*a)
  

• Go through increasing powers of ten for the value of b (use "up-arrow, HOME, Ctrl-F" to move the cursor in place for changing the value of b). The resulting lines are

          IDL> b=1e2 & print, (-b+sqrt(b^2-4.*a*c))/(2.*a)
          IDL> b=1e3 & print, (-b+sqrt(b^2-4.*a*c))/(2.*a)
          IDL> b=1e4 & print, (-b+sqrt(b^2-4.*a*c))/(2.*a)
          IDL> b=1e5 & print, (-b+sqrt(b^2-4.*a*c))/(2.*a)
          IDL> b=1e6 & print, (-b+sqrt(b^2-4.*a*c))/(2.*a)
          IDL> b=1e7 & print, (-b+sqrt(b^2-4.*a*c))/(2.*a)
  

• There is clearly a problem with this expression when b becomes large.

  How large is b (even powers of ten), when the result first vanishes? 10^4 10^5 10^6 10^7 Credits: 4/-4 OK

• An alternative expression for the root in question may be obtained by observing that the product of the two roots is equal to c/a (cf. the discussion in Chapter 5 of Numerical Recipes). Append the correct expression to the line, by typing "up-arrow, Ctrl-E, 2.*c/(...)".

          IDL> b=1e5 & print, (-b+sqrt(b^2-4.*a*c))/(2.*a), 2.*c/(...)
  

• Exactly what to replace the (...) with should be clear from either
  • The lecture notes, or ...
  • ... the discussion in Numerical Recipes. What? You haven't got it yet? Well then you can either look it up on-the-fly, using the free web access to Numerical Recipes, or you're left with the final alternative ...
  • for this to work you may have to create a directory called .fileopen in your home directory.
  • If you want it to work on your own laptop/desktop, you need to use Acrobat Reader version 7 or 8 and have a plug_in installed. You can find details
  here
  • ... derive the formula from what was just mentioned about the relation of the two roots to the coefficients a and c.
• If you have time now (if not, come back to this point some other time), you may want to check up on a couple of extra questions:
  • What if one just computes in double precision (free hint)? You may try this by simply using double precision for b:

            IDL> b=1d5 & print, (-b+sqrt(b^2-4.*a*c))/(2.*a)
    

    How large is b (even powers of ten), when the result first vanishes?? 10^12 10^18 10^8 10^6 Credits: 4/-4 OK

  • How large is b when one first gets into trouble (zero value) even with the improved expression (try with b again in single precision)

            IDL> b=1e17 & print, 2.*c/(...)
    

    How large is b in single precision for the improved formula? 10^17 10^20 10^38 Credits: 4/-4 OK

    What could one do about it, in addition to making b double precision? Try to figure it out, before taking this free hint? (No minus points - just your selfsatisfaction if you can figure it out without the hint!-:).

Next, we will look at a ready-made IDL program, where the problem with computing one of the roots of a second order polynomial is acute.

• Open the IDL program H+He.pro, in the ComputerPhysics/2_H+He/ directory. Use the Open .. menu item in the File .. menu of your idlde, and double click, first on the directory names, then on the file name.
• The IDL program appears in the upper window panel. Read through it to familiarize yourself with the program. Use the scrollbar at the right window edge to reach the bottom of the file.
  • If you prefer to see the program in a separate window, use the menus File -> Preferences -> Layout, and press the Editor layout: multiple button to change your preferences permanently. For a quick change, there is a Multiple windows menu item in the Window menu.
• Run the program, by clicking first on the Compile button (gear weels - first button in the second group), then on the Start or resume ... button (green arrow).
  • Alternatively, use the menu button in the pane showing the program, and chose Compile and then Run.
  • Keyboard short-cut fans may look in the main idlde window Run menu, to discover that Ctrl-F8 followed by F8 is what they need.
• A surface plot shows the result; the fraction of electrons out of all particles, as a function of temperature and logarithmic pressure.
  
  • At sufficiently high temperatures and low pressures, hydrogen is fully ionized, so there are 10 electrons for every 21 particle (ten electrons, ten protons, one helium atom); i.e., the value should be slightly less than one half. It is always a good idea to do such "sanity checks", using "back of the envelope" calculations.
• Note the problem with precision that occurs in the corner with low densities and high temperatures. The reason is exactly the same as above.
• Find the line where the problem originates. Turn the original line into a comment line, using ;, and add a comment about why this calculation had a problem.
• Eliminate the problem, by adding a new line that calculates the answer in the way you figured out above.
  • Use the fact that you know that b is positive and c is negative, and choose the appropriate signs in the formula (no need for a SIGN function)
• Verify that the result is OK after the change. For this, you need to save the file, using the Save button (diskette pictogram), and then press the Compile and GO buttons again.
• The program uses the IDL surface procedure to make the surface plot. Change the call, so that xsurface is used instead (note that xsurface only takes one argument, so remove the rest of the arguments).
  • Take the opportunity, and learn how to find the on-line documentation for a particular built-in procedure. Use "Help" -> Help Contents. A new window opens which contains all the documentation for IDL. Type xsurf into the top search line and select highlighted xsurface from the list of hits.
• Use the window buttons to rotate the view, and to change between a 3-D surface view and a 3-D wire frame view.
• Play around with the buttons and menus (trying, for example, the various color tables and palettes)

Questions:

1. The problem with numerical precision occurs in an expression equivalent to (in IDL syntax)

   f = -1 + sqrt(1 + epsilon)
   

   For small values of epsilon, f is approximately equal to epsilon epsilon/2 sqrt(epsilon) Credits: 4/-4 OK

2.  

   If the square root function has 7 significant digits, and epsilon is equal to 2 10^-4 (and also has 7 significant digits), how many significant digits are there in f? (compare the error in the result with the size of the result) 3 3.5 4 7 Credits: 4/-4 OK

The 2_H+He/H+He.f90 program solves the same problem, but since Fortran is a language that requires compilation (i.e., turning the program text into a binary, executable file), there are more steps to go through. Also, in order to plot the results, data must be written to a file and then read with IDL.

At this point you need to edit a Fortran text file. Most of you probably already have a favorite text editor that you can do it with. If not, you may want to try kedit, nedit, gedit, kate or possibly emacs / xemacs -- you may then find Nina Jansens 45 minute emacs tutorial useful - or vim -- here is an introduction by Bo Milvang.

• Open the file H+He.f90 in the editor, and take a look at it, comparing it with the IDL procedure (in the IDLDE window, or open it too in another window of the editor). Notice that, while the IDL program makes a plot of the results at the end, the Fortran program writes the data out into a file.
  
  
• To compile and run the Fortran program, do the following in a terminal window:

  cd
  cd Co<TAB>/2<TAB>             (where <TAB> represents the TAB-key; see tab expansion)
  ifort H+He.f90                (or use gfortran, g95, or whatever Fortran compiler you have ...)
  ./a.out
  

• To read the result from the Fortran-90 program into IDL and make a surface plot:
• • open the IDL procedure read_surf.pro in IDLDE
  • compile it (press the compile icon) and
  • run it (press the run icon
    
  • This gives a problem...
    • Where ever you are starting idlde from, it always assume you are standing in your home directory when referring to file names in programs.
    • Either change the name of the file in read_surf.pro to the full path name
    • Or get idl to change its location using this command:

      IDL> cd,"ComputerPhysics/2_H+He"
      

    
    
• Fix the problem in the Fortran-90 program, in the same way as in the IDL program.
• • Recompile (ifort command) and rerun (a.out command) the Fortran-90 program.
  • Verify that the result is now OK, by looking at the output numbers, and plotting the results again with IDL.
  • If there is a compilation problem with the new Fortran-90 line (or the comment line), or the result is not correct, fix it and try again until it works.

There are two competing softwares for doing analytical mathematics, namely Maple and Mathematica. Both are available at NBI and are used by different research groups. Here you have the possibility for trying out both packages. From the CVS download files for both packages are included. It will therefore be simple for you to compare the two competing softwares. To start Maple or Mathematica from an xterm prompt, do

    cd ~/ComputerPhysics/2_H+He  # Use the TAB key instead of typing the whole thing
    xmaple & # or mathematica &

Use the File -> Open menu to find and open the Quadratic.{mw/nb} and H+He.{mw/nb} notebooks in the 2_H+He subdirectory.

• Execute the notebook cells, one by one: Put the cursor in the first cell, and type SHIFT-ENTER. Repeatedly typing SHIFT-ENTER executes subsequent cells (the first cell takes a while, since Mathematica must first start the so called "kernel"; a background program that does all the computations).
• Note that Mathematica works with much higher precision than the hardware floating point numbers that Fortran and IDL use (this is one of the reasons why Mathematica is much slower), and that, therefore, the problem with round off errors is only reproducible at the same pressures by explicitly adding "synthetic round-off errors" (Random noise). At lower pressures, the problem is visible without adding extra noise.
• Feel free to play around with the notebooks, changing or adding new cells. If you want to save a copy, use "Save as .." to save it under a different name.

Mathematica documentation (including the full Mathematica Book) is available on-line, from the HELP menu in the upper right hand corner of the Mathematica window.

• Find the "Tour of Mathematica" help entry, and try out some of the examples (the cells are editable, and executable!).

• To save space (your account has a finite quota of disk space), remove the a.out and *.res files, and other large files that you don't need, after you are done (throw them in the file manager waste basket and empty the waste basket, or else use the Unix/Linux rm command).

If you still have time left at his point; use it to play with IDL, Maple or Mathematica (if you don't have time left, do this as part of the home work). There are a number of IDL demos, that may be accessed by typing demo at the IDL prompt.

After you are finished with this exercise, use time at home or Wednesday to study and compare the IDL, Fortran-90, Maple and Mathematica programs. A good exercise might be to trace the use of one-dimensional and two-dimensional variables in the programs, and mark them by single and double underlining on your hard copy. Use the variable trace pane to browse the variables in IDL. Use analogy reasoning to figure out the corresponding properties of the Fortran-90, Maple and Mathematica variables.

Take ample time to review relevant sections (on assignments, declarations, Parameters, etc.) in one of the Fortran book. Familiarize yourself with IDL, by accessing the on-line PDF copy of Getting Started in IDL. If / when you have installed IDL, you should also familiarize yourself with where to find the IDL help and documentation files on your own PC.

Note that you may use CVS from home, either from a terminal window under Linux, or from the WinCVS and MacCVS graphical user interface clients, to check out and maintain a consistent set of exercise files on your home computer or laptop. There will be more information about this later, but if you are eager to try this out you may want to start out with the CVS tutorial. With access to a cvs command line (e.g. from Linux), the CVS commands used in the exercises should work from home as well.