C===============================================================
C     PROGRAM FOR CALCULATION OF MESON CAPTURE
C     IN THE HARTREE-FOCK "FROZEN" CORE
C    authors: L.V.Chernysheva, N.A. Cherepkov
C--- LI-initial, LM-final, LL -hole, LM -ejected electron
C===============================================================
      SUBROUTINE MESON(Z,IS,R,NE,BET,H,IN,IL,IQ,L1,SS1,
     *     NU1,IG,GAM,KK,ISIGM1,AKAP,K1,EPS,
     *     MU,IT,NU2,GAM1,GAM2,UD,KSU1,
     *     ISIGM2,IP,S1,TETA,ALM3,ALA,C3,ALA1,ALA2,DIR,
     *     NNL,P1,SM,LAMM,
     *     R1,R2,R3,R4,R5,R6,R9,R10,R8)

      INTEGER IS,IG,SS1,LA,AS,AS1,LM,
     * IP( 1),ISIGM1( 1),ISIGM2( 1),KK( 1),
     * S1(IS),UD(IS),IN(IS),IL(IS),IQ(IS),DIR(IS),
     *  M(3),M1(3),MZ(3),NNL(K1)

      REAL*8 Z,R,BET,H,EPS,AKAP,MU,MU1,
     *   ALFA,RO1,RO2,C,D,EPS1,
     *   A1,LAMI,LAML,LAMK,RO1N,LAMM(3,K1),
     *   DAT(15),GAM( 1),GAM1( 1),
     *   TETA(IS),ALM3(IS),ALA(IS),GAM2(IS),
     *   C3(IS),ALA1(IS),ALA2(IS),P1(K1),SM(K1),
     * R1(NE+15),R2(NE+15),R3(NE+15),R4(NE+15),
     * R5(NE+15),R6(NE+15),R9(NE+15),R10(NE+15),
     *   R8(NE+15,*)

      NE5=NE+15
      EPS1=1.D-05
C     IF(ORT.EQ.0) UD(IS)=0
      IF(L1.EQ.111) GOTO 40
      DO 20 I=1,IS
20    IP(I)=1
      CH=1
      I1=1
      S1(1)=SS1
      DO 30 I=1,IS
      JJA=2*IL(I)+1
      IF(S1(1).LE.99) JJA=JJA*2
      DO 30 J=I1,IS
      JJ=2*IL(J)+1
      IF(S1(1) .LE.99)  JJ=JJ*2
      IF(I.LT.J ) CH=CH+     ( IL(I)+IL(J)-
     *IABS(IL(I)-IL(J)))/2 +1
      IF(I.EQ.J.AND.IQ(I).NE.1) CH=CH+IL(I)
      IF(KSU1.NE.1.AND.(IQ(I).EQ.1.
     * AND.IQ(J).EQ.JJ-1).OR.
     * (IQ(J).EQ.1.AND.IQ(I).EQ.JJA-1))
     * CH=CH+ MIN0(IL(I),IL(J))
 30    CONTINUE
       CALL COEF(GAM,DIR,IS,KSU1,IL,IQ,
     * L1,S1,CH,KK,ISIGM1,ISIGM2,IG,1 )
      NU1=CH
      NU=CH
      GOTO 60
 40   II=IS
      DO 50 I=1,NU1
 50   ISIGM2(I)=II
      ISIGM2(NU1+NU2+1)=IT
 60   WRITE(6,62) UD(IS),IS,NU1,IG
      WRITE(*,62) UD(IS),IS,NU1,IG
 62   FORMAT(//'          MESON CAPTURE  '/
     *         ' ORTOGONALIZATION SIGN',13X,I2/
     * 21H NUMBER OF N-L SHELLS,11X,
     *I5/20H NO OF COEFFICIENTS ,12X, I5,
     * /' NO OF EXCH. COEFFICIENTS       ',I5)
      REWIND 9
      READ(9)(DAT(JJ1),JJ1=1,15)
      IF(DAT(1).EQ.1.) GOTO 70
      WRITE(6,80)DAT(1)
 80   FORMAT(5X,' KSU#1,   KSU=  ',E8.0)
      RETURN
  70  IF(DABS(DAT(7)-Z).LT.EPS1.AND.
     *DABS(DAT(2)-R).LT.EPS1.AND.
     *DABS(DFLOAT(IS-2) -DAT(8)).LT.EPS1)  GOTO 110
      WRITE(6,100)(DAT(JJ1),JJ1=1,15),
     *R, Z,IS,EPS1
 100  FORMAT(5X,' CHECK "CORE" AND "EXCITED" PARAMETERS',
     */,5X,' PL= ',E15.8,3X,' R=  ',E15.8/
     *  3X,' N=  ',E15.8,3X,' H=  ',E15.8,3X/
     *     ' BET=',E15.8,3X,' EPS=',E15.8,3X/
     *     ' Z=  ',E15.8,3X,2(4(3X,E15.8),/),
     *                   3X,2E12.5,I3,E15.8)
       RETURN
 110  CONTINUE

         UD(IS+1)=1
        DO 194 I=1,IS-1
 194    IF(UD(I).EQ.1) GOTO 195
        GOTO 196
 195    IF(UD(IS).EQ.0) THEN
        UD(IS)=1
        UD(IS+1)=0
        END IF
 196  CALL VAR11(NE5,Z,R,BET,H,ALFA,RO1N,R1)
      CALL COMON3( Z,R,BET,H,GAM,GAM1,GAM2,EPS,ALFA,
     * NE,IN,IL,IQ,IG,NU,NU1,NU2,
     * IT,UD,DIR,IP,KK,ISIGM1,ISIGM2,IS,L1,SS1,MU)
      IF(UD(IS).EQ.0) THEN
      ISIGM2(NU+2)=0
      ELSE
      ISIGM2(NU+2)=1
      END IF
      WRITE(6,140)R1(3)
 140  FORMAT(' R MIN',22X,E15.8  )
      IRA=IS-2
      WRITE(6,190)
      WRITE(*,190)
 190  FORMAT('==========================='/
     *       '  CORE WAVE FUNCTIONS:  '/
     *       ' E                 N    L    Q  ZONE')
      DO 260 IW=1,IRA
      CALL  DISKO1(9,IW,2,NE5,R8(1,IW))
      MMOM=IFIX(SNGL(R8(2,IW)))
      NOBL=IFIX(SNGL(R8(NE+5,IW)))
      R8(NE+6,IW)=1.
      KOLE=IFIX(SNGL(R8(NE+6,IW)))
      R8(2,IW)=1.D 00
      R8(NE+3,IW)=0.D 00
      DO 270 J=1,NE
 270  R8(J+2,IW)=R8(J+2,IW)*DSQRT(R8(2,IW))
      WRITE(6,220) R8(1,IW),NOBL,MMOM,KOLE,IW
      WRITE(*,220) R8(1,IW),NOBL,MMOM,KOLE,IW
 220  FORMAT( E15.8,4I5)
 260  CONTINUE
      REWIND 9
      NH=NE
 250  ISIGM=IS
 290  I=IS
 350  CONTINUE
C      IF(K2.LT.1) GOTO 360
      RO1=-BET*(10.D 00 + DLOG(Z))
      RO2=RO1+DFLOAT(NE)*H
      ALFA=(RO2- BET*DLOG(R))/R
C ========GET WAVEFUNCTIONS OF MESON STATES==========
         NN1=0
         AS=IS+2
      DO 61 I=1,2
         LNO=11+I-1
         REWIND LNO
         READ (LNO) (DAT(J),J=1,15)
C        DAT(13)=11
         M(I)=DAT(13)+DAT(12)
         M1(I)=DAT(12)
         MZ(I)=AS+1
       WRITE(6,71)M(I),M1(I),MZ(I)
       WRITE(*,71)
  71  FORMAT(3I3, 60('=')/' MESON  STATES  '/
     *       '      ENERGY      PHASE OR N',
     *       '        IN ZONE')
      DO 81 J=1,M(I)
      LA=AS+J
      CALL DISKO1(LNO,J,2,NE5,R8(1,LA))
      R8(NE+2,LA)=0.0D 00
      IA=7
      IF(R8(1,LA).LT.0.D 00) IA=5
      WRITE(6,91)R8(1,LA),R8(NE+IA,LA),LA
      WRITE(*,91)R8(1,LA),R8(NE+IA,LA),LA
  91  FORMAT(2X,E15.8,1X,E15.8,I6)
  81  CONTINUE
      AS=LA
      WRITE(6,101) R8(2,LA)
  101 FORMAT(' WITH ELECTRON ORBITAL MOMENTUM ', F4.1)
  61  CONTINUE
C+++++++++++++++++++++++++++
C======to calculate auger-electron wave-function=====
      IF(K1.LT.1) GOTO 360
      DO 450 I1=1,K1
      SM(I1)=0.
      CALL COLUMN(NE5,R8(1,IT),R4)
      LL=IL(IT)
      LAML=R8(1,IT)
      WRITE(6, 451)LL,R8(1,IT),IT
451   FORMAT('****************************************'/
     *       ' HOLE '/
     *       ' L ',I4,' E= ',E15.8,' ZONE',I4)
C=================
      CALL COLUMN(NE5,R8(1,MZ(1)),R2)
      LAMI=R8(1,MZ(1))
      LI=IFIX(SNGL(R8(2,MZ(1))))
      WRITE(6,453)LI,LAMI,MZ(1)
453   FORMAT(' MESON INITIAL STATE '/
     *       ' L ',I4,' E= ',E15.8,' ZONE',I4)

C=================================================
      AS1=IS+2
      III=I1+MZ(2)-1
C     WRITE(6,600)III
C 600  FORMAT('III=  ' I3)
      CALL COLUMN(NE5,R8(1,III),R3)
      CALL COLUMN(NE5,R3,R8(1,IS-1))
      LM=IL(III)
      LAMM(1,I1)=R8(1,III)
      NNL(I1)=IFIX(SNGL(R8(NE+5,III)))
      WRITE(6,452)I1,LM,NNL(IT+1),R3(1),III
452   FORMAT(' MESON FINAL STATE '/
     *  I3,' L ',I3,' N ',I3,' E= ',E15.8,' ZONE',I4)
      R8(1,3)=0.
      WRITE(6,601)I1,R8(1,1),R8(1,2),R8(1,3) 
 601  FORMAT(I4,'E1,E2,E3-field',3E15.8)
C====================
      LAMK=LAMI+LAML-LAMM(1,I1)
      LAMM(3,I1)=LAMK
      LK=IL(IS)
      WRITE(6, 456)LAMI,LAMM(1,I1),LAML,LAMK
 456  FORMAT(' ENERGIES imlk ',4E15.8)
      IF (LAMK.GT.0.) THEN
      AKAP=DSQRT(LAMK)
           WRITE(6,461)LI,LM,LL,LK
461   FORMAT(' LI,LM,LL,LK ',4I4)
C===========================================

      WRITE(6,460) I1
      WRITE(*,460) I1
 460  FORMAT('  CONTINUOUS SPECTRUM WAVE ' 
     *       'FUNCTION  OF EJECTED ELECTRON    '/
     *      ' IN THE FIELD OF CORE AND FINAL MESON',I5)
      CALL ZERO( R5,Z,R8 ,TETA,IS,IN,IQ,NE5,KSU1)
      ISIGM=IS
       MU1=1.
       DO 470 IW=1,IS
  470  ALM3(IW)=0.0
C       PRINT 471,R8(2,1),R8(2,2),R8(2,3)
C 471  FORMAT(' NORMA ',3E15.8)
      CALL MAIN23(ALFA,BET,H,Z,AKAP,EPS,TETA,
     *R1,R2,R3,R4,R5,R6,R9,R10,R8,ALA,ALA1,ALA2,ALM3,C3,
     *GAM,GAM1,GAM2,MU1,IS,NE5,IL,KK,ISIGM1,ISIGM2,
     *IN,IQ,NU,IG,ISIGM,DIR,UD,1,1)
C=================
      CALL COLUMN(NE5,R8(1,IT),R4)
      CALL COLUMN(NE5,R8(1,III),R3)
      CALL COLUMN(NE5,R8(1,MZ(1)),R2)
      WRITE(6,472) R5(1)
 472  FORMAT(' ENERGY',27X, E15.8)

      L3=MAX0(ABS(LI-LM),ABS(LL-LK))
      L2=MIN0(LI+LM,LL+LK)
C      WRITE(6,483)LI,LM,LL,LK,L3,L2
C483   FORMAT(' LI,LM,LL,LK, L3,L2 ',6I4)
      DO 454 IW=1,NE
      R6(IW+1)=R3(IW+2)*R2(IW+2)
454   CONTINUE
        DO 500 I=L3, L2, 2
        A1=1./DFLOAT(2*I+1)
       CALL POT(ALFA,BET,H,R1,R6,I,NE5)
       CALL SIMPS2(ALFA,BET,H,P1(I1),R1,R4,R5,R6,R6,NE5,0)
       CALL IOT31(LI,I,LM,0,0,0,C)
       CALL IOT31(LL,I,LK,0,0,0,D)
       WRITE(6,501) P1(I1),C,D
 501   FORMAT(' p1,c,d',3E15.8)
      WRITE(6,483)LI,LM,LL,LK,L3,L2
483   FORMAT(' LI,LM,LL,LK, L3,L2 ',6I4)
       P1(I1)=P1(I1)*2.*C*D*DSQRT(DFLOAT((2*LI+1)*
     * (2*LK+1)*(2*LM+1)*(2*LL+1)))
       SM(I1)=SM(I1)+A1*P1(I1)**2
 500   CONTINUE
      IF ( LAMI.GT.0. ) THEN
       SM(I1)=SM(I1)*69.461/LAMI/MU
      ELSE
       SM(I1)=SM(I1)*1.297
      END IF
      END IF
450    CONTINUE
C       PRINT 520,II
C520    FORMAT('  MATRIX ELEMENT ',I4)
C---------------------------------------------
      WRITE(6,565)LAML,LL,LAMI,LI,MU
565   FORMAT(// '==================================='/
     *       '  E core           ',E15.8,' L ',I4/
     *       '  E initial meson  ',E15.8,' L ',I4,/
     *       '  MESON MASS       ',10X,F5.1)

      IF ( LAMI.GT.0 ) THEN
       WRITE(6, 530)
       WRITE(*, 530)
530    FORMAT (10X/'=========CROSS SECTION=========== '/
     *        '  FINAL STATE OF MESON  '/
     *        '      N   L  ',
     *        '  ERY(MESON)    E1Ry(ELECTR) ',
     *        '  MATRIX ELEMENT ',
     *        '  CROSS SECTION ')
      ELSE
       WRITE(6, 540)
       WRITE(*, 540)
 540   FORMAT('=======AUGER TRANSITI0N PROBABILITY '/
     *        '  FINAL STATE OF MESON  '/
     *        '      N    L  ',
     *        '  ERy(MESON)    E1Ry(ELECTR) ',
     *        '  MATRIX ELEMENT ',
     *        '  PROBABIL.  ')
      END IF
      DO 550 IW=1, K1
      WRITE(6,560)NNL(IW),LM,LAMM(1,IW),LAMM(3,IW),P1(IW),SM(IW)
560   FORMAT(3X,2I4,4E15.8)
550   CONTINUE
360   RETURN
      END

