function [gdat_data_out] = gdat2time_out(gdat_data_in,option,others_as_data,varargin);
%
% function [gdat_data_out] = gdat2time_out(gdat_data_in,time_out,option,gdat_timedim,varargin);
%
% Aim: Reduce data output to desired time interval or compute values at desired time_instances
%
% Inputs:
%    gdat_data_in: gdat_data structure
%    option: 1: standard .data, .t, .dim{gdat_timedim} to be changed
%            21, 22: deal with grids_1d as well: 21 with rhopol 1D, 22 with rhopol 2D
%
%    others_as_data: {'field1','field2',...} other fields gdat_data_in.(others_as_data{i}) to be treated as .data
%
% Inferred/necessary inputs:
%                 gdat_data_in.gdat_params.time_out: time interval [t1,t2] or list of times to get output
%                 gdat_data_in.mapping_for.(gdat_data_in.gdat_params.machine).gdat_timedim: dim of time
%
% Outputs:
% gdat_data_out: with reduced data
%

gdat_data_out = gdat_data_in;

if isempty(gdat_data_in) || isempty(gdat_data_in.t) || ~isnumeric(gdat_data_in.data) || ~isnumeric(gdat_data_in.t)
  return
end

if ~exist('option') || isempty(option)
  option = 1;
end

iothers_as_data = false;
others_as_data_eff = '';
if exist('others_as_data') && ~isempty(others_as_data) && iscell(others_as_data)
  others_as_data_eff = others_as_data;
  iothers_as_data = true;
end

try
  dim_time = gdat_data_in.mapping_for.(gdat_data_in.gdat_params.machine).gdat_timedim;
catch ME
  try
    getReport(ME)
    error(['expect gdat_timedim in mapping_for.' gdat_data_in.gdat_params.machine ' from data_request: ' gdat_data_in.gdat_request]);
  catch ME2
    getReport(ME2)
    error(['expect gdat_timedim in mapping_for.machine fro data_request']);
  end
end

try
  time_out = sort(gdat_data_in.gdat_params.time_out);
  time_out = reshape(time_out,1,numel(time_out)); % usual shape if input as [t1 t2] or [t1,t2];
catch ME
  getReport(ME)
  error(['expect gdat_data_in.gdat_params.time_out'])
end

gdat_data_out = gdat_data_in;

% do this special part first otherwise ref time may have changed
if any(option==[21,22])
  % add extra fields then correct
  ab_tmp_rho = gdat_data_out.grids_1d;
  ab_tmp_rho.data = gdat_data_out.grids_1d.rhotornorm; ab_tmp_rho.t = gdat_data_in.t;
  ab_tmp_rho.gdat_params = gdat_data_out.gdat_params; ab_tmp_rho.mapping_for = gdat_data_out.mapping_for;
  extra_fields = {'psi','rhotornorm','rhovolnorm','rhopolnorm_equil','rhotornorm_equil', ...
                    'rhovolnorm_equil','rhotor_edge','volume_edge','b0'};
  if option==22; extra_fields{end+1} = 'rhopolnorm'; end
  [ab_tmp_rho] = gdat2time_out(ab_tmp_rho,1,extra_fields);
  gdat_data_out.grids_1d = rmfield(ab_tmp_rho,{'data','t','gdat_params','mapping_for'});
end

if any(option == [1, 21, 22])
  if numel(time_out) == 2
    % time interval provided, take all existing time points within this interval
    ij = [gdat_data_out.t>=time_out(1) & gdat_data_out.t<=time_out(2)];
    gdat_data_out.t = gdat_data_out.t(ij);
    if isfield(gdat_data_out,'dim'); gdat_data_out.dim{dim_time} = gdat_data_out.t; end
    nbdims = numel(size(gdat_data_out.data));
    nbdims_eff = nbdims;
    if nbdims==2 && any(size(gdat_data_out.data)==1)
      nbdims_eff = 1;
    end
    switch nbdims_eff
     case 1
      gdat_data_out.data = gdat_data_out.data(ij);
      if iothers_as_data
        for j=1:length(others_as_data_eff)
          eval(['gdat_data_out.' others_as_data_eff{j} ' = gdat_data_out.' others_as_data_eff{j} '(ij);']);
        end
      end
     case 2
      if dim_time ==1
        gdat_data_out.data = gdat_data_out.data(ij,:);
        if iothers_as_data
          for j=1:length(others_as_data_eff)
            eval(['atmp = gdat_data_out.' others_as_data_eff{j} ';']);
            if ~iscell(atmp)
              eval(['gdat_data_out.' others_as_data_eff{j} ' = gdat_data_out.' others_as_data_eff{j} '(ij,:);']);
            else
              eval(['gdat_data_out.' others_as_data_eff{j} '{1} = gdat_data_out.' others_as_data_eff{j} '{1}(ij);']);
            end
          end
        end
      else
        gdat_data_out.data = gdat_data_out.data(:,ij);
        if iothers_as_data
          for j=1:length(others_as_data_eff)
            eval(['atmp = gdat_data_out.' others_as_data_eff{j} ';']);
            if ~iscell(atmp)
              eval(['gdat_data_out.' others_as_data_eff{j} ' = gdat_data_out.' others_as_data_eff{j} '(:,ij);']);
            else
              eval(['gdat_data_out.' others_as_data_eff{j} '{2} = gdat_data_out.' others_as_data_eff{j} '{2}(ij);']);
            end
          end
        end
      end
     case 3
      if dim_time == 1
        gdat_data_out.data = gdat_data_out.data(ij,:,:);
        if iothers_as_data
          for j=1:length(others_as_data_eff)
            eval(['atmp = gdat_data_out.' others_as_data_eff{j} ';']);
            if ~iscell(atmp)
              eval(['gdat_data_out.' others_as_data_eff{j} ' = gdat_data_out.' others_as_data_eff{j} '(ij,:,:);']);
            else
              eval(['gdat_data_out.' others_as_data_eff{j} '{1} = gdat_data_out.' others_as_data_eff{j} '{1}(ij);']);
            end
          end
        end
      elseif dim_time == 2
        gdat_data_out.data = gdat_data_out.data(:,ij,:);
        if iothers_as_data
          for j=1:length(others_as_data_eff)
            eval(['atmp = gdat_data_out.' others_as_data_eff{j} ';']);
            if ~iscell(atmp)
              eval(['gdat_data_out.' others_as_data_eff{j} ' = gdat_data_out.' others_as_data_eff{j} '(:,ij,:);']);
            else
              eval(['gdat_data_out.' others_as_data_eff{j} '{2} = gdat_data_out.' others_as_data_eff{j} '{2}(ij);']);
            end
          end
        end
      else
        gdat_data_out.data = gdat_data_out.data(:,:,ij);
        if iothers_as_data
          for j=1:length(others_as_data_eff)
            eval(['atmp = gdat_data_out.' others_as_data_eff{j} ';']);
            if ~iscell(atmp)
              eval(['gdat_data_out.' others_as_data_eff{j} ' = gdat_data_out.' others_as_data_eff{j} '(:,:,ij);']);
            else
              eval(['gdat_data_out.' others_as_data_eff{j} '{3} = gdat_data_out.' others_as_data_eff{j} '{3}(ij);']);
            end
          end
        end
      end
     otherwise
      warning(['nbdims_eff = ' num2str(nbdims_eff) ' not yet defined in gdat2time_out, ask O. Sauter if needed']);
    end
  else
    % Assume 1 or several (>2) time points are requested. Interpolate to get data at these time instances specifically
    % default time_out_interp='linear', will implement {'interpos','tension','option',..} later when needed
    % time interval provided, take all existing time points within this interval
    ij = iround_os(gdat_data_out.t,time_out);
    gdat_data_out.t = time_out'; % usual shape n X 1  for time array in gdat, at least when 1D array
    if isfield(gdat_data_out,'dim'); gdat_data_out.dim{dim_time} = gdat_data_out.t; end
    nbdims = numel(size(gdat_data_out.data));
    nbdims_eff = nbdims;
    if nbdims==2 && any(size(gdat_data_out.data)==1)
      nbdims_eff = 1;
    end
    switch nbdims_eff
     case 1
      ij_ok = find(isfinite(gdat_data_in.data));
      if numel(ij_ok) > 1
        gdat_data_out.data = interpos(21,gdat_data_in.t(ij_ok),gdat_data_in.data(ij_ok),gdat_data_out.t);
      elseif numel(ij_ok) == 1
        gdat_data_out.data = gdat_data_in.data(ij_ok) .* ones(size(gdat_data_out.t));
      else
        gdat_data_out.data = NaN .* ones(size(gdat_data_out.t));
      end
      if iothers_as_data
        for j=1:length(others_as_data_eff)
          ij_ok = eval(['find(isfinite(gdat_data_in.' others_as_data_eff{j} '));']);
          if numel(ij_ok) > 1
            eval(['gdat_data_out.' others_as_data_eff{j} ' = interpos(21,gdat_data_in.t(ij_ok),gdat_data_in.' ...
                  others_as_data_eff{j} '(ij_ok),gdat_data_out.t);']);
          elseif numel(ij_ok) == 1
            eval(['gdat_data_out.' others_as_data_eff{j} ' = gdat_data_in.' others_as_data_eff{j} '(ij_ok) .* ones(size(gdat_data_out.t));']);
          else
            eval(['gdat_data_out.' others_as_data_eff{j} ' = NaN .* ones(size(gdat_data_out.t));']);
          end
        end
      end
     case 2
      len_t = numel(ij);
      if dim_time ==1
        gdat_data_out.data = NaN*ones(len_t,size(gdat_data_out.data,2));
        for ix=1:size(gdat_data_out.data,2)
          ij_ok = find(isfinite(gdat_data_in.data(:,ix)));
          if numel(ij_ok) > 1
            gdat_data_out.data([1:len_t],ix) = interpos(21,gdat_data_in.t(ij_ok),gdat_data_in.data(ij_ok,ix),gdat_data_out.t);
          elseif numel(ij_ok) == 1
            gdat_data_out.data([1:len_t],ix) = gdat_data_in.data(ij_ok,ix) .* ones(size(gdat_data_out.t));
          else
            gdat_data_out.data([1:len_t],ix) = NaN .* ones(size(gdat_data_out.t));
          end
        end
        if iothers_as_data
          for j=1:length(others_as_data_eff)
            eval(['atmp = gdat_data_out.' others_as_data_eff{j} ';']);
            if ~iscell(atmp)
              ix_len = eval(['size(gdat_data_out.' others_as_data_eff{j} ',2);']);
              eval(['gdat_data_out.' others_as_data_eff{j} ' =  NaN*ones(len_t,ix_len);']);
              for ix=1:ix_len
                ij_ok = eval(['find(isfinite(gdat_data_in.' others_as_data_eff{j} '(:,ix)));']);
                if numel(ij_ok) > 1
                  eval(['gdat_data_out.' others_as_data_eff{j} '([1:len_t],ix) = interpos(21,gdat_data_in.t(ij_ok),gdat_data_in.' ...
                        others_as_data_eff{j} '(ij_ok,ix),gdat_data_out.t);']);
                elseif numel(ij_ok) == 1
                  eval(['gdat_data_out.' others_as_data_eff{j} '([1:len_t],ix) = gdat_data_in.' ...
                        others_as_data_eff{j} '(ij_ok,ix) .* ones(size(gdat_data_out.t));']);
                else
                  eval(['gdat_data_out.' others_as_data_eff{j} '([1:len_t],ix) = NaN .* ones(size(gdat_data_out.t));']);
                end
              end
            else
              eval(['gdat_data_out.' others_as_data_eff{j} '{1} = NaN*ones(len_t,1);']);
              ij_ok = eval(['find(isfinite(gdat_data_in.' others_as_data_eff{j} '{1}));']);
              if numel(ij_ok) > 1
                eval(['gdat_data_out.' others_as_data_eff{j} '{1} = interpos(21,gdat_data_in.t(ij_ok),gdat_data_in.' ...
                      others_as_data_eff{j} '{1}(ij_ok),gdat_data_out.t);']);
              elseif numel(ij_ok) == 1
                eval(['gdat_data_out.' others_as_data_eff{j} '{1} = gdat_data_in.' ...
                      others_as_data_eff{j} '{1}(ij_ok) .* ones(size(gdat_data_out.t));']);
              else
                eval(['gdat_data_out.' others_as_data_eff{j} '{1} = NaN .* ones(size(gdat_data_out.t));']);
              end
            end
          end
        end
      else
        gdat_data_out.data = NaN*ones(size(gdat_data_out.data,1),len_t);
        for ix=1:size(gdat_data_out.data,1)
          ij_ok = find(isfinite(gdat_data_in.data(ix,:)));
          if numel(ij_ok) > 1
            gdat_data_out.data(ix,[1:len_t]) = interpos(21,gdat_data_in.t(ij_ok),gdat_data_in.data(ix,ij_ok),gdat_data_out.t);
          elseif numel(ij_ok) == 1
            gdat_data_out.data(ix,[1:len_t]) = gdat_data_in.data(ix,ij_ok) .* ones(size(gdat_data_out.t));
          else
            gdat_data_out.data(ix,[1:len_t]) = NaN .* ones(size(gdat_data_out.t));
          end
        end
        if iothers_as_data
          for j=1:length(others_as_data_eff)
            eval(['atmp = gdat_data_out.' others_as_data_eff{j} ';']);
            if ~iscell(atmp)
              ix_len = eval(['size(gdat_data_out.' others_as_data_eff{j} ',1);']);
              eval(['gdat_data_out.' others_as_data_eff{j} ' =  NaN*ones(ix_len,len_t);']);
              for ix=1:ix_len
                ij_ok = eval(['find(isfinite(gdat_data_in.' others_as_data_eff{j} '(ix,:)));']);
                if numel(ij_ok) > 1
                  eval(['gdat_data_out.' others_as_data_eff{j} '(ix,[1:len_t]) = interpos(21,gdat_data_in.t(ij_ok),gdat_data_in.' ...
                        others_as_data_eff{j} '(ix,ij_ok),gdat_data_out.t);']);
                elseif numel(ij_ok) == 1
                  eval(['gdat_data_out.' others_as_data_eff{j} '(ix,[1:len_t]) = gdat_data_in.' ...
                        others_as_data_eff{j} '(ix,ij_ok) .* ones(size(gdat_data_out.t));']);
                else
                  eval(['gdat_data_out.' others_as_data_eff{j} '(ix,[1:len_t]) = NaN .* ones(size(gdat_data_out.t));']);
                end
              end
            else
              eval(['gdat_data_out.' others_as_data_eff{j} '{2} = NaN*ones(1,len_t);']);
              ij_ok = eval(['find(isfinite(gdat_data_in.' others_as_data_eff{j} '{2}));']);
              if numel(ij_ok) > 1
                eval(['gdat_data_out.' others_as_data_eff{j} '{2} = interpos(21,gdat_data_in.t(ij_ok),gdat_data_in.' ...
                      others_as_data_eff{j} '{2}(ij_ok),gdat_data_out.t);']);
              elseif numel(ij_ok) == 1
                eval(['gdat_data_out.' others_as_data_eff{j} '{2} = gdat_data_in.' ...
                      others_as_data_eff{j} '{2}(ij_ok) .* ones(size(gdat_data_out.t));']);
              else
                eval(['gdat_data_out.' others_as_data_eff{j} '{2} = NaN .* ones(size(gdat_data_out.t));']);
              end
            end
          end
        end
      end
     case 3
      % do not assume can have a cell case until need to deal with it, since not sure what dim is left
      len_t = numel(ij);
      if dim_time == 1
        gdat_data_out.data = NaN*ones(len_t,size(gdat_data_out.data,2),size(gdat_data_out.data,3));
        for ix=1:size(gdat_data_out.data,2)
          for jy=1:size(gdat_data_out.data,3)
            ij_ok = find(isfinite(gdat_data_in.data(:,ix,jy)));
            if numel(ij_ok) > 1
              gdat_data_out.data([1:len_t],ix,jy) = interpos(21,gdat_data_in.t(ij_ok),gdat_data_in.data(ij_ok,ix,jy),gdat_data_out.t);
            elseif numel(ij_ok) == 1
              gdat_data_out.data([1:len_t],ix,jy) = gdat_data_in.data(ij_ok,ix,jy) .* ones(size(gdat_data_out.t));
            else
              gdat_data_out.data([1:len_t],ix,jy) = NaN .* ones(size(gdat_data_out.t));
            end
          end
        end
        if iothers_as_data
          for j=1:length(others_as_data_eff)
            ix_len = eval(['size(gdat_data_out.' others_as_data_eff{j} ',2);']);
            jy_len = eval(['size(gdat_data_out.' others_as_data_eff{j} ',3);']);
            eval(['gdat_data_out.' others_as_data_eff{j} ' =  NaN*ones(len_t,ix_len,jy_len);']);
            for ix=1:ix_len
              for jy=1:jy_len
                ij_ok = eval(['find(isfinite(gdat_data_in.' others_as_data_eff{j} '(:,ix,jy)));']);
                if numel(ij_ok) > 1
                  eval(['gdat_data_out.' others_as_data_eff{j} '([1:len_t],ix,jy) = interpos(21,gdat_data_in.t(ij_ok),gdat_data_in.' ...
                        others_as_data_eff{j} '(ij_ok,ix,jy),gdat_data_out.t);']);
                elseif numel(ij_ok) == 1
                  eval(['gdat_data_out.' others_as_data_eff{j} '([1:len_t],ix,jy) = gdat_data_in.' ...
                        others_as_data_eff{j} '(ij_ok,ix,jy) .* ones(size(gdat_data_out.t));']);
                else
                  eval(['gdat_data_out.' others_as_data_eff{j} '([1:len_t],ix,jy) = NaN .* ones(size(gdat_data_out.t));']);
                end
              end
            end
          end
          for j=1:length(others_as_data_eff)
          end
        end
      elseif dim_time == 2
        gdat_data_out.data = NaN*ones(size(gdat_data_out.data,1),len_t,size(gdat_data_out.data,3));
        for ix=1:size(gdat_data_out.data,1)
          for jy=1:size(gdat_data_out.data,3)
            ij_ok = find(isfinite(gdat_data_in.data(ix,:,jy)));
            if numel(ij_ok) > 1
              gdat_data_out.data(ix,[1:len_t],jy) = interpos(21,gdat_data_in.t(ij_ok),gdat_data_in.data(ix,ij_ok,jy),gdat_data_out.t);
            elseif numel(ij_ok) == 1
              gdat_data_out.data(ix,[1:len_t],jy) = gdat_data_in.data(ix,ij_ok,jy) .* ones(size(gdat_data_out.t));
            else
              gdat_data_out.data(ix,[1:len_t],jy) = NaN .* ones(size(gdat_data_out.t));
            end
          end
        end
        if iothers_as_data
          for j=1:length(others_as_data_eff)
            eval(['gdat_data_out.' others_as_data_eff{j} ' =  NaN*ones(size(gdat_data_out.' others_as_data_eff{j} ...
                  ',1),len_t,size(gdat_data_out.' others_as_data_eff{j} ',3));']);
          end
          for j=1:length(others_as_data_eff)
            ix_len = eval(['size(gdat_data_out.' others_as_data_eff{j} ',1);']);
            jy_len = eval(['size(gdat_data_out.' others_as_data_eff{j} ',3);']);
            for ix=1:ix_len
              for jy=1:jy_len
                ij_ok = eval(['find(isfinite(gdat_data_in.' others_as_data_eff{j} '(ix,:,jy)));']);
                if numel(ij_ok) > 1
                  eval(['gdat_data_out.' others_as_data_eff{j} '(ix,[1:len_t],jy) = interpos(21,gdat_data_in.t(ij_ok),gdat_data_in.' ...
                        others_as_data_eff{j} '(ix,ij_ok,jy),gdat_data_out.t);']);
                elseif numel(ij_ok) == 1
                  eval(['gdat_data_out.' others_as_data_eff{j} '(ix,[1:len_t],jy) = gdat_data_in.' ...
                        others_as_data_eff{j} '(ix,ij_ok,jy) .* ones(size(gdat_data_out.t));']);
                else
                  eval(['gdat_data_out.' others_as_data_eff{j} '(ix,[1:len_t],jy) = NaN .* ones(size(gdat_data_out.t));']);
                end
              end
            end
          end
        end
      else
        gdat_data_out.data = NaN*ones(size(gdat_data_out.data,1),size(gdat_data_out.data,2),len_t);
        for ix=1:size(gdat_data_out.data,1)
          for jy=1:size(gdat_data_out.data,2)
            ij_ok = find(isfinite(gdat_data_in.data(ix,jy,:)));
            if numel(ij_ok) > 1
              gdat_data_out.data(ix,jy,[1:len_t]) = interpos(21,gdat_data_in.t(ij_ok),gdat_data_in.data(ix,jy,ij_ok),gdat_data_out.t);
            elseif numel(ij_ok) == 1
              gdat_data_out.data(ix,jy,[1:len_t]) = gdat_data_in.data(ix,jy,ij_ok) .* ones(size(gdat_data_out.t));
            else
              gdat_data_out.data(ix,jy,[1:len_t]) = NaN .* ones(size(gdat_data_out.t));
            end
          end
        end
        if iothers_as_data
          for j=1:length(others_as_data_eff)
            eval(['gdat_data_out.' others_as_data_eff{j} ' =  NaN*ones(size(gdat_data_out.' others_as_data_eff{j} ...
                  ',1),size(gdat_data_out.' others_as_data_eff{j} ',2),len_t);']);
          end
          for j=1:length(others_as_data_eff)
            ix_len = eval(['size(gdat_data_out.' others_as_data_eff{j} ',1);']);
            jy_len = eval(['size(gdat_data_out.' others_as_data_eff{j} ',2);']);
            for ix=1:ix_len
              for jy=1:jy_len
                ij_ok = eval(['find(isfinite(gdat_data_in.' others_as_data_eff{j} '(ix,jy,:)));']);
                if numel(ij_ok) > 1
                  eval(['gdat_data_out.' others_as_data_eff{j} '(ix,jy,[1:len_t]) = interpos(21,gdat_data_in.t(ij_ok),gdat_data_in.' ...
                        others_as_data_eff{j} '(ix,jy,ij_ok),gdat_data_out.t);']);
                elseif numel(ij_ok) == 1
                  eval(['gdat_data_out.' others_as_data_eff{j} '(ix,jy,[1:len_t]) = gdat_data_in.' ...
                        others_as_data_eff{j} '(ix,jy,ij_ok) .* ones(size(gdat_data_out.t));']);
                else
                  eval(['gdat_data_out.' others_as_data_eff{j} '(ix,jy,[1:len_t]) = NaN .* ones(size(gdat_data_out.t));']);
                end
              end
            end
          end
        end
      end
     otherwise
      warning(['nbdims_eff = ' num2str(nbdims_eff) ' not yet defined in gdat2time_out, ask O. Sauter if needed']);
    end
  end
else
  error(['option = ' num2str(option) ' not yet implemented'])
end