From dc00b8e7d4344988556d161d45a80cc9c3fe7d37 Mon Sep 17 00:00:00 2001
From: Olivier Sauter <olivier.sauter@epfl.ch>
Date: Thu, 18 Feb 2016 10:44:21 +0000
Subject: [PATCH] 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 | 74 ++++++++++++++++++++++++++++++++++++++++++
1 file changed, 74 insertions(+)
diff --git a/crpptbx/AUG/gdat_aug.m b/crpptbx/AUG/gdat_aug.m
index 816b5244..b2c205af 100644
--- a/crpptbx/AUG/gdat_aug.m
+++ b/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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
--
GitLab