Chapter 11 Appendix Introduction to FORTRAN by rrk61112


									Chapter 11

Appendix: Introduction to FORTRAN

Fortran is a relatively old programming language, going back to the 1950’s, and has undergone
many changes over the years. Currently there are three versions of Fortran in common use:
Fortran 77, Fortran 90 and Fortran 95. There are few differences between Fortran 90 and Fortran
95 (none which occur in this description), however, Fortran 77 is quite different. In Fortran 77 the
first several spaces in each line are reserved for special characters, for instance, lines numbers,
line continuations, comments, etc. In Fortran 90 these no longer exist and program statements
can start in the first space on a line. Additionally, a new program component called a “Module”
was introduced which will be discussed later. This summary is meant to describe the main
features in Fortran 90.
    It is also important to note that all of the versions of Fortran are not case sensitive. Hence, in
Fortran print = Print = PRINT, and variable = Variable = VARIABLE. In this summary for clarity
Fortran intrinsic functions (such as PRINT, CALL, PROGRAM, INTEGER) will be written in
upper case letters.

11.1    Sample Program 1: Solving the Quadratic Equation
Before we begin to examine the components of a program, let’s look at an example.
Step 1: Specify the Problem
   We would like to write a program which evaluates the roots of the quadratic equation ax2 +
bx + c = 0 which are given by the following formula:
                                             −b ±     b2 − 4ac
                                        x=                                                      (11.1)

Note that there will be three cases of interest, when the discriminant (D = b2 − 4ac) is 0, positive
or negative.
Step 2: Write a Pseudocode (or Algorithm)

   • read in the values of a, b, and c for which you would like to know the roots.

   • if a is zero stop - in this case we do not have a quadratic equation

   • evaluate D
   • if D = 0 then the root is   2a
                                              −b± D
   • if D > 0 then there are two real roots     2a

98                                 CHAPTER 11. APPENDIX: INTRODUCTION TO FORTRAN
                                                     −b±i −D
     • if D < 0 then there are two imaginary roots      2a

     • print the solution
Step 3: Write the Program

! PROGRAM: QuadraticEquationRoots ! DESCRIPTION: This program
calculates the roots of the Quadratic Equation ! DATE: Aug 25,
REAL    :: a, b, c, D REAL    :: realcomp, imagcomp, root1, root2

! EXECUTION PART PRINT*, ’What are a, b, and c?’ READ*, a, b, c
PRINT*, a, b, c

IF (ABS(a).GT.0.e-10) THEN
  D = b**2 - 4*a*c
  IF (ABS(D).LT.0.1e-10) THEN
       root1 = -b/(2.0*a)
       PRINT*, ’There is only one root and it is’, root1
  ELSE IF (D.GT.1e-10) THEN
       root1 = (-b + sqrt(D))/(2.0*a)
       root2 = (-b - sqrt(D))/(2.0*a)
       PRINT*, ’There are 2 roots which are’, root1, ’and’, root2
       realcomp = -b/(2.0*a)
       imagcomp = sqrt(-D)/(2.0*a)
       PRINT*, ’There are 2 roots which are’, realcomp, ’+ i’, imagcomp, &
              ’and’, realcomp, ’- i’, imagcomp
  PRINT*, ’PROBLEM: this is not a quadratic equation since a is 0’

END PROGRAM QuadraticEquationRoots

11.2     Program Components
A program usually consists of two parts:
     • an Initialization part where all the variables are defined

     • a Execution part which is the main body of the code where data is read in, manipulated,
Note that as shown in the example, programs in Fortran always start with the statement
PROGRAM program_name
and end with
11.3. INITIALIZATION PART                                                                        99

END PROGRAM program_name
Additionally, comments are inserted into Fortran 90 by proceeding the comment by a ‘!’. For
! There are three possible cases, D =,>,< 0.
D = b**2 - 4*a*c ! D = Descriminant

11.3     Initialization Part
As seen in the example, before the main execution part of the code starts, there are a number of
initializations, such as variable declarations. Here is a short list of some of the statements found
in the initialization part:
   • IMPLICIT NONE (this should always be included as it forces all variables to be defined. )

   • INTEGER :: var1 (defines “var1” as an integer such as -1, 10, 12838)

   • INTEGER, PARAMETER :: N = 10 (defines N to be a parameter ALWAYS equal to 10, you
     cannot change the value of N later on)

   • REAL :: var2 (defines “var2” as a real number such as 2.18847)

   • REAL, PARAMETER :: Pi = 3.1414 (defines ‘Pi’ to be a real number equal to 3.1415)

   • CHARACTER :: initial ( defines ”initial” as a single alphanumeric character such as “a” or
     “Z”, etc.)

   • CHARACTER(LEN = 12) :: name ( defines “name” to be a string containing 12 alphanu-
     meric characters such as “Program test”.)

   • LOGICAL :: raining ( defines “raining” to be a logical variable which takes on the possible
     values of .TRUE. and .FALSE.)
Note that PARAMETER statements can also be used in conjunction with CHARACTER and
LOGICAL statements.
    There are also both single and double precision variables in Fortran 90. All of the above
examples are in single precision. To use double precision we write REAL(8). Double precision
in a program is written with a D as in 3.411D0.

11.3.1   Implicit Typing
As mentioned previously, the best way to avoid problems in your code is to include an IM-
PLICIT NONE statement at the beginning of the program. If you do not do this, then all vari-
ables which you have not declared in the initialization part of the code have an implicitly de-
clared type which depends on the first letter in the variable name. If the first letter of a name is
i, j, k, l, m, or n then the variable is automatically of type INTEGER, otherwise it is assumed to
be REAL. For instance, if you had wanted to have
... REAL         :: epsilon epsilon = 1.0 ... PRINT*, epsilon
but instead had written
100                                CHAPTER 11. APPENDIX: INTRODUCTION TO FORTRAN

... REAL         :: epsilon espilon = 1.0 ... PRINT*, epsilon

then if you had NOT included IMPLICIT NONE the program would have printed junk instead
of 1.0 when you asked it to print. There are many ways in which this can cause problems in
your code, so the rule of thumb is to always include an IMPLICIT NONE statement.

11.4     Execution Part
11.4.1   Input and Output
To and From the screen
To set a prompt for the user, for instance as in the example program for Quadratic Roots, we use

READ*,     a, b, c

During the execution of the program the user will be prompted with the request “what are a,
b, and c” and then can enter these values into the terminal. Note that the * indicates that we
are using “standard input” for the file formating. Formatting will be discussed in a following
    To write to the screen we write:

PRINT*, ’the roots are’, root1, ’and’, root2

Note that everything in quotations is printed as is, and commas are used to separate elements
in the PRINT statement. The * in this case implies “standard output”.

To and from files
Before reading or writing to a file, the file must be first opened. To open a file called ‘file-
name.dat’ we use the following statement:


The first number (in the above example 21) gives us a way to identify the file later and can be any
positive integer except 5, 6, 100, 101, and 102 which have default meanings. It is also possible
to give information about HOW a file is to be opened, for instance, it is possible to append data
to the end of a file without writing over the file itself. To do this you can include a POSITION
comment in the OPEN statement as follows:


There are many other options available when opening files which are not described here. See
any text on Fortran 90 for a complete list of the options.
   After opening a file you can read data from the file. If “filename.dat’ contains a list of N
pieces of data to be read into the array variable1 then we use the statement

DO i=1,N
   READ(21,*) variable1(i)

    It is also possible to write to the file ‘filename.dat’. To write the contents of the array variable1
to this file use the statement:
11.4. EXECUTION PART                                                                            101

DO i=1,N
  WRITE(21,*) variable1(i)
   When finished with a file it can be closed with the statement

While in the previous sections the formatting of the file was done automatically, it is also possi-
ble to specify the formatting when reading from or writing to a file. For instance, it is possible to
control the number of decimal places included when writing a REAL number. This is possible
through the use of the FORMAT statement. As an example,
15 FORMAT (2E10.4,3I5)
implies that there are 2 exponential numbers each taking up 10 character positions with 4 num-
bers after the decimal followed by 3 integers each taking 5 of the next character positions. The
number ‘15’ preceding the word FORMAT is a label for the statement and can be any integer. It
will replace the * in the PRINT, READ, and WRITE statements to indicate which format state-
ment you wish to use. The following table lists the options in format statements:
 Format Interpretation
 Is         Outputs an integer into the next s character places.
 Fs.d       Outputs a real number into the next s character spaces
            with d decimal places.
 Es.d       Outputs a real number in exponential form in the next s
            character spaces and with d decimal places.
 As         Outputs a character string into the next s character

    There are several ways to use the FORMAT statement, however, only a few will be included
here. The following consists of a list of some of the more common applications (in all examples
it has been assumed that the following FORMAT statement
40 FORMAT(E10.4,A7)
has also been included in the code) :
   • PRINT 40, real1, character1
   • WRITE(21,FMT=40) real1, character1
   • READ(21,FMT=40) real1, character1

11.4.2   Numeric Operators
Fortran contains the following numeric operators:
 Operator Function             Example
 +          Addition           2+3
 −          Subtraction        922.18 − 9263.44
 *          Multiplication     3×6→3∗6
 /          Division           2.2 → 4.9/2.2
 **         Exponentiation 24 → 2 ∗ ∗4
102                               CHAPTER 11. APPENDIX: INTRODUCTION TO FORTRAN

11.4.3   Relational Operators
The following relational operators are   included in Fortran (for instance, for use in an IF state-
 Operator Function                        Integer or Real
 .GT.       Greater Than                  Integer and Real
 .GE.       Greater Than or Equal To      Integer and Real
 .LT.       Less Than                     Integer and Real
 .LE.       Less Than or Equal To         Integer and Real
 .EQ.       Equal To                      Integer
 .NE.       Not Equal To                  Integer

11.4.4   Logical Operators
Logical operators produce either a .TRUE. or .FALSE. result. The following table lists the logical
operators in Fortran (assume α and β are logical variables):
 Operator Usage           Explanation
 .AND.       α.AND.β      .TRUE. when both α and β are .TRUE.
                          .FALSE. otherwise
 .OR.        α.OR.β       .TRUE. when either α or β are .TRUE.
                          .F ALU E. when both are .FALSE.
 .NOT.       .NOT. α      .TRUE. when α is .FALSE.
                          .FALSE. when α is .TRUE.
 .EQV.       α.EQV.β      .TRUE. when α and β are the same (ie. both
                          are .TRUE. or .FALSE.)
                          .FALSE. otherwise
 .NEQV.      α.NEQV.β .TRUE. when α and β are different (ie. one
                          is .TRUE. and the other is .FALSE.)
                          .FALSE. otherwise

11.4.5   IF Statements
There are essentially two types of IF statements: a single line IF and a block-IF statement. An
example of a single line IF is
IF (i.LT.10) PRINT*, ’i is less than 10’
An example of a block IF statement is:
IF (i.LT.10) THEN
     PRINT*, ’i is less than 10’
ELSE IF (i.GE.10.AND.i.LT.20) THEN
     PRINT*, ’i is between 10 and 19’
   PRINT*, ’i is greater than or equal to 20’

11.4.6   DO Statements
There are two types of DO loops, ones which execute the loop a fixed number of times and ones
which are terminated on reaching some condition. An example of a “fixed” loop is given by the
11.4. EXECUTION PART                                                                            103

DO i=1,10
    PRINT*, ’i is’, i

An example of a loop terminated based on a condition is

counter = 0 DO
  counter = counter + 1
  PRINT*, ’counter is’, counter
  IF (counter.GE.10) EXIT

Both programs will print the numbers from 1 to 10, and then exit the loop. However, depending
on what you are doing, you may find one form more advantageous than the other.

11.4.7   Arrays

There are many times where it is useful to be able to store lists of data. To do this we use arrays.
Here are many examples of defining an array

INTEGER                 :: A(1:3,1:2) REAL                                              ::
B(-10:10,0:6,2:3) INTEGER                  :: C(3,2) INTEGER,
DIMENSION(3,2)         :: D REAL, DIMENSION(-10:10,0:6,2:3)                                ::

The arrays A,C, and D are 3 by 2 integer arrays, starting with the index (1,1). For instance
A would contain the elements A(1,1), A(1,2), A(2,1), A(2,2), A(3,1), A(3,2) which are all of type
INTEGER. The arrays B and E contain REAL numbers and are of the form B(-10,0,2), B(-10,0,3)....
There are three ways to reference an array: you can reference a specific element, such as A(1,1),
or the entire array such as A, or even an array section. There are several different ways to
reference an array section, however, only one will be discussed here. You can access an entire
section of an array using a ‘:’, for example, if we had an array Particles(1:10,1:3) where the
first indice referred to a particle number and the last indice referred to the particle coordinate,
we could then access the coordinates of particle 3 by writing Particles(3,:). The result is a 1
by 3 vector. This shorthand notation is often useful. For instance, if we wanted to find the
displacement vector between particle 3 and particle 8 then we could write

REAL        :: displacementvect(1:3) REAL        ::
Particles(1:10,1:3) ... displacementvect = Particles(8,:) -

Whole Array Operations

Here are some operations which can be carried out on whole arrays (and parts of arrays).
104                               CHAPTER 11. APPENDIX: INTRODUCTION TO FORTRAN

 Operation               Requirements         Explanation
 A = 0.0                 None                 sets all elements of A to 0.0
 A=B                     A,B conformable      assigns the values of B to
 A=B+C                   A,B,C         con-   preforms a “matrix” like
                         formable             addition of B and C and
                                              assigns it to A.
 A=B−C                   A,B,C         con-   preforms a “matrix” like
                         formable             subtraction of B and C and
                                              assigns it to A.
 A=B∗C                   A,B,C         con-   NOT MATRIX MULTIPLI-
                         formable             CATION, instead multi-
                                              ples the elements on a ele-
                                              ment by element basis like
                                              addition and assigns it to
 A = B/C                 A,B,C         con-   divides the elements of B
                         formable             with C on an element by
                                              element basis and assigns
                                              them to A.
 A∗α                     α scalar, A array    multiples all elements of A
                                              by α
 A/α                     α scalar, A array    divides all elements of A
                                              by α
 A+α                     α scalar, A array    adds all elements of A to α
 A−α                     α scalar, A array    subtracts all elements of A
                                              by α
 MATMUL(A,B)      N columns in A =            matrix multiples A and B
                  N rows in B
 DOT PRODUCT(A,B) A and B are 1xN             computes the dot product
                                              of A and B

Note that ‘conformable’ means that the arrays have the same size, for instance A(1:N,1:M) and

Allocatable Length Arrays
It is also possible in fortran to use arrays with a variable length, ie. arrays whose size is deter-
mined by the program while it is running. To declare such arrays the syntax for the initialization
REAL, ALLOCATABLE           :: ravg(:,:)
In the main code we could then use
N = 10 M = 5 ALLOCATE(ravg(1:N,1:M))
and then ravg is a N by M array. At the end of the program, or when the program no longer
needs this array it can be cleared by:
11.4. EXECUTION PART                                                                         105

11.4.8   Some Fortran Implicit Functions
There are many implicit functions in Fortran which manipulate integers, reals and character
variables. Some more useful ones will be described here, however, see

for a more complete list and explanations.

Type Operators
 Function       Input    Output   Description
                Type     Type
 FLOAT(i)       I        R        outputs a real approximation to i
 INT(z)         R or I   I        outputs the truncation of x
 DFLOAT(d)      I        RDP      converts d to double precision real
 CEILING(x)     R        I        outputs the smallest integer greater
                                  than or equal to the x
 FLOOR(x)       R        I        output the largest integer less than or
                                  equal to x
 MOD(a,b)       I or R   I or R   the remainder on division of a by b

Note that R, I and RDP refer to real, integer, and real double precision respectively.

Mathematical Functions
 Function       Explanation
 ABS(x)         absolute value of x
 SQRT(x)        square root of x
 EXP(x)         exponential of x
 LOG(x)         natural logarithm of x
 COS(x)         cosine of x, assumes x is in radians
 SIN(x)         sine of x, assumes x is in radians
 TAN(x)         tangent of x, assumes x is in radians
 ACOS(x)        arccosine of x, result is in radians
 ASIN(x)        arcsine of x, result is in radians
 ATAN(x)        arctangent of x, result is in radians
 ATAN2(y,x)     arctangent of y/x where the signs of signs of x and y
                are used to determine the quadient for the angle (See
                note following table)
 COSH(x)        hyperbolic cosine of x, assumes x is in radians
 SINH(x)        hyperbolic sine of x, assumes x is in radians
 TANH(x)        hyperbolic tangent of x, assumes x is in radians

Note that all of these functions can be converted to double precision by adding a D to the begin-
ning, for instance, SIN(x) → DSIN(x).
    The function ATAN2(y,x) is often useful when determining angles. Here is a better descrip-
tion of how it works taken from the link mentioned above.
106                                 CHAPTER 11. APPENDIX: INTRODUCTION TO FORTRAN

      • The result is in radians and is the principal value of the argument of the complex number
        (X, Y).

      • The result approximates arctan(Y/X) for all X.

      • The range of the result is −π < AT AN 2(Y, X) ≤ π.

      • When Y > 0 the result is positive and when Y < 0 the result is negative.

      • When Y = 0 then if X > 0 the result is zero and if X < 0 the result is π.

      • When X = 0 the absolute value of the result is π/2.

11.4.9     Functions
An example function is

FUNCTION average(N,array) IMPLICIT NONE INTEGER                           :: N REAL
:: array(1:N) INTEGER     :: i REAL        :: sum

sum = 0.0 DO i=1,N
   sum = sum + array(i)

average = sum / REAL(N)


    In the main program it is necessary to declare the function. This is done by including the
following statement in the initialization part of the code:

REAL, EXTERNAL                :: average

Now, if in the main program there is a real array A of length 100, this function can be used in
the following manner:

avergeofA = average(10,A)

11.4.10     Subroutines
A subroutine is similar to a function however, it transfers data to the main program in a slightly
different way. Here is an example of a subroutine which rotates the vector ‘r’ by ‘theta’ degrees.


REAL, INTENT(INOUT)                              :: r(1:2) REAL
:: rnew(1:2) REAL, INTENT(IN)                                 :: theta


r = rnew

11.4. EXECUTION PART                                                                            107

Note that in the declaration statements we have included the comment INTENT(IN), and IN-
TENT(INOUT) for the variables being read in from the main program. The INTENT(IN) state-
ment will produce an error message when you try to compile the code if you try and write over
that variable in the subroutine. The INOUT option implies that the variable may be changed
during the functioning of the subroutine. The option OUT can be used if the variable is the
product of the subroutine. It is not necessary to use the INTENT option in subroutines, how-
ever, sometimes it can help to prevent errors - particularly if you end up trying to write over a
variable which you did not intend to change in the subroutine.
   To use the subroutine in the main program we use the statement (assuming we have a 2x2
vector ’rvector’ which we would like to rotate by 1.53 radians) :

CALL rotate2D(rvector,1.53)

Note that subroutines are not declared in the main program.

11.4.11   Common Variables / Modules
There are times when it is convenient to have variables which are passed between the main
program and the subroutines without needing to list them all in the CALL line. In Fortran 77
this was done using a COMMON statement, but in Fortran 90, it is usually done with a Module.
Here is an example of a module:


INTEGER, PARAMETER                                        :: d = 3 REAL, PARAMETER
:: temperature = 20.0 REAL                                      :: array1, array2


In the main program and any subroutines and functions that need access to these variables the
following line is included right after the PROGRAM, SUBROUTINE, or FUNCTION statements
and before the IMPLICIT NONE

USE input

There are many other ways in which modules are used, however they will not be discussed

11.4.12   ’Save’ variables in subroutines
There are times when it is convenient in a subroutine for a variable to be remembered by the
subroutine for the next instance when you call it. For instance, if you are calculating averages, it
can be useful to know the number of times you called the subroutine function which calculates
the averages. This can be done by using the SAVE attribute in the initialization the that variable.
For instance:


INTEGER, SAVE                 :: measurements = 0 REAL, INTENT(IN)     :: r
REAL, INTENT(OUT)                 :: ravg REAL, SAVE           :: rtot =
108                               CHAPTER 11. APPENDIX: INTRODUCTION TO FORTRAN

measurements = measurements+1

rtot = rtot + r ravg = rtot/float(measurements)


This code can be used to calculate the average value of some variable r. Each time the code
enters the routine the integer ’measurements’ is increased by 1 and there is a new rtot, which
is the sum of it’s previous value and the current r. To see how this works, try using the above
subroutine with the following main program:


INTEGER                :: i REAL                   :: x, xavg

DO i=1,10 x = float(i) CALL statistics(x,xavg) print*, x, xavg END


11.5      Useful References
      • “ Fortran 90 Programming” by T.M.E. Ellis, I.R. Philips, and T.M. Lahey. This is a good
        introductory book on Fortran 90 for both people with and without previous programming


Chapter 12

Appendix: Introduction to C

12.1      Introduction
12.1.1    The language
C is a general-purpose computer programming language developed in 1972 by Dennis Ritchie
at the Bell Telephone Laboratories for use with the Unix operating system.
    Although C was designed for implementing system software, it is also widely used for de-
veloping portable application software.
    C is one of the most popular programming languages. It is widely used on many different
software platforms, and there are few computer architectures for which a C compiler does not
exist. C has greatly influenced many other popular programming languages, most notably C++,
which originally began as an extension to C.
    C is a relatively simplistic language, lacking some of the features of newer languages. How-
ever, for calculation and simulation, these features can be done without, as the needed opera-
tions there are mostly straightforward as well.

12.1.2    A first program
As is usual when starting out with a programming language, here is a simple program:
# i n c l u d e <s t d i o . h>

main ( )
  p r i n t f ( ” Hello , world ! \ n ” ) ;
    To run this program, we will have to compile it into an executable file. For this, we save the
file as hello.c and compile the code with a compiler called gcc. We can run it from a terminal
gcc h e l l o . c −o h e l l o
   The gcc command is always supplied with the file(s) to compile, and a number of compiler
options. The only option here is ”-o hello”, which tells the compiler the program it creates
should be called hello.
   If you have entered this correctly, you should now see a file called hello. This file is the
binary version of your program, and when run should display ”Hello, world!”

110                                         CHAPTER 12. APPENDIX: INTRODUCTION TO C

    Here is an example of what compiling and running looks like when using a terminal on a
unix system. ls is a common unix command that will list the files in the current directory, which
in this case is the directory ’progs’ inside the home directory (represented with the special tilde,
∼ , symbol).
˜/ progs$ l s
hello . c
˜/ progs$ gcc −o h e l l o h e l l o . c
˜/ progs$ l s
hello hello . c
˜/ progs$ . / h e l l o
Hello , world !
˜/ progs$
    We’ll take a quick look at the lines of the program individually. The first line tells the com-
piler to include the standard header file <stdio.h> There are many such header files, each
containing standard function descriptions which can be used in your program. In this case,
<stdio.h> supplies a number of standard input and output functions, such as printf.
    The line ’main()’ shows the start of the function ’main’. Each C program has a ’main’ func-
tion. Generally, the ’main’ function is where a program begins. The body of the function is
enclosed by braces ({}), which show where the function starts and ends.
    The line that does the actual printing is below this: the printf function is a function for
printing formatted strings to the screen, and is defined in <stdio.h>. It can print strings, but
also the value of any variables in your code. See the section on printing later on for details. The
”\n” in this case is an escape sequence which adds a new line at the end of the printed text.
    Every command in C is ended by a semicolon (;). This is easy to forget when starting out,
and can cause confusing error messages from the compiler. This also means that starting a new
line in your code if a statement is very long is always possible.
    Note that many of the examples given in this chapter are not full programs. Often just a few
lines or functions are given. You may have to add a main function and include the right libraries
to try them out.

12.1.3   Operators
The operators you will need for basic programming are fairly straightforward:
x   = 1;           // S e t x t o one
x   = 1 + 2;       //Addition
y   = 5 − x;       // S u b t r a c t i o n
z   = y ∗ x;       // M u l t i p l i c a t i o n
x   = x / 2;       // D i v i s i o n
y   = 10 % 2 ;     //10 mod 2
x   += 6 ;         //Add 6 t o x .
x   /= 3 ;         //Divide x by t h r e e .
x   −= 2 ;         //∗= and −= work s i m i l a r l y
x   = (1 + 2) ∗ 3;       // B r a c k e t s
Spaces in these statement are not important at all, but may make your code look nicer.

12.1.4   Coding
When coding, there are a few things to keep in mind.
12.2. HEADER FILES                                                                             111

    First of all, whitespace generally does not affect the code. This means you can put extra
spaces, tabs, and line breaks anywhere you want.
    It is generally useful to use indentation to show the structure of your program. Use spaces
to move lines within a function, loop, or if-statement to the right. This way, you can easily see
what code is grouped together.
    If a line of code contains two slashes in a row (//), then anything after this is designated as
a comment, and ignored by the compiler. This allows you to leave notes explaining what your
code is doing, or to temporarily disable a line of code. In addition, anything between /* and */
is also ignored, even if this spans multiple lines.
# i n c l u d e <s t d i o . h>

main ( )
  int i ;
  f o r ( i = 1 ; i < 2 0 ; i ++) // S t a r t loop
      i f ( i % 2 == 0 )                 // i s i d i v i s i b l e by two ?
         p r i n t f (”%d i s even ! \ n ” , i ) ;
         p r i n t f (”%d i s odd ! \ n ” , i ) ;
The following code is equivalent, but a lot less easy to read...
# i n c l u d e <s t d i o . h>

main ( ) {
int i ;
f o r ( i = 1 ; i <20; i ++)
i f ( i %2==0) p r i n t f (”%d i s even ! \ n ” , i ) ;
e l s e p r i n t f (”%d i s odd ! \ n ” , i ) ; }

12.2      Header files
As shown in the first little program, it is often necessary to include header files at the start of
your program to use some of the built-in functions in C. If there are multiple files to include,
you can simply use multiple #include lines.
    <stdio.h> includes all standard input and output functions. This includes printing to the
screen, reading from and writing to files and getting input from the keyboard. You will likely
need this in every C program you write.
    <stdlbi.h> is the standard library, including a number of standard functions. The impor-
tant ones for this course are the random number generator ’drand48()’, and memory allocating
112                                            CHAPTER 12. APPENDIX: INTRODUCTION TO C

    <string.h> includes a number of functions for manipulating strings.
    <math.h> includes a large number of mathematical functions, such as ’exp’, ’sin’, ’sqrt’, etc.
Also useful is the function ’pow(x,y)’, which returns xy . NOTE: When including <math.h> you
also have to let the compiler know to link this library. You do this by adding the compiler option
”-lm” to your compile command.
gcc −lm prog . c −o prog
      Some examples of the math functions:
# i n c l u d e <math . h>

main ( )
  double x , y , z ;
  double p i ;
  p i = atan ( 1 ) ∗ 4 ;           // C a l c u l a t e p i
  x = exp ( p i ) ;
  y = sinh ( pi ) ;
  z = pow( y , x ) ;

12.3      Variables
12.3.1    Data types
There are many types of variables in C. The most useful ones include:
    int: This is an integer, either positive or negative.
    double: A double-precision floating point number, used for any real number you might want
to use.
    float: A single-precision floating point number. For most calculations that require any accu-
racy, it is better to use doubles.
    char: A single character, such as ’a’.
    You can change the data types of a value by assigning it to a variable with the new data type,
or by directly telling the program it should treat a variable as a variable of a different type. This
is done by placing the new variable type between brackets (see below). For example, if x is a
double, ’(int) x’ instead uses x as an integer, rounding down.
double x =       3.5;
int i = x ;                // i i s now 3
double a =      i;               //a i s now 3 . 0
double k =      i / 2;           // I n t e g e r d i v i s i o n : k i s now 1
double m =      ( double ) i / 2 ;
                    // T r e a t i as a double : m i s 1 . 5
p r i n t f (”%d\n ” , ( i n t ) x ) ; // P r i n t s x as an i n t e g e r , so 3 .

12.3.2    Declaring variables
Before using a variable anywhere in your program, it must be declared. Declaring variables is
the way in which a C program shows the number of variables it needs, what they are going to
be named, and how much memory they will need. After the declaration, the variable does not
12.3. VARIABLES                                                                                    113

have a well-defined value yet. You should only use the value of a variable after it has actually
been given one. In the example below, the print statement tries to print the value of y, but since
it has not been given a value, the outcome of this is unpredictable. Do not assume the value is
# i n c l u d e <s t d i o . h>

int j = 5;                 //Global v a r i a b l e

main ( )
  int i ;              //Some d e c l a r a t i o n s
  double x , y ;
  double a = 0 . 5 ;         // D e c l a r a t i o n and i n i t i a l i z a t i o n
  i = 3;               // I n i t i a l i z a t i o n o f i
  i ++;                //Add one t o i , i i s now 4
  x = i ∗ j ∗ a;
  p r i n t f ( ” i = %d , x = %l f , y = %l f \n ” , i , x , y ) ;
    Note that variable names (and function names) are case-sensitive.
    Any variables declared outside any function, just below the #include lines, is considered
a global variable. This means that any function in your program can access and modify this
variable. If a variable is declared within a function, it is only available in that function. This also
means you can use a variable with the same name in another function.

12.3.3     Arrays
In many cases, you will need to work with lists of values, such as coordinates, rather than single
values. A list in C is called an array, and can be declared as follows:
i n t values [ 3 ] ;
values [ 0 ] = 1 ;
values [ 1 ] = 2 ;
values [ 2 ] = values [ 0 ] + values [ 1 ] ;
// v a l u e s [ 3 ] = 4 ; i s not allowed ! !
Array indices always start at zero, and reading or writing outside the array bounds can cause
your program to crash. If this happens, you will usually see the message ”Segmentation fault”.
Note that declaring an array as a global variable can only be done if the size of the array is a
# i n c l u d e <s t d i o . h>
# d e f i n e SIZE 5             //R e p l a c e s a l l o c c u r e n c e s o f SIZE
                     // i n t h e program by 5

int   s = 5;
int   values1 [ 5 ] ;             //This works
int   v a l u e s 2 [ SIZE ] ;    //This works
int   values3 [ s ] ;             //This doesn ’ t work

main ( )
114                                            CHAPTER 12. APPENDIX: INTRODUCTION TO C

    //Do s t u f f
Operations on arrays can generally not be done in one command in C. Normally, array opera-
tions are done in loops.

12.4     Loops and conditions
12.4.1    if
With if statements you can create code that is only executed if certain conditions are met. You
can use ’else’ if you also want other code to execute if the if statement turns out to be false. The
code to execute is enclosed by curly braces, though these can be omitted if you only want to
have one statement depend on the condition.
func ( i n t i )
  i f ( i == 0 )
     p r i n t f ( ” This number i s zero ! \ n ” ) ;
  else i f ( i > 0)
     p r i n t f ( ”A p o s i t i v e number\n ” ) ;
     p r i n t f ( ” Negative ! \ n ” ) ;
  i f ( i ! = 1 && i ! = 8 ) p r i n t f ( ” Not one or e i g h t . \ n ” ) ;
     // && means ”and”
     // | | means ” or ”
     // ! = means ” not equal ”
Note that to check if two numbers are equal, you need two equals signs (==).

12.4.2    for
To execute a block of code multiple times, you can use for loops. These generally look like this:
int i ;
f o r ( i = 0 ; i < 1 0 ; i ++)                // i ++ adds one t o i
   //Code t o run 10 ti me s
    p r i n t f ( ” i i s now %d\n ” , i ) ;
The for statement has three parts. The first part is executed when the loop begins, and sets the
starting value for the loop index. The second part is the condition for the loop: it keeps going
until this condition is untrue. The last part is executed at the end of every cycle in the loop.
In the code above, this causes the integer i to run from 0 to 9, and the body of the for loop is
executed for each value of i. This is very useful for doing operations on arrays.
12.5. FUNCTIONS                                                                             115

12.4.3   while
Another common type of loop is the while loop. While loops simply continue looping until a
certain condition evaluates to false.
int i = 1000;
int n = 0;
while ( i >= 1 ) // i ++ adds one t o i
  i = i / 2;
  n += 1 ;   //adds one t o n

12.4.4   Escaping loops
Sometimes you want to get out of a loop before it would normally end. In this case, the break
statement will break out of the innermost loop. The continue statement will just move on to the
next iteration of the loop.
//Search f o r s q u a r e s i n a s i l l y way :
int i ;
int notthis = 3;                 //This number we s k i p
f o r ( i = 0 ; i < 1 0 ; i ++)
    if     ( i == n o t t h i s ) c o n t i n u e ;
    f o r ( j = 0 ; j < i ; j ++)
        i f ( j ∗ j == i )
           p r i n t f (”%d squared = %d\n ” , j , i ) ;
           break ;        //Move on t o t h e n ex t i

12.5     Functions
In C programming, all executable code resides within a function. A function is a named block
of code that performs a task and then returns control to a caller. A function is often executed
(called) several times, from several different places, during a single execution of the program.
After finishing a subroutine, the program will branch back (return) to the point after the call.
    The following function, for example, calculates the square of a number:
i n t square ( i n t x )
     i n t sq ;
     sq = x ∗ x ;
     r e t u r n sq ;
116                                          CHAPTER 12. APPENDIX: INTRODUCTION TO C

printsquares ( )
  int i , isquare ;
  f o r ( i = 0 ; i < 1 0 ; i ++)
      i s q u a r e = square ( i ) ;
      p r i n t f (”%d squared i s %d\n ” , i , i s q u a r e ) ;
Functions can return a value, but they don’t necessarily have to do so. If a function returns a
value, it will state what type of value it returns at the start of the function, before the name. The
above function square returns an integer, while printsquares does not return anything.
    A function can also take any number of arguments, which can be used within the function.
The function square, for example, takes one integer as its argument. In C, a function normally
is only able to change any of its arguments locally: try the following program.
# i n c l u d e <s t d i o . h>
main ( )
    int x = 3;
    func ( x ) ;
    p r i n t f ( ” In main , x i s %d\n ” , x ) ;

func ( i n t x )
  x = 2 ∗ x;
  p r i n t f ( ” In func , x i s %d\n ” , x ) ;
The variable x is set to 3 in the main function. The function func receives this value 3, doubles it,
and prints the result. However, the value of x in the main function has not changed. If you want
to use results from a function, you have three options. The first way is to let the function return
a value, as shown above. Another option is the use of global variables. Any changes you make
to global variables in any function will be seen by any other function using those variables. The
last option is a bit more complicated, and requires the use of pointers.

12.5.1   Pointers
Pointers are a slightly more tricky concept in C, but they have many uses. We won’t need them
in many cases in C, so this section only covers the basics. A pointer is a variable that holds the
memory address of another variable. To work with pointers, two operators are important. The
operator & means ’address of’, and a * means ’value at’. The following example shows the use
of a pointer to an integer:
i n t ∗ pt ;        //pt i s a p o i n t e r t o an i n t e g e r
int i = 5;
pt = &i ;           // S e t pt equal t o t h e address o f i
∗ pt = 2 ∗ ∗ pt ; //Double t h e value a t pt ;
p r i n t f ( ” i i s now %d\n ” , i ) ;
12.6. PRINTING AND READING                                                                   117

   Of course, we could use ’i = 6’ instead, so what’s the use? It becomes a lot more useful if we
pass the address of i to a function:
# i n c l u d e <s t d i o . h>
main ( )
    int i = 5;
    t w i c e (& i ) ;
    p r i n t f ( ” i i s now %d\n ” , i ) ;

t w i c e ( i n t ∗ pt )
   ∗ pt = 2 ∗ ( ∗ pt ) ;
This way, we can change the value one or more parameters of a function, rather than just return-
ing one value and/or changing global variables. There are several built-in functions that work
the same way, such as fscanf.

12.5.2   Static variables
When leaving a function, any local variables will be forgotten. The next time you enter the
function, it will not be influenced by any previous function calls. You can make your program
remember local variables by declaring them with the static keyword:
count ( )
  s t a t i c i n t counter = 0 ;
  int i = 0;
  i ++;
  c o u n t e r ++;
  p r i n t f ( ” Function c a l l e d %d times , but i i s s t i l l %d\n ” ,
                     counter , i ) ;

12.6     Printing and reading
Most programs will at some point write something to the screen or to a file. You can write to
the screen this using a printf (’print formatted’) statement, as shown in many examples in this
chapter. The printf arguments consist of one string to print, and a number of variables to fill in
at specified points in the string. As mentioned before, ’\n’ is a line break.
double x = 1 . 0 / 3 . 0 ;
int i = 4;

p r i n t f (”%d , %l f \n ” , i , x ) ; //%d : i n t e g e r
                          //% l f : double
       // l f s t a n d s f o r ’ long f l o a t ’
118                                               CHAPTER 12. APPENDIX: INTRODUCTION TO C

You can format the printed numbers in various ways. For example ’%.2lf’ will round a double
to two digits after the decimal point, and ’%3d’ will add spaces in front of the integer to make
sure the total width printed is at least 3.
    You can also print to files, using fprintf. Files have to be opened first, and closed after you
have finished writing. File pointers have the type ’FILE*’.
FILE ∗ f i l e ;
f i l e = fopen ( ” out . dat ” , ”w” ) ;

f p r i n t f ( ” This t e x t ends up i n t h e f i l e \n ” ) ;

fclose ( f i l e ) ;
The above code opens a file named out.dat, and writes to it. The ’w’ passed to fopen tells it to
open the file for writing. This means the file is emptied before writing to it. If you replace the
’w’ with an ’a’, it will append to the file instead, leaving anything already there in place. You
can also open for reading:
# i n c l u d e <s t d i o . h>
# d e f i n e MAX 1000

double x [MAX] ;
double y [MAX] ;
i n t npart ;

main ( )
  int i ;
  readparts ( ) ;
  f o r ( i = 0 ; i < np art ; i ++)
      p r i n t f (”%2d : (% l f , %l f ) \ n ” , i , x [ i ] , y [ i ] ) ;

readparts ( )
  int i ;
  i n t numbersread ;
  FILE ∗ f i l e ;
  f i l e = fopen ( ” i n . dat ” , ” r ” ) ;
  i f ( f i l e == NULL) //fopen r e t u r n s NULL i f i t f a i l s
      p r i n t f ( ” F i l e not found ! \ n ” ) ;
      return ;
  f s c a n f ( f i l e , ”%d\n ” , &n par t ) ;
  i f ( npar t > MAX) np art = MAX; //At most MAX p a r t i c l e s
  f o r ( i = 0 ; i < np art ; i ++)
      numbersread= f s c a n f ( f i l e ,”% l f %l f \n” ,&( x [ i ] ) , & ( y [ i ] ) ) ;
12.7. EXAMPLE PROGRAM                                                                           119

       i f ( numbersread ! = 2 ) //Something wrong?
          p r i n t f ( ” Read e r r o r ! \ n ” ) ;
          break ;
    fclose ( f i l e ) ;

/∗ Example i n . dat :
1.0 1.0
2.0 2.0
The above program reads in a file called in.dat, and expects it to contain a line with the number
of particles, and then lines with xy-coordinates for each particle. If the file isn’t there, or if it
can’t read the right numbers, it will complain. It then prints the result. The function fscanf
works very similar to fprintf, except that it takes in the address of any variables it should assign
data to, using the & operator. It returns the number of values read, which can be used to check
if anything went wrong.
    In addition to printing to the screen and to files, it can also be useful to print to a string,
using the function sprintf. Strings are just arrays of characters.
char s t r [ 2 0 ] ;
int x = 3;
s p r i n t f ( s t r , ”x%d ” , x ) ;
p r i n t f ( ” S t r i n g : %s \n ” , s t r ) ;

12.7      Example program
As a conclusion, let’s look at an example.

12.7.1    Step 1: Specify the Problem
We would like to write a program which evaluates the roots of the quadratic equation ax2 +bx+
c = 0 which are given by the following formula:
                                          −b ± b2 − 4ac
                                     x=                                                  (12.1)
Note that there will be three cases of interest, when the discriminant (D = b2 − 4ac) is 0, positive
or negitive.

12.7.2    Step 2: Write a Pseusocode (or Algorithm)
    • read in the values of a, b, and c for which you would like to know the roots.

    • if a is zero stop - in this case we do not have a quadratic equation

    • evaluate D
120                                              CHAPTER 12. APPENDIX: INTRODUCTION TO C

      • if D = 0 then the root is   2a
                                                 −b± D
      • if D > 0 then there are two real roots     2a
                                                         −b±i −D
      • if D < 0 then there are two imaginary roots         2a

      • print the solution

12.7.3     Step 3: Write the Program

# i n c l u d e <s t d i o . h>
# i n c l u d e <math . h>

main ( )
  double a , b , c ;
  double D, r e a l , imag ;
  double root1 , r o o t 2 ;
  double d e l t a = 0 . 0 0 0 0 0 0 0 0 0 1 ; //rounding e r r o r
  p r i n t f ( ” What a r e a , b , and c ?\n ” ) ;
  s c a n f (”% l f \n%l f \n%l f ” , &a , &b , &c ) ;
  p r i n t f (”% l f x ˆ 2 + %l f x + %l f \n ” , a , b , c ) ;
  i f ( a == 0 )
      p r i n t f ( ” Not a q u a d r a t i c e q u a t i o n ! \ n ” ) ;
      return ;
  D = b∗b −4∗a ∗ c ;
  i f ( D > −d e l t a && D < d e l t a )
      r o o t 1 = −b / ( 2 ∗ a ) ;
      p r i n t f ( ” There i s only one r o o t : %l f \n ” , r o o t 1 ) ;
  e l s e i f (D > 0 )
      r o o t 1 = (−b − s q r t (D) ) / ( 2 ∗ a ) ;
      r o o t 2 = (−b + s q r t (D) ) / ( 2 ∗ a ) ;
      p r i n t f ( ” There a r e two r o o t s : %l f , %l f \n ” , root1 , r o o t 2 ) ;
      r e a l = −b / ( 2 ∗ a ) ;
      imag = s q r t (−D) / ( 2 ∗ a ) ;
      p r i n t f ( ” There a r e two complex r o o t s : %l f + %l f i , %l f − % l f i \n ” ,
                              r e a l , imag , r e a l , imag ) ;
Chapter 13

Appendix: Short Introduction into

13.1    Getting started
To start Mathematica, double-click on the Mathematica icon. When you start Mathematica you
will be presented with an empty window similar to the window shown here:

You can open a new window in your copy of Mathematica by selecting New from the File menu.
To enter your input, select the new window by clicking on it, and start typing. After entering
your input, the window will look something like this:

To evaluate your input, press Shift-Enter (press the enter key while holding down the shift key).
When the calculation is finished, the result is displayed just below your input:

When you are finished with your calculations you can end your Mathematica session (which
will close this window) by selecting Quit from the File menu. To stop a calculation, select inter-
rupt or abort Evaluation or Quit Kernel from the Evaluation tab.


13.2        Help
  1. A good tutorial can be found here:

  2. For help with individual functions, for example with a function whose name begins with
     Expression, type:
           ??Expression   for more help
           ?E*            for help on all functions beginning with E.

  3. For help on mathematical functions, use the Function Browser under the Help menu.

  4. With Esc you can type in special characters like by typing Esc pi Esc.

13.3        Elementary calculations
13.3.1      Arithmetic
+ is the sum of two numbers, while - is minus

      2.3 + 5.63 - 3.00


/ stands for division and ˆ for power



Spaces stand for multiplication, but you can also use *



//N gives an approximate numerical value and % denotes the previous output. %% the output
before this, et cetera. %n denotes Out[n].
13.3. ELEMENTARY CALCULATIONS                                                            123


   1 267 650 600 228 229 401 496 703 205 376



   1/3 + 1/16 //N


13.3.2   Some mathematical functions
Sqrt[x] gives the square root of x



Mod[m,n] gives the remainder on division of m by n.



Random[] gives a uniformly distributed random number in the range 0 to 1.



RandomReal[{a,b}] gives a uniformly distributed real random number in the range a to b



RandomInteger[{a,b}] gives a uniformly distributed random integer in the range a to b



IntegerPart[x] gives the integer part of x.



13.3.3      Defining Variables
When you do long calculations, it is often convenient to give names to your intermediate results.
Just as in standard mathematics, or in other computer languages, you can do this by introducing
named variables.

This sets the value of the variable x to be 5.


13.3. ELEMENTARY CALCULATIONS                                               125

Whenever x appears, Mathematica now replaces it with the value 5.



To unassign x use:

   x =. ;


   x = x

You can perform two evaluations by separating them with a semicolon ” ; ”

13.3.4      Making tables
This gives a table of the values of i2 , with i running from 1 to 6.

   x=Table[iˆ 2, {i, 1, 6}]

   {1, 4, 9, 16, 25, 36}

//TableForm displays the table in tableform

   % //TableForm

x[[4]] extracts the fourth element



This creates a 2x3 table, and gives it the name m.

      m = Table[i - j, {i, 3}, {j, 2}]


This extracts the first sublist from the list of lists that makes up the table.



This extracts the second element of that sublist.



This displays m in a ”tabular” form.


       0   -1
       1   0
       2   1
Flatten[m,1] flattens table m to level 1.



This is a 2x2 matrix.

      m = {{a, b}, {c, d}}

      {{a, b}, {c, d}}

Here is the first row.



Here is the element m12 .


13.3. ELEMENTARY CALCULATIONS                                      127

Another example:





    1    1   1   1
    1    2   3   4
    2    2   2   2
    1    2   3   4
    3    3   3   3
    1    2   3   4



    1    1
    1    2
    1    3
    1    4
    2    1
    2    2
    2    3
    2    4
    3    1
    3    2
    3    3
    3    4
Here is a xyt-table

   nmove=2; npart=3;

   xyt = Table[{x[[j]],y[[j]],i},{i,nmove},{j,npart}]



13.3.5       For, if, goto, label

For [start,test,incr,body] executes start, then repeatedly evaluates body and incr until test fails
to give True


If [condition,t] gives t if condition evaluates to true



      If [ ipart != j, Print[ipart] ]


Note that ”a = b” should be entered as ”!=” and ”a ≤ b” as ”<=”.
Goto[tag] scans for Label [tag], and transfers control to that Point


      For[i=1,i ≤npart,i++,


      Print[”i smaller than npart”]; Print [i]; Print[npart];

13.3. ELEMENTARY CALCULATIONS                                    129

     i smaller than npart
     i smaller than npart
     i smaller than npart

13.3.6   Basic Plotting
To plot a graph of sin(x) as a function of x from 0 to 2π, use

   Plot[Sin[x], {x, 0, 2Pi}]

You can give a list of functions to plot.

   Plot[{Sin[x], Sin[2x], Sin[3x]}, {x, 0, 2Pi}]

      Plot[{Sin[x], Sin[2x], Sin[3x]}, {x, 0, 2Pi}, Frame –>True]

      Plot[{Sin[x],Sin[2 x],Sin[3 x]},{x,0,2 Pi}, PlotRange –>{{0,1},{0,1.5}}]

Graphics[{Pink,Disk[{0,0},5]}] plots a disk in the color pink at position {0,0} with a radius of 5.
PlotRangeClipping allows graphics objects to spread beyond PlotRange:

      Graphics[{Pink,Disk[{0,0},5]},Frame –>True,PlotRange –>4, PlotRangeClipping –>False]
13.3. ELEMENTARY CALCULATIONS                                                     131

   Graphics[{Pink,Disk[{0,0},5]},Frame –>True,PlotRange –>4, PlotRangeClipping –>True]

AspectRatio uses numerical values to specify the height-to-width ratio:

   Graphics[Circle[ ],AspectRatio –>Automatic]

   Graphics[Circle[ ],AspectRatio –>0.5]

Graphics[primitives,options] represents a two-dimensional graphical image.


To top