bhe_dke_yp

    Differential electron-ion bremstrahlung cross-section (dsigma/dk.domega)
    Bethe-Heitler model including Elwert factor.

    Input:

        - ec: kinetic energy of the incoming electron (keV) [1,m] 
        - k: photon energy (keV)  [1,n] 
        - t: cosine of the angle between the direction of displacement of the incoming 
          electron and the photon emitted by bremsstrahlung (radian) [1,p]
          (default = 0)
        - Z: target ion charge (default = 1) [1,1]

    Output: 

        - seBHE: Bethe-Heitler + Elwert factor bremstrahlung cross-section (m^2) [p,m,n]
        - seBH: Bethe-Heitler bremstrahlung cross-section (m^2) [p,m,n]
        - Elwert: Elwert correction factor [p,m,n]
        - ec: kinetic energy of the incoming electron (mc2) [p,m,n] 
        - k: photon energy (mc2)  [p,m,n] 
        - c: cosine angle between the direction of displacement of the incoming electron and the photon emitted by bremsstrahlung (radian) [p,m,n]

    Warning: Cross-section units : m^2 but energies are in relativistic units. To get 
             cross-sections in standard m^2/keV units, seBH or seBHe must be divided 
             by 511 keV.


by Y.Peysson CEA-DRFC <yves.peysson@cea.fr> and J. Decker MIT-RLE <jodecker@mit.edu>

0001 function [seBHE,seBH,Elwert,ec,k,c] = bhe_dke_yp(ec_in,k_in,t_in,Z)
0002 %
0003 %    Differential electron-ion bremstrahlung cross-section (dsigma/dk.domega)
0004 %    Bethe-Heitler model including Elwert factor.
0005 %
0006 %    Input:
0007 %
0008 %        - ec: kinetic energy of the incoming electron (keV) [1,m]
0009 %        - k: photon energy (keV)  [1,n]
0010 %        - t: cosine of the angle between the direction of displacement of the incoming
0011 %          electron and the photon emitted by bremsstrahlung (radian) [1,p]
0012 %          (default = 0)
0013 %        - Z: target ion charge (default = 1) [1,1]
0014 %
0015 %    Output:
0016 %
0017 %        - seBHE: Bethe-Heitler + Elwert factor bremstrahlung cross-section (m^2) [p,m,n]
0018 %        - seBH: Bethe-Heitler bremstrahlung cross-section (m^2) [p,m,n]
0019 %        - Elwert: Elwert correction factor [p,m,n]
0020 %        - ec: kinetic energy of the incoming electron (mc2) [p,m,n]
0021 %        - k: photon energy (mc2)  [p,m,n]
0022 %        - c: cosine angle between the direction of displacement of the incoming electron and the photon emitted by bremsstrahlung (radian) [p,m,n]
0023 %
0024 %    Warning: Cross-section units : m^2 but energies are in relativistic units. To get
0025 %             cross-sections in standard m^2/keV units, seBH or seBHe must be divided
0026 %             by 511 keV.
0027 %
0028 %
0029 %by Y.Peysson CEA-DRFC <yves.peysson@cea.fr> and J. Decker MIT-RLE <jodecker@mit.edu>
0030 %
0031 if nargin < 4,
0032     infoyp(2,'Wrong number of input arguments for bhe_dke_yp');
0033     return;
0034 end
0035 %
0036 [qe,me,mp,mn,epsi0,mu0,re,mc2,clum,alpha] = pc_dke_yp;%Physics constant
0037 %
0038 ec = repmat(ones(length(t_in),1)*ec_in(:)'/mc2,[1,1,length(k_in)]);%Relativistic units
0039 k = shiftdim(repmat(ones(length(ec_in),1)*k_in(:)'/mc2,[1,1,length(t_in)]),2);%Relativistic units
0040 if (length(ec_in) == 1) & (length(t_in) == 1),%For the case ec_in and t_in are scalars
0041     k = reshape(k,1,1,length(k_in));
0042 end
0043 c = repmat(t_in(:)*ones(1,length(ec_in)),[1,1,length(k_in)]);
0044 %
0045 s = sqrt(1-c.^2);
0046 %
0047 mask = ec > k;%Photon are only emitted by electron of higher energies
0048 %
0049 e0 = ec+1;
0050 p0 = sqrt(e0.^2-1);
0051 v0 = p0./e0;
0052 p = sqrt((e0-k).^2-1);
0053 ep = e0-k;d0 = e0-c.*p0;
0054 ksi0 = alpha*Z*e0./p0;
0055 ksi1 = alpha*Z*ep./p;
0056 %
0057 Elwert = (ksi1./ksi0).*((1-exp(-2*pi*ksi0))./(1-exp(-2*pi*ksi1)));
0058 %
0059 q = sqrt(p0.^2+k.^2-2*p0.*k.*c);
0060 w1 = log((ep.*e0-1+p.*p0)./(ep.*e0-1-p.*p0));
0061 w2 = log((ep+p)./(ep-p));
0062 w3 = log((q+p)./(q-p));
0063 s1 = 8*(s.^2).*(2*e0.^2+1).*(p0.^(-2)).*(d0.^(-4));
0064 s2 = 2*(5*e0.^2+2*e0.*ep+3).*(p0.^(-2)).*(d0.^(-2));
0065 s3 = 2*(p0.^2-k.^2).*(q.^(-2)).*(d0.^(-2));
0066 s4 = 4*ep.*(p0.^(-2)).*(d0.^(-1));
0067 s5 = 4*e0.*(s.^2).*(3*k-(p0.^2).*ep).*(p0.^(-2)).*(d0.^(-4));
0068 s6 = 4*(e0.^2).*(e0.^2+ep.^2).*(p0.^(-2)).*(d0.^(-2));
0069 s7 = (2-2*(7*e0.^2-3*e0.*ep+ep.^2)).*(p0.^(-2)).*(d0.^(-2));
0070 s8 = 2*k.*(e0.^2+ep.*e0-1).*(p0.^(-2)).*(d0.^(-1));
0071 s9 = 4*w2.*(p.^(-1)).*(d0.^(-1));
0072 s10 = w3.*(p.^(-1)).*(q.^(-1));
0073 s11 = 4*d0.^(-2);
0074 s12 = 6*k.*d0.^(-1);
0075 s13 = 2*k.*(p0.^2-k.^2).*(q.^(-2)).*(d0.^(-1));
0076 s14 = w1.*(p.^(-1)).*(p0.^(-1));
0077 %
0078 s0 = s1-s2-s3+s4+s14.*(s5+s6+s7+s8)-s9+s10.*(s11-s12-s13);
0079 a = alpha*Z^2*re^2*p./(8*pi*k.*p0);%Cross-section in cm^2 unit
0080 %
0081 seBH = s0.*a.*mask;
0082 seBHE = seBH.*Elwert;
0083 %
0084 seBHE(isnan(seBHE)) = 0;
0085 seBHE(~isreal(seBHE)) = 0;
0086