function Spec=sedm_extract_spec_cube(varargin)
%--------------------------------------------------------------------------
% sedm_extract_spec_cube function                                     SEDM
% Description: Extract a source and background spectrum from a SEDM cube.
% Input  : * Arbitrary number of pairs of ...,key,val,.. arguments.
%            The following keywords are available:
%            'Cube'  - Cube structure as returned by
%                      sedm_segments2cube.m.
%                      This can be either a Cube structure or a mat file
%                      name containing a cube structure.
%                      This parameter must be provided.
%            'Wave'  - Vector of wavelengths corresponding to the first
%                      dimension of the cube.
%                      Default is logspace(log10(360),log10(950),200).'.
%                      However, if the Cube structure contains a 'Wave'
%                      field, then the wavelength vector will be taken from
%                      this field.
%            'DispWaveRange' - Wavelength range from which to build
%                      and present the 2D image.
%                      Default is [400 900] (i.e., in nm).
%            'Display' - Type of display in which to show the image.
%                      {'ds9'}. Default is 'ds9'.
%            'UseExistDisp' - Use existing open display (true) or open
%                      a new display even if one is already opened (false).
%                      Default is true.
%            'Pos'   - [X, Y] positions in the cube of the source to
%                      extract. Of empty then use interactive (dispaly)
%                      mode. Default is empty. Note that if several lines
%                      are provided then multiple positions will be
%                      extracted.
%                      If 'max' then will extract the source around the
%                      highest pixel in the image.
%                      If 'sex' then use SExtractor to find sources in
%                      image.
%            'Aper'  - Aperture radii of extraction aperture and inner
%                      and outer radii of the background annulus.
%                      Default is [2 8 12] in the cube pixel units.
%            'ImageName' - Name of the output 2D image.
%                      Default is 'tmpImage.fits'.
%            'Frame'  - Frame number to use in ds9. Default is 1.
%            'Centeroid' - Method for centeroiding the source.
%                       {'none'|'moment'}. Default is 'moment'.
%            'Gain'   - Instrument gain for error calculations.
%                       Default is 1.
%            'PlotSpec' - Plot extracted spectrum {false|true}.
%                       Default is true.
%            'InvTran'- Structure containing the transmission information
%                       as generated by spec_fluxcalib.m.
%                       Alternatively this can be a mat file name
%                       containing the structure.
%                       If empty then don't correct for transmission.
%                       Default is 'sedm_InvTran.mat'.
%            'PlotCalib' - Plot the calibrated spectrum (true) rather
%                       than the spectrum in ADU units (false).
%                       Default is true.
%            'AM'     - AirMass of image. The airmass is used to correct
%                       for atmospheric extinction. Default is 1.
%            'ExpTime'- Exposure time of image [seconds]. Default is 1.
%                       This is used for flux calibrate the spectrum.
%            'InterpMethod' - Interpolation method. See interp1.m for
%                       options. Default is 'linear'.
%            'SexKeys'- SExtractor parameters that will be derived from
%                       image. Default is 
%                       {'XWIN_IMAGE','YWIN_IMAGE','FLUX_AUTO','FLUXERR_AUTO'}.
%            'SexPars'- Additional parameters to pass to SExtractor.
%                       Default is {'DETECT_THRESH','2'}. 
% Output : Structure array containing the extracted spectrum of the sources
%          and background. The following fields are available:
%          .Xpix         - Source X position [pixels].
%          .Ypix         - Source Y position [pixels].
%          .Wave         - Vector of wavelengths.
%          .BS_Source    - Background subtracted source intensity.
%          .BS_SourceErr - Background subtracted source intensity error.
%          .Source       - Source intensity before background subtraction.
%          .Back         - median background intensity per pixel.
%          .BackSpex     - median background from the spexcell extraction
%                          stage (ADU/pix).
%          .SourceErr    - Error on source intensity.
%          .BackErr      - Error on background intesnity.
%          .Aper         - The input Aper vector.
%          .SrcNpix      - Number of cube pixels in source extraction
%                          region.
%          .BckNpix      - Number of cube pixels in background extraction
%                          region.
% Tested : Matlab R2011b
%     By : Eran O. Ofek                    Sep 2013
%    URL : http://weizmann.ac.il/home/eofek/matlab/
% Example: Spec=sedm_extract_spec_cube('Cube',Cube);
% Reliable: 2
%--------------------------------------------------------------------------

WaveField = 'Wave';
NM_TO_ANG = 10;

DefV.Cube           = [];
DefV.Wave           = logspace(log10(360),log10(950),200).';
%DefV.SI             = [];
DefV.DispWaveRange  = [400 900];
DefV.Display        = 'ds9';   % {'ds9'}
DefV.UseExistDisp   = true;
DefV.Pos            = [];         
DefV.Aper           = [2 8 12];   % [pixels of cube]
DefV.ImageName      = 'tmpImage.fits';   % FITS image name to save
DefV.Frame          = 1;
DefV.Centeroid      = 'moment';   % {'none'|'moment'}
DefV.Gain           = 1;
DefV.PlotSpec       = true;
DefV.InvTran        = 'sedm_InvTran.mat';
DefV.PlotCalib      = true;
DefV.AM             = 1;
DefV.ExpTime        = 1;          % [s]
DefV.Ext            = 'KPNO_atmospheric_extinction.dat';
DefV.InterpMethod   = 'linear';
DefV.SexKeys        = {'XWIN_IMAGE','YWIN_IMAGE','FLUX_AUTO','FLUXERR_AUTO'};
DefV.SexPars        = {'DETECT_THRESH','2'};

InPar = set_varargin_keyval(DefV,'y','use',varargin{:});

if (isempty(InPar.Cube)),
    error('Cube must be provided');
end
if (ischar(InPar.Cube)),
    % read Cube from mat file
    InPar.Cube = load2(InPar.Cube);
end

if (ischar(InPar.InvTran)),
    % read from mat file
    Tran = load2(InPar.InvTran);
else
    Tran = InPar.InvTran;
end

if (isfield(InPar.Cube,WaveField)),
    InPar.Wave = InPar.Cube.(WaveField);
end

Iw = InPar.Wave>=InPar.DispWaveRange(1) & InPar.Wave<=InPar.DispWaveRange(2);



% collapse cube to 2D image
Image = squeeze(sum(InPar.Cube.Cube(Iw,:,:),1));
% store 2D image
fitswrite(Image,InPar.ImageName);

% build the image coordinate system
[MatX,MatY] = meshgrid((1:1:size(Image,2)),(1:1:size(Image,1)));
Cube2d      = reshape(InPar.Cube.Cube,length(InPar.Wave),numel(Image));
MatX2d      = reshape(MatX,1,numel(Image));
MatY2d      = reshape(MatY,1,numel(Image));

% The Pos parameter options
if (ischar(InPar.Pos)),
    switch lower(InPar.Pos)
        case 'max'
            [~,MaxInd] = maxnd(Image);
            InPar.Pos  = [MaxInd(2), MaxInd(1)];
        case 'sex'
            
            [Cat,~,ColC] = run_sextractor(InPar.ImageName,[],...
                                          InPar.SexKeys,[],[],...
                                          InPar.SexPars{:});
            InPar.Pos    = Cat{1}(:,[ColC.XWIN_IMAGE, ColC.YWIN_IMAGE]);
            if (isempty(InPar.Pos)),
                warning('No sources found using SExtractor');
            end
                                      
        otherwise
            error('Unknown Pos option');
    end
end

ExtractSpec = true;
SpecInd     = 1;
Npos        = size(InPar.Pos,1);

switch lower(InPar.Display)
    case 'ds9'
        
        % loop over all positions
        % or while the user feed more positions
        while ExtractSpec
            
            if (isempty(InPar.Pos)),

                if (InPar.UseExistDisp),
                    % check if ds9 exist and if not open new
                    ds9_start;
                else
                    % open new ds9 display
                    system('ps -C ds9');
                end

                % display image
                %fitswrite(Image,InPar.ImageName);

                ds9_disp(InPar.ImageName,InPar.Frame,'Zoom',10);
                
                fprintf('\n');
                fprintf('Select source to extract using ank key/mouse\n');
                fprintf('  Use q to select and quit\n');

                [UserX,UserY,~,Key]=ds9_getcoo(1,'image','any');
                % plot source
                switch lower(InPar.Centeroid)
                    case 'none'
                        CooX = UserX;
                        CooY = UserY;
                    case 'moment'
                        M1   = moment_2d(cut_image(Image,...
                                         [round(UserX) round(UserY) ceil(InPar.Aper(1)) ceil(InPar.Aper(1))],'center'),UserX,UserY);
                        CooX = M1(1);
                        CooY = M1(2);
                    otherwise
                        error('Unknown Centeroid option');
                end
                % plot extracted source and background in ds9
                ds9_plotregion(CooX,CooY,'Coo','image','Type','circle','Size',InPar.Aper(1),'Color','red');
                ds9_plotregion([CooX;CooX],[CooY;CooY],'Coo','image','Type','circle','Size',InPar.Aper(2:3).','Color','green');
                
                switch lower(Key{1})
                    case 'q'
                        % Exit
                        ExtractSpec = false;
                    otherwise
                        ExtractSpec = true;
                         
                end
                
            else
                
                CooX = InPar.Pos(SpecInd,1);
                CooY = InPar.Pos(SpecInd,2);
                
                if (SpecInd>Npos),
                    ExtractSpec = false;
                end
            end

            % run aperture "photometry"
	    MatDist = Util.Geom.plane_dist(MatX2d,MatY2d,CooX,CooY);

            Isource = find(MatDist(:)<=InPar.Aper(1));
            Iback   = find(MatDist(:)>InPar.Aper(2) & MatDist(:)<InPar.Aper(3));

            % store info into Spec output argument
            Spec(SpecInd).Xpix        = CooX;
            Spec(SpecInd).Ypix        = CooY;
            Spec(SpecInd).Aper        = InPar.Aper;
            Spec(SpecInd).SrcNpix     = length(Isource);
            Spec(SpecInd).BckNpix     = length(Iback);
            Spec(SpecInd).Wave        = InPar.Wave;
            Spec(SpecInd).Source      = sum(Cube2d(:,Isource),2);
            Spec(SpecInd).SourceErr   = sqrt(Spec(SpecInd).Source./InPar.Gain);
            Spec(SpecInd).Back        = nanmedian(Cube2d(:,Iback),2);
            Spec(SpecInd).BackSpex    = nanmedian(InPar.Cube.Back(Iback));  % background from the spexcell extraction stage
            
            Spec(SpecInd).BackErr     = sqrt(Spec(SpecInd).Back./InPar.Gain./Spec(SpecInd).BckNpix);
            Spec(SpecInd).BS_Source   = Spec(SpecInd).Source - Spec(SpecInd).Back.* Spec(SpecInd).SrcNpix;  % back subtracted source spec
            Spec(SpecInd).BS_SourceErr= sqrt(Spec(SpecInd).Source + (Spec(SpecInd).Back + Spec(SpecInd).BackSpex).* Spec(SpecInd).SrcNpix);
            if (isempty(Tran) || ~InPar.PlotCalib),
                % spectrum is not flux calibrated
                Spec(SpecInd).FC           = false;
                Spec(SpecInd).FC_Source    = Spec(SpecInd).BS_Source;
                Spec(SpecInd).FC_SourceErr = Spec(SpecInd).BS_SourceErr;
                Spec(SpecInd).FC_Back      = Spec(SpecInd).Back;
                Spec(SpecInd).FC_BackErr   = Spec(SpecInd).BackErr;    
            else
                % flux calibrate the spectrum
                
                InvTran = interp1(Tran.Wave,Tran.InvTran,Spec(SpecInd).Wave.*NM_TO_ANG,InPar.InterpMethod);

                AirlessSpec = atmospheric_ext([Spec(SpecInd).Wave.*NM_TO_ANG,Spec(SpecInd).BS_Source],...
                                              InPar.AM,...
                                              InPar.Ext,...
                                              InPar.InterpMethod);
                
                Spec(SpecInd).FC           = true;              
                Spec(SpecInd).FC_Source    = AirlessSpec(:,2).*InvTran./InPar.ExpTime;
                Spec(SpecInd).FC_SourceErr = Spec(SpecInd).FC_Source.* (Spec(SpecInd).BS_SourceErr./Spec(SpecInd).BS_Source);
                Spec(SpecInd).FC_Back      = Spec(SpecInd).FC_Source.* (Spec(SpecInd).Back./Spec(SpecInd).BS_Source);
                Spec(SpecInd).FC_BackErr   = Spec(SpecInd).FC_Back.* (Spec(SpecInd).BackErr./Spec(SpecInd).Back);
            end
            
            % Plot the spectrum
            if (InPar.PlotSpec),
                clf;
                
                stairs(Spec(SpecInd).Wave,Spec(SpecInd).FC_Source,'k-','LineWidth',2);
                hold on;
                stairs(Spec(SpecInd).Wave,Spec(SpecInd).FC_Back.*Spec(SpecInd).SrcNpix,'k-','Color',[0.8 0.8 0.8]);
                box on;
                legend('Source','Background');
                patch_band(Spec(SpecInd).Wave,Spec(SpecInd).FC_Source,Spec(SpecInd).FC_SourceErr,[0.8 0.8 0.8]);
                plot_invchildren
                H = xlabel('Wavelength [nm]');
                set(H,'FontSize',18);
                if (Spec(SpecInd).FC),
                    H = ylabel('Flux [erg cm^{-2} s^{-1} A^{-1}]');
                else
                    H = ylabel('Counts');
                end
                set(H,'FontSize',18);
            end
            
            SpecInd = SpecInd + 1;
            if (SpecInd>Npos && Npos>0),
                ExtractSpec = false;
            end
        end
        
    otherwise
        error('Unkown Display option');
end
