下段程序，求高手给个注解，谢谢
C
C    ****************************************************************
C
      SUBROUTINE GALERKIN(NEQ)
C
C    ****************************************************************
C     FROM THE SYSTEM OF EQUATIONS A X=F THROUGH GALERKIN TECHENIQUE.
      INCLUDE 'TREFFTZ.INC'
    IMPLICIT REAL*8(A-H,O-Z)
      COMMON PI,NNOD,NITP,NTERM,NEXP,NECP
    COMMON/SYSMT/ A(MXEQ,MXEQ),AC(MXEQ,MXEQ),F(MXEQ),FC(MXEQ)
    COMMON/GENRL/ Q(MXEQ),U(MXEQ),VP(MAXI),DPX(MAXI),DPY(MAXI),
     1              RES(3,MAXE),CBC(3,MAXE),VBC(3,MAXE),IBC(3,MAXE)
    COMMON/MESHS/ XI(MAXI),YI(MAXI),X(MAXN),Y(MAXN),IEL(3,MAXE)
    COMMON/SFJAC/ H,SHAPF(3),SHAPD(3),SF(3,MXGP),SD(3,MXGP),RJAC,
     1              UN(2)
      COMMON/EXPAN/ SCL,XC,YC,XCE(MXEX),YCE(MXEX),ALFA(MXEX),BETA(MXEX),
     1              IGE(MXEX),CORNEL(MXEX,MAXE)
    COMMON/ELMNT/ XE(3),YE(3),IJK(3)
    COMMON/TRANS/ XQ,YQ,XX,YY,CB,SB,R,THETA,DPDX,DPDY,DDX,DDY,
     1              FAC,IGEN
    LOGICAL GAMMAQ,CORNEL
C
      NEQ=NEXP*NTERM
C-----SET UP ABCISSAE AND WEIGHTS FOR GAUSS INTEGRATION
      NGAUS=MIN(4*NEXP,MXGP)
    CALL GAULEG(-1.D0,+1.D0,GA,GW,NGAUS)
C-----STORE SHAPE FUNCTIONS AND LOCAL DERIVATIVES AT GAUSS POINTS
      DO 5 IG=1,NGAUS 
       H=GA(IG)
       CALL SHAPES(IG)
5     CONTINUE 
C-----ELEMENT LOOP
      DO 100 M=1,NNEL
       GAMMAU=.TURE.
       IF(IBC(2,M).EQ.1) GAMMAU=.FALSE.
       CALL ARRAYS(M,1)
C--------INTEGRATION LOOP
         DO 130 IG=1,NGAUS
C-----------JACOBIAN & UNIT OUTWARD NORMAL COMPONENTS
            CALL JACOBI(IG)
          FACT=RJAC*GW(IG)
C-----------COORDINATES AND BOUNDARY VALUE OF THE INTEGARTION POINT
            XQ=0.
          YQ=0.
          BV=0.
          DO 140 N=1,3
             XQ=XQ+SF(N,IG)*XE(N)
             YQ=YQ+SF(N,IG)*YE(N)
             BV=BV+SF(N,IG)*VBC(N,M)
140         CONTINUE
C-----------EIGENEXPANSION LOOP
            DO 300 K=1,NEXP
             NCOL=(K-1)*NTERM
             IGEN=IGE(K)
             IF (IGEN.EQ.2) FAC=PI/ALFA(K)
C--------------LOCAL X AXIS OF THE EIGENEXPANSION
               CB=DCOS(DBLE(BETA(K)))
             SB=DSIN(DBLE(BETA(K)))
C--------------GLOBAL ORIGIN COORDINATES OF THE EIGENEXPANSION
               XX=XCE(K)
             YY=YCE(K)
C--------------COMPUTATIONS IN THE LOCAL REFERENCE SYSTEM
               CALL SAMBA(0,.FALSE.)
C--------------EXPANSION-TERMS LOOP
               DO 150 N=1,NTERM
                IF (CORNEL(K,M).OR.R.LT.1.D-6) THEN
                   Q(NCOL+N)=0.
                ELSE
C--------------------TRANSFORMATIONS TO THE GLOBAL REFERENCE SYSTEM
                     CALL SAMBA(J,.FALSE.)
C--------------------NORMAL FLUX
                     Q (NCOL+N)=DDX*UN(1)+DDY*UN(2)
                  ENDIF
C-----------------POTENTIAL
                  U(NCOL+N)=POT(N,R,THE,IGEN,FAC)
150            CONTINUE
300         CONTINUE
            DO 350 K=1,NEXP
             NCOL=(K-1)*NTERM
C--------------LOOP ON ROWS
               DO 60 N=1,NTERM
                I=NCOL+N
                IF (N.EQ.1) THEN 
                   IF (GAMMAU) THEN
C-----------------------EQUATION FOR THE CONSTANT TERM
                        DO 75 J=1,NEQ
                         A(I,J)=A(I,J)+U(J)*FACT
75                      CONTINUE
                        F(I)=F(I)+BV*FACT
                   ENDIF
                ELSE
C--------------------REMAINING EQUATIONS
                     DO 70 J=1,NEQ
                      IF (GAMMAU) THEN
                           A(I,J)=A(I,J)+Q(I)*U(J)*FACT
                      ELSE
                         A(I,J)=A(I,J)-U(I)*U(J)*FACT
                      ENDIF
70                   CONTINUE
                     IF (GAMMAU) THEN
                        F(I)=F(I)+Q(I)*BV*FACT
                     ELSE
                        F(I)=F(I)-U(I)*BV*FACT
                   ENDIF
                ENDIF
60             CONTINUE
350         CONTINUE
130      CONTINUE
100   CONTINUE 
      RETURN
    END