next up previous


Electrostatics with FORTRAN and Mathematica

Aleksandar Donev - Dr. Phillip Duxbury

Due Friday 27th Oct. - two weeks

Physics 201 home

We will now start a series of worksheets dealing with electrostatic and electrodynamics problems. You will solve each problem using both Fortran 90 and Mathematica. Mathematica is very often useful in order to check the validity of results obtained via alternative methods, but for complex problems it is often inefficient.

In this first application to EM problems, you will find and plot the electrostatic potential and electric field of a set of static charges. You should write the code for an arbitrary number of charges whose positions and charges are input by the user. You will plot the potential and electric field in a given plane (input by the user). Illustrate your code (both Mathematica and FORTRAN) by plotting the potential and field near: (i) a dipole; (ii) a quadrapole; (iii) a configuration of 10 charges of your chosing.

Electrostatic potential and field of a system of charges

The potential of a single charge is,


 \begin{displaymath}V(\mathbf{r})=\frac{Kq}{r} \end{displaymath} (1)

while the electric field is given by:


 \begin{displaymath}\mathbf{E}(\mathbf{r})= - {\bf\nabla} V = \frac{Kq\mathbf{r}}{r^{3}} \end{displaymath} (2)

Here bold letters denote a vector, while regular letters are scalars. For this worksheet, just take K=1. Also remember that the potential and field of a group of charges is the sum of the individual potentials (scalar sum) and fields(vector sum).

Solution in Mathematica

Store the magnitudes of the charges in an array (lists), and store their positions in a two-dimensional array (nested list). For example for two charges: 

   q={1,-2};
   rq={{1,2,0},{0,-1,0}};

Now define a Mathematica function (operator) to evaluate the potential at a given point (vector) $\mathbf{r}=\left( x,y,z\right) $. Try to devise an operator that will work with any q, rq and r.

Just to remind you, here is an example of defining an operator in Mathematica that calculates the square of an expression x: 

   sqr[x_]:=x^2

You may also need a function to calculate the length of a vector r, $r=\left\Vert \mathbf{r}\right\Vert $. You can use, for example, a dot product ( .) to do this: 

   magnitude[r_]:= Sqrt[r.r]

Plot the electric potential both as a surface (Plot3D in Mathematica) and as a contour plot (ContourPlot) and compare the two. Which one is better and why? Try to make both plots look as nice as possible. Then plot the electric field as a vector field ( PlotVectorField loaded with the package via << Graphics`PlotField`).

Solution in FORTRAN

Plotting the potential in Fortran involves more work then doing it in Mathematica. This is your first attempt to write a complete FORTRAN application, so you will probably find it quite a challenge to design and write this code on your own. You will need to use simple arrays extensively, so review this material (sections 1.9.1, 1.9.2, 2.6.1, 2.6.2, 2.6.3, 2.6.4 and 2.6.9 in the manual. In particular, try to understand the concept of array constructors, as these can make your life much easier for this assignment). The positions of the charges will be stored in a two-dimensional array (matrix), and this is the first time you encounter these in this class. They are no different from vectors, at least at this level of usage. We have prepared some additional help in getting started with the FORTRAN code, but you should try to do it without this help first.

Methodology

The first thing to do is to make functions which calculate the potential and the electric field for a given array of charges.

In particular, you should make a new module and place a function in it that gives the electrostatic potential at a certain point $% \left( x,y\right) $. The interface of this function (the name is irrelevant) should be: 

     FUNCTION V(x,y) RESULT(U)
       REAL(KIND=sp), INTENT(IN) :: x,y
       REAL(KIND=sp) :: U
       ...Calculate U=V(x,y) here...
       ...Sum the potential contributions from each charge...
     END FUNCTION V

It is up to you to devise a way of calculating the potential inside this function. The best and recommended approach is to make arrays that will hold the magnitudes and positions of the charges, just like we did in Mathematica, and then sum the individual contributions to the potential from each charge. To insure that you never divide by zero, add a small constant $% \varepsilon $, called TINY(0.0_sp) in Fortran 90, to the denominator in equations [*] and [*].

This function V will need to know the positions and magnitudes of the charges. Since these can not be passed to the procedure (because the plotting routine described below only accepts functions with the above interface), you should make these global (module) variables in the module where you put the function V(x,y). That way the main program can work with these arrays as well.

We will do the plotting in Fortran similarly to the way we did it in Mathematica. However, now you can pass the function V as an argument to the plotting routines described below, which will calculate this function at a grid of points and invoke the DISLIN graphics library to render the plot.

After this works and you get what you consider to be a good plot of the potential, move to calculating and plotting the electric field. You will need to make a function, for example E, that returns the two components of the electric field E(x,y) for given x and y. The result should in fact be an array of length 2, in the form (/Ex,Ey /): 

     FUNCTION E(x,y) RESULT(E_)
       REAL(KIND=sp), INTENT(IN) :: x,y
       REAL(KIND=sp), DIMENSION(2) :: E_
       ...Calculate E_=E(x,y) here...
     END FUNCTION V

There are many possibilities for coding equations [*] and [*] in Fortran, some faster than others, some more elegant than others. For those that like trying new things, we would suggest reading section 2.6.5 and 2.6.6. Using array operations in Fortran 90, it is possible to code the above formulas so they almost fully resemble the Mathematica coding, but since the number of dimensions is fixed, faster solutions exist.

The FunGraphics module

In the last worksheet you learned how to make a simple 2D plot with the SimpleGraphics module. In 3D, making all the required arrays manually and assigning their values is more cumbersome, so we made a new module with pre-compiled routines that will plot a function instead of an array, thus simplifying the usage significantly. This FunGraphics module is documented in:

http://computation.pa.msu.edu/phy201/FunGraphics.html

You should now get a feeling of how these routines work and be able to use the routines FunPlot3D and/or FunPlotContour to plot the potential V(x,y) and the routine FunVectorPlot2D to plot the electric field. The example given in the above wepage, FunPlot_3D.f90, is very helpful. The executable of the example is in the file /classes/phy201/FunPlot_3D.x, and the code is also there. You only need to change some of the numbers and pass the functions you just made as an argument, as in: 

   CALL FunPlot3D(f_xy=V,xy_range=(/x_min,x_max,y_min,y_max/),...)

The same goes for the contour and vector plot as well. Play with the plots until they at least partially resemble the Mathematica plots and compare the two. If you need more time to finish this worksheet, talk to us.

Miscellaneous

Notice that we do not give a lot of pseudocode here. This is because you should already have an understanding of how to make a module and encapsulate functions in it and then USE the module in a program. The example file erf_fun_plot.f90 in /classes/phy201 solves last week's assignment using the FunGraphics module and can be very helpful to look at. The executable solution to this worksheet is in the file charges.x. The procedure for compiling and linking the program is the same as in the last worksheet.

Notice how we have used all the things we learned so far (with the error function project), only now in a physics context. Tell us how you feel about these worksheets!

About this document ...

Electrostatics with FORTRAN and Mathematica

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 worksheet7_f00.tex.

The translation was initiated by Phil Duxbury on 2000-10-16

next up previous
Phil Duxbury
2000-10-16