      PROGRAM FPRIME
C     PROGRAM TO INTERPOLATE CROSS SECTIONS
C     AND COMPUTE ANOMALOUS SCATTERING FACTORS
C
C
C     QUESTIONS CONCERNING THIS PROGRAM SHOULD BE ADDRESSED TO
C            DON T. CROMER   MST-5
C            MAIL STOP G730
C            LOS ALAMOS NATIONAL LABORATORY
C            LOS ALAMOS, NEW MEXICO 87545
C
C
C
C     Revised for VAX FORTRAN by R.B. Von Dreele
C
C
C     NORD=0 INTERPOLATION IS BY FITTING THREE CLOSEST POINTS
C       TO A QUADRATIC
C     NORD=1,2,3 ETC. USE AITKEN'S INTERPOLATION METHOD AND NORD
C      IS THE INTERPOLATION ORDER. NORD=2 IS PROBABLY THE BEST
C       VALUE TO USE
C
C
C          *********NOTE********
C
C
C     THREE CONTRIBUTIONS TO FPRIME ARE LISTED SEPARATELY
C     IN THE OUTPUT.
C         1.  P.E.      PHOTOELECTRIC CONTRIBUTION
C         2. ETERM      TERM DEPENDENT ON THE TOTAL ENERGY OF
C                       THE ATOM. CROMER AND LIBERMAN J. CHEM. PHYS.
C                       53,1891-1898(1970).
C         3. JENSEN     TERM DEPENDENT ON THE XRAY ENERGY.
C                       JENSEN, PHYSICS LETTERS,74A,41-44(1979)
C
C
C
C
C     THE CROSS SECTION FILE WILL HAVE A NUMBER OF ORBITALS
C     FOR EACH ATOM FROM ATOMIC NUMBER 3-98
C     THE FIRST MX CARDS (MX=5 FOR THIS XSECTION FILE AND IS SET IN PROGRAM)
C     FOR EACH ORBITAL WILL HAVE CROSS SECTIONS AT MX VALUES OF ENERGY
C     FROM ABOUT 1 TO 80 KEV APPROXIMATELY EQUALLY SPACED IN LOG(ENERGY)
C     THE NEXT FIVE RECORDS WILL BE CROSS SECTIONS AT ENERGIES SELECTED
C     BY GAUSS INTEGRATION SCHEME. IF THE FUNCTION TYPE IS 0 (IFBE=0)
C     ( FUNCTION TYPE GIVEN IN CROMER AND LIBERMAN
C         J. CHEM. PHYS. 53,1891-1898(1970)     )
C     A SIXTH VALUE IS READ IN FOR AN ENERGY=1.001*(BINDING ENERGY).
C     IF THE XRAY ENERGY IS LESS THAN BINDING ENERGY, FUNCTION SIGMA3
C     WILL BE USED (CROMER AND LIBERMAN ACTA CRYST. A37,267-268(1981))
C
C     FILES IN IO AND IS ARE SET AT THE BEGINNING OF THE PROGRAM
C     THESE WILL DEPEND ON LOCAL CONVENTIONS
C
C***********WARNING*************
C
C
C     IF AN XRAY ENERGY IS VERY CLOSE TO ONE OF THE ENERGIES
C     USED IN THE GAUSS INTEGRATION AN "ANOMALOUS" ANOMALOUS SCATTERING
C     FACTOR MAY RESULT. THERE IS NO EASY WAY OUT OF THIS PROBLEM. A
C     SUGGESTED WAY IS TO COMPUTE SEVERAL VALUES AT NEARBY ENERGIES
C     AND DRAW A SMOOTH CURVE. THIS METHOD SHOULD WORK PROVIDED
C     THE POINTS DO NOT PASS THROUGH AN EDGE
C
C
C
        IMPLICIT REAL*8 (A-H,O-Z),INTEGER*4 (I-N)
C
        CHARACTER*80     EVNAME,AWAVE
        CHARACTER*2     NATOM,NAT
        CHARACTER*6     NSHEL(24)
        EXTERNAL SIGMA0, SIGMA1, SIGMA2, SIGMA3, DLGNDR
        DIMENSION XK(10), XR(10), FP(10,24), FPP(10,24), ABS(10,24)
        DIMENSION EW(14), SIG(14), EL(14), SL(14)
        DIMENSION CXB(10), SUMFP(10), SUMFPP(10)
        DIMENSION T(32)
        DIMENSION XJENSN(10),ET(10)
        DIMENSION EG(5)
        DIMENSION ETERM(98),WT(98),NORB(98)
      CHARACTER*115 AT
        COMMON /GAUS/ CX, BB, SIGG(5), RX, ICOUNT
        COMMON /EDGE/ SEDGE
        DATA ETERM/0.0,0.0,0.001,0.001,0.002,0.003,0.005,0.007,0.009,
     1  0.011,0.014,0.018,0.021,0.026,0.030,0.035,0.041,0.047,0.053,
     1  0.060,0.068,0.075,0.084,0.093,0.102,0.113,0.123,0.135,0.146,
     1  0.159,0.172,0.186,0.200,0.215,0.231,0.247,0.264,0.282,0.300,
     1  0.319,0.338,0.359,0.380,0.401,0.424,0.447,0.471,0.496,0.521,
     1  0.547,0.575,0.602,0.631,0.660,0.690,0.721,0.753,0.786,0.819,
     1  0.854,0.889,0.925,0.962,1.0,1.039,1.079,1.119,1.161,1.204,
     1  1.248,1.293,1.338,1.385,1.433,1.482,1.532,1.583,1.636,1.689,
     1  1.743,1.799,1.856,1.914,1.973,2.033,2.095,2.157,2.221,2.287,
     1  2.353,2.421,2.490,2.561,2.633,2.707,2.782,2.858,2.936/
        DATA WT/1.0079,4.00260,6.941,9.01218,10.81,12.011,14.0067,
     1  15.9994,18.9984,20.179,22.98977,24.305,26.98154,28.0855,
     2  30.97376,32.06,35.453,39.948,39.0983,40.08,44.9559,47.88,
     1  50.9415,
     3  51.996,54.9380,55.847,58.9332,58.69,63.546,65.38,69.72,72.59,
     4  74.9216,78.96,79.904,83.80,85.4678,87.62,88.9059,91.22,92.9064,
     5  95.94,98.,101.07,102.9055,106.42,107.868,112.41,114.82,118.69,
     6  121.75,127.60,126.9045,131.29,132.9054,137.33,138.9055,140.12,
     7  140.9077,144.24,145.,150.36,151.96,157.25,158.9254,162.50,
     8  164.9304,167.26,168.9342,173.04,174.67,178.49,180.9479,183.85,
     9  186.207,190.2,192.22,195.08,196.9655,200.59,204.383,207.2,
     1  208.9804,209.,210.,222.,223.,226.0254,227.0278,232.0381,
     2  231.0359,238.0482,237.0482,239.,243.,247.,247.,251./
        DATA NORB/0,0,2,2,3, 3,4,4,4,4, 4,4,5,6,7, 7,7,7,7,7,
     1  7,7,8,9,9, 9,9,9,9,9, 9,9,9,9,9, 9,9,12,12,13,
     1  13,14,14,14,14, 14,14,14,14,14, 14,14,14,14,17,
     1  17,17,18,18,18, 18,18,18,19,19, 19,19,19,19,19,
     1  19,20,20,20,21, 21,21,21,21,22, 22,23,23,24,24,
     1  24,24,24,24,24, 24,24,24,24,24, 24,24,24/
        DATA C /137.0367/
        DATA PI /3.14159265/
        DATA AU /2.80022E+7/
        DATA C1 /.02721/
        DATA XKEV /12.397639/
        DATA MX/5/
        DATA EVNAME/'ATMXSECT.DAT'/

C**********  SET INPUT, OUTPUT AND X SECTION FILES*****
        IS=1
        OPEN(IS,FILE=EVNAME,ACCESS='DIRECT',STATUS='OLD',
     1       FORM='FORMATTED',RECL=120)
C**********
      IF (NARGS().EQ.3) THEN
          NW = 1
          CALL GETARG(1,NATOM,LC)
          CALL GETARG(2,AWAVE,LCW)
          READ(AWAVE(1:LCW),'(F10.0)') XR(1)
          NORD = 2               
C          XR(1)
      ELSE
        PRINT '(A)',' ENTER ATOMIC SYMBOL (LEFT JUSTIFIED - ALL CAPS)'
        READ '(A)',NATOM
        PRINT '(A)',' ENTER NUMBER OF WAVELENGTHS OR ENERGIES'
        PRINT '(A)','  TO CALCULATE (<=10) & INTERPOLATION METHOD (0-3)'
        READ *, NW,NORD
        IF ( NW.GT.10 ) NW = 10
        IF ( NW.LT.1 ) NW = 1
        IF ( NORD.GT.3 .OR. NORD.LT.0 ) NORD = 0

C*****READ NUMBER OF XRAY VALUES TO BE INTERPOLATED,
C*****READ NW VALUES OF XRAY(ANGSTROMS) AT WHICH INTERPOLATION TO BE DONE
C***** AND CHANGE ANGSTROMS TO KEV
        DO I=1,NW
          PRINT '(A,I1)',' ENTER WAVELENGTH (a) OR ENERGY (Ev) NO.',I
          READ '(F10.0)',XR(I)
        END DO
      ENDIF
        DO K=1,NW
          IF ( XR(K).LT.100.0 ) THEN
            XK(K)=XKEV/XR(K)
          ELSE
            XK(K)=XR(K)/1000.0
            XR(K) = XKEV/XK(K)
          END IF
        END DO
C        NAT = ' '
        IREC = 1
        DO WHILE ( NATOM.NE.NAT )
          READ(IS,REC=IREC,FMT='(A2,I3,A115)',IOSTAT=IERR) NAT,IZ,AT
          IREC = IREC+3
        END DO
        IF ( IERR.NE.0 ) THEN
          PRINT '(A)',' ****  error - ELEMENT NOT FOUND ****'
          STOP
        END IF
        IREC = IREC-3

        NO=NORB(IZ)
        DO I=1,NW
          CXB(I)=0.0
          DO J=1,NO
            FP(I,J)=0.0
            FPP(I,J)=0.0
          END DO
        END DO

        DO J=1,NO
          READ(IS,REC=IREC,FMT=230) NJ,NSHEL(J),IFBE,BE,
     1      (EW(K),SIG(K),K=1,3)
230       FORMAT(5X,I3,A6,I1,E15.8,1P6E15.8)
          READ(IS,FMT=231,REC=IREC+1) (EW(K),SIG(K),K=4,7)
          READ(IS,FMT=231,REC=IREC+2) (EW(K),SIG(K),K=8,11)
231       FORMAT(1P8E15.8)
          IREC = IREC+3
          DO K=1,5
            EG(K) = EW(K+MX)
            SIGG(K) = SIG(K+MX)
          END DO
          NX = MX+5
          IF ( IFBE.EQ.0 ) THEN
            NX = NX+1
            EDG = EW(NX)
            SIGEDG = SIG(NX)
          ELSE
            EDG = 0.0
            SIGEDG = 0.0
          END IF
          BB=BE/C1
          DO I=1,5
            SIGG(I)=SIGG(I)/AU
          END DO             
C*****SORT ALL X SECTIONS
          CALL FPSORT(NX,EW,SIG)
C*****SORT THE FIVE CROSS SECTIONS AT INTEGRATION POINTS
          CALL FPSORT(5,EG,SIGG)
C*****CHANGE TO LOGS FOR INTERPOLATION

          DO K=1,NX
            EL(K)=0.
            SL(K)=0.
            EL(K)=LOG(EW(K))
            IF (SIG(K).GT.0.0) SL(K)=LOG(SIG(K))
          END DO

          DO K=1,NW
            MF=0
            ZX=LOG(XK(K))                       !LOG OF XRAY(KEV) ENERGY
            CX=0.
            IF (BE.LE.XK(K)) THEN
              IF (NORD.EQ.0) THEN
                CALL XSECT (ZX,EL,SL,CX,NX)
              ELSE

                DO M=1,NX
                  N1=M
                  IF (SL(M).NE.0.0) GO TO 100
                END DO
100             MM=NX-N1+1

                CX=AKNINT(ZX,MM,NORD,EL(N1),SL(N1),T)
                CX=EXP(CX)                      !CX IS THE INTERPOLATED X SECTION IN BARNS
              END IF
              CXB(K)=CXB(K)+CX                  !CXB IS SUM TO GET MU/RHO
              CX=CX/AU                          !CHANGE CX TO ATOMIC UNITS
            END IF
            ICOUNT=6
            RX=XK(K)/C1                 !XRAY ENERGY IN ATOMIC UNITS

            IF (IFBE.EQ.0 .AND. BE.GE.XK(K) ) THEN
              SEDGE=SIGEDG/AU                   !XSECTION IN ATOMIC UNITS AT ENERGY=1.001*BE
              CX=0.0
              FP(K,J)=DGAUSS(1,SIGMA3,5,DLGNDR)*C/(2.0*PI**2)
              MF=3
            ELSE
              IF (IFBE.EQ.0) FP(K,J)=DGAUSS(1,SIGMA0,5,DLGNDR)*C/
     1  (2.0*PI**2)
              IF (IFBE.EQ.1) FP(K,J)=DGAUSS(1,SIGMA1,5,DLGNDR)*C/
     1  (2.0*PI**2)
              IF (IFBE.EQ.2) FP(K,J)=DGAUSS(1,SIGMA2,5,DLGNDR)*C/
     1  (2.0*PI**2)
            END IF
            FPP(K,J)=0.
            IF (CX.NE.0.0) FPP(K,J)=C*CX*RX/(4.0*PI)
            CORR=0.0
            IF (CX.NE.0.0) CORR=-CX*RX*0.5*LOG((RX+BB)/(RX-BB))*C/
     1  (2.0*PI**2)
            IF (MF.EQ.3) CORR=0.5*SEDGE*BB**2*LOG((-BB+RX)/(-BB-RX))/
     1   RX*C/(2.*PI**2)
            FP(K,J)=FP(K,J)+CORR
            ABS(K,J)=CX*AU

          END DO
        END DO

        PRINT 290, NATOM,IZ,NO,WT(IZ),ETERM(IZ),NW,NORD
290     FORMAT (/,' SYMBOL= ',A2,' z,NORB= ',2I3,' AT.WT.,ETERM= ',2F8.3,
     1       ' NW= ',I2,' NORD= ',I1)
        PRINT 260, (XR(I),I=1,NW)
260     FORMAT (/,' LAMBDA ',10F10.3)
        DO I=1,NW
          SUMFP(I)=0.0
          SUMFPP(I)=0.
        END DO
        DO J=1,NO
          DO K=1,NW
            SUMFP(K)=SUMFP(K)+FP(K,J)
          END DO
        END DO
        DO K=1,NW
           XJENSN(K) = 0.0
!          XJENSN(K)=-0.5*FLOAT(IZ)*(XK(K)/C1/137.0367**2)**2
        END DO
        DO K=1,NW
          ET(K)=-ETERM(IZ)
        END DO
        DO K=1,NW
          SUMFP(K)=SUMFP(K)+ET(K)+XJENSN(K)
        END DO
        PRINT 320, (SUMFP(K),K=1,NW)
320     FORMAT (/,'     FP=',10F10.3)
        DO J=1,NO
          DO K=1,NW
            SUMFPP(K)=SUMFPP(K)+FPP(K,J)
          END DO
        END DO
        PRINT 330, (SUMFPP(K),K=1,NW)
330     FORMAT ('    FPP=',10F10.3)
        DO K=1,NW
          CXB(K)=CXB(K)*0.602472/WT(IZ)
        END DO
        PRINT 340, (CXB(K),K=1,NW)
340     FORMAT (' MU/RHO=',10F10.1)
C        READ *
        STOP
        END
C
        SUBROUTINE FPSORT (N,A,B)
        IMPLICIT REAL*8 (A-H,O-Z),INTEGER*4 (I-N)
        DIMENSION A(1),B(1)
        M=N-1
        DO I=1,M
          I1=I+1
          DO J=I1,N
            IF(A(J).LE.A(I)) THEN
              X=A(J)
              Y=A(I)
              A(I)=X
              A(J)=Y
              X=B(J)
              Y=B(I)
              B(I)=X
              B(J)=Y
            END IF
          END DO
        END DO
        RETURN
        END
        SUBROUTINE XSECT (ZX,EL,SL,CX,NX)
        IMPLICIT REAL*8 (A-H,O-Z),INTEGER*4 (I-N)
        DIMENSION EL(14), SL(14)
C*****FIND EL(K) CLOSEST TO ZX
        ER=1000000.
        DO L=1,NX
          P=ABS(ZX-EL(L))
          IF (P.LE.ER) THEN
            ER=P
            LL=L
          END IF
        END DO

        LL=LL-1
        IF (LL.EQ.0) LL=1
        IF (LL.EQ.12) LL=11
        IF (SL(LL).EQ.0.) LL=LL+1
        DET=EL(LL+2)**2*(EL(LL+1)-EL(LL))+EL(LL+1)**2*(EL(LL)-
     1  EL(LL+2))+EL(LL)**2*(EL(LL+2)-EL(LL+1))
        A0=(EL(LL)**2*(SL(LL+1)*EL(LL+2)-SL(LL+2)*EL(LL+1))+
     1  EL(LL+1)**2*(SL(LL+2)*EL(LL)-SL(LL)*EL(LL+2))+EL(LL+2)**2*
     1  (SL(LL)*EL(LL+1)-SL(LL+1)*EL(LL)))/DET
        A1=(EL(LL)**2*(SL(LL+2)-SL(LL+1))+EL(LL+1)**2*(SL(LL)-
     1  SL(LL+2))+EL(LL+2)**2*(SL(LL+1)-SL(LL)))/DET
        A2=(SL(LL)*(EL(LL+2)-EL(LL+1))+SL(LL+1)*(EL(LL)-EL(LL+2))+
     1  SL(LL+2)*(EL(LL+1)-EL(LL)))/DET
        CX=EXP(A0+A1*ZX+A2*ZX**2)
        RETURN
        END
        REAL*8 FUNCTION AKNINT (XBAR,IN,IM,X,Y,T)

        IMPLICIT REAL*8 (A-H,O-Z),INTEGER*4 (I-N)
C
C AITKEN REPEATED INTERPOLATION
C
C   XBAR = ABSCISSA AT WHICH INTERPOLATION IS DESIRED
C   IABS(IN) = NO. OF VALUES IN TABLE
C              IF IN.GT.0, CHK ORDERING OF X(I).
C              IF IN.LT.0, SKIP PRECEEDING TEST.
C   IM   = DEGREE OF APPROXIMATING POLYNOMIAL
C   X    = VECTOR OF IABS(IN) VALUES OF ABSCISSA
C   Y    = VECTOR OF IABS(IN) VALUES OF ORDINATE
C   T    = TEMPORARY STORAGE VECTOR OF 4*(M+1) LOCATIONS)
C
        REAL*8 T,DXBAR
        DIMENSION T(80), X(10), Y(10)

        DXBAR = XBAR
        N = IABS(IN)
        M = IM
        IF ( M.GE.N ) GO TO 120
10      K = N-1
        IF ( N.LT.2 ) GO TO 110
        S = X(2)-X(1)
        IF ( IN.LT.0 ) GO TO 30
C CHK IF ORDER MONOTONIC
        IF ( N.EQ.2 ) GO TO 30
        DO 20 I=3,N
        Z = (X(I)-X(I-1))*S
20      IF ( Z.LE.0. ) GO TO 100
30      IF ( S.LT.0. ) GO TO 50
C INCREASING ORDER
        DO 40 J=1,N
40      IF ( XBAR.LE.X(J) ) GO TO 70
        J = N
        GO TO 70
C DECREASING ORDER
50      DO 60 J=1,N
60      IF ( XBAR.GE.X(J) ) GO TO 70
        J = N
70      K = M
        M = M+1
        J = J-M/2
        J = MAX0(J,1)
        J = MIN0(J,N-K)
        MEND = J+K
        DO I=J,MEND
          KK = I-J+1
          T(KK) = Y(I)
          T(KK+M) = X(I)-DXBAR
        END DO
        DO I = 1,K
          KK = I+1
          DO JJ=KK,M
            T(JJ) = (T(I)*T(JJ+M)-T(JJ)*T(I+M))/(X(JJ+J-1)-X(I+J-1))
          END DO
        END DO
        AKNINT = T(M)
        RETURN
100     PRINT *,' AKNINT X(I) NOT SEQUENCED PROPERLY'
110     PRINT *,' AKNINT N.LT.2 YBAR RETURNED AS Y(1)'
        AKNINT = Y(1)
        RETURN
120     PRINT *,' AKNINT WARNING ORDER OF INTERPOLATION TOO LARGE'
        M = N-1
        GO TO 10
        END
        REAL*8 FUNCTION DGAUSS (N,Y,M,LTBL)

        IMPLICIT REAL*8 (A-H,O-Z),INTEGER*4 (I-N)
        CHARACTER*50 CMT
        DIMENSION M(6),G(6),Z(6),A(6)
        EXTERNAL Y,LTBL
        DATA CMT/'GAUSS   N NOT IN (1,6).  RESULTS   SET TO ZERO.'/

        IF ( (N.GT.6) .OR. (N.LT.1) ) THEN
          CALL LABRT (1,CMT,1)
          DGAUSS = 0.
        ELSE
          NN = N
          GO TO (110,90,70,50,30,10), NN
10        J = 1
          G(6) = 0.
20        CALL LTBL (M(6),J,A(6),Z(6))
30        K = 1
          G(5) = 0.
40        CALL LTBL (M(5),K,A(5),Z(5))
50        L = 1
          G(4) = 0.
60        CALL LTBL (M(4),L,A(4),Z(4))
70        JJ = 1
          G(3) = 0.
80        CALL LTBL (M(3),JJ,A(3),Z(3))
90        KK = 1
          G(2) = 0.
100       CALL LTBL (M(2),KK,A(2),Z(2))
110       LL = 1
          G(1) = 0.
120       CALL LTBL (M(1),LL,A(1),Z(1))
          G(1) = G(1)+A(1)*Y(Z)
          LL = LL+1
          IF ( LL.LE.M(1) ) GO TO 120
          IF ( NN.EQ.1 ) GO TO 130
          G(2) = G(2)+A(2)*Y(Z)*G(1)
          KK = KK+1
          IF ( KK.LE.M(2) ) GO TO 100
          IF ( NN.EQ.2 ) GO TO 130
          G(3) = G(3)+A(3)*Y(Z)*G(2)
          JJ = JJ+1
          IF ( JJ.LE.M(3) ) GO TO 80
          IF ( NN.EQ.3 ) GO TO 130
          G(4) = G(4)+A(4)*Y(Z)*G(3)
          L = L+1
          IF ( L.LE.M(4) ) GO TO 60
          IF ( NN.EQ.4 ) GO TO 130
          G(5) = G(5)+A(5)*Y(Z)*G(4)
          K = K+1
          IF ( K.LE.M(5) ) GO TO 40
          IF ( NN.EQ.5 ) GO TO 130
          G(6) = G(6)+A(6)*Y(Z)*G(5)
          J = J+1
          IF ( J.LE.M(6) ) GO TO 20
130       DGAUSS = G(NN)
        END IF
        RETURN
        END

        REAL*8 FUNCTION SIGMA0 (X)

        IMPLICIT REAL*8 (A-H,O-Z),INTEGER*4 (I-N)
        COMMON /GAUS/ CX, BB, SG(5), RX, ICOUNT

        ICOUNT = ICOUNT-1
        SIGMA0 = SG(ICOUNT)*BB**3/X**2/(RX**2*X**2-BB**2)-
     1       BB*CX*RX**2/(RX**2*X**2-BB**2)
        RETURN
        END

        REAL*8 FUNCTION SIGMA1 (X)

        IMPLICIT REAL*8 (A-H,O-Z),INTEGER*4 (I-N)
        COMMON /GAUS/ CX, BB, SG(5), RX, ICOUNT

        ICOUNT = ICOUNT-1
        SIGMA1 = 0.5*BB**3*SG(ICOUNT)/(SQRT(X)*(RX**2*X**2-BB**2*X))
        RETURN
        END

        REAL*8 FUNCTION SIGMA2 (X)

        IMPLICIT REAL*8 (A-H,O-Z),INTEGER*4 (I-N)
        COMMON /GAUS/ CX, BB, SG(5), RX, ICOUNT

        ICOUNT = ICOUNT-1
        DENOM = X**3*RX**2-BB**2/X
        SIGMA2 = 2.0*BB*SG(ICOUNT)*BB**2/X**4/DENOM-
     1    2.0*BB*CX*RX**2/DENOM
        RETURN
        END

        REAL*8 FUNCTION SIGMA3 (X)

        IMPLICIT REAL*8 (A-H,O-Z),INTEGER*4 (I-N)
        COMMON /EDGE/ SEDGE
        COMMON /GAUS/ CX, BB, SG(5), RX, ICOUNT

        ICOUNT = ICOUNT-1
        SIGMA3 = BB**3*(SG(ICOUNT)-SEDGE*X**2)/(X**2*(X**2*RX**2-BB**2))
        RETURN
        END

        SUBROUTINE LABRT (ISW,LHOL,INX)
        IMPLICIT REAL*8 (A-H,O-Z),INTEGER*4 (I-N)
        CHARACTER*50 LHOL
        LOGICAL PS, TS
        DATA NP /10/, PS /.TRUE./, TS /.FALSE./
        IF ((ISW.EQ.0).OR.(ISW.GT.5)) RETURN
        IF ( ISW.EQ.1 ) THEN
          IF (PS.AND.(NP.GT.0)) PRINT 60, LHOL,INX
60        FORMAT (1H0,9X,A50,3X,I6)
          NP=NP-1
          IF (TS) STOP
        ELSE IF ( ISW.EQ.2) THEN
          PS=.FALSE.
        ELSE IF ( ISW.EQ.3) THEN
          PS=.TRUE.
          NP=INX
        ELSE IF ( ISW.EQ.4) THEN
          TS=.TRUE.
        ELSE IF ( ISW.EQ.5) THEN
          TS=.FALSE.
        END IF
        RETURN
        END

        SUBROUTINE DLGNDR (M,K,AA,Z)

!PURPOSE:

        IMPLICIT REAL*8 (A-H,O-Z),INTEGER*4 (I-N)

!CALLING ARGUMENTS:

        INTEGER*4       M       !
        INTEGER*4       K       !
        REAL*8          AA      !
        REAL*8          Z       !

!INCLUDE STATEMENTS:

!LOCAL VARIABLES:

        REAL*8          A(68), X(62)

!DATA STATEMENTS:

        DATA X(1)       /.06943184420297/
        DATA X(2)       /.33000947820757/
        DATA X(3)       /.04691007703067/
        DATA X(4)       /.23076534494716/
        DATA X(5)       /.03376524289992/
        DATA X(6)       /.16939530676687/
        DATA X(7)       /.38069040695840/
        DATA X(8)       /.02544604382862/
        DATA X(9)       /.12923440720030/
        DATA X(10)      /.29707742431130/
        DATA X(11)      /.01985507175123/
        DATA X(12)      /.10166676129319/
        DATA X(13)      /.23723379504184/
        DATA X(14)      /.40828267875217/
        DATA X(15)      /.01591988024619/
        DATA X(16)      /.08198444633668/
        DATA X(17)      /.19331428364971/
        DATA X(18)      /.33787328829809/
        DATA X(19)      /.01304673574141/
        DATA X(20)      /.06746831665551/
        DATA X(21)      /.16029521585049/
        DATA X(22)      /.28330230293537/
        DATA X(23)      /.42556283050918/
        DATA X(24)      /.01088567092697/
        DATA X(25)      /.05646870011595/
        DATA X(26)      /.13492399721298/
        DATA X(27)      /.24045193539659/
        DATA X(28)      /.36522842202382/
        DATA X(29)      /.00921968287664/
        DATA X(30)      /.04794137181476/
        DATA X(31)      /.11504866290285/
        DATA X(32)      /.20634102285669/
        DATA X(33)      /.31608425050091/
        DATA X(34)      /.43738329574426/
        DATA X(35)      /.00790847264071/
        DATA X(36)      /.04120080038851/
        DATA X(37)      /.09921095463335/
        DATA X(38)      /.17882533027983/
        DATA X(39)      /.27575362448178/
        DATA X(40)      /.38477084202243/
        DATA X(41)      /.00685809565159/
        DATA X(42)      /.03578255816821/
        DATA X(43)      /.08639934246512/
        DATA X(44)      /.15635354759416/
        DATA X(45)      /.24237568182092/
        DATA X(46)      /.34044381553605/
        DATA X(47)      /.44597252564632/
        DATA X(48)      /.600374098758E-2/
        DATA X(49)      /.31363303799647E-1/
        DATA X(50)      /.75896708294787E-1/
        DATA X(51)      /.13779113431991/
        DATA X(52)      /.21451391369574/
        DATA X(53)      /.30292432646121/
        DATA X(54)      /.39940295300128/
        DATA X(55)      /.00529953250417/
        DATA X(56)      /.02771248846338/
        DATA X(57)      /.06718439880608/
        DATA X(58)      /.12229779582250/
        DATA X(59)      /.19106187779868/
        DATA X(60)      /.27099161117138/
        DATA X(61)      /.35919822461038/
        DATA X(62)      /.45249374508118/
        DATA A(1)       /.17392742256873/
        DATA A(2)       /.32607257743127/
        DATA A(3)       /.11846344252810/
        DATA A(4)       /.23931433524968/
        DATA A(5)       /.28444444444444/
        DATA A(6)       /.85662246189585E-1/
        DATA A(7)       /.18038078652407/
        DATA A(8)       /.23395696728635/
        DATA A(9)       /.06474248308443/
        DATA A(10)      /.13985269574464/
        DATA A(11)      /.19091502525256/
        DATA A(12)      /.20897959183674/
        DATA A(13)      /.05061426814519/
        DATA A(14)      /.11119051722669/
        DATA A(15)      /.15685332293894/
        DATA A(16)      /.18134189168918/
        DATA A(17)      /.04063719418079/
        DATA A(18)      /.09032408034743/
        DATA A(19)      /.13030534820147/
        DATA A(20)      /.15617353852000/
        DATA A(21)      /.16511967750063/
        DATA A(22)      /.03333567215434/
        DATA A(23)      /.07472567457529/
        DATA A(24)      /.10954318125799/
        DATA A(25)      /.13463335965500/
        DATA A(26)      /.14776211235738/
        DATA A(27)      /.02783428355809/
        DATA A(28)      /.06279018473245/
        DATA A(29)      /.09314510546387/
        DATA A(30)      /.11659688229599/
        DATA A(31)      /.13140227225512/
        DATA A(32)      /.13646254338895/
        DATA A(33)      /.02358766819326/
        DATA A(34)      /.05346966299766/
        DATA A(35)      /.08003916427167/
        DATA A(36)      /.10158371336153/
        DATA A(37)      /.11674626826918/
        DATA A(38)      /.12457352290670/
        DATA A(39)      /.02024200238266/
        DATA A(40)      /.04606074991886/
        DATA A(41)      /.06943675510989/
        DATA A(42)      /.08907299038097/
        DATA A(43)      /.10390802376845/
        DATA A(44)      /.11314159013145/
        DATA A(45)      /.11627577661544/
        DATA A(46)      /.01755973016588/
        DATA A(47)      /.04007904357988/
        DATA A(48)      /.06075928534395/
        DATA A(49)      /.07860158357910/
        DATA A(50)      /.09276919873897/
        DATA A(51)      /.10259923186065/
        DATA A(52)      /.10763192673158/
        DATA A(53)      /.01537662099806/
        DATA A(54)      /.03518302374405/
        DATA A(55)      /.05357961023359/
        DATA A(56)      /.06978533896308/
        DATA A(57)      /.08313460290850/
        DATA A(58)      /.09308050000778/
        DATA A(59)      /.09921574266356/
        DATA A(60)      /.10128912096278/
        DATA A(61)      /.01357622970588/
        DATA A(62)      /.03112676196932/
        DATA A(63)      /.04757925584125/
        DATA A(64)      /.06231448562777/
        DATA A(65)      /.07479799440829/
        DATA A(66)      /.08457825969750/
        DATA A(67)      /.09130170752246/
        DATA A(68)      /.09472530522754/

!CODE:

        KK = K
        IF ((M.GT.16).OR.(M.LT.4)) KK=4
        IS = 0
        IH = (M+1)/2
        Z = .5
        IF (MOD(M,2).EQ.1) IS=-1
        IP = KK
        T = 0.
        IF (IP.GT.IH) THEN
          IP = M+1-IP
          T = -1
        END IF
        I4 = M-4
        IA = (I4*(M+4)+IS)/4+IP
        AA = A(IA)
        IF ((IP.EQ.IH).AND.(IS.LT.0)) RETURN
        IA = IA-(I4+IS)/2
        Z = -T+SIGN(X(IA),T)
        RETURN
        END
