collect_parameters_for_luke

    function dummy = collect_parameters_for_luke (shot,t1,t2,laun,filename,teflag,zeff)

    collect all TCV data needed for Luke run (with the exception of the
      equilibrium, which must be supplied separately in EQDSK format:
      use make_eqdsk(shot,times,outdir,machine,liuqetype) for this)
    input:
        shot = shot number (vector)
        t1, t2 = start and end times (vectors): data will be averaged
               from t1 to t2
        laun = launcher numbers of interest (vector)
        filename = filename for output data storage
        teflag = flag, if 0 (default) Thomson data are used for Te,
             if nonzero Thomson data are renormalized to Te(0)
             supplied by XTe
        zeff = structure array (length = shot):
           zeff(i).data = zeff profile (dimensions = rho x time)
           zeff(i).rho = rho_psi grid
           zeff(i).time = time vector
             where i is the shot number
             the output will be remapped on the Thomson grid
             if this variable is empty or undefined, zeff will be 
               read from data but will then be scalar
    rho coordinate is rho_psi

0001 function collect_parameters_for_luke (shot,t1,t2,laun,filename,teflag,zeff)
0002 
0003 %
0004 %    function dummy = collect_parameters_for_luke (shot,t1,t2,laun,filename,teflag,zeff)
0005 %
0006 %    collect all TCV data needed for Luke run (with the exception of the
0007 %      equilibrium, which must be supplied separately in EQDSK format:
0008 %      use make_eqdsk(shot,times,outdir,machine,liuqetype) for this)
0009 %    input:
0010 %        shot = shot number (vector)
0011 %        t1, t2 = start and end times (vectors): data will be averaged
0012 %               from t1 to t2
0013 %        laun = launcher numbers of interest (vector)
0014 %        filename = filename for output data storage
0015 %        teflag = flag, if 0 (default) Thomson data are used for Te,
0016 %             if nonzero Thomson data are renormalized to Te(0)
0017 %             supplied by XTe
0018 %        zeff = structure array (length = shot):
0019 %           zeff(i).data = zeff profile (dimensions = rho x time)
0020 %           zeff(i).rho = rho_psi grid
0021 %           zeff(i).time = time vector
0022 %             where i is the shot number
0023 %             the output will be remapped on the Thomson grid
0024 %             if this variable is empty or undefined, zeff will be
0025 %               read from data but will then be scalar
0026 %    rho coordinate is rho_psi
0027 
0028 %    S. Coda     30/04/07
0029 %               21/06/11
0030 %   J. Kamleitner   18/10/13
0031 %                   06/11/13
0032 %
0033 %                   19/03/2014: new version based on luke_get_trdata
0034 %
0035 
0036 if(nargin<7 || isempty(zeff))
0037     zeff.data=NaN;
0038 end
0039 if(nargin<6 || isempty(teflag))
0040     teflag=0;
0041 end
0042 if(nargin<5 || isempty(filename))
0043     filename='luke_input_parameters.mat';
0044 end
0045 
0046 laun=laun(:)';
0047 
0048 t0=0.5*(t1+t2);
0049 
0050 newversion=true;
0051 
0052 if(newversion)
0053     a.shot=shot;
0054     a.t=t0;
0055     a.dtres=[];
0056     a.dtreslim=0.5*(t2-t1);
0057     a.diag{1}.name='EQUIL';
0058     a.diag{2}.name='IP';
0059     a.diag{3}.name='VLOOP';
0060     a.diag{4}.name='NE';
0061     a.diag{5}.name='TE';
0062     a.diag{6}.name='BT';
0063     a.diag{7}.name='ZEFF';
0064     a.diag{8}.name='ZAXIS';
0065     a.diag{9}.name='TI';
0066     a.diag{10}.name='ECRH';
0067     a.diag{10}.isel=laun;
0068     
0069     % get the data
0070     d=luke_get_trdata(a);
0071     
0072     % arrange the data for saving
0073     vloop=-d.vloop.data; % vloop sign is inverse to luke_get_trdata to keep it consistent with old LUKE version (returns in LUKE sign convention directly, luke_get_trdata returns everything in TCV sign conventions)
0074 %     nel=nel';
0075     te=d.te.data';
0076 %     dte=dte';
0077     ne=d.ne.data';
0078 %     dne=dne';
0079     rho=d.te.rho';
0080     ti=d.ti.data';
0081     if(~all(isnan(zeff.data)))
0082         % interpolate the zeff from input (as in old version)
0083         [~,jt1]=min(abs(zeff.time-t1));
0084         [~,jt2]=min(abs(zeff.time-t1));
0085         zefff=nanmean(interp1(zeff.rho(:),zeff.data(:,jt1:jt2),rho),2);
0086         clear zeff;
0087         zeff=zefff';
0088         clear zefff;
0089     else
0090         % get zeff data
0091         clear zeff;
0092         zeff=d.zeff.data';
0093     end
0094     zaxis=d.zaxis.data';
0095     ipsign=d.ipsign.data';
0096     btsign=d.btsign.data';
0097     Ip=d.ip.data';
0098     for j=laun
0099         ang=d.ecrh.ang{j};
0100         ang.P0=d.ecrh.pow{j}.data;
0101         ang.id=sprintf('%d',shot);
0102         eval(['angles_',int2str(j),'=ang;']);
0103         clear ang;
0104     end
0105     clear a d j;
0106     
0107 else % old version
0108     for i=1:length(shot)
0109     
0110         mdsopen(shot(i));
0111 
0112         % get B_T sign
0113         btsign(i)=sign(mdsdata('\magnetics::iphi[$]',(t1(i)+t2(i))/2));
0114 
0115         % get I_P sign
0116         ip0=mdsdata('\ipol[[0.0,0.01],"OH_001"]');
0117         if(ip0(2)>=ip0(1))
0118             ipsign(i)=-1;
0119         else
0120             ipsign(i)=1;
0121         end
0122         clear ip0;
0123 
0124         % get loop voltage and plasma current
0125         [~,vl]=tcvget('VLOOP',[t1(i):.002:t2(i)]);
0126         vloop(i)=-nanmean(vl);
0127         clear vl;
0128 
0129         % get plasma current
0130         nipavg=10;
0131         ipabs=ipsign(i)*mdsdata('\magnetics::iplasma:trapeze');
0132         tip=mdsdata('dim_of(\magnetics::iplasma:trapeze,0)');
0133         [~,itip0]=min(abs(tip-t0(i)));  % index of acquired time point closest to t0
0134         [~,itipmin(1)]=min(ipabs(tip<t1(i)));    % position of minimum (baseline) before plasma
0135         [~,itipmin(2)]=min(ipabs(tip>t2(i)));    % position of minimum (baseline) after plasma
0136         itipmin(2)=itipmin(2)+find(tip<=t2(i),1,'last');    %correct the shortening of the ipabs vector in the command above
0137         ipoff(1)=nanmean(ipabs(max(1,itipmin(1)-nipavg):itipmin(1)));
0138         ipoff(2)=nanmean(ipabs(itipmin(2):min(length(ipabs),itipmin(2)+nipavg)));
0139         Ip(i)=ipsign(i)*(nanmean(ipabs( tip>=t1(i) & tip<=t2(i) ))-interp1(tip(itipmin),ipoff,tip(itip0)));  % subtract linearly interpolated baseline off the averaged measured Ip value and multiply with ipsign again to get the correct value for Ip
0140         clear nipavg ipabs tip itip0 itipmin ipoff;
0141 
0142         % get electron density and temperature from Thomson measurements
0143         [~,netemp]=tcvget('NEL',[t1(i):.002:t2(i)]);
0144         nel(i)=nanmean(netemp);
0145         if(shot(i)<16000)
0146             teraw=mdsdata('\results::th_prof_te');
0147             dteraw=mdsdata('\results::th_prof_te:error_bar');
0148             neraw=mdsdata('\results::th_prof_ne');
0149             dneraw=mdsdata('\results::th_prof_ne:error_bar');
0150             rho(:,i)=mdsdata('dim_of(\results::th_prof_te,1)');
0151             tth=mdsdata('dim_of(\results::th_prof_te,0)');
0152         else
0153             teraw=mdsdata('\results::proffit.avg_time:teft')';
0154             dteraw=mdsdata('\results::proffit.avg_time:teft_std')';
0155             neraw=mdsdata('\results::proffit.avg_time:neft_abs')';
0156             dneraw=mdsdata('\results::proffit.avg_time:neft_std')'.*repmat(mdsdata('\results::proffit.avg_time:neft_firat'),1,size(neraw,2));
0157             rho(:,i)=mdsdata('dim_of(\results::proffit.avg_time:teft)');
0158             tth=mdsdata('dim_of(\results::proffit.avg_time:teft,1)');
0159         end
0160         teind1=find(tth>=t1(i)&tth<=t2(i), 1 );
0161         teind2=find(tth>=t1(i)&tth<=t2(i), 1, 'last' );
0162         if(~isempty(teind1))
0163             te(:,i)=nanmean(teraw(teind1:teind2,:),1)';
0164             dte(:,i)=max(std(dteraw(teind1:teind2,:),0,1),sqrt(nanmean(dteraw(teind1:teind2,:).^2,1)/(teind2-teind1)))';
0165             ne(:,i)=nanmean(neraw(teind1:teind2,:),1)';
0166             dne(:,i)=max(std(dneraw(teind1:teind2,:),0,1),sqrt(nanmean(dneraw(teind1:teind2,:).^2,1)/(teind2-teind1)))';
0167         else
0168             te(:,i)=NaN;
0169             dte(:,i)=NaN;
0170             ne(:,i)=NaN;
0171             dne(:,i)=NaN;
0172         end
0173         if(teflag)
0174             te0(i)=mdsdata('\results::te_x_a[[$1],13]',t0(i));
0175             te(:,i)=te(:,i)/max(te(:,i))*te0(i);
0176             dte(:,i)=dte(:,i)/max(te(:,i))*te0(i);
0177         end
0178         clear netemp teraw dteraw neraw dneraw tth teind1 teind2;
0179 
0180         % get Ti data if available and reasonable
0181         tiraw=mdsdata('\results::cxrs:ti');
0182         if ~isempty(tiraw),
0183             rhotiraw=mdsdata('\results::cxrs:rho');
0184             if isempty(rhotiraw),
0185                 rhotiraw=mdsdata('dim_of(\results::cxrs:ti,0)');
0186             end
0187             ttiraw=mdsdata('dim_of(\results::cxrs:ti,1)');
0188             tiind1=find(ttiraw>=t1(i)&ttiraw<=t2(i), 1 );
0189             tiind2=find(ttiraw>=t1(i)&ttiraw<=t2(i), 1, 'last' );
0190             if ~isempty(tiind1),
0191                 tiraw=nanmean(tiraw(:,tiind1:tiind2),2);
0192                 rhotiraw=nanmean(rhotiraw(:,tiind1:tiind2),2);
0193                 titemp=interp1(rhotiraw,tiraw,rho(:,i),'spline');
0194                 if(min(rhotiraw)>0.5)
0195                     warning('DNBI:DataProblem','DNBI penetration only up to rho=%g>0.5 => ignoring Ti measurement',min(rhotiraw));
0196                     clearti=true;
0197                 elseif(min(titemp)<0)
0198                     warning('DNBI:DataProblem','Ti from DNBI spline interpolation gives negative Ti values => ignoring Ti measurement');
0199                     clearti=true;
0200                 elseif(max(titemp)>2*max(tiraw))
0201                     warning('DNBI:DataProblem','Ti from DNBI spline interpolation goes up to Ti=%geV, this is twice the maximal underlying raw data Ti=%geV => ignoring Ti measurement',max(titemp),2*max(tiraw));
0202                     clearti=true;
0203                 else
0204                     clearti=false;
0205                 end
0206                 if(clearti)
0207                     ti(:,i)=NaN;
0208                     tiind1=[];
0209                     tiind2=[];
0210                     tiraw=[];
0211                 else
0212                     ti(:,i)=titemp;
0213                     clear('titemp');
0214                 end
0215             else
0216                 ti(:,i)=NaN;
0217             end
0218         else
0219             ti(:,i)=NaN;
0220         end
0221         clear rhotiraw tiraw;
0222 
0223         % get data for X2 and X3 launchers
0224         for j=laun
0225             clear ang
0226             js=int2str(j);
0227             if(j<7) % X2 launcher
0228                 ang.theta=nanmean(mdsdata(['\ecrh::launchers:theta_l:x2_',js,'[[$1]:[$2]:*]'],t1(i),t2(i)))/180*pi;
0229                 ang.phi=nanmean(mdsdata(['\ecrh::launchers:phi_l:x2_',js,'[[$1]:[$2]:*]'],t1(i),t2(i)))/180*pi;
0230                 ang.theta_eff=nanmean(mdsdata(['-atand(tand(\ecrh::launchers:theta_l:x2_',js, ...
0231                     '[[$1]:[$2]:*])*cosd(\ecrh::launchers:phi_l:x2_',js, ...
0232                     '[[$1]:[$2]:*]))'],t1(i),t2(i)))/180*pi;
0233                 dummy=mdsdata(['_y01=tand(\ecrh::launchers:theta_l:x2_',js,'[[$1]:[$2]:*])'],t1(i),t2(i));
0234                 dummy=mdsdata(['_y02=\ecrh::launchers:phi_l:x2_',js,'[[$1]:[$2]:*]'],t1(i),t2(i));
0235                 ang.phi_eff=nanmean(mdsdata(['atand(_y01*sind(_y02)/','sqrt(1+_y01^2*(cosd(_y02))^2))']))/180*pi;
0236                 ang.freq=82.7e9;
0237                 pow=mdsdata('_p=\results::toray.input:p_gyro[*,$1]',j)*1000;
0238                 tpow=mdsdata('dim_of(_p)');
0239                 i1=find(tpow>=t1(i), 1 );
0240                 if ~isempty(i1),
0241                     i2=find(tpow<=t2(i), 1, 'last' );
0242                     ang.P0=nanmean(pow(i1:i2));
0243                 else
0244                     ang.P0=NaN;
0245                 end
0246                 ang.P1_R=1.3072;
0247                 ang.P1_theta=52.013/180*pi;
0248                 ang.P2_R=1.1906;
0249                 ang.P2_d=0.01;
0250                 if(j==1 || j==4)
0251                     ang.P1_Z=-0.05865;
0252                     ang.P2_Z=0.0645;
0253                     ang.P3_Z=-0.0025;
0254                 else
0255                     ang.P1_Z=0.39885;
0256                     ang.P2_Z=0.522;
0257                     ang.P3_Z=0.455;
0258                 end
0259                 ang.P3_tor=0.;
0260             else % X3 launcher
0261                 thetaraw=nanmean(mdsdata(['\ecrh::launchers:theta_l:x3','[[$1]:[$2]:*]'],t1(i),t2(i)))/180*pi;
0262                 phiraw=mdsdata(['\ecrh::launchers:phi_l:x3_',js])/180*pi;
0263                 %  convert from X3 theta and phi to X2 theta and phi:
0264                 %    X3: tan(theta_l)=-y/x, tan(phi_l)=z/sqrt(x^2+y^2)
0265                 %       x radially inward, y upward, z toroidal clockwise from top
0266                 %    X2: tan(theta_l)=sqrt(y^2+z^2)/x, tan(phi_l)=z/y
0267                 %    x radially inward, y upward, z toroidal counterclockwise from top
0268                 ang.theta=atan2(sqrt(sin(thetaraw)^2+tan(phiraw)^2),cos(thetaraw));
0269                 ang.phi=atan2(tan(phiraw),sin(thetaraw));
0270                 ang.theta_eff=-thetaraw;
0271                 dummy=mdsdata(['_y01=cosd(\ecrh::launchers:theta_l:x3','[[$1]:[$2]:*])'],t1(i),t2(i));
0272                 dummy=mdsdata(['_y02=tand(\ecrh::launchers:phi_l:x3_',js,')']);
0273                 ang.phi_eff=nanmean(mdsdata('atan2(_y02,abs(_y01))'));
0274                 ang.freq=118.e9;
0275                 ang.P0=0.47e6;
0276                 ang.P1_R=nanmean(mdsdata('\ecrh::launchers:R:x3[[$1]:[$2]:*]',t1(i),t2(i)));
0277                 ang.P1_Z=.9885;
0278                 ang.P1_theta=0;
0279                 ang.P2_R=ang.P1_R;
0280                 ang.P2_Z=ang.P1_Z;
0281                 ang.P2_d=0;
0282                 ang.P3_Z=ang.P1_Z;
0283                 ang.P3_tor=(ang.P1_R+0.125)*tan(phiraw);
0284             end
0285             ang.id=int2str(shot(i));
0286             eval(['angles_',int2str(j),'(1,i)=ang;']);
0287         end
0288 
0289         % get Zeff
0290         if isnan(zeff(1).data),
0291             z=mdsdata('\results::zx[[$1]:[$2]:0.05]',t1(i),t2(i));
0292             if(~isempty(z) && ~iscell(z) && ~ischar(z))
0293                 zefff(i)=nanmean(z(:,1));
0294             else
0295                 zefff(i)=NaN;
0296             end
0297         else
0298             [dummy,jt1]=min(abs(zeff(i).time-t1(i)));
0299             [dummy,jt2]=min(abs(zeff(i).time-t1(i)));
0300             zefff(:,i)=nanmean(interp1(zeff(i).rho(:),zeff(i).data(:,jt1:jt2),rho(:,i)),2);
0301         end
0302 
0303         % get z position
0304         zaxis(:,i)=nanmean(mdsdata('\results::z_axis[[$1]:[$2]:0.05]',t1(i),t2(i)));
0305 
0306         mdsclose;
0307     end
0308     
0309     % arrange data for saving in the old version:
0310     %   clear what is not needed and transpose if required
0311     clear ang dummy i i1 i2 j js pow tpow z
0312     vloop=vloop';
0313     nel=nel';
0314     te=te';
0315     dte=dte';
0316     ne=ne';
0317     dne=dne';
0318     rho=rho';
0319     ti=ti';
0320     clear zeff;
0321     zeff=zefff';
0322     clear zefff
0323     zaxis=zaxis';
0324     ipsign=ipsign';
0325     btsign=btsign';
0326     Ip=Ip';
0327 end
0328 
0329 % save the data to a mat file
0330 save(filename,'-v7.3');
0331 %save(filename,'-v7.3','vloop','te','ne','rho','ti','zeff','zaxis','ipsign','btsign','Ip');
0332 fprintf('\n\nData saved to %s\n',filename);
0333 
0334 end