VIEWS: 5 PAGES: 36 CATEGORY: Computers & Internet POSTED ON: 5/1/2010
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 ﬁrst 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 ﬁrst 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) 2a 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 −b • if D = 0 then the root is 2a √ −b± D • if D > 0 then there are two real roots 2a 97 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, 2007 !--------------------------------------------------------------------------- PROGRAM QuadraticEquationRoots ! INITIALIZATION PART IMPLICIT NONE 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 ELSE 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 END IF ELSE PRINT*, ’PROBLEM: this is not a quadratic equation since a is 0’ END IF END PROGRAM QuadraticEquationRoots 11.2 Program Components A program usually consists of two parts: • an Initialization part where all the variables are deﬁned • a Execution part which is the main body of the code where data is read in, manipulated, etc. 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 instance, ! There are three possible cases, D =,>,< 0. and 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 deﬁned. ) • INTEGER :: var1 (deﬁnes “var1” as an integer such as -1, 10, 12838) • INTEGER, PARAMETER :: N = 10 (deﬁnes N to be a parameter ALWAYS equal to 10, you cannot change the value of N later on) • REAL :: var2 (deﬁnes “var2” as a real number such as 2.18847) • REAL, PARAMETER :: Pi = 3.1414 (deﬁnes ‘Pi’ to be a real number equal to 3.1415) • CHARACTER :: initial ( deﬁnes ”initial” as a single alphanumeric character such as “a” or “Z”, etc.) • CHARACTER(LEN = 12) :: name ( deﬁnes “name” to be a string containing 12 alphanu- meric characters such as “Program test”.) • LOGICAL :: raining ( deﬁnes “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 ﬁrst letter in the variable name. If the ﬁrst 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 ﬁle formating. Formatting will be discussed in a following section. 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 ﬁles Before reading or writing to a ﬁle, the ﬁle must be ﬁrst opened. To open a ﬁle called ‘ﬁle- name.dat’ we use the following statement: OPEN(21,FILE=’filename.dat’) The ﬁrst number (in the above example 21) gives us a way to identify the ﬁle 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 ﬁle is to be opened, for instance, it is possible to append data to the end of a ﬁle without writing over the ﬁle itself. To do this you can include a POSITION comment in the OPEN statement as follows: OPEN(21,FILE=’filename.dat’,POSITION=’Append’) There are many other options available when opening ﬁles which are not described here. See any text on Fortran 90 for a complete list of the options. After opening a ﬁle you can read data from the ﬁle. If “ﬁlename.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) END DO It is also possible to write to the ﬁle ‘ﬁlename.dat’. To write the contents of the array variable1 to this ﬁle use the statement: 11.4. EXECUTION PART 101 DO i=1,N WRITE(21,*) variable1(i) END DO When ﬁnished with a ﬁle it can be closed with the statement CLOSE(21) Formatting While in the previous sections the formatting of the ﬁle was done automatically, it is also possi- ble to specify the formatting when reading from or writing to a ﬁle. 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 spaces. 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 4.9 / 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- ment.) 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’ ELSE PRINT*, ’i is greater than or equal to 20’ END IF 11.4.6 DO Statements There are two types of DO loops, ones which execute the loop a ﬁxed number of times and ones which are terminated on reaching some condition. An example of a “ﬁxed” loop is given by the following: 11.4. EXECUTION PART 103 DO i=1,10 PRINT*, ’i is’, i END DO 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 END DO Both programs will print the numbers from 1 to 10, and then exit the loop. However, depending on what you are doing, you may ﬁnd 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 deﬁning 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) :: E 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 speciﬁc 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 ﬁrst 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 ﬁnd 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,:) - Particles(3,:) 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 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. 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 B(1:N,1:M). 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 is: 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: DEALLOCATE(ravg) 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 http://publib.boulder.ibm.com/infocenter/pseries/v5r3/ index.jsp?topic=/com.ibm.xlf101a.doc/xlflr/minloc.htm 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) END DO average = sum / REAL(N) END FUNCTION 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. SUBROUTINE rotate2D(r,theta) IMPLICIT NONE REAL, INTENT(INOUT) :: r(1:2) REAL :: rnew(1:2) REAL, INTENT(IN) :: theta rnew(1)=r(1)*cos(theta)-r(2)*sin(theta) rnew(2)=r(1)*sin(theta)+r(2)*cos(theta) r = rnew END SUBROUTINE 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: MODULE input IMPLICIT NONE INTEGER, PARAMETER :: d = 3 REAL, PARAMETER :: temperature = 20.0 REAL :: array1, array2 END MODULE input 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 here. 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: SUBROUTINE statistics(r,ravg) IMPLICIT NONE INTEGER, SAVE :: measurements = 0 REAL, INTENT(IN) :: r REAL, INTENT(OUT) :: ravg REAL, SAVE :: rtot = 0.0 108 CHAPTER 11. APPENDIX: INTRODUCTION TO FORTRAN measurements = measurements+1 rtot = rtot + r ravg = rtot/float(measurements) END SUBROUTINE 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: PROGRAM test IMPLICIT NONE INTEGER :: i REAL :: x, xavg DO i=1,10 x = float(i) CALL statistics(x,xavg) print*, x, xavg END DO END PROGRAM 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 experience. • http://www.dur.ac.uk/resources/its/info/guides/138fortran90.pdf • http://www.lancs.ac.uk/iss/researchsupport/f90_intro/index.html 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 inﬂuenced 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 ﬁrst 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 ﬁle. For this, we save the ﬁle as hello.c and compile the code with a compiler called gcc. We can run it from a terminal window: gcc h e l l o . c −o h e l l o The gcc command is always supplied with the ﬁle(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 ﬁle called hello. This ﬁle is the binary version of your program, and when run should display ”Hello, world!” 109 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 ﬁles 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 ﬁrst line tells the com- piler to include the standard header ﬁle <stdio.h> There are many such header ﬁles, 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 deﬁned 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. Example: # 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 ) ; } else { 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 ﬁles As shown in the ﬁrst little program, it is often necessary to include header ﬁles at the start of your program to use some of the built-in functions in C. If there are multiple ﬁles 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 ﬁles 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 functions. 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 ﬂoating point number, used for any real number you might want to use. ﬂoat: A single-precision ﬂoating 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-deﬁned 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 zero. # 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 constant: # 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 ” ) ; } else { 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 ﬁrst 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 ﬁnishing 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 ﬁrst 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 inﬂuenced 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 ﬁle. 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 ﬁll in at speciﬁed 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 ﬁles, using fprintf. Files have to be opened ﬁrst, and closed after you have ﬁnished 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 ﬁle named out.dat, and writes to it. The ’w’ passed to fopen tells it to open the ﬁle for writing. This means the ﬁle is emptied before writing to it. If you replace the ’w’ with an ’a’, it will append to the ﬁle 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 : 2 1.0 1.0 2.0 2.0 ∗/ The above program reads in a ﬁle 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 ﬁle 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 ﬁles, 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) 2a 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 −b • 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 ) ; } else { 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 Mathematica 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 ﬁnished, the result is displayed just below your input: When you are ﬁnished 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. 121 122 CHAPTER 13. APPENDIX: SHORT INTRODUCTION INTO MATHEMATICA 13.2 Help 1. A good tutorial can be found here: http://library.wolfram.com/conferences/devconf99/withoff/index2.html. 2. For help with individual functions, for example with a function whose name begins with Expression, type: ?Expression ??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 4.93 / stands for division and ˆ for power 2.4/8.9ˆ2 0.0302992 Spaces stand for multiplication, but you can also use * 234 24 //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 2ˆ100 1 267 650 600 228 229 401 496 703 205 376 N[%] 1.26765x1030 1/3 + 1/16 //N 0.395833 13.3.2 Some mathematical functions Sqrt[x] gives the square root of x Sqrt[3.0] 1.73205 Mod[m,n] gives the remainder on division of m by n. Mod[1,2] 1 Random[] gives a uniformly distributed random number in the range 0 to 1. Random[] 0.227695 RandomReal[{a,b}] gives a uniformly distributed real random number in the range a to b RandomReal[{-0.5,0.5}] 124 CHAPTER 13. APPENDIX: SHORT INTRODUCTION INTO MATHEMATICA -0.427695 RandomInteger[{a,b}] gives a uniformly distributed random integer in the range a to b RandomInteger[{1,100}] 56 IntegerPart[x] gives the integer part of x. IntegerPart[13.45] 13 13.3.3 Deﬁning 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. x=5 5 13.3. ELEMENTARY CALCULATIONS 125 Whenever x appears, Mathematica now replaces it with the value 5. xˆ2 25 To unassign x use: x =. ; 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 1 4 9 16 25 36 x[[4]] extracts the fourth element x[[4]] 16 126 CHAPTER 13. APPENDIX: SHORT INTRODUCTION INTO MATHEMATICA This creates a 2x3 table, and gives it the name m. m = Table[i - j, {i, 3}, {j, 2}] {{0,-1},{1,0},{2,1}} This extracts the ﬁrst sublist from the list of lists that makes up the table. m[[1]] {0,-1} This extracts the second element of that sublist. m[[1,2]] -1 This displays m in a ”tabular” form. TableForm[m] 0 -1 1 0 2 1 Flatten[m,1] ﬂattens table m to level 1. Flatten[m,1] {0,-1,1,0,2,1} This is a 2x2 matrix. m = {{a, b}, {c, d}} {{a, b}, {c, d}} Here is the ﬁrst row. m[[1]] {a,b} Here is the element m12 . m[[1,2]] b 13.3. ELEMENTARY CALCULATIONS 127 Another example: xy=Table[{i,j},{i,1,3},{j,1,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}}} TableForm[xy] 1 1 1 1 1 2 3 4 2 2 2 2 1 2 3 4 3 3 3 3 1 2 3 4 xy=Flatten[xy,1] {{1,1},{1,2},{1,3},{1,4},{2,1},{2,2},{2,3},{2,4},{3,1},{3,2},{3,3},{3,4}} TableForm[xy] 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}] {{{x[[1]],y[[1]],1},{x[[2]],y[[2]],1},{x[[3]],y[[3]],1}}, {{x[[1]],y[[1]],2},{x[[2]],y[[2]],2},{x[[3]],y[[3]],2}}} 128 CHAPTER 13. APPENDIX: SHORT INTRODUCTION INTO MATHEMATICA 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 For[i=0,i<=3,i++,Print[i]] 0 1 2 3 If [condition,t] gives t if condition evaluates to true ipart=5; j=4; If [ ipart != j, Print[ipart] ] 5 Note that ”a = b” should be entered as ”!=” and ”a ≤ b” as ”<=”. Goto[tag] scans for Label [tag], and transfers control to that Point npart=5; For[i=1,i ≤npart,i++, If[i>3,Goto[end]]; Print[”i smaller than npart”]; Print [i]; Print[npart]; Label[end];] 13.3. ELEMENTARY CALCULATIONS 129 i smaller than npart 1 5 i smaller than npart 2 5 i smaller than npart 3 5 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}] 130 CHAPTER 13. APPENDIX: SHORT INTRODUCTION INTO MATHEMATICA 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. Graphics[{Green,Rectangle[{0,-1},{2,1}],Red,Disk[]}]} 132 CHAPTER 13. APPENDIX: SHORT INTRODUCTION INTO MATHEMATICA