Solution to worksheet 7
phy201
!______________________________________________________________________________
! This program plots the potential of a system of 4 charges in a plane
!______________________________________________________________________________
MODULE Charges
IMPLICIT NONE
PUBLIC
INTEGER, PARAMETER :: wp=KIND(0.0E0) ! All variables that are plotted must be single precision!
INTEGER :: n_charges ! The number of charges present
REAL(KIND=wp), ALLOCATABLE, DIMENSION(:) :: q ! The positions
REAL(KIND=wp), ALLOCATABLE, DIMENSION(:,:) :: r ! The magnitudes of the charges
! In this program, q is of shape [n_charges], while r is of shape [2,n_charges]
CONTAINS
! This function gives the potential at the point by summing over the potential
! of the individual charges. Notice the TINY(0.0_wp) added for numerical stability:
FUNCTION E(x,y) RESULT(E_)
REAL(KIND=wp), INTENT(IN) :: x,y
REAL(KIND=wp), DIMENSION(2) :: E_
INTEGER :: i
REAL(KIND=wp) :: dx,dy,dr
E_=0.0
DO i=1,n_charges
dx=x-r(1,i)
dy=y-r(2,i)
dr=SQRT(dx**2+dy**2)+EPSILON(0.0_wp)
E_(1)=E_(1)+q(i)*dx/dr**3
E_(2)=E_(2)+q(i)*dy/dr**3
END DO
END FUNCTION E
FUNCTION V(x,y) RESULT(U)
REAL(KIND=wp), INTENT(IN) :: x,y
REAL(KIND=wp) :: U
INTEGER :: i
REAL(KIND=wp) :: dx,dy,dr
U=0.0
DO i=1,n_charges
dx=x-r(1,i)
dy=y-r(2,i)
dr=SQRT(dx**2+dy**2)+EPSILON(0.0_wp)
U=U+q(i)/dr
END DO
END FUNCTION V
END MODULE Charges
PROGRAM Plot_Field
USE SimpleGraphics
USE FunGraphics
USE Charges
IMPLICIT NONE
REAL :: length, pi
CHARACTER(LEN=100), DIMENSION(3) :: title
pi=4*ATAN(1.0)
! This is an example of how an interactive input of the charges and magnitudes could be done:
! WRITE(*,*) "Enter the number of charges :"
! READ(*,*) n_charges
n_charges=4
ALLOCATE(q(n_charges),r(2,n_charges))
! WRITE(*,*) "Enter the magnitudes of the charges in the form q1,q2... :"
! READ(*,*) q
! WRITE(*,*) "Enter the coordinates of the charges in the form x1,y1,x2,y2... :"
! READ(*,*) r
q=(/1.0,-2.0,.5,-1.0/) ! The magnitudes
r=RESHAPE(SOURCE=(/1.0,2.0,-1.0,0.0,-.5,-1.0,2.0,-1.0/),SHAPE=(/2,n_charges/)) ! The charges
length=3.0 ! We will plot in the interval [-length,length] in both x and y
WRITE(title(2),"(A,10F6.2)") "q=",q
WRITE(title(3),"(A,20F6.2)") "r=",r
title(1)="The potential of a charge system (surface plot):"
CALL InitGraphics(file="potential3D.png",file_type="CONS", &
plot_title=title, &
x_label="x",y_label="y",z_label="V(x,y)")
CALL FunPlot3D(f_xy=V,xy_range=(/-length,length,-length,length/),plot_spec="3CR", &
axis=(/-length,length,-length,length,-3.0,3.0/),num_points=(/50,50/))
CALL EndGraphics()
title(1)="The potential of a charge system (contour plot):"
CALL InitGraphics(file="potential2D.png",file_type="CONS", &
plot_title=title, &
x_label="x",y_label="y",z_label="V(x,y)")
CALL FunPlotContour(f_xy=V,xyz_range=(/-length,length,-length,length,-3.0,3.0/), &
num_points=(/50,50,20/))
CALL EndGraphics()
title(1)="The electrostatic field of a charge system (vector field plot):"
CALL InitGraphics(file="field2D.png",file_type="CONS", &
plot_title=title, &
x_label="x",y_label="y",z_label="|E(x,y)|")
CALL FunVectorPlot2D(f_xy=E,xy_range=(/-length,length,-length,length/), &
num_points=(/25,25/),plot_spec="11g", &
axis=(/-length,length,-length,length,0.1,2.5/),zoom=0.3)
CALL EndGraphics()
END PROGRAM Plot_Field