C	to compute certain Airy-function values, using the SLATEC routines 
C	zairy,zbiry
C       note 1: the SLATEC routines daie,daide,dbie,dbide would be faster 
C               for real arguments, however I have found dbie and dbide 
C               to be unreliable for nonzero arguments between -1 and +1, 
C               hence the "ccc commented" statements; I have found the 
C               corresponding IMSL routines reliable, but they are not 
C               freely available
C	note 2: below the four pertinent functions are computed "separately", 
C	        it would be more efficient to lift out certain "common" 
C	        computations

C	Sven Ivansson  1990



	SUBROUTINE AIRY4E (Z, ALN,BLN, AIE,AIDE,BIE,BIDE)
C	upon exit: 
C	  Ai(z)  = dexp(-aln) * aie
C	  Ai'(z) = dexp(-aln) * aide
C	  Bi(z)  = dexp(+bln) * bie
C	  Bi'(z) = dexp(+bln) * bide

	IMPLICIT NONE

	COMPLEX*16 Z

	REAL*8 ALN,BLN
	COMPLEX*16 AIE,AIDE,BIE,BIDE

	INTEGER NZ,IERR
	REAL*8 X, DAIE,DAIDE,DBIE,DBIDE, 
     $	REZ,IMZ, REUT,IMUT
	COMPLEX*16 ZNOW, CZET,CSCALE

	ZNOW = Z

ccc	IF (DIMAG(ZNOW).EQ.0.0D0) THEN

ccc	  X = DREAL(ZNOW)
ccc	  IF (X.GT.0.0D0) THEN
ccc	    ALN = 0.6666666666666667D0 * X * DSQRT(X)
ccc	  ELSE
ccc	    ALN = 0.0D0
ccc	  END IF
ccc	  BLN = ALN
ccc	  AIE = DAIE(X)
ccc	  AIDE = DAIDE(X)
ccc	  BIE = DBIE(X)
ccc	  BIDE = DBIDE(X)

ccc	ELSE

	  REZ = DREAL(ZNOW)
	  IMZ = DIMAG(ZNOW)

	  CZET = 0.6666666666666667D0 * ZNOW*CDSQRT(ZNOW)
	  CSCALE = CDEXP ( -DCMPLX(0.0D0,DIMAG(CZET)) )
	  ALN = DREAL(CZET)
	  BLN = DABS(ALN)

	  CALL ZAIRY (REZ,IMZ,0,2, REUT,IMUT,NZ,IERR)
	  AIE = DCMPLX(REUT,IMUT)
	  IF (IERR.NE.0) THEN
	    IF (IERR.NE.3) THEN
	      WRITE(*,*) '%%% zairy 0, ierr=',IERR
	      WRITE(*,*) 'znow=',ZNOW
	      STOP 567
	    END IF
	  END IF
	  AIE = CSCALE * AIE
	  CALL ZAIRY (REZ,IMZ,1,2, REUT,IMUT,NZ,IERR)
	  AIDE = DCMPLX(REUT,IMUT)
	  IF (IERR.NE.0) THEN
	    IF (IERR.NE.3) THEN
	      WRITE(*,*) '%%% zairy 1, ierr=',IERR
	      WRITE(*,*) 'znow=',ZNOW
	      STOP 567
	    END IF
	  END IF
	  AIDE = CSCALE * AIDE
	  CALL ZBIRY (REZ,IMZ,0,2, REUT,IMUT,IERR)
	  BIE = DCMPLX(REUT,IMUT)
	  IF (IERR.NE.0) THEN
	    IF (IERR.NE.3) THEN
	      WRITE(*,*) '%%% zbiry 0, ierr=',IERR
	      WRITE(*,*) 'znow=',ZNOW
	      STOP 567
	    END IF
	  END IF
	  CALL ZBIRY (REZ,IMZ,1,2, REUT,IMUT,IERR)
	  BIDE = DCMPLX(REUT,IMUT)
	  IF (IERR.NE.0) THEN
	    IF (IERR.NE.3) THEN
	      WRITE(*,*) '%%% zbiry 1, ierr=',IERR
	      WRITE(*,*) 'znow=',ZNOW
	      STOP 567
	    END IF
	    IF (IERR.NE.3) STOP 567
	  END IF

ccc	END IF

	RETURN
	END