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