PROGRAM ELCWRE
C
C     -------------------------------------
C     AUTHOR: ROBERT EHRLICH (GEORGE MASON U.   1973)
C     Revised: Peter Signell (Mich. St. U.      3/86)
C     Revised: Gary Herzenstiel (Mich. St. U.   3/86)
C     Revised: E.J.D.Kales (Mich. St. U.       11/87)
C
C     Calculates the electric potential for charged wire
C     assumed to be a continous finite line of charge
C
C     SUBROUTINES:
C                  GETINF  - get input values of N,A,X0,Y0
C                  CALC    - calculates V1,V2,V3,V4
C                  DISPLA  - write to file V1,V2,V3,V4
C                  QUERY   - ask if user wants to do again
C
C                     E L C W R E
C                          :
C         :          :           :           :
C      GETINF      CALC        DISPLA      QUERY
C   (N,A,X0,Y0) (N,A,X0,Y0, (V1,V2,V3,V4) (AGAIN)
C               V1,V2,V3,V4)    
C
C     ****** DECLARATIONS ******
      REAL N,A,X0,Y0,V1,V2,V3,V4
      CHARACTER AGAIN
      OPEN(UNIT=3,FILE='m349p1f.out')
      WRITE(3,'(A1//)') 12
C
C     ***** M A I N   B O D Y *****
C
 10   CALL GETINF (               N,A,X0,Y0)
      CALL CALC   (N,A,X0,Y0,   V1,V2,V3,V4)
      CALL DISPLA (V1,V2,V3,V4             )
      CALL QUERY  (                   AGAIN)
      IF(AGAIN.EQ.'Y') GOTO 10
      STOP
      END
C
C     ***** G E T I N F *****
C
      SUBROUTINE GETINF(   N,A,X0,Y0)
C 
C     get input values for N,A,X0,Y0
C
C     ***** DECLARATIONS *****
      REAL N,A,X0,Y0
C
C     N = number of point charges (number of terms used in 
C         approximation)
C     A = half-length of the line from -A to +A along the
C         X-axis
C     X0,Y0 = X and Y of point at which potential is to be
C         calculated
C
      PRINT *,'Please supply N, A, X0, & Y0 values:'
      WRITE(*,20) ' N='
      READ *,N
      WRITE(*,20) ' A='
      READ *,A
      WRITE(*,20) 'X0='
      READ *,X0
      WRITE(*,20) 'Y0='
      READ *,Y0     
 20   FORMAT(A3/)
      PRINT *
      WRITE(*,24) N,A,X0,Y0
      WRITE(3,24) N,A,X0,Y0
 24   FORMAT(' m349p1f.for:  Electrical Potential of a Wire'//
     +'   N=',F10.5,'   A=',F10.5,'  X0=',F10.5,
     +'  Y0=',F10.5/)
      PRINT *
      RETURN
      END
C
C     ***** C A L C *****
C
      SUBROUTINE CALC(N,A,X0,Y0,   V1,V2,V3,V4)
C 
C     calculates V1,V2,V3,V4 from values of N,A,X0,Y0
C
C     V1 = exact potentials for finite line of (continuous)
C          charge
C     V2 = potential from summing contributions of N point
C          charges
C     V3 = Trapezoidal Rule approximation to the line of
C          charge
C     V4 = Simpson's Rule approximation to the line of
C          charge
C
C     ***** DECLARATIONS *****
      REAL N,A,X0,Y0,V1,V2,V3,V4,DX,X,C,L,F
      INTEGER J,M
C
      V1=ALOG((SQRT((X0-A)**2+Y0**2)-(X0-A))/
     +        (SQRT((X0+A)**2+Y0**2)-(X0+A)))
C
C     First redefine N to be the nearest larger odd integer-
C     Necessary for the application of Simpson's Rule
C
      M=N
      M=2*(M/2)+1
      N=M
C
C     V2,V3,V4 are computed as sums over number of charges
C     making up the line.  Set to zero initially.
C
      V2=0
      V3=0
      V4=0
C
C     DX spacing between adjacent point charges
C
      DX=2.0*A/(N-1.0)
C
C     start at one end of the line:
C
      X=-A
C
C     Add the potential due to each of the N charges:
C
      DO 40 J=1,M
C
C          L (Lambda) is charge density at position X;
C              initialize to 1.0
C
           L=1.0
C
C          F is the function to be integrated
C
           F=L/SQRT((X-X0)**2+Y0**2)
C
C          find C'S for Finite Sum
C
           C=1.0
           V2=V2+C*F
C
C          find C'S for Trapezoidal Rule
C
           C=1.0
           IF(J.EQ.1.OR.J.EQ.M) C=0.5
           V3=V3+C*F
C
C          find C'S for Simpson's Rule
C
           C=2
           JJ=2*(J/2)
           IF (J.EQ.JJ)C=4
           IF(J.EQ.1.OR.J.EQ.M) C=1.0
           V4=V4+C*F
C
           X=X+DX
C                 increment X by DX
 40   CONTINUE
      V2=V2*DX
      V3=V3*DX
      V4=V4*DX/3
      RETURN
      END
C
C     ***** D I S P L A *****
C
      SUBROUTINE DISPLA(V1,V2,V3,V4     )
C 
C     SHOWS THE VALUES OF V1,V2,V3,V4
C
C     ***** DECLARATIONS *****
      REAL V1,V2,V3,V4
C     
      WRITE(3,50)V1,V2,V3,V4
      WRITE(*,50)V1,V2,V3,V4
 50   FORMAT('  V1=',F10.5/'  V2=',F10.5/'  V3=',
     +F10.5/'  V4=',F10.5///)
      RETURN
      END 
C
C     ***** Q U E R Y *****
C
      SUBROUTINE QUERY(    AGAIN)
C 
C     ***** DECLARATIONS *****
      CHARACTER AGAIN
C
      PRINT *
      WRITE(*,62)
  62  FORMAT('Would you like to run the program again? [Y/N]
     +')
      READ(*,'(A1)')AGAIN
      PRINT *
      RETURN
      END
C     ***** P R O G R A M    E N D *****