C     PROGRAM FOR CALCULATION OF SINGLE PHOTON DECAY
C         OF DOUBLE  VACANCY STATE
C     AUTHORS L.V.CHERNYSHEVA, M.YA.AMUSIA, V.K.KILIN,
C             A.N.KOLESNIKOVA, I.S.LI
C     fortran version: L.V.Chernysheva
C     new version  July.1992 with V.A Kilin
C                 VARIABLES
C     ND -number of diagrams
C     NES - number of excited ststes
C     AA, Y - matrices
C     ORI.ORF- state momenta
C     SP spin
C    NORM = NORM init++(-2)*NORM fin**(-2)
C     EXPI=1(input KK,IND,DLAM,GAM- arrays) or EXPI=0
C     KK,IND,DLAM,GAM -in dauger.f
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 PHODEC(ND,NES,AA,ORI,ORF,SP,NORM,EK,Y,
     *    EXPI,KK,IND,DLAM,GAM,NE,S,NW1,LI,ME,
     *    ME11,ME12,ME13,ME14,ME15,ME16,ME17,ME18,ME19,
     *    LAMI,LAMM,
     *    CC,QQ,DD1,RR,DLAM1,GAM1,
     *    R1,R2,R5,R6,R7,R8)

      INTEGER NE,S,ND,LV,KSU,M,M1,I1,I2,I3,I4,
     *    A,B,C,D,ORI,ORF,SP,I,J,K,L,NW1,EXPI,KK,
     *    AA(11,ND),LI(S),LX(6),IND(KK)

      REAL*8 Z,R,BET,ALFA,H,E,NORM,B1,B2,RO1,DELE,
     *    SUM1,SUM2,SUM3,SUM4,RES1,RES2,
     *    RDIDL,RDIDV,REXDL,REXDV,IDIDL,IDIDV,IEXDL,IEXDV,
     *    OM,OM1,OM2,OM3,OM4,OME,OME1,RES3,RES4,EK,
     *    DLAM(KK),DLAM1((NES+1)*NW1),GAM(KK),GAM1((NES+1)*NW1),
     *    SUMI1,SUMI2,SUMI3,SUMI4,
     *    LAMI(S),LAMX(6),ME4(1),ME5(1),ME6(1),ME7(1),
     *    ME8(1),ME9(1),Y(3,ND),DAT(15),ME(1,NW1),
     *    ME11(NW1),ME12(NW1),ME13(NW1),ME14(NW1),
     *    ME15(NW1),ME16(NW1),ME17(NW1),ME18(NW1),
     *    ME19(NW1),LAMM(NW1),PR(20),
     *    CC(NW1),QQ(NW1),DD1(NW1),RR(NW1),
     *    R1(NE+15),R2(NE+15),R5(NE+15),
     *    R6(NE+15),R7(NE+15),R8(NE+15,S+NW1*ND)

      PRINT 10
 10   FORMAT('====AMPLITUDES OF THE  DOUBLE HOLE=== '/
     *       '======STATES  RADIATIONAL DECAY===== ')
      READ(11) (DAT(I),I=1,15)
      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,NE,R,H,BET,ALFA,EK,ORI,ORF,SP,NORM
  20  FORMAT('  CHARGE OF NUCLEUS     ',E15.8/
     *       '  NUMBER OF SHELLS      ',I10/
     *       '  NO OF EXCITED STATES  ',I10/
     *       '  NUMBER OF POINTS      ',I10/
     *       '  R MAX                 ',E15.8/
     *       '  STEP OF INTEGRATION   ',E15.8/
     *       '  BETTA                 ',E15.8/
     *       '  ALFA                  ',E15.8/
     *       '  PHOTON ENERG“Y         ',E15.8/
     *       '  INITIAL STATE MOMENTUM',I10/
     *       '  FINAL STATE MOMENTUM  ',I10/
     *       '  SPIN                  ',I10/
     *       '  NORMA                 ',E15.8)
       PRINT 30
   30  FORMAT('  INITIAL MATRICES')
       PRINT 40
   40  FORMAT('  WAVE FUNCTIONS IN ZONES    M   M1',
     *  '  LV1  KSU    N   KSU1')
       WRITE(6,41) ((AA(J,I),J=1,11),I=1,ND)
   41  FORMAT(11I5)
       IF ( AA(10,1).EQ.0 ) THEN
       PRINT 50
   50  FORMAT('  ANGLE     SPIN DI     SPIN EX')
       WRITE(6,51) ((Y(J,I),J=1,3),I=1,ND)
       END IF
   51  FORMAT(3F11.3)
      PRINT 60
  60  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 82,I,LI(I),LAMI(I),I, DLAM1(I),GAM1(I)
       END IF
  80   CONTINUE
  90   FORMAT(2I6,E15.8,I6)
  82   FORMAT(2I6,E15.8,I6,2E15.8)
C ========GET WAVEFUNCTIONS OF EXCITED STATES=======

  120  A=S
       B=D-1
      DO 125 I=1, NES
C      M=M3(I)
      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 130,I
 130  FORMAT('============================'/ 
     *       '    EXCITED STATES   CHANNEL ',I6)
      PRINT 135
 135  FORMAT('  ENERGIES          PHASE or n',
     *       '        IN ZONES',
     *       '  EXP.ENERGY    GAMMA ')
      DO 140 K=1, M
       LA=A+K
       DLAM1(LA)=0
       GAM1(LA)=0
        IF(EXPI.EQ.1) THEN
          DO 141 II=1,KK
          IF (LA.EQ.IND(II)) THEN
             DLAM1(LA)=DLAM(II)
             GAM1(LA)=GAM(II)
          END IF
 141      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,150) K,R8(1,LA),R8(NE+IA,LA),LA
       ELSE
       PRINT 142,K,R8(1,LA),R8(NE+IA,LA),LA,DLAM1(LA),GAM1(LA)
       END IF
  140  CONTINUE
       A=LA
  150  FORMAT(I4,2E15.8,I6)
  142  FORMAT(I4,2E15.8,I6,2E15.8)
       PRINT 155, R8(2,LA)
  155  FORMAT(' WITH ORBITAL  MOMENTUM=',F10.4)
  125  CONTINUE

C=======CALCULATIONS OF DIAGRAM CONTRIBUTION========
  300  RES1=0
       RES2=0
       RES3=0
       RES4=0
      DO 310 I=1, ND
       I1=AA(1,I)
       I2=AA(2,I)
       I3=AA(3,I)
       I4=AA(4,I)
       M=AA(6,I)
       M1=AA(7,I)
       PRINT 315,I
 315  FORMAT('====================================='/
     *   ' DIAGRAM NO',I6/'  ORBITALS   L   ENERGIES')
      DO 317 J=1, 5
        A=AA(J,I)
       CALL COLUMN(NE5,R8(1,A),R5)
       LX(J)=IFIX(SNGL(R5(2)))
       LAMX(J)=R5(1)
       PRINT 320,J,LX(J),LAMX(J),A
  317  CONTINUE
  320  FORMAT(2I6,F10.4,I4)
       LV=AA(8,I)
       KSU=AA(9,I)
       L=LX(5)
       PRINT 330,LV
 330   FORMAT('  TRANSFERED MOMENTA ',I5)
      IF ( AA(10,I).EQ.1 ) THEN
       CALL IOT61(ORI,1,ORF,LX(3),LX(1),L,B1)
       CALL IOT61(ORI,L,LX(1),LV,LX(2),LX(4),B2)
       Y(1,I)=B1*B2*(-1)**(ORF+LX(2))
       Y(2,I)=1
       Y(3,I)=(-1)**SP
       END IF
      IF ( AA(10,I).EQ.2 ) THEN
       CALL IOT61(ORI,1,ORF,LX(3),LX(1),L,B1)
       CALL IOT61(ORI,L,LX(1),LV,LX(2),LX(4),B2)
       Y(1,I)=B1*B2*(-1)**(ORI+LX(3)+LX(4)+LX(1))
       Y(2,I)=1
       Y(3,I)=(-1)**SP
       END IF
      IF ( AA(10,I).EQ.3 ) THEN
      CALL IOT61(ORI,1,ORF,L,LX(1),LX(3),B1)
      CALL IOT61(L,LX(1),ORF,LX(2),LX(4),LV,B2)
       Y(1,I)=B1*B2*(-1)**(LX(2)+L+LX(1)+1)
       Y(2,I)=(-1)**SP
       Y(3,I)=1
       END IF
      IF ( AA(10,I).EQ.4 ) THEN
      CALL IOT61(L,LX(1),ORF,LX(2),LX(4),LV,B1)
      CALL IOT61(ORI,1,ORF,L,LX(1),LX(3),B2)
       Y(1,I)=B1*B2*(-1)**(ORI+ORF+LX(4))
       Y(2,I)=(-1)**SP
       Y(3,I)=1
       END IF
      IF ( AA(10,I).GT.4 ) PRINT 340
 340   FORMAT('  MISTAKE')
       DO 355 J=1,3
  355  PRINT 356, J,I,Y(J,I)
  356  FORMAT('   Y(',I4,',',I4,')=',F10.5)
       Y(1,I)=Y(1,I)/DSQRT(NORM)
       SUM1=0
       SUM2=0
       SUM3=0
       SUM4=0
       SUMI1=0
       SUMI2=0
       SUMI3=0
       SUMI4=0
       OM=LAMX(2)+LAMX(4)-LAMX(1)
       DELE=0
      IF ( KSU.EQ.0 ) GOTO 350
       PRINT 360,L
 360  FORMAT('  INTERMEDIATE STATE  ORBITAL MOMENTUM ',I4)
       ME6(1)=LAMI(J)
      DO 370 J=1, S
      IF ( L.EQ.LI(J) ) THEN
       E=LAMI(J)
       OME=OM-E
       PRINT 372,DLAM1(J),GAM1(J)
  372  FORMAT(' DLAM1,GAM1 ',2E15.8)
       IF((EXPI.EQ.1).AND.((DLAM1(J).NE.0.) 
     *  .OR.(GAM1(J).NE.0.))) THEN
       DELE=DLAM1(J)-OME
      ELSE
       DELE=0
       END IF
       DELE=DELE*AA(11,I)
       PRINT 371,DELE,OME
       OME=OME+DELE
 371   FORMAT(' delE ome ',2E15.8)
      CALL MAEL1(ALFA,BET,NE5,H,J,I3,1,LX(3),
     * L,LV,1,ME4,ME5,R1,R2,R5,R6,R8)
      CALL MAEL5(ALFA,BET,NE5,H,NW1,I4,I2,I1,J,1,LX(4),LX(2),LX(1),
     *L,LV,R1,PR,ME,ME7,ME8,ME6,ME9,R8,R2,R5,R6,R7)
       OME1=E-LAMX(3)
       RDIDL=ME4(1)*ME7(1)*Y(1,I)
       RDIDV=ME5(1)*ME7(1)*Y(1,I)*2/OME1
       REXDL=ME4(1)*ME8(1)*Y(1,I)
       REXDV=ME5(1)*ME8(1)*Y(1,I)*2/OME1
       IDIDL=-RDIDL*GAM1(J)/2/(OME**2+GAM1(J)**2/4)
       IDIDV=-RDIDV*GAM1(J)/2/(OME**2+GAM1(J)**2/4)
       IEXDL=-REXDL*GAM1(J)/2/(OME**2+GAM1(J)**2/4)
       IEXDV=-REXDV*GAM1(J)/2/(OME**2+GAM1(J)**2/4)
       RDIDL=RDIDL*OME/(OME**2+GAM1(J)**2/4)
       RDIDV=RDIDV*OME/(OME**2+GAM1(J)**2/4)
       REXDL=REXDL*OME/(OME**2+GAM1(J)**2/4)
       REXDV=REXDV*OME/(OME**2+GAM1(J)**2/4)
      PRINT 380,J,E,OME,ME4,ME5,ME7,ME8,RDIDL,RDIDV,REXDL,REXDV,
     *          IDIDL,IDIDV,IEXDL,IEXDV
       SUM1=SUM1+RDIDL
       SUM2=SUM2+RDIDV
       SUM3=SUM3+REXDL
       SUM4=SUM4+REXDV
       SUMI1=SUMI1+IDIDL
       SUMI2=SUMI2+IDIDV
       SUMI3=SUMI3+IEXDL
       SUMI4=SUMI4+IEXDV
       END IF
  370  CONTINUE
  380 FORMAT('=====  M-<F ',I6,'  E ',E12.5,' OM ',E11.4/
     *' DIP LF ',E10.4,' DIP VF  ',E10.4/
     *' DIRECT ',E10.4,' EXCHAN  ',E10.4/
     *'  AMPLITUDE    '/
     *' RDIDL ',E10.4,' RDIDV ',E10.4,' REXDL ',E10.4,' REXDV ',E10.4/
     *' IDIDL ',E10.4,' IDIDV ',E10.4,' IEXDL ',E10.4,' IEXDV ',E10.4)
       PRINT 390,SUM1,SUM2,SUM3,SUM4,
     *           SUMI1,SUMI2,SUMI3,SUMI4
 390  FORMAT ('  SUM IN M-<F'/4E15.6/4E15.6)
 350   C=M-M1
       D=AA(5,I)
      CALL MAEL5(ALFA,BET,NE5,H,NW1,I4,I2,I1,D,M, LX(4),LX(2),LX(1),
     * L,LV,R1,PR,ME,ME11,ME12,ME13,ME14,R8,R2,R5,R6,R7)
      CALL MAEL1(ALFA,BET,NE5,H,D,I3,M,LX(3),
     *L,LV,1,ME14,ME15,R1,R2,R5,R6,R8)
      PRINT 391
 391  FORMAT('======  M>F  DISCRETE===========')
      DO 395 J=1 , M1
      CALL COLUMN(NE5,R8(1,D+J-1),R5)
      E=R5(1)
      OME=(OM-E)
       PRINT 397,DLAM1(D+J-1),GAM1(D+J-1)
 397   FORMAT(' DLAM1, GAM1 ',2E15.8)
      IF((EXPI.EQ.1).AND.((DLAM1(D+J-1).NE.0.).OR.
     *  (GAM1(D+J-1).NE.0))) THEN  
       DELE=DLAM1(D+J-1)-OME
       ELSE
       DELE=0
       END IF   
       DELE=DELE*AA(11,I)
       PRINT 396,DELE,OME
       OME=OME+DELE
 396   FORMAT(' del E ome  ',2E15.8)
       OME1=E-LAMX(3)
       RDIDL=ME11(J)*ME14(J)*Y(1,I)
       RDIDV=ME11(J)*ME15(J)*Y(1,I)*2/OME1
       REXDL=ME12(J)*ME14(J)*Y(1,I)
       REXDV=ME12(J)*ME15(J)*Y(1,I)*2/OME1
       IDIDL=-RDIDL*GAM1(D+J-1)/2/(OME**2+GAM1(D+J-1)**2/4)
       IDIDV=-RDIDV*GAM1(D+J-1)/2/(OME**2+GAM1(D+J-1)**2/4)
       IEXDL=-REXDL*GAM1(D+J-1)/2/(OME**2+GAM1(D+J-1)**2/4)
       IEXDV=-REXDV*GAM1(D+J-1)/2/(OME**2+GAM1(D+J-1)**2/4)
       RDIDL=RDIDL*OME/(OME**2+GAM1(D+J-1)**2/4)
       RDIDV=RDIDV*OME/(OME**2+GAM1(D+J-1)**2/4)
       REXDL=REXDL*OME/(OME**2+GAM1(D+J-1)**2/4)
       REXDV=REXDV*OME/(OME**2+GAM1(D+J-1)**2/4)
      PRINT 400,J,E,OME,ME14(J),ME15(J),ME11(J),ME12(J),
     *RDIDL,RDIDV,REXDL,REXDV,IDIDL,IDIDV,IEXDL,IEXDV
  400 FORMAT('  M>F ',I6,'  E=',E15.8,' OM ',E15.8/
     *' DIP LF  ',E10.4,' DIP VF  ',E10.4/
     *' DIRECT  ',E10.4,' EXCHAN  ',E10.4/
     *'  AMPLITUDE     '/
     *' RDIDL ',E10.4,' RDIDV ',E10.4,' REXDL ',E10.4,' REXDV ',E10.4/
     *' IDIDL ',E10.4,' IDIDV ',E10.4,' IEXDL ',E10.4,' IEXDV ',E10.4)
      SUM1=SUM1+RDIDL
      SUM2=SUM2+RDIDV
      SUM3=SUM3+REXDL
      SUM4=SUM4+REXDV
      SUMI1=SUMI1+IDIDL
      SUMI2=SUMI2+IDIDV
      SUMI3=SUMI3+IEXDL
      SUMI4=SUMI4+IEXDV
 395  CONTINUE
      PRINT 410,SUM1,SUM2,SUM3,SUM4,
     *          SUMI1,SUMI2,SUMI3,SUMI4
 410  FORMAT('  SUM IN DISCRETE  M  '/4E15.6)
      PRINT 420
 420  FORMAT('  M>F CONTINUUM')
      PRINT 450
 450  FORMAT('  COULOMB MATRIX ELEMENTS')
      DO 430 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)
      PRINT 460,R5(1),ME16(J),ME17(J),ME18(J),ME19(J)
 430  CONTINUE
 460  FORMAT('  E ',E10.4,
     *       ' DIR ',E10.4,
     *       ' EXC ',E10.4,
     *       ' DIPLF ',E10.4,
     *       ' DIPVF ',E10.4)
      DO 470 J=1,C
        B=J+M1
      ME16(J)=ME11(B)*ME14(B)*Y(1,I)
      ME17(J)=ME11(B)*ME15(B)*2/(LAMM(B)-LAMX(3))*Y(1,I)
      ME18(J)=ME12(B)*ME14(B)*Y(1,I)
      ME19(J)=ME12(B)*ME15(B)*2/(LAMM(B)-LAMX(3))*Y(1,I)
 470   CONTINUE
      CALL ONE(OM,LAMM,ME16,C,0,RDIDL,IDIDL,CC,QQ,DD1,RR)
      CALL ONE(OM,LAMM,ME17,C,0,RDIDV,IDIDV,CC,QQ,DD1,RR)
      CALL ONE(OM,LAMM,ME18,C,0,REXDL,IEXDL,CC,QQ,DD1,RR)
      CALL ONE(OM,LAMM,ME19,C,0,REXDV,IEXDV,CC,QQ,DD1,RR)
      PRINT 480,RDIDL,RDIDV,REXDL,REXDV,IDIDL,IDIDV,IEXDL,IEXDV
480   FORMAT('  INTEGRAL IN M>F CONTINUUM'/
     * '  RE   AMPL DIDL  ',E10.4,' AMPL DIDV  ',E10.4/
     * '     AMPL EXDL    ',E10.4,' AMPL EXDV  ',E10.4/
     * '  IM   AMPL DIDL  ',E10.4,' AMPL DIDV  ',E10.4/
     * '     AMPL EXDL    ',E10.4,' AMPL EXDV  ',E10.4)
        SUM1=SUM1+RDIDL
        SUM2=SUM2+RDIDV
        SUM3=SUM3+REXDL
        SUM4=SUM4+REXDV
        SUMI1=SUMI1+IDIDL
        SUMI2=SUMI2+IDIDV
        SUMI3=SUMI3+IEXDL
        SUMI4=SUMI4+IEXDV
        PRINT 490,SUM1,SUM2,SUM3,SUM4,
     *            SUMI1,SUMI2,SUMI3,SUMI4
490     FORMAT('  RADIAL AMPLITUDE'/
     * '  RE ',4E15.6/'  IM ',4E15.6)
        OM1=(Y(2,I)*SUM1+Y(3,I)*SUM3)
        OM3=(Y(2,I)*SUMI1+Y(3,I)*SUMI3)
        OM2=(Y(2,I)*SUM2+Y(3,I)*SUM4)
        OM4=(Y(2,I)*SUMI2+Y(3,I)*SUMI4)
        PRINT 500,I,OM1,OM3,I,OM2,OM4
500     FORMAT('  COMPLETE AMPLITUDE LF',I3,E10.4,'+I(',E10.4,')'/
     *         '  COMPLETE AMPLITUDE VF',I3,E10.4,'+I(',E10.4,')')
         RES1=RES1+OM1
         RES2=RES2+OM2
         RES3=RES3+OM3
         RES4=RES4+OM4
  310    CONTINUE
C======AMPLITUDE OF PROCESS=====================
      E=(RES1**2+RES3**2)*13.605*(EK/137)**3
     *  *(2*ORF+1)/3
      H=(RES2**2+RES4**2)*13.605*(EK/137)**3
     *  *(2*ORF+1)/3
      PRINT 510,RES1,RES3,RES2,RES4,E,H
 510  FORMAT('=================================='/
     *       '    AMPLITUDE OF PROCESS LF         ',E10.4,
     *       ' + I*',E10.4/
     *       '    AMPLITUDE OF PROCESS VF         ',E10.4,
     *       ' + I*',E10.4/
     *       '    DOUBLE PROBABILITY   LF (IN EV) ',E15.8/
     *       '    DOUBLE PROBABILITY   VF (IN EV) ',E15.8)
       RETURN
      END
