diff --git a/matlab/gdat.m b/matlab/gdat.m
index 11efb5931b88064d3e98894b1767e7922d8eb0d0..ae16099ae4d47f47bce0c83fd8674ca907c85825 100644
--- a/matlab/gdat.m
+++ b/matlab/gdat.m
@@ -187,10 +187,21 @@ try
   else
     args = [{shot,data_request},varargin_eff];
   end
-  [gdat_data,gdat_params,error_status,varargout] = feval(['gdat_' lower(machine_eff)],args{:});
-  % because data request can be an actual detailed tree related signal, upper and lower case need to be kept, but other places remain with lower case when case insensitive
-  % needed since some substructure have machine name like mapping_for
-  gdat_data.gdat_params.machine = lower(gdat_data.gdat_params.machine);
+  % treat multiple shot numbers
+  if exist('shot') && numel(args) > 0 && isnumeric(args{1})
+    args_in = args;
+    for i=1:numel(shot)
+      args{1} = shot(i);
+      [gdat_data_i,gdat_params,error_status(i),varargout] = feval(['gdat_' lower(machine_eff)],args{:});
+      gdat_data_i.gdat_params.machine = lower(gdat_data_i.gdat_params.machine);
+      gdat_data(i) = gdat_data_i;
+    end
+  else
+    [gdat_data,gdat_params,error_status,varargout] = feval(['gdat_' lower(machine_eff)],args{:});
+    % because data request can be an actual detailed tree related signal, upper and lower case need to be kept, but other places remain with lower case when case insensitive
+    % needed since some substructure have machine name like mapping_for
+    gdat_data.gdat_params.machine = lower(gdat_data.gdat_params.machine);
+  end
 
 catch ME_gdat
   warning(['problems calling gdat_' lower(machine_eff)]);
@@ -204,9 +215,13 @@ catch ME_gdat
   return
 end
 
-gdat_data.gdat_call = gdat_call;
+if numel(gdat_data) > 0
+  for i=1:numel(gdat_data)
+    gdat_data(i).gdat_call = gdat_call;
+  end
+end
 
-if gdat_data.gdat_params.doplot
+if numel(gdat_data) > 0 && gdat_data(1).gdat_params.doplot
   try
     % plot gdat_data versus 1st dim by default, if nb_dims<=2, otherwise do not plot
     hhh = gdat_plot(gdat_data); % return handle to figure
diff --git a/matlab/gdat_plot.m b/matlab/gdat_plot.m
index ccbf80d03b49bead91fae55fdf1432a7e4941c77..29895699b09e889c7085e485ec014d4c66f645f3 100644
--- a/matlab/gdat_plot.m
+++ b/matlab/gdat_plot.m
@@ -31,180 +31,189 @@ end
 doplot=1;
 if isfield(gdat_plot_params,'doplot') && ~isempty(gdat_plot_params.doplot)
   doplot = gdat_plot_params.doplot;
-elseif isfield(gdat_data,'gdat_params') && isfield(gdat_data.gdat_params,'doplot') && ~isempty(gdat_data.gdat_params.doplot)
-  if gdat_data.gdat_params.doplot~=0 % assume one does not call gdat_plot not to plot
-    doplot = gdat_data.gdat_params.doplot;
+elseif isfield(gdat_data(1),'gdat_params') && isfield(gdat_data(1).gdat_params,'doplot') && ~isempty(gdat_data(1).gdat_params.doplot)
+  if gdat_data(1).gdat_params.doplot~=0 % assume one does not call gdat_plot not to plot
+    doplot = gdat_data(1).gdat_params.doplot;
   end
 end
 if doplot==0; return; end
 redo_legend_from_Tags = 0;
 do_legend = 1;
 
-if all(isfield(gdat_data,{'data','t'})) && ~isempty(gdat_data.data) && ~isempty(gdat_data.t)
-  fighandle = get(0,'CurrentFigure');
-  if doplot == 1
-    fighandle = figure;
-  elseif doplot > 1
-    fighandle = figure(doplot);
-    clf;
-  elseif doplot == -1
-    hold all
-  elseif doplot < -1
-    fighandle = figure(abs(doplot));
-    hold all
-  end
-  if isfield(gdat_data,'gdat_request') && ischar(gdat_data.gdat_request) ...
-        && strcmp(gdat_data.gdat_request,'eqdsk')
-    % special case, plot contours of first equil
-    endstr = '';
-    if iscell(gdat_data.eqdsk); endstr = '{1}'; end
-    eval(['contour(gdat_data.eqdsk' endstr '.rmesh,gdat_data.eqdsk' endstr '.zmesh,gdat_data.eqdsk' endstr '.psi'',100);'])
-    hold on
-    eval(['linehandles{end+1} = plot(gdat_data.eqdsk' endstr '.rplas,gdat_data.eqdsk' endstr '.zplas,''k'');'])
-    eval(['linehandles{end+1} = plot(gdat_data.eqdsk' endstr '.rlim,gdat_data.eqdsk' endstr '.zlim,''k'');'])
-    axis equal;
-    title(eval(['gdat_data.eqdsk' endstr '.stitle']));
-    do_legend = 0;
-  elseif any(find(size(gdat_data.data)==length(gdat_data.t)))
-    try
-      if length(size(gdat_data.data))<=2
-        if isnumeric(gdat_data.t)
-          linehandles{end+1} = plot(gdat_data.t,gdat_data.data);
-        elseif iscell(gdat_data.t)
-          linehandles{end+1} = plot(str2num(cell2mat(gdat_data.t)),gdat_data.data);
+nb_linehandles_current = numel(linehandles);
+for ishot=1:numel(gdat_data)
+  nb_linehandles_prev = nb_linehandles_current;
+  if all(isfield(gdat_data(ishot),{'data','t'})) && ~isempty(gdat_data(ishot).data) && ~isempty(gdat_data(ishot).t)
+    fighandle = get(0,'CurrentFigure');
+    if ishot > 1 && gdat_data(1).gdat_params.doplot==1
+      % assume that for a series of shot, need doplot=-1 after 1st plot if doplot was 1 in request
+      doplot = -1;
+    end
+    if doplot == 1
+      fighandle = figure;
+    elseif doplot > 1
+      fighandle = figure(doplot);
+      clf;
+    elseif doplot == -1
+      hold all
+    elseif doplot < -1
+      fighandle = figure(abs(doplot));
+      hold all
+    end
+    if isfield(gdat_data(ishot),'gdat_request') && ischar(gdat_data(ishot).gdat_request) ...
+        && strcmp(gdat_data(ishot).gdat_request,'eqdsk')
+      % special case, plot contours of first equil
+      endstr = '';
+      if iscell(gdat_data(ishot).eqdsk); endstr = '{1}'; end
+      eval(['contour(gdat_data(ishot).eqdsk' endstr '.rmesh,gdat_data(ishot).eqdsk' endstr '.zmesh,gdat_data(ishot).eqdsk' endstr '.psi'',100);'])
+      hold on
+      eval(['linehandles{end+1} = plot(gdat_data(ishot).eqdsk' endstr '.rplas,gdat_data(ishot).eqdsk' endstr '.zplas,''k'');'])
+      eval(['linehandles{end+1} = plot(gdat_data(ishot).eqdsk' endstr '.rlim,gdat_data(ishot).eqdsk' endstr '.zlim,''k'');'])
+      axis equal;
+      title(eval(['gdat_data(ishot).eqdsk' endstr '.stitle']));
+      do_legend = 0;
+    elseif any(find(size(gdat_data(ishot).data)==length(gdat_data(ishot).t)))
+      try
+        if length(size(gdat_data(ishot).data))<=2
+          if isnumeric(gdat_data(ishot).t)
+            linehandles{end+1} = plot(gdat_data(ishot).t,gdat_data(ishot).data);
+          elseif iscell(gdat_data(ishot).t)
+            linehandles{end+1} = plot(str2num(cell2mat(gdat_data(ishot).t)),gdat_data(ishot).data);
+          end
+        else
+          disp('WARNING')
+          disp(['size(gdat_data(ishot).data) = ' num2str(size(gdat_data(ishot).data)) ', >2D thus do not try to plot'])
+          disp(' ')
+          return
+        end
+      catch ME
+        if exist('ME','var')
+          disp('Problem in gdat_plot')
+          rethrow(ME)
         end
-      else
-        disp('WARNING')
-        disp(['size(gdat_data.data) = ' num2str(size(gdat_data.data)) ', >2D thus do not try to plot'])
-        disp(' ')
         return
       end
-    catch ME
-      if exist('ME','var')
-        disp('Problem in gdat_plot')
-        rethrow(ME)
+      if ~isfield(gdat_data(ishot),'shot'); return; end % allows to plot if just .t and .data exist
+      if doplot < 0
+        title([gdat_data(ishot).gdat_params.machine]);
+      else
+        title([gdat_data(ishot).gdat_params.machine ' #' num2str(gdat_data(ishot).shot)]);
       end
-      return
-    end
-    if ~isfield(gdat_data,'shot'); return; end % allows to plot if just .t and .data exist
-    if doplot < 0
-      title([gdat_data.gdat_params.machine]);
-    else
-      title([gdat_data.gdat_params.machine ' #' num2str(gdat_data.shot)]);
-    end
-    if isfield(gdat_data,'mapping_for')
-      dimunits_x = gdat_data.dimunits{gdat_data.mapping_for.(gdat_data.gdat_params.machine).gdat_timedim};
-      if ischar(dimunits_x)
-        xlabel(['time [' dimunits_x ']']);
+      if isfield(gdat_data(ishot),'mapping_for')
+        dimunits_x = gdat_data(ishot).dimunits{gdat_data(ishot).mapping_for.(gdat_data(ishot).gdat_params.machine).gdat_timedim};
+        if ischar(dimunits_x)
+          xlabel(['time [' dimunits_x ']']);
+        else
+          xlabel(['time']);
+        end
       else
         xlabel(['time']);
       end
-    else
-      xlabel(['time']);
-    end
-    ylabel_eff = gdat_data.label;
-    if iscell(gdat_data.label) && length(gdat_data.label)>=2; ylabel_eff = gdat_data.label{end}; end
-    if ~isempty(gdat_data.units)
-      hylabel=ylabel([ylabel_eff '[' gdat_data.units ']'],'interpreter','none');
-    else
-      hylabel=ylabel(ylabel_eff,'interpreter','none');
-    end
-    ab=get(gca,'children');
-    if iscell(gdat_data.label) && numel(gdat_data.label) > 1;
-      if numel(ab)==numel(gdat_data.label) || mod(numel(ab),numel(gdat_data.label))==0
-        % Assumes a single shot with several lines and labels, if mod==0 then this is a new shot, only Tag present lines
-        for iab=1:numel(gdat_data.label)
-          set(ab(iab),'Tag',[num2str(gdat_data.shot) ' ' gdat_data.label{end-iab+1}]);
+      ylabel_eff = gdat_data(ishot).label;
+      if iscell(gdat_data(ishot).label) && length(gdat_data(ishot).label)>=2; ylabel_eff = gdat_data(ishot).label{end}; end
+      if ~isempty(gdat_data(ishot).units)
+        hylabel=ylabel([ylabel_eff '[' gdat_data(ishot).units ']'],'interpreter','none');
+      else
+        hylabel=ylabel(ylabel_eff,'interpreter','none');
+      end
+      ab=get(gca,'children');
+      if iscell(gdat_data(ishot).label) && numel(gdat_data(ishot).label) > 1;
+        if numel(ab)==numel(gdat_data(ishot).label) || mod(numel(ab),numel(gdat_data(ishot).label))==0
+          % Assumes a single shot with several lines and labels, if mod==0 then this is a new shot, only Tag present lines
+          for iab=1:numel(gdat_data(ishot).label)
+            set(ab(iab),'Tag',[num2str(gdat_data(ishot).shot) ' ' gdat_data(ishot).label{end-iab+1}]);
+          end
+          if numel(ab)>numel(gdat_data(ishot).label), redo_legend_from_Tags = 1; end
+        end
+        if iscell(gdat_data(ishot).label) && length(ab) < length(gdat_data(ishot).label)
+          if length(ab) == 1
+            % assume combine label is best
+            tempaaa = sprintf('%s/',gdat_data(ishot).label{:});
+            hhhleg=legend(tempaaa(1:end-1));
+            set(ab(1),'Tag',[num2str(gdat_data(ishot).shot) ' ' tempaaa(1:min(10,numel(tempaaa)-1))]);
+          else
+            % not sure why would get there
+            hhhleg=legend(gdat_data(ishot).label{1:length(ab)});
+          end
+        elseif numel(ab)==numel(gdat_data(ishot).label)
+          hhhleg=legend(gdat_data(ishot).label);
+        end
+        if exist('hhhleg','var'), set(hhhleg,'Interpreter','none'); end
+      else
+        % assume one line per call
+        if isempty(gdat_data(ishot).label)
+          set(ab(1),'Tag',[num2str(gdat_data(ishot).shot)]);
+        elseif ischar(gdat_data(ishot).label)
+          % assume one signal, take max 10 chars
+          set(ab(1),'Tag',[num2str(gdat_data(ishot).shot) ' ' gdat_data(ishot).label(1:min(10,numel(gdat_data(ishot).label)))]);
+        elseif iscell(gdat_data(ishot).label) && numel(gdat_data(ishot).label) == 1
+          if isempty(gdat_data(ishot).label{1})
+            set(ab(1),'Tag',[num2str(gdat_data(ishot).shot)]);
+          else
+            set(ab(1),'Tag',[num2str(gdat_data(ishot).shot) ' ' gdat_data(ishot).label{1}(1:min(10,numel(gdat_data(ishot).label{1})))]);
+          end
         end
-        if numel(ab)>numel(gdat_data.label), redo_legend_from_Tags = 1; end
       end
-      if iscell(gdat_data.label) && length(ab) < length(gdat_data.label)
-        if length(ab) == 1
-          % assume combine label is best
-          tempaaa = sprintf('%s/',gdat_data.label{:});
-          hhhleg=legend(tempaaa(1:end-1));
-          set(ab(1),'Tag',[num2str(gdat_data.shot) ' ' tempaaa(1:min(10,numel(tempaaa)-1))]);
-        else
-          % not sure why would get there
-          hhhleg=legend(gdat_data.label{1:length(ab)});
+      if redo_legend_from_Tags
+        for iab=1:numel(ab)
+          % Use Tag for DisplayName, when overlay plots of multiple data at this stage
+          set(ab(iab),'DisplayName',strrep(get(ab(iab),'Tag'),'_','\_'));
         end
-      elseif numel(ab)==numel(gdat_data.label)
-        hhhleg=legend(gdat_data.label);
       end
-      if exist('hhhleg','var'), set(hhhleg,'Interpreter','none'); end
+      zoom on;
+    end
+    maxnblines = 1;
+    nb_linehandles_current = numel(linehandles);
+    if ~exist('ab','var'), ab=get(gca,'children'); end
+    if do_legend==0 || redo_legend_from_Tags || any(strcmp(gdat_data(ishot).gdat_params.data_request,'powers'))  ...
+          || (numel(ab)==numel(gdat_data(ishot).label) && numel(ab)>1)
+      % keep legend as is
     else
-      % assume one line per call
-      if isempty(gdat_data.label)
-        set(ab(1),'Tag',[num2str(gdat_data.shot)]);
-      elseif ischar(gdat_data.label)
-        % assume one signal, take max 10 chars
-        set(ab(1),'Tag',[num2str(gdat_data.shot) ' ' gdat_data.label(1:min(10,numel(gdat_data.label)))]);
-      elseif iscell(gdat_data.label) && numel(gdat_data.label) == 1
-        if isempty(gdat_data.label{1})
-          set(ab(1),'Tag',[num2str(gdat_data.shot)]);
+      for i=nb_linehandles_prev+1:nb_linehandles_current
+        maxnblines = max(maxnblines,numel(linehandles{i}));
+        if numel(linehandles{i}) == 1
+          set(linehandles{i},'DisplayName',[num2str(gdat_data(ishot).shot)]);
         else
-          set(ab(1),'Tag',[num2str(gdat_data.shot) ' ' gdat_data.label{1}(1:min(10,numel(gdat_data.label{1})))]);
+          for j=1:numel(linehandles{i})
+            set(linehandles{i}(j),'DisplayName',[num2str(gdat_data(ishot).shot) ',' num2str(j)]);
+          end
         end
       end
     end
-    if redo_legend_from_Tags
-      for iab=1:numel(ab)
-        % Use Tag for DisplayName, when overlay plots of multiple data at this stage
-        set(ab(iab),'DisplayName',strrep(get(ab(iab),'Tag'),'_','\_'));
-      end
-    end
-    zoom on;
-  end
-  maxnblines = 1;
-  if ~exist('ab','var'), ab=get(gca,'children'); end
-  if do_legend==0 || redo_legend_from_Tags || any(strcmp(gdat_data.gdat_params.data_request,'powers'))  ...
-      || (numel(ab)==numel(gdat_data.label) && numel(ab)>1)
-    % keep legend as is
-  else
-    for i=1:numel(linehandles)
-      maxnblines = max(maxnblines,numel(linehandles{i}));
-      if numel(linehandles{i}) == 1
-        set(linehandles{i},'DisplayName',[num2str(gdat_data.shot)]);
+    if do_legend==1 && maxnblines==1; legend show; end
+    if strcmp(gdat_data(ishot).gdat_request,'mhd')
+      % special case, add legend instead of ylabel
+      delete(hylabel);
+      legend(gdat_data(ishot).label);
+      % add spectrogram per signal
+      mhd_sum_data = 0.;
+      if isfield(gdat_data(ishot).gdat_params,'nfft') && ~isempty(gdat_data(ishot).gdat_params.nfft)
+        nfft = gdat_data(ishot).gdat_params.nfft;
       else
-        for j=1:numel(linehandles{i})
-          set(linehandles{i}(j),'DisplayName',[num2str(gdat_data.shot) ',' num2str(j)]);
-        end
+        nfft=512;
       end
-    end
-  end
-  if do_legend==1 && maxnblines==1; legend show; end
-  if strcmp(gdat_data.gdat_request,'mhd')
-    % special case, add legend instead of ylabel
-    delete(hylabel);
-    legend(gdat_data.label);
-    % add spectrogram per signal
-    mhd_sum_data = 0.;
-    if isfield(gdat_data.gdat_params,'nfft') && ~isempty(gdat_data.gdat_params.nfft)
-      nfft = gdat_data.gdat_params.nfft;
-    else
-      nfft=512;
-    end
-    tmhdm=mean(reshape(gdat_data.t(1:nfft*fix(length(gdat_data.t)/nfft)),nfft,fix(length(gdat_data.t)/nfft)));
-    for i=1:size(gdat_data.data,2)
-      [B,F,T]=specgram(gdat_data.data(:,i),nfft,1/mean(diff(gdat_data.t)),hanning(nfft),nfft/2);
+      tmhdm=mean(reshape(gdat_data(ishot).t(1:nfft*fix(length(gdat_data(ishot).t)/nfft)),nfft,fix(length(gdat_data(ishot).t)/nfft)));
+      for i=1:size(gdat_data(ishot).data,2)
+        [B,F,T]=specgram(gdat_data(ishot).data(:,i),nfft,1/mean(diff(gdat_data(ishot).t)),hanning(nfft),nfft/2);
+        figure;
+        imagesc(T+tmhdm(1),F/1e3,20*log10(abs(B)));axis xy;colormap jet;
+        ylabel('freq [kHz]')
+        xlabel(gdat_data(ishot).dimunits{1})
+        ylabel_eff = gdat_data(ishot).label;
+        if iscell(gdat_data(ishot).label) && length(gdat_data(ishot).label)>=i; ylabel_eff = gdat_data(ishot).label{i}; end
+        title([upper(gdat_data(ishot).gdat_params.machine) '#' num2str(gdat_data(ishot).shot) ' ' ylabel_eff])
+        mhd_sum_data = mhd_sum_data + gdat_data(ishot).data(:,i);
+      end
+      [B,F,T]=specgram(mhd_sum_data./size(gdat_data(ishot).data,2),nfft,1/mean(diff(gdat_data(ishot).t)),hanning(nfft),nfft/2);
       figure;
       imagesc(T+tmhdm(1),F/1e3,20*log10(abs(B)));axis xy;colormap jet;
       ylabel('freq [kHz]')
-      xlabel(gdat_data.dimunits{1})
-      ylabel_eff = gdat_data.label;
-      if iscell(gdat_data.label) && length(gdat_data.label)>=i; ylabel_eff = gdat_data.label{i}; end
-      title([upper(gdat_data.gdat_params.machine) '#' num2str(gdat_data.shot) ' ' ylabel_eff])
-      mhd_sum_data = mhd_sum_data + gdat_data.data(:,i);
+      xlabel(gdat_data(ishot).dimunits{1})
+      ylabel_eff = gdat_data(ishot).label;
+      if iscell(gdat_data(ishot).label); ylabel_eff = sprintf('%s ',gdat_data(ishot).label{:}); end
+      title([upper(gdat_data(ishot).gdat_params.machine) '#' num2str(gdat_data(ishot).shot) ' sum of ' ylabel_eff])
     end
-    [B,F,T]=specgram(mhd_sum_data./size(gdat_data.data,2),nfft,1/mean(diff(gdat_data.t)),hanning(nfft),nfft/2);
-    figure;
-    imagesc(T+tmhdm(1),F/1e3,20*log10(abs(B)));axis xy;colormap jet;
-    ylabel('freq [kHz]')
-    xlabel(gdat_data.dimunits{1})
-    ylabel_eff = gdat_data.label;
-    if iscell(gdat_data.label); ylabel_eff = sprintf('%s ',gdat_data.label{:}); end
-    title([upper(gdat_data.gdat_params.machine) '#' num2str(gdat_data.shot) ' sum of ' ylabel_eff])
+  else
+    disp(['cannot plot gdat_data(' num2str(ishot) '), has empty data or t field'])
   end
-else
-  disp('cannot plot gdat_data, has empty data or t field')
 end