irunhxr_jd

 This function calculates the hard X-ray emisssion from LUKE results

0001 function varargout = irunhxr_jd(opt_gui,workdir,basestr,output,external,id_simul,equil,dkeparam,dkedisplay,select)
0002 %
0003 % This function calculates the hard X-ray emisssion from LUKE results
0004 %
0005 if nargin < 10,
0006     select = struct;
0007 end
0008 %
0009 quitmsg = '-----> HXR calculation aborted';
0010 warning off
0011 %
0012 % load LUKE results if needed
0013 %
0014 if nargin < 9,
0015     %
0016     % data from LUKE_RESULTS_*.mat file
0017     %
0018     [luke_file,luke_path] = uigetfile('LUKE_RESULTS_*.mat','Please select the LUKE RESULTS file',output);
0019     %
0020     if luke_file == 0,
0021         disp('-----> No LUKE RESULTS file specified.')
0022         disp(quitmsg)
0023         return
0024     end  
0025     load([luke_path,luke_file],'id_simul','equil','dkeparam','dkedisplay','dke_out','equilDKE','radialDKE','momentumDKE','mksa','Zmomcoef');
0026     %
0027 elseif any(isfield(output,{'savefile','dke_out'})),
0028     %
0029     if isfield(output,'savefile'),% output saved in a file by run_lukert
0030         lukedata = load(output.savefile);
0031         output = lukedata.lukestructs.output;
0032         clear lukedata
0033     end        
0034     %
0035     % data from output structure
0036     %
0037     dke_out = output.dke_out;
0038     equilDKE = output.equilDKE;
0039     radialDKE = output.radialDKE;
0040     momentumDKE = output.momentumDKE;
0041     mksa = output.mksa;
0042     Zmomcoef = output.Zmomcoef;
0043     %
0044     if isempty(equil),
0045         equil = dke_out.equil;
0046     end        
0047     %
0048 end
0049 %
0050 if isfield(equil,'shotnum'),
0051     shotnum = equil.shotnum;
0052 elseif isfield(equil,'info') && isfield(equil.info,'shotnum'),
0053     shotnum = equil.info.shotnum;
0054 elseif isfield(equil,'id'),
0055     idsep = find(equil.id == '_');
0056     shotnum = equil.id(idsep(1)+1:idsep(2)-1); 
0057 else
0058     shotnum = '';
0059 end
0060 %
0061 % load HXR camera spec
0062 %
0063 if isfield(external,'param') && ~isempty(external.param),% parameters from external - iluke and irunluke
0064     hxrparam = external.param;
0065     %
0066     if isfield(external,'cam'),
0067         cam = external.cam;
0068     else
0069         cam = {};
0070     end
0071     %
0072     hxrspec = cam2hxrspec_jd(cam);
0073     %
0074 elseif strcmp(basestr,'TS') && ~isempty(shotnum) && exist('tsdate','file') && input_dke_yp('Do you want to load the HXR camera spec automatically (y/n)','y',{'y','n'}) == 'y',
0075     %
0076     [hxrspec,hxrparam] = get_hxrspec_forLUKE_jd('TS',str2num(shotnum));
0077     %
0078 elseif strcmp(basestr,'WEST') && ~isempty(shotnum) && exist('tsdate','file') && input_dke_yp('Do you want to load the HXR camera spec automatically (y/n)','y',{'y','n'}) == 'y',
0079     %
0080     [hxrspec,hxrparam] = get_hxrspec_forLUKE_jd('WEST',str2num(shotnum));
0081 else
0082     %
0083     hxr_file = which(['HXR_',basestr,'.mat']);
0084     [hxr_file,hxr_path] = igetfile_jd(opt_gui,'HXR_*.mat','Please select the HXR camera spec file',hxr_file);
0085     if hxr_file == 0,
0086         disp('-----> No HXR camera spec file specified.')
0087         disp(quitmsg)
0088         return
0089     end
0090     %
0091     load([hxr_path,hxr_file],'hxr','hxrparam');
0092     %
0093     hxrspec = hxr;
0094     clear hxr %to avoid name conflict
0095     %
0096 end
0097 %
0098 if isfield(hxrspec,'beta_hxr')
0099     nc = length(hxrspec.beta_hxr);% number of chords
0100 else
0101     nc = 0;% no camera defined
0102 end
0103 %
0104 % display options
0105 %
0106 if isfield(select,'p_opt'),
0107     p_opt = select.p_opt;
0108 else
0109     p_opt = input_dke_yp('Option for output figures? -1 (nothing), 0 (print), 1 (print and save), 2 (save)',-1,-1:2,'',[1,1]);
0110 end
0111 %
0112 if nc == 0,
0113     d_opt = 'n';
0114 elseif isfield(select,'d_opt'),
0115     d_opt = select.d_opt;
0116 else
0117     d_opt = input_dke_yp('Do you want to plot the chords? [y/n]','n',{'y','n'},'',[1,1]);
0118 end
0119 %
0120 if strcmp(d_opt,'y'),
0121     %
0122     figure(1),clf
0123     %
0124     plot(equil.ptx' + equil.Rp,equil.pty');hold on
0125     %
0126     L = [-2:0.1:2];%Chord length
0127     R = zeros(nc,length(L));
0128     Z = zeros(nc,length(L));
0129     for chord = 1:nc,
0130         R(chord,:) = hxrspec.R_hxr(chord) + L*sin(hxrspec.beta_hxr(chord))*cos(hxrspec.alpha_hxr(chord));
0131         Z(chord,:) = hxrspec.Z_hxr(chord) + L*cos(hxrspec.beta_hxr(chord));
0132     end
0133     plot(R',Z'),
0134     plot(R(1),Z(1),'*')
0135     hold off
0136     axis('equal')
0137     xlabel('R(m)');ylabel('Z(m)');title(basestr)
0138     print_jd(p_opt,[basestr,'_HXR'],1)
0139 end 
0140 %
0141 % load and treat HXR experimental results
0142 %
0143 if isfield(select,'n_opt'),% iluke
0144     n_opt = select.n_opt;
0145 elseif isfield(external,'hxr'),
0146     n_opt = input_dke_yp(['Option for hxr experimental data? 0 (structure external), 1 (data file), 2 (',basestr,' database) -1 (no exp. data)'],0,-1:2,'',[1,1]);
0147 else
0148     n_opt = input_dke_yp(['Option for hxr experimental data? 1 (data file), 2 (',basestr,' database) -1 (no exp. data)'],-1,[-1,1,2],'',[1,1]);
0149 end
0150 %
0151 if n_opt == 0,% structure external with irunluke
0152     %
0153     t = external.hxr.t;
0154     energy = external.hxr.energy;
0155     hxr = external.hxr.counts;
0156     %filter_thickness = external.hxr.filter_thickness;
0157     %filter_material = external.hxr.filter_material;
0158     %
0159     if isfield(external,'beta')
0160         beta = external.beta;
0161     else
0162         beta = NaN;
0163     end
0164     %
0165     if iscell(t),%old TCV data file
0166         t = t{ishot};
0167         energy = energy{ishot};
0168         hxr = hxr{ishot};
0169     end
0170     %
0171 elseif n_opt > 0,
0172     if isempty(shotnum)
0173         disp('-----> Shot number information missing');
0174         return
0175     end
0176     %
0177     if n_opt == 1,% hxr experimental data from file
0178         %
0179         [data_file,data_path] = igetfile_jd(opt_gui,'*.*','Please select the hxr experimental data file :');
0180         %
0181         if data_file == 0,
0182             disp(quitmsg)
0183            return
0184         end
0185         %
0186         load([data_path,data_file],'shot','energy','hxr','t','beta');%,'filter_thickness','filter_material'
0187         %
0188         ishot = find(shot == str2num(shotnum));
0189         %
0190         if isempty(ishot) || ~exist('hxr','var'),
0191             disp(['-----> hxr data not found for shot :',shotnum]);
0192             disp(quitmsg)
0193             return
0194         end
0195         %
0196         if length(ishot) > 1,
0197             disp(['-----> WARNING : Several entries are found for shot :',shotnum,', first entry selected']);
0198             ishot = ishot(1);
0199         end
0200         %
0201         clear data_path data_file
0202         %
0203         if iscell(t),%old TCV data file
0204             t = t{ishot};
0205             energy = energy{ishot};
0206             hxr = hxr{ishot};
0207         end
0208         %
0209         clear ishot
0210         %
0211     elseif n_opt == 2,% hxr experimental data from database
0212         %
0213         opt.R = input_dke_yp('Pile-up rejection (TS only)? (0) no, (1) yes)',1,[0,1],'',[1,1]);
0214         opt.V = input_dke_yp('Real time diagnostic (TS only)? (0) no, (1) yes)',1,[0,1],'',[1,1]);
0215         %
0216         [energy,hxr,t,beta] = collect_hxr_forLUKE_jd(basestr,str2num(shotnum),opt);
0217         %
0218     end
0219     %
0220     if ~exist('beta','var'),
0221         beta = NaN;
0222     end
0223     %
0224 else% no data or external with iluke
0225     %
0226     if isfield(external,'exp') && ~isempty(external.exp),
0227         t = external.exp.time;
0228         energy = external.exp.energy;
0229         hxr = external.exp.data;      
0230         %
0231         beta = NaN;
0232     else
0233         energy = NaN;
0234     end
0235 end
0236 %
0237 hxrspec.kfit = hxrparam.kphot(1:2:end);
0238 %
0239 if any(isnan(energy)),% no or useless exp. data
0240     %
0241     hxrexp = NaN;
0242     dt = 1;
0243     %
0244 else
0245     %
0246     % test angles between diagnostic prescription and database read, if applicable
0247     %
0248     if nc > 0,
0249         beta(beta > pi) = 2*pi - beta(beta > pi);% correct beta definition
0250         %
0251         if ~all(isnan(beta)),
0252             if length(beta) ~= nc,
0253                 error('number of chords mismatch')
0254             end
0255             if any(abs(beta - hxrspec.beta_hxr) > 1e-15),
0256                 error('angle of chords mismatch')
0257             end
0258         end
0259     end
0260     %
0261     % time sample of experimental data
0262     %
0263     if isfield(select,'t_opt'),
0264         t_opt = select.t_opt;
0265     else
0266         t_opt = true;
0267     end
0268     %
0269     % WARNING : For now, TCV HXRS data are count rate within energy bins, while TS HXR is
0270     % integrated photon count above energy threshold. Quid of TCV HXRC?
0271     %
0272     if t_opt,
0273         if isfield(equil,'info') && isfield(equil.info,'t1') && isfield(equil.info,'t2'),
0274             t1 = equil.info.t1;
0275             t2 = equil.info.t2;
0276         elseif isfield(external,'id') && length(find(external.id == '_')) >= 5,
0277             idsep = find(external.id == '_');
0278             t1 = str2num(external.id(idsep(3)+1:idsep(4)-1));
0279             t2 = str2num(external.id(idsep(4)+1:idsep(5)-1));
0280         else
0281             t1 = NaN;
0282             t2 = NaN;
0283         end
0284         %
0285         t1 = input_dke_yp('Initial time for TCV HXR data averaging',t1,[0;Inf],'',[1,1]);
0286         t2 = input_dke_yp('Final time for TCV HXR data averaging',t2,[0;Inf],'',[1,1]);
0287         %
0288         it1 = find(t >= t1,1,'first');
0289         it2 = find(t <= t2,1,'last');
0290         %
0291         % time averaging
0292         %
0293         hxrexp = squeeze(sum(hxr(it1:it2,:,:)));
0294         dt = t(it2) - t(it1);
0295         hxrexp = hxrexp/dt;%taux de comptage
0296         %
0297         % incremental contribution of energy canal
0298         %
0299         hxrexp = hxrexp - [hxrexp(2:end,:);zeros(1,size(hxrexp,2))];% raw data is all hits above given energy
0300         %
0301     else
0302         %
0303         CRmin = min(hxr(hxr(:) > 0));
0304         %
0305         dt = 1/CRmin;% assuming a single photon count occured somewhere... to improve from database!
0306         %
0307         hxrexp = hxr;
0308     end
0309     %
0310 end        
0311 %
0312 if all(isfield(output,{'equilHXR','radialHXR','Zbremplasma'})),% iluke - local bremsstrahlung emissivity already calculated
0313     equilHXR = output.equilHXR;
0314     radialHXR = output.radialHXR;
0315     Zbremplasma = output.Zbremplasma;
0316 else
0317     %
0318     % refined LUKE equilibrium for FEB calculations
0319     %
0320     radialHXR.xpsin_S_dke = radialDKE.xpsin_S_dke;% warning : check in DKE mode
0321     radialHXR.xpsin_f = radialDKE.xpsin_f_dke;% select only FP points
0322     %
0323     [equilHXR] = equilibrium_jd(equil,radialHXR,dkedisplay,'spline',hxrparam.mfactor);%Magnetic equilibrium for HXR at points of the distribution function
0324     %
0325     % FEB emissivity on LUKE equilibrium
0326     %
0327     if iscell(dke_out.XXf0),% select only FP points
0328         XXf0 = dke_out.XXf0{end}(:,:,radialDKE.r_dke);
0329     else % older LUKE versions
0330         XXf0 = dke_out.XXf0(:,:,radialDKE.r_dke);
0331     end    
0332     %
0333     [Zbremplasma] = bremsstrahlung_dke_yp(equilHXR,mksa,momentumDKE,Zmomcoef,XXf0,hxrparam,dkedisplay);
0334     %
0335     if all(isfield(external,{'param','cam','exp'})),% iluke - only local bremsstrahlung emissivity requested
0336         varargout{1} = equilHXR;
0337         varargout{2} = radialHXR;
0338         varargout{3} = Zbremplasma;
0339         %
0340         return
0341     end
0342     %
0343 end
0344 %
0345 % HXR diagnostic lines of sight
0346 %
0347 [Zbremchord] = bremchord_dke_yp(equilHXR,radialHXR,hxrspec,hxrparam);
0348 %
0349 % FEB emission along the lines of sight
0350 %
0351 [Zbremplasma] = bremchordemission_dke_jd(Zbremplasma,Zbremchord,hxrspec);
0352 %
0353 % HXR synthetic diagnostic signal
0354 %
0355 [Zbremdiag] = bremdiag_dke_yp(Zbremplasma,hxrspec,hxrparam,dt);
0356 %
0357 %For display HXR
0358 %
0359 %
0360 if isfield(select,'d_opt'),
0361     d_opt = select.d_opt;
0362 else
0363     d_opt = input_dke_yp('Do you want to plot the results? [y/n]','n',{'y','n'},'',[1,1]);
0364 end
0365 %
0366 if strcmp(d_opt,'y'),
0367     %
0368     c_opt = input(['Focus on chord? (1-',num2str(nc),') [1] : ']);
0369     if isempty(c_opt),
0370         c_opt = 1;
0371     end
0372     if ~isnumeric(c_opt) | ~any(c_opt == 1:nc),
0373         disp('-----> Invalid chord number');
0374         disp(quitmsg)
0375         return
0376     end
0377     %
0378     r_opt = input('Local bremsstrahlung at radial position? (0<rho<1) [0.5] : ');
0379     if isempty(r_opt),
0380         r_opt = 0.5;
0381     end
0382     if ~isnumeric(r_opt) | r_opt < 0 | r_opt > 1,
0383         disp('-----> Invalid radial position');
0384         disp(quitmsg)
0385         return
0386     end
0387     %
0388     t_opt = input('Local bremsstrahlung at poloidal position? (0<theta<360) [0] : ');
0389     if isempty(t_opt),
0390         t_opt = eps;
0391     end
0392     if ~isnumeric(t_opt) || t_opt < 0 || t_opt > 360,
0393         disp('-----> Invalid poloidal position');
0394         disp(quitmsg)
0395         return
0396     end
0397     %
0398     k_opt = input(['Local bremsstrahlung at energy? (',num2str(Zbremdiag.kfit),') [',num2str(Zbremdiag.kfit(1)),'] : ']);
0399     if isempty(k_opt),
0400         k_opt = Zbremdiag.kfit(1);
0401     end
0402     if ~isnumeric(t_opt) || ~any(k_opt == Zbremdiag.kfit),
0403         disp('-----> Invalid energy');
0404         disp(quitmsg)
0405         return
0406     else
0407         ik = find(k_opt == Zbremdiag.kfit);
0408     end
0409     %
0410     bremchorddisplay_dke_yp(Zbremchord,c_opt,[],p_opt,equil.id);
0411     %
0412     if ~isempty(r_opt) && ~isempty(t_opt) && ~isempty(k_opt)
0413         bremlocaldisplay_dke_yp(dkeparam,dkedisplay,equilDKE,momentumDKE,Zmomcoef,radialDKE,dke_out,mksa,Zbremchord,Zbremplasma,r_opt,t_opt,k_opt,p_opt,equil.id)
0414     end
0415     %
0416     if any(isnan(energy)),
0417         bremdiagdisplay_dke_yp(hxrspec,Zbremchord,Zbremplasma,Zbremdiag,c_opt,k_opt,p_opt,equil.id)
0418     else
0419         bremdiagdisplay_dke_yp(hxrspec,Zbremchord,Zbremplasma,Zbremdiag,c_opt,k_opt,p_opt,equil.id,hxrexp(ik,1:nc))
0420     end
0421 end
0422 %
0423 if isfield(select,'s_opt'),
0424     s_opt = select.s_opt;
0425 else
0426     s_opt = input_dke_yp('Do you want to save the HXR output data? [y/n]','n',{'y','n'},'',[1,1]);
0427 end
0428 %
0429 if strcmp(s_opt,'y'),
0430     sd = pwd;
0431     cd(workdir);
0432     [save_file,save_path] = uiputfile(['HXR_RESULTS_',id_simul,'*.mat'],'Please provide the RESULTS file name',['HXR_RESULTS_',id_simul,'.mat']);
0433     %
0434     if save_file == 0,
0435         disp('-----> No file name specified, HXR output data not saved.')
0436     else
0437         save([save_path,save_file],'Zbremchord','Zbremplasma','Zbremdiag','equilHXR','hxrexp');
0438         disp('-----> The output data is saved')
0439     end
0440     cd(sd);
0441 end 
0442 %
0443 if all(isfield(external,{'param','cam','exp'})),% iluke - only local bremsstrahlung emissivity requested
0444     varargout{1} = Zbremchord;
0445     varargout{2} = Zbremplasma;
0446     varargout{3} = Zbremdiag;
0447 end
0448 
0449 
0450