function QL_PC_Index
%% Calculation of Quick Look PC index
% Compute Quick Look PC-index for the northern hemisphere (thl).
% Results are written to a SQL database defined in db_config.pc
%
% SQL table information
% +-----------+--------------+------+-------------------+-------------------------+
% | COLUMN    | TYPE         | NULL | DEFAULT           | COMMENTS                |
% +-----------+--------------+------+-------------------+-------------------------+
% | id        | int(11)      | No   |                   | SQL auto increment id   |
% +-----------+--------------+------+-------------------+-------------------------+
% | timestamp | timestamp    | No   | CURRENT_TIMESTAMP | SQL entry time          |
% +-----------+--------------+------+-------------------+-------------------------+
% | time      | datetime     | No   |                   | time of the measurement |
% +-----------+--------------+------+-------------------+-------------------------+
% | x         | float        | Yes  | NULL              | X component             |
% +-----------+--------------+------+-------------------+-------------------------+
% | y         | float        | Yes  | NULL              | Y component             |
% +-----------+--------------+------+-------------------+-------------------------+
% | qd_x      | float        | Yes  | NULL              | QD(x)                   |
% +-----------+--------------+------+-------------------+-------------------------+
% | qd_y      | float        | Yes  | NULL              | QD(y)                   |
% +-----------+--------------+------+-------------------+-------------------------+
% | pc        | float        | Yes  | NULL              | PC index                |
% +-----------+--------------+------+-------------------+-------------------------+
%
% $Revision: 1.0 $

% $Id: QL_PC_Index.m,v 1.0 2019/03/14 15:25:52 jnsnl Exp $

%% Input
% Path to station nrt second cdf-files
path_cdf = '/m1/projects/magobs/measurements/magdata/winsys/data/THL6/cdf';

% Longitude of station (291 for THL and 106.9 for VOS)
lon = 291;

% Margin for preliminary data +/- n days
marg = 30;

% Load load phi, alpha and beta coefficients
load('coeff.mat');

% Read config file db_cong.pc
f = fopen('db_config.pc','r');
    serv = fgetl(f);
    log = fgetl(f);
    pass = fgetl(f);
    dbname  = fgetl(f);
    table = fgetl(f);
fclose(f);

%% Initialize data structure
data  = struct('t', [], 'x', [], 'y', [], 'qd_x', [], 'qd_y', [], 'pc', []);

mysql('open', serv,log,pass);
mysql(['use ' dbname]);
% Find start time from latest timestamp from QL datebase
t1 = mysql(['select time from ' table ' where id=(select max(id) from ' ...
            table ');']) + 1/1440;

% Get timestamp from latest nrt second cdf-files
cdf_files = dir(fullfile(path_cdf, 'THL6_*.cdf'));
t2 = datenum(cdf_files(end).name(6:13),'yyyymmdd');

% Number of days
ndays = t2-floor(t1)+marg*2+1;
% Write time in minutes to data structure
[y, m, d, H, M] = datevec(floor(t1));
data.t = (datenum(y, m, d-marg*2, H, M:ndays*1440-1, 0))';
% Initialize the remaining data structure with database data and NaN's
fns = fieldnames(data);
for fn = fns(2:end)' % [x, y, qd_x, qd_y, pc]
    data.(char(fn)) = NaN(length(data.t), 1);
    i = (data.t >= t1-marg*2 & data.t < t1);
    data.(char(fn))(i) = mysql(['select ' char(fn) ' from ' table ...
        ' where time >= "' datestr(t1-marg*2,'yyyy-mm-dd HH:MM:SS') '"' ...
        ' and time <= "' datestr(t1-1/1440,'yyyy-mm-dd HH:MM:SS') '"' ...
        'order by time ;']);
end%for fn
mysql('close');

% Produce quasi-definitive minute XY data from second cdf-files from start
% date until now and write to data structure
for cdf_date = floor(t1):t2
    i = find(data.t >= cdf_date & data.t < cdf_date+1);
    [~, data.x(i), data.y(i)] = cdf2min(cdf_date, path_cdf);
end%for cdf_date

% Timestamp of latest data sample
t2 = data.t(find(~isnan(data.x+data.y), 1, 'last'));
if(t1 >= t2)
  % No new data
  return
end%if

%% Calculate QL PC index
% Initialize
qdc = struct('x', NaN(length(data.t), 1), 'y', NaN(length(data.t), 1));
% If the gap between t1 and t2 are larger than one day, fill in in one 
% day increments
n = marg*2+1;
for t = [floor(t1):floor(t2)-1, t2]
%% Find Quiet Days (Actual Days)
for fn = fns(2:3)' %[x, y]
    % Update once evrey day at midnight to calculate a series of actual QD:
    % (If interpolation and smoothing Actual QDC for every day fails, try
    %  recalculate actual QD)
    for m = 1:2
        if(t2-floor(t) > 1 || m == 2)
          % Calculate QDC for 30-days dataset and save it into
          % "Actual Day" in data structure
          start_dt = floor(t)-30;
          end_dt = floor(t)-1/1440;
          i = (data.t >= start_dt & data.t <= end_dt);
          
          % Calculate QDC for an "Actual Day"
          [qday, ActDay] = q_day(data.(char(fn))(i));
          
          % Write QDC for "Atcual Day" to data structure
          if ~isnan(ActDay) % If any
            j = data.t >= start_dt + ActDay - 1 & ...
                data.t <= start_dt + ActDay - 1/1440;
            if(any(j))
              data.(['qd_' char(fn)])(j) = qday;
            end%if
          end%if ~isnan
        end%if
        
        % Interpolation of QDC for every day from "Actual Day" array
        k = (data.t >= floor(t)-marg*2 & data.t <= floor(t)+1-1/1440);
        qd = data.(['qd_' char(fn)])(k);
        % "Actual Days"
        ActDays = find(~isnan(sum(reshape(qd, [1440 n]))));
        
        % Construct 2D array of QDC
        qdays = reshape(qd, [1440 n]);
        
        % Interpolation and smoothing Actual QDC for every day
        interp_arr = NaN(1440, n);
        try
        for i = 1:1440
            interp_arr(i,:) = smooth(interp1(ActDays, qdays(i,ActDays), 1:n, ...
                                     'nearest', 'extrap'), 60, 'lowess');
        end
        break
        catch
        if(m == 2)
          error('Faild to interpolate and smooth Actual QDC for every day')
        end%if
        end%try
    end%for m
    
    % Make 1D array of the QDC and smooth and write to data structure
    l_arr = reshape(interp_arr, [1 n*1440]);
    qdc.(char(fn))(k) = smooth(l_arr,120,'loess');
end%for fn

%% Caculate PC index
i = (data.t >= t-marg & data.t <= t);
% Calculate disturbed part of the variations
dist_x = data.x(i) - qdc.x(i);
dist_y = data.y(i) - qdc.y(i);

% H Projection and PC calculation for the period
[y, m, d, H, M] = datevec(data.t(i));
% Day of year since 1. january
doy = (datenum(y, m, d) - datenum(y, 1, 1));
% Minute in the year
min_in_year = doy*1440 + H*60 + M + 1;
% Calculate H_proj
UT = H .* 15 + M .* 0.25; % TIME to DEG
y = lon + coeff.f(min_in_year) + UT;
H_proj = dist_x.*sin(y.*pi./180) - dist_y.*cos(y.*pi./180);
% Calculate PC index and write to data structure
data.pc(i) = (H_proj - coeff.b(min_in_year))./coeff.a(min_in_year);

end%for t

%% Write data structure to SQL database
mysql('open', serv, log, pass);
mysql(['use ' dbname]);
% Index of t1
i_t1 = find(data.t >= t1-1/2880 & data.t <= t1+1/2880);
% Index of t2
i_t2 = find(data.t >= t2-1/2880 & data.t <= t2+1/2880);
% Write data to SQL database in bulks of 1000 rows (max limit of SQL
% INSERT-statement) to save time
bulk = [i_t1-marg*1440:1000:i_t2, i_t2];
for i = 1:length(bulk)-1
    % Construct SQL query string
    str = [];
    for j = bulk(i)+1:bulk(i+1)
        str = [str '('];
        % Format time for SQL database
        str = [str '''' datestr(data.t(j),'yyyy-mm-dd HH:MM:SS') ''''];
        
        % Format data for SQL database
        for fn = fns(2:end)' % [x, y, qd_x, qd_y, pc]
            if ~isnan(data.(char(fn))(j))
              % Save data with two decimals in SQL database
              str = [str ', ' num2str(round(data.(char(fn))(j)*100)/100)];
            else
              % NaN is set to NULL
              str = [str ', NULL'];
            end%if ~isnan
        end%for fn
        str = [str '), '];
    end%for j
    % SQL query
    query = ['INSERT INTO ' table ' (time, x, y, qd_x, qd_y, pc) ' ...
        'VALUES ' str(1:end-2) ' '...
        'ON DUPLICATE KEY UPDATE x = VALUES(x), y = VALUES(y), ' ...
        'qd_x = VALUES(qd_x), qd_y = VALUES(qd_y), ' ...
        'pc = VALUES(pc);'];
    mysql(query);
end%for i
mysql('close');

%% Make ASCII, CDF and CSV files
pcnq2iaga(t1, t2, marg, data);

end%function