      PROGRAM ELECPOT
C
C     ELECTRICAL POTENTIAL FOR TWO POINT CHARGES
C
C     -----------------------------------------------
C     AUTHOR: ROBERT EHRLICH (GEORGE MASON U    1973)
C     UPDATE: PETER SIGNEL (MSU PHYSICS DEPT 3/01/86)
C     UPDATE: GARY HERZENSTIEL (MSU CPS DEPT 3/25/86)
C     Updated  November 1987  E.J.D.Kales
C
C     calculate the electrical potential for 2 point charges in an
C     array (21x21) in the range X= -10 to +10 and Y= -10 to +10
C
C     SUBROUTINES: GETCORD - get coordinates for the 2 charges
C                  CALCEP  - calculate electrical potential matrix
C                  DISPLAY - write the matrix to output
C
C                        E L E C P O T
C                              :
C            ---------------------------------------
C            :             :            :          :
C         GETCORD        CALCEP       DISPLAY     QUERY
C      (A1,A2,Q1,Q2) (A1,A2,Q1,Q2,V)    (V)      (AGAIN)
C
C     ***** DECLARATIONS *****
      REAL V(-10:10,-10:10),A1,A2,Q1,Q2
      CHARACTER AGAIN
      OPEN(UNIT=3,FILE='m348p1f.out')
C                    
C     ***** M A I N   B O D Y *****
C
 10   WRITE(3,'(A1)') 12
      CALL GETCORD(A1,A2,Q1,Q2     )
      CALL CALCEP (A1,A2,Q1,Q2,   V)
      CALL DISPLAY(V               )
      CALL QUERY  (           AGAIN)
C     IF USER IS DONE, END
      IF(AGAIN.EQ.'Y') GOTO 10
      STOP
      END
C
C     ***** G E T C O R D *****
C
      SUBROUTINE GETCORD(A1,A2,Q1,Q2)
C
C     ***** DECLARATIONS *****
      REAL A1,A2,Q1,Q2
C
      PRINT *,'Please supply A1, A2, Q1, & Q2 values:'
      WRITE(*,20) 'A1='
      READ *,A1
      WRITE(*,20) 'A2='
      READ *,A2
      WRITE(*,20) 'Q1='
      READ *,Q1
      WRITE(*,20) 'Q2='
      READ *,Q2      
  20  FORMAT(A3/)
      PRINT *
      WRITE(*,34) A1,A2,Q1,Q2
      WRITE(3,34) A1,A2,Q1,Q2
 34   FORMAT('m348p1f.for:    Potential for 2 Point-Charges:'/
     +' A1=',F8.5,'  A2=',F8.5,'  Q1=',F8.5,'  Q2=',F8.5)
      RETURN
      END
C
C     ***** C A L C E P *****
C
      SUBROUTINE CALCEP(A1,A2,Q1,Q2,V    )
C
C     calculate electric potential matrix
C
C     ***** DECLARATIONS *****
      INTEGER X,Y
      REAL V(-10:10,-10:10),A1,A2,Q1,Q2
C
C     V    array holding electrical potentials
C
      DO 40,X=-10,10
          DO 40,Y=-10,10
C
C     calculate distances R1 & R2 from (X,Y) to each of 
C     the 2 charges. increment by +.000001 to insure R1 & R2 
C	  are not zero.
C
              R1=SQRT((X-A1)**2+Y**2)+.000001
              R2=SQRT((X-A2)**2+Y**2)+.000001
C
C             calculate net electrical potential for the 2 
C             point charges, store the potentials in the matrix.
C
              V(X,Y)=Q1/R1+Q2/R2
C
 40   CONTINUE
      RETURN
      END
C
C     ***** D I S P L A Y *****
C
      SUBROUTINE DISPLAY(V     )
C
C     ***** DECLARATIONS *****
      REAL V(-10:10,-10:10)
      INTEGER X,Y
C
C     display matrix V with proper printer carriage width
C
CC for a wide printer:
CC        WRITE(3,55) ((V(X,Y),X=-10,10),Y=-10,10)
CC 55     FORMAT('132 Width'/(21(F4.2,1X)//))
CC for a small printer:
      WRITE(3,55) ((V(X,Y),X=-10,5),Y=-10,10)
 55   FORMAT('80 Width'/(16(F4.2,1X)//))
      RETURN
      END
C
C     ***** Q U E R Y *****
C
      SUBROUTINE QUERY(    AGAIN)
C
C     ***** DECLARATIONS *****
      CHARACTER AGAIN
C
      WRITE(*,62)
 62   FORMAT(
     +'Would you like to do this again using new values? (Y/N)? '/)
      READ( *,'(A1)')AGAIN
      PRINT *
      RETURN
      END
C     ***** P R O G R A M   E N D *****
