C       subroutines for rpress.f
C	rcomp r1comp r1comppos

C	Sven Ivansson  1990






	SUBROUTINE RCOMP (UP, RFAKT)
C	to compute  J0(k*r) or 0.5d0*H0(1)(k*r) or CDEXP(i*k*r)  
C             resp. J1(k*r) or 0.5d0*H1(1)(k*r) or CDEXP(i*k*r)
C	for the k-value rk=rk0*up and 
C	for all li r-values in the or array

	IMPLICIT NONE

	INCLUDE 'RPPARAMETER.INC'
C       ..., li <= lidim

	INCLUDE 'INRANGE.INC'
	INCLUDE 'INPARAM.INC'	! only " kfwbess, omeg,omeg2, u0,rk0 " needed
	INCLUDE 'INSPEC.INC'

	COMPLEX*16 UP
	COMPLEX*16 RFAKT(*)

	INTEGER IR
	REAL*8 DREALRK, X, SLASK, PIG4, DBES0RAT, DBES1RAT
	COMPLEX*16 RK, CDLOGZBY2, Z, CSLASK, CSLASK1, BESSRA, 
     $  CBES0BOD,CHAN0BOD, DHAN0RAT, CHAN0SPEC, 
     $  CBES1BOD,CHAN1BOD, DHAN1RAT, CHAN1SPEC

	DATA PIG4 /0.7853981633974483D0/        ! pi/4
	REAL*8 PI
	DATA PI /3.141592653589793D0/

	RK = RK0*UP
	IF (IBESSTYP.LE.0) GOTO 100

C       here J0 resp. J1
	IF (IRKSPEC.NE.0) STOP 76542 ! should not be used in this case
	IF (IRKREAL.NE.0) THEN
	  DREALRK = DREAL(RK)
	  IF (KFWBESS.EQ.0) THEN
	    DO IR = 1, LI
	      X = DREALRK*OR(IR)
	      RFAKT(IR) = DBES0RAT (X)
	    END DO
	  ELSE
	    DO IR = 1, LI
	      X = DREALRK*OR(IR)
	      RFAKT(IR) = DBES1RAT (X)
	    END DO
	  END IF
	ELSE IF (IBESSTYP.EQ.1) THEN
	  IF (KFWBESS.NE.0) STOP 76543	! not implemented
	  DO IR = 1, LI
	    Z = RK*OR(IR)
	    RFAKT(IR) = BESSRA (0,CDLOGZBY2,Z)
	  END DO
	ELSE
	  IF (KFWBESS.EQ.0) THEN
	    DO IR = 1, LI
	      Z = RK*OR(IR)
	      RFAKT(IR) = CBES0BOD (Z)
	    END DO
	  ELSE
	    DO IR = 1, LI
	      Z = RK*OR(IR)
	      RFAKT(IR) = CBES1BOD (Z)
	    END DO
	  END IF
	END IF
	RETURN

C       here 0.5d0*H0(1) resp. 0.5d0*H1(1)
	ENTRY RCOMPPOS (UP, RFAKT)
	RK = RK0*UP
 100	IF (IBESSTYP.EQ.0) GOTO 200

	IF (IRKSPEC.NE.0) THEN
	  IF (IBESSTYP.NE.-2) STOP 76542 ! should not be used in this case
	  DLNSPEC = ORMINSPEC * DIMAG(RK)
	  IF (KFWBESS.EQ.0) THEN
	    DO IR = 1, LI
	      Z = RK*OR(IR)
	      RFAKT(IR) = 0.5D0 * CHAN0SPEC (Z)
	    END DO
	  ELSE
	    DO IR = 1, LI
	      Z = RK*OR(IR)
	      RFAKT(IR) = 0.5D0 * CHAN1SPEC (Z)
	    END DO
	  END IF
	  RETURN
	END IF

	IF (IRKSPEC.NE.0) STOP 76542 ! should not be used in this case
	IF (RK.NE.(0.0D0,0.0D0)) THEN
	  IF (IBESSTYP.EQ.-3) THEN
	    SLASK = RMIN * DREAL(RK)
	    IF (KFWBESS.EQ.0) THEN
	      CSLASK = DCMPLX(1.0D0,-1.0D0) * DEXP(-RMIN*DIMAG(RK)) *
     $        DCMPLX ( DCOS(SLASK), DSIN(SLASK) ) / CDSQRT(RK)
	    ELSE
	      CSLASK = DCMPLX(-1.0D0,-1.0D0) * DEXP(-RMIN*DIMAG(RK)) *
     $        DCMPLX ( DCOS(SLASK), DSIN(SLASK) ) / CDSQRT(RK)
	    END IF
	    RFAKT(1) = RS(1) * CSLASK
	    SLASK = DR * DREAL(RK)
	    CSLASK1 = DEXP(-DR*DIMAG(RK)) * 
     $	    DCMPLX ( DCOS(SLASK), DSIN(SLASK) )
	    DO IR = 2, LI
	      CSLASK = CSLASK * CSLASK1
	      RFAKT(IR) = RS(IR) * CSLASK
	    END DO
	  ELSE IF (IBESSTYP.EQ.-4) THEN
	    IF (KFWBESS.NE.0) STOP 76543 ! not implemented
	    DO IR = 1, LI
	      SLASK = OR(IR) * DREAL(RK)
	      CSLASK = 0.5D0*DCMPLX(1.0D0,-1.0D0) * DEXP(-OR(IR)*DIMAG(RK)) *
     $        DCMPLX ( DCOS(SLASK), DSIN(SLASK) )
	      RFAKT(IR) = CSLASK / CDSQRT(1+PI*RK*OR(IR))
	    END DO		! weinberg approx
	  ELSE IF (IRKREAL.NE.0) THEN
	    DREALRK = DREAL(RK)
	    IF (KFWBESS.EQ.0) THEN
	      DO IR = 1, LI
		X = DREALRK*OR(IR)
		RFAKT(IR) = 0.5D0 * DHAN0RAT (X)
	      END DO
	    ELSE
	      DO IR = 1, LI
		X = DREALRK*OR(IR)
		RFAKT(IR) = 0.5D0 * DHAN1RAT (X)
	      END DO
	    END IF
	  ELSE IF (IABS(IBESSTYP).EQ.1) THEN
	    IF (KFWBESS.NE.0) STOP 76543 ! not implemented
	    CSLASK = CDLOG(RK)
	    DO IR = 1, LI
	      CDLOGZBY2 = CSLASK + RL(IR)
	      Z = RK*OR(IR)
	      RFAKT(IR) = 0.5D0 * BESSRA (-1,CDLOGZBY2,Z)
	    END DO
	  ELSE
	    IF (KFWBESS.EQ.0) THEN
	      DO IR = 1, LI
		Z = RK*OR(IR)
		RFAKT(IR) = 0.5D0 * CHAN0BOD (Z)
	      END DO
	    ELSE
	      DO IR = 1, LI
		Z = RK*OR(IR)
		RFAKT(IR) = 0.5D0 * CHAN1BOD (Z)
	      END DO
	    END IF
	  END IF
	ELSE
C         "overflow", choose symmetry in "J0 = 0.5d0 * ( H0(1) + H0(2) )"
C                                  resp. "J1 = 0.5d0 * ( H1(1) + H1(2) )"
	  IF (KFWBESS.EQ.0) THEN
	    DO IR = 1, LI
	      RFAKT(IR) = (0.5D0,0.0D0)
	    END DO
	  ELSE
	    DO IR = 1, LI
	      RFAKT(IR) = (0.0D0,0.0D0)
	    END DO
	  END IF
	END IF
	RETURN

 200    IF (IRKSPEC.NE.0) STOP 76542 ! should not be used in this case
	IF (RK.NE.(0.0D0,0.0D0)) THEN
	  DO IR = 1, LI
	    Z = RK*OR(IR)
	    RFAKT(IR) = CDEXP ( DCMPLX(-DIMAG(Z),DREAL(Z)) )  ! cdexp(i*k*r)
	  END DO
	ELSE
	  DO IR = 1, LI
	    RFAKT(IR) = (1.0D0,0.0D0)
	  END DO
	END IF
	RETURN

	END








	SUBROUTINE R1COMP (UP, R1FAKT)
C	to compute "parts" of  J0(k*r) or 0.5d0*H0(1)(k*r) or CDEXP(i*k*r)  
C	for the k-value rk=rk0*up and for a certain r-value in the or array
C	let z=k*r and  x=Re(z)
C	  - for J0
C	      J0(z) = R1FAKT(1) + R1FAKT(2)
C	      so that r1fakt(1) is the term with dcos(x-pig4) and
C	              r1fakt(2) is the term with dsin(x-pig4)
C	      note: when Im(z)=0.0d0 we put r1fakt(2)=(0.0d0,0.0d0)
C	            continuity of the "parts" when Im(z)-->0.0d0 is only 
C	            achieved asymptotically so in a "call sequence" z with 
C	            Im(z)=0 and "nonzero Im(z)" should not be mixed
C	  - for 0.5d0*H0(1)
C	      0.5d0*H0(1)(z) = R1FAKT(1) + R1FAKT(2)
C	      so that r1fakt(1) is the term with dcos(x-pig4) and
C	              r1fakt(2) is the term with dsin(x-pig4)
C	  - for CDEXP
C	      CDEXP(i*z) = R1FAKT(1) + R1FAKT(2)
C	      so that r1fakt(1) is the term with dcos(x) and
C	              r1fakt(2) is the term with dsin(x)

	IMPLICIT NONE

	INCLUDE 'RPPARAMETER.INC'
C        ... , li <= lidim

	INCLUDE 'INRANGE.INC'
	INCLUDE 'INTAIL.INC'     ! only  " nhyd2, rk0or "  needed
	INCLUDE 'INPARAM.INC'	 ! only  " kfwbess, rk0 "  needed
	INCLUDE 'INSPEC.INC'

	COMPLEX*16 UP
	COMPLEX*16 R1FAKT(2)

	REAL*8 X, XCHI, XCOS,XSIN, SLASK, PIG4, DBES0RAT
	COMPLEX*16 Z, CH1,CH2, CSLASK, CSLASK1,CSLASK2, 
     $	BESSRA, CBES0BOD,CHAN0BOD, DHAN0RAT, CHAN0SPEC, 
     $	B(2), H(2)

	DATA PIG4 /0.7853981633974483D0/        ! pi/4

	IF (KFWBESS.NE.0) STOP 76543 ! not implemented

	Z = RK0OR * UP

	IF (IBESSTYP.LE.0) GOTO 100

	IF (IRKSPEC.NE.0) STOP 76542 ! should not be used in this case
	IF (IRKREAL.NE.0) THEN
	  X = DREAL(Z)
	  R1FAKT(1) = DBES0RAT (X)
	  R1FAKT(2) = (0.0D0,0.0D0)
	ELSE IF (DIMAG(Z).EQ.0.0D0) THEN
	  IF (IBESSTYP.EQ.1) THEN
	    R1FAKT(1) = BESSRA (0,CSLASK,Z)
	  ELSE
	    R1FAKT(1) = CBES0BOD (Z)
	  END IF
	  R1FAKT(2) = (0.0D0,0.0D0)
	ELSE
	  IF (IBESSTYP.EQ.1) THEN
	    B(1) = BESSRA (0,CSLASK,Z)
	    H(1) = BESSRA (1,CSLASK,Z)
	  ELSE
	    CALL CDHAN (1,Z,B,H)
	  END IF
	  XCHI = DREAL(Z) - PIG4
	  XCOS = DCOS(XCHI)
	  XSIN = DSIN(XCHI)
	  CSLASK1 = DCMPLX(XCOS,-XSIN)   ! cdexp(-i*xchi)
	  CSLASK2 = DCONJG(CSLASK1)      ! cdexp(+i*xchi)
	  CSLASK = 0.5D0 * H(1)
	  CH1 = CSLASK * CSLASK1          ! 0.5*h0(1)(z) * cdexp(-i*(x-pig4))
	  CH2 = (B(1)-CSLASK) * CSLASK2   ! 0.5*h0(2)(z) * cdexp(+i*(x-pig4))
	  R1FAKT(1) = (CH1+CH2) * XCOS
	  CSLASK = CH1 - CH2
	  R1FAKT(2) = DCMPLX(-DIMAG(CSLASK),DREAL(CSLASK)) * XSIN ! "i*"
	END IF
	RETURN

 100	IF (IBESSTYP.EQ.0) GOTO 200

	IF (IRKSPEC.NE.0) THEN
	  IF (IBESSTYP.NE.-2) STOP 76542 ! should not be used in this case
	  DLNSPEC = ORMINSPEC * DIMAG(RK0*UP)
	  CSLASK = 0.5D0 * CHAN0SPEC (Z)
	ELSE
	  IF (IRKSPEC.NE.0) STOP 76542 ! should not be used in this case
	  IF (Z.NE.(0.0D0,0.0D0)) THEN
	    IF (IRKREAL.NE.0) THEN
	      X = DREAL(Z)
	      CSLASK = 0.5D0 * DHAN0RAT(X)
	    ELSE IF (IBESSTYP.EQ.-1) THEN
	      CSLASK = 0.5D0 * BESSRA (1,CSLASK,Z)
	    ELSE
	      CSLASK = 0.5D0 * CHAN0BOD (Z)
	    END IF
	  ELSE
C	    "overflow", choose symmetry in  " J0 = 0.5d0 * ( H0(1) + H0(2) ) "
	    CSLASK = (0.5D0,0.0D0)
	  END IF
	END IF
	XCHI = DREAL(Z) - PIG4
	XCOS = DCOS(XCHI)
	XSIN = DSIN(XCHI)
	CSLASK = CSLASK * DCMPLX(XCOS,-XSIN)   ! cslask * cdexp(-i*xchi)
	R1FAKT(1) = CSLASK * XCOS
	R1FAKT(2) = DCMPLX(-DIMAG(CSLASK),DREAL(CSLASK)) * XSIN    ! i*
	RETURN

 200    IF (IRKSPEC.NE.0) STOP 76542 ! should not be used in this case
	IF (Z.NE.(0.0D0,0.0D0)) THEN
	  SLASK = DEXP(-DIMAG(Z))
	  X = DREAL(Z)
	  R1FAKT(1) = DCMPLX ( SLASK*DCOS(X) , 0.0D0 )
	  R1FAKT(2) = DCMPLX ( 0.0D0 , SLASK*DSIN(X) )
	ELSE
	  R1FAKT(1) = (1.0D0,0.0D0)
	  R1FAKT(2) = (0.0D0,0.0D0)
	END IF
	RETURN

	END






	SUBROUTINE R1COMPPOS (UP, R1FAKTPOS)
C	to compute "pos parts" of "parts" of  J0(k*r) or 0.5d0*H0(1)(k*r) or 
C	CDEXP(i*k*r)  for the k-value rk=rk0*up and for a certain r-value in 
C	the or array
C	      r1faktpos(1) is the "cdexp(i*x) part of r1fakt(1) from r1comp"
C	      r1faktpos(2) is the "cdexp(i*x) part of r1fakt(2) from r1comp"
C	      where z=k*r and x=Re(z)
C	difference from r1comp for J0 and 0.5d0*H0(1):  
C	                         dcos(xchi) --> cdexp(i*xchi)/2
C	                         dsin(xchi) --> cdexp(i*xchi)/(2*i)
C	                  ( and J0 --> 0.5d0*H0(1) for "J0 with Im(z)=0.0d0" )
C	difference from r1comp for CDEXP:
C	                         dcos(x) --> cdexp(i*x)/2
C	                         dsin(x) --> cdexp(i*x)/(2*i)

	IMPLICIT NONE

	INCLUDE 'RPPARAMETER.INC'
C        ... , li <= lidim

	INCLUDE 'INRANGE.INC'
	INCLUDE 'INTAIL.INC'     ! only  " nhyd2, rk0or "  needed
	INCLUDE 'INPARAM.INC'	 ! only  " kfwbess, rk0 "  needed
	INCLUDE 'INSPEC.INC'

	COMPLEX*16 UP
	COMPLEX*16 R1FAKTPOS(2)

	REAL*8 X, XCHI, SLASK, PIG4
	COMPLEX*16 Z, CH1,CH2, CSLASK, CSLASK2, 
     $	BESSRA, CHAN0BOD, DHAN0RAT, CHAN0SPEC, 
     $	B(2), H(2)

	DATA PIG4 /0.7853981633974483D0/        ! pi/4

	IF (KFWBESS.NE.0) STOP 76543 ! not implemented

	Z = RK0OR * UP

	IF (IBESSTYP.LE.0) GOTO 100

	IF (IRKSPEC.NE.0) STOP 76542 ! should not be used in this case
	IF (IRKREAL.NE.0) THEN
	  X = DREAL(Z)
	  R1FAKTPOS(1) = 0.5D0 * DHAN0RAT(X)
	  R1FAKTPOS(2) = (0.0D0,0.0D0)
	ELSE IF (DIMAG(Z).EQ.0.0D0) THEN
	  IF (IBESSTYP.EQ.1) THEN
	    R1FAKTPOS(1) = 0.5D0 * BESSRA (1,CSLASK,Z)
	  ELSE
	    R1FAKTPOS(1) = 0.5D0 * CHAN0BOD (Z)
	  END IF
	  R1FAKTPOS(2) = (0.0D0,0.0D0)
	ELSE
	  IF (IBESSTYP.EQ.1) THEN
	    B(1) = BESSRA (0,CSLASK,Z)
	    H(1) = BESSRA (1,CSLASK,Z)
	  ELSE
	    CALL CDHAN (1,Z,B,H)
	  END IF
	  XCHI = DREAL(Z) - PIG4
	  SLASK = 2*XCHI
	  CSLASK2 = DCMPLX ( DCOS(SLASK), DSIN(SLASK) )   ! cdexp(2*i*xchi)
	  CSLASK = 0.5D0 * H(1)
	  CH1 = CSLASK
	  CH2 = (B(1)-CSLASK) * CSLASK2
	  R1FAKTPOS(1) = (CH1+CH2) * 0.5D0
	  R1FAKTPOS(2) = (CH1-CH2) * 0.5D0
	END IF
	RETURN

 100	IF (IBESSTYP.EQ.0) GOTO 200

	IF (IRKSPEC.NE.0) THEN
	  IF (IBESSTYP.NE.-2) STOP 76542 ! should not be used in this case
	  DLNSPEC = ORMINSPEC * DIMAG(RK0*UP)
	  CSLASK = 0.5D0 * CHAN0SPEC (Z)
	ELSE
	  IF (IRKSPEC.NE.0) STOP 76542 ! should not be used in this case
	  IF (Z.NE.(0.0D0,0.0D0)) THEN
	    IF (IRKREAL.NE.0) THEN
	      X = DREAL(Z)
	      CSLASK = 0.25D0 * DHAN0RAT(X)
	    ELSE IF (IBESSTYP.EQ.-1) THEN
	      CSLASK = 0.25D0 * BESSRA (1,CSLASK,Z)
	    ELSE
	      CSLASK = 0.25D0 * CHAN0BOD (Z)
	    END IF
	  ELSE
C	    "overflow", choose symmetry in  " J0 = 0.5d0 * ( H0(1) + H0(2) ) "
	    CSLASK = (0.25D0,0.0D0)
	  END IF
	END IF
	R1FAKTPOS(1) = CSLASK
	R1FAKTPOS(2) = R1FAKTPOS(1)
	RETURN

 200    IF (IRKSPEC.NE.0) STOP 76542 ! should not be used in this case
	IF (Z.NE.(0.0D0,0.0D0)) THEN
	  CSLASK = 0.5D0 * CDEXP ( DCMPLX(-DIMAG(Z),DREAL(Z)) ) !0.5*cdexp(i*z)
	ELSE
	  CSLASK = (0.5D0,0.0D0)
	END IF
	R1FAKTPOS(1) = CSLASK
	R1FAKTPOS(2) = R1FAKTPOS(1)
	RETURN

	END






	SUBROUTINE RCOMPONE (UP, RFAKTONE)
C	dummy to set  rfaktone(1) = (1.0d0,0.0d0)

	IMPLICIT NONE

	INCLUDE 'INPARAM.INC'	 ! only  " kfwbess "  needed
	INCLUDE 'INSPEC.INC'     ! only  " irkspec "  needed

	COMPLEX*16 UP
	COMPLEX*16 RFAKTONE(1)

	IF (IRKSPEC.NE.0) STOP 76542 ! should not be used in this case

	RFAKTONE(1) = (1.0D0,0.0D0)

	RETURN
	END





	SUBROUTINE RCOMPPOSONE (UP, RFAKTONE)
C	dummy to set  rfaktone(1) = (0.5d0,0.0d0)

	IMPLICIT NONE

	INCLUDE 'INPARAM.INC'    ! only  " kfwbess "  needed
	INCLUDE 'INSPEC.INC'     ! only  " irkspec "  needed

	COMPLEX*16 UP
	COMPLEX*16 RFAKTONE(1)

	IF (IRKSPEC.NE.0) STOP 76542 ! should not be used in this case

	RFAKTONE(1) = (0.5D0,0.0D0)

	RETURN
	END