haug_dke_yp

    Differential electron-electron bremstrahlung cross-section (dsigma/dk.domega)
    Haug model including Coulomb correction factor.

    Input:

        - ec: kinetic energy of the incoming electron (keV) [1,m] 
          (default = 500)
        - e: photon energy (keV)  [1,n] 
        - t: angle between the speed direction of the incoming electron and 
          the direction of emission of the bremsstrahlung photon (radian) [1,p]
          (default = 0)

    Output: 

        - seHC: Haug + Coulomb factor bremstrahlung cross-section (m^2) [p,m,n]
        - seH: Haug bremstrahlung cross-section (m^2) [p,m,n]
        - fee: Elwert correction factor [p,m,n]
        - ec: kinetic energy of the incoming electron (mec2) [p,m,n] 
        - k: photon energy (mec2) [p,m,n] 
        - c: 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, seH or seHC must be divided 
             by 511 keV.

0001 function [seHC,seH,fee,ec,k,c] = haug_dke_yp(ec_in,k_in,t_in)
0002 %
0003 %    Differential electron-electron bremstrahlung cross-section (dsigma/dk.domega)
0004 %    Haug model including Coulomb correction factor.
0005 %
0006 %    Input:
0007 %
0008 %        - ec: kinetic energy of the incoming electron (keV) [1,m]
0009 %          (default = 500)
0010 %        - e: photon energy (keV)  [1,n]
0011 %        - t: angle between the speed direction of the incoming electron and
0012 %          the direction of emission of the bremsstrahlung photon (radian) [1,p]
0013 %          (default = 0)
0014 %
0015 %    Output:
0016 %
0017 %        - seHC: Haug + Coulomb factor bremstrahlung cross-section (m^2) [p,m,n]
0018 %        - seH: Haug bremstrahlung cross-section (m^2) [p,m,n]
0019 %        - fee: Elwert correction factor [p,m,n]
0020 %        - ec: kinetic energy of the incoming electron (mec2) [p,m,n]
0021 %        - k: photon energy (mec2) [p,m,n]
0022 %        - c: 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, seH or seHC 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 < 3,
0032     infoyp(2,'Wrong number of input arguments for haug_dke_yp');
0033     return;
0034 end
0035 %
0036 [qe,me,mp,mn,e0,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;
0054 %
0055 %Calculation parameters
0056 %
0057 x1 = k.*(e0-p0.*c);
0058 x2 = k;
0059 w2 = 2.*(e0+1);
0060 x = x1+x2;
0061 rho2 = w2-2.*x;
0062 w24 = w2-4;
0063 rho24 = rho2-4;
0064 w22 = w2-2;
0065 rho22 = rho2-2;
0066 %
0067 %Coulomb correction factor
0068 %
0069 betaw = sqrt(w2).*sqrt(w24)./w22;
0070 betarho = sqrt(rho2).*sqrt(rho24)./rho22;
0071 ap = alpha./betaw;
0072 a2 = alpha./betarho;
0073 fee = (a2./ap).*(exp(2.*pi.*ap)-1)./(exp(2.*pi.*a2)-1);
0074 %
0075 r1 = rho24+4*x1+4*(x1.^2)./rho2;
0076 r2 = rho24+2*x1;
0077 gw2 = sqrt(x2.*(rho24.*x2./4+2*x.*x1./rho2));
0078 gw4 = sqrt(w24.*(w24.*rho24./4+4*x1.*x2./rho2));
0079 xl = log((r2+sqrt(rho24.*r1)).*sqrt(rho2)./(4.*x1));
0080 xl1 = log((sqrt(rho2)+sqrt(rho24))./2);
0081 xl2 = log(1+rho2.*(rho24.*x2+2*sqrt(rho24).*gw2)./(4*x.*x1));
0082 xl3 = log(((sqrt(w2).*sqrt(rho24)+sqrt(rho2).*sqrt(w24)).^2)./(8*x));
0083 xl4 = log(1+rho2.*(w24.*rho24+2*sqrt(rho24).*gw4)./(8*x1.*x2));
0084 %
0085 %Main expression
0086 %
0087 s0 = sqrt(rho24);
0088 s1 = (w2+rho2).*(((x1-x2)./x).^2)./(4*x1.*x2);
0089 s2 = -0.25*((1./x1-1./x2).^2);
0090 s3 = -rho2./(2*(x.^2));
0091 s4 = 2*(1+1./w24).*(rho2./(w24.*x1.*x2));
0092 s5 = 4*rho2./(w2.*w24.*x1.*x1.*x1.*x1);
0093 s6 = 3*(w22.^2).*x1.*x2./w24;
0094 s7 = -2*(x.^2).*(1+6*1./(w2.*w24));
0095 s8 = -4*rho2.*x./(w24.*x1.*x1.*x1);
0096 s9 = rho2.*(4*1./w2-1.5-8*1./(w2.*(w24.^2))+w22.*x./(w2.*w24))./(x1.^2);
0097 s10 = -rho2./(w24.*x1);
0098 s11 = (w2-4*x2)./rho2;
0099 s12 = -w2.*w24./(4*x1.*x2);
0100 s13 = -4./x1;
0101 s14 = (w2-w22.*(rho2./2))./(x1.^2);
0102 %
0103 z1 = (s1+s2+s3+s4+s5.*(s6+s7)+s8+s9+s10+(s11+s12+s13+s14)./r1).*s0;
0104 %
0105 s0 = xl./sqrt(r1);
0106 s1 = 4+(10-w2./2)./x2;
0107 s2 = -2*(2.*w2-x2+4)./x1;
0108 s3 = 3*w2.*w24./(4*x1.*x2);
0109 s4 = 4.0./(x1.*r2);
0110 s5 = x2-3*x1+4+2*(rho2-3)./x1;
0111 s6 = (w24-4*x1.*x2./rho2)./(x1.*r1);
0112 s7 = w2.*r2./(4*x2)-rho22-2*x1;
0113 %
0114 z2 = s0.*(s1+s2+s3+s4.*s5+s6.*s7);
0115 %
0116 z3 = sqrt(rho2).*xl1.*((rho2+2)./(x.^2)+8*1./(x1.^2));
0117 %
0118 s0 = xl2./gw2;
0119 s1 = 2*w22./x2;
0120 s2 = -(rho22-x2)./x1;
0121 s3 = rho22.*(rho2+x1)./(2*x);
0122 s4 = w24./r2-2+rho22.*((w24+rho2).^2)./(8*x.*x1);
0123 s5 = 1./(r2.*x1);
0124 s6 = rho22.*(x2.*(w24+rho2)./2-w22)-2*(x2.^2)-4*x2;
0125 s7 = 1./(r2.*x);
0126 s8 = 0.5*rho22.*(3*rho24-w2.*(rho2-5))+2.*(x1.^2)-6*x1+w22.*x2;
0127 %
0128 z4 = s0.*(s1+s2+s3+s4+s5.*s6+s7.*s8);
0129 %
0130 s0 = 2*sqrt(rho2).*xl3./(sqrt(w2).*sqrt(w24).*x1);
0131 s1 = 4+8*w22./(w24.^2);
0132 s2 = -3*0.5*w22./x2;
0133 s3 = ((x2.^2)-2*w2-(w2-1).*x2+0.5*(w22.^2))./x1;
0134 s4 = (0.5*rho2.*rho22./x)-(x2.*w22.*0.5./x);
0135 s5 = -((w22.^2).*(rho2+w24))./(4*x.*x2);
0136 s6 = x./(4*x1.*x2);
0137 s7 = ((w22.^2)+rho22.*rho24)-8*rho22./w24;
0138 s8 = 4*(x.^2)./(w2.*w24.*x1);
0139 s9 = -(2*rho22.*x2+w24.*x-4*w22+w22.*(8-rho2)./(2*x1))./r2;
0140 s10 = (rho22.*0.5.*(3*rho24-w2.*(rho2-5))-x1.*(w2+4-2*x1))./(x.*r2);
0141 s11 = 2*(w2.*w22.*rho24-2*rho22+4*rho2./w2)./((w24.^2).*x1);
0142 s12 = 4*w22.*x./(w24.*(x1.^2));
0143 s13 = (12*x./(w2.*w24))-rho22./2;
0144 s14 = 1-x./(w2.*x1);
0145 %
0146 z5 = s0.*(s1+s2+s3+s4+s5+s6.*s7+s8+s9+s10+s11+s12.*s13.*s14);
0147 %
0148 s0 = (-1).*xl4./gw4;
0149 s1 = 1;
0150 s2 = rho22./(8.*x1.*x2);
0151 s3 = w22.^2+(rho22.^2)-6.*(w24+rho2)+16.*x./w24;
0152 s4 = 2.*(1-x1-(x1.^2))./(x1.*x2);
0153 s5 = (rho24-8./w24)./w24;
0154 %
0155 z6 = s0.*(s1+s2.*s3+s4+s5);
0156 %
0157 z7a = z1+z2+z3+z4+z5+z6;
0158 %
0159 x0 = x1;
0160 x1 = x2;
0161 x2 = x0;
0162 %
0163 r1 = rho24+4*x1+4*(x1.^2)./rho2;
0164 r2 = rho24+2*x1;
0165 gw2 = sqrt(x2.*(rho24.*x2./4+2*x.*x1./rho2));
0166 gw4 = sqrt(w24.*(w24.*rho24./4+4*x1.*x2./rho2));
0167 xl = log((r2+sqrt(rho24.*r1)).*sqrt(rho2)./(4.*x1));
0168 xl1 = log((sqrt(rho2)+sqrt(rho24))./2);
0169 xl2 = log(1+rho2.*(rho24.*x2+2*sqrt(rho24).*gw2)./(4*x.*x1));
0170 xl3 = log(((sqrt(w2).*sqrt(rho24)+sqrt(rho2).*sqrt(w24)).^2)./(8*x));
0171 xl4 = log(1+rho2.*(w24.*rho24+2*sqrt(rho24).*gw4)./(8*x1.*x2));
0172 %
0173 %Main expression...
0174 %
0175 s0 = sqrt(rho24);
0176 s1 = (w2+rho2).*(((x1-x2)./x).^2)./(4*x1.*x2);
0177 s2 = -0.25*((1./x1-1./x2).^2);
0178 s3 = -rho2./(2*(x.^2));
0179 s4 = 2*(1+1./w24).*(rho2./(w24.*x1.*x2));
0180 s5 = 4*rho2./(w2.*w24.*x1.*x1.*x1.*x1);
0181 s6 = 3*(w22.^2).*x1.*x2./w24;
0182 s7 = -2*(x.^2).*(1+6./(w2.*w24));
0183 s8 = -4*rho2.*x./(w24.*x1.*x1.*x1);
0184 s9 = rho2.*(4./w2-1.5-8./(w2.*(w24.^2))+w22.*x./(w2.*w24))./(x1.^2);
0185 s10 = -rho2./(w24.*x1);
0186 s11 = (w2-4*x2)./rho2;
0187 s12 = -w2.*w24./(4*x1.*x2);
0188 s13 = -4./x1;
0189 s14 = (w2-w22.*(rho2./2))./(x1.^2);
0190 %
0191 z1 = (s1+s2+s3+s4+s5.*(s6+s7)+s8+s9+s10+(s11+s12+s13+s14)./r1).*s0;
0192 %
0193 s0 = xl./sqrt(r1);
0194 s1 = 4+(10-w2./2)./x2;
0195 s2 = -2*(2.*w2-x2+4)./x1;
0196 s3 = 3*w2.*w24./(4*x1.*x2);
0197 s4 = 4./(x1.*r2);
0198 s5 = x2-3*x1+4+2*(rho2-3)./x1;
0199 s6 = (w24-4*x1.*x2./rho2)./(x1.*r1);
0200 s7 = w2.*r2./(4*x2)-rho22-2*x1;
0201 %
0202 z2 = s0.*(s1+s2+s3+s4.*s5+s6.*s7);
0203 %
0204 z3 = sqrt(rho2).*xl1.*((rho2+2)./(x.^2)+8./(x1.^2));
0205 %
0206 s0 = xl2./gw2;
0207 s1 = 2*w22./x2;
0208 s2 = -(rho22-x2)./x1;
0209 s3 = rho22.*(rho2+x1)./(2*x);
0210 s4 = w24./r2-2+rho22.*((w24+rho2).^2)./(8*x.*x1);
0211 s5 = 1./(r2.*x1);
0212 s6 = rho22.*(x2.*(w24+rho2)./2-w22)-2*(x2.^2)-4*x2;
0213 s7 = 1./(r2.*x);
0214 s8 = 0.5*rho22.*(3*rho24-w2.*(rho2-5))+2.*(x1.^2)-6*x1+w22.*x2;
0215 %
0216 z4 = s0.*(s1+s2+s3+s4+s5.*s6+s7.*s8);
0217 %
0218 s0 = 2*sqrt(rho2).*xl3./(sqrt(w2).*sqrt(w24).*x1);
0219 s1 = 4+8*w22./(w24.^2);
0220 s2 = -3*0.5*w22./x2;
0221 s3 = ((x2.^2)-2*w2-(w2-1).*x2+0.5*(w22.^2))./x1;
0222 s4 = (0.5*rho2.*rho22./x)-(x2.*w22.*0.5./x);
0223 s5 = -((w22.^2).*(rho2+w24))./(4*x.*x2);
0224 s6 = x./(4*x1.*x2);
0225 s7 = ((w22.^2)+rho22.*rho24)-8*rho22./w24;
0226 s8 = 4*(x.^2)./(w2.*w24.*x1);
0227 s9 = -(2*rho22.*x2+w24.*x-4*w22+w22.*(8-rho2)./(2*x1))./r2;
0228 s10 = (rho22.*0.5.*(3*rho24-w2.*(rho2-5))-x1.*(w2+4-2*x1))./(x.*r2);
0229 s11 = 2*(w2.*w22.*rho24-2*rho22+4*rho2./w2)./((w24.^2).*x1);
0230 s12 = 4*w22.*x./(w24.*(x1.^2));
0231 s13 = (12*x./(w2.*w24))-rho22./2;
0232 s14 = 1-x./(w2.*x1);
0233 %
0234 z5 = s0.*(s1+s2+s3+s4+s5+s6.*s7+s8+s9+s10+s11+s12.*s13.*s14);
0235 %
0236 s0 = (-1).*xl4./gw4;
0237 s1 = 1;
0238 s2 = rho22./(8.*x1.*x2);
0239 s3 = w22.^2+(rho22.^2)-6.*(w24+rho2)+16.*x./w24;
0240 s4 = 2.*(1-x1-(x1.^2))./(x1.*x2);
0241 s5 = (rho24-8./w24)./w24;
0242 %
0243 z6 = s0.*(s1+s2.*s3+s4+s5);
0244 %
0245 z7b = z1+z2+z3+z4+z5+z6;
0246 %
0247 z8 = z7a+z7b;
0248 seH = alpha*(re^2).*k.*z8./(pi*sqrt(w2).*sqrt(rho2).*sqrt(w24));
0249 seHC = alpha*(re^2).*k.*z8.*fee./(pi*sqrt(w2).*sqrt(rho2).*sqrt(w24));
0250 %
0251 fee(find(imag(fee)~=0)) = 0;
0252 seH(find(imag(seH)~=0)) = 0;
0253 seHC(find(imag(seHC)~=0)) = 0;
0254 seHC(isnan(seHC)) = 0;
0255 seH(isnan(seH)) = 0;
0256 seHC(isnan(seHC)) = 0;