!
! Use Internal procedures in a program: https://www.tutorialspoint.com/fortran/fortran_procedures.htm , https://en.wikibooks.org/wiki/Fortran/Fortran_procedures_and_functions
! All variables are accessible to internal procedures by "host association"
! Internal procedures are functions or subroutines that follow a CONTAINS statement in a program unit. 
! The program unit in which the internal procedure appears is called its host.
! Internal procedures can appear in the main program, in an external subprogram, or in a module subprogram.
!   . Only the host program unit can use an internal procedure.
!   . An internal procedure has access to host entities by host association; that is, names declared in the host program unit are useable within the internal procedure.
!   . An internal procedure must not contain an ENTRY statement.
!
! Compile with:
! gfortran ELines.f90 -o el
! Example Usage: 
! (echo 2;echo 1 0  1;echo -1 0 -1                       )| ./el | gnuplot -p -d -e 'r=5  ;plot [-r:r][-r:r] "-" w dots not'
! (echo 3;echo 1 0  1;echo -1 0 -1;echo 0 1 1            )| ./el | gnuplot -p -d -e 'r=2  ;plot [-r:r][-r:r] "-" w dots not'
! (echo 4;echo 1 0 -1;echo -1 0 -1;echo 0 1 1;echo 0 -1 1)| ./el | gnuplot -p -d -e 'r=2.5;plot [-r:r][-r:r] "-" w dots not'
! Exercise: modify ndraw to accommodate charges that have different magnitude.
!****************************************************
program      Electric_Fields
!****************************************************
 implicit none
 real   , allocatable, dimension(:) :: X,Y,Q
 integer                            :: N
 integer                            :: i , j , ndraw
 real                               :: x0, y0, R, theta
 real   , parameter                 :: PI       = 3.141592
 real   , parameter                 :: step     = 0.01 
 real   , parameter                 :: max_dist = 20.0
!----------------------------------------------------
!Set charge distribution:
 print  *, '# Enter number of charges: '
 read   *, N
 print  *, '# N = ',N

 ALLOCATE(X(N),Y(N),Q(N))

 do i = 1, N
  print *, '# Charge: ', i
  print *, '# Enter position and charge (X,Y,Q): '
  read  *,  X(i), Y(i), Q(i)
  print *, '# (X,Y)= ', X(i), Y(i), ' Q= ', Q(i)
 end do

!----------------------------------------------------
!Drawing lines:
 ndraw   = 12 
 R       = 0.1

 do  i   = 1, N
  do j   = 1, ndraw
   theta = (2.0 * PI / ndraw) * j
   x0    = X(i) + R  * cos(theta)
   y0    = Y(i) + R  * sin(theta)
!----------------------------------------------------
   call eline(x0, y0)
!----------------------------------------------------
  end do
 end do


!****************************************************
!----------------------------------------------------
contains
!----------------------------------------------------
!****************************************************


!****************************************************
 subroutine     eline(xin,yin)
!****************************************************
!No need for implicit none: taken from main program
  real             :: xin,yin,x0,y0        !x0,y0 shadow globals
  integer          :: i,direction          !i is local
  real             :: rmin,rmax,r,dx,dy,dl !r is local
  real             :: Ex,Ey,E
!----------------------------------------------------
  do direction = -1,1,2                    !direction= +/- 1
   dl = direction * step
   x0 = xin; y0 = yin
   dx = 0.0; dy = 0.0
!----------------------------------------------------
   call mdist(x0,y0,rmin,rmax)
!----------------------------------------------------
   do while(rmin > (2.0*step) .and. rmax < max_dist)
    print *,x0,y0
!----------------------------------------------------
!   We evaluate the E-field at the midpoint: This reduces
!   systematic errors 
!----------------------------------------------------
    call efield(x0+0.5*dx,y0+0.5*dy,Ex,Ey)
!----------------------------------------------------
    E  = sqrt(Ex*Ex+Ey*Ey)
    if( E .le. 1.0e-10 ) exit
    dx = dl*Ex/E ; dy = dl*Ey/E
    x0 = x0 + dx ; y0 = y0 + dy
!----------------------------------------------------
    call mdist(x0,y0,rmin,rmax)
!----------------------------------------------------
   end do                    !do while()

  end do                     !do direction = -1,1,2
 end subroutine eline
!****************************************************
 subroutine     efield(x0,y0,Ex,Ey)
!****************************************************
  real              :: x0,y0,dx,dy,Ex,Ey
  integer           :: i
  real              :: r3,xi,yi
!----------------------------------------------------
  Ex  = 0.0     ; Ey = 0.0
  do i= 1,N
   xi = x0-X(i) ; yi = y0-Y(i)
!  Exercise: Improve code so that xi*xi+yi*yi=0 is taken care of
   r3 = (xi*xi+yi*yi)**(-1.5)
   Ex = Ex + Q(i)*xi*r3
   Ey = Ey + Q(i)*yi*r3
  end do
 end subroutine efield
!****************************************************
 subroutine mdist(x0,y0,rmin,rmax)
!****************************************************
  implicit none
  real              :: x0,y0,rmin,rmax
  integer           :: i
  real              :: r
!----------------------------------------------------
  rmax = -HUGE(rmax) ; rmin = HUGE(rmin)
  do i = 1,N
   r   = sqrt((x0-X(i))**2 + (y0-Y(i))**2)
   if(r > rmax) rmax = r
   if(r < rmin) rmin = r
  end do
 end subroutine mdist
!****************************************************
!----------------------------------------------------
end program Electric_Fields
!----------------------------------------------------
!****************************************************
!  ---------------------------------------------------------------------
!  Copyright by Konstantinos N. Anagnostopoulos (2004-2021)
!  Physics Dept., National Technical University,
!  konstant@mail.ntua.gr, www.physics.ntua.gr/~konstant
!  
!  This program is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation, version 3 of the License.
!  
!  This program is distributed in the hope that it will be useful, but
!  WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
!  General Public License for more details.
!  
!  You should have received a copy of the GNU General Public Liense along
!  with this program.  If not, see <http://www.gnu.org/licenses/>.
!  -----------------------------------------------------------------------
