VIEWS: 0 PAGES: 6 CATEGORY: Computers & Internet POSTED ON: 5/1/2010
Tutorial introduction to Fortran 90 A simple program The study of physics usually begins with mechanics, so it makes sense to start our study of computational physics in the same place. The intent of this tutorial is to allow you to get some familiarity with the programming process — creating a code ﬁle, compiling and debugging a program — and not to teach you to program. The program examples can be typed in as is, and along the way you will get some idea of the language. Of course, there are other good reasons. Mechanics has evolved methods of extreme sophistication in the 350 years since Newton. One of the main impetuses has been the study of orbital motion in the solar system. But these sophisticated methods are there to provide a good starting point for perturbation treatments; often obscured in the beauty of the formalism is the fact that there are few exact solutions. So numerical treatments — starting already with Euler — have provided a powerful alternative. Let’s start oﬀ with something very simple: dropping a ball from, say, the leaning tower of Pisa. The equation of motion is d2 y m = −mg dt2 with y the height, t the time, m the mass of the ball and g the acceleration due to gravity. Of course, this has the solution at time t, assuming an initial height y0 and velocity zero at time t = 0, 1 y = y 0 − g t2 2 We can begin our programming in the simplest way, using the computer to evaluate this mathematical expression. We will evaluate the height for a set of times 0.1, 0.2, . . . , 1.0 s. We will do this pretty much in the same way we would do the task with a scientiﬁc calculator, so think of the program as a list of instructions for using such a tool. Here is a simple program in Fortran to calculate the height as a function of time: program drop0 2 ! drop0.f90 - rwf ! Simple Fortran program 4 ! straightforward evaluation of expressions 6 implicit none 1 8 real : : g =9.8 ! accel due to gravity (m/s^2) real : : m =1.0 ! mass of ball ( kg) 10 real : : y0 = 1 0 . 0 ! initial height 10 m real : : v0 = 0 . 0 ! initial vertical velocity component (m/s) 12 real : : tmax = 1 . 0 ! drop time (s) real : : dt = 0 . 1 ! time step 14 integer : : Nsteps , step real : : y 16 ! determine number of time steps: 18 Nsteps=tmax / dt+1 ! calculation loop 20 do step =1, Nsteps t = step ∗ dt 22 y = y0 −0.5∗ g∗t ∗∗2 print ∗ , t , y 24 end do 26 end program drop0 (Note that the line numbers printed here are just for reference below, and are not part of the program!) • Comments, which are there to provide explanations of what is going on, start from an exclamation mark and continue to the end of the line. • Note that the program is “bracketed” by program ... and end program ... state- ments. The name given in the statement is arbitrary, and does not have to be the same as the ﬁle. • implicit none is a statement that turns oﬀ some features of older versions of Fortran. • Note that all variables used in the program are declared in the ﬁrst section of the program; some are deﬁned to have initial values. • Lines 8-15 deﬁne various constants and variables. Numerical values could simply be used instead in the appropriate places, but it is usually better to give names to these and not have “magic numbers” ﬂoating about in your program. • This program uses a do loop to iterate over a range of integers which are used to calculate the time and hence the height of the ball. • In line 18, the number of steps to be taken is calculated. • Line 20 introduces the do loop: The variable step is given values in succession from 1 to Nsteps and the indented portion below is executed. Note that the statement ends with a statement end do. Do loops can be contained in other do loops, but cannot overlap. 2 • Lines 21-23 constitute the body of the loop. Note that the indentation is not required, but is useful for seeing the structure of the code. • Note the exponentiation operator is **. • The print statement will print the values of the variables in its argument list on the terminal. The asterisk * is a placeholder for a format statement, and implies a free format. You can specify the format of the data by using a proper format statement instead of the asterisk; this allows you to specify number of signiﬁcant digits, etc. Note the comma after the *. Type this into a ﬁle. The comments are unimportant from the compiler’s point of view and can be omitted — they are there for you to keep track of things. Save the ﬁle with a .f90 extension, e.g. something like drop0.f90. Now, before we can run the program, it must be compiled. We will use the program ifc (Intel Fortran compiler) for this. The command for this is ifc -o drop0 drop0.f90 This takes the ﬁle drop0.f90 and creates an executable program drop0 (that’s the -o drop0 If the command prompt returns after this with no messages all went well. On the other hand it is probable that you got some stream of error messages and warnings. You will have to see what they say and try and ﬁx your code, maybe with some help. OK, it compiles cleanly. Now run it: ./drop0 This line runs the program drop0 that exists in the current directory (that’s what the ./ indicates). (Technical note: the Unix PATH, which lists where executable programs exist, usually excludes the current directory for security reasons, so you need to explicitly set the path). You program runs! You get a list of numbers: do they seem sensible? Now make the problem a little less trivial by including the drag force, so that the equation of motion is d2 y 1 m 2 = −mg + CρAv 2 dt 2 Here, C is a constant; ρ the density of air and A the cross-sectional area of the body. 3 This actually has the solution, 1 1 + (γv0 )2 y = y0 + ln γ2g cosh γgt where γ = CρAv 2 /2mg and y0 is the initial height. This expression has a number of mathematical functions in it, including the square root. There are suitable versions of these supplied by the system. Here is a simple program in Fortran to evaluate this expression. You can enter this, but it’s quicker to copy your ﬁrst program to a new ﬁle (cp drop0.f90 drop1.f90) and modify the new ﬁle. Use what you have. program drop1 2 ! drop1.f90 - rwf ! Simple Fortran program 4 ! straightforward evaluation of expressions 6 implicit none 8 real : : g =9.8 ! accel due to gravity (m/s^2) real : : m =1.0 ! mass of ball ( kg) 10 real : : y0 = 1 0 . 0 ! initial height 10 m real : : v0 = 0 . 0 ! initial vertical velocity component (m/s) 12 real : : tmax = 2 . 0 ! drop time (s) real : : dt = 0 . 1 ! time step 14 real : : C =0.5 ! drag coefficient real : : rho = 1 . 0 ! density of air ( kg/m^3) 16 real : : R =0.1 ! radius of ball (m) real : : pi =3.1415926 18 integer : : Nsteps , step real : : y , t , gamma 20 ! determine number of time steps: 22 Nsteps=tmax / dt+1 ! calculation loop 24 do step =1, Nsteps t = step ∗ dt 26 gamma=sqrt ( C∗ rho ∗ pi∗R ∗ ∗ 2 / ( 2 ∗ m∗g ) ) y = y0 + ( 1 . 0 / ( g∗ gamma ∗ ∗ 2 ) ) ∗ log ( sqrt ( 1 . 0 + ( gamma ∗v0 ) ∗ ∗ 2 ) / ( cosh ( gamma ∗g∗t ) ) ) 28 print ∗ , t , y end do 30 end program drop1 • In many languages you will need to import the math functions from a library or module. Fortran has them built into the language. (You can of course import other things from a library or module). 4 • Note the way of calling a function — this is very similar to obvious mathematical notation. A function typically returns a value. • Note that trigonometric functions expect their argument in radians. • The system does not supply a value for π. The distance fallen in a given time is presumably smaller than in the original code, where there was no air resistance. How much smaller? Calculate it! Let’s add a line to calculate the vacuum value of free fall, say yvac = y0-0.5*g**2. Add this line to the do loop, and alter your program to print it out. Note that you have a new variable name. How about a graph? The simplest way to do this is to use a program like gnuplot. This requires the data to be plotted to be in a ﬁle. The simplest way to do this is to cut and paste from the screen to a text editor, and save the ﬁle, for example as drop1.out. The start gnuplot and plot it: gnuplot> plot ’drop1.out’ Hmm. This only plots one set of data. You have to tell the program which columns to plot: gnuplot> plot ’drop1.out’ using 1:2, ’drop1.out’ using 1:3 This can be abbreviated; you might like to use lines instead of points. gnuplot> plot ’drop1.out’ u 1:2 with lines, ’drop1.out’ u 1:3 w l Now let’s change the program so that it runs until one of the objects hits the ground at y = 0. We can do this by replacing the do loop with a do-while loop: do while ( yvac > 0 . 0 ) 2 step = step + 1 While we’re at it change the distance fallen to 50 m. Now, when we run the program, the data scrolls out of view. It’s time to write it to a ﬁle from within the program. We can do so by ﬁrst opening a ﬁle open(1, file=’drop2.out’) 5 Note that the ﬁle name is a character string enclosed in quotes. the “1” is a logical unit used to refer to the ﬁle. To write to the ﬁle we can replace the print statements with: write(1,*) t,y,yvac Here the “1” stands for the ﬁle; the “*” is still a placeholder for a format. Note that there is no comma before the t. Finally, we can close the ﬁle (after we have written everything to it. close(1) Plot the results using gnuplot. 6