function [SimMos]=sim_coadd_sub(Images,varargin)
%--------------------------------------------------------------------------
% sim_coadd_sub function                                           ImBasic
% Description: A version of sim_coadd.m that coadd images, but work on
%              section of the images at a time. This is designed to deal
%              with very large image stacks that can't fit into RAM
%              momeory.
% Input  : - Cell array of registered image names, image names with
%            wild cards and any input possible for create_list.m.
%          * Arbitrary number of pairs of arguments: ...,keyword,value,...
%            where keyword are one of the followings:
%            'SubSecSize' - sub section size [Y X] on which to work.
%                   Default is [256 256].
%            --- Additional parameters
%            Any additional key,val, that are recognized by one of the
%            following programs:
%               image2sim.m, image2sim.m, sim_coadd.m
% Output : - Coadded image.
% License: GNU general public license version 3
% Tested : Matlab R2014a
%     By : Eran O. Ofek                    Apr 2015
%    URL : http://weizmann.ac.il/home/eofek/matlab/
% Example: [Coadd]=sim_coadd_sub('*.fits','BackMethod','mode_fit','ScaleMethod','none','CombMethod','median');
% Reliable: 
%--------------------------------------------------------------------------


ImageField  = 'Im';
%HeaderField = 'Header';
%FileField   = 'ImageFileName';
%MaskField   = 'Mask';
%BackImField = 'BackIm';
%ErrImField  = 'ErrIm';
%WeightImField = 'WeightIm';

DefV.SubSecSize         = [256 256];
DefV.BackMethod         = 'none';   % {'none'|'mode_fit'|...}
DefV.ScaleMethod        = 'none';   % {'none'|'mode_fit'|...}
DefV.WeightMethod       = 'const';
DefV.WeightConst        = 1;
DefV.Filter             = {};
DefV.FiltPrep           = 'shift';
DefV.FiltNorm           = 1;
DefV.CombMethod         = 'sum';
DefV.Prctile            = 0.0;
DefV.Quantile           = 0.5;
DefV.CombMethodWeight   = 'fracnotnan';
DefV.CombMethodErr      = 'rss';   
DefV.CombMethodBack     = 'sum';
DefV.Verbose            = true;
InPar = set_varargin_keyval(DefV,'n','use',varargin{:});

Nfilt = numel(InPar.Filter);

[~,List] = create_list(Images,NaN);
Nim = numel(List);

InPar.WeightConst = InPar.WeightConst.*ones(Nim,1);

% read each image, subtract back, scale, weight, filter
ListTmp = cell(1,Nim);
for Iim=1:1:Nim,
    if (InPar.Verbose),
        fprintf('Image %d out of %d - background/flux/filter\n',Iim,Nim);
    end
    % read entire image
    Sim = image2sim(List{Iim},varargin{:});

    % remove background
    switch lower(InPar.BackMethod)
        case 'none'
            % do nothing
        otherwise
            Sim   = sim_background(Sim,varargin{:},'BackMethod',InPar.BackMethod);
    end

    % Scale images
    switch lower(InPar.ScaleMethod)
        case 'none'
            % do nothing
        otherwise
            [Sim] = sim_scaleflux(Sim,varargin{:},'ScaleMethod',InPar.ScaleMethod);
    end

    % Weight images
    switch lower(InPar.WeightMethod)
        case 'none'
            % do nothing
        otherwise
            [Sim] = sim_scaleflux(Sim,varargin{:},'ScaleMethod',InPar.WeightMethod,'ScaleConst',InPar.WeightConst(Iim));
    end
    
    % filter images
    if (~isempty(InPar.Filter)),
        Sim=sim_filter(Sim,varargin{:},'Back',false, 'Filter',InPar.Filter{min(Iim,Nfilt)}, 'FiltPrep',InPar.FiltPrep, 'FiltNorm',InPar.FiltNorm);
    end

    % save new images to disk
    ListTmp{Iim} = sim2fits(Sim,'OutName',List{Iim},'Prefix','tmpcoadd_');
end
[ImageSizeY,ImageSizeX] = size(Sim.(ImageField));
%% This section is to override a bug in fitsread.m
% it will remove the last column and row from the image
ImageSizeY = ImageSizeY - 1;
ImageSizeX = ImageSizeX - 1;
%%
Nsub = ceil([ImageSizeY,ImageSizeX]./InPar.SubSecSize);
NsubX = Nsub(2);
NsubY = Nsub(1);
% coadd each section seperatly
Ind = 0;
Npos = prod(Nsub);
CCDpos = zeros(Npos,4);
Coadd = simdef(Npos,1);
for Iy=1:1:NsubY,
    YLim = [(Iy-1), Iy].*InPar.SubSecSize(1) + [1 0];
    YLim(2) = min(YLim(2),ImageSizeY);
    for Ix=1:1:NsubX,
        XLim = [(Ix-1), Ix].*InPar.SubSecSize(1) + [1 0];
        XLim(2) = min(XLim(2),ImageSizeX);
        Ind = Ind + 1;
        if (InPar.Verbose),
             fprintf('Coadding subsection %d out of %d\n',Ind,Npos);
        end
        [Coadd(Ind)]=sim_coadd(ListTmp,varargin{:},'ImSec',[XLim YLim],...
                        'BackMethod','none',...
                        'ScaleMethod','none',...
                        'WeightMethod','none',...
                        'Filter',{});
                    
        CCDpos(Ind,:) = [XLim, YLim];
    end
end
[SimMos] = sim_mosaic(Coadd,CCDpos,varargin);
delete_cell(ListTmp);
