function [ResPhot,OutRef,Extra]=phot_calib(Sim,varargin)
% Photometrically calibrate an AstCat object for an image.
% Package: @AstCat
% Description: For each catalog in an AstCat object find matched sources in
%              a reference catalog and use the sources to estimate the
%              image photometric zero point.
% Input  : - A Sim or AstCat object.
%          * Arbitrary number of pairs of arguments: ...,keyword,value,...
%            where keyword are one of the followings:
%            'OutRef' - The matched reference catalog generated by
%                       SIM/xcat. If empty then call xcat.
%                       Default is empty.
%            'Extra'  -
%            'xcatPar'- 
%            'RefCat' -
%            'CatColRA'-
%            'CatColDec'-
%            'CatUnits'-
%            'CatColFlags'
%            'CatBadFlags'
%            'RefColRA'
%            'RefColDec'
%            'RefUnits'
%            'RefColFlags'
%            'RefBadFlags'
%            'RefColStar'
%            'RefColStarFun'
%            'MaxRefMag'
%            'MinRefMag'
%            'RefMaxMagErr'
%            'CatMaxMagErr'
%            'LocalRA'
%            'LocalDec'
%            'LocalRadius'
%            'CatColMag'            = 'MAG_PSF';
%            'CatColMagErr'         = 'MAGERR_PSF';
%            'RefRefMag'            = 'modelMag_r';
%            'RefRefMagErr'         = 'modelMagErr_r';
%            'RefAbsErr'            = 0.015;
%            'CalibEqConst'         = true;
%            'CalibEqRef'           = {'modelMag_g-modelMag_r'};
%            'CalibEqCat'           = {};
%            'Verbose'              = true;
% Output : - 
% License: GNU general public license version 3
%     By : Eran O. Ofek                    Jul 2017
%    URL : http://weizmann.ac.il/home/eofek/matlab/
% Example: ResPhot=phot_calib(Sim);
% Reliable: 
%--------------------------------------------------------------------------

CatField    = AstCat.CatField;


DefV.OutRef               = [];
DefV.Extra                = [];
DefV.xcatPar              = {};

%DefV.KeyCCDSEC            = 'CCDSEC';   % if empty use actual image size
DefV.RefCat               = 'sdss';     % AstCat or catalaog name
%DefV.SearchRad            = 2;
%DefV.SearchRadUnits       = 'arcsec';
% generate catalog
%DefV.GenCatFun            = @mextractor;
%DefV.GenCatPar            = {};
% catalaog attributes
DefV.CatColRA             = {'ALPHAWIN_J2000','RA'};
DefV.CatColDec            = {'DELTAWIN_J2000','Dec'};
DefV.CatUnits             = 'deg';
DefV.CatColFlags          = {'FLAGS'};
DefV.CatBadFlags          = {'Bit_ImSaturated','Bit_BiasNoisy','Bit_BiasNoise0','Bit_BiasAnom','Bit_BadPixel','Bit_FlatNaN','Bit_FlatLowNim','Bit_FlatHighStd','Bit_FlatLow','Bit_CR','Bit_Edge','Bit_Spike','Bit_Ghost','Bit_NonLin','Bit_BackGrad'};
% reference attributes
DefV.RefColRA             = {'RA','ALPHAWIN_J2000'};
DefV.RefColDec            = {'Dec','DELTAWIN_J2000'};
DefV.RefUnits             = 'rad';

DefV.RefColFlags          = {'FLAGS'};
DefV.RefBadFlags          = {};
DefV.RefColStar           = {'type'};
DefV.RefColStarFun        = @(Val) Val==6;

DefV.MaxRefMag            = 19;
DefV.MinRefMag            = 15;
DefV.RefMaxMagErr         = 0.02;
DefV.CatMaxMagErr         = 0.02;

DefV.LocalRA              = [];
DefV.LocalDec             = [];
DefV.LocalRadius          = [];

DefV.CatColMag            = 'MAG_PSF';
DefV.CatColMagErr         = 'MAGERR_PSF';
DefV.RefRefMag            = 'modelMag_r';
DefV.RefRefMagErr         = 'modelMagErr_r';
DefV.RefAbsErr            = 0.015;

DefV.CalibEqConst         = true;
DefV.CalibEqRef           = {'modelMag_g-modelMag_r'};
DefV.CalibEqCat           = {};

DefV.Verbose              = true;

InPar = InArg.populate_keyval(DefV,varargin,mfilename);

% cross match the catalogs with external reference catalogs
if (isempty(InPar.OutRef) || isempty(InPar.Extra))

    [OutRef,Extra]=xcat(Sim,'CatColRA',InPar.CatColRA,...
                            'CatColDec',InPar.CatColDec,...
                            'RefColRA',InPar.RefColRA,...
                            'RefColDec',InPar.RefColDec,...
                            'RefCat',InPar.RefCat,...
                            InPar.xcatPar{:});
else
    OutRef = InPar.OutRef;
    Extra  = InPar.Extra;
end


% some preparation
switch lower(InPar.RefCat)
    case 'apass'
        InPar.RefColStar = {};
end





Nsim = numel(Sim);
for Isim=1:1:Nsim
   
    %---------------
    % clean matches
    %---------------
    Nsrc = size(Sim(Isim).(CatField),1);
    
    MagCat    = col_get(Sim(Isim),InPar.CatColMag);
    MagErrCat = col_get(Sim(Isim),InPar.CatColMagErr);
    MagRef    = col_get(OutRef(Isim),InPar.RefRefMag);
    MagErrRef = col_get(OutRef(Isim),InPar.RefRefMagErr);
    
    % clean by input catalog flags
    % TBD
    Flag.CatFlags = true(Nsrc,1);
    
    % clean by reference catalog flags
    % TBD
    Flag.RefFlags = true(Nsrc,1);
    
    % clean by reference catalog object type
    if (~isempty(InPar.RefColStar))
        [~,RefColStar]  = select_exist_colnames(OutRef,InPar.RefColStar');
        Flag.RefType    = InPar.RefColStarFun(OutRef(Isim).(CatField)(:,RefColStar));
    else
        Flag.RefType    = true(Nsrc,1);
    end
    
    % clean by position
    if (~isempty(InPar.LocalRA) && ~isempty(InPar.LocalDec) && ~isempty(InPar.LocalRadius))
        [~,RefColRA]     = select_exist_colnames(OutRef,InPar.RefColRA');
        [~,RefColDec]    = select_exist_colnames(OutRef,InPar.RefColDec');
        RefRA  = col_get(OutRef(Isim),RefColRA);
        RefDec = col_get(OutRef(Isim),RefColDec);
        
        D = celestial.coo.sphere_dist_fast(InPar.LocalRA,InPar.LocalDec,RefRA,RefDec);
        Flag.FlagLocal = D<=InPar.LocalRadius;
    else
        Flag.FlagLocal = true(Nsrc,1);
    end

    % magnitide is defined
    Flag.MagFlag = ~isnan(MagCat) & ~isnan(MagRef);
    
    % magnitude range
    Flag.MagRange = MagRef>InPar.MinRefMag & MagRef<InPar.MaxRefMag;
    
    % magnitude errors
    Flag.MagErr   = MagErrRef<InPar.RefMaxMagErr & MagErrCat<InPar.CatMaxMagErr;
    
    % combine all flags
    FlagAll = Flag.CatFlags & Flag.RefFlags & Flag.RefType & Flag.FlagLocal & Flag.MagFlag & Flag.MagRange & Flag.MagErr;
        
    % check that there are enough stars
    Nc = sum(FlagAll);
    if (InPar.Verbose)
        fprintf('Use %d calibrators\n',Nc);
    end
        
    
    YY    = MagCat - MagRef;
    VarYY = MagErrCat.^2 + MagErrRef.^2 + InPar.RefAbsErr.^2;   % variance
    
    NeqRef = numel(InPar.CalibEqRef);
    NeqCat = numel(InPar.CalibEqCat);
    
    % prepare the design matrix
    CurCol = 0;
    if (InPar.CalibEqConst)
        Npar = NeqRef+NeqCat+1;
        H = zeros(Nsrc,Npar);
        CurCol = CurCol+1;
        H(:,CurCol) = 1;
    else
        Npar = NeqRef+NeqCat;
        H = zeros(Nsrc,Npar);
    end
    
    for IeqRef=1:1:NeqRef
        CurCol = CurCol + 1;
        H(:,CurCol) = col_arith(OutRef(Isim),InPar.CalibEqRef{IeqRef},'mat');
    end
    for IeqCat=1:1:NeqCat
        CurCol = CurCol + 1;
        H(:,CurCol) = col_arith(Sim(Isim),InPar.CalibEqCat{IeqCat},'mat');
    end   

    % 1st iteration
    [Par,ParErr] = lscov(H(FlagAll,:),YY(FlagAll),1/VarYY(FlagAll));
    Par = H(FlagAll,:)\YY(FlagAll);
    Resid = YY - H(:,:)*Par;
    RStd  = Util.stat.rstd(Resid(FlagAll));
    
    Flag.GoodResid = abs(Resid) < RStd.*3;
    FlagAll = FlagAll & Flag.GoodResid;
    
    % 2nd iteration
    [Par,ParErr] = lscov(H(FlagAll,:),YY(FlagAll),1/VarYY(FlagAll));
    Par = H(FlagAll,:)\YY(FlagAll);
    Resid = YY(FlagAll) - H(FlagAll,:)*Par;
    RStd  = Util.stat.rstd(Resid);
    
    Chi2  = sum(Resid.^2./VarYY(FlagAll));
    Nstar = numel(Resid);
    
    ResPhot(Isim).Par      = Par;
    ResPhot(Isim).ParErr   = ParErr;
    ResPhot(Isim).Resid    = Resid;
    ResPhot(Isim).rstd     = Util.stat.rstd(Resid);
    
end

