function [Res,CorrObsSpec]=spec_fit_telluric(ObsSpec,TelluricTemplate,varargin)
%--------------------------------------------------------------------------
% spec_fit_telluric function                                        ImSpec
% Description: Given an observed spectrum and a Telluric template
%              (e.g., generated by spec_telluric_template.m) adjust the
%              Telluric depth (multiplicative) by airmass or by minimizing
%              the rms of the spectrum in the Telluric region.
% Input  : - Observed spectrum [Wavelength, Flux].
%          - Telluric Template. If this is a column vector then this is
%            the Telluric absorption factor,
%            assumed to have the same sampling as the observed spectrum.
%            If this is a two column matrix [Wavelength, Absorption_Factor]
%            then it will be resampled to the observed spectrum
%            wavelengths.
%          * Arbitrary number of pairs or arguments: ...,keyword,value,...
%            where keyword are one of the followings:
%            'TelluricFlag' - A vector of flags (true/false) indicating if
%                             a specific wavelength is in a Telluric
%                             region or not.
%                             (e.g., generated by spec_telluric_template.m)
%                             Default is empty.
%                             Used if FitTelluric=rms.
%            'AM'           - AirMass of observation. Default is empty.
%                             Needed only of FitTelluric=am.
%            'Alt'          - Altitude of observation, will be used only
%                             if the air mass is not provided.
%                             Default is empty.
%            'TelluricAM'   - AirMass of Telluric template observation.
%                             Default is empty.
%                             Needed only of FitTelluric=am.
%            'InterpMethod' - Interpolation method. Default is 'linear'.
%            'FitTelluric'  - Method to adjust/fit Telluric template:
%                             'AM'  - Scale by the air mass ratio between
%                                     object and Telluric.
%                             'LogAM'- Scale by log10 of air mass ratio
%                                     (Beer's law) between object and
%                                     Telluric.
%                             'rms' - Find best scale by minimizing
%                                     the RMS in the Telluric regions.
%                             'const'- Adjust by a constant given by
%                                     'Const'.
%            'Const'        - Constant adjustment for FitTelluric=Const.
%                             Default is 1.
%            'Factor'       - Vector of factors of depth adjustments
%                             to test. Default is (0.7:0.02:1.3).'.
%                             Used if FitTelluric=rms.
%            'StdSize'      - Size of the stdfilt1.m window for calculting
%                             the rms in the Telluric regions.
%                             Default is 5.
%            'MeanFun'      - Handle for the "mean" function to use in
%                             order to calculate the mean RMS.
%                             Default is @nanmedian.
%            'ShiftTelluric'- Attempt to shift Telluric in order to match
%                             small shifts in wavelength calibration.
%                             'no' - No shift.
%                             'xc' - Cross correlation (default).
%            'ShiftTelluric'- Vector of shift Telluric factors to match
%                             small shifts in wavelength calibration.
%                             Default is (0.9995:0.00001:1.0005).'.
% Tested : Matlab R2011b
%     By : Eran O. Ofek                    Mar 2014
%    URL : http://weizmann.ac.il/home/eofek/matlab/
% Example: 
% Reliable: 
%-----------------------------------------------------------------------------

DefV.TelluricFlag    = [];
DefV.AM              = [];
DefV.Alt             = [];
DefV.TelluricAM      = [];
DefV.InterpMethod    = 'linear';
DefV.FitTelluric     = 'rms';
DefV.Const           = 1;
DefV.Factor          = (0.7:0.001:1.3).';
DefV.StdSize         = 10;
DefV.MeanFun         = @nanmedian;
DefV.ShiftTelluric   = 'xc';
DefV.ShiftVec        = (0.9995:0.00001:1.0005).';
InPar = set_varargin_keyval(DefV,'y','use',varargin{:});

% get AirMass
if (isempty(InPar.AM)),
    if (~isempty(InPar.Alt)),
        InPar.AM = hardie(pi./2 - InPar.Alt);
    end
end


WaveVec = ObsSpec(:,1);

Size = size(TelluricTemplate);
if (Size(2)==1),
    % same sampling
else
    % assumes two columns in TelluricTemplate
    if (~isempty(InPar.TelluricFlag)),
        InPar.TelluricFlag = logical(interp1(TelluricTemplate(:,1),InPar.TelluricFlag,WaveVec,InPar.InterpMethod,'extrap'));
    end
    TelluricTemplate = interp1(TelluricTemplate(:,1),TelluricTemplate(:,2),WaveVec,InPar.InterpMethod);
end

if (all(TelluricTemplate==1)),
    % TelluricTemplate is flat - skip correction
    AdjustedTelluric = TelluricTemplate;
else
    
    switch lower(InPar.FitTelluric)
        case 'const'

            AdjustedTelluric = 1 + (TelluricTemplate - 1).*InPar.Const;
            AdjustedTelluric(AdjustedTelluric<0) = 0;

            Min = [];
        case {'am','logam'}
            % adjust Telluric depth according to AM instead of fitting
            if (isempty(InPar.AM) || isempty(InPar.TelluricAM)),
                error('AM or Alt must supplied when FitTelluric=am');
            end
            Factor = InPar.AM./InPar.TelluricAM;
            switch lower(InPar.FitTelluric)
                case 'logam'
                    Factor = 1 + log10(Factor)
                otherwise
                    % do nothing
            end
            AdjustedTelluric = 1 + (TelluricTemplate - 1).*Factor;
            AdjustedTelluric(AdjustedTelluric<0) = 0;

            Min = [];
        case 'rms'
            % adjust Telluric depth by minimizing rms in ObsSpec
            if (isempty(InPar.TelluricFlag)),
                error('TelluricFlag must supplied when FitTelluric=rms');
            end
            Nfac = numel(InPar.Factor);
            TestObsSpec = ObsSpec;
            RMS = zeros(Nfac,1);
            ExtendedTelluricFlag = runmean1(double(InPar.TelluricFlag),5)>0;
            ExtendedTelluricFlag = runmean1(double(TelluricTemplate<0.8),11)>0;
            for Ifac=1:1:Nfac,
                %Ifac=230
                AdjVecFactor = (TelluricTemplate(InPar.TelluricFlag) - 1).*InPar.Factor(Ifac) + 1;
                AdjVecFactor(AdjVecFactor<0) = 0;
                TestObsSpec(InPar.TelluricFlag,2) = ObsSpec(InPar.TelluricFlag,2)./AdjVecFactor;

                %RMS(Ifac) = InPar.MeanFun(stdfilt1(TestObsSpec(InPar.TelluricFlag,2),InPar.StdSize));
                %RMS(Ifac) = InPar.MeanFun(stdfilt1(TestObsSpec(ExtendedTelluricFlag,2)./ObsSpec(ExtendedTelluricFlag,2),InPar.StdSize));
                RMS(Ifac) = InPar.MeanFun(stdfilt1(TestObsSpec(ExtendedTelluricFlag,2),InPar.StdSize));
    %RMS(Ifac) = range(TestObsSpec(ExtendedTelluricFlag,2));

                RMS(Ifac) = sum((TestObsSpec(InPar.TelluricFlag,2) - InPar.MeanFun(TestObsSpec(InPar.TelluricFlag,2))).^2);
                %figure(2); graph(TestObsSpec)
                %ll
            end
            [~,MinInd]       = min(RMS);
            %MinInd
            AdjustedTelluric = TelluricTemplate;
            AdjustedTelluric(InPar.TelluricFlag) = (TelluricTemplate(InPar.TelluricFlag) - 1).*InPar.Factor(MinInd) + 1;
            AdjustedTelluric(AdjustedTelluric<0) = 0;
            %TelluricTemplate(InPar.TelluricFlag).*InPar.Factor(MinInd);

            Min.RMS    = RMS;
            Min.Factor = InPar.Factor;
        otherwise
            error('Unknown FitTelluric option');
    end

    switch lower(InPar.ShiftTelluric)
        case {'no','none'}
            % do nothing
        case 'xc'
            Nshift = numel(InPar.ShiftVec);
            ShiftRMS = zeros(Nshift,1);
            for Ishift=1:1:Nshift,
                ShiftedTelluric  = interp1(ObsSpec(:,1),AdjustedTelluric,ObsSpec(:,1).*InPar.ShiftVec(Ishift),InPar.InterpMethod);
                ShiftRMS(Ishift) = InPar.MeanFun(stdfilt1(ObsSpec(:,2)./ShiftedTelluric,InPar.StdSize));
            end
            %plot(InPar.ShiftVec,ShiftRMS)
            [~,MinInd] = min(ShiftRMS);
            AdjustedTelluric = interp1(ObsSpec(:,1),AdjustedTelluric,ObsSpec(:,1).*InPar.ShiftVec(MinInd),InPar.InterpMethod);

        otherwise
            error('Unknwon ShiftTelluric option');
    end
end



% correct spectrum
CorrObsSpec = ObsSpec;
CorrObsSpec(:,2) = ObsSpec(:,2)./AdjustedTelluric;
        
Res.AdjustedTelluric = AdjustedTelluric;
Res.Min              = Min;
Res.AM               = InPar.AM;
Res.TelluricAM       = InPar.TelluricAM;

