
C	--------------------------------------------------------------
C	subroutines BESSRT and BESSRA
C	--------------------------------------------------------------

C       sven ivansson 1990





	SUBROUTINE BESSRT (LUNREAD,LUNBESS, BESSFIL,KAQ,XMAXQ,NXQ,XGLUEQ,
     $	YMAXQ,NYQ,YGLUEQ)
C	to read a table with J0-values and Y0REG-values created with 
C	besseltab.for

	IMPLICIT NONE

	INTEGER KADIM, NXDIM,NTOTDIM
	PARAMETER ( KADIM = 100, NXDIM = 16000,NTOTDIM = 100000 )

	COMMON /BESST/ KAEVEN,KAODD,AK(kadim)
	COMMON /BESST_X/ XMAX,NX,DDX,XGLUE, 
     $                   YMAX,NY,DDY,YGLUE, NX2, BESS(0:ntotdim-1,2)
	INTEGER KAEVEN,KAODD, NX,NY, NX2
	REAL*8 AK, XMAX,DDX,XGLUE,YMAX,DDY,YGLUE
	COMPLEX*16 BESS
	SAVE   /BESST/, /BESST_X/

	INTEGER LUNREAD,LUNBESS
	INTEGER KAQ, NXQ,NYQ
	REAL*8 XMAXQ,XGLUEQ,YMAXQ,YGLUEQ
	CHARACTER*32 BESSFIL, GETSTRING

	INTEGER KA, K, IX,IY, IPOS, ISLASK
	COMPLEX*16 BVAL(0:nxdim+1), Y0VAL(0:nxdim+1)

	IF (LUNBESS.GT.0) THEN

	  IF (LUNREAD.EQ.5) WRITE(*,*) 'Give BESSFIL (starting with "%" ',
     $	  'for "ka=0 asymptotics") :'
	  BESSFIL = ' '
	  DO WHILE (BESSFIL.EQ.' ')
	    READ(LUNREAD,'(A32)') BESSFIL
	  END DO
	  IPOS = 1
	  BESSFIL = GETSTRING (BESSFIL,IPOS)

	  IF (BESSFIL(1:1).NE.'%') THEN

	    OPEN(LUNBESS,FILE=BESSFIL,FORM='UNFORMATTED',STATUS='OLD')

	    READ(LUNBESS) KA
	    IF (KA.GT.kadim) STOP 313
	    IF (KA.GT.0) READ(LUNBESS) (AK(K),K=1,KA)

	    READ(LUNBESS) XMAX,NX,XGLUE
	    READ(LUNBESS) YMAX,NY,YGLUE

	    IF (XMAX.GE.0.0D0) THEN
	      IF (NX.LT.2) STOP 312       ! 4-point x-interpolation
	      IF (NX.GT.nxdim) STOP 313
	      DDX = DBLE(XMAX)/NX
	      IF (YMAX.GT.0.0) THEN
	        IF (NY.LT.2) STOP 314     ! 4-point y-interpolation
	        DDY = DBLE(YMAX)/NY
	      ELSE
	        IF (YMAX.NE.0.0D0 .OR. NY.NE.-1) STOP 415
	      END IF
	      NX2 = NX+2
	      IF ( (NX2*(NY+2)) .GT. ntotdim ) STOP 315
	      DO IY = 0, NY+1
	        READ(LUNBESS) (BVAL(IX),IX=0,NX+1)
	        READ(LUNBESS) (Y0VAL(IX),IX=0,NX+1)
	        ISLASK = IY*NX2
	        DO IX = 0, NX+1
	          BESS(ISLASK+IX,1) = BVAL(IX)
	        END DO
	        DO IX = 0, NX+1
	          BESS(ISLASK+IX,2) = Y0VAL(IX)
	        END DO
	      END DO
	    ELSE
	      IF (YMAX.GE.0.0D0) STOP 324
	    END IF

	    CLOSE(LUNBESS)

	  ELSE

	    KA = 0
	    XMAX = -1.0D0
	    YMAX = -1.0D0

	  END IF

	ELSE IF (LUNBESS.EQ.-1) THEN

	  BESSFIL = '-'
	  READ(-LUNBESS) KA
	  IF (KA.GT.0) READ(-LUNBESS) (AK(K),K=1,KA)
	  XMAX = -1.0D0
	  YMAX = -1.0D0

	ELSE IF (LUNBESS.EQ.0) THEN

	  BESSFIL = '%'
	  KA = 0
	  XMAX = -1.0D0
	  YMAX = -1.0D0

	ELSE

	  STOP 1452

	END IF

	KAEVEN = 2*(KA/2)
	IF (KAEVEN.EQ.KA) THEN
	  KAODD = KAEVEN-1
	ELSE
	  KAODD = KAEVEN+1
	END IF

	KAQ = KA
	XMAXQ = XMAX
	NXQ = NX
	XGLUEQ = XGLUE
	YMAXQ = YMAX
	NYQ = NY
	YGLUEQ = YGLUE

	RETURN
	END







	COMPLEX*16 FUNCTION BESSRA (IBH,CDLOGZBY2,Z)
C	to compute J0(Z) %% IBH=0 %% or H0(1)(Z) %% IBH=1 or IBH=-1 %% by 
C	reading table or using asymptotic expression
C	when IBH=-1 CDLOGZBY2 is used as the logarithm of 0.5*Z
C	for H0(1)(Z) a certain symmetry is demanded when Z=(0.0D0,0.0D0)
C	note: bessrt must have been called to update /besst/

	IMPLICIT NONE

	INTEGER KADIM, NXDIM,NTOTDIM
	PARAMETER ( KADIM = 100, NXDIM = 16000,NTOTDIM = 100000 )

	COMMON /BESST/ KAEVEN,KAODD,AK(kadim)
	COMMON /BESST_X/ XMAX,NX,DDX,XGLUE, 
     $                   YMAX,NY,DDY,YGLUE, NX2, BESS(0:ntotdim-1,2)
	INTEGER KAEVEN,KAODD, NX,NY, NX2
	REAL*8 AK, XMAX,DDX,XGLUE,YMAX,DDY,YGLUE
	COMPLEX*16 BESS
	SAVE   /BESST/, /BESST_X/

	INTEGER IBH
	COMPLEX*16 CDLOGZBY2,Z

	INTEGER IX,IY, IM1, I0,I1,I2, KAEVENLOK,KAODDLOK, K
	REAL*8 PIG4,AFAKT,TWOBYPI
	REAL*8 XP,YP,Y,QX,XR,QY,YR, B1,B2,B3,B4, A1,A2,A3,A4, SLASK
	COMPLEX*16 ZNOW, BHLP,Y0REG, PZ,QZ, CHI, ZNOW2INV, EIMULT, 
     $	BS(4),YS(4)
	LOGICAL KVADR23, HKVADR3

	DATA PIG4,AFAKT,TWOBYPI /0.7853981633974483D0,
     $	0.7978845608028654D0,0.6366197723675814D0/  ! pi/4, dsqrt(2/pi), 2/pi

	ZNOW = Z

	IF (ZNOW.EQ.(0.0D0,0.0D0)) THEN
	  IF (IBH.EQ.0) THEN
	    BESSRA = (1.0D0,0.0D0)
	  ELSE
C	    "overflow" but choose symmetry  " J0 = 0.5d0 * ( H0(1) + H0(2) ) "
	    BESSRA = DCMPLX(1.0D0,0.0D0)
C	    ( would have been  " dcmplx(1.0d0,0.5772156649015329d0*twobypi) " 
C	    if the logarithmic term had been neglected )
	  END IF
	  RETURN
	END IF

	XP = DREAL(ZNOW)
	IF (XP.LT.0.0D0) THEN
	  KVADR23 = .TRUE.
	  XP = -XP
	  Y = -DIMAG(ZNOW)
	ELSE
	  KVADR23 = .FALSE.
	  Y = DIMAG(ZNOW)
	END IF
	YP = DABS(Y)

	IF (XP.LE.XMAX .AND. YP.LE.YMAX) THEN
C	  grid interpolation in this case
	  QX = XP/DDX
	  IX = QX
	  XR = QX-IX
	  IF (IX.EQ.0) THEN
	    IX = 1
	    XR = XR-1.0
	  ELSE IF (IX.EQ.NX) THEN
	    IX = NX-1
	    XR = XR+1.0
	  END IF
	  IF (YP.GT.0.0D0) THEN
	    QY = YP/DDY
	    IY = QY
	    YR = QY-IY
	  ELSE
	    IY = 0
	    YR = 0.0D0
	  END IF
	  IF (YR.GT.0.0D0) THEN
	    IF (IY.EQ.0) THEN
	      IY = 1
	      YR = YR-1.0
	    ELSE IF (IY.EQ.NY) THEN
	      IY = NY-1
	      YR = YR+1.0
	    END IF
	    SLASK = (YR-1)*(YR-2)
	    B1 = -SLASK*YR/6
	    B2 = SLASK*(YR+1)/2
	    SLASK = YR*(YR+1)
	    B3 = -SLASK*(YR-2)/2
	    B4 = SLASK*(YR-1)/6
C	    Lagrange four-point from Abramovich & Stegun 25.2.13, cf. lagr4.for
	    IM1 = (IY-1)*NX2+IX
	    I0 = IM1+NX2
	    I1 = I0+NX2
	    I2 = I1+NX2
	    BS(1) = B1*BESS(IM1-1,1) + B2*BESS(I0-1,1) + 
     $	    B3*BESS(I1-1,1) + B4*BESS(I2-1,1)     
	    BS(2) = B1*BESS(IM1,1) + B2*BESS(I0,1) + 
     $	    B3*BESS(I1,1) + B4*BESS(I2,1)     
	    BS(3) = B1*BESS(IM1+1,1) + B2*BESS(I0+1,1) + 
     $	    B3*BESS(I1+1,1) + B4*BESS(I2+1,1)     
	    BS(4) = B1*BESS(IM1+2,1) + B2*BESS(I0+2,1) + 
     $	    B3*BESS(I1+2,1) + B4*BESS(I2+2,1)     
	    IF (IBH.NE.0) THEN
	      YS(1) = B1*BESS(IM1-1,2) + B2*BESS(I0-1,2) + 
     $	      B3*BESS(I1-1,2) + B4*BESS(I2-1,2)     
	      YS(2) = B1*BESS(IM1,2) + B2*BESS(I0,2) + 
     $	      B3*BESS(I1,2) + B4*BESS(I2,2)     
	      YS(3) = B1*BESS(IM1+1,2) + B2*BESS(I0+1,2) + 
     $	      B3*BESS(I1+1,2) + B4*BESS(I2+1,2)     
	      YS(4) = B1*BESS(IM1+2,2) + B2*BESS(I0+2,2) + 
     $	      B3*BESS(I1+2,2) + B4*BESS(I2+2,2)     
	    END IF
	  ELSE
	    I0 = IY*NX2+IX
	    BS(1) = BESS(I0-1,1)
	    BS(2) = BESS(I0,1)
	    BS(3) = BESS(I0+1,1)
	    BS(4) = BESS(I0+2,1)
	    IF (IBH.NE.0) THEN
	      YS(1) = BESS(I0-1,2)
	      YS(2) = BESS(I0,2)
	      YS(3) = BESS(I0+1,2)
	      YS(4) = BESS(I0+2,2)
	    END IF
	  END IF
	  SLASK = (XR-1)*(XR-2)
	  A1 = -SLASK*XR/6
	  A2 = SLASK*(XR+1)/2
	  SLASK = XR*(XR+1)
	  A3 = -SLASK*(XR-2)/2
	  A4 = SLASK*(XR-1)/6
C	  Lagrange four-point from Abramovich & Stegun 25.2.13, cf. lagr4.for
	  IF (IBH.EQ.0) THEN
	    BESSRA = A1*BS(1) + A2*BS(2) + A3*BS(3) + A4*BS(4)  
	    IF (Y.LT.0.0D0) BESSRA = DCONJG(BESSRA)
	  ELSE
	    BHLP = A1*BS(1) + A2*BS(2) + A3*BS(3) + A4*BS(4)  
	    Y0REG = A1*YS(1) + A2*YS(2) + A3*YS(3) + A4*YS(4)  
	    IF (Y.LT.0.0D0) THEN
	      BHLP = DCONJG(BHLP)
	      Y0REG = DCONJG(Y0REG)
	    END IF
	    IF (IBH.EQ.1) THEN
	      BESSRA = EIMULT(Y0REG) + BHLP*(1.0+DCMPLX(0.0D0,TWOBYPI)*
     $	      CDLOG(0.5D0*ZNOW))
	    ELSE
	      BESSRA = EIMULT(Y0REG) + BHLP*(1.0+DCMPLX(0.0D0,TWOBYPI)*
     $	      CDLOGZBY2)
	    END IF
	  END IF
	ELSE
C	  asymptotics in this case
	  K = 2 * ( 1 + IDINT(
     $	  DMAX1(DABS(DREAL(ZNOW)),DABS(DIMAG(ZNOW)))) ) !!! from "cdabs(znow)"
	  KAEVENLOK = MIN0 ( KAEVEN, K )     ! suggested by "besbod cashan"
	  KAODDLOK = MIN0 ( KAODD, 1+K )     ! suggested by "besbod cashan"
	  IF (IBH.EQ.0) THEN
	    IF (KVADR23) ZNOW = -ZNOW     ! will give better approximation
	  ELSE IF (KVADR23) THEN
	    IF (Y.GT.0.0D0) THEN
	      HKVADR3 = .TRUE.      ! note, the sign of y was changed
	      ZNOW = -ZNOW    ! will give better approximation, but "add 2*J0"
	    ELSE
	      HKVADR3 = .FALSE.
	    END IF
	  ELSE
	    HKVADR3 = .FALSE.
	  END IF
	  IF (KAEVENLOK.GT.0) THEN
	    ZNOW2INV = 1.0D0/(ZNOW*ZNOW)
	    PZ = AK(KAEVENLOK)
	    DO K = KAEVENLOK-2, 2, -2
	      PZ = PZ*ZNOW2INV + AK(K)
	    END DO
	    PZ = PZ*ZNOW2INV + (1.0D0,0.0D0)
	  ELSE
	    PZ = (1.0D0,0.0D0)
	  END IF
	  CHI = ZNOW-PIG4
	  IF (KAODDLOK.GT.0) THEN
	    QZ = AK(KAODDLOK)
	    DO K = KAODDLOK-2, 1, -2
	      QZ = QZ*ZNOW2INV + AK(K)
	    END DO
	    QZ = QZ/ZNOW
	    IF (IBH.EQ.0) THEN
	      BESSRA = PZ*CDCOS(CHI) - QZ*CDSIN(CHI)
	    ELSE IF (HKVADR3) THEN
C	      for the "first-quadrant znow" :   " H0(1) + 2*J0 " 
	      BESSRA = (3*PZ+EIMULT(QZ))*CDCOS(CHI) + 
     $	      (-3*QZ+EIMULT(PZ))*CDSIN(CHI)
	    ELSE
	      BESSRA = ( PZ + EIMULT(QZ) ) * CDEXP(EIMULT(CHI))
	    END IF
	  ELSE IF (IBH.EQ.0) THEN
	    BESSRA = PZ*CDCOS(CHI)
	  ELSE IF (HKVADR3) THEN
C	    for the "first-quadrant znow" :   " H0(1) + 2*J0 " 
	    BESSRA = PZ * ( 3*CDCOS(CHI) + EIMULT(CDSIN(CHI)) )
	  ELSE
	    BESSRA = PZ * CDEXP(EIMULT(CHI))
	  END IF
	  BESSRA = AFAKT * BESSRA / CDSQRT(ZNOW)
	END IF

	RETURN
	END
