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 & data_in{idata}.data(:,it_data{idata}(it))>0);
	  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 = [];
    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;
  if ~isempty(rhoeff)
    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
  else
    xxfit(1:nb_points,i) = NaN;
    dxxfitdrn(1:nb_points,i) = NaN;
  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};