      SUBROUTINE CCNFRC ( NUMPTS, XS, YS )                              00000010
CDC      SUBROUTINE CCNFRC ( NUMPTS, XS, YS, X, Y )                     00000020
C                                                                       00000030
C     THIS ROUTINE WILL INTERPOLATE OR EXTRAPOLATE A COMPLEX FUNCTION   00000040
C     THAT IS TABULATED ON A NONUNIFORM COMPLEX GRID.  A SPECIAL FORM OF00000050
C     CONTINUED FRACTION IS USED.  NO TWO INPUT VALUES OF THE           00000060
C     FUNCTION MAY BE EQUAL.  THIS ALGORITHM IS QUITE SUCCESSFUL AT     00000070
C     FUNCTIONS AROUND POLES AND MODERATELY SUCCESSFUL AT CONTINUING    00000080
C     AROUND CUTS AND OTHER ESSENTIAL SINGULARITIES.                    00000090
C                                                                       00000100
C     ENTRY  CCNFRC  IS USED TO SETUP THE COEFFICIENTS OF THE COMPLEX   00000110
C     CONTINUED FRACTION.  THEN THE ENTRY  CCONTF  MAY BE USED AS       00000120
C     MANY TIMES AS IS DESIRED TO GET VALUES OF THE FUNCTION AT         00000130
C     NEW VALUES OF THE ORDINATE.  IN ADDITION THE CONVERGENCE          00000140
C     TO THE INTERPOLATED/EXTRAPOLATED VALUE IS AVAILABLE FOR           00000150
C     INSPECTION.                                                       00000160
C                                                                       00000170
C                                                                       00000180
C     NUMPTS IS THE NUMBER OF INPUT VALUES OF XS AND YS.                00000190
C     XS IS THE ARRAY OF INPUT X'S AT WHICH THE FUNCTION VALUES ARE     00000200
C        TABULATED.  THIS MUST BE A COMPLEX*16 ARRAY.                   00000210
C        IT IS REARRANGED DURING THE CCNFRC CALCULATION, AND MUST BE    00000220
C        LEFT THAT WAY THROUGH ALL CCONTF CALLS.                        00000230
C     YS, ON INPUT, IS THE ARRAY OF INPUT FUNCTION VALUES:              00000240
C            YS(I) = Y(XS(I))                                           00000250
C        IT MUST BE A COMPLEX*16 ARRAY.                                 00000260
C     ON OUTPUT, YS GETS THE COEFFICENTS.                               00000270
C                                                                       00000280
C                                                                       00000290
C                                                                       00000300
C     VERSION OF 4/18/70 = 70.108                                       00000310
C     SAVED ON DISK 7/17/70                                             00000320
C     9/26/72  ARGONNE VERSION FOR CONTIN                               00000330
C     10/02/72   COMPLEX VERSION                                        00000340
C     5/1/74 - MADE BACK TO REAL VERSION                                00000350
C     5/23/74 - AND BACK TO COMPLEX AGAIN.                              00000360
C     8/11/75 - PTOLEMY VERSION BASED ON SPEAKEASY LINKULE VERSION.     00000370
C     7/23/79 - CHANGE ORDER, USE BIGGEST FOR BETTER ACCURACY, PASS     00000380
C        CCNFRC ARGUMENTS TO CCONTF - RPG                               00000390
C     1/31/80 - END STATEMENT, COMPLEX*8 VERSION FOR VAX - RPG          00000400
C                                                                       00000410
      IMPLICIT REAL*8 ( A-H, O-Z )                                      00000420
	COMPLEX*16 XS,YS,X,Y
CDC      ALLOCLEVEL XS,YS                                               00000440
C                                                                       00000450
	dimension ys(1),xs(1)
C                                                                       00000470
	COMPLEX*16 D,XJ,YJ
CFR%F                                                                   00000490
C                                                                       00000500
C                                                                       00000510
      DIMENSION DR(2)                                                   00000530
      EQUIVALENCE ( DR(1), D )                                          00000540
C     THIS IS THE SQUARE OF THE RELATIVE ACCURACY NEEDED IN THE         00000550
C     INTERPOLATED FUNCTION.                                            00000560
      DATA TINY / 1.D-14 /                                              00000570
C                                                                       00000580
C                                                                       00000590
C     REDUCE THE INPUT ARRAYS TO SINGLE PRECISION FOR THE VAX.          00000600
      NMAX = NUMPTS - 1                                                 00000630
C                                                                       00000640
C     SEARCH FOR THE LARGEST Y.                                         00000650
C                                                                       00000660
      COMP = -1                                                         00000670
      DO 109 I = 1, NUMPTS                                              00000680
      D = YS(I)                                                         00000690
      TEMP = DR(1)**2 + DR(2)**2                                        00000700
      IF ( TEMP .LE. COMP ) GO TO 109                                   00000710
      COMP = TEMP                                                       00000720
      K = I                                                             00000730
 109  CONTINUE                                                          00000740
      IF ( NMAX .LE. 0 ) RETURN                                         00000750
C                                                                       00000760
C     THE FIRST COEFFICIENT IS THE BIGGEST Y.                           00000770
C     EXCHANGE THE FIRST WITH THE BIGGEST.                              00000780
C                                                                       00000790
      YJ = YS(K)                                                        00000800
      IF ( K .EQ. 1 ) GO TO 120                                         00000810
      YS(K) = YS(1)                                                     00000820
      YS(1) = YJ                                                        00000830
      XJ = XS(K)                                                        00000840
      XS(K) = XS(1)                                                     00000850
      XS(1) = XJ                                                        00000860
C                                                                       00000870
C     CALCULATE THE CANDIDATES FOR THE SECOND COEFFICIENT,              00000880
C     AND SEARCH FOR THE LARGEST.                                       00000890
C                                                                       00000900
 120  COMP = -1                                                         00000910
      DO 139 I = 2, NUMPTS                                              00000920
      D = 1. - YJ/YS(I)                                                 00000930
      YS(I) = D                                                         00000940
      TEMP = DR(1)**2 + DR(2)**2                                        00000950
      IF ( TEMP .LE. COMP ) GO TO 139                                   00000960
      COMP = TEMP                                                       00000970
      K = I                                                             00000980
 139  CONTINUE                                                          00000990
C                                                                       00001000
C     LOOP THROUGH THE REST OF THE COEFFICIENTS.                        00001010
C                                                                       00001020
      DO 400 J = 2, NUMPTS                                              00001030
C                                                                       00001040
C        EXCHANGE THE LARGEST CANDIDATE WITH THE J-TH, AND              00001050
C        FINISH CALCULATING THIS COEFFICIENT.                           00001060
C                                                                       00001070
      XJ = XS(K)                                                        00001080
      YJ = YS(K)                                                        00001090
      IF ( K .EQ. J ) GO TO 209                                         00001100
      YS(K) = YS(J)                                                     00001110
      XS(K) = XS(J)                                                     00001120
      XS(J) = XJ                                                        00001130
 209  YJ = YJ/(XS(J-1)-XJ)                                              00001140
      YS(J) = YJ                                                        00001150
      IF ( J .EQ. NUMPTS ) GO TO 450                                    00001160
C                                                                       00001170
C        QUIT IF THIS COEFFICIENT IS TINY.                              00001180
C        NOTE:  THIS COMPARISON IS MEANINGFUL ONLY IF THE SIZE OF THE   00001190
C        DIFFERENCE XS(J-1)-XS(J) IS MORE OR LESS TYPICAL OF THE RANGE  00001200
C        OF X.  IF THIS DIFFERENCE IS MANY ORDERS OF MAGNITUDE SMALLER  00001210
C        THAN TYPICAL X DIFFERENCES, IT MAY QUIT PREMATURELY, BUT THEN  00001220
C        THE WHOLE PROBLEM IS PATHOLOGICAL ANYWAY.                      00001230
C                                                                       00001240
      IF ( COMP .LT. TINY ) GO TO 430                                   00001250
C                                                                       00001260
C        CALCULATE THE CANDIDATES FOR THE NEXT COEFFICIENT,             00001270
C        AND SEARCH FOR THE LARGEST.                                    00001280
C                                                                       00001290
      JP1 = J + 1                                                       00001300
      COMP = -1                                                         00001310
      XJ = XS(J-1)                                                      00001320
      DO 309 I = JP1, NUMPTS                                            00001330
      D = 1. + YJ*(XS(I)-XJ)/YS(I)                                      00001340
      YS(I) = D                                                         00001350
      TEMP = DR(1)**2 + DR(2)**2                                        00001360
      IF ( TEMP .LE. COMP ) GO TO 309                                   00001370
      COMP = TEMP                                                       00001380
      K = I                                                             00001390
 309  CONTINUE                                                          00001400
C                                                                       00001410
 400  CONTINUE                                                          00001420
C                                                                       00001430
C     END OF LOOP THROUGH COEFFICIENTS.                                 00001440
C     (IT NEVER GETS HERE ANYWAY.)                                      00001450
C                                                                       00001460
C     THE INPUT FUNCTION APPEARS TO BE A CONTINUED FRACTION EXACTLY     00001470
 430  NMAX = I-1                                                        00001480
      PRINT 433, NMAX, NUMPTS                                           00001490
 433  FORMAT ( '0**** WARNING:  CONTINUED FRACTION USED ONLY', I3,  ' OU00001500
     *T OF', I3, ' POINTS.' )                                           00001500
C                                                                       00001520
 450  CONTINUE                                                          00001530
      RETURN                                                            00001540
C                                                                       00001550
C    HERE WE COMPUTE ONE VALUE OF THE CONTINUED FRACTION                00001560
C                                                                       00001570
C      X IS THE  LOCATION OF THE DESIRED FUNCTION VALUE                 00001580
C      Y WILL BE THE FUNCTION                                           00001590
C     BOTH X AND Y MUST BE COMPLEX*16.                                  00001600
C                                                                       00001610
C                                                                       00001620
C                                                                       00001630
      ENTRY CCONTF ( NUMPTS, XS, YS, X, Y )                             00001640
C                                                                       00001650
 500  Y = 0                                                             00001660
      IF ( NMAX .LT. 1 ) GO TO 520                                      00001670
      DO 509 J=1, NMAX                                                  00001680
      K = NMAX - J + 1                                                  00001690
      Y = YS(1+K)*(X-XS(K))/(1.+Y)                                      00001700
 509  CONTINUE                                                          00001710
 520  Y = YS(1)/(1.+Y)                                                  00001720
      RETURN                                                            00001730
      END                                                               00001740
C                                                                       00001750
C     THIS SUBROUTINE COMPRESSES A REAL*8 ARRAY INTO A REAL*4 ARRAY.    00001760
C     THE ARGUMENTS MAY BE (AND ARE) THE SAME ARRAY.                    00001770
C                                                                       00001780
CIB      FUNCTION CLEBSH (/A/, /B/, /X/, /Y/, /CIN/, /ZIN/)             00000100
      FUNCTION CLEBSH ( A, B, ARG3, ARG4, ARG5, ZIN )
      IMPLICIT INTEGER (A-C,I-Z), REAL*8 (D-H)
      REAL*8 CLEBSH 
      REAL*8 THREEJ
C                                                                       00000300
C     COMPUTES CLEBSCH-GORDAN AND 3-J COEFICENTS.                       00000400
C                                                                       00000500
C       THIS ROUTINE COMPUTES THE CLEBSH GORDAN OR 3-J COEFICENT FOR    00000600
C     A GIVEN SET OF J'S AND M'S.  THE INPUT CONSISTS OF 6 INTEGERS THAT00000700
C     ARE TWICE THE VALUES OF THE J'S AND M'S.  IN THIS MANNER          00000800
C     HALF-INTEGER ANGULAR MOMENTA CAN BE ACCOMODATED.                  00000900
C                                                                       00001000
C       A FACTORIZATION OF THE RACAH SUM SUGGESTED BY                   00001100
C         J. G. WILLS, COMPUTER PHYSICS COMM. 2, 381 (1971)             00001200
C     IS USED.  THE REQUIRED FACTORIALS THAT MULTIPLY THE SUM ARE       00001300
C     FOUND IN ONE OF 3 WAYS:                                           00001400
C       FOR J1+J2+J3 < 96,        A  TABLE OF SQRT(FACTORIAL) IS USED;  00001500
C       FOR 95 < J1+J2+J3 < 1000, A TABLE OF LOG(FACTORIAL) IS USED;    00001600
C       FOR 1000 <= J1+J2+J3,     DLGAMA IS USED.                       00001700
C                                                                       00001800
C       THESE CHOICES REPRESENT REASONABLE TRADE OFFS BETWEEN SIZE AND  00001900
C     SPEED CONSIDERATIONS.  BY RECOMPILING THE BLOCK DATA PROGRAMS     00002000
C     FOR THE /FACTRL/ OR /LOGFAC/ COMMON BLOCKS, ONE CAN CHANGE        00002100
C     THE ABOVE BOUNDRIES TO EITHER REDUCE CORE (IF ONE IS NOT          00002200
C     INTERESTED IN LARGE J'S) OR INCREASE THE SPEED FOR LARGE J'S.     00002300
C                                                                       00002400
C       THE MAGNITUDE OF THE J'S AND M'S THAT CAN BE ACCOMODATED BY     00002500
C     THIS ROUTINE IS DETERMINED BY CANCELLATIONS IN THE                00002600
C     RACAH SUM.  AS CAN BE SEEN IN THE TABLE BELOW, IF ALL THREE J'S   00002700
C     GET LARGE, SEVERE CANCELLATION SETS IN FOR EACH J OF THE ORDER OF 00002800
C     75, WHILE OVERFLOWS BEGIN TO OCCUR FOR THE J'S OF THE ORDER OF    00002900
C     150.  HOWEVER IF ONE OF THE J'S IS HELD SMALL, THE PRECISION      00003000
C     REMAINS GOOD FOR VERY LARGE VALUES OF THE OTHER TWO J'S AND THERE 00003100
C     ARE NO OVERFLOW PROBLEMS (THE LARGEST VALUES TRIED WERE 2000 AS   00003200
C     IS SHOWN IN THE TABLE).                                           00003300
C                                                                       00003400
C       A SPECIAL CASE IS USED FOR M1 = M2 = M3 = 0.  ALL THREE J'S     00003500
C     MAY THEN BE VERY LARGE.                                           00003600
C                                                                       00003700
C       EACH LINE IN THE FOLLOWING TABLE IS THE RESULT OF SEVERAL       00003800
C     THOUSAND EVALUATIONS OF CLEBSH GORDAN COEFFICENTS FOR J'S AND     00003900
C     M'S IN THE RANGES INDICATED.  PRECISIONS WERE DETERMINED BY       00004000
C     COMPARING THE RESULTS TO A QUADRUPOLE PRECISION CLEBSH GORDAN     00004100
C     ROUTINE.  IT IS NOT SIGNIFICANT THAT J3 WAS HELD SMALL; ANY       00004200
C     OF THE J'S MAY BE KEPT SMALL WITH THE SAME TIMING AND RELATIVE    00004300
C     ERROR RESULTS.  THE EXECUTION TIMES ARE GIVEN IN MICROSECONDS.    00004400
C                                                                       00004500
C       J1, J2    J3   M1, M2     TIMES    RMS ABSOLUTE                 00004600
C        RANGE   RANGE  LIMIT   /195  /75     ERROR                     00004700
C        0   5   0   5    5      40   340     1E-16                     00004800
C        5  10   0  10   10      48   425     8E-17                     00004900
C       10  15   0  25   13      51   465     8E-17                     00005000
C       15  30   0  47   30      59   560     7E-17                     00005100
C       33  42   0  75   42      84   715     9E-13                     00005200
C       50  52   0 100   52      95   835     2E-10                     00005300
C       74  75   0 150   75     103   935     2E-5                      00005400
C        100     0 200  100     124  1175     .6                        00005500
C        150     0 300  150     205  2000     4E+8   SOME OVERFLOWS     00005600
C        100     0  10   10      63   455     2E-15                     00005700
C        400     0  10   10      63   455     1E-14                     00005800
C        500     0  10   10     465  3500     6E-14                     00005900
C        500     0  10  500     470           6E-14  NOTE LARGE M1, M2  00006000
C       2000     0  10   10     460           2E-13                     00006100
C        0  10   0  10    0      37   185     8E-17  M1=M2=M3=0         00006200
C        100      100     0      37   210     5E-15                     00006300
C        300      300     0      38   215     3E-12                     00006400
C        500      500     0     172  1320                               00006500
C       1000     1000     0     174           3E-12                     00006600
C                                                                       00006700
C     NO EVALUATION OF CLEBSCH-GORDAN COEFICENTS THAT WERE TRIVIALLY    00006800
C     ZERO WERE INCLUDED IN THE ABOVE TIMINGS.  THE LONG EXECUTION      00006900
C     TIMES FOR THE  J1, J2 > 500  ENTRIES ARE DUE TO THE USE OF        00007000
C     THE DLGAMA FUNCTION.  THEY MAY BE REDUCED BY INCREASING THE       00007100
C     SIZE OF THE  /LOGFAC/  TABLE.                                     00007200
C                                                                       00007300
C                                                                       00007400
C     THERE ARE THREE ENTRIES:                                          00007500
C                                                                       00007600
C     CLEBSH (2*J1, 2*J2, 2*M1, 2*M2, 2*J3, 2*M3)   RETURNS A           00007700
C                  DOUBLE PRECISION CLEBSCH-GORDAN COEFICENT.           00007800
C     THREEJ (2*J1, 2*J2, 2*J3, 2*M1, 2*M2, 2*M3)   RETURNS A           00007900
C                  DOUBLE PRECISION 3-J SYMBOL (M3 = -(M1+M2)).         00008000
C     COFCG (2*J1, 2*J2, 2*J3, 2*M1, 2*M2, 2*M3)    RETURNS A           00008100
C                  SINGLE PRECISION CLEBSCH-GORDAN COEFICENT.           00008200
C                                                                       00008300
C     IN ALL CASES THE TRIANGLE INEQUALITIES, THE REQUIREMENTS THAT     00008400
C     EACH (J, M) SATISFY   !M! <= J,   AND THE REQUIREMENT THAT        00008500
C     M1+M2 = M3 (FOR 3-J M1+M2+M3 = 0)  ARE CHECKED.  IF THEY ARE      00008600
C     VIOLATED, 0 IS RETURNED.  HOWEVER NO CHECK IS MADE TO VERIFY      00008700
C     THAT EACH (J, M) PAIR IS EITHER 1/2 INTEGER OR INTEGER NOR        00008800
C     THAT J1+J2+J3 IS INTEGER.  THE USER MUST INSURE THAT              00008900
C        2*J+2*M IS EVEN   AND THAT  2*J1+2*J2+2*J3  IS EVEN.           00009000
C     IF THIS IS NOT SO, MEANINGLESS RESULTS WILL BE RETURNED.          00009100
C                                                                       00009200
C       CORE REQUIREMENTS (BASE 10):                                    00009300
C                                                                       00009400
C          CLEBSH    4088 BYTES                                         00009500
C         /FACTRL/ 798     MAY BE REDUCED                               00009600
C         /LOGFAC/   8028     MAY BE REDUCED OR INCREASED               00009700
C          DLGAMA    1008                                               00009800
C          TOTAL   13,922                                               00009900
C                                                                       00010000
C     ALSO "DSQRT", "DLOG", "DEXP" AND THE FORTRAN I/O PACKAGE ARE USED.00010100
C     SINCE THESE ARE LIKELY TO BE USED BY OTHER PARTS OF YOUR          00010200
C     PROGRAM, THEY ARE NOT INCLUDED IN THE ABOVE (THE FORTRAN I/O      00010300
C     PACKAGE IS USED BY DLGAMA, DLOG, DEXP, AND DSQRT).  THE  /FACTRL/ 00010400
C     AND  /LOGFAC/  COMMON BLOCKS AND THE "DLGAMA" FUNCTION WILL BE    00010500
C     SHARED WITH "SIXJ" AND "RACAH" SHOULD YOU ALSO USE THEM.          00010600
C                                                                       00010700
C                                                                       00010800
C     S. PIEPER                                                         00010900
C                                                                       00011000
C     4/5/71                                                            00011100
C     8/21/72  ARGONNE VERSION IN DOUBLE PRECISION                      00011200
C     8/25/72  INTEGER INPUT VERSION                                    00011300
C     9/3/72  SPEEDED UP VERSION, ALSO QUITE ACCURATE                   00011400
C     9/7/72  USE SQRT OF FACTORIALS AND BINOMIAL COFS                  00011500
C     9/8/72  SINGLE AND DOUBLE PRECISION ENTRY POINTS                  00011600
C     10/22/72  ALMOST COMPLETE REWRITE FOR RACAH FORMULA               00011700
C     10/23/72  ADD SPECIAL CASE FOR M1 = M2 = M3 = 0 AND 3-J ENTRY     00011800
C     5/17/73  CORRECT SIGN ERROR IN THREEJ                             00011900
C     6/18/74 - TOTALLY NEW METHOD                                      00012000
C     7/8/74 - USE "XL" FUNCTIONS FOR SPEED, REDUCED SIZE.              00012100
C     9/11/74 - CORRECT THREEJ DOUBLE/SINGLE PRECISION ERROR            00012200
C     12/28/79 - CDC TO CNI, STANDARDIZE - RPG                          00012300
C     1/16/80 - CNI ERRORS: STMT 120, AINT TO DINT - RPG                00012400
C     3/13/80 - REMOVE NON-TBIT STUFF -- LET MACROS DO IT - SCP.        00012500
C     5/28/80 - CDC TEST OF TRIANGLE; FIX 3-J DECLARATION - S.P.        00012600
C     7/7/80 - REAL ALGAMA FOR CDC, CRA - S.P.                          00012700
C     7/9/80 - CRAY TRIANGLE TEST - S.P.                                00012800
C     7/31/80 - USE LOGFAC TABLE IN ALLOCATOR - S.P.                    00012900
C     8/25/80 REMOVE COFCG FOR PTOLEMY - S.P.                           00012950
C                                                                       00013100
C                                                                       00013200
C                                                                       00013300
      REAL*8 ALLOC
	COMMON/ALLOCS/ALLOC(1)
      REAL*8 SUM, ANSWER                                                00013600
CNV      REAL ALGAMA                                                    00013700
CIB      REAL*8  CLEBSH                                                 00013800
	CHARACTER*8  LOGNAM / 'LOGFACS ' /
CK     INTEGER*4 Z
	COMMON /LOCKT/ Z(100)
C                                                                       00014000
      COMMON /FACTRL/ MAXFAC, ISPACE,FACTBL(57)                         00014100
C                                                                       00014200
      REAL*8 LF                                                         00014300
      COMMON /LOGFAC/ MAXLF, LFBIAS, LF(1)                              00014400
CRA      COMMON /LOGFAC/  MAXLF, LFBIAS, LFTBL(1)                       00014500
C                                                                       00014600
CRA       POINTER ( LFADDR, LF(1) )                                     00014800
C                                                                       00014900
CIBC     FOLLOWING ARE THE FORTRAN H XL FUNCTIONS THAT WE USE.          00015000
CIBC     TBIT CHECKS TO SEE IF A BIT IS ON.  TBIT(I,31) CHECKS FOR      00015100
CIBC          ODD NUMBER, TBIT(I, 30) CHECKS FOR I/2 BEING ODD.         00015200
CIBC     SHFTR SHIFTS RIGHT THE INDICATED NUMBER OF BITS  = /2**N.      00015300
C                                                                       00015400
CIB      LOGICAL TBIT                                                   00015500
C                                                                       00015600
CIB      REAL*8 SHIFT /Z4F000000 00000000/,  ADDER /.5D0/               00015700
C                                                                       00015800
      LOGICAL DBLSW, THRESW                                             00015900
CXXCIB      REAL*4 COFCG                                                00016000
CXXC                                                                    00016100
CXXCIBC     THE FOLLOWING ARE FOR ROUNDING TO SINGLE PRECISION OUTPUT.  00016200
CXXCIB      REAL*8  ROUND / Z 00 000000 80000000 /                      00016300
CXXCIB      LOGICAL*1 ANSEXP, RNDEXP                                    00016400
CXXCIB      EQUIVALENCE (ANSEXP, ANSWER), (RNDEXP, ROUND)               00016500
C                                                                       00016600
C     FOLLOWING ARRAYS ARE DEFINED TO ALLOW USE OF LOOPS TO SHORTEN CODE00016700
C                                                                       00016800
      DIMENSION DRAY(9)                                                 00016900
C                                                                       00017100
      EQUIVALENCE (DRAY(1), DA1), (DRAY(2), DA2), (DRAY(3), DA3),  (DRAY00017200
     *(4), DY1), (DRAY(5), DX2), (DRAY(6), DY2), (DRAY(7), DX1),  (DRAY(00017200
     *8), DZ1), (DRAY(9), DZ2)                                          00017200
C                                                                       00017500
C     FOLLOWING DUMMY EXTERNAL CAUSES THE FACTORIAL TABLES TO BE DEFINED00017600
C                                                                       00017700
CIB      EXTERNAL  FACDUM, LOGDUM                                       00017800
C                                                                       00017900
C                                                                       00018000
C      CLEBSH  IS THE AMD COMPATABLE ENTRY POINT                        00018100
C                                                                       00018200
C                                                                       00018300
C                                                                       00018400
      THRESW = .FALSE.                                                  00018500
CXXCIB      DBLSW = .TRUE.                                              00018600
CXX      GO TO 10                                                       00018700
CXXCIBC                                                                 00018800
CXXCIBC     COFCG IS SINGLE PRECISION ENTRY                             00018900
CXXCIBC                                                                 00019000
CXXCIB      ENTRY COFCG (/A/, /B/, /CIN/, /X/, /Y/, /ZIN/)              00019100
CXXCIB      THRESW = .FALSE.                                            00019200
CXXCIB      DBLSW = .FALSE.                                             00019300
C                                                                       00019400
 10   ZZ= ZIN                                                           00019500
C                                                                       00019600
C                                                                       00019700
      X = ARG3                                                          00019800
      Y = ARG4                                                          00019900
      C = ARG5                                                          00020000
C                                                                       00020100
  20  CONTINUE                                                          00020200
C                                                                       00020300
CIB  20  C = CIN                                                        00020400
C                                                                       00020500
CRA      LFADDR = LOCF(LFTBL) + LFBIAS                                  00020600
C                                                                       00020700
      IF (ZZ .EQ. X+Y) GO TO 120                                         00020800
C                                                                       00020900
  110 CLEBSH=0.0                                                        00021000
      RETURN                                                            00021100
C                                                                       00021200
C                                                                       00021300
C                                                                       00021400
C     X1 = J1 + M1 ;  ETC.                                              00021500
C                                                                       00021600
C     THE FOLLOWING SHFTR FUNCTIONS EFFECT A DIVIDE BY 2                00021700
C     TO GET FROM THE 2*J STUFF TO J STUFF SINCE THE 3                  00021800
C     (J+M)'S MUST BE INTEGERS.                                         00021900
120      IF(A+X.LT.0)GO TO 110
	IF(B+Y.LT.0)GO TO 110
	IF(C+ZZ.LT.0)GO TO 110
	X1=JISHFT(A+X,-1)
      Y1 = JISHFT(B+Y, -1)                                                000222
      Z1 = JISHFT(C+ZZ, -1)                                                000223
C                                                                       00022400
C     X2 = J1 - M1 ;  ETC                                               00022500
C                                                                       00022600
      X2 = X1 - X                                                       00022700
      Y2 = Y1 - Y                                                       00022800
      Z2 = Z1 - ZZ                                                       00022900
	ILOGFC=NAMLOC(LOGNAM)
	LLOGFC=Z(ILOGFC)-1
C                                                                       00023000
C     THE  SHFTR  INSTRUCTION IS A LOGICAL SHIFT.  THUS IT CONVERTS     00023100
C     A NEGATIVE NUMBER INTO A LARGE POSITIVE NUMBER WHICH WILL LOUSE   00023200
C     UP THE TRIANGLE CONDITION CHECKS.  HERE WE OR ALL THREE NUMBERS   00023300
C     TOGETHER AND THEN CHECK TO SEE IF ANY OF THEM HAD SIGN BITS       00023400
C     THAT WERE SHIFTED INTO THE FIRST POSITION.  IF SO THEN THE        00023500
C     CORRESPONDING M WAS GREATER THAN J SO WE JUST RETURN 0.           00023600
C     ALL OF THIS IS TO AVOID DIVIDES BY 2 WHICH ARE VERY SLOW ON       00023700
C     THE /195.                                                         00023800
CCDC      ON THE CDC THE SHIFT PRESERVES THE SIGN BUT THE OR            00023900
CCDC      OF A NEGATIVE NUMBER WITH ANYTHING REMAINS NEGATIVE           00024000
CCDC      SO WE TEST FOR THAT.                                          00024100
C     ON CRAY WE TEST FOR A HUGH INTEGER                                00024200
C                                                                       00024300
C                                                                       00024500
C     R = J1+J2+J3 + 1                                                  00024600
C                                                                       00024700
      R = X1 + Y1 + Z2 + 1                                              00024800
C                                                                       00024900
C     CARRY OUT THE ABOVE MENTIONED CHECK                               00025000
C                                                                       00025100
CIB      IF (TBIT(TEST, 1))  GO TO 110                                  00025200
CDC      IF ( AND(TEST, -1)  .LT. 0 )  GO TO 110                        00025300
CRA      IF ( TEST .GT. 1 000 000 000 )  GO TO 110                      00025400
C                                                                       00025500
C     A(I) = J1+J2+J3 - 2J(I) ;  A1 = -J1 + J2 + J3 , ETC.              00025600
C                                                                       00025700
      A1 = R-1 - A                                                      00025800
      A2 = R-1 - B                                                      00025900
      A3 = R-1 - C                                                      00026000
C                                                                       00026100
C     B1 = J3 - J1 - M2 ;   B2 = J3 - J2 + M1                           00026200
C                                                                       00026300
      B1 = A1 - Y1                                                      00026400
	B2=A2-X2
C                                                                       00026600
C     WE WILL NEED THESE FLOATING POINT NUMBERS AT VARIOUS PLACES       00026700
C                                                                       00026800
      DA1 = A1                                                          00026900
      DA2 = A2                                                          00027000
      DA3 = A3                                                          00027100
      DY1 = Y1                                                          00027200
      DX2 = X2                                                          00027300
C                                                                       00027400
      DB1 = DA1 - DY1                                                   00027500
      DB2 = DA2 - DX2                                                   00027600
C                                                                       00027700
      SUM = 1                                                           00027800
C                                                                       00027900
      NMIN = MAX0 (0, -B1, -B2)                                         00028000
      NMAX = MIN0 (A3, X2, Y1)                                          00028100
C                                                                       00028200
C     THIS ONE CHECK MAKES SURE OF THE FOLLOWING                        00028300
C     1)  ALL 3 TRIANGLE RULES                                          00028400
C     2)  !M!  <=  J  FOR ALL 3 CASES                                   00028500
C     3)  J  >=  0    FOR ALL 3 CASES                                   00028600
C                                                                       00028700
      IF (NMAX - NMIN) 110, 300, 200                                    00028800
C                                                                       00028900
C     IS IT THE ALL M = 0 CASE                                          00029000
C                                                                       00029100
 200  IF (IABS(X)+IABS(Y) .EQ. 0) GO TO 500                             00029200
C                                                                       00029300
C     PREPARE TO DO THE SUM                                             00029400
      DX = NMAX                                                         00029500
      E1 = -DX                                                          00029600
      E2 = -DB1-DX                                                      00029700
      E3 = -DB2 - DX                                                    00029800
      E4 = DA3 - (DX-1)                                                 00029900
      E5 = DX2 - (DX-1)                                                 00030000
      E6 = DY1 - (DX-1)                                                 00030100
C                                                                       00030200
      F1 = E1*E2*E3                                                     00030300
      F4 = E4*E5*E6                                                     00030400
C                                                                       00030500
      NMAXM2 = NMAX - 2                                                 00030600
      SUM = 1 + F4/F1                                                   00030700
      IF (NMAXM2 .LT. NMIN) GO TO 300                                   00030800
C                                                                       00030900
      G1 = (E1+1)*(E2+1)*(E3+1) - F1                                    00031000
      G4 = (E4+1)*(E5+1)*(E6+1) - F4                                    00031100
      H1 = (E1+2)*(E2+2)*(E3+2) - G1 - G1 - F1                          00031200
      H4 = (E4+2)*(E5+2)*(E6+2) - G4 - G4 - F4                          00031300
C                                                                       00031400
C     COMPUTE THE SUM IN THE RACAH FORMULA                              00031500
C                                                                       00031600
      DO 259 N = NMIN, NMAXM2                                           00031700
      F1 = F1 + G1                                                      00031800
      F4 = F4 + G4                                                      00031900
      G1 = G1 + H1                                                      00032000
      G4 = G4 + H4                                                      00032100
      H1 = H1 + 6                                                       00032200
      H4 = H4 + 6                                                       00032300
      SUM = 1 + SUM*(F4/F1)                                             00032400
 259  CONTINUE                                                          00032500
C                                                                       00032600
C     NOW WE NEED THE FACTORIALS OUTSIDE THE SUM                        00032700
C                                                                       00032800
 300  IF (R .LE. MAXFAC) GO TO 400                                      00032900
C                                                                       00033000
C                                                                       00033100
C    SOME OF THE FACTORIALS WILL BE TOO LARGE TO USE THE TABLE.         00033200
C     WE USE ONE OF TWO METHODS*                                        00033300
C                                                                       00033400
C     1)  FOR MODERATELY SIZED J'S WE USE A TABLE OF                    00033500
C         LOG FACTORIALS.                                               00033600
C                                                                       00033700
C     2)  FOR TRUELY LARGE J'S WE USED THE LOG GAMMA FUNCTION           00033800
C         TO GET THE LOGS OF THE FACTORIALS.                            00033900
C                                                                       00034000
C                                                                       00034100
      IF (R .GT. MAXLF) GO TO 350                                       00034200
C                                                                       00034300
C     THIS IS THE MODERATELY LARGE REGION - USE TABLE                   00034400
C                                                                       00034500
C                                                                       00034600
	E=ALLOC(LLOGFC+1+A3)-ALLOC(LLOGFC+1+NMIN)
	1-ALLOC(LLOGFC+1+A3-NMIN)+ALLOC(LLOGFC+1+A1)
	2-ALLOC(LLOGFC+1+B1+NMIN)-ALLOC(LLOGFC+1+Y1-NMIN)
	3+ALLOC(LLOGFC+1+A2)-ALLOC(LLOGFC+1+B2+NMIN)
	4-ALLOC(LLOGFC+1+X2-NMIN)
	E2=.5*(-ALLOC(LLOGFC+1+A1)-ALLOC(LLOGFC+1+A2)
	1-ALLOC(LLOGFC+1+A3)+ALLOC(LLOGFC+1+X1)+ALLOC(LLOGFC+1+X2)
	2+ALLOC(LLOGFC+1+Y1)+ALLOC(LLOGFC+1+Y2)+ALLOC(LLOGFC+1+Z1)
	3+ALLOC(LLOGFC+1+Z2)-ALLOC(LLOGFC+1+R)+ALLOC(LLOGFC+2+C)
	4-ALLOC(LLOGFC+1+C))
C                                                                       00035300
      IF (DABS(E) .GT. 40) GO TO 330                                    00035400
C                                                                       00035500
C     SUM*E**2  MUST BE AN INTEGER.                                     00035600
C     SINCE E IS NOT VERY LARGE IS IS LIKELY TO BE MEANINGFUL TO        00035700
C     ATTEMPT TO FORCE THE INTEGER RESULT.                              00035800
C                                                                       00035900
CIB      SUM = (SUM*DEXP(E) + DSIGN(ADDER, SUM)) + SHIFT                00036000
      SUM = DINT( SUM*DEXP(E) + DSIGN(.50D0, SUM) )                     00036100
C                                                                       00036200
      ANSWER = DEXP(E2)*SUM                                             00036300
      GO TO 450                                                         00036400
C                                                                       00036500
C     HERE  E  IS LARGE AND AN ATTEMPT TO INTERGERIZE  SUM*EXP(E)  WOULD00036600
C     NOT BE SIGNIFICANT ON A 16 PLACE MACHINE.  INSTEAD WE SPREAD      00036700
C     THE E FACTOR AMOUNGST THE OTHER STUFF IN AN ATTEMPT TO            00036800
C     AVOID OVERFLOWS.                                                  00036900
C                                                                       00037000
 330  E = E + E2                                                        00037100
      ANSWER = DEXP(E)*SUM                                              00037200
      GO TO 450                                                         00037300
C                                                                       00037400
C     FOR LARGE VALUES OF J, WE USE THE LOG GAMMA FUNTION TO DIRECTLY   00037500
C     FIND THE LOG OF THE COMPLETE MESS OF OUTSIDE FACTORIALS.          00037600
C                                                                       00037700
C     HERE WE AVOID INTEGER TO DOUBLE CONVERSIONS                       00037800
C                                                                       00037900
 350  DX1 = DA2 + DA3 - DX2                                             00038000
      DC = C                                                            00038100
      DN = NMIN                                                         00038200
      DZ1 = DX1 + DY1 - DA3                                             00038300
      DY2 = DA1 + DA3 - DY1                                             00038400
      DZ2 = DC - DZ1                                                    00038500
C                                                                       00038600
      E = DLOG(DC+1) - DLGAMA(2+DX1+DY1+DZ2)                            00038700
      DO 369 I = 1, 9                                                   00038800
      E = E + DLGAMA(1+DRAY(I))                                         00038900
 369  CONTINUE                                                          00039000
      E = .5*E                                                          00039100
      DA1 = DB1 + (DN+1)                                                00039200
      DA2 = DB2 + (DN+1)                                                00039300
      DA3 = DA3 - (DN-1)                                                00039400
      DY1 = DY1 - (DN-1)                                                00039500
      DX2 = DX2 - (DN-1)                                                00039600
      DY2 = (DN+1)                                                      00039700
      DO 379 I = 1, 6                                                   00039800
      E = E - DLGAMA(DRAY(I))                                           00039900
 379  CONTINUE                                                          00040000
C                                                                       00040100
      ANSWER = DEXP(E) * SUM                                            00040200
C                                                                       00040300
      GO TO 450                                                         00040400
C                                                                       00040500
C     WE CAN USE THE FACTORIAL TABLES AND NOT WORRY                     00040600
C                                                                       00040700
 400  F1 = FACTBL(1+A1)                                                 00040800
      F2 = FACTBL(1+A2)                                                 00040900
      F3 = FACTBL(1+A3)                                                 00041000
C                                                                       00041100
C     HERE WE FORSCE THE SUM OF PRODUCTS OF 3 BINOMIAL                  00041200
C     COEFFICENTS TO BE AN INTEGER                                      00041300
C                                                                       00041400
      E = (F3/(FACTBL(1+NMIN)*FACTBL(1+A3-NMIN)))  * (F2/(FACTBL(1+X2-NM00041500
     *IN)*FACTBL(1+B2+NMIN)))  * (F1/(FACTBL(1+Y1-NMIN)*FACTBL(1+B1+NMIN00041500
     *)))                                                               00041500
CIB      SUM = (E*SUM*E + DSIGN(ADDER, SUM)) + SHIFT                    00041800
      SUM = DINT( E*SUM*E + DSIGN(.50D0, SUM) )                         00041900
C                                                                       00042000
C    FOR THE FOLLOWING NOTE THAT  R = J1+J2+J3 + 1                      00042100
      ANSWER = SUM * (FACTBL(1+X1)*FACTBL(1+X2)/(F3*F2))  * (FACTBL(1+Y100042200
     *)*FACTBL(1+Y2)/F1)  * (FACTBL(2+C)/FACTBL(1+C))  * (FACTBL(1+Z1)*F00042200
     *ACTBL(1+Z2)/FACTBL(1+R))                                          00042200
C                                                                       00042600
 450  IF (MOD(NMIN,2).NE.0) ANSWER = -ANSWER                            00042700
      GO TO 800                                                         00042800
C                                                                       00042900
C                                                                       00043000
C                                                                       00043100
C     USE SPECIAL CASE FOR M1 = M2 = M3 = 0 SINCE RACAH IS              00043200
C     NOT TRIVIAL FOR THAT CASE                                         00043300
C                                                                       00043400
C                                                                       00043500
C     RESULT IS 0 IF J1+J2+J3 IS ODD                                    00043600
 500  IF (MOD(R,2).EQ.0) GO TO 110                                      00043700
      JSUM = R-1                                                        00043800
      IF (R .GT. MAXFAC) GO TO 550                                      00043900
C                                                                       00044000
C     USE THE FACTORIAL TABLES                                          00044100
C                                                                       00044200
      ANSWER = ((FACTBL(1+A1)*FACTBL(1+A2))/FACTBL(1+R))  * ((FACTBL(2+C00044300
     *)*FACTBL(1+A3))/FACTBL(1+C))  * (FACTBL(1+JSUM/2)/  (FACTBL(1+A1/200044300
     *)*FACTBL(1+A2/2)*FACTBL(1+A3/2)) )**2                             00044300
      GO TO 570                                                         00044700
C                                                                       00044800
C     USE LOG(FACTORIAL)                                                00044900
C                                                                       00045000
C     IN THE FOLLOWING WE USE  SHFTR( ,1)  TO MAKE A /2.                00045100
C                                                                       00045200
 550  IF (R .GT. MAXLF) GO TO 560                                       00045300
	E=.5*(ALLOC(LLOGFC+1+A1)+ALLOC(LLOGFC+1+A2)+
	1ALLOC(LLOGFC+1+A3)+ALLOC(LLOGFC+2+C)-ALLOC(LLOGFC+1+C)
	2-ALLOC(LLOGFC+1+R))-ALLOC(LLOGFC+1+JISHFT(A1,-1))
	3-ALLOC(LLOGFC+1+JISHFT(A2,-1))-ALLOC(LLOGFC+1+JISHFT(A3,-1))
	4+ALLOC(LLOGFC+1+JISHFT(JSUM,-1))
C                                                                       00045800
      GO TO 565                                                         00045900
C                                                                       00046000
 560  DJSUM = DA1 + DA2 + DA3                                           00046100
      E = DLOG(DFLOAT(C+1)) - DLGAMA(DJSUM+2)                           00046200
      E2 = -DLGAMA(.5*DJSUM+1)                                          00046300
      DO 564 I = 1, 3                                                   00046400
      E = E + DLGAMA(1+DRAY(I))                                         00046500
      E2 = E2 + DLGAMA(1 + .5*DRAY(I))                                  00046600
 564  CONTINUE                                                          00046700
      E = .5*E - E2                                                     00046800
C                                                                       00046900
C                                                                       00047000
 565  ANSWER = DEXP(E)                                                  00047100
C                                                                       00047200
C     GET PHASE                                                         00047300
 570  IF (MOD(A3,4).GT.1) ANSWER = -ANSWER                              00047400
C                                                                       00047500
C     DETERMIN THE REQUESTED OUTPUT MODE                                00047600
C                                                                       00047700
C     WAS A 3-J REQUESTED                                               00047800
 800  IF (THRESW) GO TO 950                                             00047900
C                                                                       00048000
 810  CONTINUE                                                          00048100
CXXCIB      IF (DBLSW)  GO TO 850                                       00048200
CXXC                                                                    00048300
CXXCIBC   ROUND THE ANSWER TO SINGLE PRECISION                          00048400
CXXCIBC                                                                 00048500
CXXCIB      RNDEXP = ANSEXP                                             00048600
CXXCIB      IF (ANSWER .EQ. 0)  GOTO 830                                00048700
CXXCIB      COFCG = ANSWER + ROUND                                      00048800
CXXCIB      RETURN                                                      00048900
CXXCIBC    ROUND ZERO RESULTS IN UNDERFLOW                              00049000
CXXCIB 830  COFCG = ANSWER                                              00049100
CXXCIB      RETURN                                                      00049200
C                                                                       00049300
C     DOUBLE PRECISION ENTRY USED                                       00049400
C                                                                       00049500
 850  CLEBSH = ANSWER                                                   00049600
      RETURN                                                            00049700
C                                                                       00049800
C     THIS IS AVAILABLE ONLY IN DOUBLE PRECISION                        00049900
C                                                                       00050000
CIB      ENTRY THREEJ (/A/, /B/, /CIN/, /X/, /Y/, /ZIN/)                00050100
      ENTRY THREEJ ( A, B, ARG3, ARG4, ARG5, ZIN )                      00050200
C                                                                       00050300
      C = ARG3                                                          00050400
      X = ARG4                                                          00050500
      Y = ARG5                                                          00050600
C                                                                       00050700
      ZZ= -ZIN                                                          00050800
      THRESW = .TRUE.                                                   00050900
CXXCIB      DBLSW = .TRUE.                                              00051000
      GO TO 20                                                          00051100
C                                                                       00051200
C     NOW PUT IN THE 3-J FACTOR                                         00051300
C                                                                       00051400
 950  TEMP = A2-Z2                                                      00051500
      TEMP = IABS(TEMP) + 2                                             00051600
      IF (MOD(TEMP,2).NE.0) ANSWER = -ANSWER                            00051700
      IF (C .GT. MAXFAC-1) GO TO 970                                    00051800
      ANSWER = ANSWER * (FACTBL(1+C)/FACTBL(2+C))                       00051900
      GO TO 810                                                         00052000
 970  ANSWER = ANSWER / DSQRT(C+1.D0)                                   00052100
      GO TO 810                                                         00052200
C                                                                       00052300
      END                                                               00052500
      FUNCTION DGAMMA ( X )                                             00000100
C                                                                       00000200
C     COMPLETE GAMMA FUNCTION AND RELATED FUNCTIONS.                    00000300
C                                                                       00000400
C                                                                       00000500
C     ENTRY POINTS (ALL ARE DOUBLE PRECISION):                          00000600
C                                                                       00000700
C     DGAMMA:  COMPLETE GAMMA FUNCTION.                                 00000800
C     DGAMAP:  FIRST DERIVATIVE OF THE COMPLETE GAMMA FUNCTION.         00000900
C     DLGAMA, DLNGAM:  NATURAL LOGARITHM OF DABS( GAMMA(X) ).           00001000
C     DIGAMA:  DIGAMMA, OR PSI, FUNCTION.                               00001100
C                                                                       00001200
C                                                                       00001300
C          ALL OF THESE FUNCTIONS MUST BE DECLARED DOUBLE PRECISION     00001400
C     (OR REAL*8) IN THE CALLING PROGRAM, AND ALL REQUIRE A DOUBLE      00001500
C     PRECISION ARGUMENT.                                               00001600
C                                                                       00001700
C                                                                       00001800
C          REFERENCE:  "A&S" BELOW REFERS TO ABRAMOWITZ AND STEGUN,     00001900
C     HANDBOOK OF MATHEMATICAL FUNCTIONS, PUBLISHED BY THE NATIONAL     00002000
C     BUREAU OF STANDARDS (1964) AND BY DOVER (1965).                   00002100
C                                                                       00002200
C     METHODS:                                                          00002300
C          DGAMMA AND DLGAMA (OR DLNGAM) ARE FOUND USING THE SAME       00002400
C     METHODS (BELOW).  LOGS OR EXPONENTIALS ARE TAKEN AT THE END.      00002500
C          FOR X < O.5 (INCLUDING NEGATIVE X), THE FUNCTION AT 1-X      00002600
C     IS COMPUTED (AS BELOW) AND TRANSFORMED BACK USING A&S 6.1.17.     00002700
C     HENCE WE ASSUME BELOW THAT X >= 0.5.                              00002800
C          FOR X < 12, GAMMA(X-N) IS COMPUTED (AS BELOW), WHERE N IS    00002900
C     THE INTEGER SUCH THAT 0.5 <= X-N < 1.5.  THE RESULT IS TRANSFORMED00003000
C     BACK UP (IF NECESARY) USING A&S 6.1.15.                           00003100
C          FOR 0.5 <= X < 1.5, WE COMPUTE 1/GAMMA(X) USING THE FIRST    00003200
C     19 TERMS OF THE SERIES EXPANSION IN A&S 6.1.34.                   00003300
C          FOR 12 <= X < 164, WE COMPUTE LOG(GAMMA) USING 6 TERMS OF    00003400
C     THE ASYMPTOTIC EXPANSION IN A&S 6.1.40.  FOR X >= 164, ONLY 2     00003500
C     TERMS ARE NECESSARY.                                              00003600
C                                                                       00003700
C          DIGAMA(X) IS FOUND USING A&S 6.3.18, USING 8 TERMS, FOR      00003800
C     \X\ >= 9.0.  FOR \X\ < 9.0, THE ARGUMENT IS TRANSFORMED USING     00003900
C     A&S 6.3.6 TO BRING IT ABOVE 10.0.  FOR X < 0.5, DIGAMA(1-X)       00004000
C     IS FOUND AND TRANSFORMED BACK USING A&S 6.3.7.                    00004100
C                                                                       00004200
C          ACCURACY:  IN THE FOLLOWING, "ERROR" REFERS TO ABSOLUTE      00004300
C     OR RELATIVE ERROR, WHICHEVER IS SMALLER.                          00004400
C          DLGAMA/DLNGAM AND DIGAMA:  FOR X > 0.5, ERRORS IN THESE      00004500
C     FUNCTIONS ARE LESS THAN 1.E-16 ON A VAX.  ABSOLUTE ERRORS FOR     00004600
C     X < 0.5 SHOULD BE COMPARABLE TO THOSE FOR 1-X; HENCE, THEY        00004700
C     WILL BE MUCH WORSE WHEN THE FUNCTION IS MUCH SMALLER AT           00004800
C     X (<0.5) THAN AT 1-X.                                             00004900
C          DGAMMA:  FOR -11 < X < 12, ERRORS ARE LESS THAN 1.E-16.      00005000
C     FOR X >= 12 AND X <= 11, ACCURACY IS LIMITED BY THE DEXP FUNCTION,00005100
C     BECAUSE TO FIRST ORDER IN D WE HAVE                               00005200
C            EXP (1+D)*Y  = EXP(Y) * (1 + Y*D).                         00005300
C     HENCE, THE RELATIVE ERROR IN DGAMMA(X) IS EXPECTED TO BE DLGAMA(X)00005400
C     TIMES THE RELATIVE ERROR IN DLGAMA(X).  THIS EXPECTATION IS       00005500
C     VERIFIED BY THE TESTS.  RELATIVE ERRORS IN DGAMMA ARE ABOUT 3.E-1600005600
C     AT X = 12, AND CLIMB TO 8.E-16 AT X = 19.                         00005700
C          DGAMAP:  ERRORS ARE COMPARABLE TO DGAMMA.                    00005800
C                                                                       00005900
C          ERRORS:  IF THE ARGUMENT IS TOO BIG OR TOO NEAR A NEGATIVE   00006000
C     INTEGER OR ZERO, AN APPROPRIATE ERROR MESSAGE IS PRINTED, AND THEN00006100
C     EXECUTION IS ALLOWED TO CONTINUE.  IN MOST CASES, THE PROGRAM WILL00006200
C     SUBSEQUENTLY BOMB WITH A DIVIDE BY ZERO OR AN OVERFLOW IN THE DEXP00006300
C     FUNCTION.  THE MAXIMUM ARGUMENT FOR DGAMMA OR DGAMAP IS ABOUT 34  00006400
C     FOR A VAX, 57 FOR IBM, AND 300 FOR CDC.                           00006500
C                                                                       00006600
C     1/26/80 - ROBERT P. GODDARD                                       00006700
C                                                                       00006800
C                                                                       00006900
      IMPLICIT REAL*8 ( A-H, O-Z )                                      00007000
C                                                                       00007100
      LOGICAL LNEG, LEXP                                                00007200
      CHARACTER*6 NAME(4)                                               00007300
      DATA NAME / 'DGAMMA', 'DGAMAP', 'DLGAMA', 'DIGAMA' /              00007400
      DATA PI / 3.1415 92653 58979 324 0D0 /                            00007500
      DATA R2PILN / 0.91893 85332 04672 742 0D0 /                       00007600
C                                                                       00007700
C     ELIM IS THE BIGGEST NUMBER S.T. DEXP(ELIM) DOESN'T OVERFLOW.      00007800
C                                                                       00007900
      DATA ELIM / 88.029 0D0 /                                          00008000
CIB      DATA  ELIM / 176.058 0D0 /                                     00008100
CUN      DATA  ELIM / 352.116 0D0 /                                     00008200
CDC      DATA  ELIM / 704.23 /                                          00008300
CRA      DATA  ELIM / 5633.8 /                                          00008400
C                                                                       00008500
CDC      NOW                                                            00008600
CFR%F                                                                   00008700
C                                                                       00008800
C - - - - - - - - - - DGAMMA - - - - - - - - - -                        00008900
C                                                                       00009000
      NE = 1                                                            00009100
      GO TO 100                                                         00009200
C                                                                       00009300
 10   DGAMMA = G                                                        00009400
      RETURN                                                            00009500
C                                                                       00009600
C - - - - - - - - - - DGAMAP - - - - - - - - - -                        00009700
C                                                                       00009800
      ENTRY DGAMAP (X)                                                  00009900
C                                                                       00010000
      NE = 2                                                            00010100
      GO TO 100                                                         00010200
C                                                                       00010300
 20   DGAMAP = G * D                                                    00010400
      RETURN                                                            00010500
C                                                                       00010600
C - - - - - - - - DLGAMA, DLNGAM - - - - - - - -                        00010700
C                                                                       00010800
      ENTRY DLGAMA ( X )                                                00010900
      ENTRY DLNGAM ( X )                                                00011000
C                                                                       00011100
      NE = 3                                                            00011200
      GO TO 100                                                         00011300
C                                                                       00011400
 30   IF ( LEXP ) G = DLOG( G )                                         00011500
      DLGAMA = G                                                        00011600
      RETURN                                                            00011700
C                                                                       00011800
C - - - - - - - - - - DIGAMA - - - - - - - - - -                        00011900
C                                                                       00012000
      ENTRY DIGAMA ( X )                                                00012100
C                                                                       00012200
      NE = 4                                                            00012300
      GO TO 500                                                         00012400
C                                                                       00012500
 40   DIGAMA = D                                                        00012600
      RETURN                                                            00012700
C                                                                       00012800
C                                                                       00012900
C     LOG GAMMA CODE COMMON TO DGAMMA, DLGAMA, AND DLNGAM.              00013000
C                                                                       00013100
C                                                                       00013200
C     CHECK FOR ARGUMENT < 0.5 AND TRANSFORM.                           00013300
C                                                                       00013400
 100  Z = X                                                             00013500
      LEXP = .FALSE.                                                    00013600
      LNEG = Z .LT. 0.50D0                                              00013700
      IF ( LNEG ) Z = 1.D0 - Z                                          00013800
C                                                                       00013900
C     IF Z < 12, WE TRANSFORM IT DOWN INTO THE RANGE 0.5 <= Z < 1.5     00014000
C     (A&S 6.1.15) AND USE A 19-TERM SERIES EXPANSION OF 1/GAMMA(1+Z)   00014100
C     (FROM A&S 6.1.34).                                                00014200
C                                                                       00014300
      IF ( Z .GE. 12.D0 ) GO TO 200                                     00014400
      H = 1.D0                                                          00014500
      N = Z - 0.50D0                                                    00014600
      IF ( N .EQ. 0 ) GO TO 120                                         00014700
      DO 110 I = 1, N                                                   00014800
      Z = Z - 1.D0                                                      00014900
      H = H * Z                                                         00015000
 110  CONTINUE                                                          00015100
C                                                                       00015200
 120  Y = Z - 1.D0                                                      00015300
      G = -.00002 01348 547807 0D0                                      00015400
      IF ( DABS( Y ) .GT. 0.10D0 ) G = G  + Y * ( -.00000 12504 934821 000015500
     *D0  + Y * ( 0.00000 11330 272320 0D0  + Y * ( -.00000 02056 33841700015500
     * 0D0  + Y * ( 0.00000 00061 160950 0D0  + Y * ( 0.00000 00050 020000015500
     *75 0D0  + Y * ( -.00000 00011 812746 0D0  + Y * ( 0.00000 00001 0400015500
     *3427 0D0  ) ) ) ) ) ) )                                           00015500
      G = 1.0D0  + Y * ( 0.57721 56649 015329 0D0  + Y * ( -.65587 8071500016400
     * 202538 0D0  + Y * ( -.04200 26350 340952 0D0  + Y * ( 0.16653 86100016400
     *13 822915 0D0  + Y * ( -.04219 77345 555443 0D0  + Y * ( -.00962 100016400
     *9715 278770 0D0  + Y * ( 0.00721 89432 466630 0D0  + Y * ( -.0011600016400
     * 51675 918591 0D0  + Y * ( -.00021 52416 741149 0D0  + Y * ( 0.00000016400
     *12 80502 823882 0D0  + Y * G ) ) ) ) ) ) ) ) ) )                  00016400
      G = H / G                                                         00017600
      LEXP = .TRUE.                                                     00017700
      GO TO 250                                                         00017800
C                                                                       00017900
C     CALCULATE LOG GAMMA USING 6 TERMS OF THE BERNOULLI EXPANSION.     00018000
C     FOR \X\ > 164, WE CAN GET 17 DECIMAL PLACES WITH 2 TERMS.         00018100
C     THE CONSTANT EXPRESSIONS FOR THE BERNOULLI NUMBERS  DIVIDED       00018200
C     BY (2*M)*(2*M-1)  ARE EVALUATED AT COMPILE TIME.  THE BERNOULLI   00018300
C     NUMBERS COME FROM A&S TABLE 23.2.                                 00018400
C                                                                       00018500
 200  ZI = 1.D0 / Z                                                     00018600
      ZIS = ZI * ZI                                                     00018700
      G = (-1.D0/360.D0)                                                00018800
      IF ( Z .LT. 164.D0 ) G = G  + ZIS * ( (1.D0/1260.D0)  + ZIS * ( (-00018900
     *1.D0/1680.D0)  + ZIS * ( (1.D0/1188.D0)  + ZIS * ( (-691.D0/31284000018900
     *.D0)  ) ) ) )                                                     00018900
      G = R2PILN - Z + ( Z - 0.50D0 ) * DLOG( Z )  + ZI * ( (1.D0/12.D0)00019500
     *  + ZIS * G )                                                     00019500
C                                                                       00019800
C      IF X < 0.5, TRANSFORM THE RESULT USING A&S 6.1.17.               00019900
C                                                                       00020000
 250  IF ( .NOT. LNEG ) GO TO 300                                       00020100
      T = DSIN( PI*Z )                                                  00020200
      IF ( T .NE. 0.D0 ) GO TO 260                                      00020300
      PRINT 803, NAME(NE), X                                            00020400
      PRINT 813                                                         00020500
 260  IF ( LEXP ) GO TO 270                                             00020600
      G = DLOG( PI / DABS( T ) ) - G                                    00020700
      LNEG = T .LT. 0                                                   00020800
      GO TO 300                                                         00020900
 270  CONTINUE                                                          00021000
      G = PI / ( G * T )                                                00021100
      LNEG = .FALSE.                                                    00021200
C                                                                       00021300
C     WE NOW HAVE LOG GAMMA (OR GAMMA IF LEXP = .TRUE.).                00021400
C     BRANCH DEPENDS ON THE ENTRY POINT USED.                           00021500
C                                                                       00021600
 300  IF ( NE .EQ. 3 ) GO TO 30                                         00021700
      IF ( LEXP ) GO TO ( 10, 500 ), NE                                 00021800
      IF ( G .LE. ELIM ) GO TO 350                                      00021900
      PRINT 803, NAME(NE), X                                            00022000
      PRINT 913                                                         00022100
 350  G = DEXP( G )                                                     00022200
      IF ( LNEG ) G = -G                                                00022300
      IF ( NE .EQ. 1 ) GO TO 10                                         00022400
C                                                                       00022500
C                                                                       00022600
C     DIGAMMA CODE COMMON TO DGAMAP, DIGAMMA.                           00022700
C                                                                       00022800
C     CHECK FOR ARGUMENT < 0.5 AND TRANSFORM.                           00022900
C                                                                       00023000
 500  Z = X                                                             00023100
      LNEG = Z .LT. 0.50D0                                              00023200
      IF ( LNEG ) Z = 1.D0 - Z                                          00023300
C                                                                       00023400
C     TRANSFORM TO Z >= 9 IF NECESSARY.                                 00023500
C                                                                       00023600
 550  H = 0.D0                                                          00023700
      IF ( Z .GE. 9.D0 ) GO TO 600                                      00023800
      N = 10.D0 - Z                                                     00023900
      DO 580 I = 1, N                                                   00024000
      H = H + 1.D0/Z                                                    00024100
      Z = Z + 1.D0                                                      00024200
 580  CONTINUE                                                          00024300
C                                                                       00024400
C     CALCULATE DIGAMMA USING 8 TERMS OF THE BERNOULLI EXPANSION.       00024500
C     THE CONSTANT EXPRESSIONS FOR THE BERNOULLI NUMBERS (DIVIDED       00024600
C     BY 2*N) ARE EVALUATED AT COMPILE TIME.                            00024700
C                                                                       00024800
 600  ZI = 1.D0/Z                                                       00024900
      ZIS = ZI*ZI                                                       00025000
      D = DLOG(Z) - 0.50D0*ZI - H  - ZIS * ( (1.D0/12.D0)  + ZIS * ( (-100025100
     *.D0/120.D0)  + ZIS * ( (1.D0/252.D0)  + ZIS * ( (-1.D0/240.D0)  + 00025100
     *ZIS * ( (1.D0/132.D0)  + ZIS * ( (-691.D0/32760.D0)  + ZIS * ( (1.00025100
     *D0/12.D0)  + ZIS * ( (-3617.D0/8160.D0)  ) ) ) ) ) ) ) )          00025100
C                                                                       00026100
C     MODIFY THE RESULT IF THE ARGUMENT WAS NEGATIVE.                   00026200
C                                                                       00026300
      IF ( .NOT. LNEG ) GO TO 700                                       00026400
      T = DSIN( PI*Z )                                                  00026500
      IF ( T .NE. 0.D0 ) GO TO 690                                      00026600
      PRINT 803, NAME(NE), X                                            00026700
      PRINT 813                                                         00026800
 690  D = D + PI*DCOS( PI*Z )/T                                         00026900
 700  IF ( NE .EQ. 4 ) GO TO 40                                         00027000
      GO TO 20                                                          00027100
C                                                                       00027200
C     ERROR MESSAGES.                                                   00027300
C                                                                       00027400
 803  FORMAT ( '0**** ERROR IN ', A6, ':  ARGUMENT =', G25.16 ' ****' ) 00027500
 813  FORMAT ( ' **** ARGUMENT IS A NEGATIVE INTEGER OR ZERO ****' )    00027600
 913  FORMAT ( ' **** EXPONENT OVERFLOW ****' )                         00027700
      END                                                               00027800
      SUBROUTINE DSGMAL(ETA,NO,DSG)                                     00000100
C                                                                       00000110
C                                                                       00000120
C     COULOMB PHASES                                                    00000130
C                                                                       00000140
C     ?/?/? - ORIGINAL VERSION                                          00000150
C     12/28/79 - CDC TO CNI, TITLE COMMENT - RPG                        00000160
C     7/9/80 - USE CDM PREFIX - S.P.                                    00000170
C                                                                       00000180
      IMPLICIT REAL*8 ( A-H, O-Z )                                      00000200
CDC      ALLOCLEVEL DSG                                                 00000300
      DIMENSION P(26),Q(23),DSG(1)                                      00000400
      DATA P( 1) / 0.310165781012994887 0D-08 /                         00000600
CIB      DATA P( 1) /Z39D524F0F53B0885/                                 00000700
      DATA P( 2) / 0.301594747899910910 0D-05 /                         00000800
CIB      DATA P( 2) /Z3C32996552696E14/                                 00000900
      DATA P( 3) / 0.205644287153878958 0D-03 /                         00001000
CIB      DATA P( 3) /Z3DD7A237CE95149F/                                 00001100
      DATA P( 4) / 0.444558861472056296 0D-02 /                         00001200
CIB      DATA P( 4) /Z3F1235899B631A0C/                                 00001300
      DATA P( 5) / 0.424137251918122321 0D-01 /                         00001400
CIB      DATA P( 5) /Z3FADBA03A99B1830/                                 00001500
      DATA P( 6) / 0.206378868197296478 0D+00 /                         00001600
CIB      DATA P( 6) /Z4034D53ED97E2D92/                                 00001700
      DATA P( 7) / 0.546892269520120738 0D+00 /                         00001800
CIB      DATA P( 7) /Z408C0121BC062DCA/                                 00001900
      DATA P( 8) / 0.794948065779220198 0D+00 /                         00002000
CIB      DATA P( 8) /Z40CB81B7688A4B0A/                                 00002100
      DATA P( 9) / 0.593131560870811980 0D+00 /                         00002200
CIB      DATA P( 9) /Z4097D778502A62D8/                                 00002300
      DATA P(10) / 0.177060007536960065 0D+00 /                         00002400
CIB      DATA P(10) /Z402D53CDFDCCFE69/                                 00002500
      DATA P(11) / 0.241087856973594357 0D-08 /                         00002600
CIB      DATA P(11) /Z39A5AC9FD40B5BA8/                                 00002700
      DATA P(12) / 0.932656994995955490 0D-06 /                         00002800
CIB      DATA P(12) /Z3BFA5BB35F20D66A/                                 00002900
      DATA P(13) / 0.758581379314270606 0D-04 /                         00003000
CIB      DATA P(13) /Z3D4F8B038B78FE9E/                                 00003100
      DATA P(14) / 0.200706718860656989 0D-02 /                         00003200
CIB      DATA P(14) /Z3E8388FFEF99C18E/                                 00003300
      DATA P(15) / 0.188924467702797066 0D-01 /                         00003400
CIB      DATA P(15) /Z3F4D622A9050FC9E/                                 00003500
      DATA P(16) / 0.444315921800891838 0D-01 /                         00003600
CIB      DATA P(16) /Z3FB5FDE6B529BF7E/                                 00003700
      DATA P(17) / -0.159601618385062274 0D+00 /                        00003800
CIB      DATA P(17) /ZC028DBA6D35A3E69/                                 00003900
      DATA P(18) / -0.685141773275969618 0D+00 /                        00004000
CIB      DATA P(18) /ZC0AF65738557FF1A/                                 00004100
      DATA P(19) / -0.577201920703612828 0D+00 /                        00004200
CIB      DATA P(19) /ZC093C3814C9C4355/                                 00004300
      DATA P(20) / 0.139541051788112899 0D+03 /                         00004400
CIB      DATA P(20) /Z428B8A825EB76321/                                 00004500
      DATA P(21) / 0.109052269358000365 0D+04 /                         00004600
CIB      DATA P(21) /Z4344285CF3F17F1E/                                 00004700
      DATA P(22) / 0.485358524121951234 0D+03 /                         00004800
CIB      DATA P(22) /Z431E55BC83CA29B9/                                 00004900
      DATA P(23) / -0.658758053903885070 0D+03 /                        00005000
CIB      DATA P(23) /ZC3292C20FD215CAA/                                 00005100
      DATA P(24) / -0.413023138385032667 0D+02 /                        00005200
CIB      DATA P(24) /ZC2294D6470917FF0/                                 00005300
      DATA P(25) / 0.922933578234238228 0D+01 /                         00005400
CIB      DATA P(25) /Z4193AB5BFF4F68E2/                                 00005500
      DATA P(26) / -0.999999999999999944 0D+00 /                        00005600
CIB      DATA P(26) /ZC0FFFFFFFFFFFFFC/                                 00005700
      DATA Q( 1) / 0.162632093394676580 0D-05 /                         00005800
CIB      DATA Q( 1) /Z3C1B48FEC701B849/                                 00005900
      DATA Q( 2) / 0.180858599479503871 0D-03 /                         00006000
CIB      DATA Q( 2) /Z3DBDA4DC51C6A4B7/                                 00006100
      DATA Q( 3) / 0.577143957513920856 0D-02 /                         00006200
CIB      DATA Q( 3) /Z3F17A3CB039DF446/                                 00006300
      DATA Q( 4) / 0.787659072077549621 0D-01 /                         00006400
CIB      DATA Q( 4) /Z40142A00A37F4085/                                 00006500
      DATA Q( 5) / 0.545157981064667799 0D+00 /                         00006600
CIB      DATA Q( 5) /Z408B8F7933D37D9D/                                 00006700
      DATA Q( 6) / 0.208553610110275223 0D+01 /                         00006800
CIB      DATA Q( 6) /Z41215E5B1A4DD198/                                 00006900
      DATA Q( 7) / 0.457313299223146097 0D+01 /                         00007000
CIB      DATA Q( 7) /Z41492B8D801E48B1/                                 00007100
      DATA Q( 8) / 0.569717859194618481 0D+01 /                         00007200
CIB      DATA Q( 8) /Z415B27A4BD3E1545/                                 00007300
      DATA Q( 9) / 0.373731134057316638 0D+01 /                         00007400
CIB      DATA Q( 9) /Z413BCC06F9EBB4E9/                                 00007500
      DATA Q(10) / 0.841871279655082615 0D-09 /                         00007600
CIB      DATA Q(10) /Z3939DA5B2E6CDBD7/                                 00007700
      DATA Q(11) / 0.525199068334112846 0D-06 /                         00007800
CIB      DATA Q(11) /Z3B8CFB67B871E12F/                                 00007900
      DATA Q(12) / 0.638793895356789098 0D-04 /                         00008000
CIB      DATA Q(12) /Z3D42FB7E391D89CE/                                 00008100
      DATA Q(13) / 0.262822889969226652 0D-02 /                         00008200
CIB      DATA Q(13) /Z3EAC3E5D2BAB275D/                                 00008300
      DATA Q(14) / 0.448275937331078468 0D-01 /                         00008400
CIB      DATA Q(14) /Z3FB79D2390AC547B/                                 00008500
      DATA Q(15) / 0.345389220533413407 0D+00 /                         00008600
CIB      DATA Q(15) /Z40586B6D8E94FA2B/                                 00008700
      DATA Q(16) / 0.122326029540586134 0D+01 /                         00008800
CIB      DATA Q(16) /Z4113927963343730/                                 00008900
      DATA Q(17) / 0.188103263084231598 0D+01 /                         00009000
CIB      DATA Q(17) /Z411E18B5AC02D34E/                                 00009100
      DATA Q(18) / -0.141115652526220753 0D+03 /                        00009200
CIB      DATA Q(18) /ZC28D1D9B6769D163/                                 00009300
      DATA Q(19) / -0.996536251962568770 0D+03 /                        00009400
CIB      DATA Q(19) /ZC33E48947CF01A61/                                 00009500
      DATA Q(20) / -0.535816773446634613 0D+03 /                        00009600
CIB      DATA Q(20) /ZC3217D1181089897/                                 00009700
      DATA Q(21) / 0.654427162588470765 0D+03 /                         00009800
CIB      DATA Q(21) /Z4328E6D5A87038E8/                                 00009900
      DATA Q(22) / 0.421589251537016345 0D+02 /                         00010000
CIB      DATA Q(22) /Z422A28AF51A1A909/                                 00010100
      DATA Q(23) / -0.931266911567572087 0D+01 /                        00010200
CIB      DATA Q(23) /ZC19500B154A4BE4F/                                 00010300
      DATA EZ / 0.180554707160510697 0D+01 /                            00010400
CIB      DATA EZ    /Z411CE385537EEB8A/                                 00010500
      DATA EZ1 / 0.180554771423339844 0D+01 /                           00010600
CIB      DATA EZ1   /Z411CE38600000000/                                 00010700
      DATA EZ2 / 0.642628291517623633 0D-06 /                           00010800
CIB      DATA EZ2   /Z3BAC811476376D16/                                 00010900
C                                                                       00011000
      X=DABS(ETA)                                                       00011200
      XSQ=X*X                                                           00011300
      IF(X.GT.4.0) GO TO 200                                            00011400
      IF(X.GT.2.0) GO TO 100                                            00011500
C                                                                       00011600
C     ABS(ETA) IN (0,2) RANGE SO=ETA(ETA**2-EZSQ)R(ETA**2)              00011700
C                                                                       00011800
      I=7                                                               00011900
      J=1                                                               00012000
      K=1                                                               00012100
      M=1                                                               00012200
      RSQ=XSQ                                                           00012300
      GO TO 300                                                         00012400
C                                                                       00012500
C     ABS(ETA) IN (2,4) RANGE SO=ETA*R(ETA**2)                          00012600
C                                                                       00012700
  100 I=6                                                               00012800
      J=11                                                              00012900
      K=10                                                              00013000
      M=2                                                               00013100
      RSQ=XSQ                                                           00013200
      GO TO 300                                                         00013300
C                                                                       00013400
C     ABS(ETA) IN (4,-) RANGE SO=ATAN(ETA)/2+ETA(LN(1+ETA**2)/2         00013500
C                                +R(1/ETA**2))                          00013600
C                                                                       00013700
  200 I=4                                                               00013800
      J=20                                                              00013900
      K=18                                                              00014000
      M=3                                                               00014100
      RSQ=1.0D0/XSQ                                                     00014200
  300 QSUM=Q(K)                                                         00014300
      N=K+I                                                             00014400
      DO 305 L=K,N                                                      00014500
  305 QSUM=QSUM*RSQ+Q(L+1)                                              00014600
      PSUM=P(J)                                                         00014700
      N=J+I+1                                                           00014800
      DO 310 L=J,N                                                      00014900
  310 PSUM=PSUM*RSQ+P(L+1)                                              00015000
      R=PSUM/(QSUM*RSQ+1.0D0)                                           00015100
      GO TO (400,500,600),M                                             00015200
  400 DSG(1)=X*R*(X+EZ)*((X-EZ1)+EZ2)                                   00015300
      GO TO 610                                                         00015400
  500 DSG(1)=X*R                                                        00015500
      GO TO 610                                                         00015600
  600 DSG(1)=DATAN(X)/2.0D0+X*(DLOG(1.0D0+XSQ)/2.0D0+R)                 00015700
C                                                                       00015800
C     SL(ETA)=SO(ETA)+SUM(I=1TOL)DATAN(ETA/I)                           00015900
C                                                                       00016000
  610 IF(NO.EQ.0) GO TO 650                                             00016100
      R=1.0D0                                                           00016200
      DO 620 I=1,NO                                                     00016300
      RSQ=X/R                                                           00016400
      DSG(I+1)=DSG(I)+DATAN(RSQ)                                        00016500
  620 R=R+1.0D0                                                         00016600
      IF (ETA.GE.0.0D0) GOTO 660                                        00016700
      DO 640 I=1,NO                                                     00016800
  640 DSG(I+1)=-DSG(I+1)                                                00016900
  650 IF (ETA.LT.0.0D0) DSG(1)=-DSG(1)                                  00017000
  660 RETURN                                                            00017100
      END                                                               00017200
      SUBROUTINE EDTDAT(ILL,NOCHS,INPUT)                                00000010
C                                                                       00000020
C         RETURNS THE DATE IN THE FORM                                  00000030
C          January 9, 1975                                              00000040
C       INPUT IS THE 8 BYTE CHARACTER STRING  ' 9/08/75'                00000050
C         ILL IS THE CHARACTER STRING ARRAY INTO WHICH                  00000060
C             THE DATE IS PUT.  IT MUST BE ATLEAST 18 BYTES LONG.       00000070
C             THE FIRST 18 BYTES OF ILL ARE SET TO BLANKS AND THE       00000080
C             DATE IS RIGHT JUSTIFIED IN THIS FIELD.                    00000090
C          NOCHS IS THE NUMBER OF CHARACTERS IN THE REP.                00000100
C                                                                       00000110
C       S.C   1/75                                                      00000120
C     12/28/79 - CIB TO CND, PUT BACK LOWER-CASE LETTERS - RPG          00000130
C     1/16/80 - USE CHARTONUM, BL, ZRO, ONE, NINE; SHORTEN DATA - RPG   00000140
C     1/21/80 - REMOVE TRNSFR CALL - RPG                                00000150
C                                                                       00000160
      IMPLICIT CHARACTER*8 (A-H,O-Z)                                    00000170
      CHARACTER*1 ILL(1),INPUT(1)                                       00000180
      CHARACTER*1 COM                                                   00000190
        data COM / ',' /                                                00000190
      CHARACTER*1 IL(8)                                                 00000200
      CHARACTER*8 RMO(10)                                               00000210
        data RMO/'JanuaryF','ebruaryM','archApri','lMayJune' ,'JulyAugu'00000210
     *,'stSeptem','berOctob','erNovemb','erDecemb','er'/                00000210
      INTEGER*4 IBEG(12)                                                00000230
        data IBEG/1,8,16,21,26,29,33,37,43,52,59,67/                    00000230
      INTEGER*4 NOCH(12)                                                00000240
        data NOCH/7,8,5,5,3,4,4,6,9,7,8,8/                              00000240
      CHARACTER*1 BL, ZRO, ONE, NINE                                    00000250
        dataNINE/'9'/                                                   00000250
        dataONE/'1'/                                                    00000250
        dataZRO/'0'/                                                    00000250
        data BL / ' ' /                                                 00000250
      CHARACTER*8 AL                                                    00000260
      CHARACTER*1 MO(1)                                                 00000270
      EQUIVALENCE (MO(1),RMO(1))                                        00000280
      EQUIVALENCE(IL(1),AL)                                             00000290
      CHARACTER*16 EM                                                   00000300
        data EM / 'INVALID DATE    ' /                                  00000300
      CHARACTER*1 EM1(18)                                               00000310
      EQUIVALENCE ( EM, EM1(1) )                                        00000320
      DO 1 I=1,8                                                        00000330
      IL(I)=INPUT(I)                                                    00000340
1     CONTINUE                                                          00000350
      DO 2 I=1,18                                                       00000360
      ILL(I)=BL                                                         00000370
2     CONTINUE                                                          00000380
100   IMO=ICHAR(IL(2)) - 48                                             00000390
      IF(IL(1).NE.BL) IMO=IMO+10*(ICHAR(IL(1)) - 48)                    00000400
300   M=1                                                               00000410
      IF(IL(4).NE.BL) M=2                                               00000420
      IF(IL(4).EQ.ZRO) M=1                                              00000430
      IF ( IMO .GT. 12 .OR. IMO .LT. 1 ) GO TO 9000                     00000440
      NOCHS=M+NOCH(IMO)+7                                               00000450
      N=NOCH(IMO)                                                       00000460
      DO 310 I=1,N                                                      00000470
      ILL(I)=MO(IBEG(IMO)+I-1)                                          00000480
310   CONTINUE                                                          00000490
      ILL(N+1)=BL                                                       00000500
      IF(M.EQ.1) GO TO 311                                              00000510
      ILL(N+2)=IL(4)                                                    00000520
      ILL(N+3)=IL(5)                                                    00000530
      GO TO 312                                                         00000540
311   ILL(N+2)=IL(5)                                                    00000550
312   ILL(N+M+2)=COM                                                    00000560
      ILL(N+M+3)=BL                                                     00000570
      ILL(M+N+4)=ONE                                                    00000580
      ILL(M+N+5)=NINE                                                   00000590
      ILL(M+N+6)=IL(7)                                                  00000600
      ILL(M+N+7)=IL(8)                                                  00000610
      RETURN                                                            00000620
C                                                                       00000630
C     BAD INPUT DATE                                                    00000640
C                                                                       00000650
 9000 DO 9009 I = 1, 16                                                 00000660
      ILL(I) = EM1(I)                                                   00000670
 9009 CONTINUE                                                          00000680
      NOCHS = 12                                                        00000690
C                                                                       00000700
CDC      CHARACTER*1  ILL(18)                                           00000710
CDC      INTEGER INPUT                                                  00000720
CDC      DECODE ( 10, 13, INPUT )  ILL                                  00000730
CDC 13   FORMAT ( 10A1 )                                                00000740
CDC      DO 29  I = 11, 18                                              00000750
CDC  29  ILL(I) = 1H                                                    00000760
C                                                                       00000770
      RETURN                                                            00000780
      END                                                               00000790
      BLOCK DATA FACDUM                                                 00000100
CNV      BLOCK DATA FACDUM                                              00000200
C                                                                       00000300
C     BLOCK DATA TO DEFINE THE SQRT(FACTORIAL) TABLE                    00000400
C                                                                       00000500
C     THE TABLE IS USED BY THE CLEBSH, THREEJ, SIXJ AND RACAH           00000600
C     ROUTINES.                                                         00000700
C                                                                       00000800
C     IT MAY BE INCLUDED IN YOUR LOAD MODULE FOR OTHER PURPOSES         00000900
C     BY PLACING A DUMMY EXTERNAL CARD OF THE FORM                      00001000
C          EXTERNAL  FACDUM                                             00001100
C     IN YOUR SOURCE.                                                   00001200
C                                                                       00001300
C     THE FORM OF THE TABLE IS                                          00001400
C                                                                       00001500
C     COMMON / FACTRL /  MAX, ISPACE, TABLE                             00001600
C                                                                       00001700
C                                                                       00001800
C     MAX IS THE MAXIMUM NUMBER WHOSE SQRT(FACTORIAL) IS                00001900
C         IN TABLE.  TABLE HAS MAX+1 ELEMENTS.                          00002000
C         IN THE PRESENT VERSION:                                       00002100
C            MAX = 96    (DETERMINED BY EXPONENT RANGE)                 00002200
C     ISPACE IS NOT USED.                                               00002300
C     TABLE IS THE TABLE.  THE FIRST ELEMENT CORRESPONDS TO 0.          00002400
C                                                                       00002500
C     6/30/74 - FIRST VERSION - S. PIEPER                               00002600
C     12/28/79 - CDC TO CNI, LOWER LIMIT FOR VAX - RPG                  00002610
C     7/9/80 - USE CDM PREFIX - S.P.                                    00002620
C                                                                       00002700
C                                                                       00002800
      REAL*8 FACTBL                                                     00002900
C                                                                       00003000
CNV      COMMON  /FACTRL/  MAXFAC, SP1, FACTBL(97)                      00003100
      COMMON /FACTRL/ MAXFAC, SP1, FACTBL(57)                           00003110
      COMMON /XFACDUM/ I                                                 00003200
C                                                                       00003300
C                                                                       00003500
CNV      DATA  MAXFAC /96/                                              00003600
      DATA MAXFAC /56/                                                  00003610
C                                                                       00003700
C     THESE ARE THE SQUARE ROOTS OF FACTORIALS.                         00003800
C                                                                       00003900
C     THEY ARE STORED UP TO THE EXPONENT RANGE OF THE MACHINE.          00004000
C                                                                       00004100
      DATA FACTBL( 1) / 0.100000000000000000 0D+01 /                    00004200
CIB      DATA FACTBL( 1) / Z 4110000000000000 /                         00004300
      DATA FACTBL( 2) / 0.100000000000000000 0D+01 /                    00004400
CIB      DATA FACTBL( 2) / Z 4110000000000000 /                         00004500
      DATA FACTBL( 3) / 0.141421356237309515 0D+01 /                    00004600
CIB      DATA FACTBL( 3) / Z 4116A09E667F3BCD /                         00004700
      DATA FACTBL( 4) / 0.244948974278317810 0D+01 /                    00004800
CIB      DATA FACTBL( 4) / Z 4127311C2812425D /                         00004900
      DATA FACTBL( 5) / 0.489897948556635621 0D+01 /                    00005000
CIB      DATA FACTBL( 5) / Z 414E6238502484BA /                         00005100
      DATA FACTBL( 6) / 0.109544511501033224 0D+02 /                    00005200
CIB      DATA FACTBL( 6) / Z 41AF456E91B52C78 /                         00005300
      DATA FACTBL( 7) / 0.268328157299974777 0D+02 /                    00005400
CIB      DATA FACTBL( 7) / Z 421AD5336963EEFC /                         00005500
      DATA FACTBL( 8) / 0.709929573971953936 0D+02 /                    00005600
CIB      DATA FACTBL( 8) / Z 4246FE3274BB4687 /                         00005700
      DATA FACTBL( 9) / 0.200798406368178131 0D+03 /                    00005800
CIB      DATA FACTBL( 9) / Z 42C8CC645C183E43 /                         00005900
      DATA FACTBL(10) / 0.602395219104534419 0D+03 /                    00006000
CIB      DATA FACTBL(10) / Z 4325A652D1448BAD /                         00006100
      DATA FACTBL(11) / 0.190494094396650524 0D+04 /                    00006200
CIB      DATA FACTBL(11) / Z 43770F0E1B42B823 /                         00006300
      DATA FACTBL(12) / 0.631797435892232807 0D+04 /                    00006400
CIB      DATA FACTBL(12) / Z 4418ADF96F9619F7 /                         00006500
      DATA FACTBL(13) / 0.218861051811417556 0D+05 /                    00006600
CIB      DATA FACTBL(13) / Z 44557E1AED26BBFF /                         00006700
      DATA FACTBL(14) / 0.789114744508046861 0D+05 /                    00006800
CIB      DATA FACTBL(14) / Z 451343F79759BA1B /                         00006900
      DATA FACTBL(15) / 0.295259701280076479 0D+06 /                    00007000
CIB      DATA FACTBL(15) / Z 454815BB3871751D /                         00007100
      DATA FACTBL(16) / 0.114353590586391301 0D+07 /                    00007200
CIB      DATA FACTBL(16) / Z 461172EFE7E6B289 /                         00007300
      DATA FACTBL(17) / 0.457414362345565180 0D+07 /                    00007400
CIB      DATA FACTBL(17) / Z 4645CBBF9F9ACA23 /                         00007500
      DATA FACTBL(18) / 0.188596773062531464 0D+08 /                    00007600
CIB      DATA FACTBL(18) / Z 4711FC69D4E669B3 /                         00007700

      DATA FACTBL(19) / 0.800148342854498439 0D+08 /                    00007800
CIB      DATA FACTBL(19) / Z 474C4EDF249133DB /                         00007900
      DATA FACTBL(20) / 0.348776576634429395 0D+09 /                    00008000
CIB      DATA FACTBL(20) / Z 4814C9E880A269F7 /                         00008100
      DATA FACTBL(21) / 0.155977626862849790 0D+10 /                    00008200
CIB      DATA FACTBL(21) / Z 485CF84C0CA0E53D /                         00008300
      DATA FACTBL(22) / 0.714779281818586540 0D+10 /                    00008400
CIB      DATA FACTBL(22) / Z 491AA0AA9B22F94E /                         00008500
      DATA FACTBL(23) / 0.335261200823717098 0D+11 /                    00008600
CIB      DATA FACTBL(23) / Z 497CE4FFE925F286 /                         00008700
      DATA FACTBL(24) / 0.160785623545405884 0D+12 /                    00008800
CIB      DATA FACTBL(24) / Z 4A256F91E9F967E8 /                         00008900
      DATA FACTBL(25) / 0.787685471322938293 0D+12 /                    00009000
CIB      DATA FACTBL(25) / Z 4AB765B6745AF034 /                         00009100
      DATA FACTBL(26) / 0.393842735661469141 0D+13 /                    00009200
CIB      DATA FACTBL(26) / Z 4B394FC9045C6B10 /                         00009300
      DATA FACTBL(27) / 0.200821179442459609 0D+14 /                    00009400
CIB      DATA FACTBL(27) / Z 4C1243BB81A3B5F6 /                         00009500
      DATA FACTBL(28) / 0.104349745809073977 0D+15 /                    00009600
CIB      DATA FACTBL(28) / Z 4C5EE7D18806B1FA /                         00009700
      DATA FACTBL(29) / 0.552166953567228500 0D+15 /                    00009800
CIB      DATA FACTBL(29) / Z 4D1F631631FC3FC8 /                         00009900
      DATA FACTBL(30) / 0.297351004601291062 0D+16 /                    00010000
CIB      DATA FACTBL(30) / Z 4DA9064420755EEA /                         00010100
      DATA FACTBL(31) / 0.162865852716949560 0D+17 /                    00010200
CIB      DATA FACTBL(31) / Z 4E39DC9021B9366C /                         00010300
      DATA FACTBL(32) / 0.906798690679354880 0D+17 /                    00010400
CIB      DATA FACTBL(32) / Z 4F14228DC4812BF0 /                         00010500
      DATA FACTBL(33) / 0.512962802680363488 0D+18 /                    00010600
CIB      DATA FACTBL(33) / Z 4F71E68F384C079E /                         00010700
      DATA FACTBL(34) / 0.294674695534107341 0D+19 /                    00010800
CIB      DATA FACTBL(34) / Z 5028E4F2BE4C118C /                         00010900
      DATA FACTBL(35) / 0.171823397428756524 0D+20 /                    00011000
CIB      DATA FACTBL(35) / Z 50EE73EEF5D02E11 /                         00011100
      DATA FACTBL(36) / 0.101652092779175703 0D+21 /                    00011200
CIB      DATA FACTBL(36) / Z 51582B4C7CE89D17 /                         00011300
      DATA FACTBL(37) / 0.609912556675054240 0D+21 /                    00011400
CIB      DATA FACTBL(37) / Z 5221103CAED73AE9 /                         00011500
      DATA FACTBL(38) / 0.370995324650140906 0D+22 /                    00011600
CIB      DATA FACTBL(38) / Z 52C91DF1A584AFB3 /                         00011700
      DATA FACTBL(39) / 0.228696877430935011 0D+23 /                    00011800
CIB      DATA FACTBL(39) / Z 534D7C4AC798B791 /                         00011900
      DATA FACTBL(40) / 0.142821154179615300 0D+24 /                    00012000
CIB      DATA FACTBL(40) / Z 541E3E59AC137462 /                         00012100
      DATA FACTBL(41) / 0.903280290523322411 0D+24 /                    00012200
CIB      DATA FACTBL(41) / Z 54BF46EC33EEDF83 /                         00012300
      DATA FACTBL(42) / 0.578381592144527070 0D+25 /                    00012400
CIB      DATA FACTBL(42) / Z 554C8C54A13C8617 /                         00012500
      DATA FACTBL(43) / 0.374834112342097247 0D+26 /                    00012600
CIB      DATA FACTBL(43) / Z 561F016BCB6BE19A /                         00012700
      DATA FACTBL(44) / 0.245795164849461257 0D+27 /                    00012800
CIB      DATA FACTBL(44) / Z 56CB512695F7B480 /                         00012900
      DATA FACTBL(45) / 0.163042067417843078 0D+28 /                    00013000
CIB      DATA FACTBL(45) / Z 57544A701D9B8D88 /                         00013100
      DATA FACTBL(46) / 0.109371943781520225 0D+29 /                    00013200
CIB      DATA FACTBL(46) / Z 58235708FF48D0A1 /                         00013300
      DATA FACTBL(47) / 0.741796613622095804 0D+29 /                    00013400
CIB      DATA FACTBL(47) / Z 58EFAFFA4C2846FA /                         00013500
      DATA FACTBL(48) / 0.508550136674023691 0D+30 /                    00013600
CIB      DATA FACTBL(48) / Z 5966B36D045FD244 /                         00013700
      DATA FACTBL(49) / 0.352333869966202272 0D+31 /                    00013800
CIB      DATA FACTBL(49) / Z 5A2C788576BCCA76 /                         00013900
      DATA FACTBL(50) / 0.246633708976341607 0D+32 /                    00014000
CIB      DATA FACTBL(50) / Z 5B1374BA63F29894 /                         00014100
      DATA FACTBL(51) / 0.174396368086360610 0D+33 /                    00014200
CIB      DATA FACTBL(51) / Z 5B899310E94A8B18 /                         00014300
      DATA FACTBL(52) / 0.124543918088655869 0D+34 /                    00014400
CIB      DATA FACTBL(52) / Z 5C3D67A7088B14B4 /                         00014500
      DATA FACTBL(53) / 0.898098965431671557 0D+34 /                    00014600
CIB      DATA FACTBL(53) / Z 5D1BACC06AB97112 /                         00014700
      DATA FACTBL(54) / 0.653825915979171439 0D+35 /                    00014800
CIB      DATA FACTBL(54) / Z 5DC979C53884500D /                         00014900
      DATA FACTBL(55) / 0.480461962427038937 0D+36 /                    00015000
CIB      DATA FACTBL(55) / Z 5E5C889A89A4D3ED /                         00015100
      DATA FACTBL(56) / 0.356320127885841961 0D+37 /                    00015200
CIB      DATA FACTBL(56) / Z 5F2AE3F62FB2FD6D /                         00015300
      DATA FACTBL(57) / 0.266645567712059188 0D+38 /                    00015400
CIB      DATA FACTBL(57) / Z 60140F6826450518 /                         00015500
CNV      DATA FACTBL(58) /     0.201312988912482289  0D+39 /            00015600
CIB      DATA FACTBL(58) / Z 60977378AAA8F3F8 /                         00015700
CNV      DATA FACTBL(59) /     0.153315404682076178  0D+40 /            00015800
CIB      DATA FACTBL(59) / Z 614816ABA937E170 /                         00015900
CNV      DATA FACTBL(60) /     0.117763796875648436  0D+41 /            00016000
CIB      DATA FACTBL(60) / Z 62229B903A48F098 /                         00016100
CNV      DATA FACTBL(61) /     0.912194448171078779  0D+41 /            00016200
CIB      DATA FACTBL(61) / Z 6310C11E20742E8E /                         00016300
CNV      DATA FACTBL(62) /     0.712446639319201794  0D+42 /            00016400
CIB      DATA FACTBL(62) / Z 6382DB147A4F0363 /                         00016500
CNV      DATA FACTBL(63) /     0.560981044781264743  0D+43 /            00016600
CIB      DATA FACTBL(63) / Z 644065C05060BC65 /                         00016700
CNV      DATA FACTBL(64) /     0.445264900413724490  0D+44 /            00016800
CIB      DATA FACTBL(64) / Z 651FF239800E9F70 /                         00016900
CNV      DATA FACTBL(65) /     0.356211920330979607  0D+45 /            00017000
CIB      DATA FACTBL(65) / Z 65FF91CC0074FB83 /                         00017100
CNV      DATA FACTBL(66) /     0.287187231472474606  0D+46 /            00017200
CIB      DATA FACTBL(66) / Z 6680C77A34A281DB /                         00017300
CNV      DATA FACTBL(67) /     0.233312009780346082  0D+47 /            00017400
CIB      DATA FACTBL(67) / Z 67416350D5E735C9 /                         00017500
CNV      DATA FACTBL(68) /     0.190974110596668804  0D+48 /            00017600
CIB      DATA FACTBL(68) / Z 682173934666346B /                         00017700
CNV      DATA FACTBL(69) /     0.157481285949690890  0D+49 /            00017800
CIB      DATA FACTBL(69) / Z 69113D90D09DA2E9 /                         00017900
CNV      DATA FACTBL(70) /     0.130813780783272721  0D+50 /            00018000
CIB      DATA FACTBL(70) / Z 698F35D47C81C4CF /                         00018100
CNV      DATA FACTBL(71) /     0.109446661301155697  0D+51 /            00018200
CIB      DATA FACTBL(71) / Z 6A4AE2EE076DAC07 /                         00018300
CNV      DATA FACTBL(72) /     0.922213960297642774  0D+51 /            00018400
CIB      DATA FACTBL(72) / Z 6B277012033F8523 /                         00018500
CNV      DATA FACTBL(73) /     0.782524494037637662  0D+52 /            00018600
CIB      DATA FACTBL(73) / Z 6C14EA3FF6F928CC /                         00018700
CNV      DATA FACTBL(74) /     0.668589220786028248  0D+53 /            00018800
CIB      DATA FACTBL(74) / Z 6CB2B2BB1AB99ED5 /                         00018900
CNV      DATA FACTBL(75) /     0.575142194723999233  0D+54 /            00019000
CIB      DATA FACTBL(75) / Z 6D601384247864EE /                         00019100
CNV      DATA FACTBL(76) /     0.498087751419319663  0D+55 /            00019200
CIB      DATA FACTBL(76) / Z 6E3400B6B1742044 /                         00019300
CNV      DATA FACTBL(77) /     0.434222834690444404  0D+56 /            00019400
CIB      DATA FACTBL(77) / Z 6F1C5598C00D4C38 /                         00019500
CNV      DATA FACTBL(78) /     0.381028991060110636  0D+57 /            00019600
CIB      DATA FACTBL(78) / Z 6FF8A20D740788F4 /                         00019700
CNV      DATA FACTBL(79) /     0.336515693218106815  0D+58 /            00019800
CIB      DATA FACTBL(79) / Z 70893DE0B137AAFD /                         00019900
CNV      DATA FACTBL(80) /     0.299101690580026232  0D+59 /            00020000
CIB      DATA FACTBL(80) / Z 714C3D4BC54C4579 /                         00020100
CNV      DATA FACTBL(81) /     0.267524684928818863  0D+60 /            00020200
CIB      DATA FACTBL(81) / Z 722A9E8008F4CECD /                         00020300
CNV      DATA FACTBL(82) /     0.240772216435936965  0D+61 /            00020400
CIB      DATA FACTBL(82) / Z 7317F9280509B453 /                         00020500
CNV      DATA FACTBL(83) /     0.218028515039038915  0D+62 /            00020600
CIB      DATA FACTBL(83) / Z 73D916506EF03CEF /                         00020700
CNV      DATA FACTBL(84) /     0.198633430462262788  0D+63 /            00020800
CIB      DATA FACTBL(84) / Z 747B9C219210168E /                         00020900
CNV      DATA FACTBL(85) /     0.182050546128413285  0D+64 /            00021000
CIB      DATA FACTBL(85) / Z 7546CE740A0DA582 /                         00021100
CNV      DATA FACTBL(86) /     0.167842310350535579  0D+65 /            00021200
CIB      DATA FACTBL(86) / Z 7628CCDA4507ED5B /                         00021300
CNV      DATA FACTBL(87) /     0.155650555359345683  0D+66 /            00021400
CIB      DATA FACTBL(87) / Z 7717A5D94121476D /                         00021500
CNV      DATA FACTBL(88) /     0.145181172966040185  0D+67 /            00021600
CIB      DATA FACTBL(88) / Z 77DC928A84D09191 /                         00021700
CNV      DATA FACTBL(89) /     0.136192012341913224  0D+68 /            00021800
CIB      DATA FACTBL(89) / Z 78815271BC3A3F73 /                         00021900
CNV      DATA FACTBL(90) /     0.128483287477042949  0D+69 /            00022000
CIB      DATA FACTBL(90) / Z 794C405921789CB9 /                         00022100
CNV      DATA FACTBL(91) /     0.121889948908093383  0D+70 /            00022200
CIB      DATA FACTBL(91) / Z 7A2D36248645D4DA /                         00022300
CNV      DATA FACTBL(92) /     0.116275600522138910  0D+71 /            00022400
CIB      DATA FACTBL(92) / Z 7B1AF4A48C0EAE0E /                         00022500
CNV      DATA FACTBL(93) /     0.111527638075238134  0D+72 /            00022600
CIB      DATA FACTBL(93) / Z 7C1028CA3E2C98F5 /                         00022700
CNV      DATA FACTBL(94) /     0.107553359179601711  0D+73 /            00022800
CIB      DATA FACTBL(94) / Z 7C9BD5C2156FF2C5 /                         00022900
CNV      DATA FACTBL(95) /     0.104276850578483769  0D+74 /            00023000
CIB      DATA FACTBL(95) / Z 7D5E6E055E5F37AB /                         00023100
CNV      DATA FACTBL(96) /     0.101636501751285494  0D+75 /            00023200
CIB      DATA FACTBL(96) / Z 7E3986335FD24508 /                         00023300
CNV      DATA FACTBL(97) /     0.995830274128553335  0D+75 /            00023400
CIB      DATA FACTBL(97) / Z 7F2339EB5C74427A /                         00023500
C                                                                       00023600
C                                                                       00023700
C                                                                       00023800
      DATA I /0/                                                        00023900
C                                                                       00024100
      END                                                               00024200
      FUNCTION GAMMA ( XIN )                                            00000100
C                                                                       00000200
C     THIS IS THE BKY SUBROUTINE GAMLOG                                 00000300
C                                                                       00000400
C     MODIFIED TO GAMMA() AND ALGAMA() FORM                             00000500
C                                                                       00000600
C     12/28/79 - PUT IN ALGAMA ARGUMENT - RPG                           00000610
      DATA YMX,TS1,TS11,TS2,TS22,TS3,PI,CST/300.,10.,4.,11.,2.,50.,     00000800
     1 3.14159265,0.434294482/,A0,A1,A2,A3,A4,A5,A6,A7,A8,A9,           00000900
     2 A10/.999999999,.422784337,.411840252,                            00001000
     3 0.8157821878492E-01,0.7423790760629E-01,-0.21090746731E-03,      00001100
     4 0.1097369584174E-01,-0.246674798054E-02,0.153976810472E-02,      00001200
     5-0.34423420456E-03,0.6771057117E-04/,B1,B2,B3,B4/                 00001300
     6 8.3333333333333E-02,-2.7777777777777E-03,7.9365079365079E-04,    00001400
     7 -5.952380952380E-04/                                             00001500
      NIN = 0                                                           00001600
 999  CONTINUE                                                          00001700
           PIL = ALOG(PI)                                               00001800
           FLGTW = 0.5*(PIL + .693147180559945)                         00001900
      X = XIN                                                           00002000
      N = NIN                                                           00002100
      N1 = N - 1                                                        00002200
      H = 1.0                                                           00002300
      IF (X-TS1) 51,1,1                                                 00002400
   51 IF (X) 52,9,2                                                     00002500
   52 Z=AINT(X)                                                         00002600
      IF (X - Z) 9,9,53                                                 00002700
   53 IF (X + TS2) 3,9,4                                                00002800
    2 CONTINUE                                                          00002900
      IF (X - TS11) 54,31,31                                            00003000
   54 IX = X + 1.0                                                      00003100
      GO TO (32,33,34,35),IX                                            00003200
   32 H = 1.0 / (X*(X + 1.0))                                           00003300
      Y = X                                                             00003400
      GO TO 40                                                          00003500
   33 H = 1.0 / X                                                       00003600
      Y = X - 1.0                                                       00003700
      GO TO 40                                                          00003800
   34 H = 1.0                                                           00003900
      Y = X - 2.0                                                       00004000
      GO TO 40                                                          00004100
   35 H = X - 1.0                                                       00004200
      Y = X - 3.0                                                       00004300
      GO TO 40                                                          00004400
   31 Y = X                                                             00004500
    6 H = H / Y                                                         00004600
      Y = Y + 1.0                                                       00004700
      IF (Y - TS1) 6,10,10                                              00004800
    4 CONTINUE                                                          00004900
      IF (X + TS22) 36,56,56                                            00005000
   56 IF (X + 1.0) 57,57,37                                             00005100
   57 H = 1.0 / ((((X + 6.0)*X + 11.0)*X + 6.0)*X)                      00005200
      Y = X + 2.0                                                       00005300
      GO TO 40                                                          00005400
   37 H = 1.0 / (((X + 3.0)*X + 2.0)*X)                                 00005500
      Y = X + 1.0                                                       00005600
      GO TO 40                                                          00005700
   36 Y = X                                                             00005800
    5 Y = Y - 1.0                                                       00005900
      H = H * ABS(Y)                                                    00006000
      IF (Y + TS2) 58,9,5                                               00006100
   58 Y = 1.0 - Y                                                       00006200
      GO TO 10                                                          00006300
    1 Y = X                                                             00006400
      GO TO 10                                                          00006500
    3 Y = 1.0 - X                                                       00006600
      GO TO 10                                                          00006700
   40 CONTINUE                                                          00006800
      S = H * (A0 + Y*(A1 + Y*(A2 + Y*(A3 + Y*(A4 + Y*(A5 + Y*(A6 + Y*( 00006900
     1A7 + Y*(A8 + Y*(A9 + Y*A10))))))))))                              00007000
      IF (N1) 59,41,42                                                  00007100
   59 GAMMA = S                                                         00007200
      RETURN                                                            00007300
   41 GAMMA = ALOG(ABS(S))                                              00007400
      RETURN                                                            00007500
   42 GAMMA = ALOG10(ABS(S))                                            00007600
      RETURN                                                            00007700
   10 S = (Y - 0.5)*ALOG(Y) - Y + FLGTW                                 00007800
      YY = 1.0 / (Y*Y)                                                  00007900
      IF (N) 60,60,25                                                   00008000
   60 IF (Y-YMX) 25,25,15                                               00008100
   25 CONTINUE                                                          00008200
      IF (Y - TS3) 61,61,13                                             00008300
   61 S = S + (1.0/Y)*(B1 + YY*(B2 + YY*(B3 + YY*B4)))                  00008400
      GO TO 14                                                          00008500
   13 CONTINUE                                                          00008600
      S = S + (1.0/Y)*(B1 + YY*B2)                                      00008700
   14 CONTINUE                                                          00008800
      IF (X) 11,62,62                                                   00008900
   62 IF (N1) 63,7,8                                                    00009000
   63 GAMMA = EXP(S)*H                                                  00009100
      RETURN                                                            00009200
    7 IF (H - 1.0) 64,26,64                                             00009300
   64 GAMMA = ALOG(ABS(H*EXP(S)))                                       00009400
      RETURN                                                            00009500
   26 GAMMA = S                                                         00009600
      RETURN                                                            00009700
    8 IF (H-1.) 67,27,67                                                00009800
   67 GAMMA = ALOG 10 (ABS(H*EXP(S)))                                   00009900
      RETURN                                                            00010000
   27 GAMMA = CST*S                                                     00010100
      RETURN                                                            00010200
   11 SN = 1.0                                                          00010300
      Z=AINT(-X)                                                        00010400
      IF (AMOD(Z,2.)) 20,65,20                                          00010500
   65 SN = -1.0                                                         00010600
   20 CONTINUE                                                          00010700
      IF (N1) 66,19,21                                                  00010800
   66 GAMMA = (SN*PI/(EXP(S)*(ABS(SIN(PI*(X+Z))))))*H                   00010900
      RETURN                                                            00011000
   19 IF (H - 1.0) 68,28,68                                             00011100
   68 GAMMA = ALOG(ABS(SN*H*PI/(EXP(S)*SIN(PI*(X+Z)))))                 00011200
      RETURN                                                            00011300
   28 GAMMA = PIL - S - ALOG(ABS(SIN(PI*(X+Z))))                        00011400
      RETURN                                                            00011500
   21 IF (H - 1.0) 69,29,69                                             00011600
   69 GAMMA = ALOG10(ABS(SN*H*PI/(EXP(S)*SIN(PI*(X+Z)))))               00011700
      RETURN                                                            00011800
   29 GAMMA = CST*(PIL-S)-ALOG10(ABS(SIN(PI*(X+Z))))                    00011900
      RETURN                                                            00012000
    9 PRINT 12,X                                                        00012100
   12 FORMAT(  // 69HGAMMA FUNCTION CALLED WITH ZERO OR NEGATIVE INTEGER00012200
     1 ARGUMENT --- X = ,1PE21.10)                                      00012300
   17 CONTINUE                                                          00012400
   18 GAMMA = 0.0                                                       00012500
      RETURN                                                            00012600
   15 IF (X) 18,18,70                                                   00012700
   70 PRINT 16,X                                                        00012800
   16 FORMAT(  // 55HGAMMA FUNCTION CALLED WITH ARGUMENT TOO LARGE ---  00012900
     1X = ,1PE21.10)                                                    00013000
      GO TO 17                                                          00013100
C                                                                       00013200
      ENTRY ALGAMA ( XIN )                                              00013300
      NIN = 1                                                           00013400
      GO TO 999                                                         00013500
      END                                                               00013600
      SUBROUTINE GAUSSL(N,X,W)                                          00000100
C     POINTS AND WEIGHTS FOR GAUSS-LEGENDRE INTEGRATION OVER (-1 1)     00000200
C     COMPUTED BY METHOD GIVEN ON PP.88 AND 89 OF                       00000300
C     'METHODS OF NUMERICAL INTEGRATION' DAVIS AND RABINOWITZ           00000400
C     ACADEMIC PRESS (1975)                                             00000500
C                                                                       00000600
C     THIS SUBROUTINE IS AN ADAPTATION OF SUBROUTINE GRULE LISTED ON    00000700
C     P.369 (APP 2) OF DAVIS AND RABINOWITZ                             00000800
C     MHM  JUNE 29, 1975                                                00000900
C                                                                       00001000
C     METHOD: - - - - - - - - - - -                                     00001100
C     BASED ON A VERY CLOSE ESTIMATE FOR ZEROS OF PN(X)                 00001200
C     THIS ESTIMATE USES THE ASYMPTOTIC RELATION BETWEEN PN(X)          00001300
C     AND THE COSINE FUNCTION (WHICH IS REMARKABLY ACCURATE)            00001400
C     THE INITIAL ESTIMATE IS IMPROVED BY A FOURTH-ORDER                00001500
C     NEWTON-RAPHSON PROCEDURE                                          00001600
C     TO SIMPLIFY MATTERS THIS IS CARRIED OUT IN TWO STEPS              00001700
C     1. A THIRD-ORDER NEWTON-RAPHSON FROM X0 TO X0+H                   00001800
C     2. A FIRST-ORDER NEWTON-RAPHSON FROM X0+H TO X0+H+H1              00001900
C     FOR ALL ZEROS AND N UP TO ABOUT 2000,THIS                         00002000
C     SUFFICES TO GIVE ZEROS TO AT LEAST 13 SIGNIFICANT FIGURES         00002100
C     GAUSS POINTS AND WEIGHTS ARE ALSO GOOD TO 13 OR MORE FIGURES      00002200
C                                                                       00002300
C     MORE DETAILS (WITH FORMULAS) ARE GIVEN BELOW CLOSE TO THE         00002400
C     APPROPRIATE PARTS OF THE PROGRAM LISTING                          00002500
C                                                                       00002600
C     ACCURACY AND TIMING - - - - - -                                   00002700
C                                                                       00002800
C     GAUSSL EXECUTION TIMES : IN MILLISECONDS                          00002900
C                            : UNCERTAINTY ABOUT .1 MS                  00003000
C     /195   T = .035  +  .0133 * N  +  .00103 * N**2                   00003100
C     /75    T = 12.1*T(195)                                            00003200
C                                                                       00003300
C     NUMBER OF SIGNIFICANT DIGITS : (TO .5 DIGIT)                      00003400
C     DIGITS(BASE 10)=16.1-.91*LOG10(N)                                 00003500
C                                                                       00003600
C     NUMBER OF BINARY DIGITS LOST : (TO 1.5 BINARY DIGITS)             00003700
C     LOST DIGITS(BASE 2) = 2.5 + 3.2 * LOG10(N)                        00003800
C                                                                       00003900
C     SUMMARY OF TIMES AND ACCURACIES FOR VARIOUS VALUES OF N           00004000
C                                                                       00004100
C       N   DIGITS(10)  DIGITS(2)    TIME(MS)                           00004200
C           ACCURACY      LOST      /195    /75                         00004300
C     ------------------------------------------------                  00004400
C       10    15.7        3.7        .28     3.5                        00004500
C       20    15.2        5.5        .75     9.2                        00004600
C       50    14.4        8.1       3.35    40.8                        00004700
C      100    14.0        9.5      11.7    144.0                        00004800
C      200    13.8       10.1      44.5    540.                         00004900
C      500    13.7       10.5     270.0                                 00005000
C     1000    13.5       11.1    1060.0   12700.0                       00005100
C     2000    13.2       12.1    4110.0                                 00005200
C     ______________________________________________                    00005300
C                                                                       00005400
C     ABSCISSAE COME OUT IN ASCENDING ORDER FROM NEAR -1 TO NEAR +1     00005500
C -  -  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     00005600
      IMPLICIT REAL*8 ( A-H, O-Z )                                      00005700
CDC      ALLOCLEVEL X,W                                                 00005800
      DIMENSION X(N),W(N)                                               00005900
C                                                                       00006000
C     TAKE CARE OF SPECIAL CASES (N.LE.1)                               00006100
C                                                                       00006200
      IF(N.GT.0) GO TO 200                                              00006300
      PRINT 100,N                                                       00006400
  100 FORMAT(1H0,2X '0000 NUMBER OF GAUSS POINTS = ',I5,' 0000'/3X, '00000006500
     *0             ABSURD             0000')                           00006500
      GO TO 600                                                         00006700
  200 IF(N.NE.1) GO TO 300                                              00006800
      X(1)=0                                                            00006900
      W(1)=2                                                            00007000
      GO TO 600                                                         00007100
C                                                                       00007200
C     BEGIN COMPUTATION                                                 00007300
C     INITIAL ESTIMATE X0(I,N) OF I'TH ZERO OF PN(X) IS                 00007400
C     -  -  -                                                           00007500
C     X0(I,N)=(1-.125*(N**-2)+.125*(N**-3))*COS(T)  + O(N**-4)          00007600
C     T=((4*I-1)/(4*N+2))*PI                                            00007700
C     -  -  -                                                           00007800
C     THE ZEROS ARE ORDERED SO THAT I=1 IS NEAREST TO X=1               00007900
C     1 > X(1) > X(2) > X(3) > ............> X(M) GE 0                  00008000
C                                                                       00008100
  300 M=(N+1)/2                                                         00008200
      E1=N*(N+1)                                                        00008300
      E2=3.1415 92653 58979 320D0/(4*N+2)                               00008400
      E3=1.-(1.-1.0D0/N)/(8.0D0*N*N)                                    00008500
      DO 500 I=1,M                                                      00008600
      T=(4*I-1)*E2                                                      00008700
      X0=E3*DCOS(T)                                                     00008800
C                                                                       00008900
C     USE RECURSION FORMULAS AND THE 2ND ORDER DIFFERENTIAL EQUATION    00009000
C     SATISFIED BY PN TO COMPUTE PK AND ITS FIRST FOUR DERIVATIVES      00009100
C     AT X0                                                             00009200
C                                                                       00009300
      PKM1=1                                                            00009400
      PK=X0                                                             00009500
      FK=1                                                              00009600
      DO 400 K=2,N                                                      00009700
      FK=FK+1                                                           00009800
      T1=X0*PK                                                          00009900
      PKP1=T1-PKM1-(T1-PKM1)/FK+T1                                      00010000
      PKM1=PK                                                           00010100
  400 PK=PKP1                                                           00010200
      DEN=1.-X0*X0                                                      00010300
      D1=N*(PKM1-X0*PK)                                                 00010400
      DPN=D1/DEN                                                        00010500
      D2PN=(2.*X0*DPN-E1*PK)/DEN                                        00010600
      D3PN=(4.*X0*D2PN+(2.-E1)*DPN)/DEN                                 00010700
      D4PN=(6.*X0*D3PN+(6.-E1)*D2PN)/DEN                                00010800
C                                                                       00010900
C     USE THIRD-ORDER NEWTON-RAPHSON METHOD TO IMPROVE ON THE INITIAL   00011000
C     ESTIMATE X0 OF THE ZERO OF PN                                     00011100
C     I.E H IS DETERMINED SUCH THAT PN(X0+H)=0 UP TO AND INCLUDING      00011200
C     THE TERM IN D3(PN)/DX IN THE TAYLOR EXPANSION OF PN ABOUT X0      00011300
C                                                                       00011400
      U=PK/DPN                                                          00011500
      V=D2PN/DPN                                                        00011600
      H=-U*(1.+.5*U*(V+U*(V*V-U*D3PN/(3.*DPN))))                        00011700
C                                                                       00011800
C     EVALUE PN AT AND D(PN)/DX AT X0+H USING TAYLOR EXPANSION OF PN    00011900
C     AT XO UP TO TERMS OF ORDER H**3. FOR D(PN)/DX THIS REQUIRES       00012000
C     D4(PN)/DX  AT X0.                                                 00012100
C     THESE VALUES OF PN AND DPN AT X0+H ARE THEN USED TO               00012200
C     COMPUTE A FURTHER FIRST-ORDER NEWTON-RAPHSON CORRECTION           00012300
C     I.E. FIND H1 SUCH THAT PN(X0+H+H1)=0 UP TO THE FIRST-DERIVATIVE   00012400
C     TERM (TERM OF ORDER H1) IN THE TAYLOR EXPANSION OF PN ABOUT X0+H  00012500
C                                                                       00012600
      P=PK+H*(DPN+.5*H*(D2PN+H/3.*(D3PN+.25*H*D4PN)))                   00012700
      DP=DPN+H*(D2PN+.5*H*(D3PN+H*D4PN/3.))                             00012800
      H=H-P/DP                                                          00012900
      X(N+1-I)=X0+H                                                     00013000
      X(I)=-X(N+1-I)                                                    00013100
C                                                                       00013200
C     CALCULATE WEIGHT W(I) CORRESPONDING TO I'TH GAUSS POINT X0I       00013300
C     W(I)=2*(1-X0I**2)/(N*PN-1(X0I))**2                                00013400
C     PN-1(X) IS THE LEGENDRE POLYNOMIAL PN OF ORDER N-1                00013500
C                                                                       00013600
      FX=D1-H*E1*(PK+.5*H*(DPN+H/3.*(D2PN+.25*H*(D3PN+.2*H*D4PN))))     00013700
      W(N+1-I)=2.*(1-X(I)*X(I))/(FX*FX)                                 00013800
  500 W(I)=W(N+1-I)                                                     00013900
      IF(M+M.GT.N) X(M)=0                                               00014000
  600 RETURN                                                            00014100
      END                                                               00014200
      SUBROUTINE INTRPC (NCUBIC, XCUBES, AS, BS, CS, DS, NPTS, XS, YS)  00000100
C                                                                       00000200
C     INTERPOLATES TO A GRID OF XS USING INPUT CUBICS.                  00000300
C                                                                       00000400
C       THIS ROUTINE MAY BE USED TO FIND INTERPOLATED VALUES OF A       00000500
C     FUNCTION AFTER EITHER THE B88888 ROUTINE "SPLNCB" OR THE          00000600
C     AMD ROUTINE "SMOOTH" HAS FOUND THE COEFICENTS OF THE INTERPOLATING00000700
C     CUBICS.  FOR EACH DESIRE POINT, THE APPROPRIATE CUBIC IS FOUND    00000800
C     AND IS USED TO DO THE INTERPOLATIONS.  IF NECESSARY THE END       00000900
C     CUBICS WILL BE USED FOR EXTRAPOLATION, BUT THIS IS NOT RECOMMENDED00001000
C     NOTE THAT THE X'S ATWHICH THE OUTPUT IS DESIRED NEED NOT BE       00001100
C     ORDERED.                                                          00001200
C                                                                       00001300
C                                                                       00001400
C     THE ARGUMENTS ARE:                                                00001500
C                                                                       00001600
C     NCUBIC - THE NUMBER OF POINTS USED IN DETERMINING THE CUBICS.     00001700
C     XCUBES - THE X VALUES USED IN FINDING THE CUBICS.  THIS ARRAY     00001800
C              HAS NCUBIC ELEMENTS.  IT MUST BE MONOTONIC (INCREASING   00001900
C              OR DCREASING).                                           00002000
C     AS, BS, CS, AND DS - THE COEFICENTS OF THE CUBICS.  EACH ARRAY    00002100
C              IS NCUBIC-1 ELEMENTS LONG.  THE COEFICENTS ARE           00002200
C              DEFINED AS:                                              00002300
C        YS(N) = F(XS(N)) = AS(I) + DEL*(BS(I) + DEL*(CS(I) +           00002400
C                 DEL*DS(I)))                                           00002500
C        WHERE                                                          00002600
C          DEL = XS(N)-XCUBES(I);  XCUBES(I) < XS(N) < XCUBES(I+1)      00002700
C        THE AS, ..., DS MAY BE FOUND USING EITHER SMOOTH OR SPLNCB.    00002800
C     NPTS - THE NUMBER OF OUTPUT VALUES OF YS TO BE GENERATED.         00002900
C     XS  - THE ARRAY OF X-VALUES AT WHICH THE INTERPOLATED VALUES      00003000
C           ARE DESIRED.  THIS ARRAY HAS NPTS ELEMENTS.  IF ANY         00003100
C           XS(I) IS NOT IN THE RANGE DEFINED BY XCUBES,                00003200
C           THE EXTREME CUBICS WILL BE USED FOR EXTRAPOLATION.          00003300
C     YS  - THIS ARRAY WILL GET F(XS(I)).  IT WILL HAVE NPTS ELEMENTS.  00003400
C                                                                       00003500
C     S. PIEPER                                                         00003600
C                                                                       00003700
C     PDP  2/21/72                                                      00003800
C     12/10/72  /360 VERSION                                            00003900
C     6/26/73  REMOVE PRINTING OF ERRORS                                00004000
C     7/15/73  GENERAL XS                                               00004100
C     8/29/74 - FIX ERRORS FOR NON-INCREASING XS.                       00004200
C     9/4/74 - CHANGE NAME TO INTRPC FROM INTRP1 - S. P.                00004300
C     12/28/79 - INF FOR VAX, UNIVAC - RPG                              00004310
C                                                                       00004400
      IMPLICIT REAL*8 ( A-H, O-Z )                                      00004500
C                                                                       00004600
CDC      ALLOCLEVEL XCUBES,AS,BS,CS,DS,XS,YS                            00004700
      DIMENSION XCUBES(1), AS(1), BS(1), CS(1), DS(1), XS(1), YS(1)     00004800
CIB      REAL*8  INF / Z7FFFFFFF FFFFFFFF /                             00004900
CDC      REAL INF /1.E300/                                              00005000
      REAL*8 INF                                                        00005010
        data INF / 1.D38 /                                              00005010
CUN      REAL*8 INF / 1.D300 /                                          00005020
C                                                                       00005100
      XSIGN = DSIGN (1.D0, XCUBES(2)-XCUBES(1))                         00005300
C                                                                       00005400
C     THESE ARE FAKE STARTING CONDITIONS THAT WILL CAUSE DIRECT         00005500
C     BRANCH TO STMNT 200 INSIDE THE LOOP.                              00005600
      XPREV = INF                                                       00005700
      XNEXT = INF                                                       00005800
C                                                                       00005900
C                                                                       00006000
      DO 399 I = 1, NPTS                                                00006100
      X = XS(I)                                                         00006200
      XC = XSIGN*X                                                      00006300
C                                                                       00006400
C     IS X STILL INSIDE THIS SPLINE                                     00006500
C                                                                       00006600
 100  IF (XC .LT. XNEXT) GO TO 200                                      00006700
C                                                                       00006800
C     GO TO THE NEXT SPLINE.                                            00006900
C                                                                       00007000
      XPREV = XNEXT                                                     00007100
      XBASE = XNEXT*XSIGN                                               00007200
      N = N+1                                                           00007300
      XNEXT = XCUBES(N+1)*XSIGN                                         00007400
C     WHEN WE ARE AT THE LAST SPLINE, ALLOW IT TO GO TO INFINITY        00007500
      IF (N .EQ. NCUBIC-1) XNEXT = INF                                  00007600
C                                                                       00007700
 150  A = AS(N)                                                         00007800
      B = BS(N)                                                         00007900
      C = CS(N)                                                         00008000
      D = DS(N)                                                         00008100
C     SEE IF X IS INSIDE THIS NEW SPLINE                                00008200
      GO TO 100                                                         00008300
C                                                                       00008400
C     WE HAVE NOT PASSED THE END OF THIS SPLINE, ARE WE STILL           00008500
C     BEYOND ITS BEGINNING.                                             00008600
C                                                                       00008700
 200  IF (XC .GE. XPREV) GO TO 300                                      00008800
C                                                                       00008900
C     THE XS ARRAY IS NOT MONOTONIC INCREASING -- GO BACK TO START      00009000
C                                                                       00009100
C     ALSO INITIALIZE FOR FIRST TIME                                    00009200
      N = 1                                                             00009300
      XPREV = -INF                                                      00009400
      XNEXT = XCUBES(2)*XSIGN                                           00009500
      XBASE = XCUBES(1)                                                 00009600
      GO TO 150                                                         00009700
C                                                                       00009800
C     THE SET OF SPLINES TO USE HAS BEEN FOUND, USE IT                  00009900
C                                                                       00010000
 300  DEL = X-XBASE                                                     00010100
      YS(I) = A + DEL*(B + DEL*(C + DEL*D))                             00010200
 399  CONTINUE                                                          00010300
C                                                                       00010400
      RETURN                                                            00010500
      END                                                               00010600
      SUBROUTINE LAGUER ( NN, X, W, ALF, EPS, CSX, CSW, TSX, TSW,       00000020
     1   A, B, C )                                                      00000030
C                                                                       00000040
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 00000050
C                                                                       00000060
C     POINTS AND WEIGHTS FOR GAUSS-LAGUERRE INTEGRATION                 00000070
C                                                                       00000080
C     THE ARGUMENTS ARE:                                                00000090
C                                                                       00000100
C     NN - THE NUMBER OF POINTS DESIRED.                                00000110
C     X - WILL BE SET TO THE POINTS WITH THE MINIMUM POINT FIRST.       00000120
C     W - WILL BE SET TO THE WEIGHTS.                                   00000130
C         X AND W MUST BE REAL*8 ARRAYS OF DIMENSION NN.                00000140
C     ALF - THE POWER OF X TO BE INCLUDED IN THE WEIGHT FUNCTION.       00000150
C           MUST HAVE -1 < ALF.  LARGE ALF CAN RESULT IN OVERFLOWS.     00000160
C           NOTE THAT ALF DOES NOT NEED TO BE AN INTEGER.               00000170
C     EPS - THE DESIRED ACCURACY OF THE ROOTS.  EPS = 1.D-20 IS         00000180
C           REASONABLE.                                                 00000190
C     CSX, CSW, TSX, TSW - WILL BE SET TO AN INDICATION OF THE          00000200
C          ACCURACY OF X AND W AS FOLLOWS:                              00000210
C             CSX = CALC SUM X(I)        TSX = TRUE SUM X(I)            00000220
C             CSW = CALC SUM W(I)        TSW = TRUE SUM W(I)            00000230
C     A, B, C - REAL*8 SCRATCH ARRAYS OF DIMENSION NN.                  00000240
C                                                                       00000250
C     THE POINTS AND WEIGHTS THAT ARE RETURNED ARE USED AS:             00000260
C                                                                       00000270
C     INTEGRAL( 0 TO INFINITY ) DX EXP(-X) X**ALF F(X)                  00000280
C              =  SUM  W(I) F( X(I) )                                   00000290
C                                                                       00000300
C                   THIS PROGRAM IS TAKEN FROM:                         00000310
C                   STROUD AND SECREST, GAUSSIAN                        00000320
C                   QUADRATURE FORMULAS,PRENTICE                        00000330
C                   HALL (1966)  ..  PAGES 32-33                        00000340
C                                                                       00000350
C                                                                       00000360
C     THIS ROUTINE CALLS  LGROOT, LGRECR , LAGBC  AND  DGAMMA           00000370
C                                                                       00000380
C     4/22/75 - F. SERDUKE                                              00000390
C     1/27/80 - USE STANDARD (A&S) LAGUERRE POLYNOMIALS TO PREVENT      00000400
C        OVERFLOW; PASS A, B, C AS ARGUMENTS - RPG                      00000410
C                                                                       00000420
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 00000430
C                                                                       00000440
      IMPLICIT REAL*8 ( A-H, O-Z )                                      00000460
C                                                                       00000470
      DIMENSION X(NN) , W(NN) , A(NN) , B(NN) , C(NN)                   00000480
CFR%F                                                                   00000490
C                                                                       00000500
C     THE COEFFICIENTS                                                  00000510
C                                                                       00000520
C                     A(N)  =  1 / N                                    00000530
C                                                                       00000540
C                     B(N)  =  (ALF + 2N - 1)                           00000550
C                                                                       00000560
C                     C(N)  =  (ALF + N - 1) / N                        00000570
C                                                                       00000580
C      IN THE RECURRENCE RELATION                                       00000590
C                                                                       00000600
C             L(N)  =  A(N) * (B(N) - X) * L(N-1)  -  C(N)*L(N-2)       00000610
C                                                                       00000620
C     THESE ARE THE LAGUERRE POLYNOMIALS DEFINED IN                     00000630
C     ABRAMOWITZ AND STEGUN, HANDBOOK OF MATHEMATICAL FUNCTIONS.        00000640
C                                                                       00000650
      CALL LAGBC ( NN , ALF , A , B , C )                               00000660
      FN = NN                                                           00000670
      CSX = 0.D0                                                        00000680
      CSW = 0.D0                                                        00000690
      CC = DGAMMA(ALF+1.D0)                                             00000700
      TSX = FN*(FN+ALF)                                                 00000710
      TSW = CC                                                          00000720
      DO 1 J = 2,NN                                                     00000730
    1 CC = CC*C(J)                                                      00000740
      DO 7 I = 1,NN                                                     00000750
      IF ( I - 2 ) 2, 4, 5                                              00000760
C         SMALLEST ZERO                                                 00000770
    2 XT = (1.D0+ALF)*(3.D0+0.920D0*ALF) / (1.D0+2.40D0*FN + 1.80D0*ALF)00000780
      GO TO 6                                                           00000790
C        SECOND ZERO                                                    00000800
    4 XT = XT + (15.D0 + 6.250D0*ALF) / (1.D0 + 0.90D0*ALF + 2.50D0*FN) 00000810
      GO TO 6                                                           00000820
C        ALL OTHER ZEROS                                                00000830
    5 FI = I-2                                                          00000840
      R1 = (1.D0+2.550D0*FI) / (1.90D0*FI)                              00000850
      R2 = 1.260D0 * FI*ALF / (1.D0+3.50D0*FI)                          00000860
      RATIO = (R1+R2) / (1.D0 + 0.30D0*ALF)                             00000870
      XT = XT + RATIO*(XT-X(I-2))                                       00000880
C                                                                       00000890
    6 CONTINUE                                                          00000900
      CALL LGROOT(XT,NN,ALF,DPN,PN1,A,B,C,EPS)                          00000910
      X(I) = XT                                                         00000920
      W(I) = -CC / DPN / PN1                                            00000930
      CSX = CSX + XT                                                    00000940
    7 CSW = CSW + W(I)                                                  00000950
      RETURN                                                            00000960
      END                                                               00000970
      SUBROUTINE LGROOT ( X , NN , ALF , DPN , PN1 , A , B , C , EPS )  00000980
      IMPLICIT REAL*8 ( A-H, O-Z )                                      00001000
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 00001010
C        THIS SUBROUITNE IMPROVES THE APPROXIMATE ROOT S OF THE         00001020
C        NN-TH ORDER LAGUERRE POLYNOMIAL.  IN ADDITION, IT ALSO         00001030
C        RETURNS                                                        00001040
C                DPN  =  DERIVATIVE OF P(N) AT X                        00001050
C                PN1  =  VALUE OF P(N-1) AT X                           00001060
C        WHICH ARE NEEDED IN THE CALCULATION OF THE WEIGHTS FOR         00001070

C        NN-TH ORDER LAGUERRE QUADRATURE FORMULA.                       00001080
C                                                                       00001090
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 00001100
C                                                                       00001110
      DIMENSION A(NN), B(NN) , C(NN)                                    00001120
CFR%F                                                                   00001130
C                                                                       00001140
      DO 2 ITER = 1, 10                                                 00001150
      CALL LGRECR ( P , DP , PN1 , X , NN , ALF , A , B , C )           00001160
      D = P/DP                                                          00001170
      X = X - D                                                         00001180
      IF ( DABS(D/X) .LT. EPS ) GO TO 3                                 00001190
    2 CONTINUE                                                          00001200
    3 DPN = DP                                                          00001210
      RETURN                                                            00001220
      END                                                               00001230
      SUBROUTINE LGRECR ( PN , DPN , PN1 , X , NN , ALF , A , B , C )   00001250
      IMPLICIT REAL*8 ( A-H, O-Z )                                      00001270
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 00001280
C        THIS SUBROUTINE IS THE RECURRENCE RELATION FOR THE LAGUERRE    00001290
C        POLYNOMIALS                                                    00001300
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 00001310
      DIMENSION A(NN) , B(NN) , C(NN)                                   00001320
C                                                                       00001330
CFR%F                                                                   00001350
C                                                                       00001360
      P1 = 1.D0                                                         00001370
      P = ALF + 1.D0 - X                                                00001380
      DP1 = 0.D0                                                        00001390
      DP = -1.D0                                                        00001400
      DO 1 J = 2,NN                                                     00001410
      Q = A(J) * (B(J) - X) * P - C(J) * P1                             00001420
      DQ = A(J) * (B(J) - X) * DP - C(J) * DP1 - A(J)*P                 00001430
      P1 = P                                                            00001440
      P = Q                                                             00001450
      DP1 = DP                                                          00001460
    1 DP = DQ                                                           00001470
      PN = P                                                            00001480
      DPN = DP                                                          00001490
      PN1 = P1                                                          00001500
      RETURN                                                            00001510
      END                                                               00001520
      SUBROUTINE LAGBC ( NN , ALF , A , B , C )                         00001540
      IMPLICIT REAL*8 ( A-H, O-Z )                                      00001560
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 00001570
C     THIS SUBROUTINE CALCULATES THE COEFFICIENTS A(N), B(N), AND C(N)  00001580
C     REQUIRED BY THE LAGUERRE QUADRATURE W AND X GENERATION            00001590
C     SUBROUTINE.                                                       00001600
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 00001610
      DIMENSION A(NN) , B(NN) , C(NN)                                   00001620
CFR%F                                                                   00001630
C                                                                       00001640
      DO 10 N = 1,NN                                                    00001650
      EN = N                                                            00001660
      A(N) = 1.D0 / EN                                                  00001670
      B(N) = ALF + EN + EN - 1.D0                                       00001680
      C(N) = (ALF + EN - 1.D0) / EN                                     00001690
   10 CONTINUE                                                          00001700
      RETURN                                                            00001710
      END                                                               00001720
      SUBROUTINE LINSAV (ASV, NMAX, NDIM, BSV, MMAX, A, B, POWER, DET)  00000010
C                                                                       00000020
C     UNIVAC VERSION OF 7/21/70                                         00000030
C     7/14/72  ARGONNE VERSION                                          00000040
C     11/06/72  REVERSE ORDER OF SUBSCRIPTS FOR /195                    00000050
C     1/16/79 - CORRECT FORMAT ERROR - RPG                              00000060
C                                                                       00000070
C     SOLUTION OF SEVERAL SETS OF SIMULTANEOUS LINEAR EQUATIONS.        00000080
C     THE MATRIX MUST BE SYMETRIC.  ITS MODE MAY BE CHANGED BY CHANGING 00000090
C      THE MARKED CARDS BELOW.  NOTE THAT EVEN COMPLEX MATRICES ARE     00000100
C      SYMETRIC... NOT HERMEATIAN.                                      00000110
C     SOLVE AX = B                                                      00000120
C     THERE IS NO PIVOT SEARCHING.  THE DETERMINANT IS COMPUTED BUT     00000130
C     NO CHECK FOR SINGULARITY IS MADE.                                 00000140
C                                                                       00000150
C     'LINSAV'... THIS ENTRY PRESERVES THE ORIGINAL MATRIX AND          00000160
C      INHOMOGENEOUS TERM.                                              00000170
C     'ASV' IS THE MATRIX                                               00000180
C     'BSV' IS THE INHOMOGENEOUS TERM.                                  00000190
C     'A' WILL CONTAIN THE TRIANGULAR MATRIX.  THE UPPER HALF IS THE    00000200
C      TI  TRIANGLE, THE LOWER HALF IS JUNK.                            00000210
C     'B' WILL HAVE THE SOLUTIN FOR X.                                  00000220
C     'NMAX' IS THE ORDER OF THE MATRIX--NOTE 'NMAX' = 1 IS ALLOWED.    00000230
C     'NDIM' IS THE FIRST DIMENSION OF 'ASV', 'BSV', 'A', AND 'B'.      00000240
C     'MMAX' IS THE NUMBER OF SETS OF EQUATIONS (NUMBER OF COLUMNS IN   00000250
C            'B' AND 'BSV').                                            00000260
C     'DET' WILL CONTAIN                                                00000270
C       (SIGN OF DET) * @DET@**POWER                                    00000280
C     THIS IS TO HELP PREVENT UNDER/OVER FLOW IN THE DETERMINANT.       00000290
C                                                                       00000300
C     NOTE THAT THE FORTLIB ROUTINE FIHZ IS USED.  AT RETURN FIHZ WILL  00000310
C      BE ON SO THAT 10 UNDER OR OVER FLOWS WILL CAUSE  TERMINATION OF  00000320
C                                                                       00000330
C     IMPLICIT COMPLEX*16 (A-H, O-Z)                                    00000340
      IMPLICIT REAL*8 ( A-H, O-Z )                                      00000350
C                                                                       00000360
C     NOTE THAT THE SECOND DEMENSION OF B, BSV IS UNNECESSARY           00000370
      DIMENSION A(NDIM, NDIM), ASV(NDIM, NDIM), B(NDIM, MMAX),  BSV(NDIM00000380
     *, MMAX)                                                           00000380
C                                                                       00000400
C     MOVE THE STUFF OVER TO SAVE IT                                    00000410
      IF ((NMAX .LE. 0) .OR. (MMAX .LE. 0) .OR. (NDIM .LT. NMAX))GOTO90000000420
      DO 10 I = 1, NMAX                                                 00000430
      DO 10 J = 1, NMAX                                                 00000440
   10 A(J,I) = ASV(J,I)                                                 00000450
C                                                                       00000460
      ENTRY LINBSV (A, NMAX, NDIM, BSV, MMAX, B, POWER, DET)            00000470
C                                                                       00000480
C     'LINBSV'... THIS ENTRY PRESERVES ONLY THE INHOMOGENEOUS TERM.     00000490
C     'A' AT ENTRY THIS IS THE MATRIX, AT EXIT IT IS AS ABOVE           00000500
C     'BSV' AND 'B' ARE AS ABOVE                                        00000510
C     THE DIMENSIONS AND OTHER ARGS ARE AS ABOVE                        00000520
      IF ((NMAX .LE. 0) .OR. (MMAX .LE. 0) .OR. (NDIM .LT. NMAX))GOTO90000000530
  3   DO 5 M = 1, MMAX                                                  00000540
      DO 5 I = 1, NMAX                                                  00000550
  5   B(I,M) = BSV(I,M)                                                 00000560
C                                                                       00000570
      ENTRY LIN (A, NMAX, NDIM, B, MMAX, POWER, DET)                    00000580
C                                                                       00000590
C     'LIN'... THIS ENRRY SAVES NOTHING                                 00000600
C     'A' IS AS UNDER 'LINBSV'                                          00000610
C     'B' AT ENTRY THIS IS THE INHOMOGENEOUS TERM, AT EXIT IT IS AS     00000620
C         UNDER 'LINSAV'.                                               00000630
C     THE DIMENSIONS AND OTHER ARGUEMENTS ARE ALL AS ABOVE.             00000640
C                                                                       00000650
      IF ((NMAX .LE. 0) .OR. (MMAX .LE. 0) .OR. (NDIM .LT. NMAX))GOTO90000000660
      IF (NMAX .GT. 1) GO TO 100                                        00000670
      DO 50 M = 1, MMAX                                                 00000680
   50 B(1,M) = B(1,M)/A(1,1)                                            00000690
      DET = A(1,1)                                                      00000700
      RETURN                                                            00000710
C                                                                       00000720
C                                                                       00000730
  100 NMAXM1 = NMAX-1                                                   00000740
C                                                                       00000750
      DO 150 K = 1, NMAXM1                                              00000760
      KP1 = K+1                                                         00000770
      DO 150 J = KP1, NMAX                                              00000780
      PIV = -A(J,K)/A(K,K)                                              00000790
      DO 130 M = 1, MMAX                                                00000800
  130 B(J,M) = B(J,M) + PIV*B(K,M)                                      00000810
      DO 150 I = J, NMAX                                                00000820
  150 A(I,J) = A(I,J) + PIV*A(I,K)                                      00000830
C                                                                       00000840
      DO 250 M = 1, MMAX                                                00000850
      B(NMAX,M) = B(NMAX,M)/A(NMAX,NMAX)                                00000860
      DO 250 II = 1, NMAXM1                                             00000870
      I = NMAX-II                                                       00000880
      SUM = 0                                                           00000890
      IP1 = I+1                                                         00000900
      DO 230 J = IP1, NMAX                                              00000910
  230 SUM = SUM + A(J,I)*B(J,M)                                         00000920
  250 B(I,M) = (B(I,M) - SUM)/A(I,I)                                    00000930
C                                                                       00000940
C******* CALLS TO FIHZ TEMPORARILY REMOVED ************                 00000950
C     CALL FIHZ (12, 'OFF')                                             00000960
C     CALL FIHZ (13, 'OFF')                                             00000970
      IF ((POWER .NE. 1.D0) .AND. (POWER .NE. 0)) GO TO 400             00000980
      DET = A(1,1)                                                      00000990
      DO 350 I = 2, NMAX                                                00001000
  350 DET = DET*A(I,I)                                                  00001010
      GO TO 500                                                         00001020
  400 DET = 1.                                                          00001030
      DO 450 I = 1, NMAX                                                00001040
  450 DET = DET*DSIGN(DABS(A(I,I))**POWER, A(I,I))                      00001050
  500 CONTINUE                                                          00001060
C     CALL FIHZ (12, 'ON')                                              00001070
C     CALL FIHZ (13, 'ON')                                              00001080
      RETURN                                                            00001090
C                                                                       00001100
C     HERE A(I,J) IS IN A(J,I+1); I .GE. J.  BSV AND A(J,I+1)           00001110
C       ARE SAVED.                                                      00001120
C    *** NOTE ***  THIS WAS CHANGED FROM SAVING THE LOWER               00001130
C    TRIANGLE TO SAVING THE UPPER TRIANGLE ON 11/06/72 TO MAKE          00001140
C    THINGS FASTER ON THE /195                                          00001150
C                                                                       00001160
      ENTRY LINTRI (A, NMAX, NDIM, BSV, MMAX, B, POWER, DET)            00001170
      IF ((NMAX .LE. 0) .OR. (MMAX .LE. 0) .OR. (NDIM .LT. NMAX)  ) GO T00001180
     *O 900                                                             00001180
      DO 515 I = 1, NMAX                                                00001200
      DO 515 J = 1, I                                                   00001210
 515  A(I,J) = A(J,I+1)                                                 00001220
      GO TO 3                                                           00001230
C                                                                       00001240
  900 PRINT 903, NMAX, NDIM, MMAX                                       00001250
  903 FORMAT ( / '0$$$ LIN/LINSAV -- INVALID INPUT...', 3I10 )          00001260
      CALL SYSERR                                                       00001270
C                                                                       00001280
      END                                                               00001290
      BLOCK DATA LOGDUM                                                 00000100
CNV      BLOCK DATA LOGDUM                                              00000200
C                                                                       00000300
C     BLOCK DATA TO DEFINE THE LOG(FACTORIAL) TABLE.                    00000400
C                                                                       00000500
C     THIS TABLE IS USED BY THE CLEBSH, THREEJ, SIXJ AND RACAH          00000600
C     ROUTINES.                                                         00000700
C                                                                       00000800
C     IT MAY BE INCLUDED IN YOUR LOAD MODULE FOR OTHER PURPOSES         00000900
C     BY PLACING A DUMMY EXTERNAL CARD OF THE FORM                      00001000
C          EXTERNAL  LOGDUM                                             00001100
C     IN YOUR SOURCE.                                                   00001200
C                                                                       00001300
C     THE FORM OF THE TABLE IS                                          00001400
C                                                                       00001500
C     COMMON / LOGFAC /  MAX, LFBIAS, TABLE                             00001600
C                                                                       00001700
C                                                                       00001800
C     MAX IS THE MAXIMUM NUMBER WHOSE LOG(FACTORIAL) IS                 00001900
C         IN TABLE.  TABLE HAS MAX+1 ELEMENTS.                          00002000
C         IN THE PRESENT VERSION:                                       00002100
C            MAX = 1000                                                 00002200
C     ISPACE IS NOT USED.                                               00002300
C     TABLE IS THE TABLE.  THE FIRST ELEMENT CORRESPONDS TO 0.          00002400
C                                                                       00002500
C     6/30/74 - FIRST VERSION - S. PIEPER                               00002600
C     12/28/79 - CDC TO CNI - RPG                                       00002700
C     7/9/80 - USE CDM PREFIX - S.P.                                    00002800
C     7/31/80 - LFBIAS AND ONLY 100 ENTRIES - S.P.                      00002900
C     12/17/80 - ONLY 1 ENTRY SINCE IT IS A FAKE NOW - S.P.             00003000
C                                                                       00003100
C                                                                       00003200
      REAL*8 LF                                                         00003300
C                                                                       00003400
C                                                                       00003500
      COMMON /LOGFAC/ MAXLF, LFBIAS, LF(1)                              00003600
      COMMON /XLOGDUM/ I                                                 00003700
C                                                                       00003800
C                                                                       00004000
      DATA MAXLF /0/, LFBIAS /0/                                        00004100
C                                                                       00004200
      DATA LF(1) / 0.0 0D+00 /                                          00004300
CIB      DATA LF(   1) / Z 0000000000000000 /                           00004400
C                                                                       00004500
      DATA I /0/                                                        00004600
C                                                                       00004800
      END                                                               00004900
      SUBROUTINE LSQPOL (X,Y,W,RESID,NSUB,SUM,LSUB,A,B,MSUB,NMAX,MMAX)  00000010
C                                                                       00000020
C     LEAST SQUARE POLYNOMIAL FIT                                       00000030
C                                                                       00000040
C     12/28/79 - REMOVE EXTRANEOUS IPIVOT - RPG                         00000050
C     2/20/80 - SCALE X TO PREVENT OVERFLOW - RPG                       00000060
C                                                                       00000070
C     X = INPUT ARRAY OF X VALUES, DIMENSION NSUB.                      00000080
C     Y = INPUT ARRAY OF Y VALUES, DIMENSION (NMAX,LSUB).               00000090
C     W = INPUT ARRAY OF WEIGHTS, DIMENSION NSUB.                       00000100
C     RESID = OUTPUT ARRAY OF RESIDUALS (FIT - INPUT Y),                00000110
C        DIMENSION (NMAX,LSUB).                                         00000120
C     NSUB = NUMBER OF X VALUES.                                        00000130
C     SUM = OUTPUT ARRAY FOR WEIGHTED SUMS OF SQUARES,                  00000140
C        DIMENSION LSUB.                                                00000150
C     LSUB = THE NUMBER OF RIGHT-HAND-SIDE VECTORS Y TO BE FITTED.      00000160
C     A = OUTPUT ARRAY FOR THE ERROR MATRIX, DIMENSION (MMAX,MSUB).     00000170
C     B = OUTPUT ARRAY FOR THE COEFFICIENTS OF X**(I-1) IN THE FINAL    00000180
C        FIT, DIMENSION (MMAX,LSUB).                                    00000190
C     MSUB = NUMBER OF TERMS IN THE POLYNOMIAL = DEGREE - 1.            00000200
C        MSUB MUST BE <= 25 UNLESS ARRAY XPOWER (BELOW) IS EXTENDED.    00000210
C     NMAX = FIRST DIMENSION OF Y AND RESID.  MUST BE >= NSUB.          00000220
C     MMAX = FIRST DIMENSION OF A AND B.  MUST BE >= MSUB.              00000230
C                                                                       00000240
C     FOR SINGLE PRECISION- ADD C'S TO THE FOLLOWING 2 CARDS.           00000250
      IMPLICIT REAL*8 ( A-H, O-Z )                                      00000260
CDC      ALLOCLEVEL X,Y,W,RESID,SUM,A,B                                 00000270
      REAL*8 X,Y,W,RESID,SUM,A,B                                        00000280
      DIMENSION X(NSUB),Y(NMAX,LSUB),W(NSUB),RESID(NMAX,LSUB),  SUM(LSUB00000290
     *),A(MMAX,MSUB),B(MMAX,LSUB)                                       00000290
C                                                                       00000310
C     THIS COMMON LIMITS MSUB TO 25.                                    00000320
C                                                                       00000330
      COMMON/F402/XPOWER(100)                                           00000340
      EQUIVALENCE (K1,DETERM),(K2,I2),(TERM,POLY)                       00000350
      N=NSUB                                                            00000360
      L=LSUB                                                            00000370
      M=MSUB                                                            00000380
      M1=M+1                                                            00000390
      M3=M+M+M                                                          00000400
      M31=M3-1                                                          00000410
      M41=M31+M                                                         00000420
C                                                                       00000430
C     SCALE X INTO (-1,1) TO PREVENT OVERFLOW OR UNDERFLOW.             00000440
C                                                                       00000450
      XMAX = 0.                                                         00000460
      DO 20 K = 1, N                                                    00000470
      XMAX = DMAX1( XMAX, DABS(X(K)) )                                  00000480
20         continue
         xscale=1./xmax
      DO 30 K = 1, N                                                    00000510
      X(K) = XSCALE*X(K)                                                00000520
   30 CONTINUE                                                          00000530
C                                                                       00000540
C     FORMATION AND INVERSION OF SYSTEM OF NORMAL EQUATIONS             00000550
C                                                                       00000560
      DO 100 K2=M1,M41                                                  00000570
      XPOWER(K2)=0.0                                                    00000580
  100 CONTINUE                                                          00000590
      DO 200 K1=1,N                                                     00000600
      TERM=W(K1)                                                        00000610
      DO 200 K2=M1,M31                                                  00000620
      XPOWER(K2)=TERM+XPOWER(K2)                                        00000630
      TERM=X(K1)*TERM                                                   00000640
  200 CONTINUE                                                          00000650
      DO 300 I=1,M                                                      00000660
      DO 300 J=1,M                                                      00000670
      K2=I+J+M-1                                                        00000680
      A(I,J)=XPOWER(K2)                                                 00000690
  300 CONTINUE                                                          00000700
      DO 500 J=1,L                                                      00000710
      DO 400 K=1,N                                                      00000720
      TERM=W(K)*Y(K,J)                                                  00000730
      DO 400 K2=M3,M41                                                  00000740
      XPOWER(K2)=TERM+XPOWER(K2)                                        00000750
      TERM=X(K)*TERM                                                    00000760
  400 CONTINUE                                                          00000770
      DO 500 I=1,M                                                      00000780
      K2=I+M31                                                          00000790
      B(I,J)=XPOWER(K2)                                                 00000800
      IF(J.NE.L) XPOWER(K2)=0.0                                         00000810
  500 CONTINUE                                                          00000820
C                                                                       00000830
C     FOR DOUBLE PRECISION- USE DOUBLE PRECISION VERSION OF MATINV      00000840
C                                                                       00000850
      CALL MATINV(A,M,B,L,DETERM,MMAX)                                  00000860
C                                                                       00000870
C     EVALUATION OF RESIDUALS                                           00000880
C                                                                       00000890
      DO 700 J=1,L                                                      00000900
      SUM(J)=0.0                                                        00000910
      DO 700 K=1,N                                                      00000920
      POLY=0.0                                                          00000930
      DO 600 I2=1,M                                                     00000940
      I=M1-I2                                                           00000950
      POLY=X(K)*POLY+B(I,J)                                             00000960
  600 CONTINUE                                                          00000970
      RESID(K,J)=POLY-Y(K,J)                                            00000980
      SUM(J)=SUM(J)+W(K)*RESID(K,J)**2                                  00000990
  700 CONTINUE                                                          00001000
C                                                                       00001010
C     NOW UN-SCALE THE COEFFICIENTS AND X.  DO IT THE SLOW WAY TO       00001020
C     PREVENT UNNECESSARY OVERFLOW OR UNDERFLOW.                        00001030
C                                                                       00001040
      DO 800 J = 1, L                                                   00001050
      DO 800 I2 = 2, M                                                  00001060
      DO 800 I = I2, M                                                  00001070
      B(I,J) = XSCALE*B(I,J)                                            00001080
  800 CONTINUE                                                          00001090
      DO 900 K = 1, N                                                   00001100
      X(K) = XMAX*X(K)                                                  00001110
  900 CONTINUE                                                          00001120
      RETURN                                                            00001130
      END                                                               00001140
      SUBROUTINE MATINV (A,N,B,M,DETERM,NMAX)                           00000100
C                                                                       00000200
C     MATRIX INVERSION WITH ACCOMPANYING SOLUTION OF LINEAR EQUATIONS   00000300
C                                                                       00000400
C     12/28/79 - DETMAX,-MIN FOR CVA,CUN - RPG                          00000410
C                                                                       00000420
      IMPLICIT REAL*8 ( A-H, O-Z )                                      00000500
CDC      ALLOCLEVEL A,B                                                 00000600
      REAL*4 PIVOT                                                      00000700
	DIMENSION A(NMAX,1),B(NMAX,1)
      COMMON /F402/ PIVOT(50), INDEX(50)                                00000900
C     (CARDS BELOW WITH ******** IN COLUMNS 73-80 IMPLEMENT THE IMPROVED00001000
C     COMPUTATION OF THE DETERMINANT IN THE MODIFICATION OF MARCH 1975.)00001100
      REAL*8 DETMAX,DETMIN                                              00001200
CDC      DATA DETMAX,DETMIN/ 1.E+200,1.E-200 /                          00001300
CIB      DATA DETMAX,DETMIN/Z6110000000000000,Z2110000000000000/        00001400
CUN      DATA DETMAX, DETMIN / 1.D200, 1.D-200 /                        00001410
        data DETMAX, DETMIN / 1.D35, 1.D-35 /                           00001420
C                                                                       00001500
C                                                                       00001700
C     INITIALIZE DETERMINANT AND PIVOT ELEMENT ARRAY                    00001800
C                                                                       00001900
      DETERM=1.0                                                        00002000
      IDET=0                                                            00002100
      DO 20 I=1,N                                                       00002200
      PIVOT(I)=0.0                                                      00002300
   20 INDEX(I)=0.0                                                      00002400
C                                                                       00002500
C     PERFORM SUCCESSIVE PIVOT OPERATIONS (GRAND LOOP)                  00002600
C                                                                       00002700
      DO 550 I=1,N                                                      00002800
C                                                                       00002900
C     SEARCH FOR PIVOT ELEMENT AND EXTEND DETERMINANT PARTIAL PRODUCT   00003000
C                                                                       00003100
      AMAX=0.0                                                          00003200
      DO 105 J=1,N                                                      00003300
      IF (PIVOT(J).NE.0.0) GO TO 105                                    00003400
      DO 100 K=1,N                                                      00003500
      IF (PIVOT(K).NE.0.0) GO TO 100                                    00003600
      TEMP=DABS(A(J,K))                                                 00003700
      IF (TEMP.LT.AMAX) GO TO 100                                       00003800
      IROW=J                                                            00003900
      ICOLUM=K                                                          00004000
      AMAX=TEMP                                                         00004100
  100 CONTINUE                                                          00004200
  105 CONTINUE                                                          00004300
      INDEX(I)=4096*IROW+ICOLUM                                         00004400
      J=IROW                                                            00004500
      AMAX=A(J,ICOLUM)                                                  00004600
      DETERM=AMAX*DETERM                                                00004700
      IF (DABS(DETERM).LT.DETMAX) GO TO 130                             00004800
      DETERM=DETERM*DETMIN                                              00004900
      IDET=IDET+1                                                       00005000
      GO TO 140                                                         00005100
  130 IF (DABS(DETERM).GT.DETMIN) GO TO 140                             00005200
      DETERM=DETERM*DETMAX                                              00005300
      IDET=IDET-1                                                       00005400
  140 CONTINUE                                                          00005500
C                                                                       00005600
C     RETURN IF MATRIX IS SINGULAR (ZERO PIVOT) AFTER COLUMN INTERCHANGE00005700
C                                                                       00005800
      IF (DETERM.EQ.0.0) GO TO 600                                      00005900
C                                                                       00006000
      PIVOT(ICOLUM)=AMAX                                                00006100
C                                                                       00006200
C     INTERCHANGE ROWS TO PUT PIVOT ELEMENT ON DIAGONAL                 00006300
C                                                                       00006400
      IF (IROW.EQ.ICOLUM) GO TO 260                                     00006500
      DETERM=-DETERM                                                    00006600
      DO 200 K=1,N                                                      00006700
      SWAP=A(J,K)                                                       00006800
      A(J,K)=A(ICOLUM,K)                                                00006900
      A(ICOLUM,K)=SWAP                                                  00007000
  200 CONTINUE                                                          00007100
      IF (M.LE.0) GO TO 260                                             00007200
      DO 250 K=1,M                                                      00007300
      SWAP=B(J,K)                                                       00007400
      B(J,K)=B(ICOLUM,K)                                                00007500
      B(ICOLUM,K)=SWAP                                                  00007600
  250 CONTINUE                                                          00007700
C                                                                       00007800
C     DIVIDE PIVOT ROW BY PIVOT ELEMENT                                 00007900
C                                                                       00008000
  260 K=ICOLUM                                                          00008100
      A(ICOLUM,K)=1.0                                                   00008200
      DO 350 K=1,N                                                      00008300
      A(ICOLUM,K)=A(ICOLUM,K)/AMAX                                      00008400
  350 CONTINUE                                                          00008500
      IF (M.LE.0) GO TO 380                                             00008600
      DO 370 K=1,M                                                      00008700
      B(ICOLUM,K)=B(ICOLUM,K)/AMAX                                      00008800
  370 CONTINUE                                                          00008900
C                                                                       00009000
C     REDUCE NON-PIVOT ROWS                                             00009100
C                                                                       00009200
  380 DO 550 J=1,N                                                      00009300
      IF (J.EQ.ICOLUM) GO TO 550                                        00009400
      T=A( J,ICOLUM)                                                    00009500
      A( J,ICOLUM)=0.0                                                  00009600
      DO 450 K=1,N                                                      00009700
      A( J,K)=A( J,K)-A(ICOLUM,K)*T                                     00009800
  450 CONTINUE                                                          00009900
      IF (M.LE.0) GO TO 550                                             00010000
      DO 500 K=1,M                                                      00010100
      B( J,K)=B( J,K)-B(ICOLUM,K)*T                                     00010200
  500 CONTINUE                                                          00010300
  550 CONTINUE                                                          00010400
C                                                                       00010500
C     INTERCHANGE COLUMNS AFTER ALL PIVOT OPERATIONS HAVE BEEN PERFORMED00010600
C                                                                       00010700
  600 DO 710 I=1,N                                                      00010800
      I1=N+1-I                                                          00010900
      K=INDEX(I1)/4096                                                  00011000
      ICOLUM=INDEX(I1)-4096*K                                           00011100
      IF (K.EQ.ICOLUM) GO TO 710                                        00011200
      DO 705 J=1,N                                                      00011300
      SWAP=A(J,K)                                                       00011400
      A(J,K)=A(J,ICOLUM)                                                00011500
      A(J,ICOLUM)=SWAP                                                  00011600
  705 CONTINUE                                                          00011700
  710 CONTINUE                                                          00011800
C                                                                       00011900
      PIVOT(1)=IDET                                                     00012000
      RETURN                                                            00012100
      END                                                               00012200
      SUBROUTINE PLMSUB (LMAX, IMMAX, Z, PRAY)                           00000100
C                                                                       00000200
C                                                                       00000300
C                                                                       00000400
C     COMPUTES THE ASSOCIATED LEGENDRE FUNCTIONS OF THE FIRST KIND:     00000500
C            P(L,M)(X)  FOR M = 0, MMAX  ,  L = M, LMAX  AND  REAL X    00000600
C                                                                       00000700
C       THE SIGN CONVENTION OF ABRAMOWITZ AND STEIGEN (M SUPER, L SUB)  00000800
C     AND OF JACKSON ("CLASSICAL ELECTRODYNAMICS") IS USED.             00000900
C     THERE IS NOT A (-1)**M BETWEEN P(L,M) AND Y(L,M) (WHERE Y IS      00001000
C     CONDON-SHORTLEY)  SO THAT  P(L,M) AND Y(L,M)                      00001100
C     HAVE THE SAME SIGN.  NOTE THAT THIS IS DIFFERENT FROM             00001200
C     MESSIAH'S CONVENTION.                                             00001300
C                                                                       00001400
C       FOR  -1 <= X <= 1,  THE "VALUE ON THE CUT" GIVEN BY ABRAMOWITZ  00001500
C     AND STEGUN  8.3.1  IS COMPUTED.  THIS IS RELATED TO P(L,M)(Z) BY  00001600
C          P(L,M)(X)  =  EXP(I .5 M PI) * P(L,M)(Z = X + I 0)           00001700
C     FOR   X  > 1,  P(L,M)(Z=X) IS COMPUTED.  THUS THE OUTPUT IS REAL  00001800
C     FOR ALL REAL X.                                                   00001900
C                                                                       00002000
C                                                                       00002100
C                                                                       00002200
C     T__ _________ ___:                                                00002300
C                                                                       00002400
C     LMAX - THE MAXIMUM L VALUE TO BE FOUND.                           00002500
C                                                                       00002600
C     MMAX - ___ _______ M _____.  MMAX ______ _______                  00002700
C            0 <= MMAX <= LMAX  ___ __  MMAX > LMAX, __________         00002800
C            ____ ________ __ __ MMAX = LMAX  ___ __ _____ _______      00002900
C            ____ __ _______.  H______,  MMAX < 0  ____ ______ __       00003000
C            A FATAL ERROR MESSAGE.  NOTE THAT THE P(L,M) ARE GENERATED 00003100
C            ONLY FOR  M > 0.                                           00003200
C                                                                       00003300
C     Z - THE VALUE OF X ATWHICH P(L,M) IS TO BE FOUND.  SEE ABOVE      00003400
C         FOR THE ACTION TAKEN DEPENDING ON THE RANGE OF Z.             00003500
C                                                                       00003600
C     PLMRAY - A ONE DIMENSIONAL ARRAY OF DIMENSION AT LEAST            00003700
C                 ( (MMAX+1) * (2*LMAX - MMAX + 2) ) / 2                00003800
C              THAT WILL RECIEVE THE P(L, M).                           00003900
C              THEY ARE STORED IN PLMRAY AS (READ L AS LMAX)            00004000
C        P00, P10, P20,... PL0, P11, P21, ... PL1, P22, P32, ... ... PLL00004100
C              THUS                                                     00004200
C                 P(L,M) = PLMRAY(L+1+ M*(2*LMAX+1-M)/2)                00004300
C                                                                       00004400
C     LMAX AND MMAX ARE INTEGER INPUT ARGUMENTS.                        00004500
C     Z IS A REAL INPUT ARGUMENT.                                       00004600
C     PLMRAY IS A REAL OUTPUT ARRAY.                                    00004700
C                                                                       00004800
C                                                                       00004900
C     9/25/72  NEW REVISED (IMPROVED) RECURSION SCHEME - S. PIEPER      00005000
C     10/23/74 - CORRECT OUTPUT FOR  Z  > 1.                            00005100
C     8/1/80 - PROGRAM IT EFFICENTLY - S.P.                             00005200
C     8/5/80 - DIFFERENT RECURSION FOR CRAY - S.P.                      00005300
C                                                                       00005400
C                                                                       00005500
      IMPLICIT REAL*8 ( A-H, O-Z )                                      00005600
CDC      ALLOCLEVEL PRAY                                                00005700
C                                                                       00005800
C                                                                       00005900
      DIMENSION PRAY(1)                                                 00006000
C                                                                       00006100
C                                                                       00006200
CRA      INTEGER WKBIAS                                                 00006300
CRA      REAL INTTBL, INTVLS                                            00006400
CRA      COMMON /INTWRK/  MAXINT, INBIAS, WKBIAS, INTTBL(1)             00006500
CRA      POINTER ( INBASE, INTVLS(1) ), ( WRK1BS, WORK1(1) ),           00006600
CRA     1   ( WRK2BS, WORK2(1) )                                        00006700
C                                                                       00006800
      MMAX=IMMAX
      IF (MMAX .EQ. 0) GO TO 200                                        00006900
C                                                                       00007000
C                                                                       00007100
      ROOT = SQRT(ABS(1-Z**2))                                          00007200
      IF ( ABS(Z) .GT. 1 ) ROOT = -ROOT                                 00007300
CRA      ROOTIN = 0                                                     00007400
CRA      IF ( ROOT .NE. 0 )  ROOTIN = 1/ROOT                            00007500
CRA      MLIM = MMAX                                                    00007600
CRA      IF ( ABS(Z) .LT. 1 )  MLIM = -40/(ALOG(1-ABS(Z))-.01) - 1.5    00007700
C                                                                       00007800
C                                                                       00007900
 200  PRAY(1) = 1.                                                      00008000
C                                                                       00008100
CRA      IF ( LMAX .LE. MAXINT )  GO TO 250                             00008200
CRA      PRINT 223, LMAX, MMAX, MAXINT                                  00008300
CRA 223  FORMAT ( "0***INTERNAL ERROR, INTTBL TOO SMALL FOR PLM", 3I5 ) 00008400
CRA      STOP 2222                                                      00008500
CRAC                                                                    00008600
CRA 250  INBASE = LOCF(INTTBL) + INBIAS                                 00008700
CRA      WRK1BS = LOCF(INTTBL) + WKBIAS                                 00008800
C                                                                       00008900
      I = 1                                                             00009000
      J = 1                                                             00009100
      IZERO = 0                                                         00009200
      DO 499 M = IZERO, MMAX                                            00009300
      IF (M .GE. LMAX) RETURN                                           00009400
      JPRE = J                                                          00009500
      J = I                                                             00009600
C     GENERATE THE P(L=M+1, M) FOR THIS M                               00009700
      I = I+1                                                           00009800
      PRAY(I) = (2*M+1)*Z*PRAY(I-1)                                     00009900
      MP2 = M + 2                                                       00010000
      IF (LMAX .LT. MP2) GO TO 450                                      00010100
C                                                                       00010200
C     NOW DO L = M+2, M+3, ..., LMAX FOR THIS M                         00010300
C                                                                       00010400
      DL = MP2                                                          00010500
      TEMP1 = (2*DL-1)*Z                                                00010600
      TEMP2 = M-1+DL                                                    00010700
      TEMP3 = DL-M                                                      00010800
      IST = I+1                                                         00010900
      IEND = I+1 + LMAX-MP2                                             00011000
C                                                                       00011100
CRA      IF ( M .GT. 0 )  GO TO 350                                     00011200
C                                                                       00011300
      DO 399 I = IST, IEND                                              00011400
      PRAY(I) = ( TEMP1*PRAY(I-1) - TEMP2*PRAY(I-2) ) / TEMP3           00011500
      TEMP1 = TEMP1 + (2.*Z)                                            00011600
      TEMP2 = TEMP2 + 1.                                                00011700
      TEMP3 = TEMP3 + 1.                                                00011800
C                                                                       00011900
 399  CONTINUE                                                          00012000
C                                                                       00012100
CRAC     ON THE CRAY WE RECURSE ON L ONLY FOR M = 0                     00012200
CRAC                                                                    00012300
CRA      WRK2BS = LOCF(PRAY)                                            00012400
CRA      DO 329  L = MP2, LMAX                                          00012500
CRA         WORK1(L+1) = (2*INTVLS(L+1)-1.)*Z/INTVLS(L+1)               00012600
CRA         WORK2(L+1) = (INTVLS(L+1)-1.)/INTVLS(L+1)                   00012700
CRA 329  CONTINUE                                                       00012800
CRAC                                                                    00012900
CRA      DO 339  I = IST, IEND                                          00013000
CRAC                                                                    00013100
CRAC     HERE, SINCE M=0, I = L+1                                       00013200
CRAC                                                                    00013300
CRA         PRAY(I) = WORK1(I)*PRAY(I-1) - WORK2(I)*PRAY(I-2)           00013400
CRA 339  CONTINUE                                                       00013500
CRA      GO TO 400                                                      00013600
CRAC                                                                    00013700
CRAC     ON THE CRAY, FOR M > 0 WE RECURSE ON M TO GET VECTOR MODE.     00013800
CRAC     THIS IS, HOWEVER, UNSTABLE FOR !M! NEAR 1.                     00013900
CRAC                                                                    00014000
CRA 350  IF ( M .GT. MLIM )  GO TO 370                                  00014100
CRAC                                                                    00014200
CRAC     HERE                                                           00014300
CRAC                                                                    00014400
CRAC     PRAY(L+JPRE-M+1) = P(L,M-1)                                    00014500
CRAC     PRAY(L+IST-M-2)  = P(L,M)   IS COMPUTED                        00014600
CRAC                                                                    00014700
CRACDIR$ IVDEP                                                          00014800
CRA      MOFF = IST-MP2                                                 00014900
CRA      MM1OFF = JPRE-M+1                                              00015000
CRA      TDMM1 = 2*(M-1)                                                00015100
CRA      DO 359  L = MP2, LMAX                                          00015200
CRA         PRAY(L+MOFF) = ROOTIN*( INTVLS(L-M+2)*Z*PRAY(L+MM1OFF)      00015300
CRA     1      - (INTVLS(L-M+2)+TDMM1)*PRAY(L+MM1OFF-1) )               00015400
CRA 359  CONTINUE                                                       00015500
CRA      GO TO 400                                                      00015600
CRAC                                                                    00015700
CRAC     HERE EVEN ON CRAY WE MUST RECURSE ON L DUE TO                  00015800
CRAC     M-RECURSION BEING UNSTABLE                                     00015900
CRAC                                                                    00016000
CRA 370  WRK2BS = LOCF(PRAY) + IST-1 - MP2                              00016100
CRA      DM = M                                                         00016200
CRA      DO 379  L = MP2, LMAX                                          00016300
CRA         WORK1(L+1) = (2*INTVLS(L+1)-1.)*Z/(INTVLS(L+1)-DM)          00016400
CRA         WORK2(L+1) = (INTVLS(L+1)+(DM-1.))/(INTVLS(L+1)-DM)         00016500
CRA 379  CONTINUE                                                       00016600
CRA      IOFF = MP2-IST+1                                               00016700
CRA      DO 389  I = IST, IEND                                          00016800
CRA         PRAY(I) = WORK1(I+IOFF)*PRAY(I-1) - WORK2(I+IOFF)*PRAY(I-2) 00016900
CRA 389  CONTINUE                                                       00017000
C                                                                       00017100
 400  I = IEND                                                          00017200
C                                                                       00017300
 450  IF (M .EQ. MMAX) RETURN                                           00017400
C                                                                       00017500
C    NOW GENERATE THE NEXT P(L=M, M)                                    00017600
C    NOTE THAT HERE THE DESIRED M VALUE IS REALLY M+1 SINCE             00017700
C     THE DO LOOP HAS NOT YET INCREMENTED M                             00017800
C                                                                       00017900
      I = I+1                                                           00018000
C     J POINTS TO THE LAST P(M, M)                                      00018100

      PRAY(I) = - (2*M+1)*ROOT * PRAY(J)                                00018200
 499  CONTINUE                                                          00018300
      RETURN                                                            00018400
C                                                                       00018500
C     THIS ENTRY COMPUTES THE LEGENDRE POLYNOMIALS                      00018600
C                                                                       00018700
      ENTRY PLSUB (LMAX, Z, PRAY)                                       00018800
      MMAX = 0                                                          00018900
      GO TO 200                                                         00019000
C                                                                       00019100
      END                                                               00019200
      SUBROUTINE RCWFN ( RHO,ETA,MINL,MAXL,FC,FCP,GC,GCP,ACCUR, IRET )  00000010
C                                                                       00000020
C     COMPUTES REGULAR AND IRREGULAR COULOMB WAVEFUNCTIONS.             00000030
C                                                                       00000040
C       THIS SUBROUTINE RETURNS THE REGULAR AND IRREGULAR COULOMB       00000050
C     WAVEFUNCTIONS AND THEIR DERIVATIVES FOR SPECIFIC VALUES OF        00000060
C     ETA AND RHO AND FOR A RANGE OF L'S.  THE RANGE OF L'S (ORBITAL    00000070
C     ANGULAR MOMENTUM) NEED NOT BEGIN AT L = 0 AND, IN GENERAL,        00000080
C     IT WILL BE LESS TIME CONSUMING TO REQUEST ONLY THOSE L'S          00000090
C     THAT ARE ACTUALLY DESIRED.  THE CONVENTIONS OF ABRAMOWITZ         00000100
C     AND STEGUN (HANDBOOK OF MATHEMATICAL FUNCTIONS) ARE USED.         00000110
C                                                                       00000120
C       THIS SUBROUTINE IS AN ADAPTATION OF THE MANCHESTER SUBROUTINE   00000130
C     OF THE SAME NAME (BUT SLIGHTLY DIFFERENT ARGUMENT LIST) THAT      00000140
C     WAS PUBLISHED IN COMPUTER PHYSICS COMMUNICATIONS 8, 377 (1974).   00000150
C     BY BARNETT, FENG, STEED AND GOLDFARB.  THE CONTINUED FRACTIONS    00000160
C     INTRODUCED IN THAT ARTICLE ARE USED IN THE RANGE RHO >            00000170
C     RHO(TURN) BUT COMPLETELY DIFFERENT PROCEDURES ARE USED            00000180
C     FOR RHO < RHO(TURN).  THE NEW PROCEDURES ALLOW RESULTS            00000190
C     FOR RHO < RHO(TURN) OF ACCURACY COMPARABLE TO THOSE FOR           00000200
C     RHO > RHO(TURN) IN COMPARABLE TIME.                               00000210
C                                                                       00000220
C       ONE OF THE CONTINUED FRACTIONS IS NOT CONVERGENT FOR VERY       00000230
C     SMALL RHO AND THUS THIS ROUTINE WILL FIND ONLY F AND F' FOR       00000240
C     RHO < .005.  IN ADDITION THE ROUTINE BECOMES RATHER SLOW FOR      00000250
C     RHO < .1.  FOR VERY LARGE RHO ( >> 1000 ), ANOTHER CONTINUED      00000260
C     FRACTION CONVERGES VERY SLOWLY AND THE SUBROUTINE WILL FAIL       00000270
C     TO CONVERGE IN THE MAXIMUM ALLOWED NUMBER OF ITERATIONS FOR       00000280
C     RHO > 9500 WITH THE EXCEPTION THAT LARGER RHO ARE POSSIBLE        00000290
C     IF ETA IS NEAR RHO/2.                                             00000300
C                                                                       00000310
C       ASIDE FROM THE ABOVE LIMITATIONS, THE SUBROUTINE HAS A VERY     00000320
C     LARGE RANGE OF ETA,  RHO, AND L FOR WHICH IT WILL RETURN          00000330
C     RELIABLE VALUES IN REASONABLE TIMES.  IT HAS BEEN TESTED FOR      00000340
C     -2000 < ETA < 5000,  1E-6 < RHO < 10000;  0 < L < 1000 .          00000350
C     THE TESTS CONSISTED OF COMPARASONS TO EXISTING ROUTINES TO        00000360
C     LOCATE CODING ERRORS AND OF COMPARASONS TO QUADRUPLE PRECISION    00000370
C     RESULTS TO DETERMINE THE NUMERICAL STABILITY OF THE ALGORITHMS.   00000380
C     HOWEVER, NO VERIFICATION OF THE CORRECTNESS OF G OR G' FOR        00000390
C     ETA < 0 HAS BEEN MADE.                                            00000400
C                                                                       00000410
C                                                                       00000420
C     ARGUMENTS (ALL FLOATING POINT ARGUMENTS ARE DOUBLE PRECISION) -   00000430
C                                                                       00000440
C     RHO - THE VALUE OF THE RADIAL COORDINATE AT WHICH F, F', G, AND   00000450
C           G' ARE DESIRED.  0 < RHO < 9500 IS REQUIRED.  IF            00000460
C           RHO < .005, ONLY F AND F' WILL BE FOUND.                    00000470
C                                                                       00000480
C     ETA - THE VALUE OF THE CHARGE PARAMETER TO BE USED.  THE ROUTINE  00000490
C           HAS BEEN TESTED FOR -2000 < ETA < 5000.  ETA = 0 WILL       00000500
C           RESULT IN THE COMPUTATION OF THE SPHERICAL BESSEL           00000510
C           FUNCTIONS TIMES RHO.                                        00000520
C                                                                       00000530
C     MINL - THE MINIMUM VALUE OF L FOR WHICH THE FUNCTIONS ARE DESIRED.00000540
C     MAXL - THE MAXIMUM VALUE OF L.                                    00000550
C                                                                       00000560
C     FC, FCP, GC, AND GCP - THESE ONE DIMENSIONAL ARRAYS WILL BE       00000570
C           SET TO THE COMPUTED VALUES OF F, F', G, AND G' RESPECTIVELY.00000580
C           EACH ARRAY MUST BE OF LENGTH AT LEAST MAXL+1 AND WILL       00000590
C           BE SET AS                                                   00000600
C             ARRAY(L+1) = FUNCTION(L)   MINL =< L =< MAXL .            00000610
C           DEPENDING ON ETA AND RHO, THE FIRST MINL ELEMENTS OF        00000620
C           EACH ARRAY MAY BE USED FOR INTERMEDIATE COMPUTATIONS AND    00000630
C           THEIR VALUES ARE NOT PREDICTABLE UPON EXIT FROM RCWFN.      00000640
C                                                                       00000650
C     ACCUR - THE DESIRED ACCURACY OF CONVERGENCE OF THE CONTINUED      00000660
C          FRACTIONS AND OTHER SERIES.  IN GENERAL THIS WILL BE         00000670
C          THE RELATIVE ACCURACY OF THE FINAL VALUES EXCEPT NEAR        00000680
C          A ZERO OF ONE OF THE FUNCTIONS.  ACCUR SHOULD BE IN          00000690
C          THE RANGE                                                    00000700
C              1E-6 < ACCUR < 1E-15 ,                                   00000710
C          IF IT IS NOT, ONE OF THE ABOVE TWO VALUES WILL BE USED.      00000720
C                                                                       00000730
C     IRET - THIS INTEGER IS SET TO A RETURN CODE TO INDICATE THAT      00000740
C          RCWFN WAS SUCCESSFUL OR THE REASON FOR FAILURE:              00000750
C        0 - ALL O.K.; ALL FUNCTIONS FOUND.                             00000760
C        1 - SOME O.K., THOSE FOR HIGHER L'S HAVE UNDER/OVER-FLOWED     00000770
C        2 - RHO < .005; ONLY F, F' WERE COMPUTED.                      00000780
C        3 - RHO < .005; ONLY F, F' WERE COMPUTED AND SOME OF THEM      00000790
C                        UNDERFLOWED.                                   00000800
C        4 - ALL WILL UNDER/OVER-FLOW; NO RESULTS DEFINED.              00000810
C        5 - FAILED TO AVOID DIVIDE BY ZERO IN F'/F LOOP                00000820
C        6 - NONCONVERGENCE OF R                                        00000830
C        7 - NONCONVERGENCE OF P + IQ                                   00000840
C        8 - NONCONVERGENCE OF MACLAUREN SERIES                         00000850
C        9 - NONCONVERGENCE OF THE TAYLOR SERIES ABOUT RHO = 2*ETA      00000860
C     IRET'S OF 5 TO 9 SHOULD BE REPORTED TO STEVE PIEPER ( X4523 )     00000870
C                                                                       00000880
C       THE FOLLOWING TABLE GIVES SAMPLE EXECUTION TIMES ON THE         00000890
C     /75 AND /195.  TIMES ARE FOR ACCUR = 1D-14 AND FOR                00000900
C     MINL = MAXL = 0.  IN GENERAL EXECUTION TIME IS WEAKLY DEPENDANT   00000910
C     ON ACCUR AND MINL.  IF MAXL > MINL, THE TIMES REQUIRED FOR        00000920
C     EACH ADDITIONAL L ARE APPROXIMATLEY 0.01 MILLISEC ON THE /195     00000930
C     AND 0.15 MILLISEC ON THE /75.  TIMES FOR ETA < 0 ARE COMPARABLE TO00000940
C     TO THOSE FOR \ETA\.  THE FINAL COLUMN OF THE TABLE GIVES          00000950
C     THE NUMBER OF BITS LOST DUE TO ROUND OFF AND TRUNCATION           00000960
C     ERRORS AND INDICATES THE MAXIMUM PRECISION POSSIBLE ON A          00000970
C     GIVEN MACHINE.  THE /360 AND /370 HAVE 56 BITS OF PRECISION       00000980
C     (16.8 DECIMAL PLACES) AND EACH 3.3 BITS LOST REPRESENTS           00000990
C     ONE DECIMAL PLACE LOST.                                           00001000
C                                                                       00001010
C      ETA    RHO     TIME IN MILLISECONDS    BITS LOST                 00001020
C                       /195      /75                                   00001030
C                                                                       00001040
C        0.      .01     .15       1.2            3                     00001050
C        0.     1.       .13       1.2            3                     00001060
C        0.   100.      1.0       11.            11                     00001070
C        0.  1000.      7.6       81.                                   00001080
C        1.      .1     5.4       81.             9                     00001090
C        1.     1.      1.0       12.5            8                     00001100
C        1.    10.       .33       3.7            3                     00001110
C        1.   100.      1.1       11.                                   00001120
C        1.  1000.      8.1       81.                                   00001130
C       10.     1.      1.5       20.             6                     00001140
C       10.    10.       .8        8.3            6                     00001150
C      100.    75.      1.9       21.            10                     00001160
C      100.   100.      2.0       21.            10                     00001170
C      100.  1000.      7.3       74.            13                     00001180
C     1000.  5000.     29.       385.            16                     00001190
C                                                                       00001200
C                                                                       00001210
C     MAY 23, 1976 - REVISED VERSION BY S. PIEPER.                      00001220
C     NOV 19, 1976 - FIX CHOICE OF ROOT FOR RHO << ETA, ETA > 10        00001230
C     12/14/77 - CDC VERSION USING MORTRAN                              00001240
C     6/27/78 - FIX CDC VRYBIG, PRINTOUTS FOR ERROR CONDITIONS - S.P.   00001250
C     1/29/80 - CUN, CVA VERSIONS, CFR%F - RPG                          00001260
C     3/13/80 - BACK TO CFR%F - S.P.                                    00001261
C     7/9/80 - CRAY DATA STATEMENST - S.P.                              00001262
C                                                                       00001270
C                                                                       00001280
      IMPLICIT REAL*8 ( A-H, O-Z )                                      00001290
C                                                                       00001300
CDC      ALLOCLEVEL FC, FCP, GC, GCP                                    00001310
      DIMENSION FC(1),FCP(1),GC(1),GCP(1)                               00001320
C                                                                       00001330
C     VRYBIG IS THE RESULT OF AN OVERFLOW                               00001340
C     BIG IS REPRESENTATIVE OF MAX NUMBER, SMALL IS ITS INVERSE         00001350
C     SMALLN IS  LOG(SMALL)  (MUST BE ACCURATE)                         00001360
C     PRECIS IS ABOUT  100*(MACHINE PRECISION)                          00001370
C     PRECLN IS  > -\LOG(MACHINE PRECISION)\                            00001380
C     PRERT3 IS LIKE THE CUBEROOT OF THE MACHINE PRECISION              00001390
C                                                                       00001400
CIB      DATA VRYBIG / Z 7FFFFFFF FFFFFFFF /                            00001410
CIB      DATA BIG / Z 70800000 00000000 /                               00001420
CIB      DATA SMALL / Z 11200000 00000000 /                             00001430
C     REAL*16 SMALLN / -132.39111148694955409869133519851171 Q+00 /     00001440
CIB      REAL*8 SMALLN / Z C284641FE1E589CC /                           00001450
CIB      DATA PRERT3 / 1.D-7 /,  PRECIS / 1.D-15 /,  PRECLN / 44.D0 /   00001460
C                                                                       00001470
C     FOR CDC WE STOP BEFORE OVERFLOW SINCE CANNOT TRAP OVERFLOW        00001480
CDC      DATA VRYBIG / 3657 4000 0000 0000 0000 B /                     00001490
CDC      DATA BIG / 3621 4000 0000 0000 0000 B /                        00001500
CDC      DATA SMALL  / 0020 4000 0000 0000 0000 B /                     00001510
CDC      DATA SMALLN / 6046 2632 2411 4172 5514 B /                     00001520
CDC      DATA PRERT3 / 1.E-5 /,  PRECIS /1.E-12 /,                      00001522
CDC     1  PRECLN / 35.E0 /                                             00001524
C                                                                       00001526
CRA      DATA VRYBIG / 1.E+2450 /                                       00001530
CRA      DATA BIG    / 1.E+2100 /                                       00001540
CRA      DATA SMALL  / 1.E-2100 /                                       00001550
CRA      DATA SMALLN / -4835.4286 9528 749 /                            00001560
CRA      DATA PRERT3 / 1.E-5 /,  PRECIS / 1.E-12 /,                     00001570
CRA     1  PRECLN / 35.E0 /                                             00001580
C                                                                       00001590
CUN      DATA VRYBIG / 1.D290 /, BIG / 1.D250 /, SMALL / 1.D-250 /      00001600
CUN      DATA SMALLN / -575.64627 32487 61421 0D0 /                     00001610
CUN      DATA PRERT3 / 1.D-7 /, PRECIS / 1.D-16 /, PRECLN / 44.D0 /     00001620
C                                                                       00001630
        data VRYBIG / 1.D32 /, BIG / 1.D25 /, SMALL / 1.D-25 /          00001640
        data SMALLN / -57.564627 32487 61421 0D0 /                      00001650
        data PRERT3 / 1.D-7 /, PRECIS / 1.D-15 /, PRECLN / 44.D0 /      00001660
C                                                                       00001670
CXX      REAL*16 PI / 3.1415926535897932384626433832795028 0Q+00 /      00001680
      REAL*8 PI                                                         00001690
        data PI / 3.141592653589793238 0D+00 /                          00001690
CIB      DATA PI / Z 413243F6A8885A31 /                                 00001700
C                                                                       00001710
      LOGICAL FRSTSW                                                    00001720
C                                                                       00001730
CFR%F                                                                   00001750
C                                                                       00001760
C     COULOMB WAVEFUNCTIONS CALCULATED AT RHO = RHO BY THE              00001770
C     CONTINUED-FRACTION METHOD OF STEED   MINL,MAXL ARE ACTUAL L-VALUES00001780
C     SEE BARNETT FENG STEED AND GOLDFARB COMPUTER PHYSICS COMMON 1974  00001790
C                                                                       00001800
C                                                                       00001810
C     HERE WE LIMIT ACCURACY TO REASONABLE VALUES FOR THE MACHINE       00001820
C                                                                       00001830
      ACC = ACCUR                                                       00001840
      ACC = DMAX1( ACC, PRECIS )                                        00001850
      ACC = DMIN1( ACC, PRERT3 )                                        00001860
C                                                                       00001870
      LMAX = MAXL                                                       00001880
      LMIN = MINL                                                       00001890
      LMIN1= LMIN + 1                                                   00001900
      XLL1 = LMIN*LMIN1                                                 00001910
      ETA2 = ETA*ETA                                                    00001920
C                                                                       00001930
C     DETERMIN WHICH REGION WE ARE IN                                   00001940
C                                                                       00001950
C     FOR RHO < .45, Q OF P+IQ IS POORLY DETERMINED SO WE DON'T USE IT. 00001960
C     EXCEPT THAT FOR LARGE NEGATIVE ETA THE MACLAUREN SERIES ALSO      00001970
C     HAS PROBLEMS                                                      00001980
C                                                                       00001990
      IF ( RHO .GT. .45D0 ) GO TO 20                                    00002000
      IF ( ETA .GE. 0 ) GO TO 10                                        00002010
      IF ( -ETA*RHO .GT. 7 ) GO TO 20                                   00002020
C                                                                       00002030
C     FOR RHO < .005, WE ONLY RETURN F AND F' SINCE THE P+IQ RECURSION  00002040
C     IS VERY SLOWLY CONVERGENT  ( FOR VERY SMALL ETA IT IS POSSIBLE    00002050
C     TO GO TO SMALLER RHO ( RHO > 1E+4*ETA )  BUT WE IGNOR THAT HERE.  00002060
C                                                                       00002070
 10   IF ( RHO .GT. .005D0 ) GO TO 60                                   00002080
      IGOTO = 5                                                         00002090
      GO TO 70                                                          00002100
C                                                                       00002110
 20   TURN = ETA + DSQRT(ETA2 + XLL1)                                   00002120
      IGOTO = 1                                                         00002130
      IF ( RHO .GE. TURN-1.D-4 ) GO TO 100                              00002140
C                                                                       00002150
C     WE ARE INSIDE THE TURNING POINT FOR MINL, CAN WE GET OUTSIDE      00002160
C     OF IT BY REDUCING MINL. (THIS IS ALWAYS POSSIBLE FOR              00002170
C     ETA < 0).                                                         00002180
C                                                                       00002190
      IF ( RHO .LT. ETA+DABS(ETA) ) GO TO 60                            00002200
C                                                                       00002210
C     YES, DO SO - THIS IS THE SAME AS  RHO > RHO(TURN)  EXCEPT         00002220
C     WE GENERATE SOME EXTRA F(L), G(L) FOR  L < MINL                   00002230
C                                                                       00002240
      LMIN = .5*( DSQRT(1+4*((RHO-ETA)**2-ETA2)) - 1 )                  00002250
      LMIN1 = LMIN + 1                                                  00002260
      GO TO 80                                                          00002270
C                                                                       00002280
C     MUST USE A DIFFERENT METHOD TO SUPPLIMENT THE BAD I Q             00002290
C     VALUE.  ALWAYS START WITH LMIN = 0 FOR SIMPLICITY.                00002300
C                                                                       00002310
C     NOTE ONLY ETA > 0 GETS TO HERE ( EXCEPT WHEN RHO < .45 )          00002320
C                                                                       00002330
 60   IGOTO = 2                                                         00002340
 70   LMIN = 0                                                          00002350
      LMIN1 = 1                                                         00002360
      IF ( ETA .LT. 10 .OR. RHO .LE. ETA ) GO TO 80                     00002370
      IGOTO = 3                                                         00002380
C                                                                       00002390
 80   XLL1 = LMIN*LMIN1                                                 00002400
C                                                                       00002410
C     HERE WE COMPUTE  F'/F  FOR L = MAXL                               00002420
C     WE THEN RECURSE DOWN TO LMIN TO GENERATE THE UNNORMALIZED F'S     00002430
C     THIS SECTION IS USED FOR ALL RHO.                                 00002440
C                                                                       00002450
 100  PL = LMAX + 1                                                     00002460
      RHOUSE = RHO                                                      00002470
 105  PLSAVE = PL                                                       00002480
 110  FRSTSW = .TRUE.                                                   00002490
C     CONTINUED FRACTION FOR  R = FP(MAXL)/F(MAXL)                      00002500
      R = ETA/PL + PL/RHOUSE                                            00002510
      DQ = (ETA*RHOUSE)*2.0 + 6*PL**2                                   00002520
      DR = 12*PL + 6                                                    00002530
      DEL = 0.0                                                         00002540
      D = 0.0                                                           00002550
      F = 1.0                                                           00002560
      X = (PL*PL - PL + (ETA*RHOUSE))*(2.0*PL - 1.0)                    00002570
      AI = RHOUSE*PL**2                                                 00002580
      DI = (2*PL+1)*RHOUSE                                              00002590
C                                                                       00002600
C     LOOP AND CONVERGE ON R                                            00002610
C                                                                       00002620
      DO 139 I = 1, 100000                                              00002630
      H = (AI + RHOUSE*ETA2)*(RHOUSE - AI)                              00002640
      X = X + DQ                                                        00002650
      D = D*H + X                                                       00002660
C                                                                       00002670
C     IF WE PASS NEAR A ZERO OF THE DIVISOR, START OVER AT              00002680
C     LARGER LMAX                                                       00002690
C                                                                       00002700
      IF ( DABS(D) .GT. PRERT3*DABS(DR) ) GO TO 130                     00002710
      PL = PL + 1                                                       00002720
      IF ( PL .LT. PLSAVE+10 ) GO TO 110                                00002730
      IRET = 5                                                          00002740
      GO TO 990                                                         00002750
C                                                                       00002760
 130  D = 1/D                                                           00002770
      DQ = DQ + DR                                                      00002780
      DR = DR + 12                                                      00002790
      AI = AI + DI                                                      00002800
      DI = DI + 2*RHOUSE                                                00002810
      DEL = DEL*(D*X - 1.0)                                             00002820
      IF (FRSTSW) DEL = -RHOUSE*(PL*PL + ETA2)*(PL + 1.0)*D/PL          00002830
      FRSTSW = .FALSE.                                                  00002840
      R = R + DEL                                                       00002850
      IF(D.LT.0.0) F = -F                                               00002860
      IF ( DABS(DEL) .LT. DABS(R*ACC) ) GO TO 140                       00002870
 139  CONTINUE                                                          00002880
      IRET = 6                                                          00002890
      GO TO 990                                                         00002900
C                                                                       00002910
C     R HAS CONVERGED;  DID WE INCREASE LMAX                            00002920
C                                                                       00002930
 140  IF ( PL .EQ. PLSAVE ) GO TO 160                                   00002940
C                                                                       00002950
C     RECURSE DOWN ON R TO LMAX                                         00002960
C     HERE THE ONLY PART OF F THAT IS OF INTEREST IS THE SIGN           00002970
C                                                                       00002980
      PL = PL-1                                                         00002990
 150  D = ETA/PL + PL/RHOUSE                                            00003000
      F = (R+D)*F                                                       00003010
      R = D - (1+ETA2/PL**2)/(R+D)                                      00003020
      PL = PL - 1                                                       00003030
      IF ( PL .GT. PLSAVE ) GO TO 150                                   00003040
C                                                                       00003050
C     NOW HAVE R(LMAX, RHO) OR IF IGOTO=4, R(LMIN, 2*ETA)               00003060
C                                                                       00003070
 160  IF ( IGOTO .EQ. 4 ) GO TO 210                                     00003080
      FC (LMAX+1) = F                                                   00003090
      FCP(LMAX+1) = F*R                                                 00003100
      IF( LMAX.EQ.LMIN) GO TO 200                                       00003110
C     DOWNWARD RECURSION TO LMIN FOR F AND FP, ARRAYS GC,GCP ARE STORAGE00003120
      L = LMAX                                                          00003130
      PL = LMAX                                                         00003140
      AR = 1/RHO                                                        00003150
      DO 189 LP = LMIN1,LMAX                                            00003160
      GC (L+1) = ETA/PL + PL*AR                                         00003170
      GCP(L+1) = DSQRT( (ETA/PL)**2 + 1 )                               00003180
      FC (L) = (GC(L+1)*FC(L+1) + FCP(L+1))/GCP(L+1)                    00003190
      FCP(L) = GC(L+1)*FC(L) - GCP(L+1)*FC(L+1)                         00003200
      PL = PL - 1                                                       00003210
      L = L - 1                                                         00003220
C                                                                       00003230
C     IF WE ARE GETTING NEAR AN OVERFLOW, RENORMALIZE EVERYTHING DOWN   00003240
C                                                                       00003250
      IF ( DABS(FC(L+1)) .LT. BIG ) GO TO 189                           00003260
      DO 179 LL = L, LMAX                                               00003270
      FC(LL+1) = SMALL*FC(LL+1)                                         00003280
      FCP(LL+1) = SMALL*FCP(LL+1)                                       00003290
 179  CONTINUE                                                          00003300
 189  CONTINUE                                                          00003310
      F = FC (LMIN1)                                                    00003320
      R = FCP(LMIN1)/F                                                  00003330
C                                                                       00003340
C     HERE WE FIND                                                      00003350
C        P + IQ  =  (G'+IF')/(G+IF)                                     00003360
C     THIS SECTION IS USED IN ALL CASES EXCEPT WHEN                     00003370
C        15 < ETA < RHO < 2*ETA                                         00003380
C                                                                       00003390
 200  IF ( IGOTO .EQ. 3 ) GO TO 500                                     00003400
      IF ( IGOTO .EQ. 5 ) GO TO 400                                     00003410
C                                                                       00003420
C     NOW OBTAIN P + I.Q FOR LMIN FROM CONTINUED FRACTION (32)          00003430
C     REAL ARITHMETIC TO FACILITATE CONVERSION TO IBM USING REAL*8      00003440
 210  P = 0.0                                                           00003450
      Q = RHOUSE - ETA                                                  00003460
      PL = 0.0                                                          00003470
      AR = -(ETA2 + XLL1)                                               00003480
      AI = ETA                                                          00003490
      BR = Q + Q                                                        00003500
      BI = 2.0                                                          00003510
      WI = ETA + ETA                                                    00003520
      DR = BR/(BR*BR + BI*BI)                                           00003530
      DI = -BI/(BR*BR + BI*BI)                                          00003540
      DP = -(AR*DI + AI*DR)                                             00003550
      DQ = (AR*DR - AI*DI)                                              00003560
C                                                                       00003570
C     LOOP AND CONVERGE ON P + IQ                                       00003580
C                                                                       00003590
 230  P = P + DP                                                        00003600
      Q = Q + DQ                                                        00003610
      PL = PL + 2.0                                                     00003620
      AR = AR + PL                                                      00003630
      AI = AI + WI                                                      00003640
      BI = BI + 2.0                                                     00003650
      D = AR*DR - AI*DI + BR                                            00003660
      DI = AI*DR + AR*DI + BI                                           00003670
      T = 1.0/(D*D + DI*DI)                                             00003680
      DR = T*D                                                          00003690
      DI = -T*DI                                                        00003700
      H = BR*DR - BI*DI - 1.0                                           00003710
      X = BI*DR + BR*DI                                                 00003720
      T = DP*H - DQ*X                                                   00003730
      DQ = DP*X + DQ*H                                                  00003740
      DP = T                                                            00003750
      IF(PL.GT.46000.) GO TO 920                                        00003760
      IF(DABS(DP)+DABS(DQ).GE.(DABS(P)+DABS(Q))*ACC) GO TO 230          00003770
      P = P/RHOUSE                                                      00003780
      Q = Q/RHOUSE                                                      00003790
C                                                                       00003800
C     WE NOW HAVE  R  AND  P+IQ,  IS THIS ENOUGH                        00003810
C                                                                       00003820
      IF ( IGOTO .EQ. 2 ) GO TO 400                                     00003830
C                                                                       00003840
C     SOLVE FOR FP,G,GP AND NORMALISE F  AT L=LMIN                      00003850
C                                                                       00003860
C     SINCE THIS IS FOR  RHO > RHO(TURN), F AND G ARE REASONABLE        00003870
C     NUMBERS                                                           00003880
C                                                                       00003890
      X = (R-P)/Q                                                       00003900
      FMAG = DSQRT( 1/(Q*(1+X**2)) )                                    00003910
      W = FMAG/DABS(F)                                                  00003920
      F = W*F                                                           00003930
      G = F*X                                                           00003940
      GP = R*G - 1/F                                                    00003950
      IF ( IGOTO .EQ. 4 ) GO TO 600                                     00003960
      GO TO 800                                                         00003970
C                                                                       00003980
C     HERE   RHO < ETA  OR  RHO < 2*ETA < 20  OR  RHO < .45             00003990
C     WE USE THE MACLAUREN SERIES TO GET  F( L=0, ETA, RHO )            00004000
C                                                                       00004010
C     FIRST COMPUTE  RHO*C(L=0, ETA)                                    00004020
C                                                                       00004030
 400  C = 2*PI*ETA                                                      00004040
      IF ( DABS(C) .GT. .5 ) GO TO 410                                  00004050
C                                                                       00004060
C     USE MACLAURIN EXPANSION OF  X / (EXP(X)-1)                        00004070
C                                                                       00004080
      X = 0                                                             00004090
      T = 1                                                             00004100
      AR = 1                                                            00004110
      BR = C                                                            00004120
      AI = 1                                                            00004130
      C = 1                                                             00004140
 405  AI = AI + 1                                                       00004150
      AR = AR*BR/AI                                                     00004160
      C = C + AR                                                        00004170
      IF ( DABS(AR) .GE. ACC*C ) GO TO 405                              00004180
      C = 1/C                                                           00004190
      GO TO 430                                                         00004200
C                                                                       00004210
C     HERE ETA IS NOT TINY.                                             00004220
C                                                                       00004230
 410  IF ( ETA .GT. 0 ) GO TO 420                                       00004240
      C = -C                                                            00004250
      X = 0                                                             00004260
      T = 1                                                             00004270
      GO TO 425                                                         00004280
 420  X = -SMALLN - PI*ETA                                              00004290
      T = SMALL                                                         00004300
 425  IF ( C .LT. PRECLN ) C = C / (1-DEXP(-C))                         00004310
 430  C = RHO*DSQRT(C)                                                  00004320
      B1 = 1                                                            00004330
      B2 = ETA*RHO                                                      00004340
      SUM = B1 + B2                                                     00004350
      AI = 6                                                            00004360
      DI = 6                                                            00004370
      DO 449 I = 1, 10000                                               00004380
      B3 = ( (2*ETA*RHO)*B2 - (RHO**2)*B1 ) / AI                        00004390
      AI = AI + DI                                                      00004400
      DI = DI + 2                                                       00004410
      SUM = SUM + B3                                                    00004420
      STOP = DABS(B1) + DABS(B2) + DABS(B3)                             00004430
      B1 = B2                                                           00004440
      B2 = B3                                                           00004450
      IF ( DABS(SUM) .LT. BIG ) GO TO 445                               00004460
      X = X - SMALLN                                                    00004470
      SUM = SUM*SMALL                                                   00004480
      B1 = B1*SMALL                                                     00004490
      B2 = B2*SMALL                                                     00004500
 445  IF ( STOP .LT. ACC*DABS(SUM) ) GO TO 450                          00004510
 449  CONTINUE                                                          00004520
      IRET = 8                                                          00004530
      GO TO 990                                                         00004540
C                                                                       00004550
 450  SUM = ( C*DEXP(X)*SUM ) * T                                       00004560
C                                                                       00004570
C     DID IT UNDERFLOW                                                  00004580
C                                                                       00004590
      IF ( SUM .EQ. 0 ) GO TO 900                                       00004600
C                                                                       00004610
C     WE NOW HAVE  F (=SUM),  R,  AND P  ( P ONLY IF RHO > .005 )       00004620
C     USE THE WRONSKIAN AS THE 4TH CONDITION                            00004630
C                                                                       00004640
      W = SUM/F                                                         00004650
      F = SUM                                                           00004660
      IF ( IGOTO .EQ. 5 ) GO TO 850                                     00004670
      X = (R-P)*F                                                       00004680
      IF ( DABS(X) .GT. PRERT3 ) GO TO 480                              00004690
C                                                                       00004700
C     HERE F**3 AND F**4 TERMS ARE LESS THAN MACHINE PRECISION          00004710
C                                                                       00004720
      G = 1/X                                                           00004730
      GP = P*G                                                          00004740
      GO TO 800                                                         00004750
C                                                                       00004760
C     HERE WE MUST INCLUDE F**3, F**4; WE MUST ALSO WORRY ABOUT         00004770
C     WHICH SIGN OF THE ROOT TO USE.                                    00004780
C     THE POSITIVE ROOT APPLIES FOR G > F; ELSE THE NEGATIVE ROOT       00004790
C     WE USE Q IN DETERMINING WHICH IS CORRECT                          00004800
C                                                                       00004810
 480  B1 = .5/X                                                         00004820
      B2 = B1 * DSQRT(1-4*(X*F)**2)                                     00004830
      G = B1 + B2                                                       00004840
C     G > F IN ALL OF REGION 2 FOR ETA > 0                              00004850
      IF ( ETA .GE. 0 ) GO TO 490                                       00004860
      SUM = 1/Q - F**2                                                  00004870
      GP = B1 - B2                                                      00004880
      IF ( DABS(G**2-SUM) .GT. DABS(GP**2-SUM) ) G = GP                 00004890
 490  GP = P*G - X*F/G                                                  00004900
      GO TO 800                                                         00004910
C                                                                       00004920
C     ETA > 15  AND  ETA < RHO < 2*ETA                                  00004930
C                                                                       00004940
C     WE FIND G AND G' FOR LMIN, RHO=2*ETA USING THE ABOVE METHOD       00004950
C     CONSISTING OF R, P+IQ, AND W.                                     00004960
C                                                                       00004970
 500  RHOUSE = ETA+ETA                                                  00004980
      PL = LMIN+1                                                       00004990
      IGOTO = 4                                                         00005000
      GO TO 105                                                         00005010
C                                                                       00005020
C     NOW WE HAVE  G, G'  AT THE TURNING POINT, GO IN USING TAYLOR      00005030
C                                                                       00005040
C                                                                       00005050
 600  DEL = RHOUSE - RHO                                                00005060
      B1 = G                                                            00005070
      B2 = -DEL*GP                                                      00005080
      B3 = 0                                                            00005090
      G = B1+B2                                                         00005100
      ACCR = ACC/2                                                      00005110
      DELINV = -1/DEL                                                   00005120
      DFACTR = 3*DELINV                                                 00005130
      X = DEL/RHOUSE                                                    00005140
      AI = X+X                                                          00005150
      DI = AI+AI                                                        00005160
      AR = 6                                                            00005170
      DR = 6                                                            00005180
      DO 639 I = 1, 10000                                               00005190
      S = ( AI*B3 + (X*DEL**2)*B1 ) / AR                                00005200
      AR = AR + DR                                                      00005210
      DR = DR + 2                                                       00005220
      AI = AI + DI                                                      00005230
      DI = DI + 2*X                                                     00005240
      G = G + S                                                         00005250
      GP = GP + DFACTR*S                                                00005260
      IF ( G .GE. VRYBIG ) GO TO 900                                    00005270
      DFACTR = DFACTR + DELINV                                          00005280
      B1 = B2                                                           00005290
      B2 = B3                                                           00005300
      B3 = S                                                            00005310
      IF ( S .LT. ACCR*G ) GO TO 650                                    00005320
 639  CONTINUE                                                          00005330
      IRET = 9                                                          00005340
      GO TO 990                                                         00005350
C                                                                       00005360
C     HERE WE HAVE  R = F'/F,  G,  G'                                   00005370
C     USE WRONSKIAN AS THE 4TH CONDITION                                00005380
C                                                                       00005390
 650  F = FC(LMIN1)                                                     00005400
      R = FCP(LMIN1)/F                                                  00005410
      SUM = 1/(R*G-GP)                                                  00005420
      W = SUM/F                                                         00005430
      F = SUM                                                           00005440
C                                                                       00005450
C     WE NOW HAVE  F, R = F'/F, G, G'  AT LMIN                          00005460
C                                                                       00005470
C     UPWARD RECURSION FROM GC(LMIN) AND GCP(LMIN),STORED VALUES ARE RHO00005480
C     RENORMALISE FC,FCP FOR EACH L-VALUE                               00005490
C                                                                       00005500
 800  GC (LMIN1) = G                                                    00005510
      GCP(LMIN1) = GP                                                   00005520
      FC(LMIN1) = F                                                     00005530
      FCP(LMIN1) = R*F                                                  00005540
      IRET = 0                                                          00005550
      IF(LMAX.EQ.LMIN) RETURN                                           00005560
      DO 829 L = LMIN1,LMAX                                             00005570
      T = GC(L+1)                                                       00005580
      GC (L+1) = (GC(L)*GC (L+1) - GCP(L))/GCP(L+1)                     00005590
      GCP(L+1) = GC(L)*GCP(L+1) - GC(L+1)*T                             00005600
      FC (L+1) = W*FC (L+1)                                             00005610
 829  FCP(L+1) = W*FCP(L+1)                                             00005620
 840  IF ( DABS(FC(LMAX+1))+DABS(FCP(LMAX+1)) .EQ. 0 ) IRET = IRET+1    00005630
      RETURN                                                            00005640
C                                                                       00005650
C     RHO < .005;  WE CANNOT FIND P OR Q AND SO RETURN ONLY F, F'.      00005660
C                                                                       00005670
 850  FC(LMIN1) = F                                                     00005680
      FCP(LMIN1) = R*F                                                  00005690
      IRET = 2                                                          00005700
      IF ( LMAX .EQ. LMIN ) RETURN                                      00005710
      DO 859 L = LMIN1, LMAX                                            00005720
      FC(L+1) = W*FC(L+1)                                               00005730
      FCP(L+1) = W*FCP(L+1)                                             00005740
 859  CONTINUE                                                          00005750
      GO TO 840                                                         00005760
C                                                                       00005770
C     F AND G ARE OUT OF THE MACHINE EXPONENT RANGE FOR LMIN.           00005780
C     IT WILL BE EVEN WORSE FOR  L > LMIN  SO GIVE UP AND RETURN        00005790
C                                                                       00005800
 900  IRET = 4                                                          00005810
      GO TO 990                                                         00005820
C                                                                       00005830
C     P + IQ FAILED TO CONVERGE                                         00005840
 920  IRET = 7                                                          00005850
C                                                                       00005860
C     ERROR DETECTED, PRINT THE INPUT AND SOME STUFF                    00005870
C                                                                       00005880
990      continue
      RETURN                                                            00005940
C                                                                       00005950
      END                                                               00005960
      FUNCTION SIXJ (A, B, C, X, Y, ZZ)                                  00000100
CDC      REAL FUNCTION SIXJ ( A, B, ARG3, ARG4, ARG5, Z )               00000200
C                                                                       00000300
C                                                                       00000400
C     COMPUTES 6-J AND RACAH COEFICENTS.                                00000500
C                                                                       00000600
C       THIS ROUTINE COMPUTES THE 6-J OR RACAH COEFICENT FOR            00000700
C     A GIVEN SET OF J'S.  THE INPUT CONSISTS OF 6 INTEGERS THAT        00000800
C     ARE TWICE THE VALUES OF THE J'S.  IN THIS MANNER                  00000900
C     HALF-INTEGER ANGULAR MOMENTA CAN BE ACCOMODATED.                  00001000
C                                                                       00001100
C       FOR THE FOLLOWING DISCUSSION CONSIDER THE 6-J COEFFICENT:       00001200
C                                                                       00001300
C           ( J1  J2  J3 )                                              00001400
C           ( J4  J5  J6 )                                              00001500
C                                                                       00001600
C                                                                       00001700
C       A FACTORIZATION OF THE RACAH SUM SUGGESTED BY                   00001800
C         J. G. WILLS, COMPUTER PHYSICS COMM. 2, 381 (1971)             00001900
C     IS USED.  THE REQUIRED FACTORIALS THAT MULTIPLY THE SUM ARE       00002000
C     FOUND IN ONE OF 3 WAYS:                                           00002100
C       FOR J1+J2+J4+J5 < 96,        A  TABLE OF SQRT(FACTORIAL) IS USED00002200
C       FOR 95 < J1+J2+J4+J5 < 1000, A TABLE OF LOG(FACTORIAL) IS USED; 00002300
C       FOR 1000 <= J1+J2+J4+J5,     DLGAMA IS USED.                    00002400
C                                                                       00002500
C       THESE CHOICES REPRESENT REASONABLE TRADE OFFS BETWEEN SIZE AND  00002600
C     SPEED CONSIDERATIONS.  BY RECOMPILING THE BLOCK DATA PROGRAMS     00002700
C     FOR THE /FACTRL/ OR /LOGFAC/ COMMON BLOCKS, ONE CAN CHANGE        00002800
C     THE ABOVE BOUNDRIES TO EITHER REDUCE CORE (IF ONE IS NOT          00002900
C     INTERESTED IN LARGE J'S) OR INCREASE THE SPEED FOR LARGE J'S.     00003000
C                                                                       00003100
C       THE MAGNITUDE OF THE J'S AND M'S THAT CAN BE ACCOMODATED BY     00003200
C     THIS ROUTINE IS DETERMINED BY CANCELLATIONS IN THE                00003300
C     RACAH SUM.  AS CAN BE SEEN IN THE TABLE BELOW, IF ALL SIX J'S     00003400
C     GET LARGE, SEVERE CANCELLATION SETS IN FOR EACH J OF THE ORDER OF 00003500
C     100.  HOWEVER IF THREE OF THE J'S ARE HELD SMALL, THE PRECISION   00003600
C     REMAINS GOOD FOR VERY LARGE VALUES OF THE OTHER THREE J'S AND     00003700
C     THERE ARE NO OVERFLOW PROBLEMS (THE LARGEST VALUES TRIED WERE     00003800
C     500 AS SHOWN IN THE TABLE).                                       00003900
C                                                                       00004000
C                                                                       00004100
C       EACH LINE IN THE FOLLOWING TABLE IS THE RESULT OF SEVERAL       00004200
C     THOUSAND EVALUATIONS OF 6-J COEFFICENTS FOR J'S                   00004300

C     IN THE RANGES INDICATED.  PRECISIONS WERE DETERMINED BY           00004400
C     COMPARING THE RESULTS TO A QUADRUPOLE PRECISION 6-J               00004500
C     ROUTINE.  IT IS NOT SIGNIFICANT THAT J3, J4 AND J5 WERE HELD      00004600
C     SMALL; ANY THREE OF THE J'S MAY BE KEPT SMALL WITH THE SAME TIMING00004700
C     AND RELATIVE ERROR RESULTS.  THE TIMING INFORMATION               00004800
C     IS IN MICROSECONDS.                                               00004900
C                                                                       00005000
C                                                                       00005100
C       J1, J2, J6    J3, J4, J5        TIMES      RMS ABSOLUTE         00005200
C         RANGE         RANGE        /195   /75       ERROR             00005300
C         0   5         0   5         49    400       4E-17             00005400
C         5  10         0  10         64    615       7E-18             00005500
C        10  15         0  25         59    550       5E-18             00005600
C        15  30         0  47         73    650       2E-16             00005700
C        33  42         0  75         95    790       2E-15             00005800
C        50  52         0 100        112   1020       3E-13             00005900
C          75           0 150        130   1140       3E-10             00006000
C         100           0 200        146   1315       5E-6              00006100
C         150            150         575   5520       3E+2              00006200
C                                                                       00006300
C          50           0  10         68    485       9E-16             00006400
C         100           0  10         68    485       7E-16             00006500
C         400           0  10         70    478       6E-15             00006600
C         500           0  10        670   5080       9E-15             00006700
C                                                                       00006800
C                                                                       00006900
C     NO EVALUATION OF 6-J COEFICENTS THAT WERE TRIVIALLY               00007000
C     ZERO WERE INCLUDED IN THE ABOVE RESULTS.  THE LONG EXECUTION      00007100
C     TIMES IN THE LAST LINE OF THE TABLE ARE DUE TO THE USE OF THE     00007200
C     DLGAMA FUNCTION.  THEY MAY BE MADE COMPARABLE TO THE TIMES IN     00007300
C     THE PREVIOUS LINE BY INCREASING THE SIZE OF THE  /LOGFAC/         00007400
C     TABLE.                                                            00007500
C                                                                       00007600
C                                                                       00007700
C     THERE ARE THREE ENTRIES:                                          00007800
C                                                                       00007900
C     SIXJ (2*J1, 2*J2, 2*J3, 2*J4, 2*J5, 2*J6)   RETURNS THE           00008000
C                  DOUBLE PRECISION 6-J COEFICENT.                      00008100
C                        ( J1  J2  J3 )                                 00008200
C                        ( J4  J5  J6 )                                 00008300
C                                                                       00008400
C     RACAH (2*J1, 2*J2, 2*J5, 2*J4, 2*J3, 2*J6)   RETURNS THE          00008500
C                  DOUBLE PRECISION RACAH SYMBOL                        00008600
C                     W (J1, J2, J5, J4; J3, J6)                        00008700
C                                                                       00008800
C     COF6J (2*J1, 2*J2, 2*J3, 2*J4, 2*J5, 2*J6)    RETURNS A           00008900
C                  PROPERLY ROUNDED SINGLE PRECISION 6-J COEFICENT.     00009000
C                                                                       00009100
C     IN ALL CASES THE TRIANGLE INEQUALITIES ARE CHECKED.  IF THEY ARE  00009200
C     VIOLATED, 0 IS RETURNED.  HOWEVER NO CHECK IS MADE TO VERIFY      00009300
C     THAT  J1+J2+J3,  J4+J2+J6,  J1+J5+J6,  AND  J4+J5+J3  ARE         00009400
C     INTEGERS.  THE USER MUST INSURE THAT  2*J1 + 2*J2 + 2*J3,  ETC,   00009500
C     ARE EVEN.  IF THIS IS NOT SO MEANINGLESS RESULTS WILL BE RETURNED.00009600
C                                                                       00009700
C       CORE REQUIREMENTS (BASE 10):                                    00009800
C                                                                       00009900
C          SIXJ      4000 BYTES                                         00010000
C         /FACTRL/    798     MAY BE MADE SMALLER                       00010100
C         /LOGFAC/   8028     MAY BE MADE SMALLER OR LARGER             00010200
C          DLGAMA    1008                                               00010300
C          TOTAL   13,834                                               00010400
C                                                                       00010500
C     IN ADDITION "DEXP" AND THE FORTRAN I/O PACKAGE WILL BE USED.      00010600
C     SINCE THESE ARE LIKELY TO BE USED BY OTHER PARTS OF YOUR          00010700
C     PROGRAM, THEY ARE NOT INCLUDED IN THE ABOVE (THE FORTRAN I/O      00010800
C     PACKAGE IS CALLED ONLY BY DLGAMA AND DEXP).  THE  /FACTRL/        00010900
C     AND  /LOGFAC/  COMMON BLOCKS AND THE "DLGAMA" FUNCTION WILL BE    00011000
C     SHARED WITH "CLEBSH" AND "THREEJ" SHOULD YOU ALSO USE THEM.       00011100
C                                                                       00011200
C                                                                       00011300
C     S. PIEPER                                                         00011400
C                                                                       00011500
C      4/1/71  MODIVIED FROM COFJU TO GIVE 6-J                          00011600
C     8/21/72  ARGONNE DOUBLE PRECISION VERSION                         00011700
C     9/5/72  INTEGER VERSION -- S. PIEPER                              00011800
C     9/6/72  USE SQRT OF BINOMIAL COFS AND FACTORIALS                  00011900
C     9/8/72  SINGLE AND DOUBLE PRECISION ENTRY POINTS                  00012000
C     9/23/73  SIXJ AND RACAH ENTRIES.                                  00012100
C     7/6/74 - MAJOR REWRITE - TOTALLY NEW METHOD                       00012200
C     1/16/80 - STANDARD MORTRAN - RPG                                  00012300
C     3/13/80 - BACK TO TBITS, ETC.,  - S.P.                            00012400
C     5/29/80 - PUT BACK TRIANGLE TEST FOR CDC - S.P.                   00012500
C     7/9/80 - CRAY TRIANGLE TEST - S.P.                                00012600
C     7/31/80 - USE LOGFAC TABLE IN THE ALLOCATOR - S.P.                00012700
C     8/25/80 - COF6J REMOVED FOR PTOLEMY - S.P.                        00012800
C                                                                       00012900
      IMPLICIT INTEGER (A-C, I-Z), REAL*8 (D-H)                         00013000
C                                                                       00013100
CNV      REAL ALGAMA                                                    00013200
      REAL*8 SUM, ANSWER                                                00013300
      REAL*8 SIXJ, RACAH                                                00013400
CXXCND      REAL*4  COF6J                                               00013500
      LOGICAL DBLSW, RACHSW                                             00013600
C                                                                       00013700
      COMMON /FACTRL/ MAXFAC, ISPACE,FACTBL(57)                         00013800
C                                                                       00013900
      REAL*8 LF                                                         00014000
      COMMON /LOGFAC/ MAXLF, LFBIAS, LF(1)                              00014100
CRA      COMMON /LOGFAC/ MAXLF, LFBIAS, LFTBL(1)                        00014200
	REAL*8 ALLOC
	COMMON/ALLOCS/ALLOC(1)
	CHARACTER*8 LOGNAM / 'LOGFACS ' /
	INTEGER*4 Z
	COMMON /LOCKT/ Z(100)
C                                                                       00014300
CRA      POINTER ( LFBASE, LF(1) )                                      00014500
C                                                                       00014600
C                                                                       00014700
C     THE FOLLOWING FAKE EXTERNALS DRAG IN THE REQUIRED BLOCK DATAS     00014800
C     FOR THE ABOVE COMMON BLOCKS                                       00014900
C                                                                       00015000
      EXTERNAL FACDUM, LOGDUM                                           00015100
C                                                                       00015200
CCCCC**********                                                         00015300
CCCC      COMMON /NUMBLK/  NUMIT, IBIG, INOTSO, ILF, E, SUM             00015400
CCCCC**********                                                         00015500
C                                                                       00015600
C     WE USE THE FOLLOWING FORTRAN H XL FUNCTIONS:                      00015700
C     TBIT - TESTS TO SEE IF A BIT IS ON.  TBIT(I, 31) IS TRUE          00015800
C            IF  I  IS ODD.                                             00015900
C     SHFTR - SHIFTS RIGHT AND IS EQUIVALENT TO A DIVIDE BY 2**N.       00016000
C                                                                       00016100
CIB      LOGICAL TBIT                                                   00016200
CIBC                                                                    00016300
CIB      REAL*8 SHIFT / Z 4F000000 00000000 /                           00016400
CIB      REAL*8  ADDER / .5D0 /                                         00016500
C                                                                       00016600
C                                                                       00016700
C                                                                       00016800
CXXC     THE FOLLOWING ARE FOR ROUNDING TO SINGLE PRECISION             00016900
CXXC                                                                    00017000
CXXCIB      REAL*8  ROUND /Z 00 000000 80000000 /                       00017100
CXXCIB      LOGICAL*1 ANSEXP, RNDEXP                                    00017200
CXXCIB      EQUIVALENCE (ANSEXP, ANSWER), (RNDEXP, ROUND)               00017300
C                                                                       00017400
      DIMENSION DRAY(8)                                                 00017600
      EQUIVALENCE (DA1, DRAY(1)), (DA2, DRAY(2)),  (DA3, DRAY(3)), (DA4,00017700
     * DRAY(4)),  (DC, DRAY(5)), (DZ, DRAY(7))                          00017700
C                                                                       00018000
C                                                                       00018100
C    SIXJ IS THE DOUBLE PRECISON ENTRY POINT                            00018200
C                                                                       00018300
CXXCNDCNR      DBLSW = .TRUE.                                           00018400
      RACHSW = .FALSE.                                                  00018500
CDC      C = ARG3                                                       00018600
CDC      X = ARG4                                                       00018700
CDC      Y = ARG5                                                       00018800
      GO TO 100                                                         00018900
CXXCNDCNRC                                                              00019000
CXXCNDCNRC    COF6J IS THE SINGLE PRECISION ENTRY POINT                 00019100
CXXCNDCNRC                                                              00019200
CXXCNDCNR      ENTRY COF6J (A, B, C, X, Y, Z)                           00019300
CXXCNDCNR      DBLSW = .FALSE.                                          00019400
CXXCNDCNR      RACHSW = .FALSE.                                         00019500
CXXCNDCNR      GO TO 100                                                00019600
C                                                                       00019700
C     RACAH ENTRY IS DOUBLE PRECISION                                   00019800
C                                                                       00019900
      ENTRY RACAH (A, B, Y, X, C, ZZ)                                    00020000
C                                                                       00020100
      RACHSW = .TRUE.                                                   00020200
CXXCNDCNR      DBLSW = .TRUE.                                           00020300
CDC      Y = ARG3                                                       00020400
CDC      X = ARG4                                                       00020500
CDC      C = ARG5                                                       00020600
C                                                                       00020700
C                                                                       00020800
C                                                                       00020900
C                                                                       00021000
 100  SUM=1                                                             00021100
C                                                                       00021200
CRA      LFBASE = LOCF(LFTBL) + LFBIAS                                  00021300
C                                                                       00021400
CCCCC****************                                                   00021500
CCCC      NUMIT = 0                                                     00021600
CCCCC****************                                                   00021700
C                                                                       00021800
C                                                                       00021900
	IF((X+Y-C).LT.0)GO TO 900
	IF((A+Y-ZZ).LT.0)GO TO 900
	IF((B+X-ZZ).LT.0)GO TO 900
	IF((A+B-C).LT.0)GO TO 900
      A2 = JISHFT(X+Y-C, -1)                                              000220
      A3 = JISHFT(A+Y-ZZ, -1)                                              000221
      A4 = JISHFT(B+X-ZZ, -1)                                              000222
      A1 = JISHFT(A+B-C, -1)                                              000223
C     B1 = J3 +J6 - J1 - J4                                             00022400
C     B2 = J3 + J6 - J2 - J5                                            00022500
C     R  = J1 + J2 + J4 + J5 + 1                                        00022600
      B1 = Y - A2 - A3                                                  00022700
      B2 = A3 - A1 + ZZ - Y                                              00022800
C                                                                       00022900
C     THE SHFTR IS A LOGICAL SHIFT AND HENCE CONVERTS NEGATIVE          00023000
C     QUANTITIES INTO VERY LARGE POSSITIVE QUANTITIES.  IF ANY OF       00023100
C     THE AI'S ARE NEGATIVE, THEN THE CORRESPONDING TRIANGLE            00023200
C     IS VIOLATED.  WE TEST FOR THAT CASE HERE.                         00023300
CCDC     ON CDC ORING A NEGATIVE NUMBER WITH ANYTHING RESULTS IN        00023400
CCDC     A NEGATIVE NUMBER SO WE TEST THE RESULT FOR < 0.               00023500
CCRA      ON CRAY WE TEST FOR A HUGE INTEGER                            00023600
C                                                                       00023700
      TEST = IOR(A2,IOR(A3,IOR(A4,A1)))                                 00023800
C                                                                       00023900
      R = A1 + A2 + C + 1                                               00024000
C                                                                       00024100
CIB      IF (TBIT(TEST, 1))  GO TO 900                                  00024200
CDC      IF ( AND( TEST, -1 )  .LT. 0 )  GO TO 900                      00024300
CRA      IF ( TEST .GT. 1 000 000 000 )  GO TO 900                      00024400
C                                                                       00024500
      NMAX = MIN0 (A1, A2, A3, A4)                                      00024600
      NMIN = MAX0 (0, -B1, -B2)                                         00024700
C                                                                       00024800
C                                                                       00024900
CCCCC***************                                                    00025000
CCCC      NUMIT = NMAX - NMIN + 1                                       00025100
CCCCC***************                                                    00025200
C                                                                       00025300
C     CONVERT VARIOUS INTEGERS TO DOUBLE PRECISON FOR NOW AND           00025400
C     POSSIBLY LATER.                                                   00025500
C                                                                       00025600
	ILOGFC=NAMLOC(LOGNAM)
	LLOGFC=Z(ILOGFC)-1
      DA1 = A1                                                          00025700
      DA2 = A2                                                          00025800
      DA3 = A3                                                          00025900
      DA4 = A4                                                          00026000
      DY = Y                                                            00026100
      DZ = ZZ                                                            00026200
      DC = C                                                            00026300
      DB1 = DY - DA2 - DA3                                              00026400
      DB2 = DA3 - DA1 + DZ - DY                                         00026500
      DR = DA1 + DA2 + DC + 1                                           00026600
      DN = NMIN                                                         00026700
      DX = NMAX                                                         00026800
C                                                                       00026900
C     THIS ONE TEST CHECKS ALL 12 INEQUALITIES IMPLIED BY THE           00027000
C     4 TRIANGLE CONDITIONS.                                            00027100
C                                                                       00027200
      IF (NMAX - NMIN) 900, 300, 200                                    00027300
C                                                                       00027400
 200  D1 = DA1 - (DX-1)                                                 00027500
      D2 = DA2 - (DX-1)                                                 00027600
      D3 = DA3 - (DX-1)                                                 00027700
      D4 = DA4 - (DX-1)                                                 00027800
      D5 = -DX                                                          00027900
      D6 = -DB1 - DX                                                    00028000
      D7 = -DB2 - DX                                                    00028100
      D8 = DR - (DX-1)                                                  00028200
C                                                                       00028300
      E1 = D1*D2*D3*D4                                                  00028400
      E5 = D5*D6*D7*D8                                                  00028500
C                                                                       00028600
      NMAXM2 = NMAX-2                                                   00028700
      SUM = 1 + E1/E5                                                   00028800
      IF (NMAXM2 .LT. NMIN) GO TO 300                                   00028900
C                                                                       00029000
      F1 = (D1+1)*(D2+1)*(D3+1)*(D4+1) - E1                             00029100
      F5 = (D5+1)*(D6+1)*(D7+1)*(D8+1) - E5                             00029200
C                                                                       00029300
      G1 = (D1+2)*(D2+2)*(D3+2)*(D4+2) - F1 - F1 - E1                   00029400
      G5 = (D5+2)*(D6+2)*(D7+2)*(D8+2) - F5 - F5 - E5                   00029500
C                                                                       00029600
      H1 = (D1+3)*(D2+3)*(D3+3)*(D4+3) - 3*(G1+F1) - E1                 00029700
      H5 = (D5+3)*(D6+3)*(D7+3)*(D8+3) - 3*(G5+F5) - E5                 00029800
C                                                                       00029900
C     COMPUTE THE SUM IN THE RACAH FORMULA                              00030000
C                                                                       00030100
      DO 259 N = NMIN, NMAXM2                                           00030200
      E1 = E1 + F1                                                      00030300
      E5 = E5 + F5                                                      00030400
      F1 = F1 + G1                                                      00030500
      F5 = F5 + G5                                                      00030600
      G1 = G1 + H1                                                      00030700
      G5 = G5 + H5                                                      00030800
      H1 = H1 + 24                                                      00030900
      H5 = H5 + 24                                                      00031000
      SUM = 1 + SUM*(E1/E5)                                             00031100
 259  CONTINUE                                                          00031200
C                                                                       00031300
C     NOW WE NEED THE OUTSIDE FACTORIALS                                00031400
C     WE USE ONE OF 3 METHODS DEPENDING ON THE SIZE OF J1+J2+J4+J5      00031500
C        1)  A TABLE OF SQRT(FACTORIALS)                                00031600
C        2)  A TABLE OF LN(FACTORIALS)                                  00031700
C        3)  LN(GAMMA) FUNCTION                                         00031800
C                                                                       00031900
 300  IF (R .GT. MAXFAC) GO TO 400                                      00032000
C                                                                       00032100
C     FIRST THE STUFF WE FACTORED OUT OF THE SUM                        00032200
C                                                                       00032300
      E = FACTBL(1+R-NMIN) / (FACTBL(1+NMIN)  * FACTBL(1+B1+NMIN) * FACT00032400
     *BL(1+B2+NMIN) * FACTBL(1+A1-NMIN)  * FACTBL(1+A2-NMIN) * FACTBL(1+00032400
     *A3-NMIN) * FACTBL(1+A4-NMIN) )                                    00032400
C                                                                       00032700
C     SUM IS AN INTEGER NOW                                             00032800
C                                                                       00032900
CIB      SUM = (E*SUM*E + DSIGN(ADDER, SUM)) + SHIFT                    00033000
      SUM = DINT ( E*SUM*E + DSIGN(.50D0, SUM) )                        00033100
C                                                                       00033200
C     AND NOW THE OUTSIDE FACTORS                                       00033300
C                                                                       00033400
C                                                                       00033500
      DENOM = (FACTBL(2+A1+C)/(FACTBL(1+A1)*FACTBL(1+A4+B1)  *FACTBL(1+A00033600
     *3+B2)))  * (FACTBL(2+A2+C)/(FACTBL(1+A2)*FACTBL(1+A3+B1)  *FACTBL(00033600
     *1+A4+B2))) *(FACTBL(2+A3+ZZ )/(FACTBL(1+A3)*FACTBL(1+A2+B1)  *FACT00033600
     *BL(1+A1+B2)))  * (FACTBL(2+A4+ZZ)/(FACTBL(1+A4)*FACTBL(1+A1+B1) *F00033600
     *ACTBL(1+A2+B2)))                                                  00033600
      ANSWER = SUM/(DENOM)                                              00034400
      GO TO 600                                                         00034500
C                                                                       00034600
C     TRY FOR LOG FACTORIAL TABLE                                       00034700
C                                                                       00034800
 400  IF (R .GT. MAXLF) GO TO 500                                       00034900
C                                                                       00035000
CCCCC**********                                                         00035100
CCCC      IBIG = IBIG + 1                                               00035200
CCCC      ILF = ILF + 1                                                 00035300
CCCCC************                                                       00035400
C                                                                       00035500
	E=ALLOC(LLOGFC+1+R-NMIN)-ALLOC(LLOGFC+1+NMIN)-
	1ALLOC(LLOGFC+1+B1+NMIN)-ALLOC(LLOGFC+1+B2+NMIN)-
	2ALLOC(LLOGFC+1+A1-NMIN)-ALLOC(LLOGFC+1+A2-NMIN)-
	3ALLOC(LLOGFC+1+A3-NMIN)-ALLOC(LLOGFC+1+A4-NMIN)
	E2=.5*(ALLOC(LLOGFC+1+A1)+ALLOC(LLOGFC+1+A4+B1)+
	1ALLOC(LLOGFC+1+A3+B2)-ALLOC(LLOGFC+2+A1+C)+
	2ALLOC(LLOGFC+1+A2)+ALLOC(LLOGFC+1+A3+B1)
	3+ALLOC(LLOGFC+1+A4+B2)-ALLOC(LLOGFC+2+A2+C)+
	4ALLOC(LLOGFC+1+A3)+ALLOC(LLOGFC+1+A2+B1)+
	5ALLOC(LLOGFC+1+A1+B2)-ALLOC(LLOGFC+2+A3+ZZ)+
	6ALLOC(LLOGFC+1+A4)+ALLOC(LLOGFC+1+A1+B1)+
	7ALLOC(LLOGFC+1+A2+B2)-ALLOC(LLOGFC+2+A4+ZZ))
C     AT THIS POINT IT WOULD BE POSSIBLE TO TRY TO INTERGERIZE          00036500
C     SUM*EXP(E).  HOWEVER, TESTS SHOW THAT IN PRACTICE THIS WOULD      00036600
C     BE USEFUL (RESULT < 10**18) LESS THAN 1% OF THE TIME SO           00036700
C     WE DON'T BOTHER WITH IT.                                          00036800
C                                                                       00036900
 450  ANSWER = DEXP(E+E2) * SUM                                         00037000
      GO TO 600                                                         00037100
C                                                                       00037200
C     HERE THE J'S ARE VERY LARGE - WE USE LN(GAMMA)                    00037300
C                                                                       00037400
 500  DRAY(6) = DC                                                      00037500
      DRAY(8) = DZ                                                      00037600
      E2 = 0                                                            00037700
      E = DLGAMA(1-DN+DR) - DLGAMA(1+DN)  - DLGAMA(1+DN+DB1) - DLGAMA(1+00037800
     *DN+DB2)                                                           00037800
      DO 549 I = 1, 4                                                   00038000
      E2 = E2 + DLGAMA(1+DRAY(I)) + DLGAMA(1+DB1+DRAY(I))  + DLGAMA(1+DB00038100
     *2+DRAY(I)) - DLGAMA(2+DRAY(I)+DRAY(4+I))                          00038100
      E = E - DLGAMA(1-DN+DRAY(I))                                      00038300
 549  CONTINUE                                                          00038400
      E = E + .5*E2                                                     00038500
C                                                                       00038600
CCCC      IBIG = IBIG + 1                                               00038700
CCCCC***************                                                    00038800
CCCCC                                                                   00038900
      ANSWER = DEXP(E) * SUM                                            00039000
C                                                                       00039100
C     NOW GET THE SIGN CORRECT                                          00039200
C                                                                       00039300
 600  TEMP = NMIN + R                                                   00039400
      IF (MOD(TEMP,2).EQ.0) ANSWER = -ANSWER                            00039500
C                                                                       00039600
C                                                                       00039700
C     CALCULATION OF THE 6-J COMPLETED; WAS RACAH REQUESTED             00039800
C                                                                       00039900
 700  IF (.NOT. RACHSW) GO TO 800                                       00040000
C                                                                       00040100
C     YES, CONVERT THE SIGN                                             00040200
      IF (MOD(R,2).EQ.0) ANSWER = -ANSWER                               00040300
C                                                                       00040400
C     WHAT PRECISION WAS REQUESTED FOR THE RESULT                       00040500
C                                                                       00040600
 800  CONTINUE                                                          00040700
CXXCNDCNR      IF (DBLSW)  GO TO 850                                    00040800
CXXCNDCNRC                                                              00040900
CXXCNDCNRC    ROUND THE ANSWER TO SINGLE PRECISION                      00041000
CXXCNDCNRC                                                              00041100
CXXCNDCNR      RNDEXP = ANSEXP                                          00041200
CXXCNDCNR      IF (ANSWER .EQ. 0)  GO TO 830                            00041300
CXXCNDCNR      COF6J = ANSWER + ROUND                                   00041400
CXXCNDCNR      RETURN                                                   00041500
CXXCNDCNRC    ROUNDING ZERO WILL RESULT IN UNDERFLOW                    00041600
CXXCNDCNR 830  COF6J = ANSWER                                           00041700
CXXCNDCNR      RETURN                                                   00041800
C                                                                       00041900
C     DOUBLE PRECISION ENTRY POINT WAS USED                             00042000
C                                                                       00042100
 850  SIXJ = ANSWER                                                     00042200
      RETURN                                                            00042300
C                                                                       00042400
C     THE TRIANGLE RULES ARE NOT SATISFIED                              00042500
C                                                                       00042600
 900  SIXJ = 0                                                          00042700
      RETURN                                                            00042800
C                                                                       00042900
C                                                                       00043100
      END                                                               00043200
      SUBROUTINE SPLNCB (N, X, Y, B, C, D)                              00000010
C                                                                       00000020
C     INTERPOLATION BY CUBIC SPLINES.                                   00000030
C                                                                       00000040
C       THIS ROUTINE IS A DRASTIC SIMPLIFICATION OF THE AMD ROUTINE     00000050
C     "SMOOTH".  IT WILL FIND THE COEFICENTS OF INTERPOLATING CUBICS    00000060
C     FOR A FUNCTION TABULATED ON A MONOTONIC GRID OF X'S.  THE         00000070
C     ROUTINE "INTRPC" MAY THEN BE USED TO FIND INTERPOLATED VALUES     00000080
C     OF THE FUNCTION USING THESE COEFFICENTS.  THE INTERPOLATION       00000090
C     SCHEME IS THE METHOD OF CUBIC SPLINES.  THE COEFICENTS FOUND BY   00000100
C     THIS ROUTINE ARE THE SAME AS THOSE FOUND BY THE AMD "SMOOTH"      00000110
C     ROUTINE FOR THE CASE OF   S = 0   IF THE SPLNCB ENTRY IS USED.    00000120
C                                                                       00000130
C      THE INPUT IS                                                     00000140
C                                                                       00000150
C     N - THE NUMBER OF INPUT POINTS.  N MUST BE >= 2.                  00000160
C         IF N <= 1, NOTHING IS DONE (IMMEADIATE RETURN).               00000170
C     X - THE ARRAY OF INPUT INDEPENDENT VARIABLES.  X                  00000180
C         MUST BE MONOTONIC (INCREASING OR DECREASING) WITH             00000190
C         NO EQUAL POINTS.  IT MUST HAVE "N" ELEMENTS.                  00000200
C     Y - THE ARRAY OF FUNCTION VALUES CORRESPONDING TO X.              00000210
C     B, C, D - THE OUTPUT ARRAYS WITH THE FITTING CUBICS.              00000220
C         THESE MUST EACH BE "N" ELEMENT ARRAYS.  THE LAST              00000230
C         ELEMENT OF EACH ARRAY WILL BE SET TO THE COEFICENTS           00000240
C         CORRESPONDING THE THE CUBIC FOR THE N-1 INTERVAL.             00000250
C                                                                       00000260
C     THE INTERPOLATING CUBICS ARE:                                     00000270
C                                                                       00000280
C      Y(X) = Y(I) + DEL*B(I) + DEL**2*C(I) +                           00000290
C             DEL**3*D(I)                                               00000300
C     WHERE                                                             00000310
C       DEL = X - X(I)                                                  00000320
C     AND                                                               00000330
C       X(I) <= X <= X(I+1),  1 <= I <= N-1                             00000340
C                                                                       00000350
C                                                                       00000360
C     AN ALTERNATIVE ENTRY IS:                                          00000370
C                                                                       00000380
C     CALL  SPLNCX (N, X, Y, B, C, D, IOPT, Y1P, YNP)                   00000390
C                                                                       00000400
C     WHERE THE FIRST 6 ARGUMENTS ARE AS ABOVE AND THE LAST 3 ARE       00000410
C                                                                       00000420
C     IOPT - CONTROLS OPTIONS CONCERNING THE BOUNDRY CONDITIONS:        00000430
C          = 1 MEANS THAT THE 2ND DERIVATIVES ARE ZERO AT X(1) AND      00000440
C              X(N).  THIS IS THE "NATURAL" CHOICE AND IS ALSO THE      00000450
C              RESULT OF CALLING SPLNCB.                                00000460
C          = 2 MEANS THAT THE 3RD DERIVATIES ARE ZERO AT X(1) AND       00000470
C              X(N).  THIS ALLOWS AN EXACT FIT TO PURE PARABOLAS        00000480
C              WHICH IS NOT POSSIBLE WITH OPTION 1.                     00000490
C          = 3 MEANS THAT THE 1ST DERIVATIVES ARE BEING SPECIFIED AT    00000500
C              X(1) AND X(N).                                           00000510
C     Y1P, YNP - FOR IOPT=3 THESE ARE THE REQUIRED 1ST DERIVATIVS AT    00000520
C              X(1) AND X(N).  THEY ARE NOT USED FOR IOPT = 1 OR 2.     00000530
C                                                                       00000540
C                                                                       00000550
C     7/4/73 - DRASTIC CUTS TO AMD "SMOOTH" ROUTINE. - S. PIEPER        00000560
C     9/5/73 - ADD SPLNCX ENTRY.                                        00000570
C     4/4/78 - PROTECT AGAINST UNDERFLOWS AT BOTH ENDS;                 00000580
C              REMOVE OPTIONAL ENTRY FOR PTOLEMY VERSION - S.P.         00000590
C     5/19/78 - FIX COUNTS IN ABOVE FOR ALL ZERO ARRAYS - S.P.          00000600
C     6/28/78 - DON=T USE UNDEFINED D(NM1) AFTER 300 LOOP               00000610
C                                                                       00000620
C                                                                       00000630
      IMPLICIT REAL*8 ( A-H, O-Z )                                      00000640
CDC      ALLOCLEVEL X,Y,B,C,D                                           00000650
      DIMENSION X(1),Y(1),B(1),C(1),D(1)                                00000660
C                                                                       00000670
C     NATURAL CUBIC BY DEFAULT                                          00000680
C                                                                       00000690
CXX 40   ASSIGN 300 TO IGO                                              00000700
CXX      IOPT = 1                                                       00000710
C                                                                       00000720
  50  IF (N .LE. 1) RETURN                                              00000730
C                                                                       00000740
      NM1=N-1                                                           00000750
C                                                                       00000760
C     LONG STRINGS OF CONSTANT VALUES CAN CAUSE UNDERFLOW PROBLEMS;     00000770
C     WE REMOVE SUCH STRINGS FROM THE START.                            00000780
C                                                                       00000790
      DO 74 I = 1, NM1                                                  00000800
      B(I) = 0                                                          00000810
      C(I) = 0                                                          00000820
      D(I) = 0                                                          00000830
      IF ( 1.D15*DABS(Y(I+1)-Y(I)) .GT. DABS(Y(I)) ) GO TO 80           00000840
 74   CONTINUE                                                          00000850
      I = NM1                                                           00000860
C                                                                       00000870
 80   IBASE = I-1                                                       00000880
      DO 84 I1 = I, NM1                                                 00000890
      J = N+I - I1                                                      00000900
      B(J) = 0                                                          00000910
      C(J) = 0                                                          00000920
      D(J) = 0                                                          00000930
      IF ( 1.D15*DABS(Y(J-1)-Y(J)) .GT. DABS(Y(J)) ) GO TO 90           00000940
 84   CONTINUE                                                          00000950
  90  NUSE = J - IBASE                                                  00000960
      IF ( NUSE .LE. 1 ) RETURN                                         00000970
      NM1 = NUSE-1                                                      00000980
C                                                                       00000990
      H=X(IBASE+2)-X(IBASE+1)                                           00001000
      F=(Y(IBASE+2)-Y(IBASE+1))/H                                       00001010
      IF ( NUSE .EQ. 2) GO TO 800                                       00001020
      DO 399 I=2,NM1                                                    00001030
      G=H                                                               00001040
      H=X(IBASE+I+1)-X(IBASE+I)                                         00001050
      E=F                                                               00001060
      F=(Y(IBASE+I+1)-Y(IBASE+I))/H                                     00001070
      GBY3 = G/3                                                        00001080
      D(IBASE+I-1) = GBY3*B(IBASE+I-1)                                  00001090
      EPSIM1 = G+H                                                      00001100
      RIM1B3 = F-E                                                      00001110
C                                                                       00001120
CXX      GO TO  IGO, (200, 220, 300)                                    00001130
CXXC                                                                    00001140
CXXC     SPECIAL OPTIONS CHOOSEN:  END POINTS MUST BE PATCHED UP        00001150
CXXC                                                                    00001160
CXXC                                                                    00001170
CXXC     2)  Y'''(X1) = Y'''(XN) = 0                                    00001180
CXXC                                                                    00001190
CXX 200  IF (I .EQ. 2)  EPSIM1 = EPSIM1 + .5*G                          00001200
CXX      IF (I .EQ. NM1)  EPSIM1 = EPSIM1 + .5*H                        00001210
CXX      GO TO 300                                                      00001220
CXXC                                                                    00001230
CXXC     3)  Y'(X1)  AND  Y'(XN)  FED IN                                00001240
CXXC                                                                    00001250
CXX 220  IF (I .NE. 2)  GO TO 230                                       00001260
CXX      EPSIM1 = .75*G + H                                             00001270
CXX      RIM1B3 = F - 1.5*E - .5*Y1P                                    00001280
CXX      GO TO 300                                                      00001290
CXX 230  IF (I .NE. NM1)  GO TO 300                                     00001300
CXX      EPSIM1 = G  +  .75*H                                           00001310
CXX      RIM1B3 = 1.5*F - E - .5*YNP                                    00001320
C                                                                       00001330
 300  B(IBASE+I) = 1/((2.D0/3.D0)*EPSIM1 - GBY3*D(IBASE+I-1))           00001340
      C(IBASE+I) = RIM1B3 - D(IBASE+I-1)*C(IBASE+I-1)                   00001350
  399 CONTINUE                                                          00001360
      D(IBASE+NM1) = 0                                                  00001370
C                                                                       00001380
      DO 500 I1=2,NM1                                                   00001390
      I=NM1+2-I1                                                        00001400
      C(IBASE+I)=B(IBASE+I)*C(IBASE+I)-D(IBASE+I)*C(IBASE+I+1)          00001410
  500 CONTINUE                                                          00001420
C                                                                       00001430
CXX      GO TO (600, 520, 530),  IOPT                                   00001440
CXXC                                                                    00001450
CXXC     FIX UP C(IBASE+NUSE) FOR THE SPECIAL OPTIONS                   00001460
CXXC                                                                    00001470
CXX 520  C(IBASE+NUSE) = C(IBASE+NM1)                                   00001480
CXX      C(IBASE+1) = C(IBASE+2)                                        00001490
CXX      GO TO 600                                                      00001500
CXX 530  C(IBASE+NUSE) = 1.5*(YNP-F)/H - .5*C(IBASE+NM1)                00001510
CXX      C(IBASE+1) = -1.5*(Y1P - (Y(IBASE+2)-Y(IBASE+1))/(X(IBASE+2)   00001520
CXX     1  -X(IBASE+1)))/(X(IBASE+2)-X(IBASE+1))                        00001530
CXX     1   -.5*C(IBASE+2)                                              00001540
CXX      GO TO 600                                                      00001550
C                                                                       00001560
 600  DO 699 I=1,NM1                                                    00001570
      H=X(IBASE+I+1)-X(IBASE+I)                                         00001580
      D(IBASE+I)=(C(IBASE+I+1)-C(IBASE+I))/(3.0*H)                      00001590
      B(IBASE+I)=(Y(IBASE+I+1)-Y(IBASE+I))/H-(H*D(IBASE+I)+C(IBASE+I))*H00001600
 699  CONTINUE                                                          00001610
C                                                                       00001620
C     DEFINE THE N'TH CUBIC AS THE N-1'TH CUBIC                         00001630
C                                                                       00001640
      D(IBASE+NUSE) = D(IBASE+NM1)                                      00001650
      B(IBASE+NUSE) = B(IBASE+NM1) + (2*C(IBASE+NM1) + 3*D(IBASE+NM1)*H)00001660
     * * H                                                              00001660
C                                                                       00001680
      RETURN                                                            00001690
C                                                                       00001700
C     SPECIAL CASE FOR N = 2.                                           00001710
C                                                                       00001720
 800  D(IBASE+1) = 0                                                    00001730
      B(IBASE+1) = F                                                    00001740
      RETURN                                                            00001750
C                                                                       00001760
C     ENTRY FOR EXTRA OPTIONS                                           00001770
C                                                                       00001780
CDCC                                                                    00001790
CDCC     THIS ENTRY DOES NOT WORK IN CDC VERSION OF SUBROUTINE          00001800
CDCC     NOT USED IN PTOLEMY AT THE MOMENT                              00001810
CDCC                                                                    00001820
C                                                                       00001830
CXX      ENTRY  SPLNCX (N, X, Y, B, C, D, IIOPT, Y1P, YNP)              00001840
CXXC                                                                    00001850
CXX      IOPT = IIOPT                                                   00001860
CXX      IF (IOPT .LT. 1  .OR.  IOPT .GT. 3)  RETURN                    00001870
CXX      GO TO (40, 920, 930), IOPT                                     00001880
CXX 920  ASSIGN 200 TO IGO                                              00001890
CXX      GO TO 50                                                       00001900
CXX 930  ASSIGN 220 TO IGO                                              00001910
CXX      GO TO 50                                                       00001920
C                                                                       00001930
      END                                                               00001940
      SUBROUTINE SYSERR                                                 00000010
C                                                                       00000020
C      7/14/72  ARGONNE VERSION                                         00000030
C     7/9/80 - CRAY VERSION ADDED - S.P.                                00000032
C                                                                       00000040
C      GIVES A TRACE BACK AND STOPS EXECUTION                           00000050
C                                                                       00000060
      PRINT 10                                                          00000070
  10  FORMAT ('0' / '0' / '0$*$*$*$*$*$  SYSERR CALLED',  5X, 10('$*$*$*00000080
     *$*$*') )                                                          00000080
CIB      CALL ERRTRA                                                    00000100
CIB      STOP 9876                                                      00000110
CDC      X = -1                                                         00000120
CDC      X = SQRT(X)                                                    00000130
CDC      STOP 7777                                                      00000140
CRA      CALL ABORT ( ' *** EXECUTION TERMINATED BY PTOLEMY ***')       00000145
CRA      STOP 6HSYSERR                                                  00000146
CUN      CALL FTNWB                                                     00000150
CUN      STOP 7777                                                      00000160
      CALL EXIT(1)
      END                                                               00000180
         function wig9j(j1,j2,j3,j4,j5,j6,j7,j8,j9)
C                                                                       00000300
C                                                                       00000400
C   COMPUTES 9-J COEFICENTS.                                            00000500
C                                                                       00000600
C                                                                       00000700
C     THIS ROUTINE COMPUTES THE WIGNER 9-J COEFICENT                    00000800
C                                                                       00000900
C         (  J1  J2  J3  )                                              00001000
C         (  J4  J5  J6  )                                              00001100
C         (  J7  J8  J9  )                                              00001200
C                                                                       00001300
C   THE INPUT CONSISTS OF 9 INTEGERS THAT ARE TWICE THE CORRESPONDING   00001400
C   J'S.  IN THIS MANNER HALF-INTEGER ANGULAR MOMENTA CAN BE            00001500
C   ACCOMODATED.                                                        00001600
C                                                                       00001700
C     THE STANDARD METHOD (COMPUTING A SUM OF PRODUCTS OF THREE 6-J'S)  00001800
C   IS USED.  ALL OF THE FACTORIALS COMMON TO ALL TERMS IN THE SUM ARE  00001900
C   FACTORED OUT AND EVALUATED ONLY ONCE.  SIMILARLY ALL THE INTEGER TO 00002000
C   REAL CONVERSIONS AND OTHER SETUP REQUIRED FOR TH 6-J CALCULATIONS   00002100
C   ARE CARRIED OUT ONLY ONCE.  TO EVALUATE THE 6-J'S A FACTORIZATION OF00002200
C   THE RACAH SUM SUGGESTED BY                                          00002300
C       J. G. WILLS, COMPUTER PHYSICS COMM. 2, 38 (1971)                00002400
C   IS USED.  THE LOGARITHMS OF THE REQUIRED FACTORIALS ARE FOUND       00002500
C   EITHER FROM A TABLE OF  LOG(FACTORIAL)  OR BY USING THE  DLGAMA     00002600
C   FUNCTION.  IN THE PRESENT IMPLIMENTATION, THE TABLE CONTAINS THE    00002700
C   FIRST 1001 FACTORIALS.  THIS SIZE MAY BE REDUCED (TO SAVE CORE) OR  00002800
C   INCREASED (TO SAVE TIME FOR VERY LARGE J'S) BY RECOMPILING THE      00002900
C   /LOGFAC/ COMMON BLOCK.                                              00003000
C                                                                       00003100
C     THE MAGNITUDE OF THE J'S THAT CAN BE ACCOMODATED BY THIS ROUTINE  00003200
C   IS DETERMINED BY CANCELLATIONS IN THE RACAH AND 9-J SUMS.  AS CAN BE00003300
C   SEEN IN THE TABLE BELOW, FAIRLY SEVERE CANCELLATIONS OCCUR WHEN ALL 00003400
C   NINE J'S ARE OF THE ORDER OF 50.  COMPLETE LOSS OF SIGNIFICANCE AND 00003500
C   EXTENSIVE NUMERIC OVERFLOWS OCCUR WHEN ALL THE J'S ARE 75.  HOWEVER 00003600
C   IF FIVE OF THE J'S REMAIN SMALL WHILE THE OTHER FOUR GET LARGE, THEN00003700
C   THE LARGE J'S MAY GET VERY LARGE.                                   00003800
C                                                                       00003900
C     THE NUMBER OF TERMS IN THE 9-J SUM IS                             00004000
C            N = MIN ( J1+J9 , J2+J6 , J4+J8 ) -                        00004100
C                  MAX ( \J1-J9\ , \J2-J6\ , \J4-J8\ )  +  1            00004200
C   NO ATTEMPT IS MADE TO MINIMIZE N BY INTERCHANGING ROWS OR COLUMNS.  00004300
C   IT MAY BE ADVANTAGEOUS TO THE USER TO ARRANGE THE J'S SUCH THAT N   00004400
C   IS MINIMIZED (AT LEAST APPROXIMATELY).  IN PARTICULAR IT SHOULD BE  00004500
C   NOTICED THAT  J3, J5, AND J7  ARE NOT USED IN DETERMINING  N  SO    00004600
C   THAT A ZERO VALUE OF ONE OF THESE IS "WASTED".                      00004700
C                                                                       00004800
C     EACH LINE IN THE FOLLOWING TABLE IS THE RESULT OF SEVERAL THOUSAND00004900
C   COMPUTATIONS WITH J'S IN THE INDICATED RANGES.  ERRORS WERE         00005000
C   DETERMINED BY COMPARISON WITH QUADRUPLE PRECISION 9-J VALUES.       00005100
C   EXECUTION TIMES FOR THE 370/195 ARE GIVEN IN MICROSECONDS.  NONE OF 00005200
C   THE 9-J'S EVALUATED WERE TRIVIALLY ZERO.                            00005300
C                                                                       00005400
C    J1 J2 J5     J3 J6 J9    J4   J7   J8       TIME     RMS RELATIVE  00005500
C     RANGE        RANGE      FIXED VALUES                    ERROR     00005600
C                                                                       00005700

C   1/2 - 5/2    1/2 - 5/2    3/2  3/2  3/2       300         1E-14     00005800
C     5 - 10       5 - 10    15/2   8  17/2      1450         1E-12     00005900
C    15 - 30       0 - 47     25  45/2   20      3000         3E-12     00006000
C                                                                       00006100
C    50 - 54      50 - 54     51   52  107/2    34000         4E-8      00006200
C    25 - 30     1/2 -  5     55/2  5   5/2       450         2E-14     00006300
C    25 - 30       5 - 10     55/2  7  15/2      1600         2E-13     00006400
C      250       1/2 -  5    505/2  2   5/2       440         2E-13     00006500
C      250         5 - 10    505/2  7  15/2      1680         5E-13     00006600
C      500       1/2 -  5   1005/2  2   5/2      4850         2E-12     00006700
C      500         5 - 10   1005/2  7  15/2     13800         2E-12     00006800
C     2500       1/2 -  5   5005/2  2   5/2      4700         2E-12     00006900
C     2500         5 - 10   5005/2  7  15/2     13700         2E-12     00007000
C                                                                       00007100
C                                                                       00007200
C   USAGE:                                                              00007300
C                                                                       00007400
C     WIG9J (2*J1, 2*J2, 2*J3, 2*J4, 2*J5, 2*J6, 2*J7, 2*J8, 2*J9)      00007500
C                                                                       00007600
C   RETURNS THE DOUBLE PRECISION 9-J COEFICENT.  THE TRIANGLE RULES ARE 00007700
C   CHECKED AND 0 IS RETURNED IF THEY ARE VIOLATED.  HOWEVER NO CHECK IS00007800
C   MADE FOR INVALID COMBINATIONS OF HALF-INTEGER AND INTEGERS J'S.  THE00007900
C   USER MUST INSURE THAT TWICE THE SUM OF THE J'S IN EACH ROW AND      00008000
C   COLUMN IS EVEN.  IF THIS IS NOT SO, MEANINGLESS RESULTS WILL BE     00008100
C   RETURNED.                                                           00008200
C                                                                       00008300
C   CORE REQUIREMENTS (BASE 10):                                        00008400
C                                                                       00008500
C       WIG9J    4878                                                   00008600
C      /LOGFAC/  8028    (MAY BE REDUCED OR INCREASED)                  00008700
C       DLGAMA   1008                                                   00008800
C       TOTAL   13914                                                   00008900
C                                                                       00009000
C     IN ADDITION "DEXP" AND THE FORTRAN I/O PACKAGE ARE REQUIRED.      00009100
C   THESE ARE NOT INCLUDED IN THE ABOVE SINCE YOU ARE LIKELY TO USE THEM00009200
C   ANYWAY.  DLGAMA AND /LOGFAC/ WILL BE SHARED WITH THE CLEBSCH AND    00009300
C   6-J ROUTINES SHOULD YOU USE THEM.                                   00009400
C                                                                       00009500
C                                                                       00009600
C   S. PIEPER                                                           00009700
C                                                                       00009800
C   9/6/74 - FIRST VERSION                                              00009900
C     1/18/80 - STANDARD MORTRAN - RPG                                  00010000
C     3/13/80 - CONVERT BACK TO TBIT, ETC. - S.P.                       00010100
C     5/28/80 - CDC FIXES, LAND INSTEAD OF TBIT FOR XHALF - S.P.        00010200
C     7/7/80 - REAL ALGAMA FOR CDC, CRAY - S.P.                         00010300
C     7/15/80 - CRAY VERSION - S.P.                                     00010400
C     7/31/80 - USE LOGFAC TABLAE IN ALLOCATOR - S.P.                   00010500
C                                                                       00010700
C                                                                       00010800
C                                                                       00010900
      IMPLICIT INTEGER (A-C, I-R, T-V, X-Z), REAL*8 (D-H, S, W)         00011000
C                                                                       00011100
CIB      REAL*8  WIG9J                                                  00011200
CNV      REAL ALGAMA                                                    00011300
C                                                                       00011400
	REAL*8 ALLOC
	COMMON/ALLOCS/ALLOC(1)
	CHARACTER*8 LOGNAM / 'LOGFACS ' /
	INTEGER*4 Z
	COMMON /LOCKT/ Z(100)
C                                                                       00011500
      REAL*8 LF                                                         00011600
      COMMON /LOGFAC/ MAXLF, LFBIAS, LF(1)                              00011700
CRA      COMMON /LOGFAC/  MAXLF, LFBIAS, LFTBL(1)                       00011800
C                                                                       00011900
CRA      POINTER ( LFBASE, LF(1) )                                      00012200
C                                                                       00012300
C                                                                       00012400
C     THE FOLLOWING FAKE EXTERNALS DRAG IN THE REQUIRED BLOCK DATAS     00012500
C     FOR THE ABOVE COMMON BLOCK                                        00012600
C                                                                       00012700
      EXTERNAL LOGDUM                                                   00012800
C                                                                       00012900
CCCCC**********                                                         00013000
CCCC      COMMON /NUMBLK/  NUMIT, NUM6J, BIGOUT, BIGIN, SUMLOG, SUMOUT  00013100
CCCCC**********                                                         00013200
C                                                                       00013300
C     WE USE THE FOLLOWING FORTRAN H XL FUNCTIONS:                      00013400
C     TBIT - TESTS TO SEE IF A BIT IS ON.  TBIT(I, 31) IS TRUE          00013500
C            IF  I  IS ODD.                                             00013600
C     SHFTR - SHIFTS RIGHT AND IS EQUIVALENT TO A DIVIDE BY 2**N.       00013700
C                                                                       00013800
CIB      LOGICAL TBIT                                                   00013900
C                                                                       00014000
C                                                                       00014100
C                                                                       00014200
      REAL*8 ZERONE(2)                                                  00014300
        data ZERONE/0.D0, 1.D0/                                         00014300
C                                                                       00014400
      DIMENSION J1S(5), J2S(5), J3S(5)                                  00014500
      DIMENSION AS(4, 3), BS(2, 3), DAS(4, 3), DBS(2, 3),  MX2S(3), MN2S00014700
     *(3), MN1S(3), RS(3), DRS(3),  DMN1S(3), DMN2S(3), DMX2S(3), IDS(2,00014700
     * 3)                                                               00014700
      DIMENSION SUMS(3)                                                 00015000
C                                                                       00015100
C                                                                       00015500
C     FOLLOWING IS FACTORIAL(N-1)                                       00015600
C                                                                       00015700
      DLF(N) = DLGAMA(DFLOAT(N))                                        00015800
C                                                                       00015900
CCCCC****************                                                   00016000
CCCC      NUMIT = 0                                                     00016100
CCCC      NUM6J = 0                                                     00016200
CCCC      BIGOUT = 0                                                    00016300
CCCC      BIGIN = 0                                                     00016400
CCCCC****************                                                   00016500
C                                                                       00016600
C                                                                       00016700
CRA      LFBASE = LOCF(LFTBL) + LFBIAS                                  00016800
C                                                                       00016900
C     PICK UP THE J'S IN THE APPROPRIATE ORDER                          00017000
C     NOTE THAT THE JIS ARRAYS REALLY CONTAIN 2*J                       00017100
C                                                                       00017200
      J1S(1) = J1                                                       00017300
      J1S(2) = J6                                                       00017400
      J1S(3) = J8                                                       00017500
      J2S(1) = J2                                                       00017600
      J2S(2) = J4                                                       00017700
      J2S(3) = J9                                                       00017800
      J3S(1) = J3                                                       00017900
      J3S(2) = J5                                                       00018000
      J3S(3) = J7                                                       00018100
C                                                                       00018200
C     EXPAND THE ARRAYS TO ALLOW CYCLIC INDICES                         00018300
C     AFTER THIS IS DONE THE 3 6-J'S ARE                                00018400
C       (  J1S(I)    J2S(I)    J3S(I)   )                               00018500
C       ( J1S(I+1)  J2S(I+2)     X      )                               00018600
C     FOR I = 1, 2, 3.                                                  00018700
C                                                                       00018800
      DO 59 I = 1, 2                                                    00018900
      J1S(I+3) = J1S(I)                                                 00019000
      J2S(I+3) = J2S(I)                                                 00019100
      J3S(I+3) = J3S(I)                                                 00019200
 59   CONTINUE                                                          00019300
C                                                                       00019400
C     DETERMIN THE RANGE OF 2*X                                         00019500
C                                                                       00019600
      XSTART = IABS(J2S(1) - J1S(2))                                    00019700
      XEND = J2S(1) + J1S(2)                                            00019800
      XSUM = XEND                                                       00019900
      DO 129 I = 2, 3                                                   00020000
      XSTART = MAX0 ( XSTART, IABS(J2S(I)-J1S(I+1)) )                   00020100
      XPIECE = J2S(I) + J1S(I+1)                                        00020200
      XEND = MIN0 ( XEND, XPIECE)                                       00020300
      XSUM = XSUM + XPIECE                                              00020400
 129  CONTINUE                                                          00020500
C                                                                       00020600
      IF (XEND .LT. XSTART) GO TO 900                                   00020700
C                                                                       00020800
      XSUM = XSUM - XEND                                                00020900
C                                                                       00021000
C     EXTRACT 1/2 INTEGER PART OF X AND CONVERT 2*X TO INTPART(X)       00021100
C                                                                       00021200
      XHALF = IAND(XSTART,1)                                            00021300
      XSTART = JISHFT(XSTART, -1)                                         000214
      XEND = JISHFT(XEND, -1)                                             000215
      XSUM = JISHFT(XSUM, -1)                                             000216
      DXHALF = ZERONE(XHALF+1)                                          00021700
      DX = XSTART - 1                                                   00021800
      XSUM = XSUM + 1                                                   00021900
	ILOGFC=NAMLOC(LOGNAM)
	LLOGFC=Z(ILOGFC)-1
C                                                                       00022000
CCCCC******************                                                 00022100
CCCC      NUMIT = XEND - XSTART + 1                                     00022200
CCCCC******************                                                 00022300
C                                                                       00022400
C     SETUP THE J1+J2-J3, ETC                                           00022500
C                                                                       00022600
      DO 159 I = 1, 3                                                   00022700
C                                                                       00022800
C     THE AS, FOR EACH 6-J THEY ARE                                     00022900
C       A1 = J1+J2-J3                                                   00023000
C       A2 = J4+J5-J3                                                   00023100
C       A3 = J1+J5-J6                                                   00023200
C       A4 = J4+J2-J6                                                   00023300
C     IN A3 AND A4 THE  J6 = X  IS LEFT OUT.                            00023400
	JJJ1=J1S(I)+J2S(I)-J3S(I)
	JJJ2=J1S(I+1)+J2S(I+2)-J3S(I)
	JJJ3=J1S(I)+J2S(I+2)
	JJJ4=J2S(I)+J1S(I+1)
      AS(1,I) = JISHFT (JJJ1, -1)                               
      AS(2,I) = JISHFT (JJJ2, -1)                          
      AS(3,I) = JISHFT (JJJ3, -1)                                    
      AS(4,I) = JISHFT (JJJ4, -1)                                   
C                                                                       00024000
C     THE BS, FOR EACH 6-J THEY ARE                                     00024100
C        B1 = J3 + J6 - J1 - J4                                         00024200
C        B2 = J3 + J6 - J2 - J5                                         00024300
C     THE  J6 = X  IS LEFT OUT                                          00024400
C                                                                       00024500
      BS(1,I) = J2S(I+2) - AS(2,I) - AS(3,I)                            00024600
      BS(2,I) = J1S(I) - AS(1,I) - AS(3,I)                              00024700
C                                                                       00024800
 159  CONTINUE                                                          00024900
C                                                                       00025000
C     NOW MAKE SURE NONE OF THEM ARE NEGATIVE WHICH WOULD INDICATE      00025100
C     A VIOLATION OF THE TRIANGLE RULES.                                00025200
C                                                                       00025300
      TEST = IOR(AS(1,1),IOR(AS(2,1),IOR(AS(1,2),IOR(AS(2,2),IOR(AS(1,3)00025400
     *,AS(2,3))))))                                                     00025400
CIB      IF (TBIT(TEST,1))  GO TO 900                                   00025600
CDC      IF ( AND(TEST, -1) .LT. 0 )  GO TO 900                         00025700
	IF(TEST.GT.1000000000)GO TO 900
C                                                                       00025900
C     SETUP THE QUANTITIES FOR EACH OF THE 3 6-J SUMS.                  00026000
C     THESE INCLUDE THE STARTING MIN AND MAX VALUES, ETC.               00026100
C                                                                       00026200
      DO 299 I = 1, 3                                                   00026300
      MN1S(I) = MIN0(AS(1,I), AS(2,I))                                  00026400
      MN2S(I) = MIN0(AS(3,I), AS(4,I))                                  00026500
      MX2S(I) = -MIN0(BS(1,I), BS(2,I))                                 00026600
      RS(I) = 1 + XHALF + AS(3,I) + AS(4,I)                             00026700
      DO 239 J = 1, 4                                                   00026800
      DAS(J,I) = AS(J,I)                                                00026900
 239  CONTINUE                                                          00027000
      DBS(1,I) = BS(1,I)                                                00027100
      DBS(2,I) = BS(2,I)                                                00027200
      DRS(I) = 1 + DXHALF + DAS(3,I) + DAS(4,I)                         00027300
      DMN1S(I) = MN1S(I)                                                00027400
      DMN2S(I) = MN2S(I)                                                00027500
      DMX2S(I) = MX2S(I)                                                00027600
C                                                                       00027700
C     THESE ARE INDICES FOR THE 4 FACTORIALTS THAT DEPEND ON X          00027800
C                                                                       00027900
      IDS(1,I) = AS(2,I) + BS(2,I)                                      00028000
      IDS(2,I) = AS(1,I) + BS(1,I)                                      00028100
 299  CONTINUE                                                          00028200
C                                                                       00028300
C     FIND THE LOG OF THE OUTSIDE FACTORS NOW                           00028400
C     IT WILL BE USED IN EACH TERM OF THE 9-J SUM TO AVOID UNDER/OVER   00028500
C     FLOW PROBLEMS.                                                    00028600
C                                                                       00028700
C                                                                       00028800
C     THESE ARE ALL THE TRIANGLE FACTORIALS THAT COMPLTELY FACTOR OUT   00028900
C     OF THE SUM                                                        00029000
C                                                                       00029100
      SUMOUT = 0                                                        00029200
C                                                                       00029300
C     XSUM IS >= ANY OF THE REQUIRED FACTORIAL ARGUMENTS                00029400
C                                                                       00029500
      IF (XSUM .GT. MAXLF) GO TO 350                                    00029600
C                                                                       00029700
C     WE CAN USE THE LOG(FACTORIAL) TABLE                               00029800
C                                                                       00029900
      DO 339 I = 1, 3                                                   00030000
      DO 339 J = 1, 2                                                   00030100
	SUMOUT=SUMOUT+ALLOC(LLOGFC+1+AS(J,I))
	1-ALLOC(LLOGFC+2+J3S(I)+AS(J,I))+
	2ALLOC(LLOGFC+1-AS(J,I)+J1S(I-1+J))
	3+ALLOC(LLOGFC+1-AS(J,I)+J2S(I-2+2*J))
 339  CONTINUE                                                          00030600
      GO TO 400                                                         00030700
C                                                                       00030800
C     WE MUST USE LOG GAMMA                                             00030900
C                                                                       00031000
 350  DO 359 I = 1, 3                                                   00031100
      DO 359 J = 1, 2                                                   00031200
      SUMOUT = SUMOUT +  DLGAMA(1+DAS(J,I)) - DLF(2+J3S(I)+AS(J,I))  + D00031300
     *LF(1-AS(J,I)+J1S(I-1+J))  + DLF(1-AS(J,I)+J2S(I-2+2*J))           00031300
 359  CONTINUE                                                          00031700
C                                                                       00031800
CCCCC************                                                       00031900
CCCC      BIGOUT = 1                                                    00032000
CCCCC***********                                                        00032100
C                                                                       00032200
 400  SUMOUT = .5*SUMOUT                                                00032300
C                                                                       00032400
C     WE ARE NOW READY TO DO THE 9-J LOOP                               00032500
C                                                                       00032600
      WIG9J = 0                                                         00032700
C                                                                       00032800
      DO 799 X = XSTART, XEND                                           00032900
      DX = DX + 1                                                       00033000
C                                                                       00033100
C     DO THE PARTS OF EACH 6-J THAT DEPEND ON X                         00033200
C                                                                       00033300
      SUMLOG = SUMOUT                                                   00033400
      ISIGN = 0                                                         00033500
C                                                                       00033600
      DO 759 I = 1, 3                                                   00033700
C                                                                       00033800
C     GET THE COEFFICENTS FOR THIS 6-J                                  00033900
C                                                                       00034000
      MX2 = MX2S(I) - X                                                 00034100
      MN2 = MN2S(I) - X                                                 00034200
      NMAX = MN1S(I)                                                    00034300
      DNMAX = DMN1S(I)                                                  00034400
      DNMIN = 0                                                         00034500
      NMIN = 0                                                          00034600
C                                                                       00034700
C     GET THE RANGE OF THE 6-J SUM                                      00034800
C                                                                       00034900
      IF (MN2 .GE. NMAX) GO TO 520                                      00035000
      NMAX = MN2                                                        00035100
      DNMAX = DMN2S(I) - DX                                             00035200
C                                                                       00035300
 520  IF (MX2 .LE. 0) GO TO 540                                         00035400
      NMIN = MX2                                                        00035500
      DNMIN = DMX2S(I) - DX                                             00035600
C                                                                       00035700
 540  SUM = 1                                                           00035800
      ISIGN = ISIGN + NMIN                                              00035900
C                                                                       00036000
C                                                                       00036100
C     THIS ONE TEST CHECKS ALL 12 INEQUALITIES IMPLIED BY THE           00036200
C     4 TRIANGLE CONDITIONS.FOR THE GIVEN 6-J.                          00036300
C                                                                       00036400
      IF (NMAX - NMIN) 900, 700, 600                                    00036500
C                                                                       00036600
 600  D1 = DAS(1,I) - (DNMAX-1)                                         00036700
      D2 = DAS(2,I) - (DNMAX-1)                                         00036800
      D8 = DRS(I) - (DNMAX-1)                                           00036900
      D3 = DAS(3,I) - (DNMAX-1+DX)                                      00037000
      D4 = DAS(4,I) - (DNMAX-1+DX)                                      00037100
      D5 = -DNMAX                                                       00037200
      D6 = -DBS(1,I) - (DX + DNMAX)                                     00037300
      D7 = -DBS(2,I) - (DX + DNMAX)                                     00037400
C                                                                       00037500
      E1 = D1*D2*D3*D4                                                  00037600
      E5 = D5*D6*D7*D8                                                  00037700
C                                                                       00037800
      NMAXM2 = NMAX-2                                                   00037900
      SUM = 1 + E1/E5                                                   00038000
      IF (NMAXM2 .LT. NMIN) GO TO 700                                   00038100
C                                                                       00038200
      F1 = (D1+1)*(D2+1)*(D3+1)*(D4+1) - E1                             00038300
      F5 = (D5+1)*(D6+1)*(D7+1)*(D8+1) - E5                             00038400
C                                                                       00038500
      G1 = (D1+2)*(D2+2)*(D3+2)*(D4+2) - F1 - F1 - E1                   00038600
      G5 = (D5+2)*(D6+2)*(D7+2)*(D8+2) - F5 - F5 - E5                   00038700
C                                                                       00038800
      H1 = (D1+3)*(D2+3)*(D3+3)*(D4+3) - 3*(G1+F1) - E1                 00038900
      H5 = (D5+3)*(D6+3)*(D7+3)*(D8+3) - 3*(G5+F5) - E5                 00039000
C                                                                       00039100
C           COMPUTE THE SUM IN THE RACAH FORMULA                        00039200
C                                                                       00039300
      DO 659 N = NMIN, NMAXM2                                           00039400
      E1 = E1 + F1                                                      00039500
      E5 = E5 + F5                                                      00039600
      F1 = F1 + G1                                                      00039700
      F5 = F5 + G5                                                      00039800
      G1 = G1 + H1                                                      00039900
      G5 = G5 + H5                                                      00040000
      H1 = H1 + 24                                                      00040100
      H5 = H5 + 24                                                      00040200
      SUM = 1 + SUM*(E1/E5)                                             00040300
 659  CONTINUE                                                          00040400
C                                                                       00040500
C                                                                       00040600
700   SUMS(I) = SUM                                                     00040700
C                                                                       00040800
CCCCC***********                                                        00040900
CCCC            NUM6J = NUM6J + NMAX - NMIN + 1                         00041000
CCCCC************                                                       00041100
C                                                                       00041200
C     HERE ARE THE PARTS OF THE FACTORIALS FOR NMIN AND ALSO THOSE      00041300
C     FACTORIALS THAT DEPEND ON X                                       00041400
C                                                                       00041500
      IF (RS(I) .GT. MAXLF) GO TO 730                                   00041600
C                                                                       00041700
C     WE CAN USE THE LOG FACTORIAL TABLE                                00041800
C                                                                       00041900
	SUMLOG=SUMLOG+ALLOC(LLOGFC+1-X+AS(4,I))
	1+ALLOC(LLOGFC+1+X+IDS(1,I))+ALLOC(LLOGFC+1+X+IDS(2,I))-
	2ALLOC(LLOGFC+2+XHALF+X+AS(4,I))+
	3ALLOC(LLOGFC+1-NMIN+RS(I))-ALLOC(LLOGFC+1+NMIN)-
	4ALLOC(LLOGFC+1+NMIN+X+BS(1,I))
	5-ALLOC(LLOGFC+1+NMIN+X+BS(2,I))-
	6ALLOC(LLOGFC+1+AS(1,I)-NMIN)-ALLOC(LLOGFC+1-NMIN+AS(2,I))
	7-ALLOC(LLOGFC+1-NMIN-X+AS(3,I))
	8-ALLOC(LLOGFC+1-NMIN-X+AS(4,I))
      GO TO 759                                                         00042800
C                                                                       00042900
C     WE MUST USE LOG GAMMA                                             00043000
C                                                                       00043100
 730  SUMLOG = SUMLOG +  DLGAMA(1-DX+DAS(4,I)) + DLGAMA(1+DX+DAS(2,I)+DB00043200
     *S(2,I))  + DLGAMA(1+DX+DAS(1,I)+DBS(1,I))  - DLGAMA(2+DXHALF+DX+DA00043200
     *S(4,I))  + DLGAMA(1-DNMIN+DRS(I)) - DLGAMA(1+DNMIN)  - DLGAMA(1+DN00043200
     *MIN+DX+DBS(1,I))  - DLGAMA(1+DNMIN+DX+DBS(2,I))  - DLGAMA(1-DNMIN+00043200
     *DAS(1,I)) - DLGAMA(1-DNMIN+DAS(2,I))  - DLGAMA(1-DNMIN-DX+DAS(3,I)00043200
     *)  - DLGAMA(1-DNMIN-DX+DAS(4,I))                                  00043200
C                                                                       00044200
CCCCC***************                                                    00044300
CCCC      BIGIN = 1                                                     00044400
CCCCC**************                                                     00044500
C                                                                       00044600
 759  CONTINUE                                                          00044700
C                                                                       00044800
C     NOW MULTIPLY THE 3 6-J PARTS TOGEATHER AND SUM                    00044900
C                                                                       00045000
      SUMP = SUMS(1) * SUMS(2) * DEXP(SUMLOG) * SUMS(3)  * (1+DXHALF+(DX00045100
     *+DX))                                                             00045100
C                                                                       00045300
      IF (MOD(ISIGN,2).NE.0) SUMP = -SUMP                               00045400
      WIG9J = WIG9J + SUMP                                              00045500
C                                                                       00045600
 799  CONTINUE                                                          00045700
C                                                                       00045800
C                                                                       00045900
C                                                                       00046000
      RETURN                                                            00046100
C                                                                       00046200
C     THE TRIANGLE RULES ARE NOT SATISFIED                              00046300
C                                                                       00046400
 900  WIG9J = 0                                                         00046500
      RETURN                                                            00046600
C                                                                       00046700
C                                                                       00046900
      END                                                               00047000
