%% QDC calculation
% for actual day
function [qday, ActDay] = q_day(arr) 
% [qday, ActDay] = q_day(arr)
%
% The function calculates actual Qiuet Day Curve (QDC) for
% INPUT array (arr) of data.
% arr - is 1-minutes values for 30 days = 43200 elements
% OUTPUT
% qday - array of 1400 elements (1-day) QDC values
% ActDay - actual day in 30 days array
%
% Change list:
% - 'moving' added to smooth statement for "absolute deviations from trend"
% computation
%
% $Revision: 2.0 $

% $Id: q_day.m,v 2.0 2018/12/12 10:47:45 jnsnl Exp $

% Absolute deviations from trend
arr_d = abs(arr - smooth(arr, 120, 'moving'));

% Absolute gradients
arr_g = abs(gradient(arr));

% Shift 61 min period by 1 min at a time and find the maximum deviations
% from trend and maximum gradients
for i = 31:length(arr)-30
    % Maximum deviations from trend in the 61 min interval
    max_ds(i-30,1) = max(abs(arr_d(i-30:i+30)));
    % Maximum gradients in the 61 min interval
    max_g(i-30,1) = max(abs(arr_g(i-30:i+30)));
end%for i

% Initialize test parameters
%Quiet value limit
divq = 0;
% Selected quiet data day-matrix
day = [];
% Selected quiet time vector
act(1:50000) = NaN;
ac = 1;
% Quiet day covered
p = 0;
% Quiet day exists
flag = true;

% Find quiet data covering an entire "day" by iterating through quiet 
% limits from 2 nT to <40 nT in steps of 2 nT until a sufficient number of
% 1-min quiet segments is gained
while min(p) < 120 % Find atleast 120 quiet samples within 2 hour periodes
    % Increase quiet limit by 2 nT
    divq = divq + 2;
    if divq >= 40 % If quiet limit reach 40 nT
        % Quiet day not covered
        flag = false;
        break
    end%if actual
    clear p;
    
    % Find data within quiet limit
    q = NaN(length(arr), 1);
    i = find(max_ds < divq & max_g < divq);
    q(i+30) = arr(i+30);
    
    % Append quiet data to quiet data day-matrix
    day = [day, reshape(q, [1440 length(arr)/1440])];
    
    % Quiet minutes within quiet limit
    act(ac:ac+length(i)-1) = i+30;
    ac = ac + length(i);
    
    % Find number of quiet data-points within 2 hour periods to ensure that
    % quiet day is covered. Atleast 120 quiet samples within 2 hours
    for i = 1:1320
        p(i) = length(find(~isnan(day(i:i+120,:))));
    end%for i
end%while

if flag % Quiet day covered
    % Repeat quiet day and calculate mean for every quiet minute
    qqq = nanmean([day; day; day]');
    % Smooth Quiet Dayli Curve
    qqq_s  = smooth(qqq, 240, 'moving');
    qqq_s = smooth(qqq_s, 240, 'rlowess');
    % Quiet Dayli Curve covering one day
    qday = qqq_s(1440:2880-1);
    % Calculate Actual Day
    ActDay = round(nanmedian(act)/1440);
else
    % No Quiet Daily Curve found
    qday(1:1440) = NaN;
    % No Actual quiet day found
    ActDay = NaN;
end