
C	subroutine to be used in dhcompfw, related to "harksipart & solsipart"

C	Sven Ivansson  1992





	SUBROUTINE HSOLPART ( OMEG,OMEG2,U,U2, 
     $	UA2,UB2, DZ, 
     $	URA,URB, EXPTYP, 
     $	SLNREF,E1X, CP,CQ, STP,STQ,SDP,SDQ, 
     $	TLNREF,E2T )
C	an essential part of the subroutine dhcompfw

	IMPLICIT NONE

	INCLUDE '../RPPARAMETER.INC'

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

	REAL*8 SLNREF,E1X, TLNREF,E2T
	COMPLEX*16 URA,URB, CP,CQ, STP,STQ,SDP,SDQ
	LOGICAL EXPTYP

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

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

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

C	for convenience we choose cuts such that Im(p0),Im(q0) are nonnegative
	IF (DIMAG(P0).LT.0.0D0) THEN
	  URA = -URA
	  P0 = -P0
        ELSE IF (DZ.EQ.0.0D0) THEN
C         dz thought of as "+epsilon" in this case
	  IF (DIMAG(OMEG*URA).LT.0.0D0) THEN
	    URA = -URA
	    P0 = -P0  ! zero
          END IF
	END IF
	IF (DIMAG(Q0).LT.0.0D0) THEN
	  URB = -URB
	  Q0 = -Q0
        ELSE IF (DZ.EQ.0.0D0) THEN
C         dz thought of as "+epsilon" in this case
	  IF (DIMAG(OMEG*URB).LT.0.0D0) THEN
	    URB = -URB
	    Q0 = -Q0  ! zero
          END IF
	END IF

	IF ( URANORM.GT.exptyplim .AND. URBNORM.GT.exptyplim ) THEN

	  EXPTYP = .TRUE.
C	  < plnref := dimag(p0) >
C	  < qlnref := dimag(q0) >
C	  slnref := if (plnref.ge.qlnref) then plnref else -qlnref endif  %OBS%
C	  < e1p := dexp(-plnref) >
C	  < e1q := dexp(-qlnref) >
C	  < e1s := dexp(-dabs(slnref)) >
C	  e1x := dexp ( -dabs(plnref-qlnref) )
C	  cp := cdexp(+i*p0) / e1p  ,  always of magn one
C	  cq := cdexp(+i*q0) / e1q  ,  always of magn one
C	  stp  :=       cdexp(+i*p0+i*q0) / (e1p*e1q)  ,  always of magn one
C	  stq  := e1q * cdexp(+i*p0-i*q0) / e1p        ,  always of magn one
C	  sdp := e1p*e1p , imaginary part not needed
C	  sdq := e1q*e1q , imaginary part not needed
C	  tlnref := plnref + qlnref
C	  e2t := dexp(-tlnref) = e1p*e1q

C	  now cp,cq
	  SLASK = DREAL(P0)
	  CP = DCMPLX ( DCOS(SLASK), DSIN(SLASK) )
	  SLASK = DREAL(Q0)
	  CQ = DCMPLX ( DCOS(SLASK), DSIN(SLASK) )
C	  now stp,stq
	  STP = CP * CQ
	  STQ = CP * DCONJG(CQ)

	  SLASK1 = DIMAG(P0)
	  SLASK2 = DIMAG(Q0)
	  TLNREF = SLASK1 + SLASK2
	  IF (SLASK1.GE.SLASK2) THEN
	    SLNREF = SLASK1
	    E1X = DEXP(SLASK2-SLASK1)
	    SLASK2 = DEXP(-SLASK2)
	    E1S = E1X*SLASK2
	    E2T = E1S*SLASK2
	    SDP = E1S*E1S         ! imaginary part not needed
	    SDQ = SLASK2*SLASK2   ! imaginary part not needed
	  ELSE
	    SLNREF = -SLASK2
	    E1X = DEXP(SLASK1-SLASK2)
	    SLASK1 = DEXP(-SLASK1)
	    E1S = SLASK1*E1X
	    E2T = SLASK1*E1S
	    SDP = SLASK1*SLASK1   ! imaginary part not needed
	    SDQ = E1S*E1S         ! imaginary part not needed
	  END IF

	ELSE

	  EXPTYP = .FALSE.
C	  < possibly   plnref := dimag(p0),   always nonnegative >
C	  < possibly   qlnref := dimag(q0),   always nonnegative >
C	  slnref := if (plnref.ge.qlnref) then plnref else -qlnref endif  %OBS%
C	  < e1p := dexp(-plnref) >
C	  < e1q := dexp(-qlnref) >
C	  < e1s := dexp(-dabs(slnref)) >
C	  e1x := dexp ( -dabs(plnref-qlnref) )
C	  cp := e1p * cdcos(p0)   ,   cq := e1q * cdcos(q0)
C	  stp := e1p * cdsin(p0) * ura
C	  stq := e1q * cdsin(q0) * urb
C	  sdp := e1p * cdsin(p0) / ura
C	  sdq := e1q * cdsin(q0) / urb
C	  tlnref := plnref + qlnref
C	  e2t := dexp(-tlnref)

	  SLASK1 = DIMAG(P0)
	  SLASK2 = DIMAG(Q0)

	  IF ( URANORM.GT.exptyplim .OR. SLASK1.GT.exptyplim ) THEN

	    IF ( URBNORM.GT.exptyplim .OR. SLASK2.GT.exptyplim ) THEN

	      TLNREF = SLASK1 + SLASK2
	      IF (SLASK1.GE.SLASK2) THEN
	        SLNREF = SLASK1
	        SLASK1 = DEXP(SLASK2-SLASK1)
	        SLASK2 = DEXP(-SLASK2)
	        E1S = SLASK1 * SLASK2
	        E2T = E1S * SLASK2
	        SLASK = DREAL(P0)
	        CSLASK = DCMPLX ( DCOS(SLASK), DSIN(SLASK) )
	        CSLASK1 = (E1S*E1S) * 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
	        SLASK = DREAL(Q0)
	        CSLASK = DCMPLX ( DCOS(SLASK), DSIN(SLASK) )
	        CSLASK1 = (SLASK2*SLASK2) * CSLASK
	        CSLASK2 = DCONJG(CSLASK)
	        CQ = 0.5D0 * ( CSLASK1 + CSLASK2 )
	        CSLASK = CSLASK2 - CSLASK1
	        SQ = 0.5D0 * DCMPLX(-DIMAG(CSLASK),DREAL(CSLASK))
	        STQ = SQ*URB
	        SDQ = SQ/URB
	        E1X = SLASK1
	      ELSE
	        SLNREF = -SLASK2
	        SLASK2 = DEXP(SLASK1-SLASK2)
	        SLASK1 = DEXP(-SLASK1)
	        E1S = SLASK1 * SLASK2
	        E2T = E1S * SLASK1
	        SLASK = DREAL(Q0)
	        CSLASK = DCMPLX ( DCOS(SLASK), DSIN(SLASK) )
	        CSLASK1 = (E1S*E1S) * CSLASK
	        CSLASK2 = DCONJG(CSLASK)
	        CQ = 0.5D0 * ( CSLASK1 + CSLASK2 )
	        CSLASK = CSLASK2 - CSLASK1
	        SQ = 0.5D0 * DCMPLX(-DIMAG(CSLASK),DREAL(CSLASK))
	        STQ = SQ*URB
	        SDQ = SQ/URB
	        SLASK = DREAL(P0)
	        CSLASK = DCMPLX ( DCOS(SLASK), DSIN(SLASK) )
	        CSLASK1 = (SLASK1*SLASK1) * 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
	        E1X = SLASK2
	      END IF

	    ELSE

	      SLNREF = DIMAG(P0)
	      TLNREF = SLNREF
	      E1S = DEXP(-SLNREF)
	      E2T = E1S
	      SLASK = DREAL(P0)
	      CSLASK = DCMPLX ( DCOS(SLASK), DSIN(SLASK) )
	      CSLASK1 = (E1S*E1S) * 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
	      IF (URB.NE.(0.0D0,0.0D0)) THEN
	        SLASK = DREAL(Q0)
	        CSLASK = DCMPLX ( DCOS(SLASK), DSIN(SLASK) )
	        SLASK = DIMAG(Q0)
	        CALL RHYPSMALL ( SLASK, SLASK1,SLASK2 )
	        CSLASK1 = SLASK1 * CSLASK
	        CSLASK2 = SLASK2 * CSLASK
	        CQ = DCMPLX ( DREAL(CSLASK1), -DIMAG(CSLASK2) )
	        SQ = DCMPLX ( DIMAG(CSLASK1), DREAL(CSLASK2) )
	        STQ = SQ*URB
	        SDQ = SQ/URB
	      ELSE
	        CQ = DCMPLX(1.0D0,0.0D0)
	        STQ = (0.0D0,0.0D0)
	        SDQ = DZOMEG
	      END IF
	      E1X = E1S

	    END IF

	  ELSE

	    IF ( URBNORM.GT.exptyplim .OR. SLASK2.GT.exptyplim ) THEN

	      SLNREF = -DIMAG(Q0)
	      TLNREF = -SLNREF
	      E1S = DEXP(SLNREF)
	      E2T = E1S
	      SLASK = DREAL(Q0)
	      CSLASK = DCMPLX ( DCOS(SLASK), DSIN(SLASK) )
	      CSLASK1 = (E1S*E1S) * CSLASK
	      CSLASK2 = DCONJG(CSLASK)
	      CQ = 0.5D0 * ( CSLASK1 + CSLASK2 )
	      CSLASK = CSLASK2 - CSLASK1
	      SQ = 0.5D0 * DCMPLX(-DIMAG(CSLASK),DREAL(CSLASK))
	      STQ = SQ*URB
	      SDQ = SQ/URB
	      IF (URA.NE.(0.0D0,0.0D0)) THEN
	        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) )
	        STP = SP*URA
	        SDP = SP/URA
	      ELSE
	        CP = DCMPLX(1.0D0,0.0D0)
	        STP = (0.0D0,0.0D0)
	        SDP = DZOMEG
	      END IF
	      E1X = E1S

	    ELSE

	      SLNREF = 0.0D0
	      TLNREF = 0.0D0
	      E1S = 1.0D0
	      E2T = 1.0D0
	      IF (URA.NE.(0.0D0,0.0D0)) THEN
	        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) )
	        STP = SP*URA
	        SDP = SP/URA
	      ELSE
	        CP = (1.0D0,0.0D0)
	        STP = (0.0D0,0.0D0)
	        SDP = DZOMEG
	      END IF
	      IF (URB.NE.(0.0D0,0.0D0)) THEN
	        SLASK = DREAL(Q0)
	        CSLASK = DCMPLX ( DCOS(SLASK), DSIN(SLASK) )
	        SLASK = DIMAG(Q0)
	        CALL RHYPSMALL ( SLASK, SLASK1,SLASK2 )
	        CSLASK1 = SLASK1 * CSLASK
	        CSLASK2 = SLASK2 * CSLASK
	        CQ = DCMPLX ( DREAL(CSLASK1), -DIMAG(CSLASK2) )
	        SQ = DCMPLX ( DIMAG(CSLASK1), DREAL(CSLASK2) )
	        STQ = SQ*URB
	        SDQ = SQ/URB
	      ELSE
	        CQ = (1.0D0,0.0D0)
	        STQ = (0.0D0,0.0D0)
	        SDQ = DZOMEG
	      END IF
	      E1X = 1.0D0

	    END IF

	  END IF

	END IF

	RETURN
	END
