C	subroutine to be used in harksi (related to "solsipart for solsi")

C	Sven Ivansson  1990





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

	IMPLICIT NONE

	INCLUDE 'RPPARAMETER.INC'

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

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

	REAL*8 URANORM,URBNORM, 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
	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

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

	  EXPTYP = .TRUE.
C	  tlnref := dabs(dimag(p0)) + dabs(dimag(q0))
C	  e2t := dexp(-tlnref)
C	  cp,cq are not set
C	  stp := e2t * cdexp(+i*p0+i*q0)
C	  stq := e2t * cdexp(+i*p0-i*q0)
C	  sdp := e2t * cdexp(-i*p0+i*q0)
C	  sdq := e2t * cdexp(-i*p0-i*q0)

	  SLASK1 = DABS(DIMAG(P0))
	  SLASK2 = DABS(DIMAG(Q0))
	  TLNREF = SLASK1 + SLASK2
	  SLASK1 = DEXP(-SLASK1)
	  SLASK2 = DEXP(-SLASK2)
	  E2T = SLASK1 * SLASK2
	  SLASK = DREAL(P0) + DREAL(Q0)
	  CSLASK1 = DCMPLX ( DCOS(SLASK), DSIN(SLASK) )
	  SLASK = DREAL(P0) - DREAL(Q0)
	  CSLASK2 = DCMPLX ( DCOS(SLASK), DSIN(SLASK) )
	  IF (DIMAG(P0).GE.0.0D0) THEN
	    IF (DIMAG(Q0).GE.0.0D0) THEN
	      STP = (E2T*E2T) * CSLASK1
	      STQ = (SLASK1*SLASK1) * CSLASK2
	      SDP = (SLASK2*SLASK2) * DCONJG(CSLASK2)
	      SDQ = DCONJG(CSLASK1)
	    ELSE
	      STP = (SLASK1*SLASK1) * CSLASK1
	      STQ = (E2T*E2T) * CSLASK2
	      SDP = DCONJG(CSLASK2)
	      SDQ = (SLASK2*SLASK2) * DCONJG(CSLASK1)
	    END IF
	  ELSE
	    IF (DIMAG(Q0).GE.0.0D0) THEN
	      STP = (SLASK2*SLASK2) * CSLASK1
	      STQ = CSLASK2
	      SDP = (E2T*E2T) * DCONJG(CSLASK2)
	      SDQ = (SLASK1*SLASK1) * DCONJG(CSLASK1)
	    ELSE
	      STP = CSLASK1
	      STQ = (SLASK2*SLASK2) * CSLASK2
	      SDP = (SLASK1*SLASK1) * DCONJG(CSLASK2)
	      SDQ = (E2T*E2T) * DCONJG(CSLASK1)
	    END IF
	  END IF

	ELSE

	  EXPTYP = .FALSE.
C	  < possibly   plnref := dabs(dimag(p0)) >
C	  < possibly   qlnref := dabs(dimag(q0)) >
C	  < e1p := dexp(-plnref) >
C	  < e1q := dexp(-qlnref) >
C	  tlnref := plnref + qlnref
C	  e2t := dexp(-tlnref)
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

	  SLASK1 = DABS(DIMAG(P0))
	  IF ( URANORM.GT.exptyplim .OR. SLASK1.GT.exptyplim ) THEN
	    TLNREF = SLASK1                  ! temporary
	    SLASK1 = DEXP(-SLASK1)
	    E2T = SLASK1                       ! temporary
	    SLASK = DREAL(P0)
	    CSLASK = DCMPLX ( DCOS(SLASK), DSIN(SLASK) )
	    IF (DIMAG(P0).GE.0.0D0) THEN
	      CSLASK1 = (SLASK1*SLASK1) * CSLASK
	      CSLASK2 = DCONJG(CSLASK)
	    ELSE
	      CSLASK1 = CSLASK
	      CSLASK2 = (SLASK1*SLASK1) * DCONJG(CSLASK)
	    END IF
	    CP = 0.5D0 * ( CSLASK1 + CSLASK2 )
	    CSLASK = CSLASK2 - CSLASK1
	    SP = 0.5D0 * DCMPLX(-DIMAG(CSLASK),DREAL(CSLASK))
	    STP = SP*URA
	    SDP = SP/URA
	  ELSE
	    TLNREF = 0.0D0                   ! temporary
	    E2T = 1.0D0                        ! temporary
	    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
	  END IF

	  SLASK2 = DABS(DIMAG(Q0))
	  IF ( URBNORM.GT.exptyplim .OR. SLASK2.GT.exptyplim ) THEN
	    TLNREF = TLNREF + SLASK2
	    SLASK2 = DEXP(-SLASK2)
	    E2T = E2T * SLASK2
	    SLASK = DREAL(Q0)
	    CSLASK = DCMPLX ( DCOS(SLASK), DSIN(SLASK) )
	    IF (DIMAG(Q0).GE.0.0D0) THEN
	      CSLASK1 = (SLASK2*SLASK2) * CSLASK
	      CSLASK2 = DCONJG(CSLASK)
	    ELSE
	      CSLASK1 = CSLASK
	      CSLASK2 = (SLASK2*SLASK2) * DCONJG(CSLASK)
	    END IF
	    CQ = 0.5D0 * ( CSLASK1 + CSLASK2 )
	    CSLASK = CSLASK2 - CSLASK1
	    SQ = 0.5D0 * DCMPLX(-DIMAG(CSLASK),DREAL(CSLASK))
	    STQ = SQ*URB
	    SDQ = SQ/URB
	  ELSE
C	    tlnref and e2t unchanged
	    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
	  END IF

	END IF

	RETURN
	END