C	subroutines to be used in dhcompfw (parts of "liqsi")
C	liqpart and airypart

C	Sven Ivansson  1992





	SUBROUTINE LIQPART ( OMEG,OMEG2,U,U2, 
     $	UA2, DZ, 
     $	URA, EXPTYP, 
     $	PLNREF, CP, STP,SDP )
C	an essential part of the subroutine dhcompfw

	IMPLICIT NONE

	INCLUDE 'RPPARAMETER.INC'

	REAL*8 DZ
	COMPLEX*16 OMEG,OMEG2,U,U2, UA2

	REAL*8 PLNREF
	COMPLEX*16 URA, CP, STP,SDP
	LOGICAL EXPTYP

	REAL*8 URANORM, SLASK, SLASK1,SLASK2
	COMPLEX*16 DZOMEG, P0, SP, CSLASK, CSLASK1,CSLASK2

C	choice of branch-cut for ura unimportant, the 	"minor" 
C	propagator elements are entire analytic functions of u; but see below
	URA = CDSQRT(UA2-U2)
	URANORM = DMAX1 ( DABS(DREAL(URA)), DABS(DIMAG(URA)) )

	DZOMEG = DZ*OMEG                           ! "dzomeg" could be stored
	P0 = DZOMEG*URA

C	for convenience we choose cut such that Im(p0) is nonnegative
	IF (DIMAG(P0).LT.0.0D0) THEN
	  URA = -URA
	  P0 = -P0
	END IF

	IF (URANORM.GT.exptyplim) THEN

	  EXPTYP = .TRUE.
C	  plnref := dimag(p0)
C	  < e1p := dexp(-plnref) >
C	  cp := cdexp(+i*p0) / e1p  ,  always of magn one
C	  stp is not set
C	  sdp := e1p*e1p , imaginary part not needed

	  PLNREF = DIMAG(P0)
	  SLASK = DREAL(P0)
	  CP = DCMPLX ( DCOS(SLASK), DSIN(SLASK) )
	  SDP = DEXP(-2*PLNREF)    ! imaginary part not needed

	ELSE

	  EXPTYP = .FALSE.
C	  possibly   plnref := dimag(p0),   always nonnegative
C	  cp := dexp(-plnref) * cdcos(p0)
C	  stp := dexp(-plnref) * cdsin(p0) * ura
C	  sdp := dexp(-plnref) * cdsin(p0) / ura

	  SLASK = DABS(DIMAG(P0))
	  IF (SLASK.GT.exptyplim) THEN
	    PLNREF = SLASK
	    SLASK = DREAL(P0)
	    CSLASK = DCMPLX ( DCOS(SLASK), DSIN(SLASK) )
	    CSLASK1 = DEXP(-2*PLNREF) * CSLASK
	    CSLASK2 = DCONJG(CSLASK)
	    CP = 0.5D0 * ( CSLASK1 + CSLASK2 )
	    CSLASK = CSLASK2 - CSLASK1
	    SP = 0.5D0 * DCMPLX(-DIMAG(CSLASK),DREAL(CSLASK))
	    STP = SP*URA
	    SDP = SP/URA
	  ELSE
	    PLNREF = 0.0D0
	    SLASK = DREAL(P0)
	    CSLASK = DCMPLX ( DCOS(SLASK), DSIN(SLASK) )
	    SLASK = DIMAG(P0)
	    CALL RHYPSMALL ( SLASK, SLASK1,SLASK2 )
	    CSLASK1 = SLASK1 * CSLASK
	    CSLASK2 = SLASK2 * CSLASK
	    CP = DCMPLX ( DREAL(CSLASK1), -DIMAG(CSLASK2) )
	    SP = DCMPLX ( DIMAG(CSLASK1), DREAL(CSLASK2) )
	    IF (URA.NE.(0.0D0,0.0D0)) THEN
	      STP = SP*URA
	      SDP = SP/URA
	    ELSE
	      CP = (1.0D0,0.0D0)
	      STP = (0.0D0,0.0D0)
	      SDP = DZOMEG
	    END IF
	  END IF

	END IF

	RETURN
	END









	SUBROUTINE AIRYPART ( OMEG,OMEG2,U,U2,RK,RK2, 
     $	UA2,UB2, DZ, IOTYP, 
     $	CFAKT,BCFAKT,BCFAKTINV, YA,YB, 
     $	PLNREF, CP,CPX, STP,SDP )
C	an essential part of the subroutine dhcompfw
C	the "middle part" of the "downwards" propagator matrix will be
C	       dexp(plnref) * pi * ( cp    stp  )
C	            "         "    ( -sdp cpx   )
C	the "middle part" of the "upwards" propagator matrix will be
C	       dexp(plnref) * pi * ( cpx   -stp )
C	            "         "    ( sdp  cp    )
C	in both cases, the surroundings are   diag(1,rho) * ... * diag(1,1/rho)

	IMPLICIT NONE

	INTEGER IOTYP   ! 0 or +-1
	REAL*8 DZ
	COMPLEX*16 OMEG,OMEG2,U,U2,RK,RK2, UA2,UB2

	COMPLEX*16 CFAKT,BCFAKT,BCFAKTINV     ! output/input if iotyp=0/+-1
	COMPLEX*16 YA,YB      ! output; iotyp=+1/-1 --> yb/ya input as well

	REAL*8 PLNREF
	COMPLEX*16 CP,CPX, STP,SDP

	INTEGER I
	REAL*8 ALN,BLN, SLASK, 
     $	EA(4),EB(4), EE(8)
	COMPLEX*16 BFAKT, CSLASK, 
     $	ABA(4),ABB(4), CC(8)

	IF (IOTYP.EQ.0) THEN
	  BFAKT = -OMEG2 * UB2 / DZ
	  IF (DIMAG(BFAKT).EQ.0.0D0) THEN
	    SLASK = DABS(DREAL(BFAKT))                    ! "dabs" (branch)
	    CFAKT = DCMPLX ( SLASK**(-0.6666666666666667D0) , 0.0D0 )
	  ELSE
	    IF (DREAL(BFAKT).GE.0.0D0) THEN
	      CFAKT = BFAKT**(-0.6666666666666667D0)
	    ELSE
	      CFAKT = (-BFAKT)**(-0.6666666666666667D0)   ! "-" (branch)
	    END IF
	  END IF
	  BCFAKT = -BFAKT*CFAKT / OMEG
	  BCFAKTINV = 1/BCFAKT
	  YA = CFAKT * ( RK2 - OMEG2*(UA2-0.5D0*UB2) )
	  YB = CFAKT * ( RK2 - OMEG2*(UA2+0.5D0*UB2) )
	  CALL ROT3 (YA,YB,CFAKT,BCFAKT,BCFAKTINV)
	  CALL AIRY4E (YA, ALN,BLN, ABA(1),ABA(2),ABA(3),ABA(4))
	  EA(1) = -ALN
	  EA(2) = -ALN
	  EA(3) = BLN
	  EA(4) = BLN
	  CALL AIRY4E (YB, ALN,BLN, ABB(1),ABB(2),ABB(3),ABB(4))
	  EB(1) = -ALN
	  EB(2) = -ALN
	  EB(3) = BLN
	  EB(4) = BLN
	ELSE
	  IF (IOTYP.GT.0) THEN
	    YA = YB
	    YB = CFAKT * ( RK2 - OMEG2*(UA2+0.5D0*UB2) )
	    CSLASK = CFAKT
	    CALL ROT3 (YA,YB,CFAKT,BCFAKT,BCFAKTINV)
	    IF (CFAKT.EQ.CSLASK) THEN
	      DO I = 1, 4
	        EA(I) = EB(I)
	        ABA(I) = ABB(I)
	      END DO
	    ELSE
	      CALL AIRY4E (YA, ALN,BLN, ABA(1),ABA(2),ABA(3),ABA(4))
	      EA(1) = -ALN
	      EA(2) = -ALN
	      EA(3) = BLN
	      EA(4) = BLN
	    END IF
	    CALL AIRY4E (YB, ALN,BLN, ABB(1),ABB(2),ABB(3),ABB(4))
	    EB(1) = -ALN
	    EB(2) = -ALN
	    EB(3) = BLN
	    EB(4) = BLN
	  ELSE
	    YB = YA
	    YA = CFAKT * ( RK2 - OMEG2*(UA2-0.5D0*UB2) )
	    CSLASK = CFAKT
	    CALL ROT3 (YA,YB,CFAKT,BCFAKT,BCFAKTINV)
	    IF (CFAKT.EQ.CSLASK) THEN
	      DO I = 1, 4
	        EB(I) = EA(I)
	        ABB(I) = ABA(I)
	      END DO
	    ELSE
	      CALL AIRY4E (YB, ALN,BLN, ABB(1),ABB(2),ABB(3),ABB(4))
	      EB(1) = -ALN
	      EB(2) = -ALN
	      EB(3) = BLN
	      EB(4) = BLN
	    END IF
	    CALL AIRY4E (YA, ALN,BLN, ABA(1),ABA(2),ABA(3),ABA(4))
	    EA(1) = -ALN
	    EA(2) = -ALN
	    EA(3) = BLN
	    EA(4) = BLN
	  END IF
	END IF

	EE(1) = EA(1) + EB(4)
	CC(1) = ABA(1) * ABB(4)
	EE(2) = EA(3) + EB(2)
	CC(2) = ABA(3) * ABB(2)
	EE(3) = EA(4) + EB(2)
	CC(3) = ABA(4) * ABB(2)
	EE(4) = EA(2) + EB(4)
	CC(4) = ABA(2) * ABB(4)
	EE(5) = EA(1) + EB(3)
	CC(5) = ABA(1) * ABB(3)
	EE(6) = EA(3) + EB(1)
	CC(6) = ABA(3) * ABB(1)
	EE(7) = EA(4) + EB(1)
	CC(7) = ABA(4) * ABB(1)
	EE(8) = EA(2) + EB(3)
	CC(8) = ABA(2) * ABB(3)
	PLNREF = DMAX1 ( EE(1),EE(2),EE(3),EE(4),EE(5),EE(6),EE(7),EE(8) )
	IF (PLNREF.NE.0.0D0) THEN
	  DO I = 1, 8
	    EE(I) = EE(I) - PLNREF
	  END DO
	END IF
	CALL ADDEC ( EE(2),-CC(2), EE(1),CC(1), SLASK,CP )
	CP = DEXP(SLASK)*CP
	CALL ADDEC ( EE(4),-CC(4), EE(3),CC(3), SLASK,STP )
	STP = BCFAKT * ( DEXP(SLASK)*STP )
	CALL ADDEC ( EE(6),-CC(6), EE(5),CC(5), SLASK,SDP )
	SDP = -BCFAKTINV * ( DEXP(SLASK)*SDP )
	CALL ADDEC ( EE(8),-CC(8), EE(7),CC(7), SLASK,CPX )
	CPX = DEXP(SLASK)*CPX

	RETURN
	END