function [Sim,FitRes,Log,ReturnFlag]=astrometry(Sim,varargin)
% Find astrometric solution to a SIM image
% Package: class/@SIM
% Description: Find an astrometric solution to a SIM image. Given an image
%              with roughly known scale, rotation and location, match the
%              sources in the image with external catalog by pattern
%              recognition and apply astrometric solution to the image.
%              Comments: The function can apply the solution to different
%              part of the image, recognize bad solution, and merge only
%              the good solutions. Proper motions are applied to the
%              refernce catalog.
% Input  : - A Sim object.
%          * Arbitrary number of pairs of arguments: ...,keyword,value,...
%            where keyword are one of the followings:
%            --- Additional parameters
%            Any additional key,val, that are recognized by one of the
%            following programs:
% Output : - 
% License: GNU general public license version 3
% Tested : Matlab R2014a
%     By : Eran O. Ofek                    Mar 2015
%    URL : http://weizmann.ac.il/home/eofek/matlab/
% Example: 
% Reliable: 
%--------------------------------------------------------------------------

RAD = 180./pi;

DefV.RePopCat           = false;
DefV.SrcExtractorProg   = @mextractor;
DefV.SrcExtractorPar    = {};

DefV.Scale              = 1.01; % arcsec/pix
DefV.RA                 = [];   % rad,...
DefV.Dec                = [];   % rad,...
DefV.Equinox            = 2000.0;
DefV.Scale              = 1.01;
DefV.KeyScale           = 'SCALE'; %
DefV.KeyRA              = {'RA','OBJRA','OBJRAD','CRVAL1'};
DefV.KeyDec             = {'DEC','OBJDEC','OBJDECD','CRVAL2'};
DefV.KeyEquinox         = {'EQUINOX'};
DefV.UnitsRA            = [];       % [], 'deg','r'
DefV.UnitsDec           = [];       % [], 'deg','r'
DefV.RefCat             = @get_ucac4; %@wget_ucac4;   % string, function, struct
DefV.RCrad              = 1.4/RAD;   % [radian]
DefV.Shape              = 'box';
DefV.RC_ColRA           = 'RA';
DefV.RC_ColDec          = 'Dec';
DefV.RC_ColMag          = 'MagG'; %'MagModel'; %'ModelMag';
DefV.RC_ColColor        = '';
DefV.ApplyPM            = false; %true;
DefV.RC_ColPM_RA        = 'pmRA';
DefV.RC_ColPM_Dec       = 'pmDec';
DefV.RC_ColPlx          = 'Plx';
DefV.RC_ColRV           = 'RV';
DefV.RC_EpochInRA       = 'EpochRA';
DefV.RC_EpochInDec      = 'EpochDec';
DefV.RC_EpochInUnits    = 'yr';
DefV.ColXc              = 'XWIN_IMAGE';
DefV.ColYc              = 'YWIN_IMAGE';
DefV.ImSize             = [];  % [x,y]
DefV.Flip               = [1 -1]; %1 1];
DefV.VecRot             = 0; %(-2:0.2:2).';
DefV.ReFindMatches      = true;
DefV.SearchRad          = 2;   % pix
DefV.SearchMethod       = 'binms';
DefV.SearchRangeX       = [-3000 3000];
DefV.SearchRangeY       = [-3000 3000];
DefV.SearchStepX        = 2;
DefV.SearchStepY        = 2;
DefV.MaxMethod          = 'adaptive'; %'sort';       % {'max1'|'sort'|'
DefV.MinNinPeak         = 5;            % min number of matches in peak - if <1 then fraction of ListCat
DefV.MaxPeaks           = 10;           % maximum number of peaks to return
DefV.MaxPeaksCheckI     = 5;     % 10      % out of the returned peaks - this is the maximum nuber of peaks to check
DefV.MaxPeaksCheckF     = 5;     % 10      % out of the returned peaks - this is the maximum nuber of peaks to check
DefV.CooType            = 'plane';
DefV.SelectBest         = 'comb';    % {'N','std','meanerr','comb'}
DefV.FalseAlarm         = 1e-4;
DefV.RegionMaxConn      = 8;
DefV.AdjustStarDensity  = true;
DefV.AdjustRatio        = 2;           % max ratio between Cat and Ref density
DefV.CatColMag          = 'MAG_PSF';  % Mag column name in SIM
DefV.BlockSize          = 'full'; %[1024 1024]; %'full';
DefV.BufferSize         = 200;
DefV.AnalysisBlockSize  = [1024 1024];
InPar = InArg.populate_keyval(DefV,varargin,mfilename);


% Internal definitions
OrigColXY = {'dup_X','dup_Y'};
CatField  = AstCat.CatField;

% Check if catalog exist in SIM
IsCatPop = isfield_populated(Sim,CatField);
if (~all(IsCatPop))
    % some of Cat fields are not populated
    InPar.RePopCat = true;
end
if (InPar.RePopCat)
    Sim = InPar.SrcExtractorProg(Sim,InPar.SrcExtractorPar{:});
end


%-----------------------------------------------
%--- Convert user supplied coo user to J2000 ---
%-----------------------------------------------
if (~isempty(InPar.RA) && ~isempty(InPar.Dec))
    % User supplied RA/Dec
    RA  = celestial.coo.convertdms(InPar.RA,'gH','r');
    Dec = celestial.coo.convertdms(InPar.Dec,'gD','R');
    if (InPar.Equinox~=2000)
        Coo2000 = coco([RA,Dec],sprintf('j%06.1f',InPar.Equinox),'j2000.0');
        RA  = Coo2000(:,1);
        Dec = Coo2000(:,2);
    end    
else
    % Read RA/Dec from header
    Out     = getcoo(Sim,'KeyRA',InPar.KeyRA,'KeyDec',InPar.KeyDec,'KeyEquinox',InPar.KeyEquinox,'OutUnits','rad');
    Coo2000 = celestial.coo.coco([[Out.RA].',[Out.Dec].'],sprintf('j%06.1f',Out(1).Equinox),'j2000.0');
    RA  = Coo2000(:,1);
    Dec = Coo2000(:,2);
end

if (~isempty(InPar.ImSize))
    ImSize = InPar.ImSize;
else
    ImSize = imagesize(Sim);
end


%--------------------------
%--- Reference Catalaog ---
%--------------------------

if (AstCat.isastcat(InPar.RefCat) || isstruct(InPar.RefCat))
    % RefCat was provided
    RefCat = InPar.RefCat;
    
    % convert to AstCat object
    RefCat = struct2astcat(RefCat);
else
    % External catalog was not provided
    % try to retrieve
    
    % get catalog data around image coordinates
    %RefCat = VO.search.get_cat(InPar.RefCat,RA(1),Dec(1),InPar.RCrad,'OutType','astcat','Shape',InPar.Shape);
    
    %RefCat = VO.search.get_gaia(RA(1),Dec(1),InPar.RCrad,'SearcRadUnits','rad');
    RefCat = VO.search.cat_cone('GAIADR1',RA(1),Dec(1),InPar.RCrad,'RadiusUnits','rad','OutType','astcat');
end

if (isempty(RefCat.Cat))
    % No reference catalog found - astrometric solution failed
    warnning('Reference catalog is empty - no solution found');
    
    FitRes     = [];
    Log        = [];
    ReturnFlag = false;
    
end
    


%---------------------------
%--- apply proper motion ---
%---------------------------
if (InPar.ApplyPM)
    % applay proper motion, RV and parallax to star positions
    RefCat = apply_proper_motion(RefCat,'EpochInRA',InPar.RC_EpochInRA,...
                                        'EpochInDec',InPar.RC_EpochInDec,...
                                        'EpochInUnits',InPar.RC_EpochInUnits,...
                                        'ColPM_RA',InPar.RC_ColPM_RA,...
                                        'ColPM_Dec',InPar.RC_ColPM_Dec,...
                                        'ColPlx',InPar.RC_ColPlx,...
                                        'ColRV',InPar.RC_ColRV);
                                
    
end




% Generate a version of the reference catalog with only selected columns
% RA, Dec, Mag, Color
[RC_ColRA, RC_ColDec] = colname2ind(RefCat,{InPar.RC_ColRA, InPar.RC_ColDec});
RC                = col_arith(RefCat,{InPar.RC_ColRA,InPar.RC_ColDec,...
                                          InPar.RC_ColMag,InPar.RC_ColColor},...
                              'astcat',true);
RC.ColCell        = {InPar.RC_ColRA,InPar.RC_ColDec,'Mag','Color'};
RC                = colcell2col(RC);

% calculate average density of sources in RefCat (RC):
% this is used to adjust the number of stars in the catalog and reference.
DensityRef = size(RC.Cat,1)./(pi.*(InPar.RCrad.*RAD).^2);   % stars/deg^2


% project catalog coordinates from celestial sphere to plane
% X,Y in units of pixels (using first guess pixel scale).
% Note that the X=0,Y=0 position is the position defined
% by the guess RA/Dec (RA, Dec parameters).
% This are the "projection-plan coordinates"
% relative to the native long and lat of the fiducial point
[X,Y]=projection(RC,'tan',[RC_ColRA RC_ColDec],[RAD.*3600./InPar.Scale RA Dec],'rad');
% add the projected X/Y as the 1st and 2nd columns in the reference catalog
ColRef   = {'X','Y'};
RC = col_insert(RC,X,1,'X');
RC = col_insert(RC,Y,2,'Y');
% sort the reference catalog by the Y position
RC = sortrows(RC,'Y');




Nsim       = numel(Sim);
FitRes     = Util.struct.struct_def({'Par','ParErr','Resid','ResidX','ResidY','ResidT',...
                         'Rrms','Chi2','Nobs','Npar','Ndof','rms','AstT','Flag','rmsAll','RrmsAll','Wrms','MagRef','Xref','Yref'},Nsim,1);
Log        = Util.struct.struct_def({'PatternMatch','FitMerge','FitFinal','Analysis'},Nsim,1);
ReturnFlag = false(Nsim,1);
for Isim=1:1:Nsim
    % for each SIM image
    
    % check if a new reference catalog query is needed
    %...
    
    
    % adjust number of stars in Sim such that it will be roughly equal to
    % the star density in the reference image
    CatColMagInd = colname2ind(Sim(Isim),InPar.CatColMag);
    if (InPar.AdjustStarDensity)
        DensityCat = size(Sim(Isim).Cat,1)./prod(ImSize(Isim,:).*InPar.Scale./3600);
        Fraction = min(InPar.AdjustRatio.*DensityRef./DensityCat,1);
        Quant    = quantile(Sim(Isim).Cat(:,CatColMagInd),Fraction);
        Flag     = Sim(Isim).Cat(:,CatColMagInd)<Quant;
        Sim(Isim).Cat = Sim(Isim).Cat(Flag,:);
    end
    
    % get X and Y column indices in SIM
    [ColXc,ColYc] = colname2ind(Sim(Isim),{InPar.ColXc,InPar.ColYc});   

    % Construct a table with [X,Y] columns only
    
    SimCat = AstCat.sim2astcat(Sim(Isim));
    %SimCat  = col_select(SimCatA,[ColXc ColYc]);
    
    
    % Duplicate X/Y columns in SimCat
    % Original X/Y are named: 'dup_X' and 'dup_Y'
    % transformed X/Y are '*WIN_IMAGE'...
    SimCat     = col_duplicate(SimCat,[ColXc,ColYc],OrigColXY);
    
    % select sources in sub image 
    % If the solution is done in sub image then this generate
    % a catalog for each sub image (block).
    SubCat     = subcat_regional(SimCat,ImSize(Isim,:),InPar.BlockSize,InPar.BufferSize,[ColXc,ColYc]);
    Nsub       = numel(SubCat);
    MatchedCat = AstCat(Nsub,1);
    MatchedRef = AstCat(Nsub,1);
    ResBest    = Util.struct.struct_def({'MaxHistMatch','MaxHistShiftX','MaxHistShiftY',...
                        'ShiftX','ShiftY','Tran','Nmatch','IndRef','IndCat',...
                        'MatchedCat','MatchedRef','MatchedResid','StdResid',...
                        'Std','MeanErr','BestRot','BestFlip'},Nsub,1);
    for Isub=1:1:Nsub
        % for each sub region
        
        %InPar.VecRot = 1;
    
        % match ref catalog with image catalog
        % find also the correct image flip and rotation.
        % The origin of the reference catalog (RC) is [RA, Dec].
        % The origin of the image is its geometric center.
        % Note that the output MatchedCat
        % contains the columns: {'XWIN_IMAGE'  'YWIN_IMAGE'  'dup_X'  'dup_Y'}
        [MatchedCat(Isub),MatchedRef(Isub),Res,IndBest,~] = ImUtil.Im.pattern_match_shift_rot(SubCat(Isub),RC,...
                                                        'VecRot',InPar.VecRot,...
                                                        'Flip',InPar.Flip,...
                                                        'ImSize',ImSize(Isim,:),...
                                                        'MaxPeaksCheckI',InPar.MaxPeaksCheckI,...
                                                        'MaxPeaksCheckF',InPar.MaxPeaksCheckF,...
                                                        'MatchShiftPar',{'SearchRad',InPar.SearchRad,...
                                                                      'SearchMethod',InPar.SearchMethod,...
                                                                      'SearchRangeX',InPar.SearchRangeX,...
                                                                      'SearchRangeY',InPar.SearchRangeY,...
                                                                      'SearchStepX',InPar.SearchStepX,...
                                                                      'SearchStepY',InPar.SearchStepY,...
                                                                      'MaxMethod',InPar.MaxMethod,...
                                                                      'MinNinPeak',InPar.MinNinPeak,...
                                                                      'MaxPeaks',InPar.MaxPeaks,...
                                                                      'CooType','plane',...
                                                                      'SelectBest',InPar.SelectBest,...
                                                                      'FalseAlarm',InPar.FalseAlarm,...
                                                                      'RegionMaxConn',InPar.RegionMaxConn});


        if (~isempty(IndBest))
            ResBest(Isub) = Res(IndBest);
        else
            ResBest(Isub).MaxHistMatch = 0;
            ResBest(Isub).Nmatch       = 0;
        end
        % store info in log
        if (isempty(IndBest))
            Log(Isim).PatternMatch(Isub).MaxHistMatch = 0;
        else
            Log(Isim).PatternMatch(Isub).MaxHistMatch    = Res(IndBest).MaxHistMatch;
            Log(Isim).PatternMatch(Isub).MaxHistShiftX   = Res(IndBest).MaxHistShiftX;
            Log(Isim).PatternMatch(Isub).MaxHistShiftY   = Res(IndBest).MaxHistShiftY;
            Log(Isim).PatternMatch(Isub).Nmatch          = Res(IndBest).Nmatch;
            Log(Isim).PatternMatch(Isub).StdResid        = Res(IndBest).StdResid;
            Log(Isim).PatternMatch(Isub).Std             = Res(IndBest).Std;
            Log(Isim).PatternMatch(Isub).MeanErr         = Res(IndBest).MeanErr;
            Log(Isim).PatternMatch(Isub).BestRot         = Res(IndBest).BestRot;
            Log(Isim).PatternMatch(Isub).BestFlip        = Res(IndBest).BestFlip;
        end


    end
    
    
    % check if solution found
    if (any([ResBest.MaxHistMatch]>InPar.MinNinPeak))
        % possible solution found
        
        % check that BestFlip are consistent with each other
        %<------
        
        ReturnFlag(Isim) = true;
    
        % merge matched catalogs 
        MatchedCat = merge(MatchedCat);
        MatchedRef = merge(MatchedRef);

            % Deal with empty MatchedCat
            % <--------------
%OrigColXY
        
        if (prod(ResBest.BestFlip)<0),
            TranPlateTiltX = 'x_platetilt_m';
            TranPlateTiltY = 'y_platetilt_m';
        else
            TranPlateTiltX = 'x_platetilt';
            TranPlateTiltY = 'y_platetilt';
        end
            
        %AstT = cell2asttran({'x_shift',0;'y_shift',0;'x_rot',[1 0];'y_rot',[0 1];'x_platetilt_m',0;'y_platetilt_m',0});
     
        
        %<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
        
       % AstT = cell2asttran({'x_shift',0;'y_shift',0;'x_rot',[1 0];'y_rot',[0 1]}); %;TranPlateTiltX,0;TranPlateTiltY,0});
        AstT = cell2asttran({'x_shift',0;'y_shift',0;'x_rot',[1 0];'y_rot',[0 1];TranPlateTiltX,0;TranPlateTiltY,0});
        AstT = cell2asttran({'x_shift',0;'y_shift',0;'x_rot',[1 0];'y_rot',[0 1];TranPlateTiltX,0;TranPlateTiltY,0;'x_ofad',[0 0 0];'y_ofad',[0 0 0]});
     
     
        FitRes(Isim) = fit_transform(MatchedCat,MatchedRef,'ColCatXY',[ColXc, ColYc],'ColRefXY',ColRef,'AstT',AstT);
        
%XXXX
%FitRes(Isim) = fit_transform(MatchedCat,MatchedRef,'ColCatXY',OrigColXY,'ColRefXY',ColRef,'AstT',AstT);

% Verify that fit transformation works
%FitRes(Isim) = fit_transform(MatchedCat,MatchedCat,'ColCatXY',OrigColXY,'ColRefXY',[ColXc, ColYc],'AstT',AstT);

        %  store info in log
        if (nargout>1)
            Log(Isim).FitMerge.Par       = FitRes(Isim).Par;
            Log(Isim).FitMerge.ParErr    = FitRes(Isim).ParErr;
            Log(Isim).FitMerge.Rrms      = FitRes(Isim).Rrms;
            Log(Isim).FitMerge.rms       = FitRes(Isim).rms;
            Log(Isim).FitMerge.Nobs      = FitRes(Isim).Nobs;
            Log(Isim).FitMerge.Npar      = FitRes(Isim).Npar;
            Log(Isim).FitMerge.Ndof      = FitRes(Isim).Ndof;
            Log(Isim).FitFinal.MinX      = col_fun2scalar(MatchedCat,@min,OrigColXY(1));
            Log(Isim).FitFinal.MaxX      = col_fun2scalar(MatchedCat,@max,OrigColXY(1));
            Log(Isim).FitFinal.MinY      = col_fun2scalar(MatchedCat,@min,OrigColXY(2));
            Log(Isim).FitFinal.MaxY      = col_fun2scalar(MatchedCat,@max,OrigColXY(2));
        end

        % use best fit transformation to find more matches
        if (InPar.ReFindMatches),
            % re-project coordinates around new center of image
            %??
            
            [~,IndBestSub] = max([ResBest.Nmatch]);

            %CatFlip     = flip_center(SimCatA,ResBest(IndBestSub).BestFlip,ImSize(Isim,:).*0.5,[ColXc, ColYc]);
            CatFlip     = flip_center(SimCat,ResBest(IndBestSub).BestFlip,ImSize(Isim,:).*0.5,[ColXc, ColYc]);
            [CatTran,~] = transform(FitRes(Isim).AstT,CatFlip,'ColXY',[ColXc,ColYc]);
            [CatTran]   = sortrows(CatTran,ColYc);

            % Match transformed catalog to reference
            ResMatch = search_cat(CatTran,RC.Cat(:,1),RC.Cat(:,2),...
                                  'ColX',ColXc,...
                                  'ColY',ColYc,...
                                  'CooType','plane',...
                                  'SearchRad',InPar.SearchRad,...
                                  'SearchMethod',InPar.SearchMethod);

            IndRef = find([ResMatch.Nfound].'==1);
            IndCat = [ResMatch(IndRef).IndCat].';
            
            %IndCat = SI(IndCat);
            %MatchedCat = row_select(SimCatA,IndCat);
            %problem here - need to apply inverse transformation!!!
            
            MatchedCat = row_select(CatTran,IndCat);
            
            
            MatchedRef = row_select(RC,IndRef);

            % Final fit transformation
            %AstT = cell2asttran({'x_shift',0;'y_shift',0;'x_rot',[1 0];'y_rot',[0 1];'x_platetilt',0;'y_platetilt',0});
            
            %AstT = cell2asttran({'x_shift',0;'y_shift',0;'x_rot',[1 0];'y_rot',[0 1];'x_platetilt',0;'y_platetilt',0;'x_ofad',[];'y_ofad',[]});
            %FitRes(Isim) = fit_transform(MatchedCat,MatchedRef,'ColXY',ColRef,'AstT',AstT,'NormXY',1);
%            FitRes(Isim) = fit_transform(MatchedCat,MatchedRef,'ColCatXY',[ColXc, ColYc],'ColRefXY',ColRef,'AstT',AstT);
            

%            IndOrigColXY = colname2ind(MatchedCat,OrigColXY);
%            IndColRef    = colname2ind(MatchedRef,ColRef);
            
%             ParPX = polyfit(MatchedCat.Cat(:,IndOrigColXY(1)),MatchedRef.Cat(:,IndColRef(1)),1)
%             ParPY = polyfit(MatchedCat.Cat(:,IndOrigColXY(2)),MatchedRef.Cat(:,IndColRef(2)),1)
%             MatchedCat.Cat(:,IndOrigColXY(1)) + ParPX(2);
%             MatchedCat.Cat(:,IndOrigColXY(2)) + ParPY(2);
            
%            MatchedCat.Cat(:,IndOrigColXY(1)) = MatchedCat.Cat(:,IndOrigColXY(1)) - ImSize(Isim,1).*0.5;
%            MatchedCat.Cat(:,IndOrigColXY(2)) = MatchedCat.Cat(:,IndOrigColXY(2)) - ImSize(Isim,2).*0.5;

            % Return the Ref catalog to arcsec units (instead of pixels)
            ColRefInd = colname2ind(MatchedRef,ColRef);
            MatchedRef.Cat(:,ColRefInd) = MatchedRef.Cat(:,ColRefInd).*InPar.Scale;
            
            % find best fit transformation
            % This is the best fit transformation between the Reference
            % catalog with origin at [RA,Dec] and the catalog with original
            % x,y pixels.
            FitRes(Isim) = fit_transform(MatchedCat,MatchedRef,'ColCatXY',OrigColXY,'ColRefXY',ColRef,'AstT',AstT);


            % plot mag vs. total ast. resid
            %plot(MatchedRef.Cat(FitRes.Flag,5),FitRes.ResidT,'.')

            % plot X/Y vs. resid in X/Y
            %plot(MatchedRef.Cat(FitRes.Flag,1),FitRes.ResidX,'.') 
            %plot(MatchedRef.Cat(FitRes.Flag,1),FitRes.ResidY,'.') 
            %plot(MatchedRef.Cat(FitRes.Flag,2),FitRes.ResidX,'.') 
            %plot(MatchedRef.Cat(FitRes.Flag,2),FitRes.ResidY,'.') 

            % plot mag vs. mag
            % plot(MatchedCat.Cat(:,5),MatchedRef.Cat(:,5),'.')


            %[AstOut,AstUM]=match(RC,AstC

            %  store info in log
            if (nargout>1)
                Log(Isim).FitFinal.Par       = FitRes(Isim).Par;
                Log(Isim).FitFinal.ParErr    = FitRes(Isim).ParErr;
                Log(Isim).FitFinal.Rrms      = FitRes(Isim).Rrms;
                Log(Isim).FitFinal.rms       = FitRes(Isim).rms;
                Log(Isim).FitFinal.Nobs      = FitRes(Isim).Nobs;
                Log(Isim).FitFinal.Npar      = FitRes(Isim).Npar;
                Log(Isim).FitFinal.Ndof      = FitRes(Isim).Ndof;
                Log(Isim).FitFinal.MinX      = col_fun2scalar(MatchedCat,@min,OrigColXY(1));
                Log(Isim).FitFinal.MaxX      = col_fun2scalar(MatchedCat,@max,OrigColXY(1));
                Log(Isim).FitFinal.MinY      = col_fun2scalar(MatchedCat,@min,OrigColXY(2));
                Log(Isim).FitFinal.MaxY      = col_fun2scalar(MatchedCat,@max,OrigColXY(2));
                
                
                % analyze residuals
                
                %<<<<<<<<<<<<<<<< THIS IS WRONG
                % need to apply transformation...
                
                CatA = col_select(MatchedCat,[ColXc, ColYc, CatColMagInd],FitRes.Flag);
                RefA = col_select(MatchedRef,{'X','Y','Mag','Color'},FitRes.Flag);
                CatA.Cat(:,1) = CatA.Cat(:,1) + ImSize(Isim,1).*0.5;
                CatA.Cat(:,2) = CatA.Cat(:,2) + ImSize(Isim,2).*0.5;
                RefA.Cat(:,1) = RefA.Cat(:,1) + ImSize(Isim,1).*0.5;
                RefA.Cat(:,2) = RefA.Cat(:,2) + ImSize(Isim,2).*0.5;
                Log(Isim).Analysis = analyze_astrometric_resid(CatA,...
                                                    RefA,...
                                                    ImSize(Isim,:),...
                                                    InPar.AnalysisBlockSize);

            end
        
        end
        
        
        
    else
        ReturnFlag(Isim) = false;
    end
    
    
    %---------------------
    %--- Calculate WCS ---
    %---------------------
    
    FitRes(Isim).AstT(1).RefPoint = [RA,Dec];
    
    [~,HeadWCS] = asttran2wcs(FitRes(Isim).AstT);
    % delete old WCS keywords including SIP
    Sim(Isim) = delete_wcs(Sim(Isim));
    
    if (HEAD.ishead(Sim))
        Sim(Isim) = replace_key(Sim(Isim),HeadWCS);
    end
    
    % for debuging
    % Check that the WCS+SIP is consistent with the referemce catalog
    
    [TestRA,TestDec] = xy2coo(Sim(Isim),Sim(Isim).Cat(:,1),Sim(Isim).Cat(:,2));
    AstTest = AstCat;
    AstTest.Cat = [TestRA,TestDec];
    AstTest.ColCell = {'RA','Dec'};
    AstTest = colcell2col(AstTest);
    %[AstOut,AstUM]=match(Sim(Isim),[TestRA,TestDec],'CatColNames',[1 2])
    %[AstOut,AstUM]=match(Sim(Isim),AstTest,'CatColNames',[1 2],'CatName',{'A_'})
    %AstOut
    
    
    
end

                                
                                
                                
                                
                                
                                
                                
                                
                          


