      REAL FUNCTION MMBSJ1 (ARG,IER)                                    
C   IMSL ROUTINE NAME   - MMBSJ1                                                
C                                                                               
C-----------------------------------------------------------------------        
C                                                                               
C   COMPUTER            - CRAY/SINGLE                                           
C                                                                               
C   LATEST REVISION     - AUGUST 1, 1981                                        
C                                                                               
C   PURPOSE             - BESSEL FUNCTION OF THE FIRST KIND OF ORDER            
C                           ONE                                                 
C                                                                               
C   USAGE               - FUNCTION MMBSJ1 (ARG,IER)                             
C                                                                               
C   ARGUMENTS    MMBSJ1 - OUTPUT VALUE OF THE FUNCTION AT ARG. MMBSJ1           
C                           MUST BE TYPED APPROPRIATELY IN THE CALLING          
C                           PROGRAM. (SEE THE PRECISION/HARDWARE                
C                           SECTION.)                                           
C                ARG    - INPUT ARGUMENT. THE ABSOLUTE VALUE OF ARG             
C                           MUST BE LESS THAN OR EQUAL TO XMAX, WHICH           
C                           IS OF THE ORDER OF 10**8. THE EXACT VALUE OF        
C                           XMAX MAY ALLOW LARGER RANGES FOR ARG ON SOME        
C                           COMPUTERS. SEE THE PROGRAMMING NOTES IN THE         
C                           MANUAL FOR THE EXACT VALUES.                        
C                IER    - ERROR PARAMETER. (OUTPUT)                             
C                         TERMINAL ERROR                                        
C                           IER = 129 INDICATES THAT THE ABSOLUTE VALUE         
C                             OF ARG IS GREATER THAN XMAX. MMBSJ1 IS            
C                             SET TO ZERO.                                      
C                                                                               
C   PRECISION/HARDWARE  - DOUBLE/H32,H36                                        
C                       - SINGLE/H48,H60                                        
C                                                                               
C   REQD. IMSL ROUTINES - UERTST,UGETIO                                         
C                                                                               
C   NOTATION            - INFORMATION ON SPECIAL NOTATION AND                   
C                           CONVENTIONS IS AVAILABLE IN THE MANUAL              
C                           INTRODUCTION OR THROUGH IMSL ROUTINE UHELP          
C                                                                               
C   COPYRIGHT           - 1978 BY IMSL, INC. ALL RIGHTS RESERVED.               
C                                                                               
C   WARRANTY            - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN         
C                           APPLIED TO THIS CODE. NO OTHER WARRANTY,            
C                           EXPRESSED OR IMPLIED, IS APPLICABLE.                
C                                                                               
C-----------------------------------------------------------------------        
C                                  SPECIFICATIONS FOR ARGUMENTS                 
      IMPLICIT REAL (A-H,O-Z),
     +INTEGER (I-N)
      INTEGER            IER                                            
      REAL               ARG                                            
C                                  SPECIFICATIONS FOR LOCAL VARIABLES           
      REAL               AP1,AP2,AQ1,AQ2,AX,B,D,DEN,DEN2,P,Q,PTPI2,     
     1                   TWOPI1,TWOPI2,V,XC,XMAX,XMAX1,XMIN,XNUM,XNUM2, 
     2                   XSMALL,X01,X02,X11,X12,Y,Z,ZSQ                 
      DIMENSION          P(6),Q(4),B(7),D(6),AP1(5),AQ1(4),AP2(5),AQ2(4)
C                                  MACHINE DEPENDENT CONSTANTS                  
C                                     RTPI2 = SQRT(2/PI)                        
C                                     TWOPI1 + TWOPI2 = 2 * PI TO EXTRA         
C                                     PRECISION                                 
C                                     XMAX = 2**27, LARGEST ACCEPTABLE          
C                                     ARGUMENT                                  
C                                     XMAX1 = SMALLEST FLOATING-POINT           
C                                     CONSTANT WITH ENTIRELY INTEGER            
C                                     REPRESENTATION                            
C                                     XMIN = 2**(-927), ARGUMENT,               
C                                     SCALED BY XMAX1, BELOW WHICH J1           
C                                     IS ZERO                                   
C                                     XSMALL = 2**(-25), ARGUMENT BELOW         
C                                     WHICH MMBSJ0 MAY BE REPRESENTED           
C                                     BY ONE TERM IN THE ASCENDING              
C                                     SERIES.                                   
C                                     X01 + X02 = FIRST ZERO OF J-SUB-1         
C                                     TO EXTRA PRECISION                        
C                                     X11 + X12 = SECOND ZERO OF                
C                                     J-SUB-1 TO EXTRA PRECISION                
      DATA RTPI2/ 0.79788456080286E+00/,XMAX1/ 0.14073748835533E+15/,   
     1  TWOPI1/ 0.62831840515137E+01/,TWOPI2/ 0.12556659146019E-05/,    
     2  XMAX/ 0.13421772800000E+09/,XMIN/ 0.88144256634025E-279/,       
     3  X01/ 0.38317059702075E+01/,X02/ 0.28441110477497E-14/,          
     4  X11/ 0.70155866698156E+01/,X12/ 0.14116713058668E-13/,          
     5  XSMALL/ 0.29802322387695E-07/,IFCN/13/
C    6                              ,FCN/5HBESJ1/                      
C                                  COEFFICIENTS FOR RATIONAL                    
C                                    APPROXIMATION OF J-1(X) / (X *             
C                                    (X**2 - X0**2)), XSMALL .LT.               
C                                    ABS(X) .LE. 4.0                            
      DATA P/-0.10081930399287E+06, 0.61080730489571E+07,               
     1       -0.13374567015858E+09, 0.78532593732672E+03,               
     2       -0.30082193935614E+01, 0.47266558846037E-02/               
      DATA Q/ 0.22991844276313E+06, 0.44064062008235E+08,               
     1        0.39273000055590E+10, 0.67993320116315E+03/               
C                                  COEFFICIENTS FOR RATIONAL                    
C                                    APPROXIMATION OF J-1(X) / (X *             
C                                    (X**2 - X1**2)), 4.0 .LT. ABS(X)           
C                                    .LE. 8.0 NUMERATOR IN MINI-NEWTON          
C                                    FORM                                       
      DATA B/-0.12883354862941E+01, 0.13497183972739E+04,               
     &  -0.54811322816122E+06, 0.10236844351296E+09,
     &  -0.77226352881628E+10,  
     &  0.59588084601724E+11, 0.49023966813227E+12/                     
      DATA D/ 0.49261697794519E+06, 0.18994940353412E+09,               
     1        0.51928735318537E+11, 0.92617171760798E+13,               
     2        0.82402998615392E+15, 0.89947181683819E+03/               
C                                  COEFFICIENTS FOR HART                        
C                                    APPROXIMATION, ABS(X) .GT. 8.0             
      DATA AP1/ 0.17431379748379E+03, 0.31327529563551E+04,             
     &  0.13090420511035E+05, 0.12909184718962E+05,
     &  0.12285053764359E+01/  
      DATA AQ1/ 0.16904721775009E+03, 0.31092814167700E+04,             
     1          0.13066783087844E+05, 0.12909184718962E+05/             
      DATA AP2/ 0.37994453796981E+01, 0.51736532818366E+02,             
     &  0.17442916890924E+03, 0.14465282874995E+03,
     &  0.36363466476035E-01/  
      DATA AQ2/ 0.85223920643413E+02, 0.11191098527048E+04,             
     1          0.37343401060163E+04, 0.30859270133323E+04/             
C                                  FIRST EXECUTABLE STATEMENT                   
      IER = 0                                                           
      AX = ABS(ARG)                                                     
      IF (AX.GT.XMAX) GO TO 40                                          
      IF (AX*XMAX1.GE.XMIN) GO TO 5                                     
      MMBSJ1 = 0.0E0                                                    
      GO TO 9005                                                        
    5 IF (AX.GT.XSMALL) GO TO 10                                        
      MMBSJ1 = ARG/2.0E0                                                
      GO TO 9005                                                        
   10 IF (AX.GT.8.0E0) GO TO 30                                         
      Y = AX*AX                                                         
      IF (AX.GT.4.0E0) GO TO 20                                         
C                                  XSMALL .LT. ABS(X) .LE. 4.0                  
      XNUM = (P(6)*Y+P(5))*Y+P(4)                                       
      DEN = Y+Q(4)                                                      
      DO 15 I=1,3                                                       
         XNUM = XNUM*Y+P(I)                                             
         DEN = DEN*Y+Q(I)                                               
   15 CONTINUE                                                          
      Z = (AX-X01)-X02                                                  
      MMBSJ1 = (XNUM/DEN)*ARG*Z*(AX+X01)                                
      GO TO 9005                                                        
C                                  4.0 .LT. ABS(X) .LE. 8.0                     
   20 XNUM = B(1)                                                       
      DEN = (Y+D(6))*Y+D(1)                                             
      DO 25 I=2,5                                                       
         XNUM = XNUM*Y+B(I)                                             
         DEN = DEN*Y+D(I)                                               
   25 CONTINUE                                                          
      XNUM = XNUM*(AX-8.0E0)*(AX+8.0E0)+B(6)                            
      XNUM = XNUM*(AX-4.0E0)*(AX+4.0E0)+B(7)                            
      Z = (AX-X11)-X12                                                  
      MMBSJ1 = (XNUM/DEN)*ARG*Z*(AX+X11)                                
      GO TO 9005                                                        
C                                  ABS(X) .GT. 8.0                              
   30 XC = RTPI2/SQRT(AX)                                               
      IF (ARG.LT.0.0E0) XC = -XC                                        
      Z = 8.0E0/AX                                                      
      ZSQ = Z*Z                                                         
      V = ((AX/TWOPI1+XMAX1)-XMAX1)+0.375E0                             
      V = (AX-V*TWOPI1)-V*TWOPI2                                        
      XNUM = AP1(5)                                                     
      DEN = 1.0E0                                                       
      XNUM2 = AP2(5)                                                    
      DEN2 = 1.0E0                                                      
      DO 35 I=1,4                                                       
         XNUM = XNUM*ZSQ+AP1(I)                                         
         DEN = DEN*ZSQ+AQ1(I)                                           
         XNUM2 = XNUM2*ZSQ+AP2(I)                                       
         DEN2 = DEN2*ZSQ+AQ2(I)                                         
   35 CONTINUE                                                          
      MMBSJ1 = XC*((XNUM/DEN)*COS(V)-Z*(XNUM2/DEN2)*SIN(V))             
      GO TO 9005                                                        
C                                  ERROR RETURN FOR ABS(X) .GT. XMAX            
   40 MMBSJ1 = 0.0E0                                                    
      IER = 129                                                         
 9000 CONTINUE                                                          
      CALL UERTST (IER,'MMBSJ1')                                        
 9005 RETURN                                                            
      END                                                               
