Fortran IO

There are many ways to do IO in Fortran. We begin from the simplest possible.

The easiest way to read an array from standard input in Fortran is simply to say:

read(*,*) x

where x is an array. If the array is of type real and is 10 entries long, Fortran  will attempt to read 10 floating point numbers. The read command reads a single line of standard input at a time. This means that all array entries have to be typed on a single line.

Figure 2.5 shows a content of an input file, let's call it chi.dat, to be read by our program, in which we must, at this stage, replace

integer, parameter :: n = 100

with

integer, parameter :: n = 10

You will learn soon how to allocate an array of desired length in Fortran.

  

Figure: The content of file chi.dat: the first row corresponds to x_i, the second row corresponds to y_i and the last row corresponds to $\sigma _i$.

In the program itself we will place the following three lines in front of the computational part:

read(*,*) x
read(*,*) y
read(*,*) sigma

and place the following single line, which writes  the result of our computation on standard output, just before the word end:

write(*,*) a, sigma_a, b, sigma_b, chi_2

The program, after recompilation, can now be run as follows:

gustav@blanc:../src 12:57:44 !522 $ ./chi < chi.dat
 1.8082818294190344,   8.7731094559548306E-2,  \
 0.51538936959208925,  1.43532455281008345E-2,  \
15.310414091470935
gustav@blanc:../src 12:57:52 !523 $

The output will be actually written on a single line, because, as read reads a line at a time, write writes a line at a time.

The results returned by the program are as follows:

$$\begin{eqnarray*}a &=& 1.808 \pm 0.088 \\ b &=& 0.515 \pm 0.015 \\ \chi^2 &=& 15.310 \end{eqnarray*}$$

I have rounded the results so that we don't have to look at more digits than makes sense. The value of $\sigma_a$ and $\sigma_b$ determine the number of digits that need to be taken into account. But remember the rule: when you round an error, i.e., $\sigma_a$ or $\sigma_b$, you should always round it up, even if the usual rounding rules would let us round 0.0143 down to 0.014. With errors it is always prudent to assume less accuracy than more.

Before we go any further, let us restructure our program so that we can view the results of our $\chi ^2$ fit on a graph. We do that by printing on standard output our results in multiple columns. The first column contains x_i, the second column contains y_i, and the third column contains $\sigma _i$. In the fourth column we will print the values of

y'_i = a + b x_i

Here is the part of the code that does the job:

real(kind=long), dimension(n) :: y_prime
integer :: i
...
y_prime = a + b * x
do i=1,n
   write(*, '(4f8.3)') x(i), y(i), sigma(i), y_prime(i)
end do

Observe the new Fortran  element, the DO loop, which is Fortran's main iterative construct. We have used this construct to produce pretty output, but at no stage did we have to use DO during computations. The string '(4f8.3)' in the write statement, means ``print 4 floating point numbers, reserving 8 spaces for each number and with 3 digits after decimal point''. The printed numbers are right justified within the space allocated to them.

Running this program as before returns a table of numbers that is shown in Figure 2.6.

  

Figure: The results of $\chi ^2$ fitting printed in four columns, which correspond to arrays x_i, y_i, $\sigma _i$, and y'_i = a + b x_i.

Let us save this output, say, on a file chi.out. We can now use gnuplot  to view the results.

  

Figure 2.7: Ploting the output of program chi with gnuplot.

If you proceed as shown in Figure 2.7 you will not only get a plot on your X11 display. You will also get an encapsulated postscript file chi.eps, which, when included into a TeX document looks as shown in Figure 2.8.

  

Figure: $\chi ^2$ fitting of data from Figure 2.5.

Looking at the plot you can appreciate that points with smaller error bars have markedly greater weights than points with larger error bars. Consequently the $\chi ^2$ fit tries to come closer to the former.

The syntax of the gnuplot command plot is quite complicated, but you don't have to use it in its most complex form. Suffice to remember that the keyword using tells plot, which columns of the input file to use for the abscissa, the ordinate, and for error bars.

Let us sum up our accomplishments, because even at this early stage they are not inconsiderable. We have a 28 lines long Fortran program, shown in Figure 2.9, and a 6 lines long Gnuplot script shown in Figure 2.7, and with these two very simple tools we have done in just a few seconds what used to take me good few hours when I was a student myself and had to do all these computations by hand, because all that we had at that time were logarithmic rulers, whereas computers were very scarce, expensive and reserved for top science and military work only.

  

Figure: The $\chi ^2$ fitting program with simple IO