       program cylinder
       implicit real*8 (a-h,o-z)
*
*
       nterm = 75
       wave = 1.d0
       rad = 2.5d0
       open(unit=23,file='cylte.dat',status='unknown',form='formatted')
       open(unit=21,file='cyltm.dat',status='unknown',form='formatted')
       do 100 phi = 0.d0,360.d0,1.d0
         call hrcs(nterm,wave,rad,phi,rcs2)
         call ercs(nterm,wave,rad,phi,rcs)
         write(23,*)phi,10.d0*dlog10(rcs2)
         write(21,*)phi,10.d0*dlog10(rcs)
 100   continue
       close(unit=21)
       close(unit=23)
       stop
       end

	subroutine ERCS(nterm,wave,rad,phi,RCS)
c
c	Calculates exact series solution for RCS of an infinite cylinder.
c
c	E - POLARIZATION.
c
c******************************************************************************
c	INPUTS:
c	nterm = # terms in series. ( nterm-1 = highest order Bessel function.)
c	 wave = wavelength
c	  rad = radius of cylinder
c	  phi = bistatic angle (in degrees).
c
c	OUTPUT: 
c	  RCS = Radar Cross Section 
c
c******************************************************************************
	implicit real*8 (a,b,d-h,o-z)
	implicit complex*16 (c)
c
	dimension aj(1000),by(1000),ch(1000)
	dimension aa(1000),bb(1000),cc(1000)
c
c
	eps  = 1.0 d-7
	ord  = 0.0d0
	pi   = 3.141592654 d0
	cunit = cmplx(0.0d0,1.0d0)
	wn = 2 * pi * rad / wave
c
	call mmbsjn(wn,nterm,aj,ier)
	call mmbsyn(wn,ord,nterm,by,ies)
c
	do 100 n = 1,nterm
100	ch(n) = aj(n) + cunit * by(n)
c
c***	CHANGE  FOR  H - POLARIZATION.
	do 150 n = 1,nterm
	aa(n) = aj(n)
	bb(n) = by(n)
150	cc(n) = ch(n)
c***
c
	csum = aa(1) / cc(1)
c
	do 200 n =2,nterm
	tem = aa(n) * cosd( (n-1) * phi )
cout	if ( abs(tem) .lt. eps ) go to 200
	ctem = tem / cc(n)
c
	csum = csum - 2 * (-1)**n * ctem		
c
200	continue
c
	tem = abs(csum)**2
	RCS = 2 * wave * tem / pi
c
	return
	end

c.............................................................................

	subroutine HRCS(nterm,wave,rad,phi,RCS)
c
c	Calculates exact series solution for RCS of an infinite cylinder.
c
c	H - POLARIZATION.
c
c******************************************************************************
c	INPUTS:
c	nterm = # terms in series. ( nterm-1 = highest order Bessel function.)
c	 wave = wavelength
c	  rad = radius of cylinder
c	  phi = bistatic angle (in degrees).
c
c	OUTPUT: 
c	  RCS = Radar Cross Section 
c
c******************************************************************************
	implicit real*8 (a,b,d-h,o-z)
	implicit complex*16 (c)
c
	dimension aj(1000),by(1000),ch(1000)
	dimension aa(1000),bb(1000),cc(1000)
c
c
	eps  = 1.0 d-7
	ord  = 0.0d0
	pi   = 3.141592654 d0
	cunit = cmplx(0.0d0,1.0d0)
	wn = 2 * pi * rad / wave
c
	mterm = 1 + nterm
	call mmbsjn(wn,mterm,aj,ier)
	call mmbsyn(wn,ord,mterm,by,ies)
c
	do 100 n = 1,nterm
100	ch(n) = aj(n) + cunit * by(n)
c
c***	CHANGE  FOR  H - POLARIZATION.
	aa(1) = - aj(2)
	bb(1) = - by(2)
	cc(1) = - ch(2)
	do 150 n = 2,nterm
	aa(n) = 0.5d0 * ( aj(n-1) - aj(n+1) )
	bb(n) = 0.5d0 * ( by(n-1) - by(n+1) )                
150	cc(n) = 0.5d0 * ( ch(n-1) - ch(n+1) )
c***
c
	csum = aa(1) / cc(1)
c
	do 200 n =2,nterm
	tem = aa(n) * cosd( (n-1) * phi )
cout	if ( abs(tem) .lt. eps ) go to 200
	ctem = tem / cc(n)
c
	csum = csum - 2 * (-1)**n * ctem		
c
200	continue
c
	tem = abs(csum)**2
	RCS = 2 * wave * tem / pi
c
	return
	end

C   IMSL ROUTINE NAME   - MMBSJN                                        MMBA0010
C                                                                       MMBA0020
C-----------------------------------------------------------------------MMBA0030
C                                                                       MMBA0040
C   COMPUTER            - VAX/DOUBLE                                    MMBA0050
C                                                                       MMBA0060
C   LATEST REVISION     - JUNE 1, 1982                                  MMBA0070
C                                                                       MMBA0080
C   PURPOSE             - BESSEL FUNCTION OF THE FIRST KIND OF          MMBA0090
C                           NONNEGATIVE INTEGER ORDER FOR               MMBA0100
C                           REAL ARGUMENTS                              MMBA0110
C                                                                       MMBA0120
C   USAGE               - CALL MMBSJN (ARG,N,B,IER)                     MMBA0130
C                                                                       MMBA0140
C   ARGUMENTS    ARG    - INPUT ARGUMENT. THE ABSOLUTE VALUE OF ARG MUSTMMBA0150
C                           BE LESS THAN OR EQUAL TO 100000. ARG MUST BEMMBA0160
C                           TYPED APPROPRIATELY IN THE CALLING PROGRAM. MMBA0170
C                           (SEE THE PRECISION/HARDWARE SECTION.)       MMBA0180
C                N      - INPUT PARAMETER SPECIFYING THE NUMBER OF      MMBA0190
C                           FUNCTION VALUES TO BE COMPUTED.             MMBA0200
C                B      - OUTPUT VECTOR OF LENGTH N CONTAINING THE      MMBA0210
C                           COMPUTED FUNCTION VALUES. B MUST BE TYPED   MMBA0220
C                           APPROPRIATELY IN THE CALLING PROGRAM.       MMBA0230
C                           (SEE THE PRECISION/HARDWARE SECTION.)       MMBA0240
C                           B(1) WILL CONTAIN THE COMPUTED VALUE FOR    MMBA0250
C                           ORDER ZERO, B(2) WILL CONTAIN THE COMPUTED  MMBA0260
C                           VALUE FOR ORDER 1, B(3) FOR ORDER 2, ETC.   MMBA0270
C                IER    - ERROR PARAMETER. (OUTPUT)                     MMBA0280
C                         TERMINAL ERROR                                MMBA0290
C                           IER = 129 INDICATES THAT EITHER ARG OR N IS MMBA0300
C                             OUT OF RANGE. B(I), (I=1,N) IS SET TO     MMBA0310
C                             MACHINE INFINITY.                         MMBA0320
C                           IER = 129 + J INDICATES THAT B(I), (I=1,J)  MMBA0330
C                             ARE COMPUTED TO MACHINE PRECISION, BUT    MMBA0340
C                             PRECISION IS LOST FOR B(I), (I=J+1,N.)    MMBA0350
C                             SEE THE PROGRAMMING NOTES.                MMBA0360
C                                                                       MMBA0370
C   PRECISION/HARDWARE  - DOUBLE/H32,H36                                MMBA0380
C                       - SINGLE/H48,H60                                MMBA0390
C                                                                       MMBA0400
C   REQD. IMSL ROUTINES - UERTST,UGETIO                                 MMBA0410
C                                                                       MMBA0420
C   NOTATION            - INFORMATION ON SPECIAL NOTATION AND           MMBA0430
C                           CONVENTIONS IS AVAILABLE IN THE MANUAL      MMBA0440
C                           INTRODUCTION OR THROUGH IMSL ROUTINE UHELP  MMBA0450
C                                                                       MMBA0460
C   COPYRIGHT           - 1982 BY IMSL, INC. ALL RIGHTS RESERVED.       MMBA0470
C                                                                       MMBA0480
C   WARRANTY            - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN MMBA0490
C                           APPLIED TO THIS CODE. NO OTHER WARRANTY,    MMBA0500
C                           EXPRESSED OR IMPLIED, IS APPLICABLE.        MMBA0510
C                                                                       MMBA0520
C-----------------------------------------------------------------------MMBA0530
C                                                                       MMBA0540
      SUBROUTINE MMBSJN (ARG,N,B,IER)                                   MMBA0550
C                                  SPECIFICATIONS FOR ARGUMENTS         MMBA0560
      INTEGER            N,IER                                          MMBA0570
      DOUBLE PRECISION   ARG,B(1)                                       MMBA0580
C                                  SPECIFICATIONS FOR LOCAL VARIABLES   MMBA0590
      INTEGER            L,LARGEX,MAGX,NCALC,NN,NBMX,M,NSTART,NEND      MMBA0600
      DOUBLE PRECISION   TEST,TEMPA,TEMPB,TEMPC,EXPARG,P                MMBA0610
      DOUBLE PRECISION   RSIGN,SUM,TOVER,PLAST,POLD,PSAVE,PSAVEL        MMBA0620
      DOUBLE PRECISION   DSIG,RTEN,TMPA4,SMALLX,XINF                    MMBA0630
      DATA               DSIG/16.8D0/                                   MMBA0640
      DATA               RTEN/38.D0/                                    MMBA0650
      DATA               EXPARG/88.02969193111305D0/                    MMBA0660
      DATA               XINF/1.700000000D+38/                          MMBA0670
      DATA               LARGEX /100000/                                MMBA0680
C                                  FIRST EXECUTABLE STATEMENT           MMBA0690
      IER = 0                                                           MMBA0700
      TEMPA = DABS(ARG)                                                 MMBA0710
      MAGX = IFIX(SNGL(TEMPA))                                          MMBA0720
      IF(N.GT.0 .AND. MAGX.LE.LARGEX) GO TO 10                          MMBA0730
C                                  ERROR RETURN -- ARG,N IS OUT OF RANGEMMBA0740
      IER = 129                                                         MMBA0750
      B(1) = XINF                                                       MMBA0760
      IF(N.LT.2) GO TO 9000                                             MMBA0770
      DO 5 L=2,N                                                        MMBA0780
         B(L) = XINF                                                    MMBA0790
    5 CONTINUE                                                          MMBA0800
      GO TO 9000                                                        MMBA0810
   10 RSIGN = 1.D0                                                      MMBA0820
      NCALC = N                                                         MMBA0830
C                                  USE 2-TERM ASCENDING SERIES FOR      MMBA0840
C                                    SMALL ARG                          MMBA0850
      TMPA4 = TEMPA**4.D0                                               MMBA0860
      SMALLX = .1D0**DSIG                                               MMBA0870
      IF(TMPA4.GE.SMALLX) GO TO 20                                      MMBA0880
C                                  TWO-TERM ASCENDING SERIES FOR        MMBA0890
C                                    SMALL ARG                          MMBA0900
      TEMPA = 1.D0                                                      MMBA0910
      TEMPB = -.25D0*ARG*ARG*RSIGN                                      MMBA0920
      B(1) = 1.D0+TEMPB                                                 MMBA0930
      IF(N.EQ.1) GO TO 9005                                             MMBA0940
      DO 15 NN=2,N                                                      MMBA0950
         TEMPA = TEMPA*ARG/DBLE(FLOAT(2*NN-2))                          MMBA0960
         B(NN) = TEMPA*(1.D0+TEMPB/DBLE(FLOAT(NN)))                     MMBA0970
   15 CONTINUE                                                          MMBA0980
      GO TO 9005                                                        MMBA0990
C                                  INITIALIZE THE CALCULATION OF P*S    MMBA1000
   20 NBMX = N-MAGX                                                     MMBA1010
      NN = MAGX+1                                                       MMBA1020
      PLAST = 1.D0                                                      MMBA1030
      P = DBLE(FLOAT(2*NN))/TEMPA                                       MMBA1040
C                                  CALCULATE GENERAL SIGNIFICANCE TEST  MMBA1050
      TEST = 2.D0*1.D1**DSIG                                            MMBA1060
      M = 0                                                             MMBA1070
      IF(NBMX.LT.3) GO TO 30                                            MMBA1080
C                                  CALCULATE P*S UNTIL NN=N-1.          MMBA1090
C                                    CHECK FOR POSSIBLE OVERFLOW.       MMBA1100
      TOVER = 1.D1**(RTEN-DSIG)                                         MMBA1110
      NSTART = MAGX+2                                                   MMBA1120
      NEND = N-1                                                        MMBA1130
      DO 25 NN=NSTART,NEND                                              MMBA1140
         POLD = PLAST                                                   MMBA1150
         PLAST = P                                                      MMBA1160
         P = DBLE(FLOAT(2*NN))*PLAST/TEMPA-RSIGN*POLD                   MMBA1170
         IF(P-TOVER) 25, 25, 35                                         MMBA1180
   25 CONTINUE                                                          MMBA1190
      NN = NEND                                                         MMBA1200
C                                  CALCULATE SPECIAL SIGNIFICANCE TEST  MMBA1210
C                                    FOR NBMX.GT.2.                     MMBA1220
C                                                                       MMBA1230
      TEST = DMAX1(TEST,DSQRT(PLAST*1.D1**DSIG)*DSQRT(2.D0*P))          MMBA1240
C                                                                       MMBA1250
C                                  CALCULATE P*S UNTIL SIGNIFICANCE     MMBA1260
C                                    TEST PASSES                        MMBA1270
   30 NN = NN+1                                                         MMBA1280
      POLD = PLAST                                                      MMBA1290
      PLAST = P                                                         MMBA1300
      P = DBLE(FLOAT(2*NN))*PLAST/TEMPA-RSIGN*POLD                      MMBA1310
      IF(P.LT.TEST) GO TO 30                                            MMBA1320
      IF(M.EQ.1) GO TO 55                                               MMBA1330
C                                  FOR J*S, A STRONG VARIANT OF THE TESTMMBA1340
C                                    IS NECESSARY. CALCULATE IT, AND    MMBA1350
C                                    CALCULATE P*S UNTIL THIS TEST IS   MMBA1360
C                                    PASSED.                            MMBA1370
      M = 1                                                             MMBA1380
      TEMPB = P/PLAST                                                   MMBA1390
      TEMPC = DBLE(FLOAT(NN+1))/TEMPA                                   MMBA1400
      IF(TEMPB+1.D0/TEMPB.GT.2.D0*TEMPC)TEMPB=TEMPC+DSQRT(TEMPC**2-1.D0)MMBA1410
      TEST = TEST/DSQRT(TEMPB-1.D0/TEMPB)                               MMBA1420
      IF(P-TEST) 30, 55, 55                                             MMBA1430
C                                  TO AVOID OVERFLOW, DIVIDE P*S BY     MMBA1440
C                                    TOVER.  CALCULATE P*S UNTIL        MMBA1450
C                                    ABS(P).GT.1.                       MMBA1460
   35 TOVER = 1.D1**RTEN                                                MMBA1470
      P = P/TOVER                                                       MMBA1480
      PLAST = PLAST/TOVER                                               MMBA1490
      PSAVE = P                                                         MMBA1500
      PSAVEL = PLAST                                                    MMBA1510
      NSTART = NN+1                                                     MMBA1520
   40 NN = NN+1                                                         MMBA1530
      POLD = PLAST                                                      MMBA1540
      PLAST = P                                                         MMBA1550
      P = DBLE(FLOAT(2*NN))*PLAST/TEMPA-RSIGN*POLD                      MMBA1560
      IF(P.LE.1.D0) GO TO 40                                            MMBA1570
      TEMPB = DBLE(FLOAT(2*NN))/TEMPA                                   MMBA1580
      TEMPC = .5D0*TEMPB                                                MMBA1590
      TEMPB = PLAST/POLD                                                MMBA1600
      IF(TEMPB+1.D0/TEMPB.GT.2.D0*TEMPC)TEMPB=TEMPC+DSQRT(TEMPC**2-1.D0)MMBA1610
C                                                                       MMBA1620
C                                  CALCULATE BACKWARD TEST, AND FIND    MMBA1630
C                                    NCALC, THE HIGHEST NN SUCH THAT THEMMBA1640
C                                    TEST IS PASSED.                    MMBA1650
      TEST = .5D0*POLD*PLAST*(1.D0-1.D0/TEMPB**2)/1.D1**DSIG            MMBA1660
      P = PLAST*TOVER                                                   MMBA1670
      NN = NN-1                                                         MMBA1680
      NEND = MIN0(N,NN)                                                 MMBA1690
      DO 45 NCALC=NSTART,NEND                                           MMBA1700
         POLD = PSAVEL                                                  MMBA1710
         PSAVEL = PSAVE                                                 MMBA1720
         PSAVE = DBLE(FLOAT(2*NN))*PSAVEL/TEMPA-RSIGN*POLD              MMBA1730
         IF(PSAVE*PSAVEL-TEST) 45, 45, 50                               MMBA1740
   45 CONTINUE                                                          MMBA1750
      NCALC = NEND+1                                                    MMBA1760
   50 NCALC = NCALC-1                                                   MMBA1770
C                                  THE SUM B(1)+2B(3)+2B(5)... IS USED  MMBA1780
C                                    TO NORMALIZE. M, THE COEFFICIENT OFMMBA1790
C                                    B(NN), IS INITIALIZED TO 2 OR 0.   MMBA1800
   55 NN = NN+1                                                         MMBA1810
      M = 2*NN-4*(NN/2)                                                 MMBA1820
C                                  INITIALIZE THE BACKWARD RECURSION ANDMMBA1830
C                                    THE NORMALIZATION SUM              MMBA1840
      TEMPB = 0.D0                                                      MMBA1850
      TEMPA = 1.D0/P                                                    MMBA1860
      SUM = DBLE(FLOAT(M))*TEMPA                                        MMBA1870
      NEND = NN-N                                                       MMBA1880
      IF(NEND) 80, 70, 60                                               MMBA1890
C                                  RECUR BACKWARD VIA DIFFERENCE        MMBA1900
C                                    EQUATION, CALCULATING (BUT NOT     MMBA1910
C                                    STORING) B(NN), UNTIL NN=N.        MMBA1920
   60 DO 65 L=1,NEND                                                    MMBA1930
         NN = NN-1                                                      MMBA1940
         TEMPC = TEMPB                                                  MMBA1950
         TEMPB = TEMPA                                                  MMBA1960
         TEMPA = (DBLE(FLOAT(2*NN))*TEMPB)/ARG-RSIGN*TEMPC              MMBA1970
         M = 2-M                                                        MMBA1980
         SUM = SUM+DBLE(FLOAT(M))*TEMPA                                 MMBA1990
   65 CONTINUE                                                          MMBA2000
C                                  STORE B(NN)                          MMBA2010
   70 B(NN) = TEMPA                                                     MMBA2020
      IF(N.GT.1) GO TO 75                                               MMBA2030
C                                  N=1.  SINCE 2*TEMPA IS ADDED TO THE  MMBA2040
C                                    SUM, TEMPA MUST BE SUBTRACTED      MMBA2050
      SUM = SUM-TEMPA                                                   MMBA2060
      GO TO 110                                                         MMBA2070
C                                  CALCULATE AND STORE B(NN-1)          MMBA2080
   75 NN = NN-1                                                         MMBA2090
      B(NN) = (DBLE(FLOAT(2*NN))*TEMPA)/ARG-RSIGN*TEMPB                 MMBA2100
      IF(NN.EQ.1) GO TO 105                                             MMBA2110
      M = 2-M                                                           MMBA2120
      SUM = SUM+DBLE(FLOAT(M))*B(NN)                                    MMBA2130
      GO TO 90                                                          MMBA2140
C                                  NN.LT.N, SO STORE B(NN) AND SET      MMBA2150
C                                  HIGHER ORDERS TO ZERO                MMBA2160
   80 B(NN) = TEMPA                                                     MMBA2170
      NEND = -NEND                                                      MMBA2180
      DO 85 L=1,NEND                                                    MMBA2190
         ITEMP = NN+L                                                   MMBA2200
         B(ITEMP) = 0.0D0                                               MMBA2210
   85 CONTINUE                                                          MMBA2220
   90 NEND = NN-2                                                       MMBA2230
      IF(NEND.EQ.0) GO TO 100                                           MMBA2240
C                                  CALCULATE VIA DIFFERENCE EQUATION ANDMMBA2250
C                                    STORE B(NN), UNTIL NN=2            MMBA2260
      DO 95 L=1,NEND                                                    MMBA2270
         NN = NN-1                                                      MMBA2280
         B(NN) = (DBLE(FLOAT(2*NN))*B(NN+1))/ARG-RSIGN*B(NN+2)          MMBA2290
         M = 2-M                                                        MMBA2300
         SUM = SUM+DBLE(FLOAT(M))*B(NN)                                 MMBA2310
   95 CONTINUE                                                          MMBA2320
C                                  CALCULATE B(1)                       MMBA2330
  100 B(1) = 2.D0*B(2)/ARG-RSIGN*B(3)                                   MMBA2340
  105 SUM = SUM+B(1)                                                    MMBA2350
C                                  NORMALIZE--IF IZE=1, DIVIDE SUM BY   MMBA2360
C                                    COSH(ARG). DIVIDE ALL B(NN) BY SUM.MMBA2370
  110 CONTINUE                                                          MMBA2380
      DO 115 NN=1,N                                                     MMBA2390
  115 B(NN) = B(NN)/SUM                                                 MMBA2400
      IF(NCALC.EQ.N) GO TO 9005                                         MMBA2410
      IER = 129+NCALC                                                   MMBA2420
 9000 CONTINUE                                                          MMBA2430
      CALL UERTST(IER,6HMMBSJN)                                         MMBA2440
 9005 RETURN                                                            MMBA2450
      END                                                               MMBA2460
C   IMSL ROUTINE NAME   - MMBSYN                                        MMB20010
C                                                                       MMB20020
C-----------------------------------------------------------------------MMB20030
C                                                                       MMB20040
C   COMPUTER            - VAX/DOUBLE                                    MMB20050
C                                                                       MMB20060
C   LATEST REVISION     - JUNE 1, 1981                                  MMB20070
C                                                                       MMB20080
C   PURPOSE             - BESSEL FUNCTION OF THE SECOND KIND OF NON-    MMB20090
C                           NEGATIVE REAL FRACTIONAL ORDER FOR REAL     MMB20100
C                           POSITIVE ARGUMENTS                          MMB20110
C                                                                       MMB20120
C   USAGE               - CALL MMBSYN (ARG,ORDER,N,YN,IER)              MMB20130
C                                                                       MMB20140
C   ARGUMENTS    ARG    - INPUT ARGUMENT. ARG MUST BE TYPED APPRO-      MMB20150
C                           PRIATELY IN THE CALLING PROGRAM. (SEE THE   MMB20160
C                           PRECISION/HARDWARE SECTION.) ARG MUST BE    MMB20170
C                           GREATER THAN OR EQUAL TO XMIN AND LESS THAN MMB20180
C                           OR EQUAL TO XMAX. XMIN IS OF THE ORDER OF   MMB20190
C                           10**(-38) AND XMAX IS OF THE ORDER OF 10**8.MMB20200
C                           THE EXACT VALUES OF XMIN AND XMAX MAY ALLOW MMB20210
C                           LARGER RANGES FOR ARG ON SOME COMPUTERS. SEEMMB20220
C                           THE PROGRAMMING NOTES IN THE MANUAL FOR THE MMB20230
C                           EXACT VALUES.                               MMB20240
C                ORDER  - INPUT VALUE SPECIFYING THE DESIRED ORDER OF   MMB20250
C                           THE BESSEL FUNCTION. ORDER MUST BE TYPED    MMB20260
C                           APPROPRIATELY IN THE CALLING PROGRAM. (SEE  MMB20270
C                           THE PRECISION/HARDWARE SECTION.) ORDER MUST MMB20280
C                           BE GREATER THAN OR EQUAL TO ZERO AND LESS   MMB20290
C                           THAN ONE.                                   MMB20300
C                N      - INPUT INTEGER VALUE SPECIFYING THE NUMBER     MMB20310
C                           OF FUNCTION VALUES TO BE COMPUTED.          MMB20320
C                YN     - OUTPUT VECTOR OF LENGTH N CONTAINING THE      MMB20330
C                           COMPUTED FUNCTION VALUES. YN MUST BE TYPED  MMB20340
C                           APPROPRIATELY IN THE CALLING PROGRAM. (SEE  MMB20350
C                           THE PRECISION/HARDWARE SECTION.) YN(1) WILL MMB20360
C                           CONTAIN THE COMPUTED VALUE FOR THE INPUT    MMB20370
C                           ORDER, YN(2) WILL CONTAIN THE COMPUTED      MMB20380
C                           FUNCTION VALUE FOR ORDER + 1, YN(3) FOR     MMB20390
C                           ORDER + 2, ETC.                             MMB20400
C                IER    - ERROR PARAMETER. (OUTPUT)                     MMB20410
C                         TERMINAL ERROR.                               MMB20420
C                           IER = 129 INDICATES THAT ORDER IS OUT OF    MMB20430
C                             RANGE OR ARG IS LESS THAN XMIN.           MMB20440
C                             YN(I), (I= 1,N) IS SET TO NEGATIVE        MMB20450
C                             MACHINE INFINITY.                         MMB20460
C                           IER = 130 INDICATES THAT N IS LESS THAN OR  MMB20470
C                             EQUAL TO ZERO. YN(1) IS SET TO NEGATIVE   MMB20480
C                             MACHINE INFINITY.                         MMB20490
C                           IER = 131 INDICATES THAT OVERFLOW WOULD     MMB20500
C                             HAVE OCCURED IN VECTOR YN, AT ELEMENT J.  MMB20510
C                             YN(I), (I = J,N) IS SET TO NEGATIVE       MMB20520
C                             MACHINE INFINITY.                         MMB20530
C                           IER = 132 INDICATES THAT ARG IS GREATER     MMB20540
C                             THAN XMAX.  YN(I), (I=1,N) IS SET TO      MMB20550
C                             ZERO.                                     MMB20560
C                                                                       MMB20570
C   PRECISION/HARDWARE  - DOUBLE/H32,H36                                MMB20580
C                       - SINGLE/H48,H60                                MMB20590
C                                                                       MMB20600
C   REQD. IMSL ROUTINES - UERTST,UGETIO                                 MMB20610
C                                                                       MMB20620
C   NOTATION            - INFORMATION ON SPECIAL NOTATION AND           MMB20630
C                           CONVENTIONS IS AVAILABLE IN THE MANUAL      MMB20640
C                           INTRODUCTION OR THROUGH IMSL ROUTINE UHELP  MMB20650
C                                                                       MMB20660
C   COPYRIGHT           - 1978 BY IMSL, INC. ALL RIGHTS RESERVED.       MMB20670
C                                                                       MMB20680
C   WARRANTY            - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN MMB20690
C                           APPLIED TO THIS CODE. NO OTHER WARRANTY,    MMB20700
C                           EXPRESSED OR IMPLIED, IS APPLICABLE.        MMB20710
C                                                                       MMB20720
C-----------------------------------------------------------------------MMB20730
C                                                                       MMB20740
      SUBROUTINE MMBSYN (ARG,ORDER,N,YN,IER)                            MMB20750
C                                  SPECIFICATIONS FOR ARGUMENTS         MMB20760
      INTEGER            N,IER                                          MMB20770
      DOUBLE PRECISION   ARG,ORDER,YN(N)                                MMB20780
C                                  SPECIFICATIONS FOR LOCAL VARIABLES   MMB20790
      INTEGER            I,ITWPI1(2),ITWPI2(2),J,K,L,M,NMAX,NN          MMB20800
      INTEGER            IEND,MEND,M1,M2,M3,IENDP2,MENDM3               MMB20810
      DOUBLE PRECISION   CO,CUM,CUMP,C2,C2P,DELX,DEN,DENP,              MMB20820
     1                   DLNHAF,DLNX,DLNXNU,FACT(34),FOOLU,F1(17),      MMB20830
     2                   F2(16),F3(17),GAM,GNU,GNU1,GNU2,GP(8),GQ(8),   MMB20840
     3                   HALF,P(12,7),PI2,Q(11,7),                      MMB20850
     4                   RTPI2,RTXNU,S,SCRAT(6),SUM,SUMP,T,TWOPI,       MMB20860
     5                   TWOPI1,TWOPI2,T1,T2,T3(6),V,VCOS,VPI,VSIN,     MMB20870
     6                   W(19),X,XC,XIN,XINF,XK,XM,XNU,XNUM,XNUMP,      MMB20880
     7                   Z,ZH,ZMAX,ZMIN,ZSTOP,CON1,CON2                 MMB20890
      DOUBLE PRECISION   TWOPIA,TWOPIB,ZMAX2                            MMB20900
      LOGICAL            NFLAG                                          MMB20910
C                                                                       MMB20920
C                                  MACHINE DEPENDENT CONSTANTS          MMB20930
C                                                                       MMB20940
C                                     PI2 = 2 / PI                      MMB20950
C                                     RTPI2 = SQRT(2 /PI)               MMB20960
C                                     DLNHAF = LOG(.5)                  MMB20970
C                                     TWOPI = 2 * PI                    MMB20980
C                                     TWOPI1 + TWOPI2 = 2 * PI TO       MMB20990
C                                        EXTRA PRECISION                MMB21000
C                                     XINF = LARGEST POSITIVE MACHINE   MMB21010
C                                        NUMBER                         MMB21020
C                                     ZMAX = 2**32, LARGEST             MMB21030
C                                        ACCEPTABLE ARGUMENT            MMB21040
C                                     ZMAX1 = SMALLEST FLOATING--       MMB21050
C                                        POINT CONSTANT WITH ENTIRELY   MMB21060
C                                        INTEGER REPRESENTATION         MMB21070
C                                     ZMIN = 2**(-128), SMALLEST        MMB21080
C                                        ACCEPTABLE ARGUMENT            MMB21090
      DATA PI2   / .636619772367581343D+00/                             MMB21100
      DATA RTPI2 / .797884560802865356D+00/                             MMB21110
      DATA DLNHAF/-.693147180559945310D+00/                             MMB21120
      DATA TWOPI1/0.62831850051879880D 01/                              MMB21130
      DATA TWOPI2/0.30199159819567520D-06/                              MMB21140
      DATA TWOPIA/0.62812500000000000D 01/                              MMB21150
      DATA TWOPIB/0.19350051879882810D-02/                              MMB21160
      DATA TWOPI / .628318530717958647D+01/                             MMB21170
      DATA HALF  / .500000000000000000D+00/                             MMB21180
      DATA ZMAX2 / 4294967296.0D0/                                      MMB21190
      DATA ZMIN  / .294D-38/                                            MMB21200
      DATA XINF  / 1.7D+38/                                             MMB21210
      DATA GNU1  / .500488281250000000D-01/                             MMB21220
      DATA GNU2  / .996093750000000000D+00/                             MMB21230
      DATA ZMAX  / .429496729600000000D+10/                             MMB21240
C                                  COEFFICIENTS FOR APPROXIMATION OF    MMB21250
C                                    GAMMA(V) 1.0 .LE. V .LE. 2.0       MMB21260
      DATA CO    / .886226925452758014D+00/                             MMB21270
      DATA GP    /-.160243404626643007D+01, .148472957685686770D+02,    MMB21280
     *        -.280750941362157079D+03,-.417571326141158257D+03,        MMB21290
     *         .594221994599142643D+03,-.273925800230766587D+05,        MMB21300
     *        -.778326730220496473D+05,-.577828296022400365D+04/        MMB21310
      DATA GQ    /-.178681796751286802D+06,-.268406415244627174D+02,    MMB21320
     *         .214414834924538890D+03,-.224026496383528068D+03,        MMB21330
     *        -.459374814623425838D+04, .145077504178693444D+05,        MMB21340
     *         .329089709308774549D+05,-.114818917860301046D+06/        MMB21350
C                                                                       MMB21360
C                                  THE FOLLOWING DATA STATEMENTS        MMB21370
C                                    INITIALIZE COEFFICIENTS FOR        MMB21380
C                                    RATIONAL APPROXIMATIONS TO         MMB21390
C                                    Y-SUB-NU(X) FOR X AN INTEGER, 0.0  MMB21400
C                                    .LE. NU .LE. 2.0. THE ARRAYS PN    MMB21410
C                                    INITIALIZE NUMERATOR COEFFICIENTS  MMB21420
C                                    FOR X = N, WHILE THE ARRAYS QN     MMB21430
C                                    INITIALIZE THE CORRESPONDING       MMB21440
C                                    DENOMINATOR POLYNOMIALS.           MMB21450
C                                                                       MMB21460
      DATA P(1,1),P(2,1),P(3,1),P(4,1),P(5,1),P(6,1),P(7,1),            MMB21470
     *     P(8,1),P(9,1),P(10,1),P(11,1),P(12,1)                        MMB21480
     *           / .966220363059263586D+02,-.135166484200055635D+04,    MMB21490
     *         .864912264696622203D+04,-.265430829974670190D+05,        MMB21500
     *        -.258277503715144970D+04, .253958451563109432D+06,        MMB21510
     *        -.741666962661522701D+06,-.185303404753455436D+06,        MMB21520
     *         .249424577877157393D+07,-.328627552380477017D+07,        MMB21530
     *        -.588915321642776927D+07, .456049615140191762D+06/        MMB21540
      DATA P(1,2),P(2,2),P(3,2),P(4,2),P(5,2),P(6,2),P(7,2),            MMB21550
     *     P(8,2),P(9,2),P(10,2)                                        MMB21560
     *           / .827836377496675271D+02,-.126469763244904028D+04,    MMB21570
     *         .580689617323250046D+04, .212643684799830998D+05,        MMB21580
     *        -.229280692127117227D+06, .172961254163067550D+06,        MMB21590
     *         .265359659939117365D+07,-.205519713780387999D+07,        MMB21600
     *        -.834208771374013572D+07, .228340012506985874D+06/        MMB21610
      DATA P(1,3),P(2,3),P(3,3),P(4,3),P(5,3),P(6,3),P(7,3),            MMB21620
     *     P(8,3),P(9,3),P(10,3)                                        MMB21630
     *           /-.655464718390655217D+02, .164166175577358534D+04,    MMB21640
     *        -.106542820607880074D+05,-.319526554910620985D+05,        MMB21650
     *         .488095894287561287D+06,-.108912380388916879D+06,        MMB21660
     *        -.626824801218970995D+07, .138102270261905953D+07,        MMB21670
     *         .184668308157951082D+08, .100645521111378579D+07/        MMB21680
      DATA P(1,4),P(2,4),P(3,4),P(4,4),P(5,4),P(6,4),P(7,4),            MMB21690
     *     P(8,4),P(9,4),P(10,4)                                        MMB21700
     *           / .696963123814590427D+02,-.133020689714923170D+04,    MMB21710
     *         .452303293370654331D+04, .471702732719340836D+05,        MMB21720
     *        -.260746560822329284D+06,-.494295070906759632D+06,        MMB21730
     *         .333646692280230407D+07, .169164170133818384D+07,        MMB21740
     *        -.948024850300566443D+07,-.136168612925693158D+07/        MMB21750
      DATA P(1,5),P(2,5),P(3,5),P(4,5),P(5,5),P(6,5),P(7,5),            MMB21760
     *     P(8,5),P(9,5),P(10,5)                                        MMB21770
     *           /-.887596516599601438D+02, .145281203139013753D+04,    MMB21780
     *        -.250210995757854514D+04,-.600773443105863251D+05,        MMB21790
     *         .213021520919733862D+06, .770892892749782111D+06,        MMB21800
     *        -.295732940068198094D+07,-.311681822772432647D+07,        MMB21810
     *         .869420380824761184D+07, .211470077845670982D+07/        MMB21820
      DATA P(1,6),P(2,6),P(3,6),P(4,6),P(5,6),P(6,6),P(7,6),            MMB21830
     *     P(8,6),P(9,6),P(10,6)                                        MMB21840
     *           / .105141263149205264D+03,-.161405265322962398D+04,    MMB21850
     *         .154355389910810923D+04, .728265946613332423D+05,        MMB21860
     *        -.208095986994592776D+06,-.104246097148516522D+07,        MMB21870
     *         .312777482643423898D+07, .474231617027130897D+07,        MMB21880
     *        -.951912134902371432D+07,-.338528003709295307D+07/        MMB21890
      DATA P(1,7),P(2,7),P(3,7),P(4,7),P(5,7),P(6,7),P(7,7),            MMB21900
     *     P(8,7),P(9,7),P(10,7),P(11,7)                                MMB21910
     *           / .336889850414961424D+02,-.422443710176094398D+03,    MMB21920
     *        -.203831469281181514D+04, .331030626094266419D+05,        MMB21930
     *         .604252878602211945D+05,-.933840931715176502D+06,        MMB21940
     *        -.113618427457866563D+07, .959057577413757707D+07,        MMB21950
     *         .903854881191026607D+07,-.236832975247582365D+08,        MMB21960
     *        -.104340643565508042D+08/                                 MMB21970
      DATA Q(1,1),Q(2,1),Q(3,1),Q(4,1),Q(5,1),Q(6,1),Q(7,1),            MMB21980
     *     Q(8,1),Q(9,1),Q(10,1),Q(11,1)                                MMB21990
     *           /-.101241970477141194D+02, .216005595832218891D+02,    MMB22000
     *        -.102815919664190652D+03, .114516150295087211D+04,        MMB22010
     *        -.433090914101065738D+03, .460732511990055463D+03,        MMB22020
     *        -.851274820032285445D+05,-.768523228710292092D+05,        MMB22030
     *         .205560195496649931D+06, .364590804564888233D+07,        MMB22040
     *         .516729324640858565D+07/                                 MMB22050
      DATA Q(1,2),Q(2,2),Q(3,2),Q(4,2),Q(5,2),Q(6,2),Q(7,2),            MMB22060
     *     Q(8,2),Q(9,2)                                                MMB22070
     *           / .225647845389957788D+01, .711593716504863289D+02,    MMB22080
     *        -.141959532218571208D+03,-.609226713172482790D+03,        MMB22090
     *        -.331287173625895019D+05,-.179090061095354737D+06,        MMB22100
     *        -.134472863756067396D+07,-.392640605688557498D+07,        MMB22110
     *        -.134787513180809696D+08/                                 MMB22120
      DATA Q(1,3),Q(2,3),Q(3,3),Q(4,3),Q(5,3),Q(6,3),Q(7,3),            MMB22130
     *     Q(8,3),Q(9,3)                                                MMB22140
     *           /-.939317088650692602D+01, .523806373239931005D+02,    MMB22150
     *        -.121057364184354667D+04,-.258682317593112574D+04,        MMB22160
     *        -.844920979601552086D+05,-.292487285966100229D+06,        MMB22170
     *        -.301034666110382428D+07,-.713311890076026399D+07,        MMB22180
     *        -.387847838647490074D+08/                                 MMB22190
      DATA Q(1,4),Q(2,4),Q(3,4),Q(4,4),Q(5,4),Q(6,4),Q(7,4),            MMB22200
     *     Q(8,4),Q(9,4)                                                MMB22210
     *           /-.437886178371543549D+01, .108803383282850214D+03,    MMB22220
     *        -.396170247546345465D+03, .482799277343449138D+04,        MMB22230
     *        -.260956190444154348D+05, .795561217047195187D+05,        MMB22240
     *        -.121501380828845716D+07,-.560365677885192256D+06,        MMB22250
     *        -.244594499396832557D+08/                                 MMB22260
      DATA Q(1,5),Q(2,5),Q(3,5),Q(4,5),Q(5,5),Q(6,5),Q(7,5),            MMB22270
     *     Q(8,5),Q(9,5)                                                MMB22280
     *           /-.581690214590269560D+01, .121643903355730538D+03,    MMB22290
     *        -.519952651037368415D+03, .634791999150319729D+04,        MMB22300
     *        -.291287737773905651D+05, .162626855091396249D+06,        MMB22310
     *        -.124038130344249191D+07, .121060159492453767D+07,        MMB22320
     *        -.270394900976098049D+08/                                 MMB22330
      DATA Q(1,6),Q(2,6),Q(3,6),Q(4,6),Q(5,6),Q(6,6),Q(7,6),            MMB22340
     *     Q(8,6),Q(9,6)                                                MMB22350
     *           /-.788702380969431594D+01, .130637206116758225D+03,    MMB22360
     *        -.771892722230546718D+03, .715457420628860810D+04,        MMB22370
     *        -.431226753881803730D+05, .199060268477405625D+06,        MMB22380
     *        -.169463010448782541D+07, .196130269020172899D+07,        MMB22390
     *        -.353328964340194695D+08/                                 MMB22400
      DATA Q(1,7),Q(2,7),Q(3,7),Q(4,7),Q(5,7),Q(6,7),Q(7,7),            MMB22410
     *     Q(8,7),Q(9,7)                                                MMB22420
     *           / .359926636077034681D+01, .144952858594130714D+03,    MMB22430
     *         .998518492990023872D+03, .121000936522091880D+05,        MMB22440
     *         .998406746198100760D+05, .637436088274496556D+06,        MMB22450
     *         .488124457435495168D+07, .158873416545813497D+08,        MMB22460
     *         .952711263569843984D+08/                                 MMB22470
C                                                                       MMB22480
C                                  (Z .LT. 0.03125, NU .LT. .05 .OR. (1 MMB22490
C                                    - NU) .LT. 0.004) COEFFICIENTS FOR MMB22500
C                                    POWER SERIES IN NU FOR T1, T2, T3, MMB22510
C                                    RESPECTIVELY. T1, T2, T3 ARE USED  MMB22520
C                                    IN COEFFICIENTS FOR POWER SERIES   MMB22530
C                                    FOR W(Z) = (Z/2)**NU * Y-SUB-NU.   MMB22540
      DATA F1    / .526129484784613088D-05, .313533730646244067D-05,    MMB22550
     *         .468966778832002960D-04, .387395004571051596D-04,        MMB22560
     *        -.433600138305899192D-03, .233918540353778055D-02,        MMB22570
     *        -.840830427176014116D-03,-.210619747888291847D-01,        MMB22580
     *         .747920302081092203D-01,-.603883559512617772D-02,        MMB22590
     *        -.417981294786837086D+00, .668645110914854594D+00,        MMB22600
     *         .638141340718330557D+00,-.143970347809096185D+01,        MMB22610
     *        -.134576433585482691D+00, .318309886183790672D+00,        MMB22620
     *         .258709136016583669D-05/                                 MMB22630
      DATA F2    / .572867669670037890D-06, .229397311396353895D-04,    MMB22640
     *         .209873686010719157D-04, .143717132163052397D-03,        MMB22650
     *        -.510284654371587433D-03, .213837479821310760D-02,        MMB22660
     *         .104955004042030943D-02,-.209403929898152524D-01,        MMB22670
     *         .521523272704097449D-01, .732116630026647386D-01,        MMB22680
     *        -.427655556436905675D+00, .274382993322142905D+00,        MMB22690
     *         .130911899735173940D+01,-.696436028549460850D+00,        MMB22700
     *        -.157079632679489662D+01, .367466905196615962D+00/        MMB22710
      DATA F3    /-.484943071073704969D-05, .969366827458976827D-05,    MMB22720
     *        -.193509783140010526D-04, .387299755087458891D-04,        MMB22730
     *        -.766940410140475115D-04, .156095491967178727D-03,        MMB22740
     *        -.292709443974278280D-03, .669702782215668059D-03,        MMB22750
     *        -.908025366654622614D-03, .355044302302707344D-02,        MMB22760
     *        -.849830319153386275D-04, .236341941622100586D-01,        MMB22770
     *         .259667398807643121D-01, .131092848703934761D+00,        MMB22780
     *         .134576433585482691D+00, .318309886183790672D+00,        MMB22790
     *         .242608782760487999D-05/                                 MMB22800
C                                                                       MMB22810
C                                  1.0 / (FACTORIAL(N))                 MMB22820
      DATA FACT  / .100000000000000000D+01, .100000000000000000D+01,    MMB22830
     *         .500000000000000000D+00, .166666666666666667D+00,        MMB22840
     *         .416666666666666667D-01, .833333333333333333D-02,        MMB22850
     *         .138888888888888889D-02, .198412698412698413D-03,        MMB22860
     *         .248015873015873016D-04, .275573192239858906D-05,        MMB22870
     *         .275573192239858906D-06, .250521083854417187D-07,        MMB22880
     *         .208767569878680990D-08, .160590438368216146D-09,        MMB22890
     *         .114707455977297247D-10, .764716373181981647D-12,        MMB22900
     *         .477947733238738529D-13, .281145725434552076D-14,        MMB22910
     *         .156192069685862264D-15, .822063524662432970D-17,        MMB22920
     *         .411031762331216485D-18, .195729410633912612D-19,        MMB22930
     *         .889679139245057327D-21, .386817017063068404D-22,        MMB22940
     *         .161173757109611835D-23, .644695028438447338D-25,        MMB22950
     *         .247959626322479746D-26, .918368986379554613D-28,        MMB22960
     *         .327988923706983791D-29, .113099628864477169D-30,        MMB22970
     *         .376998762881590564D-32, .121612504155351795D-33,        MMB22980
     *         .380039075485474359D-35, .115163356207719503D-36/        MMB22990
C                                  FIRST EXECUTABLE STATEMENT           MMB23000
      ZSTOP = 20.5D0                                                    MMB23010
      CON1 = 25.0D0                                                     MMB23020
      CON2 = 35.0D0                                                     MMB23030
      IEND = 17                                                         MMB23040
      MEND = 16                                                         MMB23050
      M1 = 11                                                           MMB23060
      M2 = 8                                                            MMB23070
      M3 = 4                                                            MMB23080
      IENDP2 = IEND+2                                                   MMB23090
      MENDM3 = MEND-3                                                   MMB23100
      Z = ARG                                                           MMB23110
      IER = 0                                                           MMB23120
      V = ORDER                                                         MMB23130
      NFLAG = (N.EQ.1)                                                  MMB23140
      IOPT = 2                                                          MMB23150
      IF (NFLAG) IOPT = 1                                               MMB23160
      IF (Z.GE.ZMIN) GO TO 5                                            MMB23170
C                                  ERROR RETURN FOR Z .LT. ZMIN         MMB23180
      IER = 129                                                         MMB23190
    5 IF ((V.GE.0.0D0).AND.(V.LT.1.0D0)) GO TO 10                       MMB23200
C                                                                       MMB23210
C                                  ERROR RETURN FOR V .LT. 0.0          MMB23220
      IER = 129                                                         MMB23230
   10 IF (N.GT.0) GO TO 15                                              MMB23240
C                                  ERROR RETURN FOR N .LE. 0            MMB23250
      IER = 130                                                         MMB23260
   15 IF (IER.EQ.0) GO TO 25                                            MMB23270
      J = N                                                             MMB23280
      IF (IER.EQ.130) J = 1                                             MMB23290
      DO 20 I=1,J                                                       MMB23300
         YN(I) = -XINF                                                  MMB23310
   20 CONTINUE                                                          MMB23320
      GO TO 9000                                                        MMB23330
   25 IF (V.LT.ZMIN) V = 0.0D0                                          MMB23340
      IF (Z.LT..03125D0) GO TO 70                                       MMB23350
      IF (Z.GE.ZSTOP) GO TO 145                                         MMB23360
C                                  IF (1/32 .LE. Z .LT. 20.5), (X**NU)  MMB23370
C                                    * Y-SUB-NU(Z) IS APPROXIMATED BY A MMB23380
C                                    SEQUENCE OF TAYLOR SERIES          MMB23390
C                                    EXPANSIONS, THE FIRST ABOUT A      MMB23400
C                                    GIVEN X, AND THE NEXT ABOUT X +    MMB23410
C                                    DELX. THE SEQUENCE IS CONTINUED    MMB23420
C                                    UNTIL Z IS WITHIN DELX OF THE      MMB23430
C                                    POINT OF EXPANSION.                MMB23440
      J = IDINT((Z+3.5D0)/3.0D0)                                        MMB23450
      X = DBLE(FLOAT(J))*3.0D0-2.0D0                                    MMB23460
      M = 12                                                            MMB23470
      IF (J.GT.1) M = 10                                                MMB23480
C                                  APPROXIMATION AT STARTING POINT X IS MMB23490
C                                    OBTAINED AS A MINIMAX RATIONAL     MMB23500
      GNU = V+1.0D0                                                     MMB23510
      XNUM = P(1,J)*V+P(2,J)                                            MMB23520
      XNUMP = XNUM+P(1,J)                                               MMB23530
      DEN = V+Q(1,J)                                                    MMB23540
      DENP = DEN+1.0D0                                                  MMB23550
      DO 30 I=3,M                                                       MMB23560
         XNUM = XNUM*V+P(I,J)                                           MMB23570
         DEN = DEN*V+Q(I-1,J)                                           MMB23580
         XNUMP = XNUMP*GNU+P(I,J)                                       MMB23590
         DENP = DENP*GNU+Q(I-1,J)                                       MMB23600
   30 CONTINUE                                                          MMB23610
      IF (J.NE.7) GO TO 35                                              MMB23620
      XNUM = XNUM*V+P(11,J)                                             MMB23630
      XNUMP = XNUMP*GNU+P(11,J)                                         MMB23640
   35 DELX = .0625D0                                                    MMB23650
      IF (Z.GT.4.0D0) DELX = DELX+DELX                                  MMB23660
      T = Z-X                                                           MMB23670
      IF (T.LT.0.0D0) DELX = -DELX                                      MMB23680
C                                  EXPAND IN TERMS OF THE AUXILIARY     MMB23690
C                                    FUNCTION W-SUB-NU(Z) = Z**V *      MMB23700
C                                    Y-SUB-NU(Z)                        MMB23710
      M = MEND                                                          MMB23720
      IF (Z.GT.4.0D0) M = MENDM3                                        MMB23730
      XNU = 1.0D0                                                       MMB23740
      IF (X.NE.1.0D0) XNU = (X**V)*X                                    MMB23750
      W(1) = XNUMP/DENP*XNU                                             MMB23760
      W(2) = XNUM/DEN*XNU                                               MMB23770
      XNU = 2.0D0*(GNU-HALF)                                            MMB23780
      IF (Z.GE..5D0) GO TO 40                                           MMB23790
      T = Z+Z-1.5D0                                                     MMB23800
      IF (Z.GE..25D0) GO TO 40                                          MMB23810
      T = 4.0D0*Z-2.0D0                                                 MMB23820
      IF (Z.LT..125D0) T = 8.0D0*Z-2.5D0                                MMB23830
   40 I = IDINT(T/DELX)+1                                               MMB23840
      C2 = 0.0D0                                                        MMB23850
      C2P = 0.0D0                                                       MMB23860
      CUM = W(1)                                                        MMB23870
      CUMP = W(2)                                                       MMB23880
      NMAX = M+3                                                        MMB23890
      DO 65 NN=1,I                                                      MMB23900
         IF (NN.LT.I) GO TO 45                                          MMB23910
         DELX = Z-X                                                     MMB23920
         IF (DELX.EQ.0.0D0) GO TO 65                                    MMB23930
   45    IF (Z.GE.1.0D0) GO TO 50                                       MMB23940
         IF (NN.EQ.9) DELX = DMAX1(DELX,-.03125D0)                      MMB23950
         IF (NN.EQ.17) DELX = DMAX1(DELX,-.015625D0)                    MMB23960
         IF (NN.EQ.25) DELX = DMAX1(DELX,-.0078125D0)                   MMB23970
C                                                                       MMB23980
C                                  CALCULATION OF DERIVATIVES OF W      MMB23990
   50    XIN = 1.0D0/X                                                  MMB24000
         W(3) = XNU*XIN*W(2)-W(1)                                       MMB24010
         XK = 1.0D0                                                     MMB24020
         S = XNU                                                        MMB24030
         DO 55 K=4,NMAX                                                 MMB24040
            S = S-1.0D0                                                 MMB24050
            W(K) = (S*W(K-1)-XK*W(K-3))*XIN-W(K-2)                      MMB24060
            XK = XK+1.0D0                                               MMB24070
   55    CONTINUE                                                       MMB24080
         SUM = FACT(M+2)*W(M+2)*DELX                                    MMB24090
         SUMP = FACT(M+2)*W(M+3)*DELX                                   MMB24100
         J = M+1                                                        MMB24110
         DO 60 K=1,M                                                    MMB24120
            SUM = (SUM+FACT(J)*W(J))*DELX                               MMB24130
            SUMP = (SUMP+FACT(J)*W(J+1))*DELX                           MMB24140
            J = J-1                                                     MMB24150
   60    CONTINUE                                                       MMB24160
C                                  A SUMMING SCHEME SUGGESTED BY KAHAN  MMB24170
C                                    IS USED TO PRESERVE ACCURACY.      MMB24180
C                                    FOOLU IS AN EXTRA VARIABLE TO      MMB24190
C                                    NEGATE POSSIBLE OPTIMIZATION ERROR MMB24200
         C2 = C2+SUM                                                    MMB24210
         T = CUM+C2                                                     MMB24220
         FOOLU = CUM-T                                                  MMB24230
         C2 = C2+FOOLU                                                  MMB24240
         CUM = T                                                        MMB24250
         C2P = C2P+SUMP                                                 MMB24260
         T = CUMP+C2P                                                   MMB24270
         FOOLU = CUMP-T                                                 MMB24280
         C2P = C2P+FOOLU                                                MMB24290
         CUMP = T                                                       MMB24300
         X = X+DELX                                                     MMB24310
         W(1) = CUM                                                     MMB24320
         W(2) = CUMP                                                    MMB24330
   65 CONTINUE                                                          MMB24340
      XNU = (Z**V)*Z                                                    MMB24350
      YN(1) = W(2)/XNU                                                  MMB24360
      IF (NFLAG) GO TO 9005                                             MMB24370
      YN(2) = W(1)/XNU                                                  MMB24380
      GO TO 170                                                         MMB24390
C                                  FOR (Z .LT. .03125, .05 .LE. NU .LE. MMB24400
C                                    .996), Y-SUB-NU = (J-SUB-NU *      MMB24410
C                                    COS(NU * PI) - J-SUB-(-NU) ) /     MMB24420
C                                    SIN(NU * PI)                       MMB24430
   70 ZH = Z*0.5D0                                                      MMB24440
      X = ZH*ZH                                                         MMB24450
      IF ((V.LT.GNU1).OR.(V.GT.GNU2)) GO TO 95                          MMB24460
      X = -X                                                            MMB24470
      GNU = V                                                           MMB24480
      VPI = V/PI2                                                       MMB24490
      IF (V.GT.0.5D0) VPI = (1.0D0-V)/PI2                               MMB24500
      VPI = VPI+VPI                                                     MMB24510
      VSIN = DSIN(VPI)                                                  MMB24520
      VCOS = DCOS(VPI)                                                  MMB24530
      IF (V.GT.0.5D0) VCOS = -VCOS                                      MMB24540
      T1 = ZH**GNU                                                      MMB24550
      T2 = 1.0D0/T1                                                     MMB24560
      DO 90 I=1,2                                                       MMB24570
         XC = T1                                                        MMB24580
         DO 85 J=1,2                                                    MMB24590
            T = GNU-0.5D0                                               MMB24600
            IF (J.EQ.2) T = T+2.0D0                                     MMB24610
            IF (T.GT.0.5D0) T = T-1.0D0                                 MMB24620
C                                  EVALUATION OF GAMMA(GNU + 1.0)       MMB24630
            XNUM = GP(1)*T                                              MMB24640
            DEN = T                                                     MMB24650
            DO 75 L=2,8                                                 MMB24660
               XNUM = (XNUM+GP(L))*T                                    MMB24670
               DEN = (DEN+GQ(L))*T                                      MMB24680
   75       CONTINUE                                                    MMB24690
            DEN = DEN+GQ(1)                                             MMB24700
            GAM = CO+XNUM/DEN                                           MMB24710
            IF (J.LT.I) GAM = GAM*GNU                                   MMB24720
            IF (J.EQ.2) GAM = GAM/(GNU+1.0D0)                           MMB24730
            IF (I*J.EQ.4) GAM = GAM/(GNU+2.0D0)                         MMB24740
C                                                                       MMB24750
C                                  EVALUATION OF J-SUB-GNU              MMB24760
            SUM = FACT(5)                                               MMB24770
            XK = 4.0D0                                                  MMB24780
            DO 80 K=1,4                                                 MMB24790
               SUM = SUM*X/(GNU+XK)+FACT(5-K)                           MMB24800
               XK = XK-1.0D0                                            MMB24810
   80       CONTINUE                                                    MMB24820
            SCRAT(J) = XC*SUM/GAM                                       MMB24830
            GNU = -GNU                                                  MMB24840
            XC = T2                                                     MMB24850
   85    CONTINUE                                                       MMB24860
         YN(I) = (VCOS*SCRAT(1)-SCRAT(2))/VSIN                          MMB24870
         IF (NFLAG) GO TO 9005                                          MMB24880
         IF (I.EQ.2) GO TO 170                                          MMB24890
         VCOS = -VCOS                                                   MMB24900
         VSIN = -VSIN                                                   MMB24910
         GNU = V+1.0D0                                                  MMB24920
         T1 = T1*ZH                                                     MMB24930
         T2 = T2/ZH                                                     MMB24940
   90 CONTINUE                                                          MMB24950
      GO TO 170                                                         MMB24960
C                                  FOR (Z .LT. .03125, NU .LT. .05 OR   MMB24970
C                                    NU .GT. .996) USE FULLERTON        MMB24980
C                                    EXPANSIONS OF W-SUB-NU(Z) = (Z /   MMB24990
C                                    2) ** V * Y-SUB-NU(Z)              MMB25000
   95 GNU = V                                                           MMB25010
      IF (V.GT.GNU2) GNU = (V-HALF)-0.5D0                               MMB25020
      T1 = F1(17)                                                       MMB25030
      T2 = 0.0D0                                                        MMB25040
      T3(1) = F3(17)                                                    MMB25050
      DO 100 I=1,16                                                     MMB25060
         T1 = T1*GNU+F1(I)                                              MMB25070
         T2 = T2*GNU+F2(I)                                              MMB25080
         T3(1) = T3(1)*GNU+F3(I)                                        MMB25090
  100 CONTINUE                                                          MMB25100
      T1 = T1/(1.0D0-GNU)                                               MMB25110
      T2 = T2/(1.0D0-GNU*GNU)                                           MMB25120
      T3(1) = T3(1)/(1.0D0+GNU)                                         MMB25130
      DLNX = 2.0D0*DLOG(ZH)                                             MMB25140
      DLNXNU = GNU*DLNX                                                 MMB25150
      IF (DLNXNU.LT.DLNHAF) GO TO 110                                   MMB25160
C                                  IF X ** V .GT. 0.5 COMPUTE (X ** V   MMB25170
C                                    -1.0) / V USING POWER SERIES FOR E MMB25180
C                                    ** (V * LOG(X))                    MMB25190
      XNU = FACT(IENDP2)                                                MMB25200
      DO 105 I=1,IEND                                                   MMB25210
         XNU = (XNU*DLNXNU)+FACT(IENDP2-I)                              MMB25220
  105 CONTINUE                                                          MMB25230
      XNU = DLNX*XNU                                                    MMB25240
      SCRAT(1) = T2+XNU*T1                                              MMB25250
      XNU = GNU*XNU+1.0D0                                               MMB25260
      GO TO 115                                                         MMB25270
  110 XNU = ZH**(GNU+GNU)                                               MMB25280
      SCRAT(1) = (XNU*T1-T3(1))/GNU                                     MMB25290
  115 XK = 1.0D0                                                        MMB25300
C                                  SCRAT CONTAINS COEFFICIENTS FOR      MMB25310
C                                    EXPANSION W = SUM (SCRAT(K+1) *    MMB25320
C                                    (-X) ** K), K = 0 TO 4             MMB25330
      DO 120 I=2,6                                                      MMB25340
         T = XK-GNU                                                     MMB25350
         T3(I) = T3(I-1)/(T*XK)                                         MMB25360
         SCRAT(I) = (SCRAT(I-1)*(1.0D0-GNU/XK)/T-2.0D0*T3(I))/(XK+GNU)  MMB25370
         XK = XK+1.0D0                                                  MMB25380
  120 CONTINUE                                                          MMB25390
      XK = 5.0D0+GNU                                                    MMB25400
      SUM = SCRAT(5)                                                    MMB25410
      SUMP = XK*SCRAT(6)+T3(6)                                          MMB25420
      IF (X.LT.1.0D-30) X = 0.0D0                                       MMB25430
      DO 125 I=1,4                                                      MMB25440
         K = 5-I                                                        MMB25450
         XK = XK-1.0D0                                                  MMB25460
         SUMP = (-SUMP)*X+XK*SCRAT(K+1)+T3(K+1)                         MMB25470
         SUM = SUM*(-X)+SCRAT(K)                                        MMB25480
  125 CONTINUE                                                          MMB25490
      SUMP = XNU*T1/ZH-SUMP*ZH                                          MMB25500
      RTXNU = DSQRT(XNU)                                                MMB25510
      XNU = GNU/ZH                                                      MMB25520
      SUM = SUM/RTXNU                                                   MMB25530
      IF (GNU.LT.0.0D0) GO TO 130                                       MMB25540
      YN(1) = SUM                                                       MMB25550
      IF (NFLAG) GO TO 9005                                             MMB25560
      YN(2) = XNU*YN(1)-SUMP/RTXNU                                      MMB25570
      GO TO 170                                                         MMB25580
C                                  IF GNU = 1.0 - V, CALCULATE          MMB25590
C                                    Y(Z,GNU+1)                         MMB25600
  130 YN(1) = XNU*SUM-SUMP/RTXNU                                        MMB25610
      IF (NFLAG) GO TO 9005                                             MMB25620
      XNU = XNU+2.0D0/Z                                                 MMB25630
C                                  CHECK FOR POTENTIAL OVERFLOW IN      MMB25640
C                                    COMPUTING SEQUENCE OF VALUES       MMB25650
      IF (XNU.GE.-XINF/YN(1)) GO TO 135                                 MMB25660
      YN(2) = XNU*YN(1)-SUM                                             MMB25670
      GO TO 170                                                         MMB25680
  135 DO 140 I=2,N                                                      MMB25690
         YN(I) = -XINF                                                  MMB25700
  140 CONTINUE                                                          MMB25710
      IER = 131                                                         MMB25720
      GO TO 9000                                                        MMB25730
C                                  USE AN ASYMPTOTIC SERIES FOR (Z .GE. MMB25740
C                                    ZSTOP)                             MMB25750
  145 IF (Z.LE.ZMAX) GO TO 155                                          MMB25760
      IER = 132                                                         MMB25770
      DO 150 I=1,N                                                      MMB25780
         YN(I) = 0.0D0                                                  MMB25790
  150 CONTINUE                                                          MMB25800
      GO TO 9000                                                        MMB25810
  155 XC = RTPI2/DSQRT(Z)                                               MMB25820
      XIN = .015625D0/(Z*Z)                                             MMB25830
      M = 16                                                            MMB25840
      IF (Z.GE.CON1) M = M1                                             MMB25850
      IF (Z.GE.CON2) M = M2                                             MMB25860
      IF (Z.GE.130.0D0) M = M3                                          MMB25870
      XM = 4.0D0*DBLE(FLOAT(M))                                         MMB25880
C                                  REDUCTION OF ARGUMENT FOR SIN AND    MMB25890
C                                    COS ROUTINES                       MMB25900
      T = DINT(Z/TWOPI)                                                 MMB25910
      IF (T.LE.ZMAX2) GO TO 158                                         MMB25920
C                                  T .GT. 16**8                         MMB25930
      X = (((Z-T*TWOPIA)-T*TWOPIB)-T*TWOPI2)-(V+0.5D0)/PI2              MMB25940
      GO TO 159                                                         MMB25950
C                                  T .LE. 16**8                         MMB25960
  158 X = ((Z-T*TWOPI1)-T*TWOPI2)-(V+0.5D0)/PI2                         MMB25970
  159 VSIN = DSIN(X)                                                    MMB25980
      VCOS = DCOS(X)                                                    MMB25990
      GNU = V+V                                                         MMB26000
      DO 165 I=1,2                                                      MMB26010
         S = ((XM-1.0D0)-GNU)*((XM-1.0D0)+GNU)*XIN*0.5D0                MMB26020
         T = (GNU-(XM-3.0D0))*(GNU+(XM-3.0D0))                          MMB26030
         XNUM = S*T*FACT(2*M+1)                                         MMB26040
         T1 = (GNU-(XM+1.0D0))*(GNU+(XM+1.0D0))                         MMB26050
         DEN = S*T1*FACT(2*M+2)                                         MMB26060
         XK = XM                                                        MMB26070
         K = M+M                                                        MMB26080
         T1 = T                                                         MMB26090
         DO 160 J=2,M                                                   MMB26100
            XK = XK-4.0D0                                               MMB26110
            S = ((XK-1.0D0)-GNU)*((XK-1.0D0)+GNU)                       MMB26120
            T = (GNU-(XK-3.0D0))*(GNU+(XK-3.0D0))                       MMB26130
            XNUM = (XNUM+FACT(K-1))*S*T*XIN                             MMB26140
            DEN = (DEN+FACT(K))*S*T1*XIN                                MMB26150
            K = K-2                                                     MMB26160
            T1 = T                                                      MMB26170
  160    CONTINUE                                                       MMB26180
         XNUM = XNUM+1.0D0                                              MMB26190
         DEN = (DEN+1.0D0)*(GNU*GNU-1.0D0)/(8.0D0*Z)                    MMB26200
         YN(I) = XC*(XNUM*VSIN+DEN*VCOS)                                MMB26210
         IF (NFLAG) GO TO 9005                                          MMB26220
         T = VSIN                                                       MMB26230
         VSIN = -VCOS                                                   MMB26240
         VCOS = T                                                       MMB26250
         GNU = GNU+2.0D0                                                MMB26260
  165 CONTINUE                                                          MMB26270
C                                  IF N .GT. 2, COMPUTE Y(Z,ORDER+I) I  MMB26280
C                                    = 2, N-1                           MMB26290
  170 IF (N.LE.2) GO TO 9005                                            MMB26300
      GNU = V+V+2.0D0                                                   MMB26310
      NN = N                                                            MMB26320
      DO 185 J=3,NN                                                     MMB26330
         IF (DABS(YN(J-1)).LT.1.0D0) GO TO 180                          MMB26340
C                                  CHECK FOR POTENTIAL OVERFLOW IN      MMB26350
C                                    COMPUTING SEQUENCE OF VALUES       MMB26360
C                                                                       MMB26370
         IF (DABS(GNU/Z).LT.DABS(XINF/YN(J-1))) GO TO 180               MMB26380
         DO 175 L=J,N                                                   MMB26390
            YN(L) = -XINF                                               MMB26400
  175    CONTINUE                                                       MMB26410
         IER = 131                                                      MMB26420
         GO TO 9000                                                     MMB26430
  180    YN(J) = GNU*YN(J-1)/Z-YN(J-2)                                  MMB26440
         GNU = GNU+2.0D0                                                MMB26450
  185 CONTINUE                                                          MMB26460
      GO TO 9005                                                        MMB26470
 9000 CONTINUE                                                          MMB26480
      CALL UERTST (IER,6HMMBSYN)                                        MMB26490
 9005 RETURN                                                            MMB26500
      END                                                               MMB26510
C   IMSL ROUTINE NAME   - UERTST                                        UERT0010
C                                                                       UERT0020
C-----------------------------------------------------------------------UERT0030
C                                                                       UERT0040
C   COMPUTER            - VAX/SINGLE                                    UERT0050
C                                                                       UERT0060
C   LATEST REVISION     - JUNE 1, 1982                                  UERT0070
C                                                                       UERT0080
C   PURPOSE             - PRINT A MESSAGE REFLECTING AN ERROR CONDITION UERT0090
C                                                                       UERT0100
C   USAGE               - CALL UERTST (IER,NAME)                        UERT0110
C                                                                       UERT0120
C   ARGUMENTS    IER    - ERROR PARAMETER. (INPUT)                      UERT0130
C                           IER = I+J WHERE                             UERT0140
C                             I = 128 IMPLIES TERMINAL ERROR MESSAGE,   UERT0150
C                             I =  64 IMPLIES WARNING WITH FIX MESSAGE, UERT0160
C                             I =  32 IMPLIES WARNING MESSAGE.          UERT0170
C                             J = ERROR CODE RELEVANT TO CALLING        UERT0180
C                                 ROUTINE.                              UERT0190
C                NAME   - A CHARACTER STRING OF LENGTH SIX PROVIDING    UERT0200
C                           THE NAME OF THE CALLING ROUTINE. (INPUT)    UERT0210
C                                                                       UERT0220
C   PRECISION/HARDWARE  - SINGLE/ALL                                    UERT0230
C                                                                       UERT0240
C   REQD. IMSL ROUTINES - UGETIO,USPKD                                  UERT0250
C                                                                       UERT0260
C   NOTATION            - INFORMATION ON SPECIAL NOTATION AND           UERT0270
C                           CONVENTIONS IS AVAILABLE IN THE MANUAL      UERT0280
C                           INTRODUCTION OR THROUGH IMSL ROUTINE UHELP  UERT0290
C                                                                       UERT0300
C   REMARKS      THE ERROR MESSAGE PRODUCED BY UERTST IS WRITTEN        UERT0310
C                TO THE STANDARD OUTPUT UNIT. THE OUTPUT UNIT           UERT0320
C                NUMBER CAN BE DETERMINED BY CALLING UGETIO AS          UERT0330
C                FOLLOWS..   CALL UGETIO(1,NIN,NOUT).                   UERT0340
C                THE OUTPUT UNIT NUMBER CAN BE CHANGED BY CALLING       UERT0350
C                UGETIO AS FOLLOWS..                                    UERT0360
C                                NIN = 0                                UERT0370
C                                NOUT = NEW OUTPUT UNIT NUMBER          UERT0380
C                                CALL UGETIO(3,NIN,NOUT)                UERT0390
C                SEE THE UGETIO DOCUMENT FOR MORE DETAILS.              UERT0400
C                                                                       UERT0410
C   COPYRIGHT           - 1982 BY IMSL, INC. ALL RIGHTS RESERVED.       UERT0420
C                                                                       UERT0430
C   WARRANTY            - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN UERT0440
C                           APPLIED TO THIS CODE. NO OTHER WARRANTY,    UERT0450
C                           EXPRESSED OR IMPLIED, IS APPLICABLE.        UERT0460
C                                                                       UERT0470
C-----------------------------------------------------------------------UERT0480
C                                                                       UERT0490
      SUBROUTINE UERTST (IER,NAME)                                      UERT0500
C                                  SPECIFICATIONS FOR ARGUMENTS         UERT0510
      INTEGER            IER                                            UERT0520
      INTEGER            NAME(1)                                        UERT0530
C                                  SPECIFICATIONS FOR LOCAL VARIABLES   UERT0540
      INTEGER            I,IEQ,IEQDF,IOUNIT,LEVEL,LEVOLD,NAMEQ(6),      UERT0550
     *                   NAMSET(6),NAMUPK(6),NIN,NMTB                   UERT0560
      DATA               NAMSET/1HU,1HE,1HR,1HS,1HE,1HT/                UERT0570
      DATA               NAMEQ/6*1H /                                   UERT0580
      DATA               LEVEL/4/,IEQDF/0/,IEQ/1H=/                     UERT0590
C                                  UNPACK NAME INTO NAMUPK              UERT0600
C                                  FIRST EXECUTABLE STATEMENT           UERT0610
      CALL USPKD (NAME,6,NAMUPK,NMTB)                                   UERT0620
C                                  GET OUTPUT UNIT NUMBER               UERT0630
      CALL UGETIO(1,NIN,IOUNIT)                                         UERT0640
C                                  CHECK IER                            UERT0650
      IF (IER.GT.999) GO TO 25                                          UERT0660
      IF (IER.LT.-32) GO TO 55                                          UERT0670
      IF (IER.LE.128) GO TO 5                                           UERT0680
      IF (LEVEL.LT.1) GO TO 30                                          UERT0690
C                                  PRINT TERMINAL MESSAGE               UERT0700
      IF (IEQDF.EQ.1) WRITE(IOUNIT,35) IER,NAMEQ,IEQ,NAMUPK             UERT0710
      IF (IEQDF.EQ.0) WRITE(IOUNIT,35) IER,NAMUPK                       UERT0720
      GO TO 30                                                          UERT0730
    5 IF (IER.LE.64) GO TO 10                                           UERT0740
      IF (LEVEL.LT.2) GO TO 30                                          UERT0750
C                                  PRINT WARNING WITH FIX MESSAGE       UERT0760
      IF (IEQDF.EQ.1) WRITE(IOUNIT,40) IER,NAMEQ,IEQ,NAMUPK             UERT0770
      IF (IEQDF.EQ.0) WRITE(IOUNIT,40) IER,NAMUPK                       UERT0780
      GO TO 30                                                          UERT0790
   10 IF (IER.LE.32) GO TO 15                                           UERT0800
C                                  PRINT WARNING MESSAGE                UERT0810
      IF (LEVEL.LT.3) GO TO 30                                          UERT0820
      IF (IEQDF.EQ.1) WRITE(IOUNIT,45) IER,NAMEQ,IEQ,NAMUPK             UERT0830
      IF (IEQDF.EQ.0) WRITE(IOUNIT,45) IER,NAMUPK                       UERT0840
      GO TO 30                                                          UERT0850
   15 CONTINUE                                                          UERT0860
C                                  CHECK FOR UERSET CALL                UERT0870
      DO 20 I=1,6                                                       UERT0880
         IF (NAMUPK(I).NE.NAMSET(I)) GO TO 25                           UERT0890
   20 CONTINUE                                                          UERT0900
      LEVOLD = LEVEL                                                    UERT0910
      LEVEL = IER                                                       UERT0920
      IER = LEVOLD                                                      UERT0930
      IF (LEVEL.LT.0) LEVEL = 4                                         UERT0940
      IF (LEVEL.GT.4) LEVEL = 4                                         UERT0950
      GO TO 30                                                          UERT0960
   25 CONTINUE                                                          UERT0970
      IF (LEVEL.LT.4) GO TO 30                                          UERT0980
C                                  PRINT NON-DEFINED MESSAGE            UERT0990
      IF (IEQDF.EQ.1) WRITE(IOUNIT,50) IER,NAMEQ,IEQ,NAMUPK             UERT1000
      IF (IEQDF.EQ.0) WRITE(IOUNIT,50) IER,NAMUPK                       UERT1010
   30 IEQDF = 0                                                         UERT1020
      RETURN                                                            UERT1030
   35 FORMAT(19H *** TERMINAL ERROR,10X,7H(IER = ,I3,                   UERT1040
     1       20H) FROM IMSL ROUTINE ,6A1,A1,6A1)                        UERT1050
   40 FORMAT(27H *** WARNING WITH FIX ERROR,2X,7H(IER = ,I3,            UERT1060
     1       20H) FROM IMSL ROUTINE ,6A1,A1,6A1)                        UERT1070
   45 FORMAT(18H *** WARNING ERROR,11X,7H(IER = ,I3,                    UERT1080
     1       20H) FROM IMSL ROUTINE ,6A1,A1,6A1)                        UERT1090
   50 FORMAT(20H *** UNDEFINED ERROR,9X,7H(IER = ,I5,                   UERT1100
     1       20H) FROM IMSL ROUTINE ,6A1,A1,6A1)                        UERT1110
C                                                                       UERT1120
C                                  SAVE P FOR P = R CASE                UERT1130
C                                    P IS THE PAGE NAMUPK               UERT1140
C                                    R IS THE ROUTINE NAMUPK            UERT1150
   55 IEQDF = 1                                                         UERT1160
      DO 60 I=1,6                                                       UERT1170
   60 NAMEQ(I) = NAMUPK(I)                                              UERT1180
   65 RETURN                                                            UERT1190
      END                                                               UERT1200
C   IMSL ROUTINE NAME   - UGETIO                                        UGET0010
C                                                                       UGET0020
C-----------------------------------------------------------------------UGET0030
C                                                                       UGET0040
C   COMPUTER            - VAX/SINGLE                                    UGET0050
C                                                                       UGET0060
C   LATEST REVISION     - JUNE 1, 1981                                  UGET0070
C                                                                       UGET0080
C   PURPOSE             - TO RETRIEVE CURRENT VALUES AND TO SET NEW     UGET0090
C                           VALUES FOR INPUT AND OUTPUT UNIT            UGET0100
C                           IDENTIFIERS.                                UGET0110
C                                                                       UGET0120
C   USAGE               - CALL UGETIO(IOPT,NIN,NOUT)                    UGET0130
C                                                                       UGET0140
C   ARGUMENTS    IOPT   - OPTION PARAMETER. (INPUT)                     UGET0150
C                           IF IOPT=1, THE CURRENT INPUT AND OUTPUT     UGET0160
C                           UNIT IDENTIFIER VALUES ARE RETURNED IN NIN  UGET0170
C                           AND NOUT, RESPECTIVELY.                     UGET0180
C                           IF IOPT=2, THE INTERNAL VALUE OF NIN IS     UGET0190
C                           RESET FOR SUBSEQUENT USE.                   UGET0200
C                           IF IOPT=3, THE INTERNAL VALUE OF NOUT IS    UGET0210
C                           RESET FOR SUBSEQUENT USE.                   UGET0220
C                NIN    - INPUT UNIT IDENTIFIER.                        UGET0230
C                           OUTPUT IF IOPT=1, INPUT IF IOPT=2.          UGET0240
C                NOUT   - OUTPUT UNIT IDENTIFIER.                       UGET0250
C                           OUTPUT IF IOPT=1, INPUT IF IOPT=3.          UGET0260
C                                                                       UGET0270
C   PRECISION/HARDWARE  - SINGLE/ALL                                    UGET0280
C                                                                       UGET0290
C   REQD. IMSL ROUTINES - NONE REQUIRED                                 UGET0300
C                                                                       UGET0310
C   NOTATION            - INFORMATION ON SPECIAL NOTATION AND           UGET0320
C                           CONVENTIONS IS AVAILABLE IN THE MANUAL      UGET0330
C                           INTRODUCTION OR THROUGH IMSL ROUTINE UHELP  UGET0340
C                                                                       UGET0350
C   REMARKS      EACH IMSL ROUTINE THAT PERFORMS INPUT AND/OR OUTPUT    UGET0360
C                OPERATIONS CALLS UGETIO TO OBTAIN THE CURRENT UNIT     UGET0370
C                IDENTIFIER VALUES. IF UGETIO IS CALLED WITH IOPT=2 OR  UGET0380
C                IOPT=3, NEW UNIT IDENTIFIER VALUES ARE ESTABLISHED.    UGET0390
C                SUBSEQUENT INPUT/OUTPUT IS PERFORMED ON THE NEW UNITS. UGET0400
C                                                                       UGET0410
C   COPYRIGHT           - 1978 BY IMSL, INC. ALL RIGHTS RESERVED.       UGET0420
C                                                                       UGET0430
C   WARRANTY            - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN UGET0440
C                           APPLIED TO THIS CODE. NO OTHER WARRANTY,    UGET0450
C                           EXPRESSED OR IMPLIED, IS APPLICABLE.        UGET0460
C                                                                       UGET0470
C-----------------------------------------------------------------------UGET0480
C                                                                       UGET0490
      SWBROUTINE UGETIO(IOPT,NIN,NOUT)                                  UGET0500
C                                  SPECIFICATIONS FOR ARGUMENTS         UGET0510
      INTEGER            IOPT,NIN,NOUT                                  UGET0520
C                                  SPECIFICATIONS FOR LOCAL VARIABLES   UGET0530
      INTEGER            NIND,NOUTD                                     UGET0540
      DATA               NIND/1/,NOUTD/2/                               UGET0550
C                                  FIRST EXECUTABLE STATEMENT           UGET0560
      IF (IOPT.EQ.3) GO TO 10                                           UGET0570
      IF (IOPT.EQ.2) GO TO 5                                            UGET0580
      IF (IOPT.NE.1) GO TO 9005                                         UGET0590
      NIN = NIND                                                        UGET0600
      NOUT = NOUTD                                                      UGET0610
      GO TO 9005                                                        UGET0620
    5 NIND = NIN                                                        UGET0630
      GO TO 9005                                                        UGET0640
   10 NOUTD = NOUT                                                      UGET0650
 9005 RETURN                                                            UGET0660
      END                                                               UGET0670
C   IMSL ROUTINE NAME   - USPKD                                         USPK0010
C                                                                       USPK0020
C-----------------------------------------------------------------------USPK0030
C                                                                       USPK0040
C   COMPUTER            - VAX/SINGLE                                    USPK0050
C                                                                       USPK0060
C   LATEST REVISION     - JUNE 1, 1982                                  USPK0070
C                                                                       USPK0080
C   PURPOSE             - NUCLEUS CALLED BY IMSL ROUTINES THAT HAVE     USPK0090
C                           CHARACTER STRING ARGUMENTS                  USPK0100
C                                                                       USPK0110
C   USAGE               - CALL USPKD  (PACKED,NCHARS,UNPAKD,NCHMTB)     USPK0120
C                                                                       USPK0130
C   ARGUMENTS    PACKED - CHARACTER STRING TO BE UNPACKED.(INPUT)       USPK0140
C                NCHARS - LENGTH OF PACKED. (INPUT)  SEE REMARKS.       USPK0150
C                UNPAKD - INTEGER ARRAY TO RECEIVE THE UNPACKED         USPK0160
C                         REPRESENTATION OF THE STRING. (OUTPUT)        USPK0170
C                NCHMTB - NCHARS MINUS TRAILING BLANKS. (OUTPUT)        USPK0180
C                                                                       USPK0190
C   PRECISION/HARDWARE  - SINGLE/ALL                                    USPK0200
C                                                                       USPK0210
C   REQD. IMSL ROUTINES - NONE                                          USPK0220
C                                                                       USPK0230
C   REMARKS  1.  USPKD UNPACKS A CHARACTER STRING INTO AN INTEGER ARRAY USPK0240
C                IN (A1) FORMAT.                                        USPK0250
C            2.  UP TO 129 CHARACTERS MAY BE USED.  ANY IN EXCESS OF    USPK0260
C                THAT ARE IGNORED.                                      USPK0270
C                                                                       USPK0280
C   COPYRIGHT           - 1982 BY IMSL, INC. ALL RIGHTS RESERVED.       USPK0290
C                                                                       USPK0300
C   WARRANTY            - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN USPK0310
C                           APPLIED TO THIS CODE. NO OTHER WARRANTY,    USPK0320
C                           EXPRESSED OR IMPLIED, IS APPLICABLE.        USPK0330
C                                                                       USPK0340
C-----------------------------------------------------------------------USPK0350
      SUBROUTINE USPKD  (PACKED,NCHARS,UNPAKD,NCHMTB)                   USPK0360
C                                  SPECIFICATIONS FOR ARGUMENTS         USPK0370
      INTEGER            NC,NCHARS,NCHMTB                               USPK0380
C                                                                       USPK0390
      INTEGER            UNPAKD(1),IBLANK                               USPK0400
      INTEGER            PACKED(1)                                      USPK0410
      DATA               IBLANK /1H /                                   USPK0420
C                                  INITIALIZE NCHMTB                    USPK0430
      NCHMTB = 0                                                        USPK0440
C                                  RETURN IF NCHARS IS LE ZERO          USPK0450
      IF(NCHARS.LE.0) RETURN                                            USPK0460
C                                  SET NC=NUMBER OF CHARS TO BE DECODED USPK0470
      NC = MIN0 (129,NCHARS)                                            USPK0480
cdel      DECODE (NC,150,PACKED) (UNPAKD(I),I=1,NC)                         USPK0490
150   FORMAT (129A1)                                                    USPK0500
C                                  CHECK UNPAKD ARRAY AND SET NCHMTB    USPK0510
C                                  BASED ON TRAILING BLANKS FOUND       USPK0520
      DO 200 N = 1,NC                                                   USPK0530
         NN = NC - N + 1                                                USPK0540
         IF(UNPAKD(NN) .NE. IBLANK) GO TO 210                           USPK0550
  200 CONTINUE                                                          USPK0560
  210 NCHMTB = NN                                                       USPK0570
      RETURN                                                            USPK0580
      END                                                               USPK0590



