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



C	certain help routines for arithmetic operations with 
C	expressions of the form  dexp(e)*c

C       sven ivansson 1992





	SUBROUTINE SCALEEC ( alim1,alim2, E,C )

	IMPLICIT NONE

	REAL*8 alim1,alim2

	REAL*8 E
	COMPLEX*16 C

        INTEGER II
	REAL*8  SLASK1,SLASK2, SLASK

	SLASK1 = DREAL(C)
	SLASK2 = DIMAG(C)
        IF ( DABS(SLASK1) .GE. DABS(SLASK2) ) THEN
          II = 1
          SLASK = DABS(SLASK1)
        ELSE
          II = 2
          SLASK = DABS(SLASK2)
        END IF
	IF ( SLASK.LT.alim1 .OR. alim2.LT.SLASK ) THEN
	  IF (SLASK.GT.1.0D0) THEN
	    SLASK = 2/SLASK
	    E = E - DLOG(SLASK)
            IF (II.EQ.1) THEN
              IF (SLASK1.GE.0.0D0) THEN
                C =  2.0D0 * DCMPLX ( 1.0D0, SLASK2/SLASK1 )
              ELSE
                C = -2.0D0 * DCMPLX ( 1.0D0, SLASK2/SLASK1 )
              END IF
            ELSE
              IF (SLASK2.GE.0.0D0) THEN
                C =  2.0D0 * DCMPLX ( SLASK1/SLASK2, 1.0D0 )
              ELSE
                C = -2.0D0 * DCMPLX ( SLASK1/SLASK2, 1.0D0 )
              END IF
            END IF
	  ELSE IF (SLASK.GT.0.0D0) THEN
	    SLASK = 0.5D0/SLASK
	    E = E - DLOG(SLASK)
            IF (II.EQ.1) THEN
              IF (SLASK1.GE.0.0D0) THEN
                C =  0.5D0 * DCMPLX ( 1.0D0, SLASK2/SLASK1 )
              ELSE
                C = -0.5D0 * DCMPLX ( 1.0D0, SLASK2/SLASK1 )
              END IF
            ELSE
              IF (SLASK2.GE.0.0D0) THEN
                C =  0.5D0 * DCMPLX ( SLASK1/SLASK2, 1.0D0 )
              ELSE
                C = -0.5D0 * DCMPLX ( SLASK1/SLASK2, 1.0D0 )
              END IF
            END IF
	  ELSE
	    E = 0.0D0
	    C = (0.0D0,0.0D0)
	  END IF
	END IF

	RETURN
	END







	SUBROUTINE SCALE2EC ( alim1,alim2, E,C1,C2 )

	IMPLICIT NONE

	REAL*8 alim1,alim2

	REAL*8 E
	COMPLEX*16 C1,C2

	REAL*8 SLASK

	SLASK = DMAX1 ( DABS(DREAL(C1)), DABS(DIMAG(C1)), 
     $	DABS(DREAL(C2)), DABS(DIMAG(C2)) )
	IF ( SLASK.LT.alim1 .OR. alim2.LT.SLASK ) THEN
	  IF (SLASK.GT.1.0D0) THEN
	    SLASK = 2/SLASK
	    E = E - DLOG(SLASK)
	    C1 = SLASK*C1
	    C2 = SLASK*C2
	  ELSE IF (SLASK.GT.0.0D0) THEN
	    SLASK = 0.5D0/SLASK
	    E = E - DLOG(SLASK)
	    C1 = SLASK*C1
	    C2 = SLASK*C2
	  ELSE
	    E = 0.0D0
	    C1 = (0.0D0,0.0D0)
	    C2 = (0.0D0,0.0D0)
	  END IF
	END IF

	RETURN
	END





	SUBROUTINE SCALE3EC ( alim1,alim2, E,C1,C2,C3 )

	IMPLICIT NONE

	REAL*8 alim1,alim2

	REAL*8 E
	COMPLEX*16 C1,C2,C3

	REAL*8 SLASK

	SLASK = DMAX1 ( DABS(DREAL(C1)), DABS(DIMAG(C1)), 
     $	DABS(DREAL(C2)), DABS(DIMAG(C2)), DABS(DREAL(C3)), DABS(DIMAG(C3)) )
	IF ( SLASK.LT.alim1 .OR. alim2.LT.SLASK ) THEN
	  IF (SLASK.GT.1.0D0) THEN
	    SLASK = 2/SLASK
	    E = E - DLOG(SLASK)
	    C1 = SLASK*C1
	    C2 = SLASK*C2
	    C3 = SLASK*C3
	  ELSE IF (SLASK.GT.0.0D0) THEN
	    SLASK = 0.5D0/SLASK
	    E = E - DLOG(SLASK)
	    C1 = SLASK*C1
	    C2 = SLASK*C2
	    C3 = SLASK*C3
	  ELSE
	    E = 0.0D0
	    C1 = (0.0D0,0.0D0)
	    C2 = (0.0D0,0.0D0)
	    C3 = (0.0D0,0.0D0)
	  END IF
	END IF

	RETURN
	END




	SUBROUTINE SCALE4EC ( alim1,alim2, E,C1,C2,C3,C4 )

	IMPLICIT NONE

	REAL*8 alim1,alim2

	REAL*8 E
	COMPLEX*16 C1,C2,C3,C4

	REAL*8 SLASK

	SLASK = DMAX1 ( DABS(DREAL(C1)), DABS(DIMAG(C1)), 
     $	DABS(DREAL(C2)), DABS(DIMAG(C2)), DABS(DREAL(C3)), DABS(DIMAG(C3)), 
     $	DABS(DREAL(C4)), DABS(DIMAG(C4)) )
	IF ( SLASK.LT.alim1 .OR. alim2.LT.SLASK ) THEN
	  IF (SLASK.GT.1.0D0) THEN
	    SLASK = 2/SLASK
	    E = E - DLOG(SLASK)
	    C1 = SLASK*C1
	    C2 = SLASK*C2
	    C3 = SLASK*C3
	    C4 = SLASK*C4
	  ELSE IF (SLASK.GT.0.0D0) THEN
	    SLASK = 0.5D0/SLASK
	    E = E - DLOG(SLASK)
	    C1 = SLASK*C1
	    C2 = SLASK*C2
	    C3 = SLASK*C3
	    C4 = SLASK*C4
	  ELSE
	    E = 0.0D0
	    C1 = (0.0D0,0.0D0)
	    C2 = (0.0D0,0.0D0)
	    C3 = (0.0D0,0.0D0)
	    C4 = (0.0D0,0.0D0)
	  END IF
	END IF

	RETURN
	END





	SUBROUTINE SCALE5EC ( alim1,alim2, E,C1,C2,C3,C4,C6 )

	IMPLICIT NONE

	REAL*8 alim1,alim2

	REAL*8 E
	COMPLEX*16 C1,C2,C3,C4,C6

	REAL*8 SLASK

	SLASK = DMAX1 ( DABS(DREAL(C1)), DABS(DIMAG(C1)), 
     $	DABS(DREAL(C2)), DABS(DIMAG(C2)), DABS(DREAL(C3)), DABS(DIMAG(C3)), 
     $	DABS(DREAL(C4)), DABS(DIMAG(C4)), DABS(DREAL(C6)), DABS(DIMAG(C6)) )
	IF ( SLASK.LT.alim1 .OR. alim2.LT.SLASK ) THEN
	  IF (SLASK.GT.1.0D0) THEN
	    SLASK = 2/SLASK
	    E = E - DLOG(SLASK)
	    C1 = SLASK*C1
	    C2 = SLASK*C2
	    C3 = SLASK*C3
	    C4 = SLASK*C4
	    C6 = SLASK*C6
	  ELSE IF (SLASK.GT.0.0D0) THEN
	    SLASK = 0.5D0/SLASK
	    E = E - DLOG(SLASK)
	    C1 = SLASK*C1
	    C2 = SLASK*C2
	    C3 = SLASK*C3
	    C4 = SLASK*C4
	    C6 = SLASK*C6
	  ELSE
	    E = 0.0D0
	    C1 = (0.0D0,0.0D0)
	    C2 = (0.0D0,0.0D0)
	    C3 = (0.0D0,0.0D0)
	    C4 = (0.0D0,0.0D0)
	    C6 = (0.0D0,0.0D0)
	  END IF
	END IF

	RETURN
	END








	SUBROUTINE ADDEC ( EA,CA, EB,CB, E,C )
C	dexp(e)*c  will be  dexp(ea)*ca + dexp(eb)*cb
C	"overwrite" possible

	IMPLICIT NONE

	REAL*8 EA,EB
	COMPLEX*16 CA,CB

	REAL*8 E
	COMPLEX*16 C

	IF (EA.EQ.EB) THEN
	  E = EA
	  C = CA + CB
	ELSE IF (CA.EQ.(0.0D0,0.0D0)) THEN
	  E = EB
	  C = CB
	ELSE IF (CB.EQ.(0.0D0,0.0D0)) THEN
	  E = EA
	  C = CA
	ELSE IF (EB.GE.EA) THEN
	  C = DEXP(EA-EB)*CA + CB
	  E = EB
	ELSE
	  C = CA + DEXP(EB-EA)*CB
	  E = EA
	END IF

	RETURN
	END







	SUBROUTINE ADD2EC ( EA,C1A,C2A, EB,C1B,C2B, E,C1,C2 )
C	dexp(e)*c1  will be  dexp(ea)*c1a + dexp(eb)*c1b
C	dexp(e)*c2  will be  dexp(ea)*c2a + dexp(eb)*c2b
C	"overwrite" possible

	IMPLICIT NONE

	REAL*8 EA,EB
	COMPLEX*16 C1A,C2A, C1B,C2B

	REAL*8 E
	COMPLEX*16 C1,C2

	REAL*8 SLASK

	IF (EA.EQ.EB) THEN
	  E = EA
	  C1 = C1A + C1B
	  C2 = C2A + C2B
	  RETURN
	END IF

	IF (C1A.EQ.(0.0D0,0.0D0)) THEN
	  IF (C2A.EQ.(0.0D0,0.0D0)) THEN
	    E = EB
	    C1 = C1B
	    C2 = C2B
	    RETURN
	  END IF
	END IF
	IF (C1B.EQ.(0.0D0,0.0D0)) THEN
	  IF (C2B.EQ.(0.0D0,0.0D0)) THEN
	    E = EA
	    C1 = C1A
	    C2 = C2A
	    RETURN
	  END IF
	END IF

	IF (EB.GE.EA) THEN
	  SLASK = DEXP(EA-EB)
	  C1 = SLASK*C1A + C1B
	  C2 = SLASK*C2A + C2B
	  E = EB
	ELSE
	  SLASK = DEXP(EB-EA)
	  C1 = C1A + SLASK*C1B
	  C2 = C2A + SLASK*C2B
	  E = EA
	END IF

	RETURN
	END





	SUBROUTINE ADD4EC ( EA,C1A,C2A,C3A,C4A, EB,C1B,C2B,
     $	C3B,C4B, E,C1,C2,C3,C4 )
C	dexp(e)*c1  will be  dexp(ea)*c1a + dexp(eb)*c1b
C	dexp(e)*c2  will be  dexp(ea)*c2a + dexp(eb)*c2b
C	dexp(e)*c3  will be  dexp(ea)*c3a + dexp(eb)*c3b
C	dexp(e)*c4  will be  dexp(ea)*c4a + dexp(eb)*c4b
C	"overwrite" possible

	IMPLICIT NONE

	REAL*8 EA,EB
	COMPLEX*16 C1A,C2A,C3A,C4A, C1B,C2B,C3B,C4B

	REAL*8 E
	COMPLEX*16 C1,C2,C3,C4

	REAL*8 SLASK

	IF (EA.EQ.EB) THEN
	  E = EA
	  C1 = C1A + C1B
	  C2 = C2A + C2B
	  C3 = C3A + C3B
	  C4 = C4A + C4B
	  RETURN
	END IF

	IF (C1A.EQ.(0.0D0,0.0D0)) THEN
	  IF (C2A.EQ.(0.0D0,0.0D0)) THEN
	    IF (C3A.EQ.(0.0D0,0.0D0)) THEN
	      IF (C4A.EQ.(0.0D0,0.0D0)) THEN
	        E = EB
	        C1 = C1B
	        C2 = C2B
	        C3 = C3B
	        C4 = C4B
	        RETURN
	      END IF
	    END IF
	  END IF
	END IF
	IF (C1B.EQ.(0.0D0,0.0D0)) THEN
	  IF (C2B.EQ.(0.0D0,0.0D0)) THEN
	    IF (C3B.EQ.(0.0D0,0.0D0)) THEN
	      IF (C4B.EQ.(0.0D0,0.0D0)) THEN
	        E = EA
	        C1 = C1A
	        C2 = C2A
	        C3 = C3A
	        C4 = C4A
	        RETURN
	      END IF
	    END IF
	  END IF
	END IF

	IF (EB.GE.EA) THEN
	  SLASK = DEXP(EA-EB)
	  C1 = SLASK*C1A + C1B
	  C2 = SLASK*C2A + C2B
	  C3 = SLASK*C3A + C3B
	  C4 = SLASK*C4A + C4B
	  E = EB
	ELSE
	  SLASK = DEXP(EB-EA)
	  C1 = C1A + SLASK*C1B
	  C2 = C2A + SLASK*C2B
	  C3 = C3A + SLASK*C3B
	  C4 = C4A + SLASK*C4B
	  E = EA
	END IF

	RETURN
	END






	SUBROUTINE EQ2EC ( E1,C1, E2,C2 )
C	equalize the exponents e1,e2

	IMPLICIT NONE

	REAL*8 E1,E2
	COMPLEX*16 C1,C2

	IF (E1.EQ.E2) RETURN

	IF (C1.EQ.(0.0D0,0.0D0)) THEN
	  E1 = E2
	ELSE IF (C2.EQ.(0.0D0,0.0D0)) THEN
	  E2 = E1
	ELSE IF (E2.GE.E1) THEN
	  C1 = DEXP(E1-E2)*C1
	  E1 = E2
	ELSE
	  C2 = DEXP(E2-E1)*C2
	  E2 = E1
	END IF

	RETURN
	END



	SUBROUTINE EQ22EC ( E1,C11,C12, E2,C21,C22 )
C	equalize the exponents e1 (for c11,c12) and e2 (for c21,c22)

	IMPLICIT NONE

	REAL*8 E1, E2
	COMPLEX*16 C11,C12, C21,C22

	REAL*8 SLASK

	IF (E1.EQ.E2) RETURN

	IF (C11.EQ.(0.0D0,0.0D0)) THEN
	  IF (C12.EQ.(0.0D0,0.0D0)) THEN
	    E1 = E2
	    RETURN
	  END IF
	END IF
	IF (C21.EQ.(0.0D0,0.0D0)) THEN
	  IF (C22.EQ.(0.0D0,0.0D0)) THEN
	    E2 = E1
	    RETURN
	  END IF
	END IF

	IF (E2.GE.E1) THEN
	  SLASK = DEXP(E1-E2)
	  C11 = SLASK*C11
	  C12 = SLASK*C12
	  E1 = E2
	ELSE
	  SLASK = DEXP(E2-E1)
	  C21 = SLASK*C21
	  C22 = SLASK*C22
	  E2 = E1
	END IF

	RETURN
	END





	SUBROUTINE EQ3EC ( E1,C1, E2,C2, E3,C3 )
C	equalize the exponents e1,e2,e3

	IMPLICIT NONE

	REAL*8 E1,E2,E3
	COMPLEX*16 C1,C2,C3

	IF (E1.EQ.E2) THEN
	  IF (E2.EQ.E3) RETURN
	END IF

	IF (C1.EQ.(0.0D0,0.0D0)) THEN
	  CALL EQ2EC ( E2,C2, E3,C3 )
	  E1 = E2   ! .eq.e3
	ELSE IF (C2.EQ.(0.0D0,0.0D0)) THEN
	  CALL EQ2EC ( E1,C1, E3,C3 )
	  E2 = E1   ! .eq.e3
	ELSE IF (C3.EQ.(0.0D0,0.0D0)) THEN
	  CALL EQ2EC ( E1,C1, E2,C2 )
	  E3 = E1   ! .eq.e2
	ELSE IF (E1.GE.E2) THEN
	  IF (E1.GE.E3) THEN
	    GOTO 1
	  ELSE
	    GOTO 3
	  END IF
	ELSE
	  IF (E2.GE.E3) THEN
	    GOTO 2
	  ELSE
	    GOTO 3
	  END IF
	END IF

 1	C2 = DEXP(E2-E1)*C2
	E2 = E1
	C3 = DEXP(E3-E1)*C3
	E3 = E1
	RETURN

 2	C1 = DEXP(E1-E2)*C1
	E1 = E2
	C3 = DEXP(E3-E2)*C3
	E3 = E2
	RETURN

 3	C1 = DEXP(E1-E3)*C1
	E1 = E3
	C2 = DEXP(E2-E3)*C2
	E2 = E3
	RETURN

	END






	SUBROUTINE EQ4EC ( E1,C1, E2,C2, E3,C3, E4,C4 )
C	equalize the exponents e1,e2,e3,e4

	IMPLICIT NONE

	REAL*8 E1,E2,E3,E4
	COMPLEX*16 C1,C2,C3,C4

	IF (E1.EQ.E2) THEN
	  IF (E2.EQ.E3) THEN
	    IF (E3.EQ.E4) RETURN
	  END IF
	END IF

	IF (C1.EQ.(0.0D0,0.0D0)) THEN
	  CALL EQ3EC ( E2,C2, E3,C3, E4,C4 )
	  E1 = E2   ! .eq.e3.eq.e4
	ELSE IF (C2.EQ.(0.0D0,0.0D0)) THEN
	  CALL EQ3EC ( E1,C1, E3,C3, E4,C4 )
	  E2 = E1   ! .eq.e3.eq.e4
	ELSE IF (C3.EQ.(0.0D0,0.0D0)) THEN
	  CALL EQ3EC ( E1,C1, E2,C2, E4,C4 )
	  E3 = E1   ! .eq.e2.eq.e4
	ELSE IF (C4.EQ.(0.0D0,0.0D0)) THEN
	  CALL EQ3EC ( E1,C1, E2,C2, E3,C3 )
	  E4 = E1   ! .eq.e2.eq.e3
	ELSE IF (E1.GE.E2) THEN
	  IF (E1.GE.E3) THEN
	    IF (E1.GE.E4) THEN
	      GOTO 1
	    ELSE
	      GOTO 4
	    END IF
	  ELSE
	    IF (E3.GE.E4) THEN
	      GOTO 3
	    ELSE
	      GOTO 4
	    END IF
	  END IF
	ELSE
	  IF (E2.GE.E3) THEN
	    IF (E2.GE.E4) THEN
	      GOTO 2
	    ELSE
	      GOTO 4
	    END IF
	  ELSE
	    IF (E3.GE.E4) THEN
	      GOTO 3
	    ELSE
	      GOTO 4
	    END IF
	  END IF
	END IF

 1	C2 = DEXP(E2-E1)*C2
	E2 = E1
	C3 = DEXP(E3-E1)*C3
	E3 = E1
	C4 = DEXP(E4-E1)*C4
	E4 = E1
	RETURN

 2	C1 = DEXP(E1-E2)*C1
	E1 = E2
	C3 = DEXP(E3-E2)*C3
	E3 = E2
	C4 = DEXP(E4-E2)*C4
	E4 = E2
	RETURN

 3	C1 = DEXP(E1-E3)*C1
	E1 = E3
	C2 = DEXP(E2-E3)*C2
	E2 = E3
	C4 = DEXP(E4-E3)*C4
	E4 = E3
	RETURN

 4	C1 = DEXP(E1-E4)*C1
	E1 = E4
	C2 = DEXP(E2-E4)*C2
	E2 = E4
	C3 = DEXP(E3-E4)*C3
	E3 = E4
	RETURN

	END




	SUBROUTINE EQ24EC ( E1,C11,C12,C13,C14, 
     $	E2,C21,C22,C23,C24 )
C	equalize the exponents e1 (for c11,c12,c13,c14) and 
c	e2 (for c21,c22,c23,c24)

	IMPLICIT NONE

	REAL*8 E1, E2
	COMPLEX*16 C11,C12,C13,C14, C21,C22,C23,C24

	REAL*8 SLASK

	IF (E1.EQ.E2) RETURN

	IF (C11.EQ.(0.0D0,0.0D0)) THEN
	  IF (C12.EQ.(0.0D0,0.0D0)) THEN
	    IF (C13.EQ.(0.0D0,0.0D0)) THEN
	      IF (C14.EQ.(0.0D0,0.0D0)) THEN
	        E1 = E2
	        RETURN
	      END IF
	    END IF
	  END IF
	END IF
	IF (C21.EQ.(0.0D0,0.0D0)) THEN
	  IF (C22.EQ.(0.0D0,0.0D0)) THEN
	    IF (C23.EQ.(0.0D0,0.0D0)) THEN
	      IF (C24.EQ.(0.0D0,0.0D0)) THEN
	        E2 = E1
	        RETURN
	      END IF
	    END IF
	  END IF
	END IF

	IF (E2.GE.E1) THEN
	  SLASK = DEXP(E1-E2)
	  C11 = SLASK*C11
	  C12 = SLASK*C12
	  C13 = SLASK*C13
	  C14 = SLASK*C14
	  E1 = E2
	ELSE
	  SLASK = DEXP(E2-E1)
	  C21 = SLASK*C21
	  C22 = SLASK*C22
	  C23 = SLASK*C23
	  C24 = SLASK*C24
	  E2 = E1
	END IF

	RETURN
	END
