%% Calculation of definitive PC index
% Compute the definitive PC index for the northern hemisphere (thl) or the
% southern hemmisphere (vos).
% 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             |
% +-----------+--------------+------+-------------------+-------------------------+
% | ss_x      | float        | Yes  | NULL              | SS(x)                   |
% +-----------+--------------+------+-------------------+-------------------------+
% | ss_y      | float        | Yes  | NULL              | SS(y)                   |
% +-----------+--------------+------+-------------------+-------------------------+
% | qdc_x     | float        | Yes  | NULL              | QDC(x)                  |
% +-----------+--------------+------+-------------------+-------------------------+
% | qdc_y     | float        | Yes  | NULL              | QDC(y)                  |
% +-----------+--------------+------+-------------------+-------------------------+
% | pc        | float        | Yes  | NULL              | PC index                |
% +-----------+--------------+------+-------------------+-------------------------+
%
% NOTE! To compute the definitive PC index a data overlap of +/- 1 month is
% needed
%
% Change list:
% - Entire code rewritten to be more time efficient
% - x, y, ss_x, ss_y, qdc_x and qdc_y changed from mediumint(9) to float
% - Data overlap added to avoid edge effects (+/- 1 month, can be modified)
% - The actual QDC is now calculate for a 30-day period instead of a 31-day
%   period
% - 'moving' added to smooth statement for "absolute deviations from trend"
%   computation in q_day subfunction
% - "Actual Day" is now saved to correct date - instead of the day after
% - All data are now written to the SQL database at the same time in
%   batches of 1000 lines (max limit of SQL INSERT-statement) to be more
%   time efficient
%
% $Revision: 2.0 $

% $Id: Definitive_PC_Index.m,v 2.0 2018/12/12 10:42:36 jnsnl Exp $

%% Input
% Select station [Qaanaaq (thl) or Vostok (vos)]
station = 'thl';

app = true; % Append available data on existing SQL table
if app
  % Time period (automatic append available data on existing SQL table) 
  [t1, t2] = append_db(station);
else
  % Time periode (manual)
  t1 = 1975;      % start year
  t2 = 2014;      % end year
end%if app
% Data overlap +/- n days
overlap = 31; % 31 = +/- 1 month

% Path to definitive IAGA2002 files
iaga_dir = ['/nfs/r14/magobs/PC-index/Definitive_PC/data/' ...
            station '_archive/'];

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

% Longitude of station
switch station
    case 'thl'
        lon = 291;
    case 'vos'
        lon = 106.9;
end%switch

% 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);

%% Calculate definitive PC index for [station] in time period [t1] to [t2]
for year = t1:t2
    disp(['Computing definitive PC index (' station ') for ' num2str(year)])
    %% Initialize data structure
    data  = struct('t', [], 'x', [], 'y', [], 'ss_x', [], 'ss_y', [], ...
                   'qdc_x', [], 'qdc_y', [], 'pc', []);
    % Number of days (year +/- data overlap)
    ndays = sum(eomday(year,1:12)) + overlap*2;
    % Write time in minutes (year +/- data overlap) to data structure
    data.t = (datenum(year, 1, 1, 0, 0:ndays*1440-1, 0) - overlap)';
    % Fill the remaining data structure with NaN's
    fns = fieldnames(data);
    for fn = fns(2:end)' % [x, y, ss_x, ss_y, qdc_x, qdc_y, pc]
        data.(char(fn)) = NaN(length(data.t), 1);
    end%for fn
    
    %% Read definitive IAGA2002 data and add to data structure
    % Read every file in the year (+/- margin)
    for actual = datenum(year, 1, 1)-overlap:datenum(year, 12, 31)+overlap;
        % /[yyyy]/thl[yyyymmdd]dmin.min
        iaga_name = fullfile(iaga_dir, datestr(actual,'yyyy'), ...
                    [station, datestr(actual,'yyyymmdd'), 'dmin.min']);
        
        if ~exist(iaga_name, 'file')
          warning(['No ' iaga_name(end-18:end) '-file. Continue.'])
          continue
        end%if ~exist
        fid = fopen(iaga_name);
        % Find end of header (Not uniform across IAGA2002 files)
        % DATE       TIME         DOY     THLX      THLY      THLZ      THLF   
        % 2011-12-31 23:59:00.000 365      2510.90  -3229.90  56246.10  88888.80
        while ~feof(fid)
            f = fgetl(fid);
            if strcmp(f(1:4), 'DATE')
              % Last header-line
              break
            end%if strcmp
        end%while ~feof
        % Read data form IAGA2002 file
        f = textscan(fid, '%s %s %d %f %f %f %f');
        fclose(fid);
        % Find index into data structure
        [~, i] = intersect(data.t, datenum([char(f{1}) char(f{2})], ...
                                              'yyyy-mm-ddHH:MM:SS.FFF'));
        % Insert x and y into data structure
        for j = 2:3
            f{j+2}((f{j+2} >= 88888)) = NaN;
            data.(char(fns(j)))(i) = f{j+2};
        end%for j
    end%for actual
    
    %% Calculate SS (sector structure)
    % Calculate ss_x and ss_y
    for fn = fns(2:3)' % [x, y]
        % Calculate daily median value
        day_med = nanmedian(reshape(data.(char(fn)), [1440 ndays]));
        % Smooth medians
        w = smooth(day_med, 7, 'moving'); % 7 days
        w = smooth(w, 7, 'rloess');
        
        % Interpolation to minutes and add to data structure
        data.(['ss_' char(fn)]) = interp1(1:ndays, w, ...
                                          1:1/1440:ndays+1-1/1440, ...
                                          'spline', 'extrap')';
    end%for fn
    
    %% Calculate QDC (Actual Day)
    % Shift 30-day period by 10 days at a time to calculation a series of
    % actual QD:
    for fn = fns(2:3)' %[x, y]
        for start_dt = datenum(year, 1, 1)-overlap:10:datenum(year, 12, 31)
            % Calculate QDC for 30-days dataset and save it into
            % "Actual Day" in data structure
            end_dt = start_dt+30-1/1440;
            i = find(data.t >= start_dt & data.t <= end_dt);
            
            % Calculate COMP-SS
            by = data.(char(fn))(i) - data.(['ss_' char(fn)])(i);
            
            % Calculate QDC for an "Actual Day"
            [qday, ActDay] = q_day(by);
            
            % 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;
              data.(['qdc_' char(fn)])(j) = qday;
            end%if ~isnan
        end%for start_dt
        
        % Interpolation of QDC for every day from "Actual Day" array
        % "Actual Days"
        ActDays = find(~isnan(sum(reshape(data.(['qdc_' char(fn)]), ...
                                          [1440 ndays]))));
        
        % Construct 2D array of Actual QDC
        qdays = reshape(data.(['qdc_' char(fn)])(~isnan(data.(['qdc_' char(fn)]))), ...
                        [1440 length(ActDays)]);
        
        % Interpolation and smoothing QDC for every day
        interp_arr = NaN(1440, ndays);
        for i = 1:1440
            interp_arr(i,:) = smooth(interp1(ActDays, qdays(i,:), 1:ndays, ...
                                     'nearest', 'extrap'), 60, 'lowess');
        end
        
        % Make 1D array of the QDC, smooth and add to data structure
        l_arr = reshape(interp_arr, [1 ndays*1440]);
        data.(['qdc_' char(fn)]) = smooth(l_arr,120,'loess');
    end%for fn
    
    %% Caculate PC index
    % Calculate disturbed part of the variations
    dist_x = data.x - data.ss_x - data.qdc_x;
    dist_y = data.y - data.ss_y - data.qdc_y;
    
    % H Projection and PC calculation for the year (without data overlap)
    i = 1+overlap*1440:(ndays-overlap)*1440;
    [~, ~, ~, hour, minute, ~] = datevec(data.t(i));
    % Minute in the year
    min_in_year = i-overlap*1440;
    % Calculate H_proj
    UT = hour .* 15 + minute .* 0.25; % TIME to DEG
    y = lon + coeff.f(min_in_year) + UT;
    H_proj = dist_x(i).*sin(y.*pi./180)-dist_y(i).*cos(y.*pi./180);
    % Calculate PC index and add to data structure
    data.pc(i) = (H_proj - coeff.b(min_in_year))./coeff.a(min_in_year);
    
    %% Write data structure to SQL database
    mysql('open', serv, log, pass);
    mysql(['use ' dbname]);
    % Write data to SQL database in bulks of 1000 rows (max limit of SQL 
    % INSERT-statement) to save time
    bulk = [overlap*1440:1000:(ndays-overlap)*1440, (ndays-overlap)*1440];
    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, ss_x, ss_y, qdc_x, qdc_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, ss_x, ss_y, qdc_x, qdc_y, pc) ' ...
                 'VALUES ' str(1:end-2) ' '...
                 'ON DUPLICATE KEY UPDATE x = VALUES(x), y = VALUES(y), ' ...
                        'ss_x = VALUES(ss_x), ss_y = VALUES(ss_y), ' ...
                        'qdc_x = VALUES(qdc_x), qdc_y = VALUES(qdc_y), ' ...
                        'pc = VALUES(pc);'];
        mysql(query);
    end%for i
    mysql('close');
end%for year