PROGRAM MAIN
C     A PROGRAM TO ITERATE THE LAPLACIAN NUMERICAL SOLUTION TO THE FIELD
C     POTENTIAL.  THE POTENTIAL IS DEFINED BY DEFPOT AND THE CONDITION
C     FOR ANOTHER ITERATION IS DEFINED BY CHGPOT.
C
C     ORIGINALLY WRITTEN 1993 WHILE A STUDENT AT UW-RF.
C     COPYRIGHT 1993, 2007 JIM HALL, JHALL@FREEDOS.ORG
C Permission is hereby granted, free of charge, to any person obtaining a copy
C of this software and associated documentation files (the "Software"), to deal
C in the Software without restriction, including without limitation the rights
C to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
C copies of the Software, and to permit persons to whom the Software is
C furnished to do so, subject to the following conditions:
C 
C The above copyright notice and this permission notice shall be included in
C all copies or substantial portions of the Software.
C 
C THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
C IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
C FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
C AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
C LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
C OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
C THE SOFTWARE.


      PARAMETER (EPS=.001, NSIZE=15)
      REAL P(NSIZE,NSIZE), PNEW(NSIZE,NSIZE)

C     IDENTIFY THE PROGRAM AND GIVE PARAMETERS

      PRINT *, 'LAPLACIAN NUMERICAL SOLUTION TO THE FIELD POTENTIAL'
      PRINT *, 'ITERATE UNTIL CHANGE FALLS BELOW: ', EPS
      PRINT *, 'SIZE OF FIELD POTENTIAL:          ', NSIZE
      CALL DEFPOT (P, PNEW, NSIZE)
      PRINT *, 'INITIAL VALUES:'
      DO 10 J = 1, NSIZE
         PRINT 1000, (P(I,J), I = 1, NSIZE)
 1000    FORMAT (X, 40(X, F6.2))
 10   CONTINUE

C     ITERATE UNTIL THE CHANGE FALLS BELOW EPS

      PRINT *, 'WORKING..'
      N = 0
 20   N = N + 1
      DO 40 I = 1, NSIZE
         DO 30 J = 1, NSIZE
            IF ((I.EQ.1).AND.(J.GT.1).AND.(J.LT.NSIZE)) THEN
               PNEW(I,J) = (P(I+1,J) + P(I,J-1) + P(I,J+1)) / 3.
            ELSE IF ((I.EQ.NSIZE).AND.(J.GT.1).AND.(J.LT.NSIZE)) THEN
               PNEW(I,J) = (P(I-1,J) + P(I,J-1) + P(I,J+1)) / 3.
            ELSE IF ((J.EQ.1).AND.(I.GT.1).AND.(I.LT.NSIZE)) THEN
               PNEW(I,J) = (P(I-1,J) + P(I+1,J) + P(I,J+1)) / 3.
            ELSE IF ((J.EQ.NSIZE).AND.(I.GT.1).AND.(I.LT.NSIZE)) THEN
               PNEW(I,J) = (P(I-1,J) + P(I+1,J) + P(I,J-1)) / 3.
            ELSE IF ((I.EQ.1).AND.(J.EQ.1)) THEN
               PNEW(I,J) = (P(I+1,J) + P(I,J+1)) / 2.
            ELSE IF ((I.EQ.1).AND.(J.EQ.NSIZE)) THEN
               PNEW(I,J) = (P(I+1,J) + P(I,J-1)) / 2.
            ELSE IF ((I.EQ.NSIZE).AND.(J.EQ.1)) THEN
               PNEW(I,J) = (P(I-1,J) + P(I,J+1)) / 2.
            ELSE IF ((I.EQ.NSIZE).AND.(J.EQ.NSIZE)) THEN
               PNEW(I,J) = (P(I-1,J) + P(I,J-1)) / 2.
            ELSE
               PNEW(I,J) = (P(I-1,J) + P(I+1,J) + P(I,J-1) + P(I,J+1)) /
     $              4.
            ENDIF
 30      CONTINUE
 40   CONTINUE

      IF (CHGPOT (P, PNEW, NSIZE).GT.EPS) THEN
         CALL DEFPOT (P, PNEW, NSIZE)
         GOTO 20
      ENDIF

C     PRINT THE FINAL VALUES

      PRINT *, 'AFTER ', N, ' ITERATIONS, THE SOLUTION CONVERGED'
      PRINT *, 'FINAL VALUES:'
      DO 50 J = 1, NSIZE
         PRINT 1000, (P(I,J), I = 1, NSIZE)
 50   CONTINUE

      PRINT *, 'IN DATA SUITABLE FOR GNUPLOT:'
      DO 70 J = 1, NSIZE
         DO 60 I = 1, NSIZE
            PRINT *, I, J, P(I,J)
 60      CONTINUE
         PRINT *
 70   CONTINUE

      PRINT *, 'DONE'
      STOP
      END

CCCCCCCCCCCCCC

      SUBROUTINE DEFPOT (P, PNEW, NSIZE)
      REAL P(NSIZE,NSIZE), PNEW(NSIZE,NSIZE)
C     DEFINES THE FIELD POTENTIAL P, AND COPIES THE NEW POTENTIAL PNEW
C     INTO P
C
C     INPUT:
C     P - THE ARRAY OF THE FIELD POTENTIAL, BEFORE ITERATION
C
C     PNEW - THE ARRAY OF THE FIELD POTENTIAL, AFTER ITERATION
C
C     NSIZE - THE SIZE OF THE FIELD POTENTIAL, BOTH P AND PNEW
C
C     OUTPUT:
C     P - THE ARRAY OF THE FIELD POTENTIAL, ALL READY FOR A NEW
C     ITERATION.

      DO 20 I = 1, NSIZE
         DO 10 J = 1, NSIZE
            P(I,J) = PNEW(I,J)
 10      CONTINUE
 20   CONTINUE

C     DEFINE THE FIELD POTENTIAL---THIS MUST BE CHANGED WITH THE
C     APPLICATION OF THE PROGRAM.

      DO 30 I = 1, NSIZE
         P(I, 5) = 10.
 30   CONTINUE
      RETURN
      END

CCCCCCCCCCCCCCCCCC

      REAL FUNCTION CHGPOT (P, PNEW, NSIZE)
      REAL P(NSIZE,NSIZE), PNEW(NSIZE,NSIZE)
C     FINDS THE LARGEST CHANGE IN THE POTENTIAL AND RETURNS IT.
C     THIS MUST BE CHANGED TO
C     REFLECT THE PARTICULAR APPLICATION, TO IGNORE PRE-DEFINED
C     POTENTIALS FROM DEFPOT.
C
C     INPUT:
C     P - THE FIELD POTENTIAL ARRAY, BEFORE ITERATION
C
C     PNEW - THE FIELD POTENTIAL ARRAY, AFTER ITERATION
C
C     NSIZE - THE SIZE OF THE FIELD POTENTIAL ARRAY
C
C     EPS - THE LIMIT OF THE PERCENTAGE CHANGE FOR AN ITERATION

      DELTA = 0.

      DO 20 I = 1, NSIZE
         DO 10 J = 1, NSIZE
            IF (J.NE.5) THEN
               DIFF = ABS(PNEW(I,J) - P(I,J))
               DELTA = MAX(DIFF, DELTA)
            ENDIF
 10      CONTINUE
 20   CONTINUE

      CHGPOT = DELTA
      RETURN
      END