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



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
        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 (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   
C              diag(1,rho) * ... * diag(1,1/rho)

	IMPLICIT NONE

	INCLUDE 'AIRYP.INC'  ! to get ai4e

	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, YROT, CROT, C120,C120N, 
     $	ABA(4),ABB(4), CC(8)

	SAVE CROT, EA,ABA, EB,ABB

	DATA C120  /(-0.5D0,+0.8660254037844387D0)/    ! cdexp(+i*2*pi/3)
	DATA C120N /(-0.5D0,-0.8660254037844387D0)/    ! cdexp(-i*2*pi/3)

	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) )
	  IF (AI4E) THEN
	    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
	    CALL AIRY2E (YA, ALN, ABA(1),ABA(2))
	    EA(1) = -ALN
	    EA(2) = -ALN
	    IF ( DIMAG(YA) .GE. 0.0D0 ) THEN
	      CROT = C120N
	    ELSE
	      CROT = C120
	    END IF
	    YROT = YA*CROT
	    CALL AIRY2E (YROT, ALN, ABA(3),ABA(4))
	    ABA(4) = CROT*ABA(4)
	    EA(3) = -ALN
	    EA(4) = -ALN
	    CALL AIRY2E (YB, ALN, ABB(1),ABB(2))
	    EB(1) = -ALN
	    EB(2) = -ALN
	    YROT = YB*CROT
	    CALL AIRY2E (YROT, ALN, ABB(3),ABB(4))
	    ABB(4) = CROT*ABB(4)
	    EB(3) = -ALN
	    EB(4) = -ALN
	  END IF
	ELSE
	  IF (IOTYP.GT.0) THEN
	    YA = YB
	    YB = CFAKT * ( RK2 - OMEG2*(UA2+0.5D0*UB2) )
	    IF (AI4E) THEN
	      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
	      DO I = 1, 2
	        EA(I) = EB(I)
	        ABA(I) = ABB(I)
	      END DO
	      CSLASK = CROT
	      IF ( DIMAG(YA) .GE. 0.0D0 ) THEN
	        CROT = C120N
	      ELSE
	        CROT = C120
	      END IF
	      IF (CROT.EQ.CSLASK) THEN
	        DO I = 3, 4
	          EA(I) = EB(I)
	          ABA(I) = ABB(I)
	        END DO
	      ELSE
	        YROT = YA*CROT
	        CALL AIRY2E (YROT, ALN, ABA(3),ABA(4))
	        ABA(4) = CROT*ABA(4)
	        EA(3) = -ALN
	        EA(4) = -ALN
	      END IF
	      CALL AIRY2E (YB, ALN, ABB(1),ABB(2))
	      EB(1) = -ALN
	      EB(2) = -ALN
	      YROT = YB*CROT
	      CALL AIRY2E (YROT, ALN, ABB(3),ABB(4))
	      ABB(4) = CROT*ABB(4)
	      EB(3) = -ALN
	      EB(4) = -ALN
	    END IF
	  ELSE
	    YB = YA
	    YA = CFAKT * ( RK2 - OMEG2*(UA2-0.5D0*UB2) )
	    IF (AI4E) THEN
	      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
	    ELSE
	      DO I = 1, 2
	        EB(I) = EA(I)
	        ABB(I) = ABA(I)
	      END DO
	      CSLASK = CROT
	      IF ( DIMAG(YA) .GE. 0.0D0 ) THEN
	        CROT = C120N
	      ELSE
	        CROT = C120
	      END IF
	      IF (CROT.EQ.CSLASK) THEN
	        DO I = 3, 4
	          EB(I) = EA(I)
	          ABB(I) = ABA(I)
	        END DO
	      ELSE
	        YROT = YB*CROT
	        CALL AIRY2E (YROT, ALN, ABB(3),ABB(4))
	        ABB(4) = CROT*ABB(4)
	        EB(3) = -ALN
	        EB(4) = -ALN
	      END IF
	      CALL AIRY2E (YA, ALN, ABA(1),ABA(2))
	      EA(1) = -ALN
	      EA(2) = -ALN
	      YROT = YA*CROT
	      CALL AIRY2E (YROT, ALN, ABA(3),ABA(4))
	      ABA(4) = CROT*ABA(4)
	      EA(3) = -ALN
	      EA(4) = -ALN
	    END IF
	  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

	IF (.NOT.AI4E) THEN
	  IF ( DIMAG(YA) .GE. 0.0D0 ) THEN
	    CSLASK = DCMPLX ( DSQRT(3.0D0), -1.0D0 )  ! hmf 10.4.12
	  ELSE
	    CSLASK = DCMPLX ( DSQRT(3.0D0), +1.0D0 )  ! hmf 10.4.11
	  END IF
	  CP  = CSLASK * CP
	  CPX = CSLASK * CPX
	  STP = CSLASK * STP
	  SDP = CSLASK * SDP
	END IF

	RETURN
	END










	SUBROUTINE ROT3 (YA,YB,CFAKT,BCFAKT,BCFAKTINV)
C	rotate +-2*pi/3 (or 0) , in order to maximize min(alna,alnb) from 
C	" call airy4e (ya, alna,blna, aiea,aidea,biea,bidea) "
C	" call airy4e (yb, alnb,blnb, aieb,aideb,bieb,bideb) "

	IMPLICIT NONE

	COMPLEX*16 YA,YB,CFAKT,BCFAKT,BCFAKTINV

	INTEGER II
	REAL*8 ALN(2)
	COMPLEX*16 Y(2),CZET(2), CSLASK1,CSLASK2, C120,C120N

	DATA C120  /(-0.5D0,+0.8660254037844387D0)/    ! cdexp(+i*2*pi/3)
	DATA C120N /(-0.5D0,-0.8660254037844387D0)/    ! cdexp(-i*2*pi/3)

	IF ( DIMAG(YA).NE.0.0D0 .OR. DIMAG(YB).NE.0.0D0 ) THEN
	  Y(1) = YA
	  Y(2) = YB
	  CZET(1) = Y(1) * CDSQRT(Y(1))
	  CZET(2) = Y(2) * CDSQRT(Y(2))
	  ALN(1) = DREAL(CZET(1))
	  ALN(2) = DREAL(CZET(2))
	  IF ( DABS(ALN(1)) .GE. DABS(ALN(2)) ) THEN
	    II = 1
	  ELSE
	    II = 2
	  END IF
	  IF ( ALN(II) .LT. 0.0D0 ) THEN
C	    this is when we want to rotate
	    IF ( DIMAG(Y(II)) .GE. 0.0D0 ) THEN
	      CSLASK1 = C120N
	      CSLASK2 = C120
	    ELSE
	      CSLASK1 = C120
	      CSLASK2 = C120N
	    END IF
	    YA = CSLASK1 * YA
	    YB = CSLASK1 * YB
	    CFAKT = CSLASK1 * CFAKT
	    BCFAKT = CSLASK1 * BCFAKT
	    BCFAKTINV = CSLASK2 * BCFAKTINV
	  END IF
	END IF

	RETURN
	END
