luke_get_trdata

LUKE_GET_TRDATA d = luke_get_trdata( a )
   For an input structure a this function returns the output data
   structure d containing time resolved data from TCV for LUKE.

 Same input for all modes (other inputs optional):
   a.shot
   a.diag{i}.name
   a.diag{j}.isel for a.diag{j}.name=='ECRH'

 Input for mode with automatic selection of (all) Thomson points:
   do not set a.t -> all Thomson points
   a.t set, but not a.dtres -> use a.t with Thomson resolution

 Input for mode with fixed time points and limited averaging range:
   put time points in a.t and set a.dtres=[]
   set a.dtreslim to specify averaging range limit (default: Thomson time)

 Input for mode with equidistant resolution in time interval:
   set a.t as a time interval
   specify the resolution in a.dtres


   J. Kamleitner, CRPP, EPFL, Feb 2014

0001 function d = luke_get_trdata( a )
0002 %LUKE_GET_TRDATA d = luke_get_trdata( a )
0003 %   For an input structure a this function returns the output data
0004 %   structure d containing time resolved data from TCV for LUKE.
0005 %
0006 % Same input for all modes (other inputs optional):
0007 %   a.shot
0008 %   a.diag{i}.name
0009 %   a.diag{j}.isel for a.diag{j}.name=='ECRH'
0010 %
0011 % Input for mode with automatic selection of (all) Thomson points:
0012 %   do not set a.t -> all Thomson points
0013 %   a.t set, but not a.dtres -> use a.t with Thomson resolution
0014 %
0015 % Input for mode with fixed time points and limited averaging range:
0016 %   put time points in a.t and set a.dtres=[]
0017 %   set a.dtreslim to specify averaging range limit (default: Thomson time)
0018 %
0019 % Input for mode with equidistant resolution in time interval:
0020 %   set a.t as a time interval
0021 %   specify the resolution in a.dtres
0022 %
0023 %
0024 %   J. Kamleitner, CRPP, EPFL, Feb 2014
0025 
0026 % first version:
0027 %   * without conditional averaging
0028 %   * with preselected diagnostics (depending on availability and parameters shot, t, dtres)
0029 
0030 % required data can be seen in LUKE function load_externaldata_TCV.m
0031 % best to use EQUIL as first diagnostic since others are based on it.
0032 
0033 % dependencies:
0034 % CRPP mds functions, tcvget, TCV database access and liuqe (latter on lac5)
0035 % gti_chdtres.m
0036 % gti_filter.m
0037 % make_eqdsk.m
0038 
0039 % some constants
0040 if isfield(a,'flagTPG')
0041     flagTPG=a.flagTPG;   % true to use TPG ox_fraction function
0042 else
0043     flagTPG=true;   % true to use TPG ox_fraction function
0044 end
0045 if isfield(a,'flagFF')
0046     flagFF=a.flagFF;    % true to use FF ox_fraction function
0047 else
0048     flagFF=true;    % true to use FF ox_fraction function
0049 end
0050 % dtres0_pgyro=1e-4;
0051 try
0052     outputpath_default=sprintf('/tmp/LUKE_%s_%s',getenv('USER'),timeid_jd(clock));
0053 catch
0054     outputpath_default=sprintf('/tmp/LUKE_%s',getenv('USER'));
0055 end
0056 
0057 % keyboard;
0058 
0059 % common time base if not defined
0060 if(~isfield(a,'t') || ~isfield(a,'dtres'))
0061     mdsopen(a.shot);
0062     tthom=mdsvalue('\results::thomson:times');
0063     dtres_default=round(60*min(diff(tthom)))/60;
0064     dtreslim_default=dtres_default;
0065     [tip,ip]=tcvget('IP');
0066     mdsclose;
0067     itip50=find(abs(ip)>=max(abs(ip))/2);
0068     t_default=tip(itip50([1 end]));
0069     t_default=t_default(:)';    % row vector
0070 elseif(~isfield(a,'dtreslim') || isempty(a.dtreslim))
0071     mdsopen(a.shot);
0072     tthom=mdsvalue('\results::thomson:times');
0073     dtreslim_default=round(60*min(diff(tthom)))/60;
0074     mdsclose;
0075 end
0076 
0077 % go through the list of requested diagnostics
0078 for i=1:length(a.diag)
0079     % set the parameters (priority for parameters in the diag struct)
0080     diag=a.diag{i};
0081     % shot
0082     shot=a.shot;
0083     if(isfield(diag,'shot'))
0084         shot=diag.shot;
0085     end
0086     % t
0087     if(isfield(a,'t'))
0088         t=a.t;
0089     else
0090         t=t_default;
0091     end
0092     if(isfield(diag,'t'))
0093         t=diag.t;
0094     end
0095     % dtres
0096     if(isfield(a,'dtres'))
0097         if(~isfield(a,'t') && isempty(a.dtres))
0098             dtres=dtres_default;
0099         else
0100             dtres=a.dtres;
0101         end
0102     else
0103         dtres=dtres_default;
0104     end
0105     if(isfield(diag,'dtres'))
0106         dtres=diag.dtres;
0107     end
0108     % dtreslim
0109     if(isfield(a,'dtreslim') && ~isempty(a.dtreslim))
0110         dtreslim=a.dtreslim;
0111     else
0112         dtreslim=dtreslim_default;
0113     end
0114     if(isfield(diag,'dtreslim'))
0115         dtreslim=diag.dtreslim;
0116     end
0117     % outputpath
0118     if(isfield(a,'outputpath'))
0119         outputpath=a.outputpath;
0120     else
0121         outputpath=outputpath_default;
0122     end
0123     if(isfield(diag,'outputpath'))
0124         outputpath=diag.outputpath;
0125     end
0126     if(~exist(outputpath,'dir'))
0127         system(sprintf('mkdir -p %s',outputpath));
0128     end
0129     % isel
0130     if(isfield(diag,'isel'))
0131         isel=diag.isel;
0132     end
0133     % tca
0134     if(isfield(a,'tca'))
0135         tca=a.tca;
0136     else
0137         tca=NaN;
0138     end
0139     if(isfield(diag,'tca'))
0140         tca=diag.tca;
0141     end
0142     % tcr
0143     if(isfield(a,'tcr'))
0144         tcr=a.tcr;
0145     else
0146         tcr=NaN;
0147     end
0148     if(isfield(diag,'tcr'))
0149         tcr=diag.tcr;
0150     end
0151     if(~all(isnan(tca)) && ~all(isnan(tcr)))
0152         % conditional averaging on
0153         docrca=true;
0154     else
0155         % no conditional averaging
0156         docrca=false;
0157     end
0158     
0159     % read the data of the diagnostic
0160     fprintf('- Get %s data for shot %d, t=[%g,%g]s.\n',diag.name,shot,t(1),t(end));
0161     switch diag.name
0162         case 'BT'
0163             d.btsign.shot=shot;
0164             d.btsign.rawdiag={'MagneticProbes'};
0165             d.btsign.time=mean(t([1 end]));
0166             mdsopen(shot);
0167             d.btsign.data=sign(mdsdata('\magnetics::iphi[$]',d.btsign.time));
0168             mdsclose;
0169         case 'IP'
0170             d.ipsign.shot=shot;
0171             d.ipsign.rawdiag={'CoilCurrents'};
0172             d.ip.shot=shot;
0173             d.ip.rawdiag={'MagneticProbes'};
0174             d.ipsign.time=mean(t([1 end]));
0175             [Ip,tip,d.ipsign.data]=ip_trapeze_woos(shot);
0176             if(docrca)
0177                 d.ip.atime=tca;
0178                 d.ip.otime=tip;
0179                 d.ip.rtime=tcr;
0180                 d.ip.time=d.ip.rtime;
0181                 d.ip.data=gti_crca(tip,Ip',tca,tcr,dtreslim);
0182             else % if no coherent/conditional averaging:
0183                 [d.ip.time,d.ip.data]=gti_chdtres(tip,Ip',t,dtres,dtreslim);
0184             end
0185         case 'L7'
0186             d.l7.shot=shot;
0187             d.l7.rawdiag={'L7'};
0188             mdsopen(shot);
0189             l7raw=tdi('\base::tr2412_prefl:channel_001');
0190             mdsclose;
0191             if(docrca)
0192                 d.l7.atime=tca;
0193                 d.l7.otime=l7raw.dim{1};
0194                 d.l7.rtime=tcr;
0195                 d.l7.time=d.l7.rtime;
0196                 d.l7.data=gti_crca(l7raw.dim{1},-l7raw.data',tca,tcr,dtreslim);
0197             else % if no coherent/conditional averaging:
0198                 [d.l7.time,d.l7.data]=gti_chdtres(l7raw.dim{1},-l7raw.data',t,dtres,dtreslim);
0199             end
0200         case 'VLOOP'
0201             d.vloop.shot=shot;
0202             d.vloop.rawdiag={'FluxLoops'};
0203             mdsopen(shot);
0204             vloop=tdi('\magnetics::vloop[*,"001"]');
0205             mdsclose;
0206             if(docrca)
0207                 d.vloop.atime=tca;
0208                 d.vloop.otime=vloop.dim{1};
0209                 d.vloop.rtime=tcr;
0210                 d.vloop.time=d.vloop.rtime;
0211                 d.vloop.data=gti_crca(vloop.dim{1},vloop.data',tca,tcr,dtreslim);
0212             else % if no coherent/conditional averaging:
0213                 if(dtres_bt(2e-3,t,dtres,dtreslim))
0214                     % if requested time resolution better than 2ms filter out
0215                     % the strong 1.3kHz noise using a 630Hz low-pass filter
0216                     [tvl,vl]=gti_filter(vloop.dim{1},vloop.data',[-Inf 630]);
0217                     [d.vloop.time,d.vloop.data]=gti_chdtres(tvl,vl,t,dtres,dtreslim);
0218                 else
0219                     % otherwise no filter is applied, just averaging
0220                     [d.vloop.time,d.vloop.data]=gti_chdtres(vloop.dim{1},vloop.data',t,dtres,dtreslim);
0221                 end
0222             end
0223         case 'NE'
0224             d.ne.shot=shot;
0225             if(docrca)
0226                 % take FIR measurements in any case for the moment
0227                 d.ne.rawdiag={'FIR','LIUQE'};
0228                 firpar=gti_gds(shot,t([1 end]),'FIR','Flux1D','MinFishReg','Unique','Unique',0);
0229                 firpar.diag.frequency.filter.hlb='no';
0230                 firpar.diag.frequency.crca.tca=tca;
0231                 firpar.diag.frequency.crca.tcr=tcr;
0232                 firpar.diag.frequency.crca.dtreslim=dtreslim;
0233                 firout=gti_expert(firpar);
0234                 d.ne.source='GTI';
0235                 fprintf('  . Using %s profiles from %s for NE.\n',d.ne.rawdiag{1},d.ne.source);
0236                 d.ne.atime=tca;
0237                 d.ne.otime=firout.diag.otime;
0238                 d.ne.rtime=tcr;
0239                 d.ne.time=firout.res.time;
0240                 d.ne.data=firout.res.inverted;
0241                 d.ne.rho=firout.res.rho;
0242             else % if no coherent/conditional averaging:
0243                 if(dtres_bt(0.016,t,dtres,dtreslim))
0244                     d.ne.rawdiag={'FIR','LIUQE'};
0245                     firpar=gti_gds(shot,t([1 end]),'FIR','Flux1D','MinFishReg','Unique','Unique',0);
0246                     firpar.diag.frequency.filter.hlb='no';
0247                     firout=gti_expert(firpar);
0248                     d.ne.source='GTI';
0249                     fprintf('  . Using %s profiles from %s for NE.\n',d.ne.rawdiag{1},d.ne.source);
0250                     [d.ne.time,d.ne.data]=gti_chdtres(firout.res.time,firout.res.inverted,t,dtres,dtreslim);
0251                     d.ne.rho=firout.res.rho;
0252                 else
0253                     d.ne.rawdiag={'Thomson','LIUQE','FIR'};
0254                     neraw=luke_get_nete_raw(diag.name,shot);
0255                     d.ne.source=neraw.source;
0256                     fprintf('  . Using %s profiles from %s for NE.\n',d.ne.rawdiag{1},d.ne.source);
0257                     jnotnan=~any(isnan(neraw.data));
0258                     if(all(~jnotnan)) % if all data contains some NaN, select those where only a few radial points are NaN
0259                         nrnan=sum(isnan(neraw.data));
0260                         frnan=nrnan/size(neraw.data,1);
0261                         fprintf('NaN values in NE profiles at all points in time, NaN fractions: %s\n',mat2str(frnan));
0262                         fprintf('Using profiles with less than 50%% NaN, removing radial points with NaN values.\n');
0263                         jnotnan=(frnan<0.5); % less than 50% NaN
0264                         rnotnan=~any(isnan(neraw.data(:,jnotnan)'));
0265                     else % standard case, can take all radial points
0266                         rnotnan=1:length(neraw.rho);
0267                     end
0268                     [d.ne.time,d.ne.data]=gti_chdtres(neraw.time(jnotnan),neraw.data(rnotnan,jnotnan),t,dtres,dtreslim);
0269                     d.ne.rho=neraw.rho(rnotnan);
0270                 end
0271             end
0272         case 'TE'
0273             d.te.shot=shot;
0274             if(docrca)
0275                 d.te.rawdiag={'Thomson','LIUQE'};
0276                 teraw=luke_get_nete_raw(diag.name,shot);
0277                 d.te.source=teraw.source;
0278                 fprintf('  . Using %s profiles from %s for TE.\n',d.te.rawdiag{1},d.te.source);
0279                 jnotnan=~any(isnan(teraw.data));
0280                 d.te.atime=tca;
0281                 d.te.otime=teraw.time(jnotnan);
0282                 d.te.rtime=tcr;
0283                 d.te.time=d.te.rtime;
0284                 d.te.data=gti_crcs(teraw.time(jnotnan),teraw.data(:,jnotnan),tca,tcr,dtreslim);
0285                 d.te.rho=teraw.rho;
0286             else % if no coherent/conditional averaging:
0287 %             if(dtres_bt(0.016,t,dtres,dtreslim))
0288 %                 d.te.rawdiag={'DMPX','LIUQE','Thomson'};
0289 %                 par=
0290 %             else
0291                 d.te.rawdiag={'Thomson','LIUQE'};
0292                 teraw=luke_get_nete_raw(diag.name,shot);
0293                 d.te.source=teraw.source;
0294                 fprintf('  . Using %s profiles from %s for TE.\n',d.te.rawdiag{1},d.te.source);
0295                 jnotnan=~any(isnan(teraw.data));
0296                 if(all(~jnotnan)) % if all data contains some NaN, select those where only a few radial points are NaN
0297                     nrnan=sum(isnan(teraw.data));
0298                     frnan=nrnan/size(teraw.data,1);
0299                     fprintf('NaN values in TE profiles at all points in time, NaN fractions: %s\n',mat2str(frnan));
0300                     fprintf('Using profiles with less than 50%% NaN, removing radial points with NaN values.\n');
0301                     jnotnan=(frnan<0.5); % less than 50% NaN
0302                     rnotnan=~any(isnan(teraw.data(:,jnotnan)'));
0303                 else % standard case, can take all radial points
0304                     rnotnan=1:length(teraw.rho);
0305                 end
0306                 [d.te.time,d.te.data]=gti_chdtres(teraw.time(jnotnan),teraw.data(rnotnan,jnotnan),t,dtres,dtreslim);
0307                 d.te.rho=teraw.rho(rnotnan);
0308 %             end
0309             end
0310         case 'TI'
0311             d.ti.shot=shot;
0312             % for the moment, do not set T_i (will lead to usage of T_e normalization)
0313             d.ti.rawdiag='None';
0314             d.ti.time=mean(t);
0315             d.ti.rho=linspace(0,1,41);
0316             d.ti.data=NaN(41,1);
0317         case 'ZEFF'
0318             d.zeff.shot=shot;
0319             % for the moments use the single value from proffit
0320             % (ZX if proffit not available), should to be replaced by a
0321             % multi-diagnostic Z_eff profile
0322             d.zeff.rawdiag={'XTOMO','Thomson','LIUQE','FIR'};
0323             mdsopen(shot);
0324             zeffproffit=tdi('\results::proffit:z_eff');
0325             mdsclose;
0326             if(mod(zeffproffit.status,2) && ~isempty(zeffproffit.data) && ~all(all(isnan(zeffproffit.data))));
0327                 d.zeff.source='proffit';
0328                 zeff_tdi=zeffproffit;
0329             else
0330                 mdsopen(shot);
0331                 zeff_zx=tdi('\results::zx');
0332                 mdsclose;
0333                 d.zeff.source='zx';
0334                 zeff_tdi=zeff_zx;
0335             end
0336             fprintf('  . Using Z_eff value from %s.\n',d.zeff.source);
0337             if(docrca)
0338                 d.zeff.atime=tca;
0339                 d.zeff.otime=zeff_tdi.dim{1};
0340                 d.zeff.rtime=tcr;
0341                 d.zeff.time=d.zeff.rtime;
0342                 d.zeff.data=gti_crca(zeff_tdi.dim{1},zeff_tdi.data(:,1)',tca,tcr,dtreslim);
0343             else % if no coherent/conditional averaging:
0344                 [d.zeff.time,d.zeff.data]=gti_chdtres(zeff_tdi.dim{1},zeff_tdi.data(:,1)',t,dtres,dtreslim);
0345             end
0346         case 'ZAXIS'
0347             d.zaxis.shot=shot;
0348             d.zaxis.rawdiag={'LIUQE'};
0349             mdsopen(shot);
0350             zraw=tdi('\results::z_axis');
0351             mdsclose;
0352             if(docrca)
0353                 d.zaxis.atime=tca;
0354                 d.zaxis.otime=zraw.dim{1};
0355                 d.zaxis.rtime=tcr;
0356                 d.zaxis.time=d.zaxis.rtime;
0357                 d.zaxis.data=gti_crca(zraw.dim{1},zraw.data',tca,tcr,dtreslim);
0358             else % if no coherent/conditional averaging:
0359                 [d.zaxis.time,d.zaxis.data]=gti_chdtres(zraw.dim{1},zraw.data',t,dtres,dtreslim);
0360             end
0361         case 'ECRH'
0362             d.ecrh.shot=shot;
0363             d.ecrh.rawdiag={'ECRH'};
0364             d.ecrh.isel=isel;
0365             mdsopen(shot);
0366             praw=tdi('\results::toray.input:p_gyro');
0367             mdsclose;
0368             % get the polarisation
0369             if(docrca)
0370             else % if no coherent/conditional averaging:
0371                 ox_t=gti_chdtres([],[],t,dtres,dtreslim);
0372                 if(flagTPG)
0373                     try
0374                         ecpol_addpathTPG;
0375                         [ox_xfrac_TPG,ox_times_TPG,~]=ecpol_ox_fraction(shot,0,ox_t);
0376                         ecpol_rmpathTPG;
0377                     catch ex
0378                         fprintf('Problem with ecpol_ox_fraction of TPG, writing NaN to X-mode fraction\nError message:\n%s\n',ex.message);
0379                         ecpol_rmpathTPG;
0380                         ox_xfrac_TPG=NaN;
0381                         ox_times_TPG=NaN;
0382                     end
0383                 end
0384                 if(flagFF)
0385                     ox_xfrac_FF=nan(length(ox_t),length(isel));
0386                     try
0387                         jnotallnan=false(1,length(isel));
0388                         for j=1:length(isel)
0389                             jnotallnan(j)=~all(isnan(praw.data(:,isel(j))));
0390                         end
0391                         juse=find(jnotallnan);
0392                         ox_xfrac_FF(:,juse)=ecpol_ox_fraction(shot,ox_t,isel(juse));
0393                     catch ex
0394                         fprintf('Problem with ecpol_ox_fraction of FF, writing NaN to X-mode fraction\nError message:\n%s\n',ex.message);
0395                     end
0396                 end
0397             end
0398             for j=1:length(isel)    % do a loop over all selected launchers
0399                 k=isel(j);
0400                 % get power: independent of harmonic, kW->W
0401                 itknotnan=~isnan(praw.data(:,k));
0402                 if(docrca)
0403                     d.ecrh.pow{k}.atime=tca;
0404                     if(~any(itknotnan))
0405                         d.ecrh.pow{k}.otime=mean(t);
0406                     else
0407                         d.ecrh.pow{k}.otime=praw.dim{1}(itknotnan);
0408                     end
0409                     d.ecrh.pow{k}.rtime=tcr;
0410                     d.ecrh.pow{k}.time=d.ecrh.pow{k}.rtime;
0411                     if(~any(itknotnan))
0412                         d.ecrh.pow{k}.data=gti_crca(mean(t),NaN,tca,tcr,dtreslim);
0413                     else
0414                         d.ecrh.pow{k}.data=gti_crca(praw.dim{1}(itknotnan),1e3*praw.data(itknotnan,k)',tca,tcr,dtreslim);
0415                     end
0416                 else % if no coherent/conditional averaging:
0417                     if(~any(itknotnan))
0418                         [d.ecrh.pow{k}.time,d.ecrh.pow{k}.data]=gti_chdtres(mean(t),NaN,t,dtres,dtreslim);
0419                     else
0420                         [d.ecrh.pow{k}.time,d.ecrh.pow{k}.data]=gti_chdtres(praw.dim{1}(itknotnan),1e3*praw.data(itknotnan,k)',t,dtres,dtreslim);
0421                     end
0422                 end
0423                 if(k<7) % X2 launchers 1-6 ang and stray
0424                     % constant data (from collect_parameters_for_luke.m)
0425                     ang.freq=82.7e9;
0426                     ang.P1_R=1.3072;
0427                     ang.P1_theta=52.013/180*pi;
0428                     ang.P2_R=1.1906;
0429                     ang.P2_d=0.01;
0430                     if(j==1 || j==4)
0431                         ang.P1_Z=-0.05865;
0432                         ang.P2_Z=0.0645;
0433                         ang.P3_Z=-0.0025;
0434                     else
0435                         ang.P1_Z=0.39885;
0436                         ang.P2_Z=0.522;
0437                         ang.P3_Z=0.455;
0438                     end
0439                     ang.P3_tor=0.0;
0440                     % data from tree
0441                     mdsopen(shot);
0442                     thetaraw=mdsvalue(sprintf('_theta=\\ecrh::launchers:theta_l:x2_%d',k))'*pi/180;
0443                     timeraw=mdsvalue('dim_of(_theta)');
0444                     phiraw=mdsvalue(sprintf('\\ecrh::launchers:phi_l:x2_%d',k))'*pi/180;
0445                     if(ischar(timeraw))
0446                         fprintf('Setting angles of launcher %d to NaN because there is no data\n',k);
0447                         thetaraw=NaN;
0448                         timeraw=NaN;
0449                         phiraw=NaN;
0450                     end
0451                     try
0452                         nodetmp=sprintf('\\ecrh::ecrh_power:stray:x2_%d',k);
0453                         strayraw=tdi(nodetmp);
0454                         if(mod(strayraw.status,2)==0)
0455                             error('No data in %s',nodetmp);
0456                         end
0457                     catch ex
0458                         fprintf('Setting stray signal of launcher %d to NaN because there is no data due to the error:\n%s\n',k,ex.message);
0459                         strayraw.data=NaN;
0460                         strayraw.dim{1}=NaN;
0461                     end
0462                     mdsclose;
0463                     if(docrca)
0464                         ang.atime=tca;
0465                         ang.otime=timeraw;
0466                         ang.rtime=tcr;
0467                         ang.time=ang.rtime;
0468                         allangdata=gti_crca(timeraw,[thetaraw;phiraw],tca,tcr,dtreslim);
0469                         stray.atime=tca;
0470                         stray.otime=strayraw.dim{1};
0471                         stray.rtime=tcr;
0472                         stray.time=stray.rtime;
0473                         stray.data=gti_crca(strayraw.dim{1},strayraw.data',tca,tcr,dtreslim);
0474                     else % if no coherent/conditional averaging:
0475                         [ang.time,allangdata]=gti_chdtres(timeraw,[thetaraw;phiraw],t,dtres,dtreslim);
0476                         [stray.time,stray.data]=gti_chdtres(strayraw.dim{1},strayraw.data',t,dtres,dtreslim);
0477                     end
0478                     ang.theta=allangdata(1,:);
0479                     ang.phi=allangdata(2,:);
0480                     ang.theta_eff=-atan(ang.theta).*cos(ang.phi);
0481                     ang.phi_eff=atan(tan(ang.theta).*sin(ang.phi))./sqrt(1+tan(ang.theta).^2.*cos(ang.phi).^2);
0482                     % ox fraction data
0483                     if(docrca)
0484                         % assume 100% X-mode for the moment (O-mode fraction only on request since takes a long time: need to specify parameter or shot range if it should be read)
0485                         d.ecrh.xfrac{k}.time=tcr;
0486                         d.ecrh.xfrac{k}.data=ones(1,numel(tcr));
0487                     else % if no coherent/conditional averaging:
0488                         % put into d.ecrh.xfrac{k}.time and .data
0489                         if(flagTPG)
0490                             d.ecrh.xfrac_TPG{k}.time=ox_t;
0491                             d.ecrh.xfrac_TPG{k}.data=NaN(1,numel(ox_t));
0492                             try
0493                                 for ell=1:length(ox_t)
0494                                     if(any(ox_times_TPG==ox_t(ell)))
0495                                         d.ecrh.xfrac_TPG{k}.data(1,ell)=ox_xfrac_TPG(k,ox_times_TPG==ox_t(ell));
0496                                     end
0497                                 end
0498                             catch ex
0499                                 fprintf('Problem with output of ecpol_ox_fraction of TPG, writing NaN to X-mode fraction\nError message:\n%s\n',ex.message);
0500                             end
0501                             if(~flagFF)
0502                                 d.ecrh.xfrac{k}=d.ecrh.xfrac_TPG{k};
0503                             end
0504                         end
0505                         if(flagFF)
0506                             d.ecrh.xfrac_FF{k}.time=ox_t;
0507                             d.ecrh.xfrac_FF{k}.data=ox_xfrac_FF(:,j)';
0508                             if(~flagTPG)
0509                                 d.ecrh.xfrac{k}=d.ecrh.xfrac_FF{k};
0510                             end
0511                         end
0512                     end
0513                 else    % X3 launchers 7-9 ang and stray
0514                     % constant data (from collect_parameters_for_luke.m)
0515                     ang.freq=118.e9;
0516                     ang.P1_Z=.9885;
0517                     ang.P1_theta=0;
0518                     ang.P2_Z=ang.P1_Z;
0519                     ang.P2_d=0;
0520                     ang.P3_Z=ang.P1_Z;
0521                     % data from tree
0522                     mdsopen(shot);
0523                     datatemp=mdsvalue('_theta=\ecrh::launchers:theta_l:x3');
0524                     if(ischar(datatemp))
0525                         thetaraw=NaN;
0526                         timeraw=mean(t);
0527                     else
0528                         thetaraw=datatemp'*pi/180;
0529                         timeraw=mdsvalue('dim_of(_theta)');
0530                     end
0531                     datatemp=mdsvalue(sprintf('\\ecrh::launchers:phi_l:x3_%d',k));
0532                     if(ischar(datatemp))
0533                         phiraw=NaN(size(timeraw));
0534                     else
0535                         phiraw=datatemp*pi/180*ones(size(thetaraw));
0536                     end
0537                     datatemp=mdsvalue('\ecrh::launchers:R:x3');
0538                     if(ischar(datatemp))
0539                         rraw=NaN(size(timeraw));
0540                     else
0541                         rraw=datatemp';
0542                     end
0543                     try
0544                         nodetmp=sprintf('\\ecrh::ecrh_power:stray:x3_%d',k);
0545                         strayraw=tdi(nodetmp);
0546                         if(mod(strayraw.status,2)==0)
0547                             error('No data in %s',nodetmp);
0548                         end
0549                     catch ex
0550                         fprintf('Setting stray signal of launcher %d to NaN because there is no data due to the error:\n%s\n',k,ex.message);
0551                         strayraw.data=NaN;
0552                         strayraw.dim{1}=NaN;
0553                     end
0554                     mdsclose;
0555                     if(docrca)
0556                         ang.atime=tca;
0557                         ang.otime=timeraw;
0558                         ang.rtime=tcr;
0559                         ang.time=ang.rtime;
0560                         allangdata=gti_crca(timeraw,[thetaraw;phiraw;rraw],tca,tcr,dtreslim);
0561                         stray.atime=tca;
0562                         stray.otime=strayraw.dim{1};
0563                         stray.rtime=tcr;
0564                         stray.time=stray.rtime;
0565                         stray.data=gti_crca(strayraw.dim{1},strayraw.data',tca,tcr,dtreslim);
0566                     else % if no coherent/conditional averaging:
0567                         [ang.time,allangdata]=gti_chdtres(timeraw,[thetaraw;phiraw;rraw],t,dtres,dtreslim);
0568                         [stray.time,stray.data]=gti_chdtres(strayraw.dim{1},strayraw.data',t,dtres,dtreslim);
0569                     end
0570                     thetaraw=allangdata(1,:);
0571                     phiraw=allangdata(2,:);
0572                     ang.P1_R=allangdata(3,:);
0573                     ang.P2_R=ang.P1_R;
0574                     ang.P3_tor=(ang.P1_R+0.125).*tan(phiraw);
0575                     ang.theta=atan2(sqrt(sin(thetaraw).^2+tan(phiraw).^2),cos(thetaraw));
0576                     ang.phi=atan2(tan(phiraw),sin(thetaraw));
0577                     ang.theta_eff=-thetaraw;
0578                     ang.phi_eff=atan2(tan(phiraw),abs(cos(thetaraw)));
0579                 end
0580                 d.ecrh.ang{k}=ang;
0581                 d.ecrh.stray{k}=stray;
0582                 clear ang stray;
0583             end
0584         case 'EQUIL'
0585             d.equil.shot=shot;
0586             d.equil.rawdiag={'LIUQE'};
0587             if(docrca)
0588                 % Preliminary solution: use equilibria around the middle
0589                 % tca point
0590                 d.equil.atime=tca;
0591                 d.equil.otime=tca(ceil(length(tca)/2))+tcr;
0592                 d.equil.rtime=tcr;
0593                 d.equil.time=d.equil.rtime;
0594                 % get times where equilibria exist already
0595                 mdsopen(shot);
0596                 t_liuqe=mdsvalue('dim_of(\results::i_p)');
0597                 mdsclose;
0598                 if(length(d.equil.otime)==1 || ~((max(diff(d.equil.otime))/min(diff(d.equil.otime))-1)<5*eps) ) % one point or not equidistant d.equil.otime
0599                     for it=1:length(d.equil.otime)
0600                         if(min(d.equil.otime(it)-t_liuqe)>eps(single(1)))
0601                             fprintf('  . LIUQE not available in MDS tree for t=%gs.\n',d.equil.otime(it));
0602                             job = matstandardjob_jd('runmat');
0603                             job.funcopy=false;
0604                             tempcommand=sprintf('/home/codes/liuqe/liuqe %d %g %g %d',shot,d.equil.otime(it),d.equil.otime(it),1);
0605                             sprintf('  . running LIUQE on lac5 for requested time point %gs.\n',d.equil.otime(it));
0606                             if(~exist('matremote_jd','file'))
0607                                 pwd_save=pwd;
0608                                 cd('/home/decker');
0609                                 startup;
0610                                 cd(pwd_save);
0611                             end
0612                             [err,res,text]=matremote_jd('system',{tempcommand},feval('remote_CRPP_LAC5','','local'),job);
0613                             % check if LIUQE ran on lac5 (warning if not)
0614                             if(~isempty(err))
0615                                 warning('matremote_jd:LIUQE:failed','The LIUQE run on lac5 failed!');
0616                                 disperr_jd(err);
0617                             elseif(res)
0618                                 disp(text);
0619                                 warning('matremote_jd:LIUQE:failed','The LIUQE run on lac5 failed, see output above!');
0620                             end
0621                         else
0622                             fprintf('  . LIUQE available in MDS tree for t=%gs.\n',d.equil.otime(it));
0623                         end
0624                     end
0625                 else % more than one point and equidistant d.equil.otime
0626                     t1=d.equil.otime(1);
0627                     t2=d.equil.otime(end);
0628                     n12=length(d.equil.otime);
0629                     % check if equilibria exist already
0630                     flag_run=false(size(d.equil.otime));
0631                     for it=length(d.equil.otime)
0632                         if(min(d.equil.otime(it)-t_liuqe)>eps(single(1)))
0633                             fprintf('  . LIUQE not available in MDS tree for t=%gs.\n',d.equil.otime(it));
0634                             flag_run(it)=true;
0635                         end
0636                     end
0637                     if(any(flag_run))
0638                         % run LIUQE if not all equilibria available (using matremote_jd)
0639                         job = matstandardjob_jd('runmat');
0640                         job.funcopy=false;
0641                         if(sum(flag_run)==1)
0642                             tempcommand=sprintf('/home/codes/liuqe/liuqe %d %g %g 1',shot,d.equil.otime(flag_run),d.equil.otime(flag_run));
0643                             sprintf('  . running LIUQE on lac5 for requested time point %gs.\n',d.equil.otime(flag_run));
0644                         elseif(sum(flag_run)==2)
0645                             t_temp=d.equil.otime(flag_run);
0646                             tempcommand=sprintf('/home/codes/liuqe/liuqe %d %g %g 2',shot,t_temp(1),t_temp(2));
0647                             sprintf('  . running LIUQE on lac5 for requested time points %gs and %gs.\n',t_temp(1),t_temp(2));
0648                         else
0649                             tempcommand=sprintf('/home/codes/liuqe/liuqe %d %g %g %d',shot,t1,t2,n12);
0650                             sprintf('  . running LIUQE on lac5 for all requested times linspace(%g,%g,%d)s.\n',t1,t2,n12);
0651                         end
0652                         if(~exist('matremote_jd','file'))
0653                             pwd_save=pwd;
0654                             cd('/home/decker');
0655                             startup;
0656                             cd(pwd_save);
0657                         end
0658                         [err,res,text]=matremote_jd('system',{tempcommand},feval('remote_CRPP_LAC5','','local'),job);
0659                         % check if LIUQE ran on lac5 (warning if not)
0660                         if(~isempty(err))
0661                             warning('matremote_jd:LIUQE:failed','The LIUQE run on lac5 failed!');
0662                             disperr_jd(err);
0663                         elseif(res)
0664                             disp(text);
0665                             warning('matremote_jd:LIUQE:failed','The LIUQE run on lac5 failed, see output above!');
0666                         end
0667                     else
0668                         fprintf('  . LIUQE already available in MDS tree for all requested times: no LIUQE run.\n');
0669                     end
0670                 end
0671                 % create and read the equilibrium files
0672                 for it=1:length(d.equil.otime)
0673                     % create (need to call for single times, otherwise all equils but the first are corrupted)
0674                     make_eqdsk(shot,d.equil.otime(it),outputpath,'TCV');
0675                     % read and put into eqdsk field
0676                     eqdsk_file = [outputpath,filesep,'eqdsk.',num2str(d.equil.shot),'t',num2str(d.equil.otime(it),'%6.4f')];
0677                     if exist(eqdsk_file,'file'),
0678                         d.equil.eqdsk{it} = fileread(eqdsk_file);% magnetic equilibrium
0679                     else
0680                         warning(['File ',eqdsk_file,' could not be found.']);
0681                     end
0682                 end
0683                 %TODO: better solution using LIUQE in matlab with
0684                 %conditionally averaged experimental data instead
0685             else % if no coherent/conditional averaging:
0686                 fprintf('- Checking availability of LIUQE equilibria for shot %d at requested time points\n',shot);
0687                 % get times where equilibria exist already
0688                 mdsopen(shot);
0689                 t_liuqe=mdsvalue('dim_of(\results::i_p)');
0690                 mdsclose;
0691                 if( isempty(dtres) && (length(t)==1 || ~((max(diff(t))/min(diff(t))-1)<5*eps)) )
0692                     % specify times
0693                     d.equil.time=t;
0694                     for it=1:length(d.equil.time)
0695                         if(min(d.equil.time(it)-t_liuqe)>eps(single(1)))
0696                             fprintf('  . LIUQE not available in MDS tree for t=%gs.\n',d.equil.time(it));
0697                             job = matstandardjob_jd('runmat');
0698                             job.funcopy=false;
0699                             tempcommand=sprintf('/home/codes/liuqe/liuqe %d %g %g %d',shot,d.equil.time(it),d.equil.time(it),1);
0700                             sprintf('  . running LIUQE on lac5 for requested time point %gs.\n',d.equil.time(it));
0701                             if(~exist('matremote_jd','file'))
0702                                 pwd_save=pwd;
0703                                 cd('/home/decker');
0704                                 startup;
0705                                 cd(pwd_save);
0706                             end
0707                             [err,res,text]=matremote_jd('system',{tempcommand},feval('remote_CRPP_LAC5','','local'),job);
0708                             % check if LIUQE ran on lac5 (warning if not)
0709                             if(~isempty(err))
0710                                 warning('matremote_jd:LIUQE:failed','The LIUQE run on lac5 failed!');
0711                                 disperr_jd(err);
0712                             elseif(res)
0713                                 disp(text);
0714                                 warning('matremote_jd:LIUQE:failed','The LIUQE run on lac5 failed, see output above!');
0715                             end
0716                         else
0717                             fprintf('  . LIUQE available in MDS tree for t=%gs.\n',d.equil.time(it));
0718                         end
0719                     end
0720                 else
0721                     % specify times
0722                     if(isempty(dtres) || length(t)==1) % equidistant with dtres=[]
0723                         d.equil.time=t;
0724                         t1=d.equil.time(1);
0725                         t2=d.equil.time(end);
0726                         n12=length(d.equil.time);
0727                     else % equidistant with dtres specified
0728                         z1=ceil(t(1)/dtres-5*eps);
0729                         z2=floor(t(end)/dtres+5*eps);
0730                         t1=z1*dtres;
0731                         t2=z2*dtres;
0732                         n12=z2-z1+1;
0733                         d.equil.time=linspace(t1,t2,n12);
0734                     end
0735                     % check if equilibria exist already
0736                     flag_run=false(size(d.equil.time));
0737                     for it=length(d.equil.time)
0738                         if(min(d.equil.time(it)-t_liuqe)>eps(single(1)))
0739                             fprintf('  . LIUQE not available in MDS tree for t=%gs.\n',d.equil.time(it));
0740                             flag_run(it)=true;
0741                         end
0742                     end
0743                     if(any(flag_run))
0744                         % run LIUQE if not all equilibria available (using matremote_jd)
0745                         job = matstandardjob_jd('runmat');
0746                         job.funcopy=false;
0747                         if(sum(flag_run)==1)
0748                             tempcommand=sprintf('/home/codes/liuqe/liuqe %d %g %g 1',shot,d.equil.time(flag_run),d.equil.time(flag_run));
0749                             sprintf('  . running LIUQE on lac5 for requested time point %gs.\n',d.equil.time(flag_run));
0750                         elseif(sum(flag_run)==2)
0751                             t_temp=d.equil.time(flag_run);
0752                             tempcommand=sprintf('/home/codes/liuqe/liuqe %d %g %g 2',shot,t_temp(1),t_temp(2));
0753                             sprintf('  . running LIUQE on lac5 for requested time points %gs and %gs.\n',t_temp(1),t_temp(2));
0754                         else
0755                             tempcommand=sprintf('/home/codes/liuqe/liuqe %d %g %g %d',shot,t1,t2,n12);
0756                             sprintf('  . running LIUQE on lac5 for all requested times linspace(%g,%g,%d)s.\n',t1,t2,n12);
0757                         end
0758                         if(~exist('matremote_jd','file'))
0759                             pwd_save=pwd;
0760                             cd('/home/decker');
0761                             startup;
0762                             cd(pwd_save);
0763                         end
0764                         [err,res,text]=matremote_jd('system',{tempcommand},feval('remote_CRPP_LAC5','','local'),job);
0765                         % check if LIUQE ran on lac5 (warning if not)
0766                         if(~isempty(err))
0767                             warning('matremote_jd:LIUQE:failed','The LIUQE run on lac5 failed!');
0768                             disperr_jd(err);
0769                         elseif(res)
0770                             disp(text);
0771                             warning('matremote_jd:LIUQE:failed','The LIUQE run on lac5 failed, see output above!');
0772                         end
0773                     else
0774                         fprintf('  . LIUQE already available in MDS tree for all requested times: no LIUQE run.\n');
0775                     end
0776                 end
0777                 % create and read the equilibrium files
0778                 for it=1:length(d.equil.time)
0779                     % create (need to call for single times, otherwise all equils but the first are corrupted)
0780                     make_eqdsk(shot,d.equil.time(it),outputpath,'TCV');
0781                     % read and put into eqdsk field
0782                     eqdsk_file = [outputpath,filesep,'eqdsk.',num2str(d.equil.shot),'t',num2str(d.equil.time(it),'%6.4f')];
0783                     if exist(eqdsk_file,'file'),
0784                         d.equil.eqdsk{it} = fileread(eqdsk_file);% magnetic equilibrium
0785                     else
0786                         warning(['File ',eqdsk_file,' could not be found.']);
0787                     end
0788                 end
0789             end
0790         case 'HXR'
0791             [d.hxr.hxr,d.hxr.hxrparam,d.hxr.ap]=luke_get_hxr_diag_tcv(shot);
0792             Etphot=d.hxr.hxr.kdiag_hxr(:,end); % the interval for photon temperature should be kept
0793             if(shot<40000) % HXRC
0794                 [tcbn,cbn,Ebin]=hxrc_getcbn(shot,d.hxr.ap.chords{1});
0795                 if(docrca) % use gti_crca for conditional averaging. Best possible solution, but for most cases the time resolution of the underlying data is too low.
0796                     d.hxr.atime=tca;
0797                     d.hxr.otime=tcbn;
0798                     d.hxr.rtime=tcr;
0799                     d.hxr.time=d.hxr.rtime;
0800                     d.hxr.data=gti_crca(tcbn,cbn,tca,tcr,dtreslim);
0801                 else % if no coherent/conditional averaging:
0802                     [d.hxr.time,d.hxr.data]=gti_chdtres(tcbn,cbn,t,dtres,dtreslim);
0803                 end
0804             else % HXRS
0805                 Ebin=d.hxr.hxr.kdiag_hxr(1,1:(end-1));
0806                 % use the vrca version of hxrs_getcbn
0807                 if(docrca)
0808                     d.hxr.atime=tca;
0809                     d.hxr.rtime=tcr;
0810                     [d.hxr.time,d.hxr.data]=hxrs_getcbn(shot,d.hxr.ap.cams,d.hxr.ap.chords,[],Ebin,[],tca,tcr);
0811                 else % if no coherent/conditional averaging:
0812                     tcr=gti_chdtres([],[],t,dtres,dtreslim);
0813                     [d.hxr.time,d.hxr.data]=hxrs_getcbn(shot,d.hxr.ap.cams,d.hxr.ap.chords,[],Ebin,[],0,tcr);
0814                 end
0815 %                 % for the moment, use high time resolution and chdtres
0816 %                 [tcbn,cbn]=hxrs_getcbn(shot,d.hxr.ap.cams,d.hxr.ap.chords,[],1e-4,Ebin);
0817 %                 [d.hxr.time,d.hxr.data]=gti_chdtres(tcbn,cbn,t,dtres,dtreslim);
0818             end
0819             % redefine kdiag_hxr parameter according to measurement
0820             d.hxr.hxr.kdiag_hxr=[Ebin,Etphot(1);Ebin(2:end),2*Ebin(end)-Ebin(end-1),Etphot(2)];
0821     end
0822 end
0823 
0824 
0825 end
0826 
0827 function neteraw=luke_get_nete_raw(nete,shot)
0828     if(strcmp(nete,'NE'))
0829         proffitnode='\results::proffit.avg_time:neft_abs';
0830         autonode='\results::thomson.profiles.auto:ne:foo';
0831     elseif(strcmp(nete,'TE'))
0832         proffitnode='\results::proffit.avg_time:teft';
0833         autonode='\results::thomson.profiles.auto:te:foo';
0834     end
0835     % special shots: no Thomson data -> take data from very similar shot
0836     if(shot==49115)
0837         shot=49111;
0838         warning('No Thomson data for shot 49115, taking data from similar shot 49111 instead!');
0839     end
0840     % get proffit data
0841     mdsopen(shot);
0842     neteproffit=tdi(proffitnode);
0843     mdsclose;
0844     % check proffit data
0845     takeproffit=(mod(neteproffit.status,2) && ~isempty(neteproffit.data) && ~all(all(isnan(neteproffit.data))));
0846     if(takeproffit)
0847         % put time coordinate in second dimension, radial in
0848         % first dimension
0849         it=find(strcmp('s',neteproffit.dimunits) | strcmp('time [s]',neteproffit.dimunits));
0850         neteraw.time=neteproffit.dim{it};
0851         neteraw.rho=neteproffit.dim{3-it};
0852         if(it==1)
0853             neteraw.data=neteproffit.data';
0854         else
0855             neteraw.data=neteproffit.data;
0856         end
0857         neteraw.source='proffit';
0858     else
0859         % get Thomson autofits
0860         mdsopen(shot);
0861         neteauto=tdi(autonode);
0862         mdsclose;
0863         % check autofits data
0864         takeauto=(mod(neteauto.status,2) & ~isempty(neteauto.data) & ~all(all(isnan(neteauto.data))));
0865         if(takeauto)
0866             % put time coordinate in second dimension, radial in
0867             % first dimension
0868             it=find(strcmp('s',neteauto.dimunits) | strcmp('time [s]',neteauto.dimunits));
0869             neteraw.time=neteauto.dim{it};
0870             neteraw.rho=neteauto.dim{3-it};
0871             if(it==1)
0872                 neteraw.data=neteauto.data';
0873             else
0874                 neteraw.data=neteauto.data;
0875             end
0876             neteraw.source='autofits';
0877         else
0878             error('No Thomson data available - try to run proffit or to reset autofits');
0879             % take FIR data (to be implemented later)
0880         end
0881     end
0882 end
0883 
0884 
0885 function flag=dtres_bt(val,t,dtres,dtreslim)
0886     flag=false;
0887     if(~isempty(dtres)) % specification by dtres
0888         if(~isnan(dtres(1)) && (dtres(1)<val)) % dtres(1)
0889             flag=true;
0890         elseif((length(dtres)==2) && ~isnan(dtres(2)) && (dtres(2)<val) ) % dtres(2)
0891             flag=true;
0892         elseif((length(dtres)==3) && ~all(isnan(dtres(2:3))) && (nanmean(dtres(2:3))<val)) % dtres(2:3)
0893             flag=true;
0894         end
0895     else % specification by dtreslim
0896         if(~isempty(dtreslim))
0897             if(all(size(dtreslim)==[1,1]))
0898                 if(~isnan(dtreslim) && (dtreslim<val))
0899                     flag=true;
0900                 end
0901             elseif(all(size(dtreslim)==[1,2]))
0902                 if(~all(isnan(dtreslim)) && (nanmean(dtreslim)<val))
0903                     flag=true;
0904                 end
0905             elseif(all(size(dtreslim)==[2,1]))
0906                 if(~isnan(dtreslim(1)) && dtreslim(1)<val)
0907                     flag=true;
0908                 end
0909             elseif(all(size(dtreslim)==[2,2]))
0910                 if(~all(isnan(dtreslim(1,:))) && (nanmean(dtreslim(1,:))<val))
0911                     flag=true;
0912                 end
0913             end
0914         end
0915     end
0916     if(length(t)>1 && (min(diff(t))<val)) % specification by t spacing
0917         flag=true;
0918     end
0919 end