Aleksandar Donev - Dr. Phillip Duxbury
October 2000
Since some of you had problems organizing your Fortran solution to this problem, here are some helpful comments. However, always try to solve the problem by yourself first. Some of you had very nice ideas about how to organize the program, and it is by experimenting with these ideas that you learn most. So please read this section after you have thought about the problem yourself! It is our wish to have everyone working at their own pace and level of experience so that they are challanged, but not frustrated. This is more important than all of you finishing the same worksheet and submitting the same (boring) solution.
In the next three worksheets we will not introduce any new Fortran material. We will in fact use all of the things that we learned with the error function project to solve some numerical problems in electrostatics. Some of the things you learned so far are essentials in Fortran 90 and you need to understand them well. These include variable declaration and basic I/O, the concept of a main program and what it does, what modules are and how we use them to encapsulate (organize and protect) global (module) variables and procedures, how to declare procedures and their arguments (this particularly encompases understanding which data are private to a procedure and which are global and why), how to declare and use simple one and two dimensional (rank-1 and rank-2) arrays, and how to read the documentation (interfaces) for a ready-made library and use the library to solve a problem. You may find that additional Fortran 90 features (un)documented in the manual may make life much easier, but it won't be necessary to use them.
The physics applications that we will explore in these final worksheets will all have the same global structure for the Fortran code. The main computational routines, be this plotting or numerical integration or numerical solution of differential equations, will be ready-made and provided for you as a very easy-to-use Fortran 90 module with good documentation and examples. These mini libraries know nothing about physics, they are computational libraries! If you actually want to learn how these routines work and how to write them, take Numerical Analysis or PHY480.
Now, all of these ready-made routines will accept some input argument from
the user that contains the physics of the process. We have chosen to make
this a procedure (usually a FUNCTION) that calculates some physical
quantity (this is a frequent situation). For example, in today's worksheet
you will need to write a function V that calculates the electric potential
at a given point
,
next time you will need to write a
function that evaluates the Bio-Savart contribution to the magnetic field at
a given point, after that you will need to write a function that evaluates
the force on a particle moving in an electromagnetic field, etc.
Please review your basic electrostatics and magnetostatics for this before
coming to lab.
The main point is that this function that you will write will contain all the relevant physics to the problem. As such, this function should be placed in a separate module which should also contain in it the relevant physical variables (such as charges, currents, fields, etc.) that this procedure will need. Then the main program will USE this module and assign values to these physical variables (say the user will enter them). You should always make only those variables that both the procedure and the main program access global (module) variables--the rest should be private to the function or the main program. Finally, the main program will call the calculation/plotting ready-made routines and print the results if necessary.
We hope this gives you an idea of how this works in Fortran. Mathematica is a great tool, but most people use it without ever understanding its inner workings and thus use it inefficently. But with Fortran everything is transparent and you need to have a clear picture of what and why you are doing. Finally, we should emphasize that we chose to stick with the methodology outlined in this section in all of the remaining worksheets not because it will always be the best choice, but simply so that we do not overload you with many different strategies.
The first and most challenging task in this worksheet is to write the module that will contain the function V(x,y). First, you need to decide what variables this procedure will need other than its input arguments x and yand decide on how you are going to represent these variables. Surely the program will need to know the number of charges, their magnitudes and their positions. The number of charges n_charges is a simple integer, the magnitudes can be easily represented by a simple one-dimensional array of real numbers q(n_charges), but the the positions of the charges give you more room to experiment. For example, you can have two (three) arrays holding the x and y coordinates of each charge, x(n_charges) and y(n_charges), or you can merge these two into a two-dimensional array r(2,n_charges), where the first row will be the x and the second row the y coordinates.
You can make all the arrays either allocatable (in which case the main program can do the allocation). If you want to make them of fixed size, make sure that simply changing n_charges will be enough to work with a different number of charges:
INTEGER, PARAMETER :: n_charges=2 ! Change this later REAL(KIND=wp), DIMENSION(n_charges) :: q ! This stays the same ...
Here is an outline of a possible module Charges:
MODULE Charges
...
PUBLIC
INTEGER :: n_charges ! The number of charges present
! In this program, q is of shape [n_charges],
! while r is of shape [2,n_charges]
REAL(KIND=wp), ALLOCATABLE, DIMENSION(:) :: q ! The magnitudes
REAL(KIND=wp), ALLOCATABLE, DIMENSION(:,:) :: r ! The positions
CONTAINS
! This function gives the potential at the point <x,y> by
! summing over the contributions from the individual charges:
FUNCTION V(x,y) RESULT(U)
REAL(KIND=wp), INTENT(IN) :: x,y
REAL(KIND=wp) :: U
...
END FUNCTION V
! This function gives the electric field at the point <x,y>:
FUNCTION E(x,y) RESULT(E_)
REAL(KIND=wp), INTENT(IN) :: x,y
REAL(KIND=wp), DIMENSION(2) :: E_
...
END FUNCTION E
...
END MODULE Charges
We have not given the body of the function V(x,y) and
here.
There are a lot of options. You know that you need to perform a sum over the
contributions from different charges (probably using a DO loop).
You can calculate the magnitude (distance) r appearing in
in different ways. For example, here is the distance from the observation
point (x,y) to the
charge:
dr = SQRT( (x-r(1,i))**2 + (y-r(2,i))**2 ) + TINY(0.0_wp)
In the main program, you will need to assign values to the magnitudes and positions of the charges. For the magnitudes, you can use element-by-element assignment (see below) or array-constructors. For the positions of the charges, you can first set all of them to zero, and then do an element-by-element assignment to the non-zero entries (see below). You can also use array constructors for the positions as well. However, array constructors always return rank-1 arrays. You can use the intrinsic function RESHAPE to make this a rank-2 array:
r = RESHAPE( SOURCE = (/ (/1.0,2.0/), (/-1.0,0.0/), &
(/-.5,-1.0/), (/2.0,-1.0/) /), &
SHAPE = (/2,n_charges/) )
Other approaches include letting the user enter these with READ (in column-wise order for rank-2 arrays) or assigning random numbers (for 10 charges this may be a choice):
CALL RANDOM_NUMBER(r) ! Intrinsic function for random numbers
For example, one possible strategy is:
PROGRAM Electrostatics USE Charges ... n_charges=2 ALLOCATE(q(n_charges),r(2,n_charges)) ! A dipole: q(1)=1.0_wp q(2)=-1.0_wp ! Positions: r=0.0_wp r(1,1)=-1.0_wp r(1,2)=1.0_wp ! What are the positions of the charges? Picture? ... CALL FunPlot3D(f_xy=V,...) ... END PROGRAM Electrostatics
After these values have been assigned, the main program needs to call the
plotting routines from the FunGraphics module and pass the
functions V and 
from the module Charges as arguments.
The examples in the documentation are sufficient help for this.
Please ask your TA or instructor if this is unclear. The remaining two worksheets will follow the exact same strategy, only with different physics. Understanding this one one will save us all time and frustrations, even if you need an extra week to finish it.
This document was generated using the LaTeX2HTML translator Version 98.1p1 release (March 2nd, 1998)
Copyright © 1993, 1994, 1995, 1996, 1997, Nikos Drakos, Computer Based Learning Unit, University of Leeds.
The command line arguments were:
latex2html -split 0 worksheet7extra_f00.tex.
The translation was initiated by Phil Duxbury on 2000-10-24