
C==== PROGRAMM FOR CALCULATION OF BREMSSTRAHLUNG SPECTRA=====
C =======OF SLOW ELECTRONS ABOVE THRESHOLD=================
C       authors: L.V.CHERNYSHEVA,M.Ya.AMUSIA,N.B.AVDONINA,
C               M.Yu.KUCHIEV
C      FORTRAN VERSION: L.V.CHERNYSHEVA
C============================================================
       SUBROUTINE BREM(NES,HI,IT,MU,NU1,NU2,EXPI,LAMI1,P,NE,S,NW1,
     *     LI,NNL,U1,LAMI,IZ,M,M1,MZ,LM,LAMM,XM0,DELXM,
     *     ME,FMEA2,RR,FME,FM,FMEA,FME1,FMEA1,P1,P2,
     *     OM,OMT,SIGM,SIGM1,INT,RE,IME,Q1,Q2,FQR,
     *     CC,QQ,D1,RM,
     *     R1,R2,R5,R6,R8)
      INTEGER NES,NU1,NU2,EXPI,AS,
     *        A,LV,L,S,IT,HI,TAU,
     *        NNL(S),LI(S),U1(S),IZ(S),
     *        M(NES),M1(NES),MZ(NES),LM(NES)
       REAL*8  ALFA,BET,H,Z,P,OM0,RO1,R,
     *    DOM,A2,MU,A1,RE1,IM1,SUM,SIMPS1,
     *    OM,OM1,B,DQ2,LAMI1,Y1,XMAXM,Y,DELY,
     *    LAMI(S),LAMM(NES*NW1),XM0(NES),DELXM(NES),
     *    ME(S),FMEA2(HI,S),RR(NW1),FME(NES*NW1),FM(NW1),
     *    FMEA(NES*NW1),FME1(NES*NW1),FMEA1(NES*NW1),P1(HI,NW1),
     *    P2(HI,NW1),OMT(NW1),SIGM(NW1),SIGM1(NW1),
     *    INT(HI),RE(NW1),IME(NW1),Q1(HI),
     *    Q2(HI),FQR(HI),DAT(15),CC(NW1),QQ(NW1),D1(NW1),RM(NW1),
     *    R1(NE+15),R2(NE+15),R5(NE+15),R6(NE+15),R8(NE+15,*)

       WRITE(6,10 )
       WRITE(*,10 )
  10   FORMAT('  BREMSSTRAHLUNG OF FAST ELECTRONS  '/
     *        '         ABOVE THRESHOLD  ')
      READ(10) (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)
      WRITE(6,30)Z,S,NE,R,H,BET,ALFA,P,NU1,NU2,MU
  30  FORMAT(' LS-COUPLING HF'/
     *       ' CHARGE OF NUCLEUS           ',F7.1/
     *       ' NUMBER OF SHELLS            ',I7/
     *       ' NUMBER OF POINTS            ',I7/
     *       ' R MAX                       ',F7.1/
     *       ' STEP OF INTEGRATION         ',E15.6/
     *       ' BETA                        ',E15.6/
     *       ' ALPHA                       ',E15.8/
     *       ' IMPULSE OF INCIDENT ELECTRON',E15.6/
     *       ' NU1                         ',I7/
     *       ' NU2                         ',I7/
     *       ' MASS OF INCIDENT PARTICLE   ',E15.8)
      WRITE(6,40)
      WRITE(*,40)
  40  FORMAT('==============================='/
     *       ' GROUND STATE CONFIGURATION'/' ORBITAL   N ',
     *'     L     Q      ENERGY          ZONE')
      DO 50 I=1, S
      CALL DISKO1(10,I,2,NE5,R8(1,I))
          LAMI(I)=R8(1,I)
         LI(I)=R8(2,I)
         U1(I)=R8(NE+5,I)
         NNL(I)=R8(NE+6,I)
      WRITE(6,51)I,U1(I),LI(I),NNL(I),LAMI(I),I
      WRITE(*,51)I,U1(I),LI(I),NNL(I),LAMI(I),I
   51 FORMAT(' ',4I6,1X,E15.8,I6)
   50 CONTINUE

          IF(EXPI.NE.0) THEN
            LAMI(IT)=LAMI1
            WRITE(6, 55)
            ELSE
            WRITE(6, 57)
          END IF
      WRITE(6, 56)U1(IT),LI(IT),NNL(IT),LAMI(IT)
 57   FORMAT('==============THEORETICAL ENERGY  ')
 55   FORMAT('==============EXPERIMENTAL ENERGY ')
   56 FORMAT(' HOLE N ',I3,' L ',I3,' Q ',I3,' E= ',E15.6)

C ========GET WAVEFUNCTIONS OF EXCITED STATES==========
         NN1=0
         AS=S
         DO 60 I=1,NES
         LNO=11+I-1
         REWIND LNO
         READ (LNO) (DAT(J),J=1,15)
C        DAT(13)=11
         M(I)=DAT(13)+DAT(12)
         NN1=NN1+M(I)
         M1(I)=DAT(12)
         MZ(I)=AS+1
         XM0(I)=DAT(14)
         DELXM(I)=DAT(15)
      WRITE(6,65)I
      WRITE(*,65)I
  65  FORMAT( 60('=')/' EXCITED STATE   CHANNEL ',I6/
     *       '      ENERGIES           PHASE',
     *       ' OR N    IN ZONES')
      DO 80 J=1,M(I)
      LA=AS+J
      CALL DISKO1(LNO,J,2,NE5,R8(1,LA))
      R8(NE+2,LA)=0.0D 00
      LAMM(J)=R8(1,LA)
C      DK1(I,J)=R8(NE+7,LA)
      IA=7
      IF(R8(1,LA).LT.0.D 00) IA=5
      WRITE(6,95)J,R8(1,LA),R8(NE+IA,LA),LA
      WRITE(*,95)J,R8(1,LA),R8(NE+IA,LA),LA
  95  FORMAT('  ',I2,2X,E15.8,1X,E15.8,I6)
  80  CONTINUE
      LM(I)=R8(2,LA)
      AS=LA
      WRITE(6,105) R8(2,LA)
      WRITE(*,105) R8(2,LA)
  105 FORMAT(' WITH ELECTRON ORBITAL MOMENTUM ', F4.1)
   60  CONTINUE
C=======================================================
      DO 70 J=1 , M(1)
        OM=LAMM(J)-LAMI(IT)
       OMT(J)=OM
  70   CONTINUE
       Y=DLOG(P-DSQRT(P**2-OMT(M1(1)+1)*MU))
       DELY=DLOG(2.*P/(P-DSQRT(P**2-OMT(M1(1)+1)*MU)))/
     * DFLOAT(HI-1)
C       PRINT 85,Y,DELY
C  85   FORMAT('  Y,DY ',2E15.8)
      DO 90  I=1, HI
        Q2(I)=Y
        Y=Y+DELY
  90  CONTINUE
      N1=M(1)
      N4=M1(1)
      N2=LI(IT)
      N3=LM(1)
      N6=MZ(1)-S
      WRITE(6, 91)
      WRITE(*, 91)
  91  FORMAT('    LN MOM.TRANS.  MOM. TRANS. FORMFACTOR     ')
      DO 96 IW=1, HI
        Q1(IW)=DEXP(Q2(IW))
      CALL FORMFA(ALFA,BET,NE5,H,S,Q1(IW),FQR(IW),
     *            LI,NNL,ME,R1,R5,R6,R8)
      WRITE(6,120)IW,Q2(IW),Q1(IW),FQR(IW)
 120  FORMAT(1X,I4,3E15.8)
  96  CONTINUE
      DO 100 IW=1, HI
       WRITE(6,110)Q2(IW),Q1(IW)
       WRITE(*,110)Q2(IW),Q1(IW)
  110  FORMAT('================================'/
     *        '  LN MOMENTUM TRANSFERED TO ATOM',E15.8/
     *        '     MOMENTUM TRANSFERED TO ATOM',E15.8)
       WRITE(6,122)
       WRITE(*,122)
 122   FORMAT('   BESSEL MATRIX ELEMENTS ')
      CALL MAEL3(ALFA,BET,NE5,H,N1,N2,N3,1,Q1(IW),S+1,IT,
     *       FME1(N6),FMEA1(N6),R1,R2,R5,R6,R8)
      DO 121 J=1,M(1)
      WRITE(6,130)FME1(J)
 121  CONTINUE
 130   FORMAT(2X,E15.8)
C  100  CONTINUE
      IF ( NU2.EQ.0 ) GOTO 190
      IF (IW.EQ.1) THEN
      WRITE(6,102)
 102   FORMAT('  DIPOLE  MATRIX ELEMENTS ')
      CALL MAEL1(ALFA,BET,NE5,H,S+1,IT,N1,N2,N3,1,
     *     FME(N6),FMEA(N6),R1,R2,R5,R6,R8)
      DO 101 J=1,M(1)
      WRITE(6,103)FME(J),FMEA(J)
      WRITE(6,103)FME(J),FMEA(J)
 103   FORMAT('  R-FORM ',E15.8,'  V-FORM ',E15.8)
 101  CONTINUE
      END IF
      WRITE(6,139)
      WRITE(*,139)
 139  FORMAT(' POLARIZ.AMPLITUDE  ')
      DO 140 J=M1(1)+1 , M(1)
       B=OMT(J)**2
      DO 160 I=1, M(1)
       FM(I)=FME(I)*FME1(I)
       RR(I)=FM(I)/(OMT(J)-LAMI(IT)+LAMM(I))
      IF( I.GT.M1(1) ) THEN
       A1=DSQRT(LAMM(I))*2.
       RR(I)=RR(I)*A1
      END IF
 160  CONTINUE
       OM1=OMT(J)+LAMI(IT)
C       DO 162 I=1,M(1)
C       PRINT 163,LAMM(I),FME(I),FME1(I),FM(I),RR(I)
C  162  CONTINUE
C  163  FORMAT('  LAMM FME  FME1  FM   RR      '/5F10.5)
       CALL ONE1(OM1,LAMM,FM,N1,N4,RE1,IM1,CC,QQ,D1,RM)
       P1(IW,J)=RE1
       P2(IW,J)=IM1
       SUM=SIMPS1(N1-N4,DELXM(1),RR,1)
      DO 170 I=1, M1(1)
        SUM=SUM+RR(I)
 170   CONTINUE
C        PRINT 171,RE1,IM1,SUM
C  171  FORMAT(' RE1,IM1,SUM  '/3E15.8)
       P1(IW,J)=(P1(IW,J)-SUM)*B
       P2(IW,J)=P2(IW,J)*B
       WRITE(6,180)OM1,P1(IW,J),P2(IW,J)
 180   FORMAT(' OM= ',E15.8,' RE= ',E15.8,
     *        '  IM= ',E15.8)
 140  CONTINUE
 100  CONTINUE
 190   A1=32./3./137.04**3/P**2

       B=0.
      DO 200 I=M1(1)+1, M(1)
C      IW=0
       Y=DLOG(P-DSQRT(P**2-OMT(I)*MU))
      DO 210 J=1, HI
       DELY=DLOG(2.*P/(P-DSQRT(P**2-OMT(I)*MU)))/
     * DFLOAT(HI-1)
C      IF ( (DABS(Q2(J)-Y).LT.0.00001).OR. ((Q2(J)-Y.LT.DELY/2)
C     * .AND.(Q2(J)-Y.GT.0)).OR.  ((Y-Q2(J).LT.DELY/2).AND.(Y
C     * -Q2(J).GT.0)) ) IW=J
C      PRINT 201,NU1,NU2,Z,FQR(J),P1(J,I),P2(J,I)
C201  FORMAT('  nu1,nu2,z,fqr,p1,p2  '/
C    *        2X,2I2,4E15.8)
      INT(J)=((DFLOAT(NU1)*(Z-FQR(J))+MU*DFLOAT(NU2)*P1(J,I))**2+
     *  MU**2*DFLOAT(NU2)**2*P2(J,I)**2)*A1
C      PRINT 191,I,J,Q2(J),Y,DELY,INT(J)
C 191  FORMAT(' I,J,Q2(J),Y,DELY,INT(J)  '/
C    *               X,2I2,3F10.4,E15.8)
 210  CONTINUE
C      IF (IW.EQ.0 ) PRINT 220,J
C 220   FORMAT(' ERROR',I5)
      SUM= SIMPS1(HI,DELY,INT,1)
C      IF (DABS(Q2(IW)-Y).GT.0.00001) B=((DFLOAT(NU1)*(Z-FQR(IW))
C     *    +MU*DFLOAT(NU2)*P1(IW,I))**2+MU**2*DFLOAT(NU2)**2*P2(IW,I)
C     * **2)*(Y-Q2(IW))*A1
C       PRINT 221,SUM
C 221   FORMAT(' SUM',E15.8)
C       SUM=(SUM-B)/OMT(I)
       SIGM(I)=SUM/OMT(I)
C       PRINT 240,SIGM(I)
C  240  FORMAT(E15.8)
  200  CONTINUE
       B=0
       A1=A1/2.
C       PRINT 249
C 249   FORMAT('  ANGULAR DISTRIBUTION')
      DO 250 I=M1(1)+1, M(1)
C       IW=0
       Y=DLOG(P-DSQRT(P**2-OMT(I)*MU))
       DELY=DLOG(2.*P/(P-DSQRT(P**2-OMT(I)*MU)))/
     * DFLOAT(HI-1)
      DO 260 J=1, HI
        A2=3.*(MU*OMT(I)/2./P/Q1(J)+Q1(J)/2./P)**2-1.
C      IF ( (DABS(Q2(J)-Y).LT.0.00001).OR. ((Q2(J)-Y.LT.DELY/2.)
C     * .AND.(Q2(J)-Y.GT.0.)).OR.  ((Y-Q2(J).LT.DELY/2.).AND.(Y
C     * -Q2(J).GT.0.)) ) IW=J
      INT(J)=((DFLOAT(NU1)*(Z-FQR(J))+MU*DFLOAT(NU2)*P1(J,I))**2+
     *  MU**2*DFLOAT(NU2)**2*P2(J,I)**2)*A1*A2
 260  CONTINUE
C      IF ( IW.EQ.0 ) PRINT 270
C 270   FORMAT(' ERROR  ' )
       SUM=SIMPS1(HI,DELY,INT,1)
C      IF ( DABS(Q2(IW)-Y).GT.0.00001) B=((DFLOAT(NU1)*(Z-FQR(IW))+
C     *  MU*DFLOAT(NU2)*P1(IW,I))**2+MU**2*DFLOAT(NU2)**2*P2(IW,I)
C     *   **2)*(Y-Q2(IW))*A1*A2
C       SUM=(SUM-B)/OMT(I)
       SUM=SUM/OMT(I)
       SIGM1(I)=SUM/SIGM(I)
C       PRINT 280,SIGM1(I)
C 280   FORMAT(E15.8)
 250   CONTINUE
C =======    TABLE OF RESULTS===================
       WRITE(6,300)Z,S,P,LAMI(IT),NU1,NU2,MU,HI
  300  FORMAT( '=========================='/
     *        '  CHARGE OF NUCLEUS        ',E15.8/
     *        '  NUMBER OF SHELLS         ',I7/
     *        '  IMPULSE OF INC. PARTICLE ',E15.8/
     *        '  HOLE ENERGY              ',E15.8/
     *        '  NU1                      ',I7/
     *        '  NU2                      ',I7/
     *        '  MASS OF INCIDENT PARTICLE',E15.8/
     *        '  NUMBER OF POINTS Q       ',I7/
     * '=PHOTON ENERGY  SPECTRUM  ANG.DISTRIBUTION',
     * ' CROSS SECTION (A.E.) ')
      DO 310 I=M1(1)+1, M(1)
      A2=SIGM(I)*OMT(I)/2.
       WRITE(6,320)OMT(I),SIGM(I),SIGM1(I),A2
       WRITE(*,320)OMT(I),SIGM(I),SIGM1(I),A2
 320  FORMAT(4E15.8)
 310  CONTINUE
      RETURN
       END

