
      PROGRAM DRIVER
C
C     TITLE: CALCULATE GALERKIN COEFFICIENTS USING GNRT.
C     --------------------------------------------------
C
C     THIS IS A DRIVER PROGRAM FOR SUBROUTINE 'GNRT'.  THE USER MUST
C     SUPPLY THE INPUT DATA IN A FILE, WHICH WILL BE DENOTED BY
C     'NAMEIN' IN THIS DRIVER PROGRAM.  THE ACTUAL NAME OF THE
C     INPUT FILE WILL BE READ BY THIS PROGRAM FROM THE STANDARD
C     INPUT UNIT, AND THE FILE NAME MUST BE A CHARACTER STRING.
C     THE INPUT PARAMETERS WILL THEN BE READ FROM FILE 'NAMEIN'.
C
C     THE PARAMETERS THAT MUST BE SUPPLIED IN THE FILE 'NAMEIN' ARE
C     GIVEN BELOW, IN THE ORDER IN WHICH THEY MUST BE GIVEN IN THE FILE.
C
C     NUMSUR        THE INDEX OF THE SURFACE S IN SUBROUTINE 'SURFAC'.
C     A,B,C         THE SURFACE PARAMETERS THAT REFER TO THE ELLIPSOID
C                   PART IN THE DEFINITION OF S.  SEE SUBROUTINE
C                   'SURFAC' FOR MORE DETAILS.
C     ALPHA         AN AUXILIARY PARAMETER THAT MAY BE NEEDED IN
C                   DEFINING THE SURFACE S.
C     DEGREE        THE MAXIMUM DEGREE OF THE SPHERICAL HARMONICS TO
C                   BE USED IN PRODUCING THE GALERKIN COEFFICIENTS.
C                   THERE IS AN UPPER LIMIT ON 'DEGREE'.  IT IS CALLED
C                   'MAXDEG', AND IT IS GIVEN BELOW IN THE PARAMETER
C                   STATEMENT.
C     NINNER        THE INTEGRATION PARAMETER USED IN APPROXIMATING THE
C                   INNER SURFACE INTEGRAL FOR THE GALERKIN COEFFICIENTS.
C     NOUTER        THE INTEGRATION PARAMETER USED IN EVALUATING THE
C                   OUTER SURFACE INTEGRAL.  THERE IS AN UPPER LIMIT
C                   OF 20 FOR 'NOUTER', DUE TO THE PARTICULAR GAUSSIAN
C                   QUADRATURE ROUTINE 'ZEROLG' BEING USED BY 'GNRT'.
C     NAMEDC        THE NAME OF THE OUTPUT FILE THAT WILL CONTAIN THE
C                   DECIMAL VERSION OF THE GALERKIN COEFFICIENTS.
C     NAMEFF        THE NAME OF THE OUTPUT FILE THAT WILL CONTAIN THE
C                   FORMATFREE VERSION OF THE GALERKIN COEFFICIENTS.
C                   THEY ARE OUTPUT IN THIS FORM FOR MORE RAPID INPUT
C                   TO THE PROGRAM 'LAPLAC' FOR SOLVING LAPLACE'S
C                   EQUATION.
C
C     THE FILE NUMBERS FOR THE INPUT AND OUTPUT FILES HAVE BEEN SET TO
C     DEFAULT VALUES, GIVEN IN A DATA STATEMENT BELOW.
C
C     THE SUBROUTINE 'GNRT' REQUIRES A WORKING STORAGE ARRAY 'WORK'
C     TO BE SUPPLIED TO IT.  THE DIMENSION OF THIS ARRAY MUST BE AT
C     LEAST THE MAXIMUM OF THE FOLLOWING QUANTITIES:
C                   IT1=(DEGREE+1)**4 + 2*DEGREE**2 + 7*DEGREE +3
C                       + (3+4*DEGREE)*NOUTER + 7*NINNER
C                   IT2=2*(DEGREE+1)**4
C     THE DIMENSION OF WORK HAS BEEN DEFAULTED TO IT2, USING A MAXIMUM
C     DEGREE SPECIFIED IN THE PARAMETER STATEMENT BELOW.  THIS WAS
C     CHOSEN BECAUSE IT2 IS USUALLY LARGER THAN IT1.  IN ADDITION,
C     THERE IS A CHECK AT THE TIME OF INPUT AS TO WHETHER THE
C     DIMENSION OF 'WORK' IS LARGE ENOUGH.
C
C     MACHINE DEPENDENT CONSTANTS ARE SET THROUGH THE FUNCTION
C     'D1MACH', WHICH IS ONE OF THE SUBPROGRAMS GIVEN IN THE
C     'GNRT' PACKAGE.  THE ROUTINE 'D1MACH' SHOULD BE RESET
C     APPROPRIATELY FOR THE USER'S COMPUTER.
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      INTEGER DEGREE,FILEDC,FILEFF,FILEIN
      CHARACTER NAMEDC*50,NAMEFF*50,NAMEIN*50
      PARAMETER(MAXDEG=10,IWORK=2*(MAXDEG+1)**4)
      DIMENSION WORK(IWORK)
      COMMON/SURPAR/A,B,C,ALPHA,NUMSUR
      COMMON/BLKWRK/WORK
      EXTERNAL SURFAC
C
C     ******************************************
C     THESE ARE THE INPUT/OUTPUT UNIT NUMBERS.
C
      DATA FILEIN/8/,FILEDC/9/,FILEFF/10/
C
C     ******************************************
C
C     OBTAIN THE NAME OF THE INPUT FILE.
      PRINT *, ' GIVE THE NAME OF THE FILE CONTAINING THE'
      PRINT *, ' INPUT DATA.'
      READ(*,'(A)') NAMEIN
C
C     OBTAIN INPUT DATA.
      L=LEN(NAMEIN)
      OPEN(FILEIN,FILE=NAMEIN(1:L))
      READ(FILEIN,*)NUMSUR,A,B,C,ALPHA,DEGREE,NINNER,NOUTER
      READ(FILEIN,'(A)') NAMEDC,NAMEFF
      CLOSE(FILEIN)
C
C     OPEN OUTPUT FILES AND PRINT INPUT PARAMETERS.
      L=LEN(NAMEDC)
      OPEN(FILEDC,FILE=NAMEDC(1:L))
      L=LEN(NAMEFF)
      OPEN(FILEFF,FILE=NAMEFF(1:L),FORM='UNFORMATTED')
      WRITE(FILEDC,1000)NUMSUR,A,B,C,ALPHA,DEGREE,NINNER,NOUTER
1000  FORMAT('1SURFACE NUMBER ',I1,/,' SURFACE PARAMETERS (A,B,C)=',
     *     1PD17.10,2D21.10,/,' SURFACE PARAMETER ALPHA=',D17.10,/,
     *     ' UPPER LIMIT ON DEGREE=',I2,/,' INTEGRATION PARAMETERS ',
     *     'NINNER, NOUTER = ',I2,I5,///)
C
C     PRODUCE NEEDED SIZE FOR WORKING ARRAY 'WORK'.
      IT1=(DEGREE+1)**4+7*NINNER+(3+4*DEGREE)*NOUTER+2*DEGREE**2
     *     +7*DEGREE+3
      IT2=2*(DEGREE+1)**4
      NCHECK=MAX(IT1,IT2)
C
      IF(NCHECK .LE. IWORK) THEN
C         THERE IS ENOUGH TEMPORARY STORAGE IN ARRAY 'WORK'.
C         CALL GNRT TO GENERATE THE GALERKIN COEFFICIENTS.
          CALL GNRT(SURFAC,DEGREE,NINNER,NOUTER,FILEDC,FILEFF,WORK)
      ELSE
C         THERE IS NOT ENOUGH WORKING STORAGE IN ARRAY 'WORK'.
          WRITE(FILEDC,1001) NCHECK
1001      FORMAT(' THERE IS NOT ENOUGH STORAGE IN THE ARRAY WORK.',/,
     *           ' YOU NEED AT LEAST ',I5,' POSITIONS.')
      END IF
      CLOSE(FILEDC)
      CLOSE(FILEFF)
      STOP
      END

      SUBROUTINE SURFAC(Q,QCROSS,ST,CT,SP,CP,IFLAG)
C     -----------------
C
C     THIS CALCULATES A POINT Q ON THE SURFACE S, BASED ON S BEING THE
C     CONTINUOUS ONE-TO-ONE IMAGE OF THE UNIT SPHERE U.  LET Q=(X,Y,Z)
C     DENOTE THE DESIRED POINT ON S.  THEN IT IS TO CORRESPOND TO
C     THE POINT
C          UQ = (UX,UY,UZ)
C     ON U, WITH
C          UX = SIN(THETA)*COS(PHI) = ST*CP
C          UY = SIN(THETA)*SIN(PHI) = ST*SP
C          UZ = COS(THETA)          = CT
C
C     INPUT:
C     ST,CT,SP,CP   THESE ARE THE GIVEN TRIGONOMETRIC FUNCTION VALUES
C           FOR THE POINT UQ.
C     IFLAG = 0   COMPUTE Q.
C           = 1   COMPUTE Q AND THE CROSS-PRODUCT
C                 QCROSS = DPHI(Q) X DTHETA(Q)/SIN(THETA)
C           THE NOTATION DPHI(Q) MEANS THE DERIVATIVE OF Q REGARDED AS
C           AS A FUNCTION OF PHI, AND SIMILARLY FOR DTHETA(Q).  THE QUANTITY
C           QCROSS IS USED IN COMPUTING THE NORMAL TO S AT Q; AND ITS
C           MAGNITUDE IS THE JACOBIAN IN THE CHANGE OF THE SURFACE
C           DIFFERENTIAL FOR CHANGING AN INTEGRAL OVER S TO AN INTEGRAL
C           OVER U.  THE DIVISION BY SIN(THETA) REMOVES A SINGULARITY
C           IN THE MAPPING DUE TO THE USE OF SPHERICAL COORDINATES
C           IN REPRESENTING POINTS OF U.
C     OUTPUT:
C     Q       THE POINT ON S CORRESPONDING TO UQ ON THE SPHERE U.
C     QCROSS  THE CROSS-PRODUCT REFERRED TO ABOVE.
C
C     THIS PARTICULAR SUBROUTINE ASSUMES THAT THE SURFACE S IS GIVEN BY
C            Q = (X,Y,Z)
C              = R*(A*UX,B*UY,C*UZ)
C     WITH R = R(UQ) A POSITIVE SMOOTH CONTINUOUS FUNCTION ON U.  THE
C     CONSTANTS A,B,C ARE ASSUMED POSITIVE; AND FOR R=1, THE SURFACE
C     WILL BE AN ELLIPSOID WITH SEMI-AXES A,B,C.  THIS GENERAL FORM
C     IS EQUIVALENT TO ASSUMING THAT THE REGION ENCLOSED BY S IS
C     STARLIKE WITH RESPECT TO THE ORIGIN.  IT IS ALSO EQUIVALENT TO
C     THE STANDARD PARAMETRIC WAY OF DEFINING A SURFACE S IN TERMS
C     OF VARIABLES PHI AND THETA VARYING OVER  (O,2*PI) AND (0,PI),
C     RESPECTIVELY.
C
C     IT IS ASSUMED THAT R CAN BE EXTENDED AS A SMOOTH FUNCTION TO A
C     NEIGHBORHOOD OF U, SO THAT THE FIRST DERIVATIVES OF R(UX,UY,UZ)
C     CAN BE COMPUTED AT ALL POINTS OF U.  THE USER MUST SUPPLY R AND
C     ITS DERIVATIVES WITH RESPECT TO UX,UY, AND UZ, DENOTED BY DXR,
C     DYR, AND DZR, RESPECTIVELY.  EXAMPLES ARE GIVEN BELOW.
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION Q(3),QCROSS(3)
      COMMON/MACHCN/U100,U100SQ
C     THE VARIABLE U100 IS 100 TIMES THE MACHINE UNIT ROUND.
      COMMON/SURPAR/A,B,C,ALPHA,NUMSUR
C
C     COMPUTE UQ.
      UX = ST*CP
      UY = ST*SP
      UZ = CT
C
C     ***** NOTE TO USER *****************************************************
C     THIS BEGINS THE SECTION THAT THE USER MUST SPECIFY.
C
      GO TO (10,20,30), NUMSUR
C
C     SURFACE #1.
C     S IS AN ELLIPSOID.
10    R = 1.0D0
      IF(IFLAG .EQ. 0) GO TO 100
      DXR = 0.0D0
      DYR = 0.0D0
      DZR = 0.0D0
      GO TO 100
C
C     SURFACE #2.
C     S IS A PEANUT-SHAPED REGION, BASED ON THE OVALS OF CASSINI.
C     THE PARAMETER ALPHA IS A NON-ZERO POSITIVE CONSTANT.  AS IT
C     APPROACHES ZERO, THE SURFACE S BECOMES MORE CONSTRICTED IN
C     THE XY-PLANE.
20    C2T = 2.0D0*UZ*UZ - 1.0D0
      TEMP = SQRT(ALPHA + C2T*C2T)
      R = SQRT(C2T + TEMP)
      IF(IFLAG .EQ. 0) GO TO 100
      DXR = 0.0D0
      DYR = 0.0D0
      DZR=2.0D0*UZ*(1.0D0+C2T/TEMP)/R
      GO TO 100
C
C     SURFACE #3.
C     S IS A HEART-SHAPED REGION, INCREASINGLY SO AS THE POSITIVE
C     CONSTANT ALPHA INCREASES.
30    TEMP=1.0D0+ALPHA*(UZ+1.0D0)**2
      R=2.0D0-1.0D0/TEMP
      IF(IFLAG .EQ. 0) GO TO 100
      DXR=0.0D0
      DYR=0.0D0
      DZR=2.0D0*ALPHA*(UZ+1.0D0)/TEMP**2
      GO TO 100
C
C     THIS ENDS THE SECTION THAT THE USER MUST SUPPLY.
C     ************************************************************************
C
C     COMPUTE THE POINT Q.
100   Q(1) = A*R*UX
      Q(2) = B*R*UY
      Q(3) = C*R*UZ
      IF(IFLAG .EQ. 0) RETURN
C
C     COMPUTE CROSS-PRODUCT.
C
      IF(ST .LE. U100) GO TO 200
C     BEGIN BY COMPUTING DERIVATIVES OF Q WITH RESPECT TO PHI AND THETA.
      DPHIR = -UY*DXR + UX*DYR
      DTHER = CT*(CP*DXR+SP*DYR) - ST*DZR
      DPHIQX = A*ST*(-SP*R + CP*DPHIR)
      DTHEQX = A*CP*(CT*R + ST*DTHER)
      DPHIQY = B*ST*(CP*R + SP*DPHIR)
      DTHEQY = B*SP*(CT*R + ST*DTHER)
      DPHIQZ = C*CT*DPHIR
      DTHEQZ = C*(-ST*R + CT*DTHER)
      QCROSS(1) = (DPHIQY*DTHEQZ-DTHEQY*DPHIQZ)/ST
      QCROSS(2) = (DPHIQZ*DTHEQX-DTHEQZ*DPHIQX)/ST
      QCROSS(3) = (DPHIQX*DTHEQY-DTHEQX*DPHIQY)/ST
      RETURN
C
C     COMPUTE CROSS-PRODUCT FOR Q AT A POLE.
200   QCROSS(1) = B*C*R*DXR
      QCROSS(2) = A*C*R*DYR
      QCROSS(3) = -A*B*UZ*R*R
      RETURN
      END

      SUBROUTINE GNRT(SURFAC,DEGREE,NINNER,NOUTER,FILEDC,FILEFF,WORK)
C     ---------------
C
C        THIS ROUTINE CALCULATES THE GALERKIN COEFFICIENTS FOR SOLVING
C     THE DOUBLE LAYER POTENTIAL INTEGRAL EQUATION
C             (2*PI + K)*RHO = F .
C     IN THIS FORMULA, IT IS ASSUMED THAT THE VARIABLES HAVE BEEN
C     CHANGED SO THAT ALL FUNCTIONS ARE DEFINED ON THE UNIT SPHERE U.
C     THIS IS DONE AUTOMATICALLY IN SUBROUTINE 'GNRT', USING THE
C     SURFACE S AS SPECIFIED IN THE USER SUPPLIED SUBROUTINE 'SURFAC'.
C
C        LET H(I) DENOTE THE SPHERICAL HARMONIC WITH INDEX I, BASED
C     ON THE ORDERING OF THEM AS DEFINED IN FUNCTION 'LOCATN' AND
C     SUBROUTINE 'INDX'.  THEN 'GNRT' PRODUCES THE GALERKIN COEFFICIENTS
C             (H(I),K*H(J)) ,
C     FOR 0 .LE. DEG(H(I)),DEG(H(J)) .LE. N, N='DEGREE' AN INPUT
C     VARIABLE.  THESE ARE OUTPUT IN N TABLES.  TABLE #L CONTAINS
C     THE GALERKIN COEFFICIENTS FOR ALL HARMONICS OF DEGREE .LE. L;
C     AND 1 .LE. L .LE. N.
C
C     INPUT/OUTPUT
C          SURFAC: A USER SUPPLIED SUBROUTINE SPECIFYING THE SURFACE S.
C              SEE THE EXAMPLE DRIVER PROGRAM FOR THE SPECIFICATION OF
C              'SURFAC'.
C          DEGREE: THE MAXIMUM DEGREE OF THE SPHERICAL HARMONICS TO BE
C              USED IN CALCULATING THE GALERKIN COEFFICIENTS.
C          NINNER: THE INNER INTEGRAL K*H(J) IS TO BE INTEGRATED BY THE
C              RULE SPECIFIED IN THE ACCOMPANYING ARTICLE, WITH THE
C              INTEGRATION PARAMETER 'NINNER'.
C          NOUTER: THE OUTER INTEGRAL, THE INNER PRODUCT (H(I),K*H(J)),
C              IS TO BE INTEGRATED WITH THE SPHERICAL GAUSSIAN QUADRATURE
C              FORMULA WITH INTEGRATION PARAMETER 'NOUTER'.  THERE IS A
C              LIMIT OF 20 FOR 'NOUTER', DUE TO THE SUBROUTINE 'ZEROLG'
C              THAT IS BEING USED TO CALCULATE THE GAUSS-LEGENDRE NODES
C              AND WEIGHTS.  USE ANOTHER SUCH ROUTINE WITH A LARGER UPPER
C              LIMIT TO REMOVE THE LIMIT ON 'NOUTER'.
C          FILEDC: THE FILE NUMBER FOR THE OUTPUT FILE CONTAINING THE
C              DECIMAL VERSION OF THE GALERKIN COEFFICIENTS.
C          FILEFF: THE FILE NUMBER FOR THE OUTPUT FILE CONTAINING THE
C              FORMATFREE VERSION OF THE GALERKIN COEFFICIENTS.
C          WORK: A TEMPORARY WORKING ARRAY FOR USE IN GNRT.  ITS ORDER
C              MUST BE AT LEAST EQUAL TO THE MAXIMUM OF THE FOLLOWING
C              TWO QUANTITIES:
C                 IT1 = (DEGREE+1)**4 + 2*DEGREE**2 + 7*DEGREE + 3
C                       + (3+4*DEGREE)*NOUTER + 7*NINNER
C                 IT2 = 2*(DEGREE+1)**4
C              IN THE PRESENT ROUTINE, THE ARRAY 'WORK' IS BROKEN APART
C              INTO SMALLER ONE AND TWO DIMENSIONAL ARRAYS, FOR USE IN THE
C              FOLLOWING SUBROUTINE 'GNRT2'.
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      INTEGER DEGREE,DEG2,DEGSQ,FILEDC,FILEFF
      DIMENSION WORK(*)
      COMMON/MACHCN/U100,U100SQ
C     THE VARIABLE U100 IS 100 TIMES THE MACHINE UNIT ROUND.
      EXTERNAL SURFAC
C
C     SET MACHINE DEPENDENT CONSTANT.
      U100=100*D1MACH(4)
      U100SQ=U100**2
C
C     INITIALIZE VARIABLES TO BE USED IN GNRT2.
      DEG2=2*DEGREE+1
      DEGSQ=(DEGREE+1)**2
      NINNR2=2*NINNER
      NOUTR2=2*NOUTER
      LGDORD=1+(DEGREE*(DEGREE+3))/2
C
C     SET UP INITIAL ADDRESSES FOR BREAKING WORK INTO SMALLER
C     WORKING ARRAYS.
      I1=1
      I2=I1+DEGSQ**2
      I3=I2+NINNR2
      I4=I3+NINNR2
      I5=I4+DEGREE
      I6=I5+DEGREE
      I7=I6+NINNER
      I8=I7+NINNER
      I9=I8+NINNER
      I10=I9+LGDORD
      I11=I10+DEGREE*NOUTR2
      I12=I11+DEGREE*NOUTR2
      I13=I12+NOUTER
      I14=I13+NOUTER
      I15=I14+NOUTER
      I16=I15+LGDORD
C
C     GO TO GNRT2 FOR ACTUAL COMPUTATION OF GALERKIN COEFFICIENTS.
      CALL GNRT2(SURFAC,DEGREE,NINNER,NOUTER,FILEDC,WORK(I1),WORK(I2),
     *     WORK(I3),WORK(I4),WORK(I5),WORK(I6),WORK(I7),WORK(I8),
     *     WORK(I9),WORK(I10),WORK(I11),WORK(I12),WORK(I13),WORK(I14),
     *     WORK(I15),WORK(I16),NINNR2,NOUTR2,LGDORD,DEG2,DEGSQ)
C
C     EXPAND AND WRITE THE FILE OF GALERKIN COEFFICIENTS.
      CALL EXPAND(WORK(I1),WORK(I2),DEGREE,DEGSQ,NINNER,NOUTER,
     *     FILEFF)
      RETURN
      END

      SUBROUTINE GNRT2(SURFAC,DEGREE,NINNER,NOUTER,FILEDC,
     *     KERMAT,SNPHI,CSPHI,SNPHI2,CSPHI2,SNTHI,CSTHI,WI,VALGDI,
     *     SNPHE,CSPHE,SNTHE,CSTHE,WE,VALGDE,TEMPKH,NINNR2,NOUTR2,
     *     LGDORD,DEG2,DEGSQ)
C     ----------------
C
C     THIS ROUTINE IS WHERE THE GALERKIN COEFFICIENTS ARE
C     ACTUALLY CALCULATED.
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION KERMAT,KERNEL
      INTEGER DEGREE,DEG2,DEGSQ,FILEDC
      COMMON/BLOCKR/V,STI,CTI,SPI,CPI,UPX,UPY,UPZ
      DIMENSION KERMAT(DEGSQ,DEGSQ),SNPHI(NINNR2),CSPHI(NINNR2),
     *     SNPHI2(DEGREE),CSPHI2(DEGREE),SNTHI(NINNER),CSTHI(NINNER),
     *     WI(NINNER),VALGDI(LGDORD),SNPHE(DEGREE,NOUTR2),
     *     CSPHE(DEGREE,NOUTR2),SNTHE(NOUTER),CSTHE(NOUTER),WE(NOUTER),
     *     VALGDE(LGDORD),TEMPKH(DEGSQ)
      DIMENSION P(3),Q(3),QCROSS(3),V(3)
      PARAMETER(ZERO=0.0D0,P5=0.5D0,ONE=1.0D0,TWO=2.0D0,
     *          FOUR=4.0D0)
C
C     CALCULATE NODES,WEIGHTS, AND TRIG FUNCTIONS FOR INTERIOR AND
C     EXTERIOR INTEGRATIONS.
      PI=FOUR*ATAN(ONE)
      TWOPI=TWO*PI
      HI=PI/NINNER
      HE=PI/NOUTER
C     CALCULATE TRIG VALUES FOR INNER INTEGRATION.
      CALL INTMI(-ONE,ONE,NINNER,CSTHI,WI)
      DO 10 J=1,NINNER
10      SNTHI(J)=SIN(ACOS(CSTHI(J)))
      DO 20 L=1,NINNR2
        PHI=(L-P5)*HI
        SNPHI(L)=SIN(PHI)
20      CSPHI(L)=COS(PHI)
C
C     CALCULATE TRIG FUNCTIONS FOR EXTERIOR INTEGRATION.
      CALL ZEROLG(NOUTER,CSTHE,WE,-ONE,ONE)
      DO 30 J=1,NOUTER
        THE=ACOS(CSTHE(J))
30      SNTHE(J)=SIN(THE)
      DO 40 L=1,NOUTR2
        PHE=L*HE
        SNP=SIN(PHE)
        CSP=COS(PHE)
        SNPHE(1,L)=SNP
        CSPHE(1,L)=CSP
        DO 40 M=2,DEGREE
          SNPHE(M,L)=SNPHE(M-1,L)*CSP+CSPHE(M-1,L)*SNP
40        CSPHE(M,L)=CSPHE(M-1,L)*CSP-SNPHE(M-1,L)*SNP
C
C     INITIALIZE FOR BEGINNING MAIN LOOP TO COMPUTE ALL (TNM,K(TRS)).
C     SET KERMAT TO ZERO.
      DO 50 I=1,DEGSQ
      DO 50 J=1,DEGSQ
50      KERMAT(I,J)=ZERO
C
C     *********
C     MAIN LOOP
C     *********
C
      DO 150 I=1,NOUTER
        CALL LEGEND(CSTHE(I),DEGREE,VALGDE)
        DO 150 J=1,NOUTR2
C         PRODUCE SURFACE COORDINATES OF EXTERNAL VARIABLE POINT P=(X,Y,Z).
          ST=SNTHE(I)
          CT=CSTHE(I)
          SP=SNPHE(1,J)
          CP=CSPHE(1,J)
          CALL SURFAC(P,DUMMY,ST,CT,SP,CP,0)
C
C         PRODUCE NEEDED DATA FOR ROTATION OF COORDINATES.
          UPX=CP*ST
          UPY=SP*ST
          UPZ=CT
          ALPHA=SIGN(ONE,-UPZ)
          V(3)=SQRT((ONE-UPZ/ALPHA)/TWO)
          PV=-TWO*ALPHA*V(3)
          V(2)=UPY/PV
          V(1)=UPX/PV
C
C         INITIALIZE TEMPKH TO ZERO.
          DO 60 II=1,DEGSQ
60          TEMPKH(II)=ZERO
C
C         LOOP TO INTEGRATE K(TRS) EVALUATED AT (THETA(I),PHI(J)).
          DO 100 II=1,NINNER
          DO 100 JJ=1,NINNR2
            UPX=CSPHI(JJ)*SNTHI(II)
            UPY=SNPHI(JJ)*SNTHI(II)
            UPZ=CSTHI(II)
C           PRODUCE ROTATED POINT
            CALL ROTATE
C
C           PRODUCE NEEDED TRIG AND LEGENDRE FUNCTIONS.
            CALL LEGEND(CTI,DEGREE,VALGDI)
            CSPHI2(1)=CPI
            SNPHI2(1)=SPI
            DO 70 M=2,DEGREE
              CSPHI2(M)=CPI*CSPHI2(M-1)-SPI*SNPHI2(M-1)
70            SNPHI2(M)=SPI*CSPHI2(M-1)+CPI*SNPHI2(M-1)
C
C           CALCULATE SURFACE POSITION
            CALL SURFAC(Q,QCROSS,STI,CTI,SPI,CPI,1)
            TEMPK=KERNEL(P,Q,QCROSS)
C
            DO 80 NR=1,DEGREE+1
80            TEMPKH(NR)=TEMPKH(NR)+WI(II)*TEMPK*(VALGDI(NR)-VALGDE(NR))
            IBRS=DEGREE
            IB=DEGREE+1
            DO 90 MS=1,DEGREE
            DO 90 NR=MS,DEGREE
              IB=IB+1
              IBRS=IBRS+2
              TEMPKH(IBRS)=TEMPKH(IBRS)+WI(II)*TEMPK*(VALGDI(IB)
     *              *CSPHI2(MS)-VALGDE(IB)*CSPHE(MS,J))
90            TEMPKH(IBRS+1)=TEMPKH(IBRS+1)+WI(II)*TEMPK*
     *               (VALGDI(IB)*SNPHI2(MS)-VALGDE(IB)*SNPHE(MS,J))
100         CONTINUE
C
C         ADD BACK IN THE CONTRIBUTION SUBTRACTED EARLIER.
          DO 110 NR=1,DEGREE+1
110         TEMPKH(NR)=HI*TEMPKH(NR)+TWOPI*VALGDE(NR)
          IB=DEGREE+1
          IBRS=DEGREE
          DO 120 MS=1,DEGREE
          DO 120 NR=MS,DEGREE
            IB=IB+1
            IBRS=IBRS+2
            TEMPKH(IBRS)=HI*TEMPKH(IBRS)+TWOPI*VALGDE(IB)*CSPHE(MS,J)
120         TEMPKH(IBRS+1)=HI*TEMPKH(IBRS+1)
     *         +TWOPI*VALGDE(IB)*SNPHE(MS,J)
C
C         PRODUCE THE NEW CONTRIBUTION, AT (TH(I),PH(J)), TO THE
C         ELEMENTS (TNM,K(TRS)) STORED IN KERMAT.
          IBNM=0
          IB=0
          TEMP1=WE(I)*HE
          DO 130 N=1,DEGREE+1
            IBNM=IBNM+1
            IB=IB+1
            TEMP=TEMP1*VALGDE(N)
            DO 130 IBRS=1,DEGSQ
130           KERMAT(IBNM,IBRS)=KERMAT(IBNM,IBRS)+TEMP*TEMPKH(IBRS)
          IBNM=IBNM-1
          DO 140 M=1,DEGREE
          DO 140 N=M,DEGREE
            IBNM=IBNM+2
            IB=IB+1
            TEMP=TEMP1*VALGDE(IB)
            DO 140 IBRS=1,DEGSQ
              KERMAT(IBNM,IBRS)=KERMAT(IBNM,IBRS)
     *            +TEMP*CSPHE(M,J)*TEMPKH(IBRS)
140           KERMAT(IBNM+1,IBRS)=KERMAT(IBNM+1,IBRS)
     *           +TEMP*SNPHE(M,J)*TEMPKH(IBRS)
150       CONTINUE
C
C     OUTPUT VALUES OF KERMAT, AND ALSO STORE THEM IN COEFF.
      DO 200 I=1,DEGSQ
        CALL INDX(I,DEGREE,N,M,L)
        WRITE(FILEDC,1003)I,N,M,L
1003    FORMAT(' ROW=',I2,5X,'N,M,L=',3I3)
        WRITE(FILEDC,1002)(KERMAT(I,J),J=1,DEGSQ)
1002    FORMAT(1PD20.10,5D20.10)
200     CONTINUE
      RETURN
      END

      SUBROUTINE ROTATE
C     -----------------
C
C     CALCULATE (TX,TY,TZ)=(I-2*V*VT)*(X,Y,Z), THE
C     HOUSEHOLDER TRANSFORMATION OF (X,Y,Z).
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON/BLOCKR/V,ST,CT,SP,CP,X,Y,Z
      DIMENSION V(3)
      COMMON/MACHCN/U100,U100SQ
C     THE VARIABLE U100 IS 100 TIMES THE MACHINE EPSILON.
      PARAMETER(ZERO=0.0D0,ONE=1.0D0,TWO=2.0D0)
C
C     CALCULATE ROTATED POINT T.
      VM=TWO*(V(1)*X+V(2)*Y+V(3)*Z)
      TX=X-VM*V(1)
      TY=Y-VM*V(2)
      TZ=Z-VM*V(3)
C
C     CALCULATE TRIG FUNCTIONS FOR ANGLES DETERMINED BY T.
      CT=TZ
      ST=SQRT(TX*TX+TY*TY)
      IF(ST .LE. U100) THEN
          CP=ONE
          SP=ZERO
      ELSE
          SP=TY/ST
          CP=TX/ST
      END IF
      RETURN
      END

      FUNCTION KERNEL(P,Q,QCROSS)
C     ---------------
C
C     EVALUATE THE NORMAL DERIVATIVE OF 1/ABS(P-Q) WITH RESPECT TO THE
C     POINT Q ON THE SURFACE S.  ALSO INCLUDE THE EFFECT OF THE CHANGE
C     OF SURFACE DIFFERENTIAL IN TRANSFORMING TO AN INTEGRAL OVER THE
C     UNIT SPHERE U.
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION KERNEL
      DIMENSION P(3),Q(3),QCROSS(3)
      COMMON/MACHCN/U100,U100SQ
C     THE VARIABLE U100 IS 100 TIMES THE MACHINE EPSILON.
      PARAMETER(ZERO=0.0D0)
C
      RS=(P(1)-Q(1))**2+(P(2)-Q(2))**2+(P(3)-Q(3))**2
      IF(RS .LE. U100SQ) THEN
          KERNEL=ZERO
      ELSE
          DENOM=RS*SQRT(RS)
          TEMP=ZERO
          DO 10 J=1,3
10          TEMP=TEMP+(P(J)-Q(J))*QCROSS(J)
          KERNEL=TEMP/DENOM
      END IF
      RETURN
      END

      SUBROUTINE LEGEND(X,NMAX,VALUES)
C     -----------------
C
C        THIS ROUTINE WILL EVALUATE THE NORMALIZED ASSOCIATED LEGENDRE
C     FUNCTIONS OF THE FIRST KIND AT X,
C                                 J
C                         Y(N,J)=P (X)
C                                 N
C
C     THE DEFINITION IS TAKEN FROM CHAPTER 8 OF ABRAMOWITZ AND STEGUN,
C     HANDBOOK OF MATHEMATICAL FUNCTIONS. THESE WILL BE EVALUATED AND
C     STORED IN THE VECTOR 'VALUES' FOR 0 .LE. J .LE. N .LE. NMAX, WITH
C     NMAX .GE. 1.  THE VALUE Y(N,J) IS STORED IN JB+N-J, WITH
C                 JB=(2+(2*NMAX+3)*J-J**2)/2
C     JB IS THE INDEX FOR Y(J,J).  THE NORMALIZATION IS DONE SO THAT
C     THE SUBSEQUENT CONSTRUCTION OF SPHERICAL HARMONICS WILL RESULT
C     IN FUNCTIONS OF NORM 1 OVER THE UNIT SPHERE.
C        THE EVALUATION IS DONE USING THE TRIPLE RECURSION RELATION
C     (8.5.3) ON PAGE 334 OF THE ABOVE REFERENCE. THIS IS PRECEDED BY
C     A SPECIAL LOOP TO SET Y(N,N) AND Y(N,N-1).  FOLLOWING THAT, THE
C     FUNCTIONS ARE NORMALIZED APPROPRIATELY.
C
C     NOTE: THERE IS AN UPPER LIMIT ON NMAX, SPECIFIED BY THE
C     VARIABLE 'MAXDEG' IN THE PARAMETER STATEMENT BELOW.  TO
C     INCREASE THE UPPER LIMIT, MERELY INCREASE 'MAXDEG'.
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(MAXDEG=20,LFORD=2*MAXDEG)
      DIMENSION VALUES(*),FACT(0:LFORD)
      PARAMETER(ONE=1.0D0,TWO=2.0D0,FOUR=4.0D0)
C
C     SET INITIAL VALUES Y(N,N) AND Y(N,N-1), 0 .LE. N .LE. NMAX.
      RT=SQRT(ONE-X*X)
      VALUES(1)=ONE
      VALUES(2)=X
      IF(NMAX .EQ. 1) THEN
          VALUES(3)=-RT
          GO TO 30
      END IF
C
      JB=1
      DO 10 J=1,NMAX-1
        JBN=JB+NMAX-J+2
        VALUES(JBN)=-(2*J-1)*RT*VALUES(JB)
        VALUES(JBN+1)=-(2*J+1)*RT*VALUES(JB+1)
10      JB=JBN
      VALUES(JB+2)=-(2*NMAX-1)*RT*VALUES(JB)
C
C     PRODUCE THE REMAINDER OF THE TABLE VALUES Y(N,J).
      JB=1
      DO 20 J=0,NMAX-2
        JB=JB+2
        DO 20 N=J+2,NMAX
          VALUES(JB)=((2*N-1)*X*VALUES(JB-1)-(N+J-1)*VALUES(JB-2))/(N-J)
20        JB=JB+1
C
C     NORMALIZE THE LEGENDRE FUNCTIONS.
30    FACT(0)=ONE
      DO 40 J=1,2*NMAX
40      FACT(J)=J*FACT(J-1)
      IB=0
      PI=4*ATAN(ONE)
      TWOPI=TWO*PI
      FOURPI=FOUR*PI
      DO 50 N=0,NMAX
50      VALUES(N+1)=VALUES(N+1)*SQRT((2*N+1)/FOURPI)
      JB=NMAX+1
      DO 60 M=1,NMAX
      DO 60 N=M,NMAX
        JB=JB+1
        COEFF=SQRT(((2*N+1)/TWOPI)*FACT(N-M)/FACT(N+M))
60      VALUES(JB)=COEFF*VALUES(JB)
      RETURN
      END

      FUNCTION LOCATN(N,M,L,NDEG)
C     ---------------
C
C     GIVES THE INDEX OF THE SPHERICAL HARMONIC BASIS FUNCTION
C     WITH INDICES (N,M,L).  FOR Q=(X,Y,Z) ON THE UNIT SPHERE, WRITE
C        Q = (SIN(THETA)*COS(PHI),SIN(THETA)*SIN(PHI),COS(THETA))
C     FOR M=0, THE SPHERICAL HARMONIC CORRESPONDING TO (N,M,L) IS
C
C        H(Q) = P (Z).
C                N
C
C     FOR M>0 AND L=0, THE SPHERICAL HARMONIC IS
C
C                M
C        H(Q) = P (Z)*COS(M*PHI)
C                N
C
C     FOR M>0 AND L=1, THE SPHERICAL HARMONIC IS
C
C                M
C        H(Q) = P (Z)*SIN(M*PHI)
C                N
C
C                              M
C     THE FUNCTIONS P (Z) AND P (Z) ARE THE LEGENDRE POLYNOMIALS AND
C                    N         N
C     THE ASSOCIATED LEGENDRE FUNCTIONS, DEFINED IN SUBROUTINE 'LEGEND'.
C
      IF(M .EQ. 0) THEN
          LOCATN=N+1
      ELSE
          LOCATN=NDEG*(2*M-1)-M*(M-3)+2*(N-M)+L
      END IF
      RETURN
      END

      SUBROUTINE INDX(I,NDEG,N,M,L)
C     ---------------
C
C     THIS PROGRAM TAKES THE INDEX I OF A SPHERICAL HARMONIC,
C              1 .LE. I .LE. (NDEG+1)**2,
C     AND IT FINDS THE INDICES (N,M,L) TO WHICH I CORRESPONDS.
C     SEE FUNCTION 'INDX' FOR THE DEFINITION OF (N,M,L).
C
      KOL(MM)=NDEG*(MM-1)-(MM*(MM-3))/2-1
C
      IF(I .LE. NDEG+1) THEN
          N=I-1
          M=0
          L=0
          RETURN
      END IF
C
      J=I-NDEG-2
      IF(J .EQ. 2*(J/2)) THEN
          L=0
      ELSE
          L=1
      END IF
      J=(J-L)/2
      DO 40 M=1,NDEG
        IF(J .LT. KOL(M+1)) THEN
            N=J-KOL(M)+M
            RETURN
        END IF
40      CONTINUE
      END

      SUBROUTINE EXPAND(KERMAT,DUMMY,DEGREE,DEGSQ,NINNER,NOUTER,
     *     FILEFF)
C     -----------------
C
C     THIS PROGRAM EXPANDS A TABLE OF GALERKIN COEFFICIENTS FROM THOSE FOR
C     DEGREES .LE. N:=DEGREE, AS INPUT.  ON OUTPUT, THERE WILL BE N TABLES,
C     FOR GALERKIN COEFFICIENTS OF DEGREES .LE. L, WITH 1 .LE. L .LE. N.
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      INTEGER DEGREE,DEGSQ,FILEFF
      DOUBLE PRECISION KERMAT
      DIMENSION KERMAT(DEGSQ,DEGSQ),DUMMY(DEGSQ,DEGSQ)
C
      REWIND FILEFF
C
C     LOOP TO CONSTRUCT TABLES FOR SMALLER VALUES OF DEGREE.
      DO 20 ND=1,DEGREE
        NB=(ND+1)**2
        DO 10 MR=0,ND
          DO 10 NR=MR,ND
            DO 10 LR=0,1
              INEW=LOCATN(NR,MR,LR,ND)
              IOLD=LOCATN(NR,MR,LR,DEGREE)
              DO 10 MC=0,ND
                DO 10 NC=MC,ND
                  DO 10 LC=0,1
                    JNEW=LOCATN(NC,MC,LC,ND)
                    JOLD=LOCATN(NC,MC,LC,DEGREE)
                    DUMMY(INEW,JNEW)=KERMAT(IOLD,JOLD)
10                  CONTINUE
        WRITE(FILEFF)ND,NINNER,NOUTER
20      WRITE(FILEFF)((DUMMY(I,J),J=1,NB),I=1,NB)
      RETURN
      END

      SUBROUTINE INTMI(A,B,N,X,W)
C     ----------------
C
C     CALCULATE THE NODES AND WEIGHTS FOR THE INTEGRATION
C
C           B                       N
C       INT   F(X)*DX  APPROX=  SUM   W(I)*F(X(I))
C           A                      I=1
C
C     THE METHOD IS DUE TO IRI, MORIGUTI, AND TAKASAWA.
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(ONE=1.0D0,TWO=2.0D0)
      DIMENSION X(*),W(*)
C
C     INITIALIZE
      H=TWO/(N+1)
      IF(N .EQ. 2*(N/2)) THEN
          M=N/2
      ELSE
          M=(N+1)/2
      END IF
      FAC=B-A
      FAC2=H*(B-A)
C
C     CALCULATE NODES AND WEIGHTS
      DO 10 I=1,M
        CALL TRNSFM(I*H-ONE,T,S)
        X(I)=A+T*FAC
        X(N+1-I)=B-T*FAC
        W(I)=FAC2*S
        W(N+1-I)=W(I)
10      CONTINUE
      RETURN
      END


      SUBROUTINE TRNSFM(Z,FI,SI)
C     -----------------
C
C..........................................................................
C
C     FOR GIVEN Z, CALCULATE
C
C        SI(Z) = CON*EXP(-4/(1-Z*Z))
C
C     AND
C                    Z
C        FI(Z) = INT   SI(T)*DT
C                   -1
C
C     WHERE CON IS CHOSEN TO NORMALIZE FI(1)=1.
C     FI(Z) IS COMPUTED USING A TRUNCATED CHEBYSHEV SERIES.
C
C     THIS ROUTINE IS TAKEN FROM THE PAPER:
C         IAN ROBINSON AND ELISE DE DONCKER, ALGORITHM 45: AUTOMATIC
C         COMPUTATION OF IMPROPER INTEGRALS OVER A BOUNDED OR
C         UNBOUNDED PLANAR REGION, COMPUTING 27(1981), PP. 253-284.
C
C..........................................................................
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      INTEGER ARG(5)
      DIMENSION A1(4),A2(4),Q(68),Q0300(18),Q0603(16),Q0806(15),
     *          Q1008(19)
      PARAMETER(C40D3=40.0D0/3.0D0)
      PARAMETER(ZERO=0.0D0,ONE=1.0D0,TWO=2.0D0,FOUR=4.0D0)
      PARAMETER(DIVSOR=7.0298584066096D-3,DV1=TWO*DIVSOR,
     *          DV2=FOUR*DIVSOR)
      EQUIVALENCE(Q(1),Q1008(1)),(Q(20),Q0806(1)),(Q(35),Q0603(1)),
     *           (Q(51),Q0300(1))
      DATA Q1008/0.1307592530313668D-01,  0.8580903946529252D-02,
     *           0.1969565496511163D-02, -0.6785908047534174D-04,
     *           0.5157989084218512D-05, -0.3254283578286669D-06,
     *           0.3009165432294816D-07, -0.2853960792992456D-08,
     *           0.3125316663674964D-09, -0.369071915636499D-10,
     *           0.47226104267733D-11,   -0.6452090127750D-12,
     *           0.935240772542D-13,     -0.142880027692D-13,
     *           0.22888743282D-14,      -0.3827976860D-15,
     *           0.665889416D-16,        -0.120097537D-16,
     *           0.22395647D-17/
      DATA Q0806/0.7530091699129213D-01,  0.2215896190278894D-01,
     *           0.1517662771681320D-02, -0.1374204580631322D-04,
     *           0.2501181491115358D-05, -0.2491206397236787D-07,
     *           0.5430948732810670D-08, -0.1406070367942514D-09,
     *           0.160397857131033D-10,  -0.7706914621139D-12,
     *           0.656131517151D-13,     -0.43257167630D-14,
     *           0.3459964512D-15,       -0.264283638D-16,
     *           0.21607480D-17/
      DATA Q0603/0.2285305746758029D+00,  0.5670149715650661D-01,
     *           0.3881611773394649D-02,  0.1483758828946990D-03,
     *           0.1986416462810431D-04,  0.1166284710859293D-05,
     *           0.1048168134503124D-06,  0.6572636994171403D-08,
     *           0.5344588684897204D-09,  0.335115172128537D-10,
     *           0.26225503158527D-11,    0.1606252080762D-12,
     *           0.124685606769D-13,      0.73251000D-15,
     *           0.579829277D-16,         0.31817034D-17/
      DATA Q0300/0.5404466499651320D+00,  0.1034761954005156D+00,
     *           0.9090551216566596D-02,  0.9130257654553327D-03,
     *           0.1026685176623610D-03,  0.1035244509501565D-04,
     *           0.1034925154934048D-05,  0.1002050044392230D-06,
     *           0.9552837592432324D-08,  0.8941334072231475D-09,
     *           0.825727197051832D-10,   0.75255600668576D-11,
     *           0.6783483380205D-12,     0.605164239084D-13,
     *           0.53497536298D-14,       0.4689365198D-15,
     *           0.407909118D-16,         0.35230149D-17/
      DATA A1/20.0D0,20.0D0,C40D3,C40D3/,
     *     A2/18.0D0,14.0D0,6.0D0,2.0D0/,
     *     ARG/1,20,35,51,69/
      DATA P1/-.3D0/,P2/-.6D0/,P3/-.8D0/
C
C     SET THE MACHINE DEPENDENT CONSTANTS.
      UFLOW=2*D1MACH(1)
      EOFLOW=LOG(D1MACH(2)/2)
C
C     'UFLOW'   IS ABOUT EQUAL TO, BUT MUST NOT BE LESS THAN THE
C               SMALLEST POSITIVE REAL NUMBER REPRESENTABLE BY
C               THE MACHINE.
C     'EOFLOW'  IS ABOUT EQUAL TO, BUT MUST NOT EXCEED THE NAPERIAN
C               LOGARITHM OF THE MACHINE'S LARGEST REAL NUMBER.
C
C     Q CONTAINS THE COEFFICIENTS OF THE CHEBYSHEV EXPANSION OF FI.
C     TO ACHIEVE GREATER THAN 16-DIGIT ACCURACY, IT IS NECESSARY TO
C     INCREASE THE NUMBER OF DIGITS IN THESE COEFFICIENTS AND ALSO
C     THE NUMBER OF COEFFICIENTS USED IN THE CHEBYSHEV SERIES. THIS
C     REQUIRES CHANGES TO THE DIMENSION OF Q, Q0300, Q0603, Q0806,
C     AND Q1008, TO THE VALUES OF ARG, AND TO THE CORRESPONDING
C     SUBSCRIPTS OF Q IN THE EQUIVALENCE STATEMENT. REFER TO
C     APPENDIX 3 OF THE FOLLOWING FOR DETAILS:
C         A HAEGEMANS, ALGORITHM 34: AN ALGORITHM FOR THE AUTOMATIC
C         INTEGRATION OVER A TRIANGLE, COMPUTING 19(1977),PP. 179-187.
C
C     SET THE ARGUMENT X FOR EVALUATION.
C
      IF(Z .LE. ZERO) THEN
          X=Z
      ELSE
          X=-Z
      END IF
C
C     COMPUTE SI.
C
      IF(X .LE. -ONE) THEN
          FI=ZERO
          SI=ZERO
          GO TO 20
      ELSE
          EXPONT=FOUR/(ONE-X*X)
          IF(EXPONT .LE. EOFLOW) THEN
              SI=EXP(-EXPONT)
          ELSE
              FI=ZERO
              SI=ZERO
              GO TO 20
          END IF
      END IF
C
C     DETERMINE COEFFICIENTS TO BE USED IN CHEBYSHEV SERIES
C
      IF(X .GE. P1) THEN
          INT=4
      ELSE
          IF(X .GE. P2) THEN
              INT=3
          ELSE
              IF(X .GE. P3) THEN
                  INT=2
              ELSE
                  INT=1
              END IF
          END IF
      END IF
      Y=A1(INT)*X+A2(INT)
      KFIRST=ARG(INT)
      KLAST=ARG(INT+1)-2
C
C     COMPUTE FI.
C
      T2=ZERO
      T3=Q(KLAST+1)
      KSUM=KFIRST+KLAST
      DO 10 K=KFIRST,KLAST
        KREV=KSUM-K
        T1=T2
        T2=T3
        T3=Y*T2-T1+Q(KREV)
10      CONTINUE
      FAC=(T3-T1)/DV2
      IF(FAC .GT. UFLOW/SI) THEN
          FI=FAC*SI
      ELSE
          FI=ZERO
      END IF
20    CONTINUE
      SI=SI/DV1
      IF(Z .GT. ZERO) FI=ONE-FI
      RETURN
      END

      SUBROUTINE ZEROLG(N,ZZ,WW,A,B)
C     -----------------
C
C     PRODUCE THE NODES AND WEIGHTS FOR GAUSS-LEGENDRE QUADRATURE
C     ON (A,B), OF ORDER N. THE NUMBER OF NODES N IS RESTRICTED TO
C     BE .LE. 20.
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION Z(109),W(109),ZZ(N),WW(N)
      DATA ONE/1.0D0/,TWO/2.0D0/
      DATA Z(1),Z(2),Z(3),Z(4),Z(5),Z(6),Z(7),Z(8),Z(9),Z(10)/
     *.577350269189626D0,.774596669241483D0,0.0D0,.861136311594053D0,
     *.339981043584856D0,.906179845938664D0,.538469310105683D0,0.0D0,
     *.932469514203152D0,.661209386466265D0/
      DATA Z(11),Z(12),Z(13),Z(14),Z(15),Z(16),Z(17),Z(18),Z(19),Z(20)/
     *.238619186083197D0,.949107912342759D0,.741531185599394D0,
     *.405845151377397D0,0.0D0,.960289856497536D0,.796666477413627D0,
     *.525532409916329D0,.183434642495650D0,.968160239507626D0/
      DATA Z(21),Z(22),Z(23),Z(24),Z(25),Z(26),Z(27),Z(28),Z(29)/
     *.836031107326636D0,.613371432700590D0,.324253423403809D0,
     *0.0D0,.973906528517172D0,.865063366688985D0,.679409568299024D0,
     *.433395394129247D0,.148874338981631D0/
      DATA Z(30),Z(31),Z(32),Z(33),Z(34),Z(35),Z(36),Z(37),Z(38)/
     *.978228658146057D0,.887062599768095D0,.730152005574049D0,
     *.519096129206812D0,.269543155952345D0,0.0D0,
     *.981560634246719D0,.904117256370475D0,.769902674194305D0/
      DATA Z(39),Z(40),Z(41),Z(42),Z(43),Z(44),Z(45),Z(46),Z(47)/
     *.587317954286617D0,.367831498998180D0,.125233408511469D0,
     *.984183054718588D0,.917598399222978D0,.801578090733310D0,
     *.642349339440340D0,.448492751036447D0,.230458315955135D0/
      DATA Z(48),Z(49),Z(50),Z(51),Z(52),Z(53),Z(54),Z(55),Z(56)/
     *0.0D0,.986283808696812D0,.928434883663574D0,
     *.827201315069765D0,.687292904811685D0,.515248636358154D0,
     *.319112368927890D0,.108054948707344D0,.987992518020485D0/
      DATA Z(57),Z(58),Z(59),Z(60),Z(61),Z(62),Z(63),Z(64),Z(65)/
     *.937273392400706D0,.848206583410427D0,.724417731360170D0,
     *.570972172608539D0,.394151347077563D0,.201194093997435D0,
     *0.0D0,.989400934991650D0,.944575023073233D0/
      DATA Z(66),Z(67),Z(68),Z(69),Z(70),Z(71),Z(72),Z(73),Z(74)/
     *.865631202387832D0,.755404408355003D0,.617876244402644D0,
     *.458016777657227D0,.281603550779259D0,.0950125098376374D0,
     *.990575475314417D0,.950675521768768D0,.880239153726986D0/
      DATA Z(75),Z(76),Z(77),Z(78),Z(79),Z(80),Z(81),Z(82),Z(83)/
     *.781514003896801D0,.657671159216691D0,.512690537086477D0,
     *.351231763453876D0,.178484181495848D0,0.0D0,
     *.991565168420931D0,.955823949571398D0,.892602466497556D0/
      DATA Z(84),Z(85),Z(86),Z(87),Z(88),Z(89),Z(90),Z(91),Z(92)/
     *.803704958972523D0,.691687043060353D0,.559770831073948D0,
     *.411751161462843D0,.251886225691506D0,.0847750130417353D0,
     *.992406843843584D0,.960208152134830D0,.903155903614818D0/
      DATA Z(93),Z(94),Z(95),Z(96),Z(97),Z(98),Z(99),Z(100),Z(101)/
     *.822714656537143D0,.720966177335229D0,.600545304661681D0,
     *.464570741375961D0,.316564099963630D0,.160358645640225D0,
     *0.0D0,.993128599185095D0,.963971927277914D0/
      DATA Z(102),Z(103),Z(104),Z(105),Z(106),Z(107),Z(108),Z(109)/
     *.912234428251326D0,.839116971822219D0,.746331906460151D0,
     *.636053680726515D0,.510867001950827D0,.373706088715420D0,
     *.227785851141645D0,.0765265211334973D0/
      DATA W(1),W(2),W(3),W(4),W(5),W(6),W(7),W(8),W(9),W(10)/
     *1.0D0,.555555555555556D0,.888888888888889D0,.347854845137454D0,
     *.652145154862546D0,.236926885056189D0,.478628670499366D0,
     *.568888888888889D0,.171324492379170D0,.360761573048139D0/
      DATA W(11),W(12),W(13),W(14),W(15),W(16),W(17),W(18),W(19),W(20)/
     *.467913934572691D0,.129484966168870D0,.279705391489277D0,
     *.381830050505119D0,.417959183673469D0,.101228536290376D0,
     *.222381034453374D0,.313706645877887D0,.362683783378362D0,
     *.0812743883615744D0/
      DATA W(21),W(22),W(23),W(24),W(25),W(26),W(27),W(28),W(29)/
     *.180648160694857D0,.260610696402935D0,.312347077040003D0,
     *.330239355001260D0,.0666713443086881D0,.149451349150581D0,
     *.219086362515982D0,.269266719309996D0,.295524224714753D0/
      DATA W(30),W(31),W(32),W(33),W(34),W(35),W(36),W(37),W(38)/
     *.0556685671161737D0,.125580369464905D0,.186290210927734D0,
     *.233193764591990D0,.262804544510247D0,.272925086777901D0,
     *.0471753363865118D0,.106939325995318D0,.160078328543346D0/
      DATA W(39),W(40),W(41),W(42),W(43),W(44),W(45),W(46),W(47)/
     *.203167426723066D0,.233492536538355D0,.249147045813403D0,
     *.0404840047653159D0,.0921214998377284D0,.138873510219787D0,
     *.178145980761946D0,.207816047536889D0,.226283180262897D0/
      DATA W(48),W(49),W(50),W(51),W(52),W(53),W(54),W(55),W(56)/
     *.232551553230874D0,.0351194603317519D0,.0801580871597602D0,
     *.121518570687903D0,.157203167158194D0,.185538397477938D0,
     *.205198463721296D0,.215263853463158D0,.0307532419961173D0/
      DATA W(57),W(58),W(59),W(60),W(61),W(62),W(63),W(64),W(65)/
     *.0703660474881081D0,.107159220467172D0,.139570677926154D0,
     *.166269205816994D0,.186161000015562D0,.198431485327112D0,
     *.202578241925561D0,.0271524594117541D0,.0622535239386478D0/
      DATA W(66),W(67),W(68),W(69),W(70),W(71),W(72),W(73),W(74)/
     *.0951585116824928D0,.124628971255534D0,.149595988816577D0,
     *.169156519395003D0,.182603415044924D0,.189450610455068D0,
     *.0241483028685479D0,.0554595293739872D0,.0850361483171792D0/
      DATA W(75),W(76),W(77),W(78),W(79),W(80),W(81),W(82),W(83)/
     *.111883847193404D0,.135136368468525D0,.154045761076810D0,
     *.168004102156450D0,.176562705366993D0,.179446470356207D0,
     *.0216160135264833D0,.0497145488949698D0,.0764257302548891D0/
      DATA W(84),W(85),W(86),W(87),W(88),W(89),W(90),W(91),W(92)/
     *.100942044106287D0,.122555206711478D0,.140642914670651D0,
     *.154684675126265D0,.164276483745833D0,.169142382963144D0,
     *.0194617882297265D0,.0448142267656996D0,.0690445427376412D0/
      DATA W(93),W(94),W(95),W(96),W(97),W(98),W(99),W(100),W(101)/
     *.0914900216224500D0,.111566645547334D0,.128753962539336D0,
     *.142606702173607D0,.152766042065860D0,.158968843393954D0,
     *.161054449848784D0,.0176140071391521D0,.0406014298003869D0/
      DATA W(102),W(103),W(104),W(105),W(106),W(107),W(108),W(109)/
     *.0626720483341091D0,.0832767415767047D0,.101930119817240D0,
     *.118194531961518D0,.131688638449177D0,.142096109318382D0,
     *.149172986472604D0,.152753387130726D0/
C
C     CALCULATE THE WEIGHTS AND NODES.
C
      SCALE=(B-A)/TWO
      IF((N/2)*2 .EQ. N) THEN
          M=N/2
          IBASE=M*M
      ELSE
          IBASE=(N*N-1)/4
          M=(N-1)/2
          ZZ(M+1)=(A+B)/TWO
          WW(M+1)=W(IBASE+M)*SCALE
      END IF
      DO 100 I=1,M
        T=Z(IBASE+I-1)
        ZZ(I)=(A*(ONE+T)+(ONE-T)*B)/TWO
        ZZ(N+1-I)=(A*(ONE-T)+(ONE+T)*B)/TWO
        WW(I)=W(IBASE+I-1)*SCALE
100     WW(N+1-I)=WW(I)
      RETURN
      END

      DOUBLE PRECISION FUNCTION D1MACH (I) 
C     --------------------------------
C
C***BEGIN PROLOGUE  D1MACH
C***DATE WRITTEN   750101   (YYMMDD)
C***REVISION DATE  860501   (YYMMDD)
C***CATEGORY NO.  R1
C***KEYWORDS  MACHINE CONSTANTS
C***AUTHOR  FOX, P. A., (BELL LABS)
C           HALL, A. D., (BELL LABS)
C           SCHRYER, N. L., (BELL LABS)
C***PURPOSE  RETURN DOUBLE PRECISION MACHINE DEPENDENT CONSTANTS.
C***DESCRIPTION
C
C     D1MACH CAN BE USED TO OBTAIN MACHINE-DEPENDENT PARAMETERS
C     FOR THE LOCAL MACHINE ENVIRONMENT.  IT IS A FUNCTION
C     SUBPROGRAM WITH ONE (INPUT) ARGUMENT, AND CAN BE CALLED
C     AS FOLLOWS, FOR EXAMPLE
C
C          D = D1MACH(I)
C
C     WHERE I=1,...,5.  THE (OUTPUT) VALUE OF D ABOVE IS
C     DETERMINED BY THE (INPUT) VALUE OF I.  THE RESULTS FOR
C     VARIOUS VALUES OF I ARE DISCUSSED BELOW.
C
C  DOUBLE-PRECISION MACHINE CONSTANTS
C  D1MACH( 1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE.
C  D1MACH( 2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE.
C  D1MACH( 3) = B**(-T), THE SMALLEST RELATIVE SPACING.
C  D1MACH( 4) = B**(1-T), THE LARGEST RELATIVE SPACING.
C  D1MACH( 5) = LOG10(B)
C***REFERENCES  FOX P.A., HALL A.D., SCHRYER N.L.,*FRAMEWORK FOR A
C                 PORTABLE LIBRARY*, ACM TRANSACTIONS ON MATHEMATICAL
C                 SOFTWARE, VOL. 4, NO. 2, JUNE 1978, PP. 177-188.
C***ROUTINES CALLED  XERROR
C***END PROLOGUE  D1MACH
C
      INTEGER SMALL(4)
      INTEGER LARGE(4)
      INTEGER RIGHT(4)
      INTEGER DIVER(4)
      INTEGER LOG10(4)
C
      DOUBLE PRECISION DMACH(5)
C
      EQUIVALENCE (DMACH(1),SMALL(1))
      EQUIVALENCE (DMACH(2),LARGE(1))
      EQUIVALENCE (DMACH(3),RIGHT(1))
      EQUIVALENCE (DMACH(4),DIVER(1))
      EQUIVALENCE (DMACH(5),LOG10(1))
C
C     MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM.
C
C     DATA SMALL(1) / ZC00800000 /
C     DATA SMALL(2) / Z000000000 /
C
C     DATA LARGE(1) / ZDFFFFFFFF /
C     DATA LARGE(2) / ZFFFFFFFFF /
C
C     DATA RIGHT(1) / ZCC5800000 /
C     DATA RIGHT(2) / Z000000000 /
C
C     DATA DIVER(1) / ZCC6800000 /
C     DATA DIVER(2) / Z000000000 /
C
C     DATA LOG10(1) / ZD00E730E7 /
C     DATA LOG10(2) / ZC77800DC0 /
C
C     MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM.
C
C     DATA SMALL(1) / O1771000000000000 /
C     DATA SMALL(2) / O0000000000000000 /
C
C     DATA LARGE(1) / O0777777777777777 /
C     DATA LARGE(2) / O0007777777777777 /
C
C     DATA RIGHT(1) / O1461000000000000 /
C     DATA RIGHT(2) / O0000000000000000 /
C
C     DATA DIVER(1) / O1451000000000000 /
C     DATA DIVER(2) / O0000000000000000 /
C
C     DATA LOG10(1) / O1157163034761674 /
C     DATA LOG10(2) / O0006677466732724 /
C
C     MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS.
C
C     DATA SMALL(1) / O1771000000000000 /
C     DATA SMALL(2) / O7770000000000000 /
C
C     DATA LARGE(1) / O0777777777777777 /
C     DATA LARGE(2) / O7777777777777777 /
C
C     DATA RIGHT(1) / O1461000000000000 /
C     DATA RIGHT(2) / O0000000000000000 /
C
C     DATA DIVER(1) / O1451000000000000 /
C     DATA DIVER(2) / O0000000000000000 /
C
C     DATA LOG10(1) / O1157163034761674 /
C     DATA LOG10(2) / O0006677466732724 /
C
C     MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES.
C     FOR FTN4
C
C     DATA SMALL(1) / 00564000000000000000B /
C     DATA SMALL(2) / 00000000000000000000B /
C
C     DATA LARGE(1) / 37757777777777777777B /
C     DATA LARGE(2) / 37157777777777777777B /
C
C     DATA RIGHT(1) / 15624000000000000000B /
C     DATA RIGHT(2) / 00000000000000000000B /
C
C     DATA DIVER(1) / 15634000000000000000B /
C     DATA DIVER(2) / 00000000000000000000B /
C
C     DATA LOG10(1) / 17164642023241175717B /
C     DATA LOG10(2) / 16367571421742254654B /
C
C     MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES.
C     FOR FTN5
C
C     DATA SMALL(1) / O"00564000000000000000" /
C     DATA SMALL(2) / O"00000000000000000000" /
C
C     DATA LARGE(1) / O"37757777777777777777" /
C     DATA LARGE(2) / O"37157777777777777777" /
C
C     DATA RIGHT(1) / O"15624000000000000000" /
C     DATA RIGHT(2) / O"00000000000000000000" /
C
C     DATA DIVER(1) / O"15634000000000000000" /
C     DATA DIVER(2) / O"00000000000000000000" /
C
C     DATA LOG10(1) / O"17164642023241175717" /
C     DATA LOG10(2) / O"16367571421742254654" /
C
C     MACHINE CONSTANTS FOR THE CRAY 1
C
C     DATA SMALL(1) / 201354000000000000000B /
C     DATA SMALL(2) / 000000000000000000000B /
C
C     DATA LARGE(1) / 577767777777777777777B /
C     DATA LARGE(2) / 000007777777777777774B /
C
C     DATA RIGHT(1) / 376434000000000000000B /
C     DATA RIGHT(2) / 000000000000000000000B /
C
C     DATA DIVER(1) / 376444000000000000000B /
C     DATA DIVER(2) / 000000000000000000000B /
C
C     DATA LOG10(1) / 377774642023241175717B /
C     DATA LOG10(2) / 000007571421742254654B /
C
C     MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200
C
C     NOTE - IT MAY BE APPROPRIATE TO INCLUDE THE FOLLOWING CARD -
C     STATIC DMACH(5)
C
C     DATA SMALL/20K,3*0/,LARGE/77777K,3*177777K/
C     DATA RIGHT/31420K,3*0/,DIVER/32020K,3*0/
C     DATA LOG10/40423K,42023K,50237K,74776K/
C
C     MACHINE CONSTANTS FOR THE HARRIS 220
C
C     DATA SMALL(1), SMALL(2) / '20000000, '00000201 /
C     DATA LARGE(1), LARGE(2) / '37777777, '37777577 /
C     DATA RIGHT(1), RIGHT(2) / '20000000, '00000333 /
C     DATA DIVER(1), DIVER(2) / '20000000, '00000334 /
C     DATA LOG10(1), LOG10(2) / '23210115, '10237777 /
C
C     MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 SERIES.
C
C     DATA SMALL(1), SMALL(2) / O402400000000, O000000000000 /
C     DATA LARGE(1), LARGE(2) / O376777777777, O777777777777 /
C     DATA RIGHT(1), RIGHT(2) / O604400000000, O000000000000 /
C     DATA DIVER(1), DIVER(2) / O606400000000, O000000000000 /
C     DATA LOG10(1), LOG10(2) / O776464202324, O117571775714 /
C
C     MACHINE CONSTANTS FOR THE HP 2100
C     THREE WORD DOUBLE PRECISION OPTION WITH FTN4
C
C     DATA SMALL(1), SMALL(2), SMALL(3) / 40000B,       0,       1 /
C     DATA LARGE(1), LARGE(2), LARGE(3) / 77777B, 177777B, 177776B /
C     DATA RIGHT(1), RIGHT(2), RIGHT(3) / 40000B,       0,    265B /
C     DATA DIVER(1), DIVER(2), DIVER(3) / 40000B,       0,    276B /
C     DATA LOG10(1), LOG10(2), LOG10(3) / 46420B,  46502B,  77777B /
C
C     MACHINE CONSTANTS FOR THE HP 2100
C     FOUR WORD DOUBLE PRECISION OPTION WITH FTN4
C
C     DATA SMALL(1), SMALL(2) /  40000B,       0 /
C     DATA SMALL(3), SMALL(4) /       0,       1 /
C     DATA LARGE(1), LARGE(2) /  77777B, 177777B /
C     DATA LARGE(3), LARGE(4) / 177777B, 177776B /
C     DATA RIGHT(1), RIGHT(2) /  40000B,       0 /
C     DATA RIGHT(3), RIGHT(4) /       0,    225B /
C     DATA DIVER(1), DIVER(2) /  40000B,       0 /
C     DATA DIVER(3), DIVER(4) /       0,    227B /
C     DATA LOG10(1), LOG10(2) /  46420B,  46502B /
C     DATA LOG10(3), LOG10(4) /  76747B, 176377B /
C
C     MACHINE CONSTANTS FOR THE HP 9000
C
C     D1MACH(1) = 2.8480954D-306
C     D1MACH(2) = 1.40444776D+306
C     D1MACH(3) = 2.22044605D-16
C     D1MACH(4) = 4.44089210D-16
C     D1MACH(5) = 3.01029996D-1
C
C     DATA SMALL(1), SMALL(2) / 00040000000B, 00000000000B /
C     DATA LARGE(1), LARGE(2) / 17737777777B, 37777777777B /
C     DATA RIGHT(1), RIGHT(2) / 07454000000B, 00000000000B /
C     DATA DIVER(1), DIVER(2) / 07460000000B, 00000000000B /
C     DATA LOG10(1), LOG10(2) / 07764642023B, 12047674777B /
C
C     MACHINE CONSTANTS FOR THE IBM 360/370 SERIES,
C     THE XEROX SIGMA 5/7/9, THE SEL SYSTEMS 85/86, AND
C     THE PERKIN ELMER (INTERDATA) 7/32.
C
C     DATA SMALL(1), SMALL(2) / Z00100000, Z00000000 /
C     DATA LARGE(1), LARGE(2) / Z7FFFFFFF, ZFFFFFFFF /
C     DATA RIGHT(1), RIGHT(2) / Z33100000, Z00000000 /
C     DATA DIVER(1), DIVER(2) / Z34100000, Z00000000 /
C     DATA LOG10(1), LOG10(2) / Z41134413, Z509F79FF /
C
C     MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR).
C
C     DATA SMALL(1), SMALL(2) / "033400000000, "000000000000 /
C     DATA LARGE(1), LARGE(2) / "377777777777, "344777777777 /
C     DATA RIGHT(1), RIGHT(2) / "113400000000, "000000000000 /
C     DATA DIVER(1), DIVER(2) / "114400000000, "000000000000 /
C     DATA LOG10(1), LOG10(2) / "177464202324, "144117571776 /
C
C     MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR).
C
C     DATA SMALL(1), SMALL(2) / "000400000000, "000000000000 /
C     DATA LARGE(1), LARGE(2) / "377777777777, "377777777777 /
C     DATA RIGHT(1), RIGHT(2) / "103400000000, "000000000000 /
C     DATA DIVER(1), DIVER(2) / "104400000000, "000000000000 /
C     DATA LOG10(1), LOG10(2) / "177464202324, "476747767461 /
C
C     MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING
C     32-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL).
C
C     DATA SMALL(1), SMALL(2) /    8388608,           0 /
C     DATA LARGE(1), LARGE(2) / 2147483647,          -1 /
C     DATA RIGHT(1), RIGHT(2) /  612368384,           0 /
C     DATA DIVER(1), DIVER(2) /  620756992,           0 /
C     DATA LOG10(1), LOG10(2) / 1067065498, -2063872008 /
C
C     DATA SMALL(1), SMALL(2) / O00040000000, O00000000000 /
C     DATA LARGE(1), LARGE(2) / O17777777777, O37777777777 /
C     DATA RIGHT(1), RIGHT(2) / O04440000000, O00000000000 /
C     DATA DIVER(1), DIVER(2) / O04500000000, O00000000000 /
C     DATA LOG10(1), LOG10(2) / O07746420232, O20476747770 /
C
C     MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING
C     16-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL).
C
C     DATA SMALL(1), SMALL(2) /    128,      0 /
C     DATA SMALL(3), SMALL(4) /      0,      0 /
C
C     DATA LARGE(1), LARGE(2) /  32767,     -1 /
C     DATA LARGE(3), LARGE(4) /     -1,     -1 /
C
C     DATA RIGHT(1), RIGHT(2) /   9344,      0 /
C     DATA RIGHT(3), RIGHT(4) /      0,      0 /
C
C     DATA DIVER(1), DIVER(2) /   9472,      0 /
C     DATA DIVER(3), DIVER(4) /      0,      0 /
C
C     DATA LOG10(1), LOG10(2) /  16282,   8346 /
C     DATA LOG10(3), LOG10(4) / -31493, -12296 /
C
C     DATA SMALL(1), SMALL(2) / O000200, O000000 /
C     DATA SMALL(3), SMALL(4) / O000000, O000000 /
C
C     DATA LARGE(1), LARGE(2) / O077777, O177777 /
C     DATA LARGE(3), LARGE(4) / O177777, O177777 /
C
C     DATA RIGHT(1), RIGHT(2) / O022200, O000000 /
C     DATA RIGHT(3), RIGHT(4) / O000000, O000000 /
C
C     DATA DIVER(1), DIVER(2) / O022400, O000000 /
C     DATA DIVER(3), DIVER(4) / O000000, O000000 /
C
C     DATA LOG10(1), LOG10(2) / O037632, O020232 /
C     DATA LOG10(3), LOG10(4) / O102373, O147770 /
C
C     MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. FTN COMPILER
C
C     DATA SMALL(1), SMALL(2) / O000040000000, O000000000000 /
C     DATA LARGE(1), LARGE(2) / O377777777777, O777777777777 /
C     DATA RIGHT(1), RIGHT(2) / O170540000000, O000000000000 /
C     DATA DIVER(1), DIVER(2) / O170640000000, O000000000000 /
C     DATA LOG10(1), LOG10(2) / O177746420232, O411757177572 /
C
C     MACHINE CONSTANTS FOR VAX 11/780
C     (EXPRESSED IN INTEGER AND HEXADECIMAL)
C     ***THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSYEMS***
C     *** THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS***
C
C     DATA SMALL(1), SMALL(2) /        128,           0 /
C     DATA LARGE(1), LARGE(2) /     -32769,          -1 /
C     DATA RIGHT(1), RIGHT(2) /       9344,           0 /
C     DATA DIVER(1), DIVER(2) /       9472,           0 /
C     DATA LOG10(1), LOG10(2) /  546979738,  -805796613 /
C
C     DATA SMALL(1), SMALL(2) / Z00000080, Z00000000 /
C     DATA LARGE(1), LARGE(2) / ZFFFF7FFF, ZFFFFFFFF /
C     DATA RIGHT(1), RIGHT(2) / Z00002480, Z00000000 /
C     DATA DIVER(1), DIVER(2) / Z00002500, Z00000000 /
C     DATA LOG10(1), LOG10(2) / Z209A3F9A, ZCFF884FB /
C
C     MACHINE CONSTANTS FOR VAX 11/780 (G-FLOATING)
C     (EXPRESSED IN INTEGER AND HEXADECIMAL)
C     ***THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSYEMS***
C     *** THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS***
C
C     DATA SMALL(1), SMALL(2) /         16,           0 /
C     DATA LARGE(1), LARGE(2) /     -32769,          -1 /
C     DATA RIGHT(1), RIGHT(2) /      15552,           0 /
C     DATA DIVER(1), DIVER(2) /      15568,           0 /
C     DATA LOG10(1), LOG10(2) /  1142112243, 2046775455 /
C
C     DATA SMALL(1), SMALL(2) / Z00000010, Z00000000 /
C     DATA LARGE(1), LARGE(2) / ZFFFF7FFF, ZFFFFFFFF /
C     DATA RIGHT(1), RIGHT(2) / Z00003CC0, Z00000000 /
C     DATA DIVER(1), DIVER(2) / Z00003CD0, Z00000000 /
C     DATA LOG10(1), LOG10(2) / Z44133FF3, Z79FF509F /
C
C     MACHINE CONSTANTS FOR THE ELXSI 6400
C     (ASSUMING REAL*8 IS THE DEFAULT DOUBLE PRECISION)
C
C     DATA SMALL(1), SMALL(2) / '00100000'X,'00000000'X /
C     DATA LARGE(1), LARGE(2) / '7FEFFFFF'X,'FFFFFFFF'X /
C     DATA RIGHT(1), RIGHT(2) / '3CB00000'X,'00000000'X /
C     DATA DIVER(1), DIVER(2) / '3CC00000'X,'00000000'X /
C     DATA LOG10(1), DIVER(2) / '3FD34413'X,'509F79FF'X /
C
C     MACHINE CONSTANTS FOR THE IBM PC - MICROSOFT FORTRAN
C
C     DATA SMALL(1), SMALL(2) / #00000000, #00100000 /
C     DATA LARGE(1), LARGE(2) / #FFFFFFFF, #7FEFFFFF /
C     DATA RIGHT(1), RIGHT(2) / #00000000, #3CA00000 /
C     DATA DIVER(1), DIVER(2) / #00000000, #3CB00000 /
C     DATA LOG10(1), LOG10(2) / #509F79FF, #3FD34413 /
C
C     MACHINE CONSTANTS FOR THE IBM PC - PROFESSIONAL FORTRAN
C     AND LAHEY FORTRAN
C
C     DATA SMALL(1), SMALL(2) / Z'00000000', Z'00100000' /
C     DATA LARGE(1), LARGE(2) / Z'FFFFFFFF', Z'7FEFFFFF' /
C     DATA RIGHT(1), RIGHT(2) / Z'00000000', Z'3CA00000' /
C     DATA DIVER(1), DIVER(2) / Z'00000000', Z'3CB00000' /
C     DATA LOG10(1), LOG10(2) / Z'509F79FF', Z'3FD34413' /
C
C     MACHINE CONSTANTS FOR THE PRIME 750/850.   {9/23/86 TSENG}
C
C     DATA SMALL(1), SMALL(2) / :10000000000, :00000100000 /
C     DATA LARGE(1), LARGE(2) / :17777777777, :37777677602 /
C     DATA RIGHT(1), RIGHT(2) / :10000000000, :00000000122 /
C     DATA DIVER(1), DIVER(2) / :10000000000, :00000000123 /
C     DATA LOG10(1), LOG10(2) / :11504046502, :17572000177 /
C
C     MACHINE CONSTANTS FOR THE APOLLO.   {9/23/86 TSENG}
C
      DATA SMALL(1), SMALL(2) / 8#00004000000, 8#00000000000 /
      DATA LARGE(1), LARGE(2) / 8#17773777777, 8#37777777777 /
      DATA RIGHT(1), RIGHT(2) / 8#07450000000, 8#00000000000 /
      DATA DIVER(1), DIVER(2) / 8#07454000000, 8#00000000000 /
      DATA LOG10(1), LOG10(2) / 8#07764642023, 8#12047674777 /
C
C
C***FIRST EXECUTABLE STATEMENT  D1MACH
C     IF (I .LT. 1  .OR.  I .GT. 5)
C    1   CALL XERROR( 'D1MACH -- I OUT OF BOUNDS',25,1,2)
C
      D1MACH = DMACH(I)
      RETURN
C
      END

