    function [clearSamples, csGHI, alpha] = pvl_detect_clear_times(GHI, Time, UTCoffset, Location, win_length, sample_interval)
    % PVL_DETECT_CLEAR_TIMES identifies times with irradiance consistent with
    % clear sky conditions in a time series of GHI.
    %
    % Syntax
    %   [clearSamples, csGHI, alpha] = pvl_detect_clear_times(GHI, Time, UTCoffset, Location, win_length, sample_interval)
    %
    % Description
    %   Detects clear sky times by comparing statistics for a regular GHI
    %   time series to the Ineichen clear sky model.  Statistics are calculated
    %   using a sliding time window (e.g., 10 minutes). An iterative algorithm
    %   identifies clear periods, uses the identified periods to estimate bias
    %   in the clear sky model and/or GHI data, adjusts the clear sky model and
    %   repeats.  Code handles GHI data with some irregularities, i.e., missing
    %   values or unequal data spacing, but execution is significantly faster
    %   if data are equally spaced and without missing values.
    %
    %   If non-daylight times are included the corresponding GHI values should
    %   be NaNs.
    %
    %   Clear times are identified by meeting 5 criteria, thresholds
    %   for which are hardcoded in this version.  Values for these thresholds
    %   are appropriate for 10 minute windows of 1 minute GHI data.
    %
    % Inputs:
    %   GHI - a Nx1 vector of GHI values (W/m2)
    %   Time - a Nx1 vector of datenum with equal time spacing
    %   UTCoffset - scalar UTC offset (e.g. EST = -5)
    %   Location - a struct with the following elements, note that all
    %      elements must be scalars in this application.
    %   Location.latitude - scalar latitude in decimal degrees (positive is
    %      northern hemisphere)
    %   Location.longitude - scalar longitude in decimal degrees (positive is
    %      east of prime meridian)
    %   Location.altitude - scalar height above sea level in meters.
    %      While altitude is optional in many uses, it is required in this
    %      model implementation.
    %   win_length - length of sliding time window in minutes, integer
    %   sample_interval - nominal minutes between each GHI sample time, integer
    %
    % Output:
    %   clearSamples - column vector of logical values with True indicating
    %    a clear sample in the input GHI time series
    %   csGHI - column vector with scaled clear sky GHI derived from the
    %   Ineichen model
    %   alpha - scaling factor applied to Ineichen model to obtain output csGHI
    %
    % References
    %   [1] Reno, M.J. and C.W. Hansen, "Identification of periods of clear sky
    %   irradiance in time series of GHI measurements" Renewable Energy, v90,
    %   p. 520-531, 2016.
    %
    % Notes:
    %   Initial implementation by Matthew Reno. Modifications for computational
    %   efficiency by Joshua Patrick and Curtis Martin.
    %

    % Threshold values for each criterion.  Values are appropriate for 10
    % minute windows with 1 minute intervals between GHI data.

    % MeanDiff: threshold value in W/m2 for first criterion, i.e., agreement between mean
    % values of GHI in each interval, see Eq. 6 in [1]
    MeanDiff = 75;

    % MaxDiff: threshold value in W/m2 for second criterion, i.e., agreement
    % between maxima of GHI values in each interval, see Eq. 7 in [1]
    MaxDiff = 75;

    % LineLengthLower and LineLengthUpper: threshold values (unit is weird) for
    % third criterion, i.e., agreement between line lengths of GHI data and
    % clear sky GHI, see Eq. 8 in [1].  Criterion is satisfied when line length
    % LL meets LineLengthLower <= LL <= LineLengthUpper
    LineLengthLower = -5;
    LineLengthUpper = 10;

    % VarDiff: threshold value in 1/seconds for the fourth criterion, i.e.,
    % agreement between normalized standard deviations of rate of change in
    % irradiance, see Eq. 9 through Eq. 11 in [1]
    VarDiff = 0.005;

    % SlopeDev: threshold value in W/m2 for the fifth criterion, i.e.,
    % agreement between largest magnitude of change in successive GHI values,
    % see Eq. 12 through 14 in [1]
    SlopeDev = 8;

    dv = datevec(Time);
    dates = datenum(dv(:,1:3));
    days = unique(dates);
    ivec = (1:length(Time))';

    % check for equally spaced data, code execution is much faster with equal
    % time spacing

    sTime = diff(round(Time*1440)); % minutes between samples
    if length(unique(sTime)) == 1 && mod(win_length, sample_interval)==0
        % equally spaced data, none missing
        equal_spacing = true;
    else
        equal_spacing = false;
    end

    % get sun angles for determining sunrise/sunset
    [~, SunEl, ~, ~] = pvl_ephemeris(pvl_maketimestruct(Time,UTCoffset), Location);
    daylight = SunEl > 1;

    % Get standard Ineichen clear sky GHI
    csGHI0 = pvl_clearsky_ineichen(pvl_maketimestruct(Time,UTCoffset),Location);
    % replace NaNs with 0
    csNaNs = isnan(csGHI0);
    csGHI0(csNaNs) = 0;

    % Create index matrix for windows

    max_samples_per_window = 2 * win_length / sample_interval;  % upper bound on number of data samples in an interval
    min_samples_per_window = floor(0.8 * win_length / sample_interval) ; % lower bound on number of data samples in an interval

    if not(equal_spacing)
        % pre-allocate some arrays
        clearMean = NaN(size(GHI));
        clearMax = clearMean;
        clearMaxSlope = clearMean;
        timeDiff = zeros(max_samples_per_window, length(clearMean));
        clearGHIDiff = timeDiff;

        measuredMean  = clearMean;
        measuredMax  = clearMean;
        measuredMaxSlope  = clearMean;
        measuredSqSlope  = clearMean;
        measuredStdSlope  = clearMean;
        measuredLineLength = clearMean;

        k = 0; % counts number of intervals meeting data size criteria, i.e., samples between min_samples_per_window and max_samples_per_window
        for j = 1:length(Time)
            tu = Time>=Time(j) & Time<(Time(j)+ (win_length + 1)/1440);
            tidx = ivec(tu);
            tmp_GHI = GHI(tidx);
            if sum(daylight(tu)) < min_samples_per_window
                % sun is not above horizon, skip interval
            elseif sum(~isnan(tmp_GHI)) >= min_samples_per_window && ...
                    sum(~isnan(tmp_GHI)) <= max_samples_per_window  % e.g., at least 8 measurements within a window but not more than 20
                k = k + 1;
                tmp_csGHI0 = csGHI0(tidx);
                w = diff(Time(tidx)*1440);  % time intervals in minutes between samples

                % Calculate parameters (all but line length) for alpha=1 case
                wx = [w; 0]/sum(w); % leave last point out of mean
                clearMean(j) = sum(wx.*tmp_csGHI0);
                clearMax(j) = max(tmp_csGHI0);
                clearMaxSlope(j) = max(abs(diff(tmp_csGHI0,1,1)./w));
                tmp_diffGHI = diff(tmp_csGHI0);
                clearGHIDiff(1:length(tmp_diffGHI),j) = tmp_diffGHI;
                timeDiff(1:length(tmp_diffGHI),j) = w;

                % Calculate parameters for measured GHI
                measuredMean(j) = sum(wx.*tmp_GHI);
                measuredMax(j) = max(tmp_GHI);
                tmpSlope = diff(tmp_GHI,1,1)./w;
                measuredMaxSlope(j) = max(abs(tmpSlope));
                measuredSqSlope(j) = sum(tmpSlope.^2);
                measuredStdSlope(j) = std(tmpSlope);
                measuredLineLength(j) = sum(sqrt(diff(tmp_GHI,1,1).^2 + w.^2));
            else
                % too much missing data or too many data samples
                disp(['Data problem between ' datestr(min(Time(tu))) ' and ' datestr(max(Time(tu))) ...
                    ' have ' num2str(sum(tu)) ' values']);
            end
        end

    else % equal data spacing

        samples_per_window = win_length / sample_interval;
        H = hankel(1:samples_per_window,samples_per_window:length(Time));
        clearMean = mean(csGHI0(H));
        clearMax = max(csGHI0(H));
        clearSlope = diff(csGHI0(H),1,1);

        % Calculate parameters for measured GHI
        measuredMean = mean(GHI(H));
        measuredMax = max(GHI(H));
        measuredSlope = diff(GHI(H),1,1);
        measuredLineLength = sum(sqrt(measuredSlope.^2 + (sample_interval*ones(size(measuredSlope))).^2 ));

    end

    % Normalized variance of measured GHI < Vardeiff criterion doesn't depend on
    % clear sky GHI, compute outside of iteration on CS model adjustment alpha
    if not(equal_spacing)
        c4 = measuredStdSlope ./ measuredMean < VarDiff;
    else
        c4 = std(measuredSlope) ./ measuredMean < VarDiff;
    end

    % initialize state for iterative process
    alpha = 1;      %start off with Ineichen clear sky model

    % Start Iterative process of finding clear days, fitting clear sky model, finding clear days, ....
    % Continue iterations until not finding any more clear days or 20 iterations
    for ii=1:20

        % Update clear sky model by scaling Ineichen
        csGHI = alpha * csGHI0;

        % Adjust Clear Sky Model parameters for scaling by alpha
        if not(equal_spacing)
            clearLineLength = sum(sqrt((alpha*clearGHIDiff).^2 + timeDiff.^2 ));
            clearLineLength = clearLineLength(:);
        else
            clearLineLength = sum(sqrt((alpha*clearSlope).^2 + (sample_interval*ones(size(measuredSlope))).^2 ));
            measuredMaxSlope = max(abs(measuredSlope));
            clearMaxSlope = max(abs(clearSlope));
        end

        % Evaluate comparison criteria
        c1 = abs(measuredMean - alpha * clearMean) < MeanDiff;
        c2 = abs(measuredMax - alpha * clearMax) < MaxDiff;
        c3 = (LineLengthLower < measuredLineLength - clearLineLength) & (measuredLineLength - clearLineLength < LineLengthUpper) ;
        c5 = (measuredMaxSlope - alpha * clearMaxSlope) < SlopeDev;
        c6 = clearMean ~= 0 & ~isnan(clearMean);            % window includes some daylight

        % Identify windows meeting all five (six) criteria
        clearWindows = c1 & c2 & c3 & c4 & c5 & c6;

        % Daily clearness is proportion of samples that are clear
        dailyClearness = histcounts(Time(clearWindows),[days]) ./ ...
            histcounts(Time(csGHI>25),[days]);

        % Save current state
        previousAlpha = alpha;

        % Get new alpha by adjusting clear sky model to match clear times
        % Need to minimize RMSE between GHI and alpha * csGHI0

        clearSamples = false(size(GHI));

        if not(equal_spacing)
            for j = 1:length(Time)
                if clearWindows(j)
                    tu = Time>=Time(j) & Time<(Time(j)+win_length/1440);
                    clearSamples(tu) = true;
                end
            end
        else
            clearSamples = ismember(ivec,H(:,clearWindows));
        end

        RMSE = @(x)sqrt(mean((GHI(clearSamples) - x*csGHI0(clearSamples)).^2));
        alpha = fminsearch(RMSE, alpha);  % update scaling factor on clear sky model

        if round(alpha*10000)==round(previousAlpha*10000)
            previousAlpha = alpha; %#ok<NASGU>
            break
        end

    end  % for loop

    csGHI = alpha * csGHI0;

    end  % function

    Not enough input arguments.

    Error in pvl_detect_clear_times (line 87)
    dv = datevec(Time);