for generic gui as well

git-svn-id: https://spcsvn.epfl.ch/repos/TCV/gdat/trunk@5865 d63d8f72-b253-0410-a779-e742ad2e26cf

crpptbx/generic_fit.m

+function [gen_signals_data] = generic_fit(shot,data_in,dataname,time_array_in,tension_xx_in,coeff_xx_in,in_option);
+%
+%  [gen_signals_data,gen_signals_fit,filename_withfits] = generic_fit(shot,data_in,dataname,time_interval_in,tension_xx_in,coeff_xx_in,filename_withfits_in,in_option);
+%
+% Required inputs:
+%
+%    shot: shot number
+%    data_in{i}: array of structure from each "diagnostic" to combine, containing the following fields:
+%                        data_in{i}.data(:,it), data_in{i}.t, data_in{i}.error_bar(:,it)
+%                and a subfield .x for the radial coordinate with:
+%                        data_in{i}.x(irho,it) or data_in{i}.x(irho) if the rho values do not change with time
+%                and a subfield data_in{i}.xname which is used to describe if it is rhopol, rhotornorm, etc
+%    dataname: name of data and of subfield to fill in in the gen_signals_data structure:
+%
+%    It uses data_in{1}.t as the main time base thus as gen_signals_data.(dataname).t
+%
+%    in_option: 0 (default) computes the fits from the provided data
+%               1 (combined data already performed thus uses data_in{1}.combined_data or data_in.combined_data to compute the fits
+%
+%    output: gen_signals_data.(dataname).data, gen_signals_fit.(dataname).data
+%
+
+in_option_eff = 0;
+if exist('in_option') && ~isempty(in_option)
+  in_option_eff = in_option;
+end
+
+if in_option_eff == 0
+  % assumes receives just raw data in data_in{}
+  gen_signals_data.(dataname).raw_data = data_in;
+end
+if in_option_eff == 1
+  % assume one receives already gen_signals_data structure
+  if iscell(data_in)
+    gen_signals_data = data_in{1};
+  else
+    gen_signals_data = data_in;
+  end
+end
+
+if ~exist('coeff_xx_in') || length(coeff_xx_in)~=length(data_in)
+  coeff_xx_eff = ones(length(data_in),1);
+else
+  coeff_xx_eff = coeff_xx_in;
+end
+if exist('tension_xx_in') && ~isempty(tension_xx_in)
+  tension_xx_eff = tension_xx_in;
+else
+  tension_xx_eff=-3;
+end
+
+% combine the data to be able to do the fits
+if in_option_eff == 0
+
+  if isempty(data_in) || length(data_in)<1 || ~isfield(data_in{1},'t') || ~isfield(data_in{1},'data') || ~isfield(data_in{1},'x') || isempty(data_in{1}.data)
+    disp('no main data')
+    return
+  end
+  
+  main_time_base = data_in{1}.t;
+  
+  time_interval_in = [main_time_base(1) main_time_base(end)];
+  ij=find(main_time_base>=time_interval_in(1) & main_time_base<=time_interval_in(2));
+  
+  gen_signals_data.(dataname).combined_data.t = main_time_base(ij);
+  gen_signals_data.(dataname).t = gen_signals_data.(dataname).combined_data.t;
+  for i=1:length(ij)
+    if i==1
+      t1(i) = main_time_base(ij(i)) - 0.5*diff(main_time_base(ij([i:i+1])));
+    else
+      t1(i) = mean(main_time_base(ij([i-1:i])));
+    end
+    if i==length(ij)
+      t2(i) = main_time_base(ij(i)) + 0.5*diff(main_time_base(ij([i-1:i])));
+    else
+      t2(i) = mean(main_time_base(ij([i:i+1])));
+    end
+    gen_signals_data.(dataname).combined_data.time_interval(i,1:2) = [t1(i) t2(i)];
+  end
+  
+  for idata=1:length(data_in)
+    gen_signals_data.(dataname).combined_data.rholabel{idata} = data_in{idata}.xname;
+  end
+
+  for i=1:length(ij)
+    % for each time interval, collect the data
+    rhotor_data_tofit = []; 
+    xxdata_tofit = [];
+    xxerr_tofit = [];
+    data_provenance = [];
+    for idata=1:length(data_in)
+      it_data{idata} = find(data_in{idata}.t>=t1(i) & data_in{idata}.t<t2(i));
+      % construct 1D array with data from data_in{idata}
+      if ~isempty(it_data{idata})
+	for it=1:length(it_data{idata})
+	  jrho_data = find(data_in{idata}.x(:,it_data{idata}(it))>0 & data_in{idata}.x(:,it_data{idata}(it))<1.01);
+	  if length(jrho_data)>0
+	    rhotor_data_tofit(end+1:end+length(jrho_data)) = data_in{idata}.x(jrho_data,it_data{idata}(it));
+	    xxdata_tofit(end+1:end+length(jrho_data)) = data_in{idata}.data(jrho_data,it_data{idata}(it));
+	    xxerr_tofit(end+1:end+length(jrho_data)) = data_in{idata}.error_bar(jrho_data,it_data{idata}(it))./coeff_xx_eff(idata);
+	    data_provenance(end+1:end+length(jrho_data)) = idata.*ones(size(jrho_data));
+	  end
+	end
+      end
+    end
+    
+    if isempty(xxdata_tofit)
+      disp('xxdata_tofit empty')
+      gen_signals_data.(dataname).combined_data.perDt{i}.rho = [];
+      gen_signals_data.(dataname).combined_data.perDt{i}.data = [];
+      gen_signals_data.(dataname).combined_data.perDt{i}.error_bar = [];
+      gen_signals_data.(dataname).combined_data.perDt{i}.provenance = [];
+      return
+    else
+      [rhoeff,irhoeff] = sort(rhotor_data_tofit);
+      xxdata_tofit_eff = xxdata_tofit(irhoeff);
+      xxerr_tofit_eff = xxerr_tofit(irhoeff);
+      data_provenance_eff = data_provenance(irhoeff);
+      gen_signals_data.(dataname).combined_data.perDt{i}.rho = rhoeff;
+      gen_signals_data.(dataname).combined_data.perDt{i}.data = xxdata_tofit_eff;
+      gen_signals_data.(dataname).combined_data.perDt{i}.error_bar = xxerr_tofit_eff;
+      gen_signals_data.(dataname).combined_data.perDt{i}.provenance = data_provenance_eff;
+    end
+  end
+end
+
+% do the fit
+doplot=0; % for debugging
+nb_points = 201;
+main_rho = linspace(0,1,nb_points);
+xxfit = [];
+dxxfitdrn = [];
+for i=1:length(gen_signals_data.(dataname).combined_data.perDt)
+  rhoeff = gen_signals_data.(dataname).combined_data.perDt{i}.rho;
+  xxdata_tofit_eff = gen_signals_data.(dataname).combined_data.perDt{i}.data;
+  xxerr_tofit_eff = gen_signals_data.(dataname).combined_data.perDt{i}.error_bar;
+  xeff = [0. rhoeff];
+  yeffxx = [xxdata_tofit_eff(1) xxdata_tofit_eff];
+  erreffxx = [100.*max(xxerr_tofit_eff) xxerr_tofit_eff];
+  [xxfit(:,i),dxxfitdrn(:,i)]=interpos(xeff,yeffxx,main_rho,tension_xx_eff,[1 0],[0 0],erreffxx);
+  if doplot
+    figure(11);clf
+    errorbar(rhoeff,xxdata_tofit_eff,xxerr_tofit_eff,'*');
+    hold all
+    plot(main_rho,xxfit(:,i))
+    pause(0.01)
+  end
+end
+
+gen_signals_data.(dataname).fit.rhofit = main_rho;
+gen_signals_data.(dataname).fit.data = xxfit;
+gen_signals_data.(dataname).fit.dydrho = dxxfitdrn;
+gen_signals_data.(dataname).shot = shot;
+gen_signals_data.(dataname).fit.tension = tension_xx_eff;
+gen_signals_data.(dataname).fit.coeff_datas = coeff_xx_eff;
+gen_signals_data.(dataname).fit.rholabel = gen_signals_data.(dataname).combined_data.rholabel{1};
+
+

crpptbx/line_average.m

+function [nel_out,s_out,nes_out,varargout] = line_average(shot,time,Rdiag,Zdiag,ne_in,rhopol_in,Rmesh,Zmesh,psi_RZmesh,psi_axis,psi_edge);
+%
+% compute s and ne(s) in s_out and nes_out; where s is distance from first point into the plasma along Rdiag,Zdiag chord
+% compute int[ne(s)ds]/int[ds] in nel_out
+%
+% can get Rmesh, Zmesh and psi_RZmesh from eqdsk at shot,time from gdat if not given, but for many times, should use equil output from gdat
+%
+
+nel_out = [];
+s_out = [];
+nes_out = [];
+
+time_eff = time;
+
+if ~exist('Rmesh') || isempty(Rmesh)
+  % do only 1st time
+  time_eff = time(1);
+  ab=gdat(shot,'eqdsk','time',time_eff);
+  eqdsk=ab.eqdsk;
+  if isempty(Rdiag)
+    % assume H-1 line
+    par=sf2par('DCN',shot,'H1LFS_R','DCNgeo');
+    Rlfs=par.value/1e3;
+    par=sf2par('DCN',shot,'H1LFS_z','DCNgeo');
+    Zlfs=par.value/1e3;
+    par=sf2par('DCN',shot,'H1HFS_R','DCNgeo');
+    Rhfs=par.value/1e3;
+    par=sf2par('DCN',shot,'H1HFS_z','DCNgeo');
+    Zhfs=par.value/1e3;
+    Rdiag=[Rlfs Rhfs];
+    Zdiag=[Zlfs Zhfs];
+  end
+  if length(Rdiag)==2
+    nbpoints=50;
+    Rout=linspace(Rdiag(1),Rdiag(2),nbpoints);
+    Zout=linspace(Zdiag(1),Zdiag(2),nbpoints);
+  elseif length(Rdiag)>10
+    Rout=Rdiag;
+    Zout=Zdiag;
+  else
+    disp(['unexpected nb of diag points: length(Rdiag) = ',length(Rdiag)]);
+    return
+  end
+  Rmesh = eqdsk.rmesh;
+  Zmesh = eqdsk.zmesh;
+  psi_RZmesh(:,:,1) = eqdsk.psirz;
+  psi_axis(1) = eqdsk.psiaxis;
+  psi_edge(1) = eqdsk.psiedge;
+end
+
+if prod(size(Rmesh)) == max(size(Rmesh))
+  Rmesh_eff = repmat(reshape(Rmesh,length(Rmesh),1),1,length(time_eff));
+  Zmesh_eff = repmat(reshape(Zmesh,length(Zmesh),1),1,length(time_eff));
+end
+
+for it=1:length(time_eff)
+  psi_at_routzout(:,it) = interpos2Dcartesian(Rmesh_eff(:,it),Zmesh(:,it),psi_RZmesh(:,:,it),Rout,Zout);
+  ij_in=find((psi_at_routzout(:,it)-psi_axis(it)).*(psi_at_routzout(:,it)-psi_edge(it))<=0);
+  if ~isempty(ij_in)
+    if ij_in(1)==1
+      Rstart=Rout(ij_in(1));
+      Zstart=Zout(ij_in(1));
+    else
+      Rstart = (psi_edge(it)-psi_at_routzout(ij_in(1)-1),it)./diff(psi_at_routzout(ij_in(1)-1:ij_in(1),it)).*diff(Rout(ij_in(1)-1:ij_in(1))) + Rout(ij_in(1)-1);
+      Zstart = (psi_edge(it)-psi_at_routzout(ij_in(1)-1),it)./diff(psi_at_routzout(ij_in(1)-1:ij_in(1),it)).*diff(Zout(ij_in(1)-1:ij_in(1))) + Zout(ij_in(1)-1);
+    end
+    if ij_in(end)==length(Rout)
+      Rend=Rout(ij_in(end));
+      Zend=Zout(ij_in(end));
+    else
+      Rend = (psi_edge(it)-psi_at_routzout(ij_in(end),it))./diff(psi_at_routzout(ij_in(end):ij_in(end)+1,it)).*diff(Rout(ij_in(end):ij_in(end)+1)) + Rout(ij_in(end));
+      Zend = (psi_edge(it)-psi_at_routzout(ij_in(end),it))./diff(psi_at_routzout(ij_in(end):ij_in(end)+1,it)).*diff(Zout(ij_in(end):ij_in(end)+1)) + Zout(ij_in(end));
+    end
+    s_eff = sqrt((Rout(ij_in)-Rstart).^2 + (Zout(ij_in)-Zstart).^2);
+    rhopol_eff=sqrt((psi_at_routzout(ij_in,it)-psi_axis(it))./(psi_edge(it)-psi_axis(it)));
+    if ij_in(1)==1
+      rhopol_eff = reshape(rhopol_eff,length(rhopol_eff),1);
+      s_eff = reshape(s_eff,length(s_eff),1);
+    else
+      rhopol_eff=[1 ;reshape(rhopol_eff,length(rhopol_eff),1); 1];
+      s_eff=[0 ;reshape(s_eff,length(s_eff),1); sqrt((Rend-Rstart).^2 + (Zend-Zstart).^2)];
+    end
+    [rho_in_eff,isort]=sort(rhopol_in(:,it));
+    if rho_in_eff(1)<1e-3
+      % assume 1st point is center
+      nes = interpos(rho_in_eff,ne_in(isort,it),rhopol_eff,-0.1,[1 0],[0 0]);
+    else
+      rho_in_eff = [0 ; rho_in_eff];
+      ne_in_eff = [ne_in(isort(1),it); ne_in(isort,it)];
+      sigma_in = ones(size(ne_in_eff));
+      sigma_in(1)=100;
+      nes = interpos(rho_in_eff,ne_in_eff,rhopol_eff,-1,[1 0],[0 0],sigma_in);
+    end
+    nes_out(:,it) = nes;
+    s_out(:,it) = s_eff;
+    [dum1,~,~,nel_out_prof] = interpos(s_eff,nes,-0.1);
+    nel_out(it) = nel_out_prof(end)./(s_eff(end)-s_eff(1));
+  end
+end