C     PROGRAM FOR CALCULATION OF DOUBLE AUGER
C      EFFECT PROBABILITY
C      L.V.CHERNYSHEVA,M.YA.AMUSIA, V.A.KILIN, I.S.LI
C      *NEW VERSION  July 1992 with V.A.Kilin
C============================================
C                VARIABLES
C      ND- NUMBER OF DIAGRAMS
C      NES- NUMBER OF EXCITED STATES
C      AA - INITIAL MATRIX 1
C      Y -  INITIAL MATRIX 2
C      NORM- NORM initial**(-2) * NORM final**( -2)
C      *EXPI=1 (INPUT KK,IND,DLAM,GAM)
C      *KK- NUMBER OF POINTS IN IND,DLAM,GAM -ARRAYS
C      *IND- NUMBER OF FUNCTIONS HAVING DLAM
C      *DLAM- ENERGY DENOMINATOR
C      *GAM-  INTERMEDIATE STATE WIDTH
C      NE5= NE+15, NE - NUMBER OF POINTS IN RADIAL MESH
C      S - NUMBER OF SHELLS
C      NW1 - MAX NUMBER OF FUNCTIONS IN FILE
C     ======================================
      SUBROUTINE  DAUGER (ND,NES,AA,Y,NORM,EXPI,KK,IND,DLAM,GAM,
     * NE,S,NW1,LI,LAMM,LAMI,
     * ME,ME11,ME12,ME13,ME14,ME15,
     * ME16,ME17,ME18,ME19,
     * CC,QQ,DD1,RR,DLAM1,GAM1,R1,R2,R5,R6,R7,R8)

      INTEGER NE,NE5,S,ND,LV1,LV2,M,M1,I1,I2,I3,I4,I5,I6,
     *I,J,K,L,A,B,C,D,NW1,EXPI,KK,IND(KK),
     *AA(12,ND),LI(S),LX(7)

      REAL*8 Z,R,BET,ALFA,H,OM,OM1,OM2,E,SUM1,SUM2,SUM3,SUM4,
     *SUMI1,SUMI2,SUMI3,SUMI4,DELE,
     *RES1,RES2,RO1,NORM,
     *RDD,RDE,RED,REE,IDD,IDE,IED,IEE,
     *DLAM(KK),DLAM1((NES+1)*NW1),GAM(KK),GAM1((NES+1)*NW1),
     *Y(5,ND),R1(NE+15),R2(NE+15),R5(NE+15),R6(NE+15),
     * R7(NE+15),R8(NE+15,S+NW1*ND),
     *ME4(1),ME5(1),ME6(1),ME7(1),ME8(1),ME(NW1,*),
     *ME11(NW1),ME12(NW1),ME13(NW1),
     *ME14(NW1),ME15(NW1),ME16(NW1),
     *ME17(NW1),ME18(NW1),ME19(NW1),
     *PR(20),CC(NW1),QQ(NW1),DD1(NW1),RR(NW1),
     *LAMX(7),LAMM(NW1),LAMI(S),DAT(15)

      PRINT 10
  10  FORMAT ('   AMPLITUDES OF THE DOUBLE AUGER PROCESS'/)
      READ(11) (DAT(I),I=1,15)
C       WRITE(6,12)(DAT(J),J=1,15)
C 12  FORMAT(5E15.7)
      Z=DAT(7)
      NE5=NE+15
      R=DAT(2)
      H=DAT(4)
      BET=DAT(5)
       CALL VAR11(NE5,Z,R,BET,H,ALFA,RO1,R1)
      PRINT 20,Z,S,NES,NORM,NE,R,H,BET,ALFA
  20  FORMAT('  CHARGE OF NUCLEUS   ',E15.8/
     *       '  NUMBER OF SHELLS    ',I5/
     *       '  NO OF EXCITED STATES',I5/
     *       '  NORMA               ',E15.8/
     *       '  NUMBER OF POINTS    ',I5/
     *       '  R MAX               ',E15.8/
     *       '  STEP OF INTEGRATION ',E15.8/
     *       '  BETTA               ',E15.8/
     *       '  ALFA                ',E15.8)
        IF ( (DAT(8)-S).GT.1.0 ) GOTO 1000
       PRINT 30
  30   FORMAT('  INITIAL MATRICES')
       PRINT 40
  40  FORMAT('  WAVE FUNCTIONS IN ZONES      M  M1 LV1 LV2 KSU1 ')
        WRITE(6,41)((AA(J,I),J=1,12),I=1,ND)
  41  FORMAT(12I4)
      PRINT 50
  50  FORMAT('     ANGLE     SPIN DD     SPIN DE     SPIN ED',
     *'    SPIN EE  ')
        WRITE(6,60)((Y(J,I),J=1,5),I=1,ND)
  60  FORMAT (5F11.5)
      DO 61 I=1,ND
      Y(1,I)=Y(1,I)/DSQRT(NORM)
  61  CONTINUE
C==========================================
      PRINT 70
  70  FORMAT(' == GROUND STATE CONFIGURATION'/
     *       '  ORBITAL  L      ENERGY ',
     *  ' IN ZONE   EXP ENERGY     GAMMA  ')
      DO 80 I=1, S
         DLAM1(I)=0
         GAM1(I)=0
         IF (EXPI.EQ.1) THEN
         DO 81 II=1,KK
         IF (I.EQ.IND(II)) THEN
            DLAM1(I)=DLAM(II)
            GAM1(I)=GAM(II)
         END IF
 81   CONTINUE
         END IF
      CALL DISKO1( 9,I,2,NE5,R8(1,I))
       LAMI(I)=R8(1,I)
       LI(I)=R8(2,I)
        IF (EXPI.EQ.0) THEN
        WRITE(6,90) I,LI(I),LAMI(I),I
        ELSE
        PRINT 91,I,LI(I),LAMI(I),I,DLAM1(I),GAM1(I)
        END IF
  80  CONTINUE
  90  FORMAT(I6,I6,E15.8,I6)
  91  FORMAT(2I6,E15.8,I6,2E15.8)
C============================================
C     GET WAVEFUNCTIONS OF EXCITED STATES
       A=S
       DO 120 I=1,NES
       LNO=11+I-1
       REWIND LNO
      READ(LNO) (DAT(J),J=1,15)
      IF (DAT(1).EQ.3) THEN
      M=DAT(13)+DAT(12)
      ELSE
      M=DAT(8)
      END IF
      PRINT 121,I
  121 FORMAT('======================='/
     *    '   EXCITED STATES',I6)
      PRINT 126
  126 FORMAT('      ENERGIES       PHASE OR N ',
     *       '  IN ZONES   EXP.ENERGY      GAMMA ')
      DO 130 K=1, M
       LA=A+K
          DLAM1(LA)=0
          GAM1(LA)=0
      IF(EXPI.EQ.1) THEN
       DO 131 II=1,KK
C      PRINT 132,II,DLAM(II)
C132   FORMAT(' II,dLAm  ',I6,E15.8)
       IF (LA.EQ.IND(II)) THEN
          DLAM1(LA)=DLAM(II)
          GAM1(LA)=GAM(II)
       END IF
 131  CONTINUE
      END IF
      CALL DISKO1(LNO,K,2,NE5,R8(1,LA))
       IA=7
       IF(R8(1,LA).LT.0.)IA=5
          IF (EXPI.EQ.0) THEN
       WRITE(6,140) K,R8(1,LA),R8(NE+IA,LA),LA
          ELSE
       WRITE(6,141) K,R8(1,LA),R8(NE+IA,LA),LA,DLAM1(LA),GAM1(LA)
          END IF
  130  CONTINUE
       A=LA
  140  FORMAT(I4,2E15.8,I6)
  141  FORMAT(I4,2E15.8,I6,2E15.8)
       PRINT 150,R8(2,LA)
  150  FORMAT('  WITH ORBITAL  MOMENTUM ',F4.1)
  120  CONTINUE
C===========================================
C      calculations OF DIAGRAMS" CONTRIBUTIONS
       RES1=0
       RES2=0
      DO 300 I=1, ND
       I1=AA(1,I)
       I2=AA(2,I)
       I3=AA(3,I)
       I4=AA(4,I)
       I5=AA(5,I)
       I6=AA(6,I)
       M=AA(8,I)
       M1=AA(9,I)
       PRINT 310,I
  310  FORMAT('=========================================='/
     * ' DIAGRAM NO',I6/'  ORBITALS   L    ENERGIES',
     *'   ZONES')
      DO 320 J=1, 7
        A=AA(J,I)
       CALL COLUMN(NE5,R8(1,A),R5)
       LX(J)=R5(2)
       LAMX(J)=R5(1)
       PRINT 330,J,LX(J),LAMX(J),A
  320  CONTINUE
  330  FORMAT(2I7,E15.8,I7)
       LV1=AA(10,I)
       LV2=AA(11,I)
      PRINT 340
  340 FORMAT('  TRANSFERED MOMENTA   LV1     LV2')
      PRINT 350,LV1,LV2
  350 FORMAT('                  ',2I6)
       L=LX(7)
       PRINT 360,L
  360 FORMAT('  INTERMEDIATE STATE  ORBITAL MOMENTUM=',I4)
       OM1=LAMX(1)+LAMX(2)-LAMX(4)
       SUM1=0
       SUM2=0
       SUM3=0
       SUM4=0
       SUMI1=0
       SUMI2=0
       SUMI3=0
       SUMI4=0
      DO 380 J=1, S
      IF ( L.EQ.LI(J) )  THEN
       E=LAMI(J)
       OM=(OM1-E)
       IF(EXPI.EQ.1 .AND.( (DLAM1(J).NE.0.).OR.GAM1(J).NE.0.)) THEN
       DELE=DLAM1(J)-OM
       ELSE
       DELE=0
       END IF
       DELE=DELE*AA(12,I)
       OM=OM+DELE
       PRINT 381,DELE
 381   FORMAT('  delE = ',E15.8)  
      CALL MAEL5(ALFA,BET,NE5,H,NW1,I2,I1,I4,J,1,LX(2),LX(1),LX(4),
     *L,LV1,R1,PR,ME,ME4,ME5,ME6,ME7,R8,R2,R5,R6,R7)
      CALL MAEL5(ALFA,BET,NE5,H,NW1,I5,I6,I3,J,1,LX(5),LX(6),LX(3),
     *L,LV2,R1,PR,ME,ME7,ME8,ME6,ME6,R8,R2,R5,R6,R7)
       RDD=ME4(1)*ME7(1)*Y(1,I)
       RDE=ME4(1)*ME8(1)*Y(1,I)
       RED=ME5(1)*ME7(1)*Y(1,I)
       REE=ME5(1)*ME8(1)*Y(1,I)
       IDD=-RDD*GAM1(J)/2/(OM**2+GAM1(J)**2/4)
       IDE=-RDE*GAM1(J)/2/(OM**2+GAM1(J)**2/4)
       IED=-RED*GAM1(J)/2/(OM**2+GAM1(J)**2/4)
       IEE=-REE*GAM1(J)/2/(OM**2+GAM1(J)**2/4)
       RDD=RDD*OM/(OM**2+GAM1(J)**2/4)
       RDE=RDE*OM/(OM**2+GAM1(J)**2/4)
       RED=RED*OM/(OM**2+GAM1(J)**2/4)
       REE=REE*OM/(OM**2+GAM1(J)**2/4)
       PRINT 385,J,E,OM,ME4,ME5,ME7,ME8,RDD,RDE,RED,REE,
     *  IDD,IDE,IED,IEE
       SUM1=SUM1+RDD
       SUM2=SUM2+RDE
       SUM3=SUM3+RED
       SUM4=SUM4+REE
       SUMI1=SUMI1+IDD
       SUMI2=SUMI2+IDE
       SUMI3=SUMI3+IED
       SUMI4=SUMI4+IEE
       END IF
  380  CONTINUE
  385 FORMAT('======= M-<F===== '/I6,'  ENERGY  ',E15.8,
     *       '  OM ',E15.8/
     *'  DIRECT1   ',E15.4,'  EXCHAN1   ',E15.4/
     *'  DIRECT2   ',E15.4,'  EXCHAN2   ',E15.4/
     *'  AMPLITUDES    '/
     *' RDD ',E11.4,' RDE ',E11.4,' RED ',E11.4,' REE ',E11.4/
     *' IDD ',E11.4,' IDE ',E11.4,' IED ',E11.4,' IEE ',E11.4)
       PRINT 390,SUM1,SUM2,SUM3,SUM4,
     *           SUMI1,SUMI2,SUMI3,SUMI4
  390 FORMAT('  SUM IN M<F  RE     '/4E15.4/
     *       '  SUM IN M<F  IM     '/4E15.4)

       M=AA(8,I)
       M1=AA(9,I)
       C=M-M1
       D=AA(7,I)
      CALL MAEL5(ALFA,BET,NE5,H,NW1,I2,I1,I4,D,M, LX(2),LX(1),LX(4),
     *L,LV1,R1,PR,ME,ME11,ME12,ME13,ME16,R8,R2,R5,R6,R7)
      CALL MAEL5(ALFA,BET,NE5,H,NW1,I5,I6,I3,D,M, LX(5),LX(6),LX(3),
     *L,LV2,R1,PR,ME,ME14,ME15,ME13,ME16,R8,R2,R5,R6,R7)
       PRINT 395
  395  FORMAT('=====M>F  DISCRETE===========')
      DO 400 J=1, M1
       CALL COLUMN(NE5,R8(1,D+J-1),R5)
       E=R5(1)
       OM=(OM1-E)
       IF(EXPI.EQ.1.AND.((DLAM1(D+J-1).NE.0.).OR. 
     *   GAM1(D+J-1).NE.0.)) THEN
       DELE=DLAM1(D+J-1)-OM
       ELSE
       DELE=0
       END IF
       DELE=DELE*AA(12,I)
       OM=OM+DELE
       PRINT 401,DELE
 401   FORMAT(' delE = ',E15.8)
       RDD=ME11(J)*ME14(J)*Y(1,I)
       RDE=ME11(J)*ME15(J)*Y(1,I)
       RED=ME12(J)*ME14(J)*Y(1,I)
       REE=ME12(J)*ME15(J)*Y(1,I)
       IDD=-GAM1(D+J-1)/2*RDD/(OM**2+GAM1(D+J-1)**2/4)
       IDE=-GAM1(D+J-1)/2*RDE/(OM**2+GAM1(D+J-1)**2/4)
       IED=-GAM1(D+J-1)/2*RED/(OM**2+GAM1(D+J-1)**2/4)
       IEE=-GAM1(D+J-1)/2*REE/(OM**2+GAM1(D+J-1)**2/4)
       RDD=RDD*OM/(OM**2+GAM1(D+J-1)**2/4)
       RDE=RDE*OM/(OM**2+GAM1(D+J-1)**2/4)
       RED=RED*OM/(OM**2+GAM1(D+J-1)**2/4)
       REE=REE*OM/(OM**2+GAM1(D+J-1)**2/4)
       PRINT 410,J,E,OM,ME11(J),ME12(J),ME14(J),
     * ME15(J),RDD,RDE,RED,REE,IDD,IDE,IED,IEE
  410 FORMAT(I4,' ENERGY ',F10.4,' DEL E ',E15.4/
     *          ' MATR.EL. ',4E11.4/
     *          ' RE       ',4E11.4/
     *          ' IM       ',4E11.4)
       SUM1=SUM1+RDD
       SUM2=SUM2+RDE
       SUM3=SUM3+RED
       SUM4=SUM4+REE
       SUMI1=SUMI1+IDD
       SUMI2=SUMI2+IDE
       SUMI3=SUMI3+IED
       SUMI4=SUMI4+IEE
  400  CONTINUE
       PRINT 420,SUM1,SUM2,SUM3,SUM4,
     *           SUMI1,SUMI2,SUMI3,SUMI4
  420 FORMAT('  SUM IN DISCRETE M  RE '/4E15.4/
     *       '  SUM IN DISCRETE M  IM '/4E15.4)
       PRINT 430
  430  FORMAT('  M>F CONTINUUM  ',
     * ' COULOMB MATRIX ELEMENTS  '/'    ENERGY  ',
     * '   DIRECT1   ',' EXCHANGE1   ',
     * '  DIRECT2   ',' EXCHANGE2 ')
      DO 440 J=1, C
        ME16(J)=ME11(J+M1)
       ME17(J)=ME12(J+M1)
       ME18(J)=ME14(J+M1)
       ME19(J)=ME15(J+M1)
       CALL COLUMN(NE5,R8(1,D+J-1+M1),R5)
       LAMM(J)=R5(1)
C      PRINT 445,R5(1)
C 445  FORMAT(E15.8)
       PRINT 450,LAMM(J),ME16(J),ME17(J),
     * ME18(J),ME19(J)
  440  CONTINUE
  450  FORMAT(5E12.4)
      DO 500 J=1, C
       B=J+M1
       ME16(J)=ME11(B)*ME14(B)*Y(1,I)
       ME17(J)=ME11(B)*ME15(B)*Y(1,I)
       ME18(J)=ME12(B)*ME14(B)*Y(1,I)
       ME19(J)=ME12(B)*ME15(B)*Y(1,I)
  500  CONTINUE
       CALL ONE(OM1,LAMM,ME16,C,0,RDD,IDD,CC,QQ,DD1,RR)
       CALL ONE(OM1,LAMM,ME17,C,0,RDE,IDE,CC,QQ,DD1,RR)
       CALL ONE(OM1,LAMM,ME18,C,0,RED,IED,CC,QQ,DD1,RR)
       CALL ONE(OM1,LAMM,ME19,C,0,REE,IEE,CC,QQ,DD1,RR)
       PRINT 510
  510  FORMAT('  INTEGRAL IN M>F CONTINUUM'/
     *'       AMPL DD         AMPL DE       AMPL ED        AMPL EE')
      PRINT 520,RDD,RDE,RED,REE
  520  FORMAT('  RE',4E15.4)
      PRINT 530,IDD,IDE,IED,IEE
  530 FORMAT('  IM',4E15.4)
      SUM1=SUM1+RDD
      SUM2=SUM2+RDE
      SUM3=SUM3+RED
      SUM4=SUM4+REE
      SUMI1=SUMI1+IDD
      SUMI2=SUMI2+IDE
      SUMI3=SUMI3+IED
      SUMI4=SUMI4+IEE
      PRINT 540
 540  FORMAT('  RADIAL AMPLITUDES   ')
      PRINT 550,SUM1,SUM2,SUM3,SUM4
  550 FORMAT('  RE   ',4E15.4)
      PRINT 560,SUMI1,SUMI2,SUMI3,SUMI4
  560 FORMAT('  IM   ',4E15.4)
      OM1=Y(2,I)*SUM1+Y(3,I)*SUM2+
     *Y(4,I)*SUM3+Y(5,I)*SUM4
      OM2=Y(2,I)*SUMI1+Y(3,I)*SUMI2+
     *Y(4,I)*SUMI3+Y(5,I)*SUMI4
      PRINT 570,OM1,OM2
  570 FORMAT('  COMPLETE AMPLITUDE'/2E15.4)
      RES1=RES1+OM1
      RES2=RES2+OM2
  300 CONTINUE
       E=2*3.1416*(RES1**2+RES2**2)*13.605
      PRINT 580,RES1,RES2
  580 FORMAT('============================'/
     *       '  AMPLITUDE OF PROCESS '/2E15.4)
      PRINT 590,E
  590 FORMAT('  DOUBLE AUGER PROBABILITY(IN EV)'/E15.4)
 1000 CONTINUE
      RETURN
      END
