!___________________________________________________________________________________________________
! Aleksandar Donev
! PHY201 Solutions to Worksheet
! December, 2000, MSU
!___________________________________________________________________________________________________
!
!______________________________________________________________________________
! 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 <x,y> 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="PNG", 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="PNG", 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="PNG", 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
!