      REAL FUNCTION MMBSJ0 (ARG,IER)                                    
C   IMSL ROUTINE NAME   - MMBSJ0                                                
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                           ZERO                                                
C                                                                               
C   USAGE               - FUNCTION MMBSJ0 (ARG,IER)                             
C                                                                               
C   ARGUMENTS    MMBSJ0 - OUTPUT VALUE OF THE FUNCTION AT ARG. MMBSJ0           
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 MUST        
C                           BE LESS THAN OR EQUAL TO XMAX, WHICH IS OF          
C                           THE ORDER OF 10**8. THE EXACT VALUE OF XMAX         
C                           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. MMBSJ0 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               AX,CONST,DOWN,PI2,PROD,PRP,PRQ,P0,P1,          
     1                   QRP,QRQ,Q0,Q1,R0,R1,TWOPI1,TWOPI2,U,UP,W,      
     2                   XDEN,XMAX,XMAX1,XNUM,XSMALL,X01,X02,X11,X12,   
     3                   Z,ZSQ                                          
      DIMENSION          PRP(6),PRQ(7),P0(5),P1(5),QRP(4),QRQ(6),       
     1                   Q0(4),Q1(4)                                    
C                                                                               
C                                  MACHINE DEPENDENT CONSTANTS                  
C                                     PI2 = 2 / PI                              
C                                     TWOPI1 + TWOPI2 = 2 * PI TO               
C                                     EXTRA PRECISION                           
C                                     XMAX = 2**27, LARGEST ACCEPTABLE          
C                                     ARGUMENT                                  
C                                     XMAX1 = SMALLEST FLOATING-POINT           
C                                     CONSTANT WITH ENTIRELY INTEGER            
C                                     REPRESENTATION                            
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-0         
C                                     TO EXTRA PRECISION                        
C                                     X11 + X12 = SECOND ZERO OF                
C                                     J-SUB-0 TO EXTRA PRECISION                
C                                                                               
      DATA PI2/ 0.63661977236758E+00/,XMAX1/ 0.14073748835533E+15/,     
     1  XMAX/ 0.13421772800000E+09/,XSMALL/ 0.29802322387695E-07/,      
     2   X01/ 0.24048255576958E+01/,X02/ 0.14093185550049E-13/,         
     3   X11/ 0.55200781102863E+01/,X12/ 0.45217780699621E-14/,         
     4  TWOPI1/ 0.62831840515137E+01/,TWOPI2/ 0.12556659146019E-05/     
C                                                                               
C                                  COEFFICIENTS FOR RATIONAL                    
C                                    APPROXIMATION OF J-0(X) / (X**2 -          
C                                    X0**2), XSMALL .LT. ABS(X) .LE.            
C                                    4.0                                        
C                                                                               
      DATA PRP/-0.56872211842749E+06, 0.26129260233700E+08,             
     &  -0.40458140629737E+09, 0.47924620550636E-01,
     &  -0.25854452494242E+02, 0.55660507153209E+04/
      DATA QRP/ 0.17203745045066E+06, 0.29250600138949E+08,             
     1  0.23397695097682E+10, 0.58207145592288E+03/                     
C                                  COEFFICIENTS FOR RATIONAL                    
C                                    APPROXIMATION OF J-0(X) / (X**2 -          
C                                    X1**2), 4.0 .LT. ABS(X) .LE. 8.0           
C                                                                               
      DATA PRQ/-0.35889459418734E+03,-0.37096308273188E+03,             
     &  0.20461990476581E+03, 0.40620030105152E+03,
     &  0.61636412010778E+02,  
     & -0.11924550016636E+02,-0.11039583656827E+03/                     
      DATA QRQ/ 0.15583306932129E+03,-0.91601874750017E+03,             
     &  0.36801141942428E+04,-0.94667237302538E+04,
     &  0.12039506966788E+05,-0.17372426501636E+02/
C                                  COEFFICIENTS FOR HART                        
C                                    APPROXIMATION, ABS(X) .GT. 8.0             
C                                                                               
      DATA P0/ 0.22045010439652E+04, 0.88944375329606E+04,              
     &  0.85548225415067E+04, 0.90047934748029E+00,
     &  0.12867758574871E+03/  
      DATA Q0/ 0.22140488519147E+04, 0.89038361417096E+04,              
     1  0.85548225415067E+04, 0.13088490049992E+03/                     
      DATA P1/-0.13990976865961E+02,-0.46093826814625E+02,              
     &  -0.37510534954957E+02,-0.93525953294032E-02,
     &  -0.10497327982346E+01/  
      DATA Q1/ 0.92156697552653E+03, 0.29719837452085E+04,              
     1  0.24006742371173E+04, 0.74428389741411E+02/                     
C                                  FIRST EXECUTABLE STATEMENT                   
      IER = 0                                                           
      AX = ABS(ARG)                                                     
      IF (AX.GT.XMAX) GO TO 35                                          
      IF (AX.GT.XSMALL) GO TO 5                                         
      MMBSJ0 = 1.0E0                                                    
      GO TO 9005                                                        
    5 IF (AX.GT.8.0E0) GO TO 25                                         
      IF (AX.GT.4.0E0) GO TO 15                                         
C                                  XSMALL .LT. ABS(ARG) .LE. 4.0                
      ZSQ = AX*AX                                                       
      XNUM = (PRP(4)*ZSQ+PRP(5))*ZSQ+PRP(6)                             
      XDEN = ZSQ+QRP(4)                                                 
      DO 10 I=1,3                                                       
         XNUM = XNUM*ZSQ+PRP(I)                                         
         XDEN = XDEN*ZSQ+QRP(I)                                         
   10 CONTINUE                                                          
C                                  CALCULATION TO PRESERVE ACCURACY             
C                                    NEAR THE FIRST ZERO OF J-0                 
      PROD = ((AX-X01)-X02)*(AX+X01)                                    
      MMBSJ0 = PROD*XNUM/XDEN                                           
      GO TO 9005                                                        
C                                  4.0 .LT. ABS(ARG) .LE. 8.0                   
   15 ZSQ = 1.0E0-(AX*AX)/64.0E0                                        
      XNUM = PRQ(6)*ZSQ+PRQ(7)                                          
      XDEN = ZSQ+QRQ(6)                                                 
      DO 20 I=1,5                                                       
         XNUM = XNUM*ZSQ+PRQ(I)                                         
         XDEN = XDEN*ZSQ+QRQ(I)                                         
   20 CONTINUE                                                          
C                                  CALCULATION TO PRESERVE ACCURACY             
C                                    NEAR THE SECOND ZERO OF J-0                
      PROD = (AX+X11)*((AX-X11)-X12)                                    
      MMBSJ0 = PROD*XNUM/XDEN                                           
      GO TO 9005                                                        
C                                  ABS(ARG) .GT. 8.0                            
   25 Z = 8.0E0/AX                                                      
      W = AX/TWOPI1                                                     
      W = ((W+XMAX1)-XMAX1)+0.125E0                                     
      U = (AX-W*TWOPI1)-W*TWOPI2                                        
      ZSQ = Z*Z                                                         
      XNUM = P0(4)*ZSQ+P0(5)                                            
      XDEN = ZSQ+Q0(4)                                                  
      UP = P1(4)*ZSQ+P1(5)                                              
      DOWN = ZSQ+Q1(4)                                                  
      DO 30 I=1,3                                                       
         XNUM = XNUM*ZSQ+P0(I)                                          
         XDEN = XDEN*ZSQ+Q0(I)                                          
         UP = UP*ZSQ+P1(I)                                              
         DOWN = DOWN*ZSQ+Q1(I)                                          
   30 CONTINUE                                                          
      R0 = XNUM/XDEN                                                    
      R1 = UP/DOWN                                                      
      MMBSJ0 = SQRT(PI2/AX)*(R0*COS(U)-Z*R1*SIN(U))                     
      GO TO 9005                                                        
C                                  ERROR RETURN FOR ABS(ARG) .GT. XMAX          
   35 MMBSJ0 = 0.0E0                                                    
      IER = 129                                                         
 9000 CONTINUE                                                          
      CALL UERTST (IER,'MMBSJ0')                                        
 9005 RETURN                                                            
      END                                                               
