! ! 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 . ! -----------------------------------------------------------------------