From ce59087e9f5957ccd470bf292c4a81920386c47f Mon Sep 17 00:00:00 2001
From: Olivier Sauter <olivier.sauter@epfl.ch>
Date: Tue, 16 Apr 2019 21:47:55 +0000
Subject: [PATCH] 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
---
crpptbx/AUG/aug_help_parameters.m | 4 +-
crpptbx/AUG/gdat_aug.m | 518 ++++++++++++++++++------------
2 files changed, 315 insertions(+), 207 deletions(-)
diff --git a/crpptbx/AUG/aug_help_parameters.m b/crpptbx/AUG/aug_help_parameters.m
index b38507f0..6fae2053 100644
--- a/crpptbx/AUG/aug_help_parameters.m
+++ b/crpptbx/AUG/aug_help_parameters.m
@@ -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. = '';
diff --git a/crpptbx/AUG/gdat_aug.m b/crpptbx/AUG/gdat_aug.m
index ddc2416c..111f53c0 100644
--- a/crpptbx/AUG/gdat_aug.m
+++ b/crpptbx/AUG/gdat_aug.m
@@ -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)
--
GitLab