! copyright (c) 2008 Swedish Defence Research Agency (FOI)
! may not be copied or distributed without written permission from FOI



C	to compute certain Airy-function values, using the gil et al. 
C	routines aiz,biz
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, 2006



	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*32 Z

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

	INTEGER IERR, IARG3, ISLASK
	REAL*16 X, DAIE,DAIDE,DBIE,DBIDE, 
     $	REZ,IMZ, REUT,IMUT, SLASK
	COMPLEX*32 ZNOW, CZET,CSCALE, QCMPLX
	INTEGER IARG3UT
	LOGICAL EXTRA

	EXTRA = .FALSE.
	GOTO 100

	ENTRY AIRY4E_IARG3 (Z, ALN,BLN, AIE,AIDE,BIE,BIDE, IARG3UT)
	EXTRA = .TRUE.

 100	ZNOW = Z

ccc	IF (AIMAG(ZNOW).EQ.0.0Q0) THEN

ccc	  X = REAL(ZNOW)
ccc	  IF (X.GT.0.0Q0) THEN
ccc	    ALN = (2.0Q0/3.0Q0) * X * SQRT(X)
ccc	  ELSE
ccc	    ALN = 0.0Q0
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 = REAL(ZNOW)
	  IMZ = AIMAG(ZNOW)

	  CZET = (2.0Q0/3.0Q0) * ZNOW*SQRT(ZNOW)
	  CSCALE = EXP ( -QCMPLX(0.0Q0,AIMAG(CZET)) )
	  ALN = REAL(CZET)
	  BLN = ABS(ALN)

	  CALL AIZ ( 1,2, REZ,IMZ, REUT,IMUT, IERR )
	  AIE = QCMPLX(REUT,IMUT)
	  IF (IERR.NE.0) THEN
	    WRITE(*,*) '%%% aiz 1, ierr=',IERR
	    WRITE(*,*) 'znow=',ZNOW
	    STOP 567
	  END IF
	  AIE = CSCALE * AIE
	  CALL AIZ ( 2,2, REZ,IMZ, REUT,IMUT, IERR )
	  AIDE = QCMPLX(REUT,IMUT)
	  IF (IERR.NE.0) THEN
	    WRITE(*,*) '%%% aiz 2, ierr=',IERR
	    WRITE(*,*) 'znow=',ZNOW
	    STOP 567
	  END IF
	  AIDE = CSCALE * AIDE
	  CALL BIZ_IARG3 ( 1,2, REZ,IMZ, REUT,IMUT, IERR, IARG3 )
	  BIE = QCMPLX(REUT,IMUT)
C         scaling correction of bie remains, see below
	  IF (IERR.NE.0) THEN
	    WRITE(*,*) '%%% biz 1, ierr=',IERR
	    WRITE(*,*) 'znow=',ZNOW
	    STOP 567
	  END IF
	  CALL BIZ_IARG3 ( 2,2, REZ,IMZ, REUT,IMUT, IERR, ISLASK )
	  IF (ISLASK.NE.IARG3) STOP '%%% airy4e, biz iarg3 inconsist'
	  BIDE = QCMPLX(REUT,IMUT)
	  IF (IERR.NE.0) THEN
	    WRITE(*,*) '%%% biz 2, ierr=',IERR
	    WRITE(*,*) 'znow=',ZNOW
	    STOP 567
	  END IF
	  IF (IARG3.EQ.2) THEN
	    BIE = CSCALE * BIE
	    BIDE = CSCALE * BIDE
	    IF (ALN.GT.0.0Q0) THEN
C	      can only happen numerically
	      SLASK = EXP(-2*ALN)
	      BIE = SLASK*BIE	    
	      BIDE = SLASK*BIDE	    
	    END IF
	  ELSE IF (IARG3.EQ.1) THEN
	    BIE = CONJG(CSCALE) * BIE
	    BIDE = CONJG(CSCALE) * BIDE
	    IF (ALN.LT.0.0Q0) THEN
C	      can only happen numerically
	      SLASK = EXP(2*ALN)
	      BIE = SLASK*BIE	    
	      BIDE = SLASK*BIDE	    
	    END IF
	  ELSE
	    STOP '%%%, airy4e, biz iarg3'
	  END IF

ccc	END IF

	IF (EXTRA) IARG3UT = IARG3

	RETURN
	END
