function [OutImageFileCell,OutputMatrixCell]=imarith_fits1(Output,Input1,Operator,Input2,varargin)
%-----------------------------------------------------------------------------
% imarith_fits1 function                                              ImBasic
% Description: Perform arithmatic operation between images two sets of
%              images.
% Input  : - List of output images (see create_list.m for details).
%            If empty matrix, then set it to be equal to the Input1 list.
%          - First list of images (see create_list.m for details).
%            Alternatively this could be a scalar of vector (of scalars).
%            If the list contains only one image and the Input2 argument
%            contains multiple images then the single image will be
%            duplicated to a list with the same length as Input2.
%          - String containing any MATLAB operator.
%            Note that '.*' is multiplication, 
%            while '*' is matrix multiplication.
%          - Second list of images (see create_list.m for details).
%            Alternatively thiscould be a scalar of vector (of scalars).
%            If the list contains only one image and the Input1 argument
%            contains multiple images then the single image will be
%            duplicated to a list with the same length as Input1.
%          * Arbitrary number of pairs of input parameters: keyword,value,...
%            Available keywords are:
%            'Divide0'  - Set division by zero to specific value, default is 0.
%            'Scalar'   - Scalare operator {'y' | 'n'}, default is 'y'.
%                         If 'n' then the operator will be assume to be
%                         vectorize instead of scalar.
%            'OutPrefix'- Add prefix before output image names,
%                         default is empty string (i.e., '').
%            'OutDir'   - Directory in which to write the output images,
%                         default is empty string (i.e., '').
%            'Header'   - From which input to copy header {1 | 2} or to
%                         use minimal header {0}.
%                         Default is 1, however, if Input1 is scalar than
%                         default will be 2.
%            'AddHead'  - Cell array with 3 columns containing additional
%                         keywords to be add to the header.
%                         See cell_fitshead_addkey.m for header structure
%                         information. Default is empty matrix.
%            'DataType' - Output data type (see fitswrite.m for options), 
%                         default is float32.
%            'CCDSEC'   - Image sction for both images on which to do the
%                         arithmatic operation.
%                         This could be string containing image keyword 
%                         name (e.g., 'CCDSEC'), or a vector of 
%                         [Xmin, Xmax, Ymin, Ymax].
%                         If keyword than it will be obtained from
%                         the first input image.
%                         The output image size will be equal to that specified
%                         in the CCDSEC keyword, so in practice the output 
%                         image will be trimmed only if Trim='y'.
%                         If empty matrix (default; i.e., []) than do not use CCDSEC.
%            'Trim'     - Trim the resultant image {'y' | 'n'} according to CCDSEC).
%                         Default is 'y'. NOT AVAILABLE YET (behaves like default).
%             'Save'    - Save output FITS image {'y' | 'n'}.
%                         Default is 'y'.
% Output  : - Cell array containing output image names.
% Tested : Matlab 7.10
%     By : Eran O. Ofek                      June 2010
%    URL : http://wise-obs.tau.ac.il/~eran/matlab.html
% Example: Out=imarith_fits('try.fits','lred0015.fits','+',2);
%          Out=imarith_fits('@list1','@list2','.*','@list1','OutPrefix','out');
%          Out=imarith_fits('lred0[023-115].*','b0[023-115].*','.^','lred0[023-115].*','OutPrefix','out');
%          Out=imarith_fits('@list1','@list2','.*','lred0010.fits','OutPrefix','out');
%-----------------------------------------------------------------------------

DefV.Divide0     = 0;
DefV.Scalar      = 'y';
DefV.OutPrefix   = '';
DefV.OutDir      = '';
DefV.Header      = 1;
DefV.AddHead     = [];
DefV.DataType    = 'float32';
DefV.CCDSEC      = [];
DefV.Trim        = 'y';
DefV.Save        = 'y';
InPar = set_varargin_keyval(DefV,'y','use',varargin{:});



if (isnumeric(Input1)),
   InputCell1 = Input1;
else
   [~,InputCell1] = create_list(Input1,NaN);
end
if (isnumeric(Input2)),
   InputCell2 = Input2;
else
   [~,InputCell2] = create_list(Input2,NaN);
end
if (isempty(Output)),
   Output = Input1;
end
[~,OutputCell] = create_list(Output,NaN);

Ni1 = length(InputCell1);
Ni2 = length(InputCell2);
if (Ni1==0 || Ni2==0),
   % do nothing
else

   No  = length(OutputCell);
   if (Ni1==1 && Ni2>1),
      if (iscell(InputCell1)),
         [InputCell1{1:Ni2}] = deal(InputCell1{1});
      else
         [InputCell1(1:Ni2)] = deal(InputCell1(1));
      end
      Ni1 = Ni2;
   elseif (Ni2==1 && Ni1>1),
      if (iscell(InputCell2)),
         [InputCell2{1:Ni1}] = deal(InputCell2{1});
      else
         [InputCell2(1:Ni1)] = deal(InputCell2(1));
      end
      Ni2 = Ni1;
   elseif (Ni1==Ni2 && Ni1==No),
      % input length verified - do nothing
   else
      error('Illegal input/Ouput length');
   end
   
   switch lower(InPar.Scalar)
    case 'y'
       Operator = vectorize(Operator);
    case 'n'
       % do nothing
    otherwise
       error('Unknown Scalar option');
   end
   
   %--- Go over all images ---
   for I=1:1:Ni1,
       
      %--- read Image Input1 ---
      if (iscell(InputCell1)),
         % read image
         InImage1 = fitsread(InputCell1{I});
         StrIm1   = InputCell1{I};
         ImIsScalar1 = 0;
      else
         % assume image is a scalar
         InImage1 = InputCell1(I);
         StrIm1   = sprintf('%f',InputCell1(I));;
         ImIsScalar1 = 1;
         if (InPar.Header==1),
            % in case Header is request from the first image
            % than set it to image 2.
   	 InPar.Header = 2;
         end
      end
   
      %--- read Image Input2 ---
      if (iscell(InputCell2)),
         % read image
         InImage2 = fitsread(InputCell2{I});
         StrIm2   = InputCell2{I};
         ImIsScalar2 = 0;
      else
         % assume image is a scalar
         InImage2 = InputCell2(I);
         StrIm2   = sprintf('%f',InputCell2(I));;
         ImIsScalar2 = 1;
      end
   
      %--- CCDSEC ---
      if (isempty(InPar.CCDSEC)),
         % use entire image
         % do nothing
      elseif (ischar(InPar.CCDSEC)),
         [InPar.CCDSEC]=get_ccdsec_fits({InputCell1{I}},InPar.CCDSEC);
   
         [InImage1] = cut_image(InImage1,InPar.CCDSEC,'boundry');
         [InImage2] = cut_image(InImage2,InPar.CCDSEC,'boundry');
      
      elseif (length(InPar.CCDSEC)==4),
         [InImage1] = cut_image(InImage1,InPar.CCDSEC,'boundry');
         [InImage2] = cut_image(InImage2,InPar.CCDSEC,'boundry');
   
      else
         error('Illegal CCDSEC input');
      end
   
      OpStr = sprintf('OutImage = InImage1 %s InImage2;',Operator);
      eval(OpStr);
   
      switch Operator
       case '/'
          % deal with devision by zero
          Iinf = find(OutImage==Inf);
          OutImage(Iinf) = InPar.Divide0;
       otherwise
          % do nothing
      end
   
      OutImageFileName = sprintf('%s%s%s',InPar.OutDir,InPar.OutPrefix,OutputCell{I});
      OutImageFileCell{I} = OutImageFileName;
   
      %--- Prepare Header for output image ---
      switch InPar.Header
       case 0
          % use minimal header
          HeaderInfo = [];
       case 1
          Info = fitsinfo(InputCell1{I});
          HeaderInfo = Info.PrimaryData.Keywords;
       case 2
          Info = fitsinfo(InputCell2{I});
          HeaderInfo = Info.PrimaryData.Keywords;
       otherwise
          error('Unknown Header option');
      end
   
      %--- Add to header comments regarding file creation ---
      [HeaderInfo] = cell_fitshead_addkey(HeaderInfo,...
                                          Inf,'COMMENT','','Created by imarith_fits.m written by Eran Ofek',...
                                          Inf,'HISTORY','',OpStr,...
                                          Inf,'HISTORY','',sprintf('InImage1=%s',StrIm1),...
                                          Inf,'HISTORY','',sprintf('InImage2=%s',StrIm2));
      if (~isempty(InPar.AddHead)),
         %--- Add additional header keywords ---
         HeaderInfo = [HeaderInfo; InPar.AddHead];
      end
   
      %--- Write fits file ---
      switch lower(InPar.Save)
       case 'y'
          fitswrite(OutImage,OutImageFileName,HeaderInfo,InPar.DataType);
       case 'n'
          % do not save FITS image
       otherwise
          error('Unknown Save option');
      end
   
      if (nargout>1),
         OutputMatrixCell{I} = OutImage;
      end
   
   end
   
   
end
