Description of bremsstrahlung_dke

0001 function  [Zbremplasma] = bremsstrahlung_dke_yp(equilHXR,mksa,momentumDKE,Zmomcoef,XXf0,hxrparam,dkedisplay,fieldside)
0002 %
0003 % R5-X2 - Calculation of  the plasma bremsstrahlung emissivity
0004 %
0005 % Calculates the plasma bremsstrahlung emissivity.
0006 %
0007 %   INPUTS:
0008 %
0009 %       - equilHXR: LUKE magnetic equilibrium structure for HXR
0010 %       - mksa: MKSA structure
0011 %       - momentumDKE: LUKE momentum structure
0012 %       - Zmomcoef: additional LUKE momentum structure
0013 %       - XXf0: LUKE distribution function
0014 %       - hxrparam : Bremsstrahlung calculation and display parameter structure
0015 %       - dkedisplay: DKE display structure
0016 %       - fieldside: HFS or LFS limit
0017 %
0018 %   OUTPUTS:
0019 %
0020 %      - Zbremplasma: Plasma bremsstrahlung data structure
0021 %         * Zbremplasma.kphot : photon energies (keV) [1,nk]
0022 %         * Zbremplasma.brem_bd_out : photon number per unit energy (keV) and time (s) for each chord [nchord,nk]
0023 %         * Zbremplasma.lchord : chord length (m) [1,nchord]
0024 %
0025 % by Joan Decker (CEA/IRFM,joan.decker@cea.fr) and Yves Peysson(CEA/IRFM,yves.peysson@cea.fr)
0026 %
0027 if nargin < 8,
0028     fieldside = '';
0029 end
0030 if nargin < 7,
0031     dkedisplay.display_mode = 0;
0032 end
0033 if nargin < 6,
0034    error('Not enough input arguments in bremshtrahlung_dke_yp')
0035 end
0036 %
0037 if ~isfield(Zmomcoef,'Ec_ref'),% for older LUKE simulations
0038     if isfield(mksa,'Ec_ref'),
0039         Zmomcoef.Ec_ref = mksa.Ec_ref;
0040     else
0041         error('Kinetic energy data missing. Please correct manually.')
0042     end
0043 end    
0044 %
0045 ne_ref = mksa.ne_ref;
0046 betath_ref = mksa.betath_ref;
0047 %
0048 mhu = momentumDKE.mhu;
0049 dpn = momentumDKE.dpn;
0050 pn = momentumDKE.pn;
0051 %
0052 pn2 = pn.*pn;
0053 np = length(pn);
0054 %
0055 Ec = Zmomcoef.Ec_ref;
0056 gamma = Zmomcoef.gamma;
0057 %
0058 xpsin = equilHXR.xpsin;
0059 theta = equilHXR.mtheta;
0060 xrho = equilHXR.xrho;
0061 XB = sqrt(equilHXR.XBx.^2 + equilHXR.XBy.^2 + equilHXR.XBphi.^2);
0062 xB0 = equilHXR.xB0;%XB(:,1).';% old value, true only for circular FS
0063 xne = equilHXR.xne;
0064 xzni = equilHXR.xzni;
0065 zZi = equilHXR.zZi;
0066 nZ = length(zZi);
0067 %
0068 nr = length(xpsin);
0069 %
0070 xzni_norm = xzni/ne_ref;
0071 xne_norm = xne/ne_ref;
0072 %
0073 [~,~,~,~,~,~,~,mc2,clum] = pc_dke_yp;%Physics constant
0074 %
0075 method = 'spline';
0076 n_gauleg = hxrparam.n_gauleg;
0077 n_pol = n_gauleg;%Necessary for accurate reconstruction at ksi = -1
0078 epsi_gauleg = hxrparam.epsi_gauleg;
0079 kphot = hxrparam.kphot;
0080 nk = length(kphot);
0081 %
0082 theta_hxr = ones(nr,1)*(theta(1:end-1) + theta(2:end))/2;% poloidal location for bremsstrahlung calculation
0083 nt = size(theta_hxr,2);
0084 psin_hxr = xpsin.'*ones(1,nt);% radial location for bremsstrahlung calculation
0085 %
0086 Xr = sqrt(equilHXR.Xx.^2 + equilHXR.Xy.^2);
0087 r_hxr = interp1(theta,Xr.',theta_hxr(1,:)).';
0088 % %
0089 PSI_hxr = griddata(xpsin.',theta,XB.',psin_hxr,theta_hxr)./(xB0.'*ones(1,nt));%B/Bmin
0090 %
0091 ZZIB = zeros(nr,nt,n_pol+1,nk);%Legendre coefficients for the bremsstrahlung at all radii at poloidal positions (NaN if not used)
0092 %
0093 % Legendre polynomials and related coefficients
0094 %
0095 [x_leg,w_leg] = gauleg_dke_yp(n_gauleg,epsi_gauleg);%Gauss-Legendre abscissa and weigth calculations between -1 and + 1
0096 [Pm,cL] = leg_dke_yp(x_leg(:).',n_pol);%Legendre polynomial calculation
0097 %
0098 x_forward = 1;%Forward emission
0099 [Pm_forward] = leg_dke_yp(x_forward(:).',n_pol);%Legendre polynomial for the forward emission
0100 x_backward = -1;%Backward emission
0101 [Pm_backward] = leg_dke_yp(x_backward(:).',n_pol);%Legendre polynomial for the backward emission
0102 x_perp = 0;%Perpendicular emission
0103 [Pm_perp] = leg_dke_yp(x_perp(:).',n_pol);%Legendre polynomial for the perpendicular emission
0104 %
0105 Zbremplasma.xrho = xrho;
0106 Zbremplasma.r_hxr = r_hxr;
0107 Zbremplasma.theta_hxr = theta_hxr;
0108 Zbremplasma.x_leg = x_leg;
0109 Zbremplasma.w_leg = w_leg;
0110 Zbremplasma.Pm = Pm;
0111 Zbremplasma.cL = cL;
0112 %
0113 %Bremsstrahlung cross-section calculations
0114 %
0115 seBHE = zeros(n_gauleg,np,nk,nZ);
0116 for iZ = 1:nZ,
0117     seBHE(:,:,:,iZ) = bhe_dke_yp(Ec,kphot,x_leg,zZi(iZ));%Electron-ion bremsstrahlung (Bethe-Heitler-Elwert) (mec2 -> keV units)
0118 end
0119 seHC = haug_dke_yp(Ec,kphot,x_leg);%Electron-electron bremsstrahlung (Haug-Coulomb) (mec2 -> keV units)
0120 %
0121 Zbremplasma.seBHE = seBHE;
0122 Zbremplasma.seHC = seHC;
0123 %
0124 %Gauss-Legendre projection for the bremsstrahlung cross-sections
0125 %
0126 seBHE_c=(reshape(seBHE,n_gauleg,np*nk*nZ).'.*(w_leg*ones(1,np*nk*nZ)).')*Pm.';
0127 seHC_c=(reshape(seHC,n_gauleg,np*nk).'.*(w_leg*ones(1,np*nk)).')*Pm.';
0128 %
0129 %%%%%%%%%% For test only (beginning section) %%%%%%%%%%%%%%
0130 %
0131 Zbremplasma.seBHE1 = reshape((Pm.'*(seBHE_c.*(ones(np*nk*nZ,1)*cL.')).'),n_gauleg,np,nk,nZ);
0132 Zbremplasma.seHC1 = reshape((Pm.'*(seHC_c.*(ones(np*nk,1)*cL.')).'),n_gauleg,np,nk);
0133 Zbremplasma.epsiBHE = (seBHE-Zbremplasma.seBHE1)./seBHE;
0134 Zbremplasma.epsiBHE(isnan(Zbremplasma.epsiBHE)) = 0;
0135 Zbremplasma.epsiHC = (seHC-Zbremplasma.seHC1)./seHC;
0136 Zbremplasma.epsiHC(isnan(Zbremplasma.epsiHC)) = 0;
0137 %
0138 %%%%%%%%%% For test only (end section) %%%%%%%%%%%%%%
0139 %
0140 seBHE_c1 = repmat(reshape(seBHE_c.',n_pol+1,np,nk,nZ),[1,1,1,1,nr]);
0141 ni_c1 = shiftdim(repmat(xzni_norm,[1,1,n_pol+1,np,nk]),2);
0142 BHE_c1 = permute(sum(seBHE_c1.*ni_c1,4),[1,2,3,5,4]);%Sum over different species
0143 %
0144 seHC_c1 = repmat(reshape(seHC_c.',n_pol+1,np,nk),[1,1,1,nr]);
0145 ne_c1 = shiftdim(repmat(xne_norm',[1,n_pol+1,np,nk]),1);
0146 HC_c1 = seHC_c1.*ne_c1;
0147 %
0148 BREM_c1 = BHE_c1 + HC_c1;
0149 %
0150 pn_c1 = permute(repmat(pn(:),[1,n_pol+1,length(kphot),1]),[2,1,3,4]);
0151 gamma_c1 = permute(repmat(gamma(:),[1,n_pol+1,length(kphot),1]),[2,1,3,4]);
0152 %
0153 for ir = 1:nr,%radial position
0154     for it = 1:nt,%poloidal position
0155         if strcmp(fieldside,'HFS')
0156             PSI_hxr(ir,it) = max(PSI_hxr(ir,:));%high field side only
0157         elseif strcmp(fieldside,'LFS')
0158             PSI_hxr(ir,it) = 1;%low field side only
0159         end
0160         mhu_hxr = sign(mhu.').*sqrt(1-PSI_hxr(ir,it)*(1 - mhu.^2).');%mhu at a poloidal position that is different from Bmin
0161         mhu_hxr = mhu_hxr(1-PSI_hxr(ir,it)*(1 - mhu.^2) > 0);
0162         Zbremplasma.mhu_hxr{ir,it} = mhu_hxr;
0163         Zbremplasma.mhubounce2_hxr(ir,it) = 1 - PSI_hxr(ir,it)/max(PSI_hxr(ir,:));
0164         Xf0_hxr = XXf0(:,1-PSI_hxr(ir,it)*(1 - mhu.^2) > 0,ir);%Non-thermal distribution function at poloidal location
0165         %
0166         Xf0_i = interp1(mhu_hxr,Xf0_hxr.',x_leg,method).';%Interpolation of the distribution function on the Legendre polynomial pitch-angle grid
0167         %
0168         Xf0_c = (Xf0_i.*(w_leg*ones(1,np)).')*Pm.';%Calculation of the Legendre coefficients for local values of the distribution function
0169         %
0170         Xf0_c1 = permute(repmat(reshape(Xf0_c.',n_pol+1,np,1),[1,1,1,nk]),[1,2,4,3]);%Same as Xf0_c but expanded for all photon energies
0171         %
0172         dIB_c1 = pn_c1.*pn_c1.*pn_c1.*BREM_c1(:,:,:,ir).*Xf0_c1./gamma_c1;
0173         IB_c1 = 2*pi*integral_dke_jd(dpn,dIB_c1,2);%Momentum integration
0174         ZZIB(ir,it,:,:) = reshape(IB_c1,n_pol+1,nk)*ne_ref^2*betath_ref*clum/mc2;%MKSA units(s-1*m-3*kev-1*str-1)
0175         %
0176         if ~isfield(dkedisplay,'fhxr') || dkedisplay.fhxr,% calculate and save HXR distribution
0177             %
0178             Xf0_i1 = (Pm.'*(Xf0_c.*(ones(np,1)*cL.')).').';%Distribution function on the Legendre polynomials pitch-angle grid back calculated
0179             norm0_hxr = 2*pi*integral_dke_jd(dpn.',pn2.'.*trapz_dke_yp([mhu_hxr,Xf0_hxr.']).',1);
0180             normi_hxr = 2*pi*integral_dke_jd(dpn.',pn2.'.*trapz_dke_yp([x_leg,Xf0_i.']).',1);
0181             %
0182             Zbremplasma.Znorm0_hxr{ir,it} = norm0_hxr;
0183             Zbremplasma.Znormi_hxr{ir,it} = normi_hxr;
0184             Zbremplasma.ZXf0_hxr{ir,it} = Xf0_hxr;
0185             Zbremplasma.ZXf0_i{ir,it} = Xf0_i;
0186             Zbremplasma.ZXf0_i1{ir,it} = Xf0_i1;
0187         end
0188         %
0189         % Local bremsstrahlung properties (forward and backward with respect to B direction or p// since p// is taken positive)
0190         %
0191         Zbremplasma.brem_forward{ir,it} = permute(reshape(Pm_forward.'*(squeeze(ZZIB(ir,it,:,:)).*(cL*ones(1,nk))),1,nk),[2,3,1]);%Forward emission
0192         Zbremplasma.brem_backward{ir,it} = permute(reshape(Pm_backward.'*(squeeze(ZZIB(ir,it,:,:)).*(cL*ones(1,nk))),1,nk),[2,3,1]);%Backward emission
0193         Zbremplasma.brem_perp{ir,it} = permute(reshape(Pm_perp.'*(squeeze(ZZIB(ir,it,:,:)).*(cL*ones(1,nk))),1,nk),[2,3,1]);%Perpendicular emission
0194         Zbremplasma.brem_4pi{ir,it} = permute(reshape(squeeze(ZZIB(ir,it,1,:))/2,1,nk),[2,3,1]);%Mean 4pi emission per steradian
0195         %
0196         if dkedisplay.display_mode,
0197             info_dke_yp(2,['Bremsstrahlung integral #',int2str((ir-1)*nt + it),'/',int2str(nr*nt),' done']);
0198         end
0199     end
0200 end
0201 %
0202 Zbremplasma.ZZIB = ZZIB;
0203 Zbremplasma.kphot = kphot;
0204 %
bremsstrahlung_dke_yp

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
