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);
catch ME
  getReport(ME)
  error(['expect gdat_data_in.gdat_params.time_out'])
end

gdat_data_out = gdat_data_in;

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;
    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 = interp1(gdat_data_in.t,gdat_data_in.data,gdat_data_out.t);
      if iothers_as_data
        for j=1:length(others_as_data_eff)
          eval(['gdat_data_out.' others_as_data_eff{j} ' = interp1(gdat_data_in.t,gdat_data_in.' others_as_data_eff{j} ',gdat_data_out.t);']);
        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)
          gdat_data_out.data([1:len_t],ix) = interp1(gdat_data_in.t,gdat_data_in.data(:,ix),gdat_data_out.t);
        end
        if iothers_as_data
          for j=1:length(others_as_data_eff)
            eval(['gdat_data_out.' others_as_data_eff{j} ' =  NaN*ones(len_t,size(gdat_data_out.' others_as_data_eff{j} ',2));']);
          end
          for j=1:length(others_as_data_eff)
            ix_len = eval(['size(gdat_data_out.' others_as_data_eff{j} ',2);']);
            for ix=1:ix_len
              eval(['gdat_data_out.' others_as_data_eff{j} '([1:len_t],ix) = interp1(gdat_data_in.t,gdat_data_in.' others_as_data_eff{j} '(:,ix),gdat_data_out.t);']);
            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)
          gdat_data_out.data(ix,[1:len_t]) = interp1(gdat_data_in.t,gdat_data_in.data(ix,:),gdat_data_out.t);
        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);']);
          end
          for j=1:length(others_as_data_eff)
            ix_len = eval(['size(gdat_data_out.' others_as_data_eff{j} ',1);']);
            for ix=1:ix_len
              eval(['gdat_data_out.' others_as_data_eff{j} '(ix,[1:len_t]) = interp1(gdat_data_in.t,gdat_data_in.' others_as_data_eff{j} '(ix,:),gdat_data_out.t);']);
            end
          end
        end
      end
     case 3
      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)
            gdat_data_out.data([1:len_t],ix,jy) = interp1(gdat_data_in.t,gdat_data_in.data(:,ix,jy),gdat_data_out.t);
          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(len_t,size(gdat_data_out.' others_as_data_eff{j} ...
                  ',2),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} ',2);']);
            jy_len = eval(['size(gdat_data_out.' others_as_data_eff{j} ',3);']);
            for ix=1:ix_len
              for jy=1:jy_len
                eval(['gdat_data_out.' others_as_data_eff{j} '([1:len_t],ix,jy) = interp1(gdat_data_in.t,gdat_data_in.' others_as_data_eff{j} '(:,ix,jy),gdat_data_out.t);']);
              end
            end
          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)
            gdat_data_out.data(ix,[1:len_t],jy) = interp1(gdat_data_in.t,gdat_data_in.data(ix,:,jy),gdat_data_out.t);
          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
                eval(['gdat_data_out.' others_as_data_eff{j} '(ix,[1:len_t],jy) = interp1(gdat_data_in.t,gdat_data_in.' others_as_data_eff{j} '(ix,:,jy),gdat_data_out.t);']);
              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)
            gdat_data_out.data(ix,jy,[1:len_t]) = interp1(gdat_data_in.t,gdat_data_in.data(ix,jy,:),gdat_data_out.t);
          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
                eval(['gdat_data_out.' others_as_data_eff{j} '(ix,jy,[1:len_t]) = interp1(gdat_data_in.t,gdat_data_in.' others_as_data_eff{j} '(ix,jy,:),gdat_data_out.t);']);
              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

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_out.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
