function varargout = interp2RAPTOR_tgrid(tgridin,tgrid,rhoin,varargin)
% function varargout = interp2RAPTOR_tgrid(tgridin,tgrid,varargin)
% tgridin: time grid of the input data
% tgrid: RAPTOR time grid
% Interpolation from the input to RAPTOR time grid

% define number of input arguments
narg = numel(varargin);
% define number of elements in the input time grid
nbtin = numel(tgridin);
% define number of elements in the RAPTOR time grid
nbt = numel(tgrid);
% define a cell array of output arguments 
varargout = cell(1,narg+1);

% Sort input time grid from the smallest to the biggest value
[tgridhat,tind] = sort(tgridin);

% repmat radial grid
varargout{1} = repmat(rhoin,1,nbt);

% for each input argument
for ij=1:narg
    d0 = varargin{ij};
    Dhat = zeros(size(d0));
    % sort data in the correct time order
    for ii=1:nbtin
        Dhat(:,tind(ii)) = d0(:,ii);
    end
    
    ng = size(Dhat,1);
    D = zeros(ng,nbt);

    if numel(tgridhat)==1
        % if input time grid contains only one element
        D = repmat(Dhat,1,nbt);
    else
        % define if the first point in input and RAPTOR time grids is different
        it1 = find(tgrid<=tgridhat(1));
        if isempty(it1)
            it1 = 0;
        else
            % parameters values before the first time step in the input
            % time grid are the same as in it first time step
            D(:,it1) = repmat(Dhat(:,1),1,numel(it1));
        end
        % define if the last point in input and RAPTOR time grids is different
        it2 = find(tgrid>=tgridhat(end));
        if isempty(it2)
            it2 = numel(tgrid)+1;
        else
            % parameters values after the last time step in the input time
            % grid are the same as in it last time step
            D(:,it2) = repmat(Dhat(:,end),1,numel(it2));
        end         
        % define indexes for an overlapping part of the input and RAPTOR time
        % grids
        it3 = it1(end)+1:it2(1)-1;  
        % interpolate data form input time grid to RAPTOR time grid
        D(:,it3) = interp1(tgridhat,Dhat',tgrid(it3))'; %,'v5cubic'
    end
    varargout{ij+1} = D;
end
return