add fit to cxrs

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

crpptbx/AUG/gdat_aug.m

@@ -392,6 +392,12 @@ elseif strcmp(mapping_for_aug.method,'expression')
 elseif strcmp(mapping_for_aug.method,'switchcase')
   switch data_request_eff % not lower(...) since data_request_eff should be lower case already at this stage
    case {'cxrs', 'cxrs_rho'}
+    % if 'fit' option is added: 'fit',1, then the fitted profiles are returned which means "_rho" is effectively called
+    if ~isfield(gdat_data.gdat_params,'fit') || isempty(gdat_data.gdat_params.fit) || ~isnumeric(gdat_data.gdat_params.fit)
+      gdat_data.gdat_params.fit = 0;
+    elseif gdat_data.gdat_params.fit==1
+      data_request_eff = 'cxrs_rho';
+    end
     exp_name_eff = gdat_data.gdat_params.exp_name;
     if ~isfield(gdat_data.gdat_params,'source') || isempty(gdat_data.gdat_params.source) || strcmp(upper(gdat_data.gdat_params.source),'CEZ')
       gdat_data.gdat_params.source = 'CEZ';
@@ -507,6 +513,74 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
       gdat_data.rhopolnorm = rhopolnorm_out;
       gdat_data.rhotornorm = rhotornorm_out;
       gdat_data.rhovolnorm = rhovolnorm_out;
+      % defaults for fits, so user always gets std structure
+      gdat_data.fit.rhotornorm = []; % same for both ti and vrot
+      gdat_data.fit.rhopolnorm = [];
+      gdat_data.fit.t = [];
+      gdat_data.fit.ti.data = [];
+      gdat_data.fit.ti.drhotornorm = [];
+      gdat_data.fit.vrot.data = [];
+      gdat_data.fit.vrot.drhotornorm = [];
+      gdat_data.fit.raw.rhotornorm = [];
+      gdat_data.fit.raw.ti.data = [];
+      gdat_data.fit.raw.vrot.data = [];
+      fit_tension_default = -1.;
+      if isfield(gdat_data.gdat_params,'fit_tension')
+	fit_tension = gdat_data.gdat_params.fit_tension;
+      else
+	fit_tension = fit_tension_default;
+      end
+      if ~isstruct(fit_tension)
+	fit_tension_eff.ti = fit_tension;
+	fit_tension_eff.vrot = fit_tension;
+	fit_tension = fit_tension_eff;
+      else
+	if ~isfield(fit_tension,'ti'); fit_tension.ti = fit_tension_default; end
+	if ~isfield(fit_tension,'vrot '); fit_tension.vi = fit_tension_default; end
+      end
+      gdat_data.gdat_params.fit_tension = fit_tension;
+      %
+      if gdat_data.gdat_params.fit==1
+	% add fits
+	gdat_data.fit.raw.rhotornorm = NaN*ones(size(gdat_data.ti.data));
+	gdat_data.fit.raw.ti.data = NaN*ones(size(gdat_data.ti.data));
+	gdat_data.fit.raw.vrot.data = NaN*ones(size(gdat_data.vrot.data));
+	rhotornormfit = linspace(0,1,161)';
+	gdat_data.fit.rhotornorm = rhotornormfit;
+	for it=1:length(gdat_data.t)
+	  % make rhotor->rhopol transformation for each time since equilibrium might have changed
+	  irhook=find(gdat_data.rhotornorm(:,it)>0 & gdat_data.rhotornorm(:,it)<1); % no need for ~isnan
+	  [rhoeff isort]=sort(gdat_data.rhotornorm(irhook,it));
+	  gdat_data.fit.rhopolnorm(:,it)=interpos([0; rhoeff; 1],[0; gdat_data.rhopolnorm(irhook(isort),it); 1],rhotornormfit,-0.1,[2 2],[0 1]);
+	  idata = find(gdat_data.ti.data(:,it)>0 & gdat_data.rhotornorm(:,it)<1.01);
+	  if length(idata)>0
+	    gdat_data.fit.ti.raw.rhotornorm(idata,it) = gdat_data.rhotornorm(idata,it);
+	    gdat_data.fit.ti.raw.data(idata,it) = gdat_data.ti.data(idata,it);
+	    gdat_data.fit.ti.raw.error_bar(idata,it) = gdat_data.ti.error_bar(idata,it);
+	    gdat_data.fit.vrot.raw.rhotornorm(idata,it) = gdat_data.rhotornorm(idata,it);
+	    gdat_data.fit.vrot.raw.data(idata,it) = gdat_data.vrot.data(idata,it);
+	    gdat_data.fit.vrot.raw.error_bar(idata,it) = gdat_data.vrot.error_bar(idata,it);
+	    [rhoeff,irhoeff] = sort(gdat_data.rhotornorm(idata,it));
+	    rhoeff = [0; rhoeff];
+	    % they are some strange low error_bars, so remove these by max(Ti/error_bar)<=10; and changing it to large error bar
+	    tieff = gdat_data.ti.data(idata(irhoeff),it);
+	    ti_err_eff = gdat_data.ti.error_bar(idata(irhoeff),it);
+	    ij=find(tieff./ti_err_eff>10.);
+	    if ~isempty(ij); ti_err_eff(ij) = tieff(ij)./0.1; end
+	    vroteff = gdat_data.vrot.data(idata(irhoeff),it);
+	    vrot_err_eff = gdat_data.vrot.error_bar(idata(irhoeff),it);
+	    ij=find(vroteff./vrot_err_eff>10.);
+	    if ~isempty(ij); vrot_err_eff(ij) = vroteff(ij)./0.1; end
+	    %
+	    tieff =  [tieff(1); tieff];
+	    ti_err_eff =  [1e4; ti_err_eff];
+	    vroteff =  [vroteff(1); vroteff];
+	    vrot_err_eff =  [1e5; vrot_err_eff];
+	    [gdat_data.fit.ti.data(:,it), gdat_data.fit.ti.drhotornorm(:,it)] = interpos(rhoeff,tieff,rhotornormfit,fit_tension.ti,[1 0],[0 0],ti_err_eff);
+	    [gdat_data.fit.vrot.data(:,it), gdat_data.fit.vrot.drhotornorm(:,it)] = interpos(rhoeff,vroteff,rhotornormfit,fit_tension.vrot,[1 0],[0 0],vrot_err_eff);
+	  end
+	end
+      end
     end
 	
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%