runhxr_jd

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

0001 function hxrdata = runhxr_jd(data,basestr,hxr_file,exp_file,workdir,texp,opt_save)
0002 %
0003 % This function calculates the hard X-ray emisssion from LUKE results
0004 %
0005 if nargin < 7,
0006     opt_save = 0;
0007 end
0008 if nargin < 6,
0009     texp = NaN;
0010 end
0011 if nargin < 5,
0012     workdir = '';
0013 end
0014 %
0015 if isstruct(data),
0016     %
0017     % data from input structure
0018     %
0019     if isfield(data,'lukestructs'),
0020         data = data.lukestructs;
0021     end
0022     %
0023     if isfield(data,'luke_input'),
0024         equil = data.luke_input.equil;
0025         dkeparam = data.luke_input.dkeparam;
0026         dkedisplay = data.luke_input.dkedisplay;
0027     else
0028         equil = data.equil;
0029         dkeparam = data.dkeparam;
0030         dkedisplay = data.dkedisplay;
0031     end
0032     %
0033     if isfield(data,'output'),% case of LUKE_DATA file
0034         Zcurr = data.output.Zcurr;
0035         dke_out = data.output.dke_out;
0036         equilDKE = data.output.equilDKE;
0037         radialDKE = data.output.radialDKE;
0038         momentumDKE = data.output.momentumDKE;
0039         mksa = data.output.mksa;
0040         Zmomcoef = data.output.Zmomcoef;
0041         if isfield(data,'external'),
0042             external = data.external;
0043         else
0044             external = '';
0045         end
0046     else % case of LUKE_RESULTS file
0047         Zcurr = data.Zcurr;
0048         dke_out = data.dke_out;
0049         equilDKE = data.equilDKE;
0050         radialDKE = data.radialDKE;
0051         momentumDKE = data.momentumDKE;
0052         mksa = data.mksa;
0053         Zmomcoef = data.Zmomcoef;
0054         external = '';
0055     end        
0056     %
0057 else
0058     %
0059     % data from file
0060     %
0061     external = '';
0062     %
0063     load(data,'equil','dkeparam','dkedisplay','Zcurr','dke_out','equilDKE','radialDKE','momentumDKE','mksa','Zmomcoef','external');
0064     %
0065 end    
0066 %
0067 clear data
0068 %
0069 if ~isfield(Zmomcoef,'Ec_ref'),
0070     if isfield(mksa,'Ec_ref'),
0071         Zmomcoef.Ec_ref = mksa.Ec_ref;
0072     else
0073         error('Kinetic energy data missing. Please correct manually.')
0074     end
0075 end    
0076 %
0077 % load HXR camera spec
0078 %
0079 load(hxr_file,'hxrspec','hxrparam');
0080 %
0081 if ~exist('hxrspec','var'),%old spec file
0082     %
0083     load(hxr_file,'hxr','hxrparam');
0084     hxrspec = hxr;
0085     clear hxr %to avoid name conflict
0086     %
0087 end
0088 %
0089 nc = length(hxrspec.beta_hxr);% number of chords
0090 %
0091 % load and treat HXR experimental results
0092 %
0093 if all(ischar(exp_file)),
0094     %
0095     if isfield(external,'hxr'),
0096         %
0097         t = external.hxr.t;
0098         energy = external.hxr.energy;
0099         hxr = external.hxr.counts;
0100         %filter_thickness = external.hxr.filter_thickness;
0101         %filter_material = external.hxr.filter_material;
0102         %
0103         if isfield(external,'beta')
0104             beta = external.beta;
0105         else
0106             beta = NaN(1,size(hxr,3));
0107         end
0108         %
0109         if iscell(t),%old TCV data file
0110             t = t{ishot};
0111             energy = energy{ishot};
0112             hxr = hxr{ishot};
0113             %filter_thickness = filter_thickness{ishot};
0114             %filter_material = filter_material{ishot};
0115         end
0116         %
0117     else
0118         %
0119         % load HXR experimental results
0120         %
0121         if isfield(equil,'info') && isfield(equil.info,'shotnum'),
0122             shotnum = equil.info.shotnum;
0123         elseif isfield(equil,'id'),
0124             idsep = find(equil.id == '_');
0125             shotnum = equil.id(idsep(1)+1:idsep(2)-1); 
0126         else
0127             disp('-----> Shot number information missing');
0128             return
0129         end
0130         %
0131         if ~isempty(exp_file),% hxr experimental data from file
0132             %
0133             load(exp_file,'shot','energy','hxr','t','beta');%,'filter_thickness','filter_material'
0134             %
0135             ishot = find(shot == str2num(shotnum));
0136             %
0137             if isempty(ishot) || ~exist('hxr','var'),
0138                 disp(['-----> hxr data not found for shot :',shotnum]);
0139                 disp(quitmsg)
0140                 return
0141             end
0142             %
0143             if length(ishot) > 1,
0144                 disp(['-----> WARNING : Several entries are found for shot :',shotnum,', first entry selected']);
0145                 ishot = ishot(1);
0146             end
0147             %
0148             clear exp_file
0149             %
0150             if iscell(t),%old TCV data file
0151                 t = t{ishot};
0152                 energy = energy{ishot};
0153                 hxr = hxr{ishot};
0154                 %filter_thickness = filter_thickness{ishot};
0155                 %filter_material = filter_material{ishot};
0156             end
0157             %
0158             clear ishot
0159             %
0160         else % hxr experimental data fromdatabase
0161             %
0162             [energy,hxr,t,beta] = collect_hxr_forLUKE_jd(basestr,str2num(shotnum));%,opt,data_path,'filter_thickness','filter_material'
0163             %
0164         end
0165         %
0166         clear shotnum
0167         %
0168         if ~exist('beta','var') || (isscalar(beta) && isnan(beta)),
0169             beta = NaN(1,size(hxr,3));
0170         end
0171         %
0172     end
0173     %
0174     % test angles
0175     %
0176     beta(beta > pi) = 2*pi - beta(beta > pi);% correct beta definition
0177     %
0178     if length(beta) ~= nc,
0179         error('number of chords mismatch')
0180     end
0181     if ~all(isnan(beta)) && any(abs(beta - hxrspec.beta_hxr) > 1e-15),
0182         error('angle of chords mismatch')
0183     end
0184     %
0185     energy = energy.';
0186     %
0187     % time sample of experimental data
0188     %
0189     if isnan(texp),
0190         if isfield(equil,'info') && isfield(equil.info,'t1') && isfield(equil.info,'t2'),
0191             t1 = equil.info.t1;
0192             t2 = equil.info.t2;
0193         elseif isfield(equil,'external') && isfield(equil.external,'id'),
0194             idsep = find(equil.external.id == '_');
0195             t1 = str2num(equil.external.id(idsep(3)+1:idsep(4)-1));
0196             t2 = str2num(equil.external.id(idsep(4)+1:idsep(5)-1));
0197         else
0198             disp(['-----> hxr averaging time span not found']);
0199             disp(quitmsg)
0200             return
0201         end
0202     else
0203         t1 = texp(1);
0204         t2 = texp(2);
0205     end
0206     %
0207     it1 = find(t >= t1,1,'first');
0208     it2 = find(t <= t2,1,'last');
0209     %
0210     if isempty(it1) || isempty(it2),% no data available in averaging range
0211         %
0212         hxrexp = NaN*ones(length(energy),nc);
0213         dt = 1;
0214         %
0215     else        
0216         %
0217         % time averaging
0218         %
0219         hxrexp = squeeze(sum(hxr(it1:it2,:,:)));
0220         dt = t(it2) - t(it1);
0221         hxrexp = hxrexp/dt;%taux de comptage
0222         %
0223         % incremental contribution of energy canal
0224         %
0225         hxrexp = hxrexp - [hxrexp(2:end,:);zeros(1,size(hxrexp,2))];% raw data is all hits above given energy
0226         %
0227         hxrexp(:,any(hxrexp < 0)) = NaN;% hxrexp cannot be negative - probably malfunctioning chords
0228         %
0229     end
0230     %
0231     hxrspec.kdiag_hxr = [energy(1:end-1),hxrspec.kdiag_hxr(1,end);energy(2:end),hxrspec.kdiag_hxr(2,end)];%redefine kdiag according to measurements
0232     %
0233 else
0234     %
0235     filter_thickness = NaN;
0236     filter_material = NaN;
0237     %
0238     energy = [hxrspec.kdiag_hxr(1,1:end-1),hxrspec.kdiag_hxr(2,end-1)];
0239     hxrexp = NaN*ones(length(energy),nc);
0240     dt = 1;
0241     %
0242 end
0243 %
0244 [Zbremchord,equilHXR] = bremchord_dke_yp(dkeparam,dkedisplay,equil,radialDKE,Zcurr,hxrspec,hxrparam);
0245 %
0246 Zbremchord.mask_hxr(:) = 1;%so that the bremmstrahlung is calculated at all positions, for display purposes
0247 %
0248 [Zbremplasma] = bremsstrahlung_dke_yp(dkeparam,dkedisplay,equilHXR,radialDKE,mksa,momentumDKE,Zmomcoef,dke_out,Zbremchord,hxrspec,hxrparam);
0249 %
0250 [Zbremdiag] = bremdiag_dke_yp(Zbremplasma,hxrspec,hxrparam,dt);
0251 %
0252 % Save results
0253 %
0254 if opt_save,
0255     id_simul = dke_out.simul.id;
0256     filename = [workdir,'HXR_RESULTS_',id_simul,'.mat'];
0257     save(filename,'Zbremchord','Zbremplasma','Zbremdiag','equilHXR','hxrexp');
0258     disp(['HXR data saved in file ',filename])
0259 end
0260 %
0261 hxrdata.Zbremchord = Zbremchord;
0262 hxrdata.Zbremplasma = Zbremplasma;
0263 hxrdata.Zbremdiag = Zbremdiag;
0264 hxrdata.equilHXR = equilHXR;
0265 hxrdata.hxrexp = hxrexp;
0266 %
0267 
0268 
0269 
0270