C******************************************************************************
C                                                                             *
C                        FINITE ELEMENT METHOD                                *
C                                                                             *
C******************************************************************************
C
C     TO APPROXIMATE THE SOLUTION TO AN ELLIPTIC PARTIAL-DIFFERENTIAL
C     EQUATION SUBJECT TO DIRICHLET, MIXED, OR NEUMANN BOUNDARY
C     CONDITIONS:
C
C     INPUT   SEE STEP 0
C
C     OUTPUT  DESCRIPTION OF TRIANGLES, NODES, LINE INTEGRALS, BASIS
C             FUNCTIONS, LINEAR SYSTEM TO BE SOLVED, AND THE
C             COEFFICIENTS OF THE BASIS FUNCTIONS
C
C
C     STEP 0
C     GENERAL OUTLINE
C
C        1. TRIANGLES NUMBERED: 1 TO K FOR TRIANGLES WITH NO EDGES ON
C           SCRIPT-S-1 OR SCRIPT-S-2, K+1 TO N FOR TRIANGLES WITH
C           EDGES ON SCRIPT-S-2, N+1 TO M FOR REMAINING TRIANGLES.
C           NOTE: K=0 IMPLIES THAT NO TRIANGLE IS INTERIOR TO D.
C           NOTE: M=N IMPLIES THAT ALL TRIANGLES HAVE EDGES ON
C           SCRIPT-S-2.
C
C        2. NODES NUMBERED: 1 TO LN FOR INTERIOR NODES AND NODES ON
C           SCRIPT-S-2, LN+1 TO LM FOR NODES ON SCRIPT-S-1.
C           NOTE: LM AND LN REPRESENT LOWER CASE M AND N RESP.
C           NOTE: LN=LM IMPLIES THAT SCRIPT-S-1 CONTAINS NO NODES.
C           NOTE: IF A NODE IS ON BOTH SCRIPT-S-1 AND SCRIPT-S-2, THEN
C           IT SHOULD BE TREATED AS IF IT WERE ONLY ON SCRIPT-S-1.
C
C        3. NL=NUMBER OF LINE SEGMENTS ON SCRIPT-S-2
C           LINE(I,J) IS AN NL BY 2 ARRAY WHERE
C           LINE(I,1)= FIRST NODE ON LINE I AND
C           LINE(I,2)= SECOND NODE ON LINE I TAKEN
C           IN POSITIVE DIRECTION ALONG LINE I
C
C        4. FOR THE NODE LABELLED KK,KK=1,...,LM WE HAVE:
C           A) COORDINATES XX(KK),YY(KK)
C           B) NUMBER OF TRIANGLES IN WHICH KK IS A VERTEX= LL(KK)
C           C) II(KK,J) LABELS THE TRIANGLES KK IS IN AND
C              NV(KK,J) LABELS WHICH VERTEX NODE KK IS FOR
C              EACH J=1,...,LL(KK)
C
C        5. NTR IS AN M BY 3 ARRAY WHERE
C           NTR(I,1)=NODE NUMBER OF VERTEX 1 IN TRIANGLE I
C           NTR(I,2)=NODE NUMBER OF VERTEX 2 IN TRIANGLE I
C           NTR(I,3)=NODE NUMBER OF VERTEX 3 IN TRIANGLE I
C
C        6. FUNCTION SUBPROGRAMS:
C           A) P,Q,R,F,G,G1,G2 ARE THE FUNCTIONS SPECIFIED BY
C              THE PARTICULAR DIFFERENTIAL EQUATION
C           B) RR IS THE INTEGRAND
C              R*N(J)*N(K) ON TRIANGLE I IN STEP 4
C           C) FFF IS THE INTEGRAND F*N(J) ON TRIANGLE I IN STEP 4
C           D) GG1=G1*N(J)*N(K)
C              GG2=G2*N(J)
C              GG3=G2*N(K)
C              GG4=G1*N(J)*N(J)
C              GG5=G1*N(K)*N(K)
C              INTEGRANDS IN STEP 5
C           E) QQ(FF) COMPUTES THE DOUBLE INTEGRAL OF ANY
C              INTEGRAND FF OVER A TRIANGLE WITH VERTICES GIVEN BY
C              NODES J1,J2,J3 - THE METHOD IS AN O(H**2) APPROXIMATION
C              FOR TRIANGLES
C           F) SQ(PP) COMPUTES THE LINE INTEGRAL OF ANY INTEGRAND PP
C              ALONG THE LINE FROM (XX(J1),YY(J1)) TO (XX(J2),YY(J2))
C              BY USING A PARAMETERIZATION GIVEN BY:
C                X=XX(J1)+(XX(J2)-XX(J1))*T
C                Y=YY(J1)+(YY(J2)-YY(J1))*T
C              FOR 0 <= T <= 1
C              AND APPLYING SIMPSON'S COMPOSITE METHOD WITH H=.01
C
C        7. ARRAYS:
C           A) A,B,C ARE M BY 3 ARRAYS WHERE THE BASIS FUNCTION N
C              FOR THE ITH TRIANGLE, JTH VERTEX IS
C              N(X,Y)=A(I,J)+B(I,J)*X+C(I,J)*Y
C              FOR J=1,2,3 AND I=1,2,...,M
C           B) XX,YY ARE LM BY 1 ARRAYS TO HOLD COORDINATES OF NODES
C           C) LINE,LL,II,NV,NTR HAVE BEEN EXPLAINED ABOVE
C           D) GAMMA AND ALPHA ARE CLEAR
C
C        8. NOTE THAT A,B,C,XX,YY,I,I1,I2,J1,J2,J3,DELTA ARE LABELLED
C           COMMON STORAGE SO THAT IN ANY SUBPROGRAM WE KNOW THAT
C           TRIANGLE I HAS VERTICES (XX(J1),YY(J1)),(XX(J2),YY(J2)),
C           (XX(J3),YY(J3)). THAT IS, VERTEX 1 IS NODE J1, VERTEX 2 IS
C           NODE J2, VERTEX 3 IS NODE J3 UNLESS A LINE INTEGRAL IS
C           INVOLVED IN WHICH CASE THE LINE INTEGRAL GOES FROM NODE J1
C           TO NODE J2 IN TRIANGLE I OR UNLESS VERTEX I1 IS NODE J1
C           AND VERTEX I2 IS NODE J2 - THE BASIS FUNCTIONS INVOLVE
C           A(I,I1)+B(I,I1)*X+C(I,I1)**Y FOR VERTEX I1 IN TRIANGLE I
C           AND A(I,I2)+B(I,I2)*X+C(I,I2)*Y FOR VERTEX I2 IN TRIANGLE I
C           DELTA IS 1/2 THE AREA OF TRIANGLE I
C
C     TO CHANGE PROBLEMS:
C        1) CHANGE FUNCTION SUBPROGRAMS P,Q,R,F,G,G1,G2
C        2) CHANGE DATA INPUT FOR K,N,M,LN,LM,NL.
C        3) CHANGE DATA INPUT FOR XX,YY,LLL,II,NV.
C        4) CHANGE DATA INPUT FOR LINE.
C        5) CHANGE DIMENSION STATEMENTS TO READ :
C           DIMENSION A(M,3),B(M,3),C(M,3),XX(LM),YY(LM)
C                   THERE ARE 10 PLACES TO CHANGE.
C           DIMENSION LINE(NL,2),LL(LM),II(LM,MAX LL(LM)),
C                NV(LM,MAX LL(LM))--1 PLACE TO CHANGE
C           DIMENSION NTR(M,3),GAMMA(LM),ALPHA(LN,LN+1)
C                   THERE IS 1 PLACE TO CHANGE
C
      EXTERNAL P,Q,RR,FFF,GG1,GG2,GG3,GG4,GG5
      COMMON A(10,3),B(10,3),C(10,3),XX(11),YY(11),I,I1,I2,J1,J2,J3
      COMMON DELTA
      DIMENSION LINE(5,2),LL(11),II(11,5),NV(11,5)
      DIMENSION NTR(10,3),GAMMA(11),ALPHA(5,6)
      CHARACTER NAME*14,NAME1*14,AA*1
      INTEGER INP,OUP,FLAG
      LOGICAL OK
      OPEN(UNIT=5,FILE='CON',ACCESS='SEQUENTIAL')
      OPEN(UNIT=6,FILE='CON',ACCESS='SEQUENTIAL')
      WRITE(6,*) 'This is the Finite Element Method.'
      OK = .FALSE.
      WRITE(6,*) 'The input will be from a text file in the following'
      WRITE(6,*) ' form:'
      WRITE(6,*) 'K     N     M     n     m     nl'
      WRITE(6,*) 'Each of the above is an integer -'
      WRITE(6,*) 'separate with at least one blank.'
      WRITE(6,*) 'Follow with the input for each node in the form:'
      WRITE(6,*) 'x-coor., y-coord., number of triangles in which the'
      WRITE(6,*) ' node is a vertex.'
      WRITE(6,*) 'Continue with the triangle number and vertex number'
      WRITE(6,*) 'for each triangle in which the node is a vertex.'
      WRITE(6,*) 'Separate each entry with at least one blank.  After'
      WRITE(6,*) 'all nodes have been entered follow with information'
      WRITE(6,*) 'on the lines over which line integrals must be'
      WRITE(6,*) 'computed. The format'
      WRITE(6,*) 'of this data will be the node number of the starting'
      WRITE(6,*) 'node, followed by the node number of the ending node'
      WRITE(6,*) 'for each line, taken in the positive direction.'
      WRITE(6,*) 'There should be 2 * nl such entries, each an integer'
      WRITE(6,*) 'separated by a blank.  Have all'
      WRITE(6,*) 'the functions P,Q,R,F,G,G1,G2 been created and'
      WRITE(6,*) 'has the input file been created.  Answer Y or N.'
      WRITE(6,*) ' '
      READ(5,*) AA
      IF((AA .EQ. 'y') .OR. (AA .EQ. 'Y')) THEN
         OK = .TRUE.
         WRITE(6,*) 'Input the file name in the form - drive: name.ext'
         WRITE(6,*) 'for example:  A:DATA.DTA'
         WRITE(6,*) ' '
         READ(5,*) NAME
         INP = 4
         OPEN(UNIT=INP,FILE=NAME,ACCESS='SEQUENTIAL')
         READ(4,*) K,N,M,LN,LM,NL
C        INPUT COORDINATES OF EACH NODE, NUMBER OF TRIANGLES NODE
C        IS IN (LLL) , WHICH TRIANGLES THE NODE IS IN, AND WHICH
C        VERTEX OF EACH TRIANGLE THE NODE IS IN
         DO 2 KK=1,LM
            READ(4,*) XX(KK),YY(KK),LLL,(II(KK,J),NV(KK,J),J=1,LLL)
3           FORMAT(2E15.8,1X,I2,1X,5(I2,1X,I2))
            LL(KK)=LLL
C           COMPUTE ENTRIES FOR NTR
            DO 4 J=1,LLL
                N1=II(KK,J)
                N2=NV(KK,J)
4               NTR(N1,N2)=KK
2        CONTINUE
C        INPUT LINE INFORMATION
         IF(NL.GT.0) THEN
           READ(4,*) ((LINE(I,J),J=1,2),I=1,NL)
         END IF
         CLOSE(UNIT=4)
      ELSE
         WRITE(6,*) 'The program will end so that the input file'
         WRITE(6,*) 'can be created or so the functions'
         WRITE(6,*) 'can be created.'
         OK = .FALSE.
      ENDIF
      IF(.NOT.OK) GOTO 400
      WRITE(6,*) 'Select output destinations: '
      WRITE(6,*) '1. Screen '
      WRITE(6,*) '2. Text file '
      WRITE(6,*) 'Enter 1 or 2 '
      WRITE(6,*) ' '
      READ(5,*) FLAG
      IF ( FLAG .EQ. 2 ) THEN
         WRITE(6,*) 'Input the file name in the form - '
         WRITE(6,*) 'drive:name.ext'
         WRITE(6,*) 'with the name contained within quotes'
         WRITE(6,*) 'as example:   ''A:OUTPUT.DTA'' '
         WRITE(6,*) ' '
         READ(5,*) NAME1
         OUP = 3
         OPEN(UNIT=OUP,FILE=NAME1,STATUS='NEW')
      ELSE
         OUP = 6
      ENDIF
      WRITE(OUP,*) 'FINITE ELEMENT METHOD'
      K1=K+1
      N1=LN+1
C     OUTPUT INFORMATION ON TRIANGLES, VERTICES AND NODES
      WRITE(OUP,100)
100   FORMAT(1X,'TRIANGLES ARE AS FOLLOWS')
      WRITE(OUP,101)
101   FORMAT(1X,'TRIANGLE',20X,'VERTEX 1',20X,'VERTEX 2',20X,
     1'VERTEX 3',/)
      WRITE(OUP,102) (I,(NTR(I,J),J=1,3),I=1,M)
102   FORMAT(1X,I3,21X,'NODE ',I3,20X,'NODE ',I3,20X,'NODE ',I3)
      WRITE(OUP,103)
103   FORMAT(1X,'NODE',25X,'X-COORD.',25X,'Y-COORD.',/)
      WRITE(OUP,104) (KK,XX(KK),YY(KK),KK=1,LM)
104   FORMAT(1X,I3,22X,E15.8,18X,E15.8)
C     OUTPUT INFO ON LINES
      WRITE(OUP,109)
109   FORMAT(1X,'DESCRIPTION OF LINES FOR LINE INTEGRALS',//)
      WRITE(OUP,110) (L1,LINE(L1,1),LINE(L1,2),L1=1,NL)
110   FORMAT(1X,'LINE',I3,'  FROM NODE',I3,' TO NODE',I3)
C     STEP 1
      IF(LM.NE.LN) THEN
           DO 5 L=N1,LM
5          GAMMA(L)=G(XX(L),YY(L))
      END IF
C     STEP 2 - INITIALIZATION OF ALPHA AND NOTE THAT
C     ALPHA(I,LN+1)=BETA(I)
      DO 6 I=1,LN
            DO 6 J=1,N1
6           ALPHA(I,J)=0.0
C     STEPS 3,4 AND 6 THROUGH 12 ARE WITHIN THE NEXT LOOP
C     FOR EACH TRIANGLE I LET NODE J1 BE VERTEX 1,NODE J2 BE VERTEX 2,
C     AND NODE J3 BE VERTEX 3
C     STEP 3
      DO 7 I=1,M
           J1=NTR(I,1)
           J2=NTR(I,2)
           J3=NTR(I,3)
           DELTA=XX(J2)*YY(J3)-XX(J3)*YY(J2)-XX(J1)*(YY(J3)-YY(J2))
     *     +YY(J1)*(XX(J3)-XX(J2))
           A(I,1)=(XX(J2)*YY(J3)-YY(J2)*XX(J3))/DELTA
           B(I,1)=(YY(J2)-YY(J3))/DELTA
           C(I,1)=(XX(J3)-XX(J2))/DELTA
           A(I,2)=(XX(J3)*YY(J1)-YY(J3)*XX(J1))/DELTA
           B(I,2)=(YY(J3)-YY(J1))/DELTA
           C(I,2)=(XX(J1)-XX(J3))/DELTA
           A(I,3)=(XX(J1)*YY(J2)-YY(J1)*XX(J2))/DELTA
           B(I,3)=(YY(J1)-YY(J2))/DELTA
           C(I,3)=(XX(J2)-XX(J1))/DELTA
C          STEP 4
C          I1=J FOR STEP 4 AND I1=K FOR STEP 7
           DO 8 I1=1,3
C          STEP 8
                JJ1=NTR(I,I1)
C               I2=K FOR STEP 4 AND I2=J FOR STEP 9
                DO 9 I2=1,I1
C                    STEP 10 AND STEP 4
                     JJ2=NTR(I,I2)
C                    ZZ=Z(I1,I2)**I
                     ZZ=B(I,I1)*B(I,I2)*QQ(P)+C(I,I1)*C(I,I2)*QQ(Q)
     *                    -QQ(RR)
C                        STEP 11 AND 12
                     IF(JJ1.LE.LN) THEN
                         IF(JJ2.LE.LN) THEN
                              ALPHA(JJ1,JJ2)=ALPHA(JJ1,JJ2) + ZZ
                              IF(JJ1.NE.JJ2) ALPHA(JJ2,JJ1) =
     *                         ALPHA(JJ2,JJ1)+ZZ
                         ELSE
                              ALPHA(JJ1,N1)=ALPHA(JJ1,N1)-ZZ*GAMMA(JJ2)
                         END IF
                     ELSE
                         IF(JJ2.LE.LN) ALPHA(JJ2,N1)=ALPHA(JJ2,N1)
     *                       -ZZ*GAMMA(JJ1)
                     END IF
9               CONTINUE
           HH=-QQ(FFF)
           IF(JJ1.LE.LN) ALPHA(JJ1,N1)=ALPHA(JJ1,N1)+HH
8          CONTINUE
7     CONTINUE
C     OUTPUT BASIS FUNCTIONS
      WRITE(OUP,105)
105   FORMAT(1X,'BASIS FUNCTIONS ON EACH TRIANGLE')
      WRITE(OUP,106)
106   FORMAT(1X,'TRIANGLE',6X,'VERTEX',9X,'NODE',33X,'FUNCTION',/)
      DO 108 I=1,M
      DO 108 J=1,3
108        WRITE(OUP,107) I,J,NTR(I,J),A(I,J),B(I,J),C(I,J)
107        FORMAT(1X,I3,10X,I3,10X,I3,10X,E15.8,' + ',E15.8,' * X + '
     *     ,E15.8,' * Y')
C     STEP 5
C     FOR EACH LINE SEGMENT JI=1,...,NL AND FOR EACH TRIANGLE I,
C     I=K1,...,N WHICH MAY CONTAIN LINE JI. SEARCH ALL 3 VERTICES
C     FOR POSSIBLE CORRESPONDENCE
C     STEP 5 AND STEPS 13 THROUGH 19
      IF(NL.NE.0.AND.N.NE.K) THEN
           DO 11 JI=1,NL
                DO 12 I=K1,N
                     DO 13 I1=1,3
C                    I1=J IN STEP 5 AND I1=K IN STEP 14
C                    STEP 15
                     J1=NTR(I,I1)
                     IF(LINE(JI,1).EQ.J1) THEN
                         DO 14 I2=1,3
C                        I2=K IN STEP 5 AND I2=J IN STEP 16
C                        STEP 17
                         J2=NTR(I,I2)
                         IF(LINE(JI,2).EQ.J2) THEN
C                        WE HAVE CORRESPONDENCE OF VERTEX I1
C                        IN TRIANGLE I WITH NODE J1 AS START OF LINE JI
C                        AND VERTEX I2 WITH NODE J2 AS END OF LINE JI
C                            STEP 5
                             XJ=SQ(GG1)
                             XJ1=SQ(GG4)
                             XJ2=SQ(GG5)
                             XI1=SQ(GG2)
                             XI2=SQ(GG3)
C                            STEP 8 AND 19
                             IF(J1.LE.LN) THEN
                                    IF(J2.LE.LN) THEN
                                         ALPHA(J1,J1)=ALPHA(J1,J1)+XJ1
                                         ALPHA(J1,J2)=ALPHA(J1,J2)+XJ
                                         ALPHA(J2,J2)=ALPHA(J2,J2)+XJ2
                                         ALPHA(J2,J1)=ALPHA(J2,J1)+XJ
                                         ALPHA(J1,N1)=ALPHA(J1,N1)+XI1
                                         ALPHA(J2,N1)=ALPHA(J2,N1)+XI2
                                    ELSE
                                         ALPHA(J1,N1)=ALPHA(J1,N1)-
     *                                      XJ*GAMMA(J2)+XI1
                                         ALPHA(J1,J1)=ALPHA(J1,J1)+XJ1
                                    END IF
                             ELSE
                                    IF(J2.LE.LN) THEN
                                         ALPHA(J2,N1)=ALPHA(J2,N1)-
     *                                       XJ*GAMMA(J1)+XI2
                                         ALPHA(J2,J2)=ALPHA(J2,J2)+XJ2
                                    END IF
                             END IF
                         END IF
14                       CONTINUE
                     END IF
13                   CONTINUE
12              CONTINUE
11         CONTINUE
      END IF
C     OUTPUT ALPHA
      WRITE(OUP,111)
111   FORMAT(1X,'ENTRIES OF ALPHA',/)
      WRITE(OUP,112) ((I,J,ALPHA(I,J),J=1,N1),I=1,LN)
112   FORMAT(3(1X,'ALPHA(',I3,',',I3,') = ',E15.8))
C     STEP 20
      IF(LN.GT.1) THEN
      INN=LN-1
      DO 30 I=1,INN
           I1=I+1
           DO 31 J=I1,LN
                CCC=ALPHA(J,I)/ALPHA(I,I)
                DO 32 J1=I1,N1
32              ALPHA(J,J1)=ALPHA(J,J1)-CCC*ALPHA(I,J1)
31         ALPHA(J,I)=0.0
30    CONTINUE
      END IF
      GAMMA(LN)=ALPHA(LN,N1)/ALPHA(LN,LN)
      IF(LN.GT.1) THEN
      DO 33 I=1,INN
           J=LN-I
           CCC=ALPHA(J,N1)
           J1=J+1
           DO 34 KK=J1,LN
34              CCC=CCC-ALPHA(J,KK)*GAMMA(KK)
33         GAMMA(J)=CCC/ALPHA(J,J)
      END IF
C     STEP 21
C     OUTPUT GAMMA
      WRITE(OUP,113)
113   FORMAT(1X,'COEFFICIENTS GAMMA',/)
      DO 115 I=1,LM
      LLL=LL(I)
115   WRITE(OUP,114) I,GAMMA(I),I,(II(I,J),J=1,LLL)
114   FORMAT(1X,'GAMMA(',I3,') = ',E15.8,'  USED WITH NODE',I3,' IN TRI
     1ANGLES ',5(I3,2X))
C     STEP 23
400   CLOSE(UNIT=5)
      CLOSE(UNIT=OUP)
      IF(OUP.NE.6) CLOSE(UNIT=6)
      STOP
      END
      FUNCTION P(X,Y)
      P=1
      RETURN
      END
      FUNCTION Q(X,Y)
      Q=1
      RETURN
      END
      FUNCTION R(X,Y)
      R=0
      RETURN
      END
      FUNCTION F(X,Y)
      F=0
      RETURN
      END
      FUNCTION G(X,Y)
      G=4
      RETURN
      END
      FUNCTION G1(X,Y)
      G1=0.0
      RETURN
      END
      FUNCTION G2(X,Y)
      G2=0
      T=1.0E-05
      IF((.2-T.LE.X.AND.X.LE.4+T).AND.ABS(Y-.2).LE.T) G2=X
      IF((.5-T.LE.X.AND.X.LE..6+T).AND.ABS(Y-.1).LE.T) G2=X
      IF((-T.LE.Y.AND.Y.LE..1+T).AND.ABS(X-.6).LE.T) G2=Y
      IF((-T.LE.X.AND.X.LE..2+T).AND.ABS(Y+X-.4).LE.T)G2=(X+Y)/SQRT(
     * 2.)
      IF((.4-T.LE.X.AND.X.LE..5+T).AND.ABS(Y+X-.6).LE.T)G2=(X+Y)/SQRT(
     * 2.)
      RETURN
      END
      FUNCTION RR(X,Y)
      COMMON A(10,3),B(10,3),C(10,3),XX(11),YY(11),I,I1,I2,J1,J2,J3
      COMMON DELTA
      RR=R(X,Y)*(A(I,I1)+B(I,I1)*X+C(I,I1)*Y)*(A(I,I2)+B(I,I2)*X+C(I,I2)
     1*Y)
      RETURN
      END
      FUNCTION FFF(X,Y)
      COMMON A(10,3),B(10,3),C(10,3),XX(11),YY(11),I,I1,I2,J1,J2,J3
      COMMON DELTA
      FFF=F(X,Y)*(A(I,I1)+B(I,I1)*X+C(I,I1)*Y)
      RETURN
      END
      FUNCTION GG1(X,Y)
      COMMON A(10,3),B(10,3),C(10,3),XX(11),YY(11),I,I1,I2,J1,J2,J3
      COMMON DELTA
      GG1=G1(X,Y)*(A(I,I1)+B(I,I1)*X+C(I,I1)*Y)*(A(I,I2)+B(I,I2)*X+
     1C(I,I2)*Y)
      RETURN
      END
      FUNCTION GG2(X,Y)
      COMMON A(10,3),B(10,3),C(10,3),XX(11),YY(11),I,I1,I2,J1,J2,J3
      COMMON DELTA
      GG2=G2(X,Y)*(A(I,I1)+B(I,I1)*X+C(I,I1)*Y)
      RETURN
      END
      FUNCTION GG3(X,Y)
      COMMON A(10,3),B(10,3),C(10,3),XX(11),YY(11),I,I1,I2,J1,J2,J3
      COMMON DELTA
      GG3=G2(X,Y)*(A(I,I2)+B(I,I2)*X+C(I,I2)*Y)
      RETURN
      END
      FUNCTION GG4(X,Y)
      COMMON A(10,3),B(10,3),C(10,3),XX(11),YY(11),I,I1,I2,J1,J2,J3
      COMMON DELTA
      GG4=G1(X,Y)*(A(I,I1)+B(I,I1)*X+C(I,I1)*Y)**2
      RETURN
      END
      FUNCTION GG5(X,Y)
      COMMON A(10,3),B(10,3),C(10,3),XX(11),YY(11),I,I1,I2,J1,J2,J3
      COMMON DELTA
      GG5=G1(X,Y)*(A(I,I2)+B(I,I2)*X+C(I,I2)*Y)**2
      RETURN
      END
      FUNCTION QQ(FF)
      COMMON A(10,3),B(10,3),C(10,3),XX(11),YY(11),I,I1,I2,J1,J2,J3
      COMMON DELTA
      A1=(XX(J1)+XX(J2)+XX(J3))/3
      A2=(YY(J1)+YY(J2)+YY(J3))/3
      QQ=3*(FF(XX(J1),YY(J1))+FF(XX(J2),YY(J2))+FF(XX(J3),YY(J3)))/60
     1+8*(FF((XX(J1)+XX(J2))/2,(YY(J1)+YY(J2))/2)+FF((XX(J1)+XX(J3))/2,
     1(YY(J1)+YY(J3))/2)+FF((XX(J2)+XX(J3))/2,(YY(J2)+YY(J3))/2))/60
     1+27*FF(A1,A2)/60
      QQ=QQ*ABS(DELTA)/2
      RETURN
      END
      FUNCTION SQ(PP)
      COMMON A(10,3),B(10,3),C(10,3),XX(11),YY(11),I,I1,I2,J1,J2,J3
      COMMON DELTA
      DIMENSION X(101)
      F1(T)=PP(T1*T+X1,T2*T+Y1)*SQRT(T1*T1+T2*T2)
      X1=XX(J1)
      Y1=YY(J1)
      X2=XX(J2)
      Y2=YY(J2)
      T1=X2-X1
      T2=Y2-Y1
      H=.01
      DO 2 II=1,101
2     X(II)=(II-1)*H
      SQ=(F1(X(1))+F1(X(101)))*H/3
      Q1=0.0
      Q2=0.0
      DO 3 II=1,49
      Q1=Q1+2*H*F1(X(2*II+1))/3
3     Q2=Q2+4*H*F1(X(2*II))/3
      Q2=Q2+4*H*F1(X(100))/3
      SQ=SQ+Q1+Q2
      RETURN
      END
