Skip to content
Snippets Groups Projects
Commit ce59087e authored by Olivier Sauter's avatar Olivier Sauter
Browse files

add scrolling through cez, cuz, coz to find first non empty ti vrot 

git-svn-id: https://spcsvn.epfl.ch/repos/TCV/gdat/trunk@11763 d63d8f72-b253-0410-a779-e742ad2e26cf
parent 48517234
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
crpptbx/AUG/aug_help_parameters.m
+ 3
− 1
View file @ ce59087e
......@@ -39,7 +39,9 @@ help_struct_all.cocos = ['cocos value desired in output, uses eqdsk_cocos_transf
help_struct_all.fit = '0, no fit profiles, 1 (default) if fit profiles desired as well, relevant for _rho profiles. See also fit_type';
help_struct_all.fit_type = 'type of fits ''std'' (default) uses diagnostic error bars, ''pedestal'', uses manual error bars with smaller values outside 0.8';
help_struct_all.fit_nb_rho_points = 'nb of points for the radial mesh over which the fits are evaluated for the fitted profiles, it uses an equidistant mesh at this stage';
help_struct_all.source = 'sxr: ''G'' (default, with ssx), camera name ''J'', ''G'', ...[[F-M], case insensitive;; cxrs: ''CEZ'' (default), ''CMZ'';; raptor: ''observer'', ''predictive'' (or ''obs'', ''pre'') to restrict the ouput to these signals';
help_struct_all.source = ['sxr: ''G'' (default, with ssx), camera name ''J'', ''G'', ...[F-M], case insensitive;' char(10) ...
'cxrs: ''CEZ'' (default), ''CMZ'',''CUZ'',''COZ'',''all'';' char(10) ...
'raptor: ''observer'', ''predictive'' (or ''obs'', ''pre'') to restrict the ouput to these signals'];
help_struct_all.camera = ['[] (default, all), [i1 i2 ...] chord nbs ([1 3 5] if only chords 1, 3 and 5 are desired), ''central'' uses J_049'];
help_struct_all.freq = '''slow'', default (means ssx, 500kHz), lower sampling; ''fast'' full samples 2MHz; integer value nn for downsampling every nn''th points';
%help_struct_all. = '';
......
crpptbx/AUG/gdat_aug.m
+ 312
− 206
View file @ ce59087e
......@@ -463,231 +463,337 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
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;
if strcmp(data_request_eff,'cxrs')
gdat_data.gdat_params.fit = 0;
else
gdat_data.gdat_params.fit = 1; % add fit as default if cxrs_rho requested (since equil takes most time)
end
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';
r_node = 'R_time';
z_node = 'z_time';
elseif strcmp(upper(gdat_data.gdat_params.source),'CMZ')
gdat_data.gdat_params.source = 'CMZ';
r_node = 'R';
z_node = 'z';
elseif strcmp(upper(gdat_data.gdat_params.source),'CUZ')
gdat_data.gdat_params.source = 'CUZ';
r_node = 'R_time';
z_node = 'z_time';
elseif strcmp(upper(gdat_data.gdat_params.source),'COZ')
gdat_data.gdat_params.source = 'COZ';
r_node = 'R_time';
z_node = 'z_time';
sources_available = {'CEZ', 'CMZ', 'CUZ', 'COZ'};
r_node_available = {'R_time', 'R', 'R_time', 'R_time'};
z_node_available = {'z_time', 'z', 'z_time', 'z_time'};
% scroll through list up to first available when no sources provided
sources_available_for_scroll = {'CEZ', 'CUZ', 'COZ'};
scroll_through_list = 0;
if ~isfield(gdat_data.gdat_params,'source') || isempty(gdat_data.gdat_params.source)
gdat_data.gdat_params.source = sources_available_for_scroll(1);
scroll_through_list = 1; % means scroll through sources until one is available
elseif ischar(gdat_data.gdat_params.source)
gdat_data.gdat_params.source = {gdat_data.gdat_params.source};
end
if length(gdat_data.gdat_params.source)==1
if strcmp(upper(gdat_data.gdat_params.source{1}),'CEZ')
gdat_data.gdat_params.source = {'CEZ'};
elseif strcmp(upper(gdat_data.gdat_params.source{1}),'CMZ')
gdat_data.gdat_params.source = {'CMZ'};
elseif strcmp(upper(gdat_data.gdat_params.source{1}),'CUZ')
gdat_data.gdat_params.source = {'CUZ'};
elseif strcmp(upper(gdat_data.gdat_params.source{1}),'COZ')
gdat_data.gdat_params.source = {'COZ'};
elseif strcmp(upper(gdat_data.gdat_params.source{1}),'ALL')
gdat_data.gdat_params.source = sources_available;
scroll_through_list = 2; % means scroll through all sources and load all sources
else
warning(['source = ' gdat_data.gdat_params.source ' not expected with data_request= ' data_request_eff]);
return
end
else
warning(['source = ' gdat_data.gdat_params.source ' not expected with data_request= ' data_request_eff]);
return
sources_in = intersect(sources_available,upper(gdat_data.gdat_params.source));
if length(sources_in) ~= length(gdat_data.gdat_params.source)
disp('following sources not yet available, check with O. Sauter if need be')
setdiff(upper(gdat_data.gdat_params.source),sources_available)
end
gdat_data.gdat_params.source = sources_in;
end
diag_name = gdat_data.gdat_params.source;
extra_arg_sf2sig_eff_string = '';
if ~strcmp(gdat_data.gdat_params.extra_arg_sf2sig,'[]')
extra_arg_sf2sig_eff_string = [',' gdat_data.gdat_params.extra_arg_sf2sig];
end
% R, Z positions of measurements
try
% eval(['[r_time]=sf2ab(diag_name,shot,r_node,''-exp'',exp_name_eff' extra_arg_sf2sig_eff_string ');']);
% both for CEZ and CMZ, and.. Ti:1 is ok, otherwise introduce string above
[r_time,err]=rdaAUG_eff(shot,diag_name,r_node,exp_name_eff,[],extra_arg_sf2sig_eff_string,'area-base:Ti:1');
catch ME_R_time
% assume no shotfile
disp(getReport(ME_R_time))
return
end
gdat_data.r = r_time.data;
inotok=find(gdat_data.r<=0);
gdat_data.r(inotok) = NaN;
try
% eval(['[z_time]=sf2ab(diag_name,shot,z_node,''-exp'',exp_name_eff' extra_arg_sf2sig_eff_string ');']);
[z_time,err]=rdaAUG_eff(shot,diag_name,z_node,exp_name_eff,[],extra_arg_sf2sig_eff_string,'area-base:Ti:1');
gdat_data.z = z_time.data;
inotok=find(gdat_data.z<=0);
gdat_data.z(inotok) = NaN;
catch ME_R_time
disp(getReport(ME_R_time))
return
end
try
% eval(['[time]=sf2tb(diag_name,shot,''time'',''-exp'',exp_name_eff' extra_arg_sf2sig_eff_string ');']);
[time,err]=rdaAUG_eff(shot,diag_name,'time',exp_name_eff,[],extra_arg_sf2sig_eff_string,'time-base:Ti:0');
%[time,err]=rdaAUG_eff(shot,diag_name,'time',exp_name_eff,[],extra_arg_sf2sig_eff_string,'area-base:Ti:0');
gdat_data.t = time.data;
catch ME_R_time
disp(getReport(ME_R_time))
return
end
gdat_data.dim{1} = {gdat_data.r , gdat_data.z};
gdat_data.dimunits{1} = 'R, Z [m]';
gdat_data.dim{2} = gdat_data.t;
gdat_data.dimunits{2} = 't [s]';
gdat_data.x = gdat_data.dim{1};
% vrot
[a,error_status]=rdaAUG_eff(shot,diag_name,'vrot',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
a.name = 'vrot';
if isempty(a.data) || isempty(a.t) || error_status>0
if gdat_params.nverbose>=3;
a
disp(['with data_request= ' data_request_eff])
% set starting source
i_count = 1;
diag_name = gdat_data.gdat_params.source{i_count};
sources_tried{i_count} = diag_name;
iload = 1;
iequil = 0;
while iload==1
ishotfile_ok = 1;
i_source = strmatch(diag_name,sources_available);
r_node = r_node_available{i_source};
z_node = z_node_available{i_source};
% R, Z positions of measurements
try
% eval(['[r_time]=sf2ab(diag_name,shot,r_node,''-exp'',exp_name_eff' extra_arg_sf2sig_eff_string ');']);
% both for CEZ and CMZ, and.. Ti:1 is ok, otherwise introduce string above
[r_time,err]=rdaAUG_eff(shot,diag_name,r_node,exp_name_eff,[],extra_arg_sf2sig_eff_string,'area-base:Ti:1');
catch ME_R_time
% assume no shotfile
disp(getReport(ME_R_time))
ishotfile_ok = 0;
end
if ishotfile_ok == 1
gdat_data.r = r_time.data;
inotok=find(gdat_data.r<=0);
gdat_data.r(inotok) = NaN;
try
% eval(['[z_time]=sf2ab(diag_name,shot,z_node,''-exp'',exp_name_eff' extra_arg_sf2sig_eff_string ');']);
[z_time,err]=rdaAUG_eff(shot,diag_name,z_node,exp_name_eff,[],extra_arg_sf2sig_eff_string,'area-base:Ti:1');
gdat_data.z = z_time.data;
inotok=find(gdat_data.z<=0);
gdat_data.z(inotok) = NaN;
catch ME_R_time
disp(getReport(ME_R_time))
ishotfile_ok = 0;
end
else
gdat_data.r = [];
gdat_data.z = [];
end
end
gdat_data.vrot.data = a.data;
if any(strcmp(fieldnames(a),'units')); gdat_data.vrot.units=a.units; end
if any(strcmp(fieldnames(a),'name')); gdat_data.vrot.name=a.name; end
gdat_data.vrot.label = 'vrot_tor';
[aerr,e]=rdaAUG_eff(shot,diag_name,'err_vrot',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
gdat_data.vrot.error_bar = aerr.data;
% Ti
% [a,e]=rdaAUG_eff(shot,diag_name,'Ti',exp_name_eff);
% [aerr,e]=rdaAUG_eff(shot,diag_name,'err_Ti',exp_name_eff);
[a,e]=rdaAUG_eff(shot,diag_name,'Ti_c',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
a.name = 'Ti_c';
[aerr,e]=rdaAUG_eff(shot,diag_name,'err_Ti_c',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
gdat_data.ti.data = a.data;
gdat_data.data = a.data;
gdat_data.label = [gdat_data.label '/Ti'];
if any(strcmp(fieldnames(a),'units')); gdat_data.ti.units=a.units; end
if any(strcmp(fieldnames(a),'units')); gdat_data.units=a.units; end
if any(strcmp(fieldnames(a),'name')); gdat_data.ti.name=a.name; end
gdat_data.ti.label = 'Ti_c';
gdat_data.ti.error_bar = aerr.data;
gdat_data.error_bar = aerr.data;
gdat_data.data_fullpath=[exp_name_eff '/' diag_name '/' a.name ';vrot, Ti in data, {r;z} in .x'];
%
if strcmp(data_request_eff,'cxrs_rho')
params_equil = gdat_data.gdat_params;
params_equil.data_request = 'equil';
[equil,params_equil,error_status] = gdat_aug(shot,params_equil);
if error_status>0
if gdat_params.nverbose>=3; disp(['problems with ' params_equil.data_request]); end
return
if ishotfile_ok == 1
try
% eval(['[time]=sf2tb(diag_name,shot,''time'',''-exp'',exp_name_eff' extra_arg_sf2sig_eff_string ');']);
[time,err]=rdaAUG_eff(shot,diag_name,'time',exp_name_eff,[],extra_arg_sf2sig_eff_string,'time-base:Ti:0');
%[time,err]=rdaAUG_eff(shot,diag_name,'time',exp_name_eff,[],extra_arg_sf2sig_eff_string,'area-base:Ti:0');
gdat_data.t = time.data;
catch ME_R_time
disp(getReport(ME_R_time))
ishotfile_ok = 0;
end
else
gdat_data.t = [];
end
gdat_data.gdat_params.equil = params_equil.equil;
gdat_data.equil = equil;
inb_chord_cxrs=size(gdat_data.data,1);
inb_time_cxrs=size(gdat_data.data,2);
psi_out = NaN*ones(inb_chord_cxrs,inb_time_cxrs);
rhopolnorm_out = NaN*ones(inb_chord_cxrs,inb_time_cxrs);
rhotornorm_out = NaN*ones(inb_chord_cxrs,inb_time_cxrs);
rhovolnorm_out = NaN*ones(inb_chord_cxrs,inb_time_cxrs);
% constructs intervals within which a given equil is used: [time_equil(i),time_equil(i+1)]
time_equil=[min(gdat_data.t(1)-0.1,equil.t(1)-0.1) 0.5.*(equil.t(1:end-1)+equil.t(2:end)) max(equil.t(end)+0.1,gdat_data.t(end)+0.1)];
iok=find(~isnan(gdat_data.r(:,1)));
for itequil=1:length(time_equil)-1
rr=equil.Rmesh(:,itequil);
zz=equil.Zmesh(:,itequil);
psirz_in = equil.psi2D(:,:,itequil);
it_cxrs_inequil = find(gdat_data.t>time_equil(itequil) & gdat_data.t<=time_equil(itequil+1));
if ~isempty(it_cxrs_inequil)
if ~strcmp(upper(gdat_data.gdat_params.source),'CMZ')
rout=gdat_data.r(iok);
zout=gdat_data.z(iok);
gdat_data.dim{1} = {gdat_data.r , gdat_data.z};
gdat_data.dimunits{1} = 'R, Z [m]';
gdat_data.dim{2} = gdat_data.t;
gdat_data.dimunits{2} = 't [s]';
gdat_data.x = gdat_data.dim{1};
% vrot
if ishotfile_ok == 1
try
[a,error_status]=rdaAUG_eff(shot,diag_name,'vrot',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
if isempty(a.data) || isempty(a.t) || error_status>0
if gdat_params.nverbose>=3;
a
disp(['with data_request= ' data_request_eff])
end
end
[aerr,e]=rdaAUG_eff(shot,diag_name,'err_vrot',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
catch ME_vrot
disp(getReport(ME_vrot))
ishotfile_ok = 0;
end
gdat_data.vrot.data = a.data;
gdat_data.vrot.error_bar = aerr.data;
else
gdat_data.vrot.data = [];
gdat_data.vrot.error_bar = [];
end
a.name = 'vrot';
if any(strcmp(fieldnames(a),'units')); gdat_data.vrot.units=a.units; end
if any(strcmp(fieldnames(a),'name')); gdat_data.vrot.name=a.name; end
gdat_data.vrot.label = 'vrot_tor';
% Ti
% [a,e]=rdaAUG_eff(shot,diag_name,'Ti',exp_name_eff);
% [aerr,e]=rdaAUG_eff(shot,diag_name,'err_Ti',exp_name_eff);
if ishotfile_ok == 1
try
[a,e]=rdaAUG_eff(shot,diag_name,'Ti_c',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
[aerr,e]=rdaAUG_eff(shot,diag_name,'err_Ti_c',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
catch ME_ti
disp(getReport(ME_ti))
ishotfile_ok = 0;
end
gdat_data.ti.data = a.data;
gdat_data.ti.error_bar = aerr.data;
else
gdat_data.ti.data = [];
gdat_data.ti.error_bar = [];
end
a.name = 'Ti_c';
if any(strcmp(fieldnames(a),'units')); gdat_data.ti.units=a.units; end
if any(strcmp(fieldnames(gdat_data.ti),'units')); gdat_data.units=gdat_data.ti.units; end
if any(strcmp(fieldnames(a),'name')); gdat_data.ti.name=a.name; end
gdat_data.ti.label = 'Ti_c';
% main node ti a this stage
gdat_data.data = gdat_data.ti.data;
gdat_data.label = [gdat_data.label '/Ti'];
gdat_data.error_bar = gdat_data.ti.error_bar;
gdat_data.data_fullpath=[exp_name_eff '/' diag_name '/' a.name ';vrot, Ti in data, {r;z} in .x'];
%
if strcmp(data_request_eff,'cxrs_rho')
% defaults
if iequil == 0
gdat_data.equil.data = [];
gdat_data.psi = [];
gdat_data.rhopolnorm = [];
gdat_data.rhotornorm = [];
gdat_data.rhovolnorm = [];
% 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
rout=gdat_data.r(iok,it_cxrs_inequil);
zout=gdat_data.z(iok,it_cxrs_inequil);
fit_tension = fit_tension_default;
end
psi_at_routzout = interpos2Dcartesian(rr,zz,psirz_in,rout,zout);
if ~strcmp(upper(gdat_data.gdat_params.source),'CMZ')
psi_out(iok,it_cxrs_inequil) = repmat(psi_at_routzout,[1,length(it_cxrs_inequil)]);
if ~isstruct(fit_tension)
fit_tension_eff.ti = fit_tension;
fit_tension_eff.vrot = fit_tension;
fit_tension = fit_tension_eff;
else
psi_out(iok,it_cxrs_inequil) = reshape(psi_at_routzout,length(iok),length(it_cxrs_inequil));
if ~isfield(fit_tension,'ti'); fit_tension.ti = fit_tension_default; end
if ~isfield(fit_tension,'vrot '); fit_tension.vrot = fit_tension_default; end
end
rhopolnorm_out(iok,it_cxrs_inequil) = sqrt((psi_out(iok,it_cxrs_inequil)-equil.psi_axis(itequil))./(equil.psi_lcfs(itequil)-equil.psi_axis(itequil)));
for it_cx=1:length(it_cxrs_inequil)
rhotornorm_out(iok,it_cxrs_inequil(it_cx)) = interpos(equil.rhopolnorm(:,itequil),equil.rhotornorm(:,itequil),rhopolnorm_out(iok,it_cxrs_inequil(it_cx)),-3,[2 2],[0 1]);
rhovolnorm_out(iok,it_cxrs_inequil(it_cx)) = interpos(equil.rhopolnorm(:,itequil),equil.rhovolnorm(:,itequil),rhopolnorm_out(iok,it_cxrs_inequil(it_cx)),-3,[2 2],[0 1]);
gdat_data.gdat_params.fit_tension = fit_tension;
if isfield(gdat_data.gdat_params,'fit_nb_rho_points')
fit_nb_rho_points = gdat_data.gdat_params.fit_nb_rho_points;
else
fit_nb_rho_points = 201;
end
gdat_data.gdat_params.fit_nb_rho_points = fit_nb_rho_points;
end
if ishotfile_ok == 1 && iequil == 0
params_equil = gdat_data.gdat_params;
params_equil.data_request = 'equil';
[equil,params_equil,error_status] = gdat_aug(shot,params_equil);
if error_status>0
if gdat_params.nverbose>=3; disp(['problems with ' params_equil.data_request]); end
return
end
iequil = 1;
gdat_data.gdat_params.equil = params_equil.equil;
gdat_data.equil = equil;
inb_chord_cxrs=size(gdat_data.data,1);
inb_time_cxrs=size(gdat_data.data,2);
psi_out = NaN*ones(inb_chord_cxrs,inb_time_cxrs);
rhopolnorm_out = NaN*ones(inb_chord_cxrs,inb_time_cxrs);
rhotornorm_out = NaN*ones(inb_chord_cxrs,inb_time_cxrs);
rhovolnorm_out = NaN*ones(inb_chord_cxrs,inb_time_cxrs);
% constructs intervals within which a given equil is used: [time_equil(i),time_equil(i+1)]
time_equil=[min(gdat_data.t(1)-0.1,equil.t(1)-0.1) 0.5.*(equil.t(1:end-1)+equil.t(2:end)) max(equil.t(end)+0.1,gdat_data.t(end)+0.1)];
iok=find(~isnan(gdat_data.r(:,1)));
for itequil=1:length(time_equil)-1
rr=equil.Rmesh(:,itequil);
zz=equil.Zmesh(:,itequil);
psirz_in = equil.psi2D(:,:,itequil);
it_cxrs_inequil = find(gdat_data.t>time_equil(itequil) & gdat_data.t<=time_equil(itequil+1));
if ~isempty(it_cxrs_inequil)
if ~strcmp(upper(gdat_data.gdat_params.source),'CMZ')
rout=gdat_data.r(iok);
zout=gdat_data.z(iok);
else
rout=gdat_data.r(iok,it_cxrs_inequil);
zout=gdat_data.z(iok,it_cxrs_inequil);
end
psi_at_routzout = interpos2Dcartesian(rr,zz,psirz_in,rout,zout);
if ~strcmp(upper(gdat_data.gdat_params.source),'CMZ')
psi_out(iok,it_cxrs_inequil) = repmat(psi_at_routzout,[1,length(it_cxrs_inequil)]);
else
psi_out(iok,it_cxrs_inequil) = reshape(psi_at_routzout,length(iok),length(it_cxrs_inequil));
end
rhopolnorm_out(iok,it_cxrs_inequil) = sqrt((psi_out(iok,it_cxrs_inequil)-equil.psi_axis(itequil))./(equil.psi_lcfs(itequil)-equil.psi_axis(itequil)));
for it_cx=1:length(it_cxrs_inequil)
rhotornorm_out(iok,it_cxrs_inequil(it_cx)) = interpos(equil.rhopolnorm(:,itequil),equil.rhotornorm(:,itequil),rhopolnorm_out(iok,it_cxrs_inequil(it_cx)),-3,[2 2],[0 1]);
rhovolnorm_out(iok,it_cxrs_inequil(it_cx)) = interpos(equil.rhopolnorm(:,itequil),equil.rhovolnorm(:,itequil),rhopolnorm_out(iok,it_cxrs_inequil(it_cx)),-3,[2 2],[0 1]);
end
end
end
gdat_data.psi = psi_out;
gdat_data.rhopolnorm = rhopolnorm_out;
gdat_data.rhotornorm = rhotornorm_out;
gdat_data.rhovolnorm = rhovolnorm_out;
%
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,fit_nb_rho_points)';
gdat_data.fit.rhotornorm = rhotornormfit;
gdat_data.fit.t = gdat_data.t;
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
end
gdat_data.psi = psi_out;
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.vrot = fit_tension_default; end
end
gdat_data.gdat_params.fit_tension = fit_tension;
if isfield(gdat_data.gdat_params,'fit_nb_rho_points')
fit_nb_rho_points = gdat_data.gdat_params.fit_nb_rho_points;
else
fit_nb_rho_points = 201;
end
gdat_data.gdat_params.fit_nb_rho_points = fit_nb_rho_points;
%
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,fit_nb_rho_points)';
gdat_data.fit.rhotornorm = rhotornormfit;
gdat_data.fit.t = gdat_data.t;
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);
if scroll_through_list == 0 || scroll_through_list == 2
if scroll_through_list == 2
tmp.(diag_name) = gdat_data;
end
if length(gdat_data.gdat_params.source) > i_count
i_count = i_count + 1;
diag_name = gdat_data.gdat_params.source{i_count};
else
iload = 0;
end
elseif scroll_through_list == 1
if ishotfile_ok == 1
iload = 0;
else
sources_remaining = setdiff(sources_available_for_scroll,sources_tried,'stable');
if ~isempty(sources_remaining)
i_count = i_count + 1;
diag_name = sources_remaining{1};
sources_tried{i_count} = diag_name;
else
iload = 0;
end
end
else
disp('should not arrive here, check value of scroll_through_list')
scroll_through_list
iload = 0;
end
end
if scroll_through_list == 2
tmp_field=fieldnames(tmp);
for i=1:length(tmp_field)
gdat_data.(tmp_field{i}) = tmp.(tmp_field{i});
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'ece', 'eced', 'ece_rho', 'eced_rho'}
nth_points = 13;
......@@ -1132,8 +1238,8 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
gdat_data.(ids_top_name) = ids_top;
gdat_data.([ids_top_name '_description']) = ids_top_description;
else
gdat_data.(ids_top_name) = equil_empty;
gdat_data.([ids_top_name '_description']) = ['shot empty so return default empty structure for ' ids_top_name];
gdat_data.(ids_top_name) = equil_empty;
gdat_data.([ids_top_name '_description']) = ['shot empty so return default empty structure for ' ids_top_name];
end
catch ME_aug_get_ids
getReport(ME_aug_get_ids)
......@@ -1746,8 +1852,8 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
try
ic=gdat_aug(shot,params_eff);
catch
ic.data = [];
ic.dim = [];
ic.data = [];
ic.dim = [];
end
if ~isempty(ic.data) && ~isempty(ic.dim)
for i=1:length(fields_to_copy)
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment