display_hxr_yp

R5-X2 - Display HXR results

 Calculate and display the HXR bremsstrahlung from results of the 
 Electron Drift Kinetic solver (main_dke_yp.m)

 by Joan Decker <joan.decker@cea.fr> (CEA/IRFM) and Yves Peysson <yves.peysson@cea.fr> (CEA/IRFM)

0001 function display_hxr_yp
0002 %R5-X2 - Display HXR results
0003 %
0004 % Calculate and display the HXR bremsstrahlung from results of the
0005 % Electron Drift Kinetic solver (main_dke_yp.m)
0006 %
0007 % by Joan Decker <joan.decker@cea.fr> (CEA/IRFM) and Yves Peysson <yves.peysson@cea.fr> (CEA/IRFM)
0008 %
0009 close all
0010 warning off;
0011 %
0012 if nargin == 0,    
0013     dir('RESULTS_*.mat');
0014     filename = [];
0015     filename = input_dke_yp('Kinetic calculations filename ',filename);
0016 %
0017     dir('HXR_*.mat');
0018     filenamehxr = [];
0019     filenamehxr = input_dke_yp('HXR tomographic system filename ',filenamehxr);
0020 %
0021     auto_mode_in = input_dke_yp('Automatic display mode (yes: 1, no: 0) ',0);
0022     igor_mode = input_dke_yp('Igor software output (yes: 1, no: 0) ',0);%For figure by figure data output (Easy data transfer to Igor software for final output)
0023 elseif nargin == 1,
0024     auto_mode_in = input_dke_yp('Automatic display mode (yes: 1, no: 0) ',0);
0025     igor_mode = input_dke_yp('Igor software output (yes: 1, no: 0) ',0);%For figure by figure data output (Easy data transfer to Igor software for final output)
0026 elseif nargin == 2,
0027     igor_mode = 0;%Automatic mode when test_dke_1yp.m is running
0028 end
0029 %
0030 load(filename);%Load DKE output data
0031 load(filenamehxr);%Load HXR structure
0032 %
0033 display_mode = input_dke_yp('Display graphics (yes: 1, no: 0)',1) + 1;
0034 %
0035 if nargin >= 2,
0036     auto_mode = auto_mode_in;
0037 end
0038 %
0039 filenameroot = filename(1:findstr(filename,'_out.mat')-1);
0040 %
0041 disp('-------------------------------------------------------------------------------------------------------------------');
0042 disp('Display results of Drift Kinetic Equation solver (see test_dke_1yp.m)');
0043 disp('-------------------------------------------------------------------------------------------------------------------');
0044 %
0045 nchords = length(hxr.alpha_hxr);
0046 theta_display = 180;
0047 chord_display = 1;
0048 kphot_display = 70;%keV
0049 IPBT = 1;%Cosine between the plasma current and the toroidal B-field directions
0050 while 1,
0051     ir_display_dke = find(r_dke == ir_display);
0052     if auto_mode_in == 0,
0053         disp(' ')
0054         disp('*****************************************************************************************************************************')
0055         disp('*************************************WARNING: from tokamak topview, BT is clockwise !****************************************')
0056         IPBT = input_dke_yp(['*********************************** Sign of the scalar product BT.Ip (-1 or +1) (exit: <=0)'],IPBT);
0057         disp('************************************************************************** For local output only ****************************')
0058         ir_display_dke = input_dke_yp(['*********************************** Radial index value (1-',int2str(nr_dke),') (exit: <=0)'],ir_display_dke);
0059         theta_display = input_dke_yp(['*********************************** Poloidal angle value (0-180 deg) (exit: <=0)'],theta_display);
0060         chord_display = input_dke_yp(['*********************************** Chord index value (0-',int2str(nchords),') (exit: <=0)'],chord_display);
0061         kphot_display = input_dke_yp(['*********************************** Photon energy (20-200 keV) (exit: <=0)'],kphot_display);
0062         disp('*****************************************************************************************************************************')
0063         disp(' ')
0064         if ir_display_dke < 1,
0065             break
0066         end
0067     end
0068 %
0069     [vschords,psin_hxr,theta_hxr,PSI_hxr,mask_hxr] = bremchord_dke_yp(equil,IPBT,hxr,xpsin_dke,xpsin_S_dke,1,'spline',display_mode);%Calculation of the chord length and angle between magnetic field line and line-of-sight
0070 %
0071     [kphot,brem_bd_out,brem_forward,brem_backward,brem_perp,brem_4pi,lchord,rho_out,theta_out] = bremsstrahlung_dke_yp(equil,hxr,vschords,psin_hxr,theta_hxr,PSI_hxr,mask_hxr,XXheaviside,zZi,xzni_norm(:,r_dke),xne_norm(r_dke),ne_ref,betath_ref,xrho(r_dke),dmhu,mhu,dpn_ref,pn_ref,pn2_ref,gamma_ref,Ec_ref,dp_cyl,pnmax_S,XXfinit_ref(:,:,r_dke),XXf0_ref(:,:,r_dke),XXf0_tp_ref(:,:,r_dke),XXf0_g_ref(:,:,r_dke),nlevel,bounce_mode,display_mode,xrho(ir_display_dke),theta_display);%Local and line-integrated bremsstrahlung calculation
0072 %
0073     [kinterp,gkEinterp,kfit,dkfit,brem_bd_plasma,brem_bd_diag,tphot_plasma,afit_plasma,tphot_plasma_exp,afit_plasma_exp,meanabsx] = bremdiag_dke_yp(equil,hxr,kphot,brem_bd_out,display_mode,chord_display);%Convolution with the diagnostic system response function for each chord
0074 %
0075 %  Display results
0076 %
0077     figure('Name','Plasma bremsstrahlung spectrum'),
0078     graph1D_jd(kfit,brem_bd_plasma(chord_display,:),0,1,'Photon energy k (keV)','dN_{phot}/dk.dt [keV^{-1}.s^{-1}]',['Non-thermal Bremsstrahlung, chord #',num2str(chord_display)],NaN,'','','-','non','b',2);
0079 
0080 %
0081     figure('Name','Pulse bremsstrahlung spectrum'),
0082     graph1D_jd(kfit,brem_bd_diag(chord_display,:),0,1,'Pulse energy channel E (keV)','dN_{pulse}/dE.dt [keV^{-1}.s^{-1}]',['HXR diagnostic count rate, chord #',num2str(chord_display)],NaN,'','','-','non','b',2);
0083     
0084 %
0085     ikphot_display = max(find(kfit <= kphot_display));
0086 %
0087     figure('Name','Plasma bremsstrahlung profile'),
0088     graph1D_jd([1:length(tphot_plasma)],brem_bd_plasma(:,ikphot_display),0,0,'Chord number','dN_{phot}/dk.dt [keV^{-1}.s^{-1}]',['HXR diagnostic count rate, k=',num2str(kfit(ikphot_display)),' (keV)'],NaN,'','','-','non','b',2);
0089 
0090 %
0091     figure('Name',''Pulse bremsstrahlung profile'),
0092     graph1D_jd([1:length(tphot_plasma)],brem_bd_diag(:,ikphot_display),0,0,'Chord number','dN_{pulse}/dE.dt [keV^{-1}.s^{-1}]',['HXR diagnostic count rate, E=',num2str(kfit(ikphot_display)),' (keV)'],NaN,'','','-','non','b',2);
0093 %
0094     figure('Name','Plasma photon temperature'),
0095     graph1D_jd([1:length(tphot_plasma)],tphot_plasma,0,0,'Chord number','T_{phot} [keV]','k.dN_{phot}/dk.dt = A.exp(-k/T_{phot}) [50-110 keV]',NaN,'','','-','non','b',2);
0096 %
0097     figure('Name','Measured photon temperature'),
0098     graph1D_jd([1:length(tphot_plasma)],tphot_plasma_exp,0,0,'Chord number','T_{phot}^{*} [keV]','E.dN_{pulse}/dE.dt.{\eta}_{E} = A.exp(-E/T_{phot}^{*}) [50-110 keV]',NaN,'','','-','non','r',2);
0099 %
0100 %    figure(10015),
0101 %    graph1D_jd([1:length(tphot_plasma)],afit_plasma,0,0,'Chord number','A(k=0) [s^{-1}]','k.dN_{phot}/dk.dt = A.exp(-k/T_{phot})@k=0 [50-110 keV]',NaN,'','','-','non','b',2);
0102 %
0103 %    figure(10016),
0104 %    graph1D_jd([1:length(tphot_plasma)],afit_plasma_exp,0,0,'Chord number','A^{*}(E=0) [s^{-1}]','E.dN_{pulse}/dE.dt.{\eta}_{E} = A.exp(-E/T_{phot}^{*})@E=0 [50-110 keV]',NaN,'','','-','non','r',2);
0105 %
0106     figure('Name','Mean HXR absorption efficiency'),
0107     graph1D_jd(kfit,meanabsx,0,0,'Photon energy channel (keV)','{\eta}_{E}','Mean stopping efficiency {\eta}_{E}',NaN,'','','-','non','r',2);
0108 %
0109     figure('Name','Detector HXR response'),
0110     ikinterp_display = max(find(kinterp <= kphot_display));
0111     graph1D_jd(kinterp,gkEinterp(ikinterp_display,:),0,0,'Photon energy k','G_{kE}',['Detector response function @k=',num2str(kinterp(ikinterp_display)),' (keV), chord#',int2str(chord_display)],NaN,'','','-','non','r',2);
0112     faxis = axis;dfaxis_x = (faxis(2) - faxis(1))/100;dfaxis_y = (faxis(4) - faxis(3))/100; 
0113     text(faxis(1) + dfaxis_x*5,faxis(4) - dfaxis_y*5,['Norm. = ',num2str(trapz(kinterp,gkEinterp(ikinterp_display,:)))]);
0114     
0115 %
0116 %*********************************** Saving results ***************************************
0117 %
0118     filename = ['HXR_',tokname,'_S',int2str(snumber),'.mat'];
0119     if diary_mode == 1,%Automatic data saving when diary mode is called
0120         filename = ['HXR_',tokname,'_S',int2str(snumber),'.mat'];        
0121         snumber=snumber+1;
0122         del = 1;
0123     else,%Manual data saving
0124         while 1,
0125             if exist(filename) == 2 & auto_mode == 0,%Only in manual mode
0126                 del = input_dke_yp(['File ',filename,' already exists ! Do you want to overwrite it (yes:1, no:0)'],0);
0127                 if del == 1,
0128                     break;
0129                 else
0130                     snumber = input_dke_yp('New simulation number',snumber+1);
0131                     filename = ['HXR_',tokname,'_S',int2str(snumber),'.mat'];    
0132                 end
0133             else
0134                 del = 1;
0135                 break;
0136             end
0137         end
0138     end
0139     %
0140     if del == 1    ,
0141         delete(filename);
0142         save(filename);
0143         if exist(filename) == 2,
0144             disp('-------------------------------------------------------------------------------------------------------------------');
0145             info_dke_yp(2,['The HXR data file : ',filename,' has been successfully saved !']);
0146         else
0147             info_dke_yp(2,['The HXR data file : ',filename,' has not been saved, an unknown error occurs on the disk !']);    
0148         end
0149     end
0150 
0151 end