diff --git a/crpptbx/AUG/aug_help_parameters.m b/crpptbx/AUG/aug_help_parameters.m
index b38507f0e6f8cd6b93c152454b469a81ee7f1cd1..6fae2053ce8e0ff79a220abff60b118b3c370092 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 ddc2416c9974ee29eaca49562c12bd578d662df7..111f53c0cbefb8564b4a6f0776bd5f5baf2dd1fb 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)