iluke_hxrexpdisp

       data: [40x19x110 double]
       time: [1x40 double]
     energy: [15 25 35 45 55 65 75 85 95 105 115 125 135 145 155 165 175 185 195]

0001 function [nobj,graphids] = iluke_hxrexpdisp(hxr,sdiags,shotimes,ic,ichord,ibin,axs,style)
0002 %
0003 %       data: [40x19x110 double]
0004 %       time: [1x40 double]
0005 %     energy: [15 25 35 45 55 65 75 85 95 105 115 125 135 145 155 165 175 185 195]
0006 %
0007 camstr = ['Cam #',num2str(hxr.cam{ic}.camnum)];
0008 %
0009 chordsind = hxr.cam{ic}.chordsind;% indices for chords of camera # ic
0010 %
0011 nshotimes = length(shotimes);
0012 %
0013 if iscell(sdiags) && isfield(sdiags{1},'Zbremdiag'),% HXR calculations by LUKE performed,
0014     %
0015     lukedata = NaN(nshotimes,length(sdiags{1}.Zbremdiag.kfit),length(chordsind));
0016     %
0017     for ishotime = 1:nshotimes,
0018         lukedata(ishotime,:,:) = reshape(sdiags{ishotime}.Zbremdiag.brem_bd_diag(chordsind,:).',[1,length(sdiags{1}.Zbremdiag.kfit),length(chordsind)]);
0019     end
0020     flag_luke = true;
0021 else
0022     flag_luke = false;
0023 end
0024 %
0025 if isfield(hxr,'exp') && ~isempty(hxr.exp),% experimental data exist
0026     hxrdata = hxr.exp.data(:,:,chordsind);
0027     flag_exp = true;
0028     energy = hxr.exp.energy;
0029 else
0030     hxrdata = NaN(size(lukedata));
0031     flag_exp = false;
0032     energy = sdiags{1}.Zbremdiag.kfit;
0033 end
0034 %
0035 if ~flag_luke,
0036     lukedata = NaN(size(hxrdata));
0037 end
0038 %
0039 if flag_luke && flag_exp,
0040     leg = {'Exp.','LUKE'};
0041 else
0042     leg = NaN;
0043 end
0044 %
0045 if ichord > 0,% display for one chord
0046     hxrdata = hxrdata(:,:,ichord);
0047     lukedata = lukedata(:,:,ichord);
0048     chordstr = ['Chord #',num2str(hxr.cam{ic}.chordsnum(ichord))];
0049     nchords = 1;
0050 else
0051     chordsnum = hxr.cam{ic}.chordsnum;
0052     nchords = length(chordsnum);
0053 end
0054 %
0055 if ibin > 0,% display for one energy
0056     nbins = 1;
0057     hxrdata = hxrdata(:,ibin,:);
0058     lukedata = lukedata(:,ibin,:);
0059     if ibin < length(energy),
0060         energystr = [num2str(energy(ibin)),' < k (keV) < ',num2str(energy(ibin + 1))];
0061     else    
0062         energystr = ['k (keV) > ',num2str(energy(ibin))];
0063     end
0064 else
0065     energyind = find(energy(1:end-1) >= hxr.param.kmin & energy(2:end) <= hxr.param.kmax);
0066     %
0067     energysep = [3*energy(1) - energy(2),energy(2:end) + energy(1:end-1),3*energy(end) - energy(end-1)]/2;
0068     denergy = diff(energysep);
0069     nbins = length(energy);
0070     %
0071     energystr = [num2str(hxr.param.kmin),' < k (keV) < ',num2str(hxr.param.kmax)];
0072 end
0073 %
0074 nobj = [0,0,0];
0075 graphids = {'','',''};
0076 %
0077 if nshotimes > 1,% display time dependence
0078     %
0079     timesind = [];
0080     shotimesmask = false(1,nshotimes);
0081     timesep = [3*hxr.exp.time(1) - hxr.exp.time(2),hxr.exp.time(2:end) + hxr.exp.time(1:end-1),3*hxr.exp.time(end) - hxr.exp.time(end-1)]/2;
0082     dt = diff(timesep);
0083     xlim = [timesep(1),timesep(end)];
0084     %
0085     for ishotime = 1:nshotimes,
0086         it = find(hxr.exp.time == shotimes(ishotime));
0087         %
0088         if length(it) == 1,
0089             timesind = [timesind,it];
0090             shotimesmask(ishotime) = true;
0091         end
0092     end
0093     %
0094     shotimes = shotimes(shotimesmask);
0095     %
0096     hxrdata = hxrdata(timesind,:,:);
0097     dt = dt(timesind);
0098     %
0099     nt = length(timesind);
0100     %
0101     if ichord == 0,% display for all chords
0102         %
0103         % TODO
0104         %
0105     else% display for one chord
0106         %
0107         if ibin == 0,% display for all energies - show global count rate and photon temperature
0108             %
0109             Nphotons = repmat(dt,[nbins,1]).*hxrdata.';% number of photons measured per energy bin
0110             hxrsum = sum(hxrdata(:,energyind),2).';
0111             ehxrsum = sqrt(hxrsum./dt);
0112             %
0113             [Tphexp,eTphexp] = tphfit_dke_yp(repmat(energy(energyind).',[1,nt]),Nphotons(energyind,:),sqrt(Nphotons(energyind,:)));
0114             Tphmax = 2*Tphexp(find(hxrsum == max(hxrsum),1,'first'));
0115             if isempty(Tphmax),
0116                 Tphmax = NaN;
0117             end
0118             Tphmask = Tphexp > 0 & Tphexp < Tphmax & eTphexp < Tphexp;
0119             Tphexp(~Tphmask) = NaN;
0120             eTphexp(~Tphmask) = NaN;
0121             %
0122             Nphotons_luke = repmat(dt,[nbins,1]).*lukedata.';% number of photons measured per energy bin
0123             hxrsum_luke = sum(lukedata(:,energyind),2).';
0124             ehxrsum_luke = sqrt(hxrsum_luke./dt);
0125             %
0126             [Tphluke,eTphluke] = tphfit_dke_yp(repmat(energy(energyind).',[1,nt]),Nphotons_luke(energyind,:),sqrt(Nphotons_luke(energyind,:)));
0127             Tphmax_luke = 2*Tphluke(find(hxrsum_luke == max(hxrsum_luke),1,'first'));
0128             if isempty(Tphmax_luke),
0129                 Tphmax_luke = NaN;
0130             end
0131             Tphmask_luke = Tphluke > 0 & Tphluke < Tphmax_luke & eTphluke < Tphluke;
0132             Tphluke(~Tphmask_luke) = NaN;
0133             eTphluke(~Tphmask_luke) = NaN;
0134             %
0135             graph1D_jd(shotimes,Tphexp,0,0,'t (s)','T_{ph.} (keV)','',NaN,xlim,'','-','o','r',2,style,axs(1),1,NaN,NaN,eTphexp);%Poloidal plot
0136             graph1D_jd(shotimes,Tphluke,0,0,'','','',leg,xlim,'','-','s','b',2,style,axs(1),1,NaN,NaN,eTphluke);%Poloidal plot
0137             ylim = get(axs(1),'ylim');
0138             set(axs(1),'ylim',[0,min([ylim(2),Tphmax,Tphmax_luke])]);drawnow
0139             %
0140             nobj(1) = 1 + iscell(leg);
0141             graphids{1} = ['Tph - ',camstr,'; ',chordstr,'; ',energystr];  
0142             %
0143             graph1D_jd(shotimes,hxrsum,0,0,'t (s)','CR (Counts/s)','',NaN,xlim,'','-','o','r',2,style,axs(2),1,NaN,NaN,ehxrsum);%Poloidal plot
0144             graph1D_jd(shotimes,hxrsum_luke,0,0,'','','',leg,xlim,'','-','s','b',2,style,axs(2),1,NaN,NaN,ehxrsum_luke);%Poloidal plot
0145             %
0146             nobj(2) = 1 + iscell(leg);
0147             graphids{2} = ['CR - ',camstr,'; ',chordstr,'; ',energystr];                  
0148             %
0149         else% display for one energy - display specific count rate
0150             %
0151             ehxrdata = sqrt(hxrdata./dt.');
0152             elukedata = sqrt(lukedata./dt.');
0153             %
0154             graph1D_jd(shotimes,hxrdata,0,0,'t (s)','CR (Counts/s)','',NaN,xlim,'','-','o','r',2,style,axs(1),1,NaN,NaN,ehxrdata);%Poloidal plot
0155             graph1D_jd(shotimes,lukedata,0,0,'','','',leg,xlim,'','-','s','b',2,style,axs(1),1,NaN,NaN,elukedata);%Poloidal plot
0156             %
0157             nobj(1) = 1 + iscell(leg);
0158             graphids{1} = ['CR - ',camstr,'; ',chordstr,'; ',energystr];  
0159             %
0160         end
0161         %
0162     end
0163     %
0164 else% display for one time step
0165     %
0166     if flag_exp,
0167         it = find(hxr.exp.time == shotimes);
0168         if isempty(it),
0169             return
0170         end
0171     else
0172         it = 1;
0173     end
0174     %
0175     timestr = ['t = ',num2str(shotimes,'%6.4f')];
0176     %
0177     hxrdata = reshape(hxrdata(it,:,:),[nbins,nchords]);
0178     if flag_luke,
0179         lukedata = reshape(lukedata,[nbins,nchords]);
0180     else
0181         lukedata = reshape(lukedata(it,:,:),[nbins,nchords]);
0182     end
0183     %
0184     if flag_exp,
0185         CRmin = min(hxrdata(hxrdata(:) > 0));
0186         %
0187         dt = 1/CRmin;% assuming a single photon count occured somewhere... to improve from database!
0188     else
0189         dt = sdiags{1}.Zbremdiag.dt;
0190     end
0191     %
0192     Nphotons = dt*hxrdata;% number of photons measured per energy bin
0193     Nphotons_luke = dt*lukedata;% number of photons measured per energy bin
0194     %
0195     if ichord == 0,% display for all chords
0196         %
0197         xlim = [chordsnum(1) - 1,chordsnum(end)+1];
0198         %
0199         if ibin == 0,% display for all energies - show global count rate and photon temperature
0200             %
0201             hxrsum = sum(hxrdata(energyind,:),1);
0202             ehxrsum = sqrt(hxrsum/dt);
0203             %
0204             [Tphexp,eTphexp] = tphfit_dke_yp(repmat(energy(energyind).',[1,nchords]),Nphotons(energyind,:),sqrt(Nphotons(energyind,:)));
0205             Tphmax = 2*Tphexp(find(hxrsum == max(hxrsum),1,'first'));
0206             if isempty(Tphmax),
0207                 Tphmax = NaN;
0208             end
0209             Tphmask = Tphexp > 0 & Tphexp < Tphmax & eTphexp < Tphexp;
0210             Tphexp(~Tphmask) = NaN;
0211             eTphexp(~Tphmask) = NaN;
0212             %
0213             hxrsum_luke = sum(lukedata(energyind,:),1);
0214             ehxrsum_luke = sqrt(hxrsum_luke/dt);
0215             %
0216             [Tphluke,eTphluke] = tphfit_dke_yp(repmat(energy(energyind).',[1,nchords]),Nphotons_luke(energyind,:),sqrt(Nphotons_luke(energyind,:)));
0217             Tphmax_luke = 2*Tphluke(find(hxrsum_luke == max(hxrsum_luke),1,'first'));
0218             if isempty(Tphmax_luke),
0219                 Tphmax_luke = NaN;
0220             end
0221             Tphmask_luke = Tphluke > 0 & Tphluke < Tphmax_luke & eTphluke < Tphluke;
0222             Tphluke(~Tphmask_luke) = NaN;
0223             eTphluke(~Tphmask_luke) = NaN;
0224             %
0225             graph1D_jd(chordsnum,Tphexp,0,0,'chord #','T_{ph.} (keV)','',NaN,xlim,'','-','o','r',2,style,axs(1),1,NaN,NaN,eTphexp);%Poloidal plot
0226             graph1D_jd(chordsnum,Tphluke,0,0,'','','',leg,xlim,'','-','s','b',2,style,axs(1),1,NaN,NaN,eTphluke);%Poloidal plot
0227             ylim = get(axs(1),'ylim');
0228             set(axs(1),'ylim',[0,min([ylim(2),Tphmax,Tphmax_luke])]);drawnow
0229             %
0230             nobj(1) = 1 + iscell(leg);
0231             graphids{1} = ['Tph - ',camstr,'; ',timestr,'; ',energystr];  
0232             %
0233             graph1D_jd(chordsnum,hxrsum,0,0,'chord #','CR (Counts/s)','',NaN,xlim,'','-','o','r',2,style,axs(2),1,NaN,NaN,ehxrsum);%Poloidal plot
0234             graph1D_jd(chordsnum,hxrsum_luke,0,0,'','','',leg,xlim,'','-','s','b',2,style,axs(2),1,NaN,NaN,ehxrsum_luke);%Poloidal plot
0235             %
0236             nobj(2) = 1 + iscell(leg);
0237             graphids{2} = ['CR - ',camstr,'; ',timestr,'; ',energystr];                  
0238             %
0239         else% display for one energy - display specific count rate
0240             %
0241             ehxrdata = sqrt(hxrdata/dt);
0242             elukedata = sqrt(lukedata/dt);
0243             %
0244             graph1D_jd(chordsnum,hxrdata,0,0,'chord #','CR (Counts/s)','',NaN,xlim,'','-','o','r',2,style,axs(1),1,NaN,NaN,ehxrdata);%Poloidal plot
0245             graph1D_jd(chordsnum,lukedata,0,0,'','','',leg,xlim,'','-','s','b',2,style,axs(1),1,NaN,NaN,elukedata);%Poloidal plot
0246             %
0247             nobj(1) = 1 + iscell(leg);
0248             graphids{1} = ['CR - ',camstr,'; ',timestr,'; ',energystr];  
0249             %
0250         end
0251         %
0252     else% display spectrum
0253         %
0254         if ibin == 0,% display for all energies - show global count rate and photon temperature
0255             %
0256             xlim = [0,max(energysep)];
0257             %
0258             [Tphexp,~,Aexp] = tphfit_dke_yp(energy(energyind).',Nphotons(energyind),sqrt(Nphotons(energyind)));
0259             %
0260             [Tphluke,~,Aluke] = tphfit_dke_yp(energy(energyind).',Nphotons_luke(energyind),sqrt(Nphotons_luke(energyind)));
0261             %
0262             %log(Ei*Ni/dE) = a+b*Ei -> Tph = -1/b
0263             %
0264             ehxrdata = sqrt(hxrdata/dt);
0265             hxrdatafit = Aexp*exp(-energy/Tphexp).*denergy./energy/dt;
0266             %
0267             elukedata = sqrt(lukedata/dt);
0268             lukedatafit = Aluke*exp(-energy/Tphluke).*denergy./energy/dt;
0269             %
0270             graph1D_jd(energy,hxrdata,0,1,'k (keV)','CR (Counts/s)','',NaN,xlim,'','none','o','r',2,style,axs(1),1,NaN,NaN,ehxrdata);%Poloidal plot
0271             graph1D_jd(energy,lukedata,0,1,'','','',leg,xlim,'','none','s','b',2,style,axs(1),1,NaN,NaN,elukedata);%Poloidal plot
0272             graph1D_jd(energy,hxrdatafit,0,1,'','','',NaN,xlim,'','--','none','r',1,style,axs(1));%Poloidal plot
0273             graph1D_jd(energy,lukedatafit,0,1,'','','',NaN,xlim,'','--','none','b',1,style,axs(1));%Poloidal plot
0274             graph1D_jd(hxr.param.kmin*[1,1],get(axs(1),'ylim'),0,1,'','','',NaN,xlim,'','--','none','k',2,style,axs(1));%Poloidal plot
0275             graph1D_jd(hxr.param.kmax*[1,1],get(axs(1),'ylim'),0,1,'','','',NaN,xlim,'','--','none','k',2,style,axs(1));%Poloidal plot
0276             %
0277             nobj(1) = 1 + iscell(leg);
0278             graphids{1} = ['dCRdE - ',camstr,'; ',timestr,'; ',chordstr];  
0279             %
0280         else% display for one energy - display specific count rate
0281             %
0282             % TODO
0283             %
0284         end
0285         %
0286     end
0287     %
0288 end
0289 %