C     PROGRAM FOR CALCULATION OF THE ELECTRON
C     SHIFTS WITH POLARIZATION CORRECTIONS
C      f before m integration
C      L.V.CHERNYSHEVA,M.YA.AMUSIA, D.DAVIDOVIC,
C      N.A.CHEREPKOV
C                VARIABLES
C      ND- number of diagrams
C      NI  number of i-holes
C      NES- NUMBER OF EXCITED STATES
C      K,KZ-NUMBER OF FUNCTIONS AND THE BEGINNING ON DISC
C      F,FZ- JUST THE SAME
C      M,MZ- JUST THE SAME
C      F1-NUMBER OF DISCRETE FUNCTIONS AMONG F(M)
C      IT,IL- THE BEGINNING OF VACANCIES
C      LV-TRANSFERED ORBITAL MOMENTUM
C      EXPI=0,IPOS=0 -Positr-ENERGY   EXPI=1,IPOS=0.5
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      KK- NO of discrete core functions with l=lk
C====    IB(nd) =0,1( O - no reverse diagram)
C     ============================================
      SUBROUTINE  SCATSP(NES,ND,NI,MU,K,KK,KZ,F,F1,FZ,M,M1,MZ,
     *         IT,IL,LV,EXPI,IPOS,NE,S,NW1,EMIN,EMAX,
     *         M3,M4,LLI,LAMM,LAMF,LAMK,PHASE,
     *         PHRESH,PHIMSH,PHRESC,PHIMSC,LAM,ME,
     *         ME1,ME2,ME3,ME4,ME5,ME6,ME7,RP,RM,
     *         RE,IM,MME1,MME2,MME4,MME5,KK1,RRR,RRM,
     *         SR,SIGM1,SIGM2,SIGM3,SIGM4,SIGM5,
     *         CC,QQ,DD1,R1,R2,R5,R6,R7,R8)
C =================================================
       INTEGER NE,NES,NI,NW1,S,NK,NM,NF,EXPI,
     *        W1,W2,IWK,IW,U1,JN,I,J,AS,NA,
     *        DIRK,DIRM,SPI,UN,NNL,DIRF,DIRI,
     *        M,M1,MZ(ND),K,K1,KK,KZ,F,F1,FZ(ND),
     *        IT(ND),IL(ND),LV(ND),LI,LL,LK,LM,LF,
     *        M3(NES),M4(NES),LLI(S),KK1(S)
      REAL*8 Z,R,BET,ALFA,H,RO1,OM,PI,MU,OM1,DOP,
     *       DELXF,DELXM,IPOS,LAMI,LAML,EMIN,EMAX,
     *       RE1,RE2,RE3,RE4,RE5,IM1,IM2,IM3,IM5,
     *       SUM,SUM1,SUM2,SUM3,Y1,SIMPS1,FA,
     *       A,B,A1,A2,A3,A4,SUMRE3,SUMRE4,
     *       RESI2,IMSI2,RESI3,IMSI3,
     *       SIGM1(K),SIGM2(K),SIGM3(K),
     *       SIGM4(K), SIGM5(K),
     *       PR(20),CC(NW1),QQ(NW1),DD1(NW1),
     *       RP(NW1),RM(NW1),RE(NW1),IM(NW1),
     *       LAM(S),LAMF(F),LAMM(M),LAMK(K+KK),
     *       PHASE(K+KK),PHRESH(K),PHIMSH(K),
     *       PHRESC(K),PHIMSC(K),
     *       ME(M,F),ME1(M,F),ME2(M,F),ME3(M,F),
     *       ME4(NW1),ME5(NW1),ME6(NW1),ME7(NW1),
     *       MME1(M,F,K+KK),MME2(M,F,K+KK),
     *       MME4(M,K+KK),MME5(M,K+KK),
     *       RRR(K+KK,K+KK),RRM(K+KK,K+KK),SR(2*(K+KK),2*(K+KK)),
     *       R1(NE+15),R2(NE+15),R5(NE+15),R6(NE+15),
     *       R7(NE+15),R8(NE+15,*),DAT(15)
      PI=3.141 592
      FA=0.
C---mu=1 for electron, mu=-1 for positron-----
      IF(MU.EQ.1.D 00) THEN
      WRITE(6, 5)
      WRITE(*, 5)
      ELSE
      WRITE(6, 6)
      WRITE(*, 6)
      END IF
  5   FORMAT('         ELECTRON ')
  6   FORMAT('         POSITRON ')
      WRITE(6, 10)
      WRITE(*, 10)
  10  FORMAT('     SCATTERING PHASE SHIFTS'/
     *       '  WITH POLARIZATION CORRECTIONS'/
     *       '  f before m'/)
      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)
      WRITE(6, 20)Z,S,NE,R,H,BET,ALFA
  20  FORMAT('  CHARGE OF NUCLEUS   ',E15.8/
     *       '  NUMBER OF SHELLS    ',I5/
     *       '  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 600
      WRITE(6, 70)
      WRITE(*, 70)
  70  FORMAT('============================='/
     *       '   GROUND STATE CONFIGURATION'/
     * '    NO     N     L     Q   SPI    ENERGY     IN ZONES')
       DO 80 I=1, S
      CALL DISKO1( 9,I,2,NE5,R8(1,I))
       LAM(I)=R8(1,I)
        LLI(I)=IDINT(R8(2,I))
       UN=IDINT(R8(NE+5,I))
       SPI=IDINT(R8(NE+8,I))
       NNL=IDINT(R8(NE+6,I))
      WRITE(6,90) I,UN,LLI(I),NNL,SPI,LAM(I),I
      WRITE(*,90) I,UN,LLI(I),NNL,SPI,LAM(I),I
  90  FORMAT(5I6,X,E15.8,I6)
  80  CONTINUE
C     ====GET WAVE FUNCTIONS OF EXCITED STATES====
      NN1=0
      AS=S
      DO 95 I=1, NES
      LNO=11+I-1
      REWIND LNO
      READ (LNO) (DAT(J),J=1,15)
      M3(I)=DAT(13)+DAT(12)
      M4(I)=DAT(12)
      WRITE(6, 121)I
      WRITE(*, 121)I
  121 FORMAT('=========================='/
     *       '  EXCITED  STATES CHANNEL ',I6)
      WRITE(6, 126)
  126 FORMAT('      ENERGY   PHASE OR N ',
     *        '  IN ZONE')
      DO 130 J1=1, M3(I)
      LA=AS+J1
       CALL DISKO1(LNO,J1,2,NE5,R8(1,LA))
      DIRM=IDINT(R8(NE+8,LA))
      IA=7
      IF(R8(1,LA).LT.0.D 00) IA=5
       WRITE(6, 140)J1,R8(1,LA),R8(NE+IA,LA),LA
       WRITE(*, 140)J1,R8(1,LA),R8(NE+IA,LA),LA
  140  FORMAT(' ',I3,2F11.5,I6)
       IF(I.EQ.1) LK=R8(2,LA)
C------to correct phase-----
        IF(R8(1,LA).GT.1.)THEN
        FA=R8(NE+IA,LA)-FA
         IF(FA.GT.3.0)THEN 
          R8(NE+IA,LA)=R8(NE+IA,LA)-PI 
          WRITE(6, 141)J1,R8(1,LA),R8(NE+IA,LA),FA                          
  141  FORMAT('    ',I3,3F10.5)
         END IF
        END IF
        FA=R8(NE+IA,LA)
C----------------------
  130  CONTINUE
       AS=LA
       WRITE(6, 150)R8(2,LA)
       WRITE(*, 150)R8(2,LA)
  150  FORMAT(' WITH ELECTRON ORBITAL MOMENTUM ',F4.1)
       WRITE(6, 151)R8(NE+9,LA),R8(NE+10,LA)
       WRITE(*, 151)R8(NE+9,LA),R8(NE+10,LA)
  151  FORMAT(' AND TERM 2S+1,L ',2F6.1)
       WRITE(6, 153)DIRM
       WRITE(*, 153)DIRM
  153  FORMAT(' SPIN ',I3)

   95  CONTINUE
C      STOP                       	
C-----k1 = core functions with Li=lk---------
       K1=0
       IF(MU.EQ.1.D 00)THEN
       J=1
       DO 94 I=1,S
       IF( LLI(I).EQ.LK) THEN
C     WRITE(6,92) LLI(I),LK,K1
C     WRITE(*,92) LLI(I),LK,K1
C 92  FORMAT(' li,lk k1 ',3I4)
       K1=K1+1
       KK1(J)=I
       LAMK(J)=LAM(I)
       J=J+1
       END IF
  94  CONTINUE
       END IF
C ---k1 is equal kk from input data----
        K1=KK
      WRITE(6,91) K1,LK
      WRITE(*,91) K1,LK
  91  FORMAT(' No OF CORE FUNCTIONS ',I4,' WITH L=Lk ',I5)
       K=K+K1
C   now k=k+kk  for electron scat only!-----
C=====================K-Particles====
C      WRITE(*, 198)
      WRITE(6,198)
  198 FORMAT('    K-PARTICLE '/
     *       '     l   ENERGY     PHASE   SPIN ',
     *        '  IN ZONE')
       DO 200 IWK=1,K
       IF (IWK.LE.K1) THEN
        NK=KK1(IWK)
       ELSE
        NK=KZ+IWK-K1-1
       END IF
       CALL COLUMN(NE5,R8(1,NK),R5)
       LK=IDINT(R8(2,NK))
       LAMK(IWK)=R8(1,NK)
       DIRK=IDINT(R8(NE+8,NK))
       PHASE(IWK)=R8(NE+7,NK)
       IA=7
       IF(R5(1).LT.0.D 00) IA=5
       WRITE(6,201) IWK,LK,LAMK(IWK),R5(NE+IA),DIRK,NK
       WRITE(*,201) IWK,LK,LAMK(IWK),R5(NE+IA),DIRK,NK
  201  FORMAT(2I3,2F11.5,2I4)
  200  CONTINUE
C      STOP
       NA=0
C========cycle on k============================
       DO 202 IW=K1+1,K
       RE2=0.D 00
       IM2=0.D 00
       RE3=0.D 00
       IM3=0.D 00
       RE4=0.D 00
       IM4=0.D 00
       RE5=0.D 00
       IM5=0.D 00
        IF(LAMK(IW).LT.EMIN) GOTO 202
        IF(LAMK(IW).GT.EMAX) GOTO 550
       WRITE(6, 101)IW,DIRK,LAMK(IW)
       WRITE(*, 101)IW,DIRK,LAMK(IW)
 101   FORMAT('=== K - SCATTERING PARTICLE ',I3,
     *  ' SPIN ',I3,' Ek ',E15.8)
C      STOP
       DO 203 I=1,K
       DO 203 U1=1,K
       RRR(U1,I)=0.D 00
       RRM(U1,I)=0.D 00
 203   CONTINUE
C ======cycle on number of diagrams========
       DO 96 JN=1,ND
       WRITE(6, 152)JN,IW,LAMK(IW),LK,LV(JN)
       WRITE(*, 152)JN,IW,LAMK(IW),LK,LV(JN)
  152  FORMAT(' =',I3,' K',I3, ' Ek ',E15.8,' Lk',I3,
     *       ' TRANSFERED ORBITAL MOMENTUM ',I3)
      IF( EXPI.EQ.0) THEN
      LAML=LAM(IL(JN))
      ELSE
      LAML=LAM(IL(JN))-IPOS
C      LAML=LAML1
      WRITE(6, 154)
      WRITE(*, 154)
 154  FORMAT (' CORRECTED ')
      END IF
      LL=LLI(IL(JN))
      SPI=IDINT(R8(NE+8,IL(JN)))
      WRITE(6, 157)
      WRITE(*, 157)
 157  FORMAT(' TIME FORWARD TERM: F-PARTICLE ')
      DO 160 I=1,F
      NF=FZ(JN)+I-1
      CALL COLUMN(NE5,R8(1,NF),R5)
      LF=IDINT(R8(2,NF))
      LAMF(I)=R8(1,NF)
      DIRF=IDINT(R8(NE+8,NF))
      IF(I.EQ.1)
     * WRITE(6, 161)LF,LAMF(I),R5(NE+7),DIRF,NF
      IF(I.EQ.1)
     * WRITE(*, 161)LF,LAMF(I),R5(NE+7),DIRF,NF
 161  FORMAT(4X,' L',I4,' E',F10.5,' PHASE',F10.5,
     *    ' SPIN ',I4, ' ZONE',I5)
 160  CONTINUE
      IF (F.GT.F1) DELXF=
     *DSQRT(LAMF(F))-DSQRT(LAMF(F-1))
C-------------------------------
      WRITE(6, 155)IL(JN),LL,SPI,LAML
      WRITE(*, 155)IL(JN),LL,SPI,LAML
 155  FORMAT(' CORE ORBITAL L-HOLE',I3,' L=',I3,
     *    ' SPI ', I3,' E',F10.5)
      WRITE(6, 162)
      WRITE(*, 162)
 162  FORMAT(' EXCITED STATES ORBITALS: M-ELECTRON  ')
      DO 170 I=1,M
      NM=MZ(JN)+I-1
C      WRITE(6,172)NM
C 172  FORMAT('NM ',I4)
       CALL COLUMN(NE5,R8(1,NM),R5)
      LM=IDINT(R8(2,NM))
      LAMM(I)=R8(1,NM)
      DIRM=IDINT(R8(NE+8,NM))
      IF(I.EQ.1)
     * WRITE(6, 173)LM,LAMM(I),R5(NE+7),DIRM,NM
      IF(I.EQ.1)
     * WRITE(*, 173)LM,LAMM(I),R5(NE+7),DIRM,NM
 173  FORMAT(4X,' L',I4,' E',F10.5,' PHASE',F10.5,
     *    ' SPIN ',I4,' ZONE',I5)
 170  CONTINUE
      IF (M.GT.M1) DELXM=
     *DSQRT(LAMM(M))-DSQRT(LAMM(M-1))
C+++++++++++++++++++++++++++++++++++++++++++
       DO 214 IN=1,NI 
       LAMI=LAM(IT(IN))
       LI=LLI(IT(IN))
       DIRI=IDINT(R8(NE+8,IT(IN)))
       IF(LI.EQ.LF.AND.DIRK.EQ.DIRI)THEN
        WRITE(6, 335)IT(IN),LI,DIRI,LAMI
        WRITE(*, 335)IT(IN),LI,DIRI,LAMI
 335   FORMAT(' TIME REVERSE I-HOLE',I3,' L=',I3,
     *  ' SPIN ',I3,' E ', F10.5)
	  END IF
 214    CONTINUE
C++++++++++++++++++++++++++++++++++++++++++++++++
C==================================
       DO 205 IWK=1,K
       IF (IWK.LE.K1) THEN
        NK=KK1(IWK)
       ELSE
        NK=KZ+IWK-K1-1
       END IF
      CALL MAEL4(ALFA,BET,NE5,H,NK,FZ(JN),F,MZ(JN),M,
     *   IL(JN),LL,LF,LM,LK,
     *   LV(JN),R1,PR,ME3,ME1,ME2,ME3,R8,R2,R5,R6,R7)
      DO 180 I=1,F
      DO 181 J=1,M
      MME1(J,I,IWK)=ME1(J,I)
C--if spin dirk non=dirm--------------------
      IF(DIRM.NE.DIRK) ME2(J,I)=0.
      MME2(J,I,IWK)=ME2(J,I)
  181  CONTINUE
  180  CONTINUE
C-------------------------------------
C       IF (W1.EQ.IW.AND.W2.EQ.IW)THEN
C       WRITE(*, 199)IW
C       WRITE(6, 199)IW
C 199  FORMAT('================='/I4,
C     *  '  COULOMB MATRIX ELEMENTS (M*F)  '/
C     *    '=====TIME FORWARD===========')
C       END IF
C       DO 210 I=1,F
C        DO 212 J=1,M
C       IF(I.EQ.1.AND.J.EQ.1)
C     *  WRITE(6, 211)J,I,ME1(J,I),ME2(J,I),ME3(J,I)
C  212   CONTINUE
C  210   CONTINUE
C  211  FORMAT(2I3,' DIR.F',F11.5,' EXCH.F',F11.5,' COMB.F',F11.5)
       DO 213 I=1,M
        ME4(I)=0.
        ME5(I)=0.
 213  CONTINUE
C+++++++++++++++++++++++++++++++=
       DO 212 IN=1,NI 
       LAMI=LAM(IT(IN))
       LI=LLI(IT(IN))
       DIRI=IDINT(R8(NE+8,IT(IN)))
C++++++++++++++++++++++++++++++++++=
        IF(DIRK.EQ.DIRI) THEN
      IF (LV(JN).LE.LI+LK.AND.LV(JN).GE.IABS(LK-LI).
     *    AND.MU.EQ.1.D 00) THEN
C         WRITE(6, 214)IWK,IW,JN,LI,LF
C 214   FORMAT('===IWK,IW,JN=LI,LF=',5I3)
C+++++++++++++++++++++++++++++++++
        IF(LI.EQ.LF)THEN
C        IF (IWK.EQ.JN) THEN
C        WRITE(6, 335)IT(IN),LI,DIRI,LAMI
C        WRITE(*, 335)IT(IN),LI,DIRI,LAMI
C 335   FORMAT(' TIME REVERSE I-HOLE',I3,' L=',I3,
C     *  ' SPIN ',I3,' E ', F10.5)
C        END IF
C+++++++++++++++++++++++++++++++++++++++++++++
       CALL MAEL5(ALFA,BET,NE5,H,IL(JN),IT(IN),NK,MZ(JN),M,
     *  LL,LI,LK,LM,
     * LV(JN),R1,PR,ME,ME4,ME5,ME6,ME7,R8,R2,R5,R6,R7)
       END IF
C++++++below  <end li=lf>   ------
       END IF
C        WRITE(6, 221)
C 221   FORMAT('=====TIME REVERSE (M)======')
  223     DO 220 I=1,M
       MME4(I,IWK)=ME4(I)
C---------if dirk non=dirm------------
       IF(DIRM.NE.DIRK) ME5(I)=0.   
       MME5(I,IWK)=ME5(I)
C       IF(I.EQ.1)
C     *  WRITE(6, 223)I,ME4(I),ME5(I),ME6(I)
C  223   FORMAT(I6,' DIR.R',F11.5,' EXCH.R',F11.5,' COMB.R',F11.5)
 220   CONTINUE
C       END IF
C+++++++++++++++++++++++++++++++++++
        END IF
  212   CONTINUE
 205   CONTINUE
C======================================
C       STOP
C=====================================
       DO 204 W2=1,K
       DO 206 W1=1,K
       B=1./DFLOAT(2*LV(JN)+1)/DFLOAT(2*LK+1)
C------time forward - direct---------
       DO 320 I=1, M
       OM=LAMK(IW)-LAMM(I)+LAML
       DO 310 U1=1, F
       RP(U1)=MME1(I,U1,W2)*MME1(I,U1,W1)*4.0
       IF(DIRK.EQ.0) RP(U1)=RP(U1)*2.
C       WRITE(6, 316)OM,LAMM(I),LAMK(IW),LAMF(U1),LAML,
C     *       U1,RP(U1),MME1(I,U1,W2),MME1(I,U1,W1)
  316  FORMAT(' OM,LAMMKFL ',5E15.8/
     *     ' RP,MME1,MME2 ',I3,3E15.8)
 310   CONTINUE
      CALL ONE2(OM,F,F1,LAMF,CC,QQ,DD1,IMIN,RP,RM,RE1,IM1)
C       IF (W1.EQ.IW.AND.W2.EQ.IW.AND.LV(JN).EQ.0)THEN
C       DO 315 III=1,F
C       WRITE(6, 317)CC(III),QQ(III),DD1(III),RM(III),IMIN
C  317  FORMAT('ONE C',E15.8,'Q ',E15.8,' D',E15.8, 'INT',E15.8
C     *          ,'IMIN ',I3)
C  315  CONTINUE     
C       END IF
C------- no pole in discrete states for positron-----
        IF(MU.EQ.-1.D 00.AND.I.LE.M1)IM1=0.0
       RE(I)=RE1
       IM(I)=IM1
       IF ( I.GT.M1 ) THEN
       A=2.0*DSQRT(LAMM(I))
       RE(I)=RE(I)*A
       IM(I)=IM(I)*A
       END IF
 320   CONTINUE
       IF (W1.EQ.IW.AND.W2.EQ.IW)THEN
C       WRITE(6, 324)JN,I,RE(I),IM(I)
 324   FORMAT(' RES1 ND',I3,' M ',I3,' RE ',E15.8,' IM ',E15.8)
       END IF
       SUM=  SIMPS1(M1+1,M,DELXM,RE)
       SUM1= SIMPS1(M1+1,M,DELXM,IM)
      DO 330 I=1,M1
       SUM=SUM+RE(I)
       SUM1=SUM1+IM(I)
 330  CONTINUE
C       IF (W1.EQ.IW.AND.W2.EQ.IW.AND.LV(JN).EQ.1)THEN
C       WRITE(6, 323)B,SUM,SUM1,IMIN
C 323   FORMAT('Simps B,Re,Im   '/3E15.8,' IMIN',I3)
C       END IF
       RE1=SUM*B
       IM1=SUM1*B
       IF(MU.EQ.1.D 00)THEN
C ----time forward exchange--only for electron-----
      DO 350 I=1, M
       OM=LAMK(IW)-LAMM(I)+LAML
      DO 340 U1=1, F
       RP(U1)=MME1(I,U1,W2)*MME2(I,U1,W1)*4.0
 340  CONTINUE
      CALL ONE2(OM,F,F1,LAMF,CC,QQ,DD1,IMIN,RP,RM,RE2,IM2)
       RE(I)=RE2
       IM(I)=IM2
      IF ( I.GT.M1 ) THEN
       A=DSQRT(LAMM(I))*2.
       RE(I)=RE(I)*A
       IM(I)=IM(I)*A
       END IF
 350   CONTINUE
       IF (W1.EQ.IW.AND.W2.EQ.IW) THEN
C     DO 346 I=1,M
C     IF(I.EQ.1)
C     WRITE(6, 347)JN,I,RE(I),IM(I)
C      WRITE(*, 347)JN,I,RE(I),IM(I)
 347   FORMAT('  RES2 ND',I3,' F ',I3,' RE ',E15.8,' IM ',E15.8)
C346  CONTINUE
       END IF
        SUM= SIMPS1(M1+1,M,DELXM,RE)
        SUM1=SIMPS1(M1+1,M,DELXM,IM)
      DO 360 I=1, M1
       SUM= SUM+RE(I)
       SUM1=SUM1+IM(I)
 360   CONTINUE
       RE2=SUM*B
       IM2=SUM1*B
       SUMRE3=0.
       SUMRE4=0.
C---- begin of time reverse---only for electron-----
       DO 361 IN=1,NI 
       LAMI=LAM(IT(IN))
       LI=LLI(IT(IN))
       DIRI=IDINT(R8(NE+8,IT(IN)))
C       IF (W1.EQ.IW.AND.W2.EQ.IW) THEN
C       WRITE(6, 369)IN,LAMI
C       WRITE(*, 369)IN,LAMI
C 369   FORMAT(' REVERSE,IN,LAMI ',I3,E15.8)
C       end if
C-------spin condition---
        IF(DIRK.EQ.DIRI) THEN
       IF(LV(JN).LE.(LI+LK).AND.LV(JN).GE.IABS(LK-LI))THEN
C------time reverse -direct---------
      DO 370 I=1, M
       RP(I)=MME4(I,W2)*MME4(I,W1)*2.0
     * /(LAMK(IW)-LAMI-LAML+LAMM(I))
       IF(DIRK.EQ.0)  RP(I)=RP(I)*2.
      IF ( I.GT.M1 ) THEN
       A=DSQRT(LAMM(I))*2.
       RP(I)=RP(I)*A
       END IF
 370  CONTINUE
       RE3=  SIMPS1(M1+1,M,DELXM,RP)
       SUMRE3=SUMRE3-RE3*B
C---time reverse  exchange------------
      DO 380 I=1, M
       RP(I)=MME4(I,W2)*MME5(I,W1)*2.0
     * /(LAMK(IW)-LAMI-LAML+LAMM(I))
      IF ( I.GT.M1 ) THEN
       A=DSQRT(LAMM(I))*2.
       RP(I)=RP(I)*A
       END IF
 380  CONTINUE
       RE4=  SIMPS1(M1+1,M,DELXM,RP)
       SUMRE4=SUMRE4-RE4*B
        END IF 
        END IF
 361   CONTINUE
C---end of time reverse-----------
        END IF
C------end mu=1-----------------------------
        RE3=SUMRE3
        RE4=SUMRE4
        IF (W1.EQ.IW.AND.W2.EQ.IW) THEN
        WRITE(6, 365)RE1,IM1
        WRITE(*, 365)RE1,IM1
 365   FORMAT(
     *   ' TIME FORWARD '/
     *   ' DIRECT   RE1 ',E15.8,' IM1 ',E15.8)
        WRITE(6, 366)RE2,IM2,RE3,RE4
        WRITE(*, 366)RE2,IM2,RE3,RE4
 366   FORMAT(
     *   ' EXCHANGE RE2 ',E15.8,' IM2 ',E15.8/
     *   ' TIME REVERSE   '/
     *   ' DIRECT   RE3 ',E15.8/
     *   ' EXCHANGE RE4 ',E15.8)
       END IF
 395   SUM=RE1-RE2+RE3-RE4
       SUM1=IM1-IM2
        A=-PI*SUM
        B=-PI*SUM1
       IF(IW.EQ.W1.AND.IW.EQ.W2) THEN
       RE5=RE5+A
       IM5=IM5+B
       WRITE(6, 389) RE5,IM5
       WRITE(*, 389) RE5,IM5
       END IF
  389  FORMAT(
     *   ' COMBINED RE5 ',E15.8,' IM5 ',E15.8)
      RRR(W2,W1)=RRR(W2,W1)+SUM
      RRM(W2,W1)=RRM(W2,W1)+SUM1
 206  CONTINUE
 204  CONTINUE
C===end w2,w1=========================
 96   CONTINUE
       WRITE(6, 390)RE5,IM5
  390  FORMAT(
     *   ' A and B for SHIFT in n-diagrams '/
     *   ' RE ',E15.8,' IM ',E15.8)
C     STOP
C=========end JN=-nD-diagrams===========
         OM1=LAMK(IW)
         DOP=0.
	 IF (IW.GT.1)DOP=LAMK(IW-1) 
          OM1=OM1+(-1)**NA*0.001*DABS(LAMK(IW)-DOP)
      CALL CQD1(OM1,K,K1,LAMK,CC,QQ,DD1)
        NA=NA+1
C       WRITE(6, 399)IW,LAMK(IW),K,K1
C 399   FORMAT(' W,LAMK(W),K,K1 ',I3,E15.8,2I3/
C     *       ' I     CC       QQ       DD1   RR   RM  ')
C       DO 397 IJ=1,K
C       WRITE(6,398)IJ,CC(IJ),QQ(IJ),DD1(IJ),RRR(1,IJ),RRM(1,IJ)
C 397   CONTINUE
C 398   FORMAT(I3,2X,5(2X,E12.6))
       CALL INVMA1(IW,K,SR,Y1,RRR,RRM,CC,QQ,DD1,RE1,IM1)
      RE1=-RE1*PI
      IM1=-IM1*PI
       WRITE(6, 491)RE1,IM1
       WRITE(*, 491)RE1,IM1
 491   FORMAT('A and B for SCAT'/
     *       E15.8,' Re ',E15.8,' IM ',E15.8)
C--phase shift  for SHIFT--with DATAN----------------
       A3=RE5**2+IM5**2
       A2=0.5*DATAN(2.*RE5/(1.-A3))
       A4=DSIGN(1.D 00,RE5)
        WRITE(6, 398)A2,RE5,A3
 398    FORMAT('shift DATAN/2,A,A^2+B^2'/3E15.8)
       IF (DABS(RE5).LT.1..AND.A3.LT.1.) RE3=A2
       IF (DABS(RE5).EQ.1..AND.A3.EQ.1.) RE3=PI/2.*A4
       IF (DABS(RE5).GT.1..AND.A3.GT.1.) RE3=A2+PI/2.*A4
       IM3=DABS(0.5*DLOG(((1.+IM5)**2+RE5**2)/
     *   DSQRT((1.-A3)**2+4.*RE5**2)))
C-----phase shift for SCAT---------------------
 	 A3=RE1**2+IM1**2
       A1=0.5*DATAN(2.*RE1/(1.-A3))
	 A4=DSIGN(1.D 00,RE1)
        WRITE(6, 399)A1,RE1,A3
 399    FORMAT('scat DATAN/2,A,A^2+B^2'/3E15.8) 
       IF (DABS(RE1).LT.1..AND.A3.LT.1.) RE2=A1
       IF (DABS(RE1).EQ.1..AND.A3.EQ.1.) RE2=PI/2.*A4
       IF (DABS(RE1).GT.1..AND.A3.GT.1.) RE2=A1+PI/2.*A4   
       IM2=DABS(0.5*DLOG(((1.+IM1)**2+RE1**2)/
     *   DSQRT((1.-A3)**2+4.*RE1**2)))
C-------------------------------------
      SUM2=PHASE(IW)+RE3
      SUM3=PHASE(IW)+RE2
      PHRESH(IW)=SUM2
      PHIMSH(IW)=IM3
      PHRESC(IW)=SUM3
      PHIMSC(IW)=IM2
C----elastic cros sect in HF, SHIFT, SCAT-------------------
      SIGM1(IW)=4.*PI/LAMK(IW)*DFLOAT(2*LK+1)*DSIN(PHASE(IW))**2
      RESI2=2.*PI/LAMK(IW)*DFLOAT(2*LK+1)*(1.-DCOS(2.*SUM2))
      RESI3=2.*PI/LAMK(IW)*DFLOAT(2*LK+1)*(1.-DCOS(2.*SUM3))
      IMSI2=2.*PI/LAMK(IW)*DFLOAT(2*LK+1)*(1.-DEXP(-4.*IM3))
      IMSI3=2.*PI/LAMK(IW)*DFLOAT(2*LK+1)*(1.-DEXP(-4.*IM2))
      SIGM2(IW)=2.*PI/LAMK(IW)*DFLOAT(2*LK+1)*(DCOSH(2.*IM3)
     *    -DCOS(2.*SUM2))*DEXP(-2.*IM3)
      SIGM3(IW)=2.*PI/LAMK(IW)*DFLOAT(2*LK+1)*(DCOSH(2.*IM2)
     *    -DCOS(2.*SUM3))*DEXP(-2.*IM2)
C-- inelastic cros section in SHIFRT, SCAT------------------
      SIGM4(IW)=PI/LAMK(IW)*DFLOAT(2*LK+1)*(1.-DEXP(-4.*IM3))
      SIGM5(IW)=PI/LAMK(IW)*DFLOAT(2*LK+1)*(1.-DEXP(-4.*IM2))
C-----------------------------------------------------------
       WRITE(6, 490) IW,Z,LAMK(IW),LK,
     *    PHASE(IW),SIGM1(IW),
     *    RE3,IM3,SUM2,IM3,RESI2,IMSI2,SIGM2(IW),
     *    RE2,IM2,SUM3,IM2,RESI3,IMSI3,SIGM3(IW),SIGM4(IW),SIGM5(IW)
       WRITE(*, 490) IW,Z,LAMK(IW),LK,
     *    PHASE(IW),SIGM1(IW),
     *    RE3,IM3,SUM2,IM3,RESI2,IMSI2,SIGM2(IW),
     *    RE2,IM2,SUM3,IM2,RESI3,IMSI3,SIGM3(IW),SIGM4(IW),SIGM5(IW)
 490   FORMAT('*************K-Particle************',I4/
     *  ' Z=',F7.3,'  Ek=',F11.5,' lk=',I3/
     *  ' HF PHASE  ',E15.8/
     *  ' HF CROSS SECTION  ',E15.8,' a02'/
     *   '-------------------------------------'/
     *  ' IN THE SECOND-ORDER APPROXIMATION '/
     *  ' PHASE SHIFT '/
     *  ' RE ',E15.8,' IM ',E15.8/
     *  ' HF PHASE+ PHASE SHIFT '/
     *  ' RE ',E15.8,' IM ',E15.8/
     *  ' CROSS SECTION   Re       Im       Total (a02) '/
     *  '     ',3E15.8/
     *   '-------------------------------------'/
     *  ' WITH POLARIZATION CORRECTIONS '/
     *  ' PHASE SHIFT '/
     *  ' RE ',E15.8,' IM ',E15.8/
     *  ' HF PHASE+ PHASE SHIFT '/
     *  ' RE ',E15.8,' IM ',E15.8/
     *  ' CROSS SECTION   Re       Im       Total (a02) '/
     *  '     ',3E15.8/
     *   '-------------------------------------'/
     *  '          CROSS SECTIONS  (ao2) (inelastic) '/
     *  '    SHIFT              SCAT    '/2E15.8/
     *   '*************************************')
 202  CONTINUE
 550   WRITE(6, 500)
 500   FORMAT('====================================='/
     *  '                 PHASE SHIFTS  '/
     *  ' PARTICLE           2-ORDER      POLARISATIONAL'/
     *  '  ENERGY            APPROXIMATION    CORRECTIONS'/
     *  '   ERy       HF   SHIFTRE    SHIFIM    SCATRE    SCATIM')
      DO 510 I=K1+1, K       
       IF(LAMK(I).LT.EMIN.OR.LAMK(I).GT.EMAX ) GOTO 510
       WRITE(6, 505) LAMK(I),PHASE(I),PHRESH(I),PHIMSH(I),
     *	    PHRESC(I),PHIMSC(I)
 505   FORMAT(F7.4,5F10.4)
 510   CONTINUE
       WRITE(6, 520) LK,Z
 520  FORMAT('====================================='/
     *  ' PARTICLE WITH L ',I3,' --> ATOM WITH Z= ',F7.3/
     *  '          CROSS SECTIONS  (ao2) (elastic) '/
     *  ' PARTICLE              2-ORDER       POLARISATIONAL'/
     *  '  ENERGY            APPROXIMATION     CORRECTIONS'/
     *  '     ERy        HF          SHIFT              SCAT')
      DO 530 I=K1+1, K
       IF(LAMK(I).LT.EMIN.OR.LAMK(I).GT.EMAX ) GOTO 530
       WRITE(6, 535) LAMK(I),SIGM1(I),SIGM2(I),SIGM3(I)
 535   FORMAT(F7.4,3E15.8)
 530   CONTINUE
       WRITE(6, 540)
 540   FORMAT('====================================='/
     *  '          CROSS SECTIONS  (ao2) (inelastic) '/
     *  ' PARTICLE       2-ORDER       POLARISATIONAL '/
     *  '  ENERGY     APPROXIMATION     CORRECTIONS  '/
     *  '     ERy         SHIFT              SCAT    ')
      DO 545 I=K1+1, K       
       IF(LAMK(I).LT.EMIN.OR.LAMK(I).GT.EMAX ) GOTO 545
      WRITE(6, 555) LAMK(I),SIGM4(I),SIGM5(I)
 555   FORMAT(F7.4,2F15.8)
 545   CONTINUE
 600  RETURN
      END