C       sven ivansson 1990




	SUBROUTINE J0ASCOMP (Z, BFAKT)
C	to compute J0(Z) = BFAKT(1) + BFAKT(2)
C	using asymptotic expression
C	so that bfakt(1) is the term with dcos(x-pig4) and
C	        bfakt(2) is the term with dsin(x-pig4)
C	where x=Re(z)
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)           ! abbreviated
	INTEGER KAEVEN,KAODD
	REAL*8 AK

	COMPLEX*16 Z
	COMPLEX*16 BFAKT(2)

	INTEGER KAEVENLOK,KAODDLOK, K
	REAL*8 XCHI, PIG4,AFAKT,TWOBYPI
	COMPLEX*16 ZNOW, PZ,QZ, ZNOW2INV, EIY, CAFAKTZ, CSLASK, 
     $	CEIY,SEIY

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

	ZNOW=Z

	K = 2 * ( 1 + JIDINT(
     $	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 (DREAL(ZNOW).LT.0.0D0) ZNOW=-ZNOW    ! for better approximation

	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

	CAFAKTZ = AFAKT / CDSQRT(ZNOW)
	XCHI=DREAL(ZNOW)-PIG4
	EIY=DCMPLX(0.0D0,DIMAG(ZNOW))

	IF (KAODDLOK.GT.0) THEN
	  QZ = AK(KAODDLOK)
	  DO K=KAODDLOK-2,1,-2
	    QZ = QZ*ZNOW2INV + AK(K)
	  END DO
	  QZ = QZ/ZNOW
	  CEIY=CDCOS(EIY)
	  SEIY=CDSIN(EIY)
	  BFAKT(1) = CAFAKTZ * ( PZ*CEIY - QZ*SEIY ) * DCOS(XCHI)
	  BFAKT(2) = CAFAKTZ * ( -QZ*CEIY - PZ*SEIY ) * DSIN(XCHI)
	ELSE
	  CSLASK = CAFAKTZ*PZ
	  BFAKT(1) = CSLASK*CDCOS(EIY) * DCOS(XCHI)
	  BFAKT(2) = -CSLASK*CDSIN(EIY) * DSIN(XCHI)
	END IF

	RETURN
	END






	SUBROUTINE H0ASCOMP (Z, HFAKT)
C	to compute H0(1)(Z) = HFAKT(1) + HFAKT(2)
C	using asymptotic expression
C	so that hfakt(1) is the term with dcos(x-pig4) and
C	        hfakt(2) is the term with dsin(x-pig4)
C	where x=Re(z)
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)           ! abbreviated
	INTEGER KAEVEN,KAODD
	REAL*8 AK

	COMPLEX*16 Z
	COMPLEX*16 HFAKT(2)

	INTEGER KAEVENLOK,KAODDLOK, K
	REAL*8 XCHI,Y, PIG4,AFAKT,TWOBYPI
	COMPLEX*16 ZNOW, PZ,QZ, ZNOW2INV, EIY, CAFAKTZ, CSLASK, 
     $	CSLASK1,CSLASK2, CEIY,SEIY, EIMULT
	LOGICAL HKVADR3

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

	ZNOW=Z

	K = 2 * ( 1 + JIDINT(
     $	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 (DREAL(ZNOW).LT.0.0D0) THEN
	  IF (DIMAG(ZNOW).LT.0.0D0) THEN
	    HKVADR3=.TRUE.
	    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

	CAFAKTZ = AFAKT / CDSQRT(ZNOW)
	XCHI=DREAL(ZNOW)-PIG4
	Y=DIMAG(ZNOW)
	EIY=DCMPLX(0.0D0,Y)

	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 (HKVADR3) THEN
C	    for the "first-quadrant znow" :   " H0(1) + 2*J0 " 
	    CEIY=CDCOS(EIY)
	    SEIY=CDSIN(EIY)
	    CSLASK1=3*PZ+EIMULT(QZ)
	    CSLASK2=-3*QZ+EIMULT(PZ)
	    HFAKT(1) = CAFAKTZ * ( CSLASK1*CEIY + CSLASK2*SEIY ) * DCOS(XCHI)
	    HFAKT(2) = CAFAKTZ * ( CSLASK2*CEIY - CSLASK1*SEIY ) * DSIN(XCHI)
	  ELSE
	    CSLASK = CAFAKTZ * (PZ+EIMULT(QZ)) * DEXP(-Y)
	    HFAKT(1) = CSLASK * DCOS(XCHI)
	    HFAKT(2) = EIMULT(CSLASK) * DSIN(XCHI)
	  END IF
	ELSE IF (HKVADR3) THEN
C	  for the "first-quadrant znow" :   " H0(1) + 2*J0 " 
	  CEIY=CDCOS(EIY)
	  SEIY=CDSIN(EIY)
	  CSLASK1=3*PZ
	  CSLASK2=EIMULT(PZ)
	  HFAKT(1) = CAFAKTZ * ( CSLASK1*CEIY + CSLASK2*SEIY ) * DCOS(XCHI)
	  HFAKT(2) = CAFAKTZ * ( CSLASK2*CEIY - CSLASK1*SEIY ) * DSIN(XCHI)
	ELSE
	  CSLASK = CAFAKTZ * PZ * DEXP(-Y)
	  HFAKT(1) = CSLASK * DCOS(XCHI)
	  HFAKT(2) = EIMULT(CSLASK) * DSIN(XCHI)
	END IF

	RETURN
	END







	SUBROUTINE J0ASCOMPPOS (Z, BFAKTPOS)
C	using asymptotic expression
C	so that bfaktpos(1) is the "cdexp(i*x) part of bfakt(1) from j0ascomp"
C	        bfaktpos(2) is the "cdexp(i*x) part of bfakt(2) from j0ascomp"
C	where x=Re(z)
C	note: bessrt must have been called to update /besst/
C	difference from j0ascomp:  dcos(xchi) --> cdexp(i*xchi)/2
C	                           dsin(xchi) --> cdexp(i*xchi)/(2*i)

	IMPLICIT NONE

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

	COMMON /BESST/ KAEVEN,KAODD,AK(kadim)           ! abbreviated
	INTEGER KAEVEN,KAODD
	REAL*8 AK

	COMPLEX*16 Z
	COMPLEX*16 BFAKTPOS(2)

	INTEGER KAEVENLOK,KAODDLOK, K
	REAL*8 PIG4,AFAKT,TWOBYPI
	COMPLEX*16 ZNOW, PZ,QZ, ZNOW2INV, EIY, CAFAKTZ, CSLASK, 
     $	CEIY,SEIY, EIXCHI,CXXCHIH, EIMULT

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

	ZNOW=Z

	K = 2 * ( 1 + JIDINT(
     $	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 (DREAL(ZNOW).LT.0.0D0) ZNOW=-ZNOW    ! for better approximation

	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

	CAFAKTZ = AFAKT / CDSQRT(ZNOW)
	EIXCHI = DCMPLX(0.0D0,DREAL(ZNOW)-PIG4)
	CXXCHIH=0.5D0*CDEXP(EIXCHI)
	EIY=DCMPLX(0.0D0,DIMAG(ZNOW))

	IF (KAODDLOK.GT.0) THEN
	  QZ = AK(KAODDLOK)
	  DO K=KAODDLOK-2,1,-2
	    QZ = QZ*ZNOW2INV + AK(K)
	  END DO
	  QZ = QZ/ZNOW
	  CEIY=CDCOS(EIY)
	  SEIY=CDSIN(EIY)
	  BFAKTPOS(1) = CAFAKTZ * ( PZ*CEIY - QZ*SEIY ) * CXXCHIH
	  BFAKTPOS(2) = CAFAKTZ * ( -QZ*CEIY - PZ*SEIY ) * (-EIMULT(CXXCHIH))
	ELSE
	  CSLASK = CAFAKTZ*PZ
	  BFAKTPOS(1) = CSLASK*CDCOS(EIY) * CXXCHIH
	  BFAKTPOS(2) = -CSLASK*CDSIN(EIY) * (-EIMULT(CXXCHIH))
	END IF

	RETURN
	END






	SUBROUTINE H0ASCOMPPOS (Z, HFAKTPOS)
C	using asymptotic expression
C	so that hfaktpos(1) is the "cdexp(i*x) part of hfakt(1) from h0ascomp"
C	        hfaktpos(2) is the "cdexp(i*x) part of hfakt(2) from h0ascomp"
C	where x=Re(z)
C	note: bessrt must have been called to update /besst/
C	difference from h0ascomp:  dcos(xchi) --> cdexp(i*xchi)/2
C	                           dsin(xchi) --> cdexp(i*xchi)/(2*i)

	IMPLICIT NONE

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

	COMMON /BESST/ KAEVEN,KAODD,AK(kadim)           ! abbreviated
	INTEGER KAEVEN,KAODD
	REAL*8 AK

	COMPLEX*16 Z
	COMPLEX*16 HFAKTPOS(2)

	INTEGER KAEVENLOK,KAODDLOK, K
	REAL*8 Y, PIG4,AFAKT,TWOBYPI
	COMPLEX*16 ZNOW, PZ,QZ, ZNOW2INV, EIY, CAFAKTZ, CSLASK, 
     $	CSLASK1,CSLASK2, CEIY,SEIY, EIXCHI,CXXCHIH, EIMULT
	LOGICAL HKVADR3

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

	ZNOW=Z

	K = 2 * ( 1 + JIDINT(
     $	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 (DREAL(ZNOW).LT.0.0D0) THEN
	  IF (DIMAG(ZNOW).LT.0.0D0) THEN
	    HKVADR3=.TRUE.
	    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

	CAFAKTZ = AFAKT / CDSQRT(ZNOW)
	EIXCHI = DCMPLX(0.0D0,DREAL(ZNOW)-PIG4)
	CXXCHIH=0.5D0*CDEXP(EIXCHI)
	Y=DIMAG(ZNOW)
	EIY=DCMPLX(0.0D0,Y)

	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 (HKVADR3) THEN
C	    for the "first-quadrant znow" :   " H0(1) + 2*J0 " 
	    CEIY=CDCOS(EIY)
	    SEIY=CDSIN(EIY)
	    CSLASK1=3*PZ+EIMULT(QZ)
	    CSLASK2=-3*QZ+EIMULT(PZ)
	    HFAKTPOS(1) = CAFAKTZ * ( CSLASK1*CEIY + CSLASK2*SEIY ) * 
     $	    CXXCHIH
	    HFAKTPOS(2) = CAFAKTZ * ( CSLASK2*CEIY - CSLASK1*SEIY ) * 
     $	    (-EIMULT(CXXCHIH))
	  ELSE
	    CSLASK = CAFAKTZ * (PZ+EIMULT(QZ)) * DEXP(-Y)
	    HFAKTPOS(1) = CSLASK * CXXCHIH
	    HFAKTPOS(2) = HFAKTPOS(1)   ! eimult(cslask) * (-eimult(cxxchih))
	  END IF
	ELSE IF (HKVADR3) THEN
C	  for the "first-quadrant znow" :   " H0(1) + 2*J0 " 
	  CEIY=CDCOS(EIY)
	  SEIY=CDSIN(EIY)
	  CSLASK1=3*PZ
	  CSLASK2=EIMULT(PZ)
	  HFAKTPOS(1) = CAFAKTZ * ( CSLASK1*CEIY + CSLASK2*SEIY ) * 
     $	  CXXCHIH
	  HFAKTPOS(2) = CAFAKTZ * ( CSLASK2*CEIY - CSLASK1*SEIY ) * 
     $	  (-EIMULT(CXXCHIH))
	ELSE
	  CSLASK = CAFAKTZ * PZ * DEXP(-Y)
	  HFAKTPOS(1) = CSLASK * CXXCHIH
	  HFAKTPOS(2) = HFAKTPOS(1)   ! eimult(cslask) * (-eimult(cxxchih))
	END IF

	RETURN
	END