! copyright (c) 2008 Swedish Defence Research Agency (FOI) 
! may not be copied or distributed without written permission from FOI



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*16 DZ
	COMPLEX*32 OMEG,OMEG2,U,U2, UA2,UB2

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

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

C	choice of branch-cuts for ura and urb unimportant, the 	"minor" 
C	propagator elements are entire analytic functions of u
	URA = SQRT(UA2-U2)
	URB = SQRT(UB2-U2)
	URANORM = MAX ( ABS(REAL(URA)), ABS(AIMAG(URA)) )
	URBNORM = MAX ( ABS(REAL(URB)), ABS(AIMAG(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 := abs(aimag(p0)) + abs(aimag(q0))
C	  e2t := exp(-tlnref)
C	  cp,cq are not set
C	  stp := e2t * exp(+i*p0+i*q0)
C	  stq := e2t * exp(+i*p0-i*q0)
C	  sdp := e2t * exp(-i*p0+i*q0)
C	  sdq := e2t * exp(-i*p0-i*q0)

	  SLASK1 = ABS(AIMAG(P0))
	  SLASK2 = ABS(AIMAG(Q0))
	  TLNREF = SLASK1 + SLASK2
	  SLASK1 = EXP(-SLASK1)
	  SLASK2 = EXP(-SLASK2)
	  E2T = SLASK1 * SLASK2
	  SLASK = REAL(P0) + REAL(Q0)
	  CSLASK1 = QCMPLX ( COS(SLASK), SIN(SLASK) )
	  SLASK = REAL(P0) - REAL(Q0)
	  CSLASK2 = QCMPLX ( COS(SLASK), SIN(SLASK) )
	  IF (AIMAG(P0).GE.0.0Q0) THEN
	    IF (AIMAG(Q0).GE.0.0Q0) THEN
	      STP = (E2T*E2T) * CSLASK1
	      STQ = (SLASK1*SLASK1) * CSLASK2
	      SDP = (SLASK2*SLASK2) * CONJG(CSLASK2)
	      SDQ = CONJG(CSLASK1)
	    ELSE
	      STP = (SLASK1*SLASK1) * CSLASK1
	      STQ = (E2T*E2T) * CSLASK2
	      SDP = CONJG(CSLASK2)
	      SDQ = (SLASK2*SLASK2) * CONJG(CSLASK1)
	    END IF
	  ELSE
	    IF (AIMAG(Q0).GE.0.0Q0) THEN
	      STP = (SLASK2*SLASK2) * CSLASK1
	      STQ = CSLASK2
	      SDP = (E2T*E2T) * CONJG(CSLASK2)
	      SDQ = (SLASK1*SLASK1) * CONJG(CSLASK1)
	    ELSE
	      STP = CSLASK1
	      STQ = (SLASK2*SLASK2) * CSLASK2
	      SDP = (SLASK1*SLASK1) * CONJG(CSLASK2)
	      SDQ = (E2T*E2T) * CONJG(CSLASK1)
	    END IF
	  END IF

	ELSE

	  EXPTYP = .FALSE.
C	  < possibly   plnref := dabs(aimag(p0)) >
C	  < possibly   qlnref := dabs(aimag(q0)) >
C	  < e1p := exp(-plnref) >
C	  < e1q := exp(-qlnref) >
C	  tlnref := plnref + qlnref
C	  e2t := exp(-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 = ABS(AIMAG(P0))
	  IF ( URANORM.GT.exptyplim .OR. SLASK1.GT.exptyplim ) THEN
	    TLNREF = SLASK1                  ! temporary
	    SLASK1 = EXP(-SLASK1)
	    E2T = SLASK1                       ! temporary
	    SLASK = REAL(P0)
	    CSLASK = QCMPLX ( COS(SLASK), SIN(SLASK) )
	    IF (AIMAG(P0).GE.0.0Q0) THEN
	      CSLASK1 = (SLASK1*SLASK1) * CSLASK
	      CSLASK2 = CONJG(CSLASK)
	    ELSE
	      CSLASK1 = CSLASK
	      CSLASK2 = (SLASK1*SLASK1) * CONJG(CSLASK)
	    END IF
	    CP = 0.5Q0 * ( CSLASK1 + CSLASK2 )
	    CSLASK = CSLASK2 - CSLASK1
	    SP = 0.5Q0 * QCMPLX(-AIMAG(CSLASK),REAL(CSLASK))
	    STP = SP*URA
	    SDP = SP/URA
	  ELSE
	    TLNREF = 0.0Q0                   ! temporary
	    E2T = 1.0Q0                        ! temporary
	    IF (URA.NE.(0.0Q0,0.0Q0)) THEN
	      SLASK = REAL(P0)
	      CSLASK = QCMPLX ( COS(SLASK), SIN(SLASK) )
	      SLASK = AIMAG(P0)
	      CALL RHYPSMALL ( SLASK, SLASK1,SLASK2 )
	      CSLASK1 = SLASK1 * CSLASK
	      CSLASK2 = SLASK2 * CSLASK
	      CP = QCMPLX ( REAL(CSLASK1), -AIMAG(CSLASK2) )
	      SP = QCMPLX ( AIMAG(CSLASK1), REAL(CSLASK2) )
	      STP = SP*URA
	      SDP = SP/URA
	    ELSE
	      CP = (1.0Q0,0.0Q0)
	      STP = (0.0Q0,0.0Q0)
	      SDP = DZOMEG
	    END IF
	  END IF

	  SLASK2 = ABS(AIMAG(Q0))
	  IF ( URBNORM.GT.exptyplim .OR. SLASK2.GT.exptyplim ) THEN
	    TLNREF = TLNREF + SLASK2
	    SLASK2 = EXP(-SLASK2)
	    E2T = E2T * SLASK2
	    SLASK = REAL(Q0)
	    CSLASK = QCMPLX ( COS(SLASK), SIN(SLASK) )
	    IF (AIMAG(Q0).GE.0.0Q0) THEN
	      CSLASK1 = (SLASK2*SLASK2) * CSLASK
	      CSLASK2 = CONJG(CSLASK)
	    ELSE
	      CSLASK1 = CSLASK
	      CSLASK2 = (SLASK2*SLASK2) * CONJG(CSLASK)
	    END IF
	    CQ = 0.5Q0 * ( CSLASK1 + CSLASK2 )
	    CSLASK = CSLASK2 - CSLASK1
	    SQ = 0.5Q0 * QCMPLX(-AIMAG(CSLASK),REAL(CSLASK))
	    STQ = SQ*URB
	    SDQ = SQ/URB
	  ELSE
C	    tlnref and e2t unchanged
	    IF (URB.NE.(0.0Q0,0.0Q0)) THEN
	      SLASK = REAL(Q0)
	      CSLASK = QCMPLX ( COS(SLASK), SIN(SLASK) )
	      SLASK = AIMAG(Q0)
	      CALL RHYPSMALL ( SLASK, SLASK1,SLASK2 )
	      CSLASK1 = SLASK1 * CSLASK
	      CSLASK2 = SLASK2 * CSLASK
	      CQ = QCMPLX ( REAL(CSLASK1), -AIMAG(CSLASK2) )
	      SQ = QCMPLX ( AIMAG(CSLASK1), REAL(CSLASK2) )
	      STQ = SQ*URB
	      SDQ = SQ/URB
	    ELSE
	      CQ = (1.0Q0,0.0Q0)
	      STQ = (0.0Q0,0.0Q0)
	      SDQ = DZOMEG
	    END IF
	  END IF

	END IF

	RETURN
	END
