%% EEJ (binned) QD vs LT
clear all

load EEJ_A_extended.mat
tmp = load('EEJ_B_extended');
t_all = [t; tmp.t]; phi_all = [phi; tmp.phi]; J_all = [J; tmp.J];


idx = t > md2000(2023, 2, 11);
t(idx) = []; phi(idx) = []; J(idx,:) = [];

T = mod(t+phi/360,1)*24;
phi = mod(phi, 360);


dT = .5; T_vector = [dT/2:dT:24]';
EEJ_mean = NaN*ones(length(T_vector), 81);
for k = 1:length(T_vector)
    for k_QD = 1:81
        idx = abs(T-T_vector(k)) < dT;
        EEJ_mean(k,k_QD) = nanmedian(J(idx,k_QD));
    end % n
end % k

clf
QD_lat = [-20:.5:20';]
imagesc(T_vector, QD_lat, EEJ_mean')
hold on
contour(T_vector, QD_lat, EEJ_mean', [-15:5:150], '-w')
hold off
axis([0 24 -20 20])
set(gca, 'XTick', [0:3:24], 'YTick', [-20:4:20])
ylabel 'QD-latitude [deg]';
xlabel 'LT [hours]';
clim(100*[-1 1])
hc = colorbar;
hc.Label.String = '[A/km]';
% hc.Layout.Tile = 'east';
cmap = brewermap(128, 'RdBu');
colormap(flipud(cmap));

set(gcf, 'PaperOrientation', 'Portrait', 'PaperType', 'A4', 'PaperUnits', ...
    'centimeters', 'PaperPosition', [0.2 0.2 [25 12]]);
filename_out = 'EEJ_binned_LT_QD.png';
print('-dpng', '-r200', filename_out); crop_png(filename_out);


%% extract EEJ at dip-equator
if size(J,2) > 1; J = J(:,41); end

% remove outliers
q_1 = prctile(J,1);
q_99 = prctile(J,99);
dq = q_99 - q_1;
idx = q_1-1.5*dq < J & J < q_99 + 1.5*dq & ~isnan(J);
t = t(idx); phi = phi(idx); J = J(idx); T = T(idx);


%% T vs season and T vs nu at QD = 0
dT = .5; T_vector = [dT/2:dT:24]';
dnu = 1; nu_vector = [dnu/2:dnu:24]';
ds = 1; s_vector = [ds/2:ds:24]';
dphi = 5; phi_vector = [dphi/2:dphi:360]';

s = mod(t/365.25, 1)*24;
nu = moon_phase(t);

hl = tiledlayout(1,3, 'Padding','compact', 'TileSpacing','tight');
EEJ_mean = NaN*ones(length(T_vector), length(phi_vector));
for k = 1:length(T_vector)
    for n = 1:length(phi_vector)
        idx = abs(T-T_vector(k)) < dT & abs(phi-phi_vector(n)) < dphi;
        EEJ_mean(k,n) = nanmedian(J(idx));
    end % n
end % k
EEJ_mean = imboxfilt(EEJ_mean, 3);

nexttile
imagesc(T_vector, phi_vector, EEJ_mean')
hold on
contour(T_vector, phi_vector, EEJ_mean', [-150:10:150], '-w')
hold off
axis([0 24 0 360])
set(gca, 'XTick', [0:3:24], 'YTick', [0:60:360])
ylabel 'longitude \phi [deg]';
xlabel 'LT [hours]';
clim(100*[-1 1])

nexttile
EEJ_mean = NaN*ones(length(T_vector), length(s_vector));
for k = 1:length(T_vector)
    for n = 1:length(s_vector)
        idx = abs(T-T_vector(k)) < dT & abs(s-s_vector(n)) < dnu;
        EEJ_mean(k,n) = nanmedian(J(idx));
    end % n
end % k
EEJ_mean = imboxfilt(EEJ_mean, 3);
imagesc(T_vector, s_vector/2, EEJ_mean')
hold on
contour(T_vector, s_vector/2, EEJ_mean', [-150:10:150], '-w')
hold off
axis([0 24 0 12])
set(gca, 'XTick', [0:3:24], 'YTick', [0:3:12])
ylabel 'season s [months]';
xlabel 'LT [hours]';
clim(100*[-1 1])

nexttile
EEJ_mean = NaN*ones(length(T_vector), length(nu_vector));
for k = 1:length(T_vector)
    for n = 1:length(nu_vector)
        idx = abs(T-T_vector(k)) < dT & abs(nu-nu_vector(n)) < dnu;
        EEJ_mean(k,n) = nanmedian(J(idx));
    end % n
end % k
H = ones(2,2);
EEJ_mean = imboxfilt(EEJ_mean, 3);
imagesc(T_vector, nu_vector, EEJ_mean')
hold on
contour(T_vector, nu_vector, EEJ_mean', [-150:10:150], '-w')
hold off
axis([0 24 0 24])
set(gca, 'XTick', [0:3:24], 'YTick', [0:3:24])
ylabel 'Lunar phase \nu [hrs]';
xlabel 'LT [hours]';
clim(100*[-1 1])

hc = colorbar;
hc.Label.String = '[A/km]';
hc.Layout.Tile = 'south';
cmap = brewermap(128, 'RdBu');
colormap(flipud(cmap));

set(gcf, 'PaperOrientation', 'Portrait', 'PaperType', 'A4', 'PaperUnits', ...
    'centimeters', 'PaperPosition', [0.2 0.2 [25 12]]);
filename_out = 'EEJ_binned_LT_phi_season_nu.png';
print('-dpng', '-r200', filename_out); crop_png(filename_out);

%% maps of model values
clear all
rad = pi/180;

% read the model parameters
[ab, R, QD_lat, params] = read_EEJM('EEJ_SW-AB_01.txt');


%% maps at dip-equator for different seasons
F107 = 100;

k_QD = 41; % dip equator
m_i = ab(:,k_QD);

T = [0:1/100:1]';
phi = [0:5:360]';

[T_grid, phi_grid] = ndgrid(T, phi);

season = char('March','June','September','December');

hl = tiledlayout(2,2, 'Padding','compact', 'TileSpacing','tight');
for i_season = 1:4
    nexttile
    t_grid = T_grid-phi_grid/360 + round(i_season*365.25/4);
    nu_grid = zeros(size(t_grid));

    s_0 = md2000(2000, 3, 22);
    J = zeros(size(t_grid));
    n = 1;
    for p = 0:params.n_t
        for m = -params.n_phi:params.n_phi
            for k = -params.n_s:params.n_s
                if params.n_l == 1
                    l_sel = 0;
                elseif params.n_l == 2
                    l_sel = [-2 0];
                end
                for l = l_sel
                    % coefficients a
                    if ~( (p == 0 && l == 0 && m < 0) || (p == 0 && l == 0 && m == 0 && k < 0) )
                        J = J + m_i(n)*cos(p*t_grid*360*rad + m*phi_grid*rad + k*(t_grid-s_0)*360*rad/365.25 + l*nu_grid*15*rad);
                        n = n + 1;
                    end
                    % coefficients b
                    if ~( (m < 0 && p == 0 && l == 0) || (p == 0 && m == 0 && l == 0 && k < 0) || (p == 0 && m == 0 && l == 0 && k == 0) )
                        J = J + m_i(n)*sin(p*t_grid*360*rad + m*phi_grid*rad + k*(t_grid-s_0)*360*rad/365.25 + l*nu_grid*15*rad);
                        n = n + 1;
                    end
                end % l
            end % k
        end % m
    end % p

    J = (1+R(k_QD)*F107)*J;

    imagesc(T*24, phi, J')
    hold on
    contour(T*24, phi, J', [-150:10:150], '-w')
    hold off
    set(gca, 'YTick', [0:60:360], 'XTick', [0:3:24])
    if mod(i_season, 2) == 1; ylabel 'longitude \phi [deg]'; end
    if i_season > 2; xlabel 'LT [hours]'; end
    clim(100*[-1 1])
    title(season(i_season,:))
end
hc = colorbar;
hc.Label.String = '[A/km]';
hc.Layout.Tile = 'south';
% blue_red_colormap
cmap = brewermap(128, 'RdBu');
colormap(flipud(cmap));

set(gcf, 'PaperOrientation', 'Portrait', 'PaperType', 'A4', 'PaperUnits', ...
    'centimeters', 'PaperPosition', [0.2 0.2 [25 15]]);
filename_out = 'EEJ_phi_LT_season.png';
print('-dpng', '-r200', filename_out); crop_png(filename_out);

%% maps at dip-equator for different QD-latitudes
F107 = 100;

hl = tiledlayout(3,3, 'Padding','compact', 'TileSpacing','tight');
i_season = 1;
% for k_QD = 33:2:49
i_plt = 1;
for k_QD = 29:3:49+6
    m_i = ab(:,k_QD);

    %
    T = [0:1/100:1]';
    phi = [0:5:360]';

    [T_grid, phi_grid] = ndgrid(T, phi);

    nexttile
    t_grid = T_grid-phi_grid/360 + round(i_season*365.25/4);
    nu_grid = zeros(size(t_grid));

    s_0 = md2000(2000, 3, 22);
    J = zeros(size(t_grid));
    n = 1;
    for p = 0:params.n_t
        for m = -params.n_phi:params.n_phi
            for k = -params.n_s:params.n_s
                if params.n_l == 1
                    l_sel = 0;
                elseif params.n_l == 2
                    l_sel = [-2 0];
                end
                for l = l_sel
                    % coefficients a
                    if ~( (p == 0 && l == 0 && m < 0) || (p == 0 && l == 0 && m == 0 && k < 0) )
                        J = J + m_i(n)*cos(p*t_grid*360*rad + m*phi_grid*rad + k*(t_grid-s_0)*360*rad/365.25 + l*nu_grid*15*rad);
                        n = n + 1;
                    end
                    % coefficients b
                    if ~( (m < 0 && p == 0 && l == 0) || (p == 0 && m == 0 && l == 0 && k < 0) || (p == 0 && m == 0 && l == 0 && k == 0) )
                        J = J + m_i(n)*sin(p*t_grid*360*rad + m*phi_grid*rad + k*(t_grid-s_0)*360*rad/365.25 + l*nu_grid*15*rad);
                        n = n + 1;
                    end
                end % l
            end % k
        end % m
    end % p

    J = (1+R(k_QD)*F107)*J;

    imagesc(T*24, phi, J')
    hold on
    contour(T*24, phi, J', [-150:10:150], '-w')
    hold off
    set(gca, 'YTick', [0:60:360], 'XTick', [0:3:24])
    if mod(i_plt, 3) == 1; ylabel 'longitude \phi [deg]'; end
    if i_plt > 6; xlabel 'LT [hours]'; end
    clim(100*[-1 1])
    title(['QD-latitude: ' num2str(QD_lat(k_QD), '%234.1f') '^\circ'])
    i_plt = i_plt + 1;
end
hc = colorbar;
hc.Label.String = '[A/km]';
hc.Layout.Tile = 'south';
cmap = brewermap(128, 'RdBu');
colormap(flipud(cmap));

set(gcf, 'PaperOrientation', 'Portrait', 'PaperType', 'A4', 'PaperUnits', ...
    'centimeters', 'PaperPosition', [0.2 0.2 [25 15]]);
filename_out = 'EEJ_phi_LT_QD_lat.png';
print('-dpng', '-r200', filename_out); crop_png(filename_out);

%% rms statistics vs LT
clear all
rad = pi/180;

[ab, R, QD_lat, params] = read_EEJM('EEJ_SW-AB_01.txt');
% [ab, ab_EUVAC, QD_lat, params] = read_EEJM_EUVAC('EEJ_SW-AB_EUVAC_01.txt');

k_QD = 41;
m_i = ab(:,k_QD);
% m_i_EUVAC = ab_EUVAC(:,k_QD);

% define some colors
co = [        0    0.4470    0.7410
    0.8500    0.3250    0.0980
    0.9290    0.6940    0.1250
    0.4940    0.1840    0.5560
    0.4660    0.6740    0.1880
    0.3010    0.7450    0.9330
    0.6350    0.0780    0.1840];

% read daily values of F10.7
filename = 'SW_OPER_AUX_F10_2__20060101T000000_20231015T000000_0001.DBL';
[t_F107_all, F107_all]  = textread(filename, '%f %f', 'commentstyle', 'shell');

for i_sat = 1:5
    if i_sat == 1
        load EEJ_A_extended.mat
    elseif i_sat == 2
        load('EEJ_B_extended.mat');
    elseif i_sat == 3
        load EEJ_C.mat
    elseif i_sat == 4
        load EEJ_CSES.mat
    elseif i_sat == 5
        load EEJ_C.mat
        tmp = load('EEJ_CSES.mat');
        J = interp1(t, J, tmp.t, 'nearest');
        phi = interp1(t, phi, tmp.t, 'nearest');
        t = interp1(t, t, tmp.t, 'nearest');
    elseif i_sat == 6
        load EEJ_A_extended.mat
        idx = find(EUVAC > 81);
        t(idx) = []; phi(idx) = []; J(idx) = []; nu(idx) = []; EUVAC(idx) = [];
    elseif i_sat == 7
        load EEJ_A_extended.mat
        idx = find(EUVAC < 81);
        t(idx) = []; phi(idx) = []; J(idx) = []; nu(idx) = []; EUVAC(idx) = [];
    end

    if size(J,2) > 1; J = J(:,41); end
    idx = t > md2000(2023, 2, 11);
    t(idx) = []; phi(idx) = []; J(idx) = [];

    % remove outliers
    q_1 = prctile(J,1);
    q_99 = prctile(J,99);
    dq = q_99 - q_1;
    idx = q_1-1.5*dq < J & J < q_99 + 1.5*dq & ~isnan(J);
    t = t(idx); phi = phi(idx); J = J(idx); % nu = nu(idx); EUVAC = EUVAC(idx);
    disp([(length(idx)-sum(idx))/length(idx)*100])

    idx = ~isnan(F107_all);
    F107 = interp1(t_F107_all(idx), zp_filter(F107_all(idx), 'boxcar', 81), t, 'linear'); % interpolate F107 to time instants of "synthetic" data
    EUVAC = (F107 + interp1(t_F107_all(idx), F107_all(idx), t, 'linear'))/2;
    nu = moon_phase(t);
    T = mod(t+phi/360, 1)*24;

    %
    J_mod = zeros(size(t));
    n = 1;
    for p = 0:params.n_t
        for m = -params.n_phi:params.n_phi
            for k = -params.n_s:params.n_s
                if params.n_l == 1
                    l_sel = 0;
                elseif params.n_l == 2
                    l_sel = [-2 0];
                end
                for l = l_sel
                    % coefficients a
                    if ~( (p == 0 && l == 0 && m < 0) || (p == 0 && l == 0 && m == 0 && k < 0) )
                        J_mod = J_mod + m_i(n)*cosd(p*t*360 + m*phi + k*t*360/365.25 + l*nu*15);
                        n = n + 1;
                    end
                    % coefficients b
                    if ~( (m < 0 && p == 0 && l == 0) || (p == 0 && m == 0 && l == 0 && k < 0) || (p == 0 && m == 0 && l == 0 && k == 0) )
                        J_mod = J_mod + m_i(n)*sind(p*t*360 + m*phi + k*t*360/365.25 + l*nu*15);
                        n = n + 1;
                    end
                end % l
            end % k
        end % m
    end % p

    J_mod = (1+R(k_QD)*EUVAC).*J_mod;

    % statistics
    delta_d = [J - J_mod];

    w = ones(size(delta_d));
    for LT = 1:24
        idx = abs(LT - T - .5) < .5;
        for k = 1:50
            % rms = sqrt(((w'.*delta_d') * (delta_d)) / sum(w));
            % only dayside data to calæculate weights
            % idx = abs(T-12) < 6;
            rms = sqrt(((w(idx)'.*delta_d(idx)') * (delta_d(idx))) / sum(w(idx)));
            w(idx) = min(1.5./abs(delta_d(idx)./rms),1);  % robust Huber weights
            % c = 6.5; residuals = delta_d./rms;
            % w = (1-(residuals/c).^2).^2; w(find(abs(residuals) >= c)) = 0; %Tuckey weights
        end
    end % LT
    % sqrt(delta_d'*(w.*delta_d)/length(t))

    rms_LT = zeros(24,1);
    ave_LT = zeros(24,1);
    w_LT = zeros(24,1);
    for LT = 1:24
        idx = abs(LT - T - .5) < .5;
        rms_LT(LT) = sqrt(delta_d(idx)'*(w(idx).*delta_d(idx))/sum(w(idx)));
        rms_LT(LT) = mad(delta_d(idx));
        ave_LT(LT) = sum(w(idx).*delta_d(idx))/sum(w(idx));
        w_LT(LT) = sum(w(idx))/sum(idx);
    end

    if i_sat == 1
        rms_LT_A = rms_LT;
        ave_LT_A = ave_LT;
        w_LT_A = w_LT;
    elseif i_sat == 2
        rms_LT_B = rms_LT;
        ave_LT_B = ave_LT;
        w_LT_B = w_LT;
    elseif i_sat == 3
        rms_LT_C = rms_LT;
        ave_LT_C = ave_LT;
        w_LT_C = w_LT;
    elseif i_sat == 4
        rms_LT_CSES = rms_LT;
        ave_LT_CSES = ave_LT;
        w_LT_CSES = w_LT;
    elseif i_sat == 5
        rms_LT_C_CSES = rms_LT;
        ave_LT_C_CSES = ave_LT;
        w_LT_C_CSES = w_LT;
    elseif i_sat == 6
        rms_LT_A_F1 = rms_LT;
    elseif i_sat == 7
        rms_LT_A_F2 = rms_LT;
    end
end

%
clf
h1 = plot([.5:24], rms_LT_A, '.', 'MarkerSize', 30)
hold on
h2 = plot([.5:24], rms_LT_B, '.', 'MarkerSize', 30)
h3 = plot([.5:24], rms_LT_C, '.', 'MarkerSize', 30)
h4 = plot(13.5, nanmean(rms_LT_CSES), '.', 'MarkerSize', 30, 'Color', co(4,:))
h5 = plot(13.7, nanmean(rms_LT_C_CSES(13:15)), '.', 'MarkerSize', 30, 'Color', co(5,:))
% plot([.5:24], rms_LT_A_F1, '--', 'LineWidth', 1, 'Color', co(1,:))
% plot([.5:24], rms_LT_A_F2, '-', 'LineWidth', 1, 'Color', co(1,:))
plot([.5:24], rms_LT_A, 'LineWidth', 2, 'Color', co(1,:))
plot([.5:24], rms_LT_B, 'LineWidth', 2, 'Color', co(2,:))
plot([.5:24], rms_LT_C, 'LineWidth', 2, 'Color', co(3,:))
% plot([.5:24], ave_LT_A, 'LineWidth', 1, 'Color', co(1,:))
% plot([.5:24], ave_LT_B, 'LineWidth', 1, 'Color', co(2,:))
% plot([.5:24], ave_LT_C, 'LineWidth', 1, 'Color', co(3,:))
hold off
axis([0 24 0 32])
grid on
set(gca, 'XTick', [0:3:24])
xlabel 'LT [hours]'
ylabel 'MAD [A/km]'
legend('Swarm A (EEF extended)', 'Swarm B (EEF extended)', 'Swarm C (EEF operational)', 'CSES', ...
    'Swarm C (CSES timespan)', 'Location','south')

set(gcf, 'PaperOrientation', 'Portrait', 'PaperType', 'A4', 'PaperUnits', ...
    'centimeters', 'PaperPosition', [0.2 0.2 [25 12]]);
filename_out = 'rms_EEJ_LT.eps';
% print_pdf(filename_out)

%%
tmp = load('EEJ_A_extended_EEJM2.txt');
idx = tmp(:,5) < -9e3;
tmp(idx,:) = [];
EEJM2.t = tmp(:,1);
EEJM2.phi = tmp(:,2);
EEJM2.EUVAC = tmp(:,3);
EEJM2.nu = tmp(:,4);
EEJM2.J = tmp(:,5)*1e3;
EEJM2.T = mod(EEJM2.t+EEJM2.phi/360,1)*24;

load EEJ_A_extended.mat
EEJM2.J_obs = interp1(t, J(:,41), EEJM2.t);

% remove outliers
q_1 = prctile(EEJM2.J_obs,1);
q_99 = prctile(EEJM2.J_obs,99);
dq = q_99 - q_1;
idx = q_1-1.5*dq < EEJM2.J_obs & EEJM2.J_obs < q_99 + 1.5*dq & ~isnan(EEJM2.J_obs);
EEJM2.t = t(idx); EEJM2.phi = phi(idx); EEJM2.T = EEJM2.T(idx);
EEJM2.J = EEJM2.J(idx); EEJM2.J_obs = EEJM2.J_obs(idx); 

%
rms_LT = NaN*ones(24,1);
delta_d = EEJM2.J_obs - EEJM2.J;
for LT = 1:24
    idx = abs(LT - EEJM2.T - .5) < .5;
    if sum(idx) > 10
    rms_LT(LT) = mad(delta_d(idx));
    end
end

hold on
h6 = plot([.5:24], rms_LT, '.', 'Markersize', 30, 'Color', co(6,:))
plot([.5:24], rms_LT, 'LineWidth', 2, 'Color', co(6,:))
hold off

legend([h1 h2 h3 h4 h5 h6], 'Swarm A (EEF extended)', 'Swarm B (EEF extended)', 'Swarm C (EEF operational)', 'CSES', ...
    'Swarm C (CSES timespan)', 'EEJM-2.0 for Swarm A', 'Location','south')


set(gcf, 'PaperOrientation', 'Portrait', 'PaperType', 'A4', 'PaperUnits', ...
    'centimeters', 'PaperPosition', [0.2 0.2 [25 12]]);
filename_out = 'rms_EEJ_LT.eps';
print_pdf(filename_out)


%% maps of observed (Swarm A operational, i.e. only daytime) and modeled EEJ profiles
clear all

[ab, R, QD_lat, params] = read_EEJM('EEJ_SW-AB_01.txt');

% read daily values of F10.7
filename = 'SW_OPER_AUX_F10_2__20060101T000000_20231015T000000_0001.DBL';
[t_F107_all, F107_all]  = textread(filename, '%f %f', 'commentstyle', 'shell');

load('EEJ_A_1.mat'); % read Swarm A daytime EEJ values (from EEF_OPER)

% first 6 months of 2018
t_start = datetime(2018, 1, 1);
t_end =  datetime(2018, 7, 1);
filename_out = 'J_observed_modelled_2018.png';

% first 6 months of 2015
% t_start = datetime(2015, 1, 1);
% t_end =  datetime(2015, 7, 1);
% filename_out = 'J_observed_modelled_2015.png';


idx = ~isnan(F107_all);
F107 = interp1(t_F107_all(idx), zp_filter(F107_all(idx), 'boxcar', 81), t, 'linear'); % interpolate F107 to time instants of "synthetic" data
EUVAC = (F107 + interp1(t_F107_all(idx), F107_all(idx), t, 'linear'))/2;
nu = moon_phase(t);
t_dt = datetime(t+datenum(2000,1,1),'ConvertFrom','datenum');
T = mod(t+phi/360,1)*24;

QD_lat = [-20:.5:20]';

% synthetic values
s_0 = md2000(2000, 3, 22);
J_mod = zeros(size(J));
n = 1;
for p = 0:params.n_t
    for m = -params.n_phi:params.n_phi
        for k = -params.n_s:params.n_s
            if params.n_l == 1
                l_sel = 0;
            elseif params.n_l == 2
                l_sel = [-2 0];
            end
            for l = l_sel
                % coefficients a
                if ~( (p == 0 && l == 0 && m < 0) || (p == 0 && l == 0 && m == 0 && k < 0) )
                    tmp = cosd(p*t*360 + m*phi + k*t*360/365.25 + l*nu*15);
                    for k_QD = 1:81
                        J_mod(:,k_QD) = J_mod(:,k_QD) + ab(n,k_QD)*tmp;
                    end
                    n = n + 1;
                end
                % coefficients b
                if ~( (m < 0 && p == 0 && l == 0) || (p == 0 && m == 0 && l == 0 && k < 0) || (p == 0 && m == 0 && l == 0 && k == 0) )
                    tmp = sind(p*t*360 + m*phi + k*t*360/365.25 + l*nu*15);
                    for k_QD = 1:81
                        J_mod(:,k_QD) = J_mod(:,k_QD) + ab(n,k_QD)*tmp;
                    end
                    n = n + 1;
                end
            end % l
        end % k
    end % m
end % p

J_mod = (1+R(k_QD)*EUVAC).*J_mod;

%
idx_QD = 21:61;


idx = find(t_start < t_dt & t_dt < t_end);

hl = tiledlayout(1,7, 'Padding','compact', 'TileSpacing','tight');
nexttile
plot(T(idx), t_dt(idx), '.');
set(gca, 'XLim', [0 24], 'XTick', [0:6:24])
grid on
xlabel 'LT [hrs]'

nexttile([1 3])
imagesc(QD_lat(idx_QD), t_dt(idx), J(idx,idx_QD))
caxis(100*[-1 1])
set(gca, 'YDir', 'normal', 'YLim', [t_start t_end])
xlabel 'QD-latitude [deg]'
disp(['Mean EUVAC: ' num2str(mean(EUVAC(idx)), '%3.0f s.f.u.')])
title 'Observed J'

nexttile([1 3])
imagesc(QD_lat(idx_QD), t_dt(idx), J_mod(idx,idx_QD))
hc = colorbar;
hc.Label.String = '[A/km]';
hc.Layout.Tile = 'east';
cmap = brewermap(128, 'RdBu');
colormap(flipud(cmap));
caxis(100*[-1 1])
set(gca, 'YDir', 'normal', 'YLim', [t_start t_end])
xlabel 'QD-latitude [deg]'
title 'Model J'

set(gcf, 'PaperOrientation', 'Portrait', 'PaperType', 'A4', 'PaperUnits', ...
    'centimeters', 'PaperPosition', [0.2 0.2 [25 10]]);
print('-dpng', '-r200', filename_out); crop_png(filename_out);


%% maps LT vs nu
F107 = 100;

k_QD = 41;
m_i = ab(:,k_QD);

T = [0:1/100:1]';
nu = [0:1/100:1]*24';

[T_grid, nu_grid] = ndgrid(T, nu);


i_season = 1;
phi_grid = zeros(size(T_grid));
t_grid = T_grid-phi_grid/360 + round(i_season*365.25/4);
t_grid = -phi_grid/360 + round(i_season*365.25/4);
phi_grid = zeros(size(t_grid));

s_0 = md2000(2000, 3, 22);
J = zeros(size(t_grid));
n = 1;
for p = 0:params.n_t
    for m = -params.n_phi:params.n_phi
        for k = -params.n_s:params.n_s
            if params.n_l == 1
                l_sel = 0;
            elseif params.n_l == 2
                l_sel = [-2 0];
            end
            for l = l_sel
                % coefficients a
                if ~( (p == 0 && l == 0 && m < 0) || (p == 0 && l == 0 && m == 0 && k < 0) )
                    J = J + m_i(n)*cos(p*t_grid*360*rad + m*phi_grid*rad + k*(t_grid-s_0)*360*rad/365.25 + l*nu_grid*15*rad);
                    n = n + 1;
                end
                % coefficients b
                if ~( (m < 0 && p == 0 && l == 0) || (p == 0 && m == 0 && l == 0 && k < 0) || (p == 0 && m == 0 && l == 0 && k == 0) )
                    J = J + m_i(n)*sin(p*t_grid*360*rad + m*phi_grid*rad + k*(t_grid-s_0)*360*rad/365.25 + l*nu_grid*15*rad);
                    n = n + 1;
                end
            end % l
        end % k
    end % m
end % p

J = (1+R(k_QD)*F107)*J;

imagesc(T*24, nu, J')
hold on
contour(T*24, nu, J', [-150:10:150], '-w')
hold off
set(gca, 'YTick', [0:60:360], 'XTick', [0:3:24])
if mod(i_season, 2) == 1; ylabel 'longitude \phi [deg]'; end
if i_season > 2; xlabel 'LT [hours]'; end
clim(100*[-1 1])
hc = colorbar;
hc.Label.String = '[A/km]';
% hc.Layout.Tile = 'east';
cmap = brewermap(128, 'RdBu');
colormap(flipud(cmap));

set(gcf, 'PaperOrientation', 'Portrait', 'PaperType', 'A4', 'PaperUnits', ...
    'centimeters', 'PaperPosition', [0.2 0.2 [25 15]]);
filename_out = 'EEJ_phi_LT_season.png';
% print('-dpng', '-r200', filename_out); crop_png(filename_out);

%% Correlation between Swarm A and B
clear all

A = load('EEJ_A_extended.mat');
B = load('EEJ_B_extended.mat');

k_QD = 41;
t = A.t;
phi = A.phi;
nu = A.nu;
EUVAC = A.EUVAC;
J_A = A.J(:,k_QD);
J_B = interp1(B.t, B.J(:,k_QD), t, 'nearest');
T_A = mod(t+phi/360,1)*24;
T_B = interp1(B.t, mod(B.t+B.phi/360,1)*24, t, 'nearest');

% remove outliers SW-A
q_1 = prctile(J_A,1);
q_99 = prctile(J_A,99);
dq = q_99 - q_1;
idx = q_1-1.5*dq < J_A & J_A < q_99 + 1.5*dq & ~isnan(J_A);
t = t(idx); phi = phi(idx); J_A = J_A(idx); J_B = J_B(idx); T_A = T_A(idx); T_B = T_B(idx); nu = nu(idx); EUVAC = EUVAC(idx);

% remove outliers SW-B
q_1 = prctile(J_B,1);
q_99 = prctile(J_B,99);
dq = q_99 - q_1;
idx = q_1-1.5*dq < J_B & J_B < q_99 + 1.5*dq & ~isnan(J_B);
t = t(idx); phi = phi(idx); J_A = J_A(idx); J_B = J_B(idx); T_A = T_A(idx); T_B = T_B(idx); nu = nu(idx); EUVAC = EUVAC(idx);

LT = [0:24]';
rho = zeros(length(LT));
for i_LT_A = 1:length(LT)
    for i_LT_B = 1:length(LT)
        idx = abs(T_A - LT(i_LT_A)) < .5 & abs(T_B - LT(i_LT_B)) < .5;
        rho(i_LT_A, i_LT_B) = corr(J_A(idx), J_B(idx));
    end
end

imagesc(LT, LT, rho)
caxis([-1 1])
colorbar
xlabel 'LT of SW-A [hours]'
ylabel 'LT of SW-B [hours]'
set(gca, 'XTick', [0:6:24], 'YTick', [0:6:24])
axis([0 24 0 24])
grid on

% plot(t(idx), J_A(idx), t(idx), J_B(idx))
cmap = brewermap(128, 'RdBu');
colormap(flipud(cmap));
axis square

set(gcf, 'PaperOrientation', 'Portrait', 'PaperType', 'A4', 'PaperUnits', ...
    'centimeters', 'PaperPosition', [0.2 0.2 [15 15]]);
filename_out = 'correlation_EEJ_A_B.png';
% print('-dpng', '-r200', filename_out); crop_png(filename_out);

%% model values SW-A
[ab, R, QD_lat, params] = read_EEJM('EEJ_SW-AB_03.txt');
m_i = ab(:,k_QD);
J_mod = zeros(size(t));
n = 1;
for p = 0:params.n_t
    for m = -params.n_phi:params.n_phi
        for k = -params.n_s:params.n_s
            if params.n_l == 1
                l_sel = 0;
            elseif params.n_l == 2
                l_sel = [-2 0];
            end
            for l = l_sel
                % coefficients a
                if ~( (p == 0 && l == 0 && m < 0) || (p == 0 && l == 0 && m == 0 && k < 0) )
                    J_mod = J_mod + m_i(n)*cosd(p*t*360 + m*phi + k*t*360/365.25 + l*nu*15);
                    n = n + 1;
                end
                % coefficients b
                if ~( (m < 0 && p == 0 && l == 0) || (p == 0 && m == 0 && l == 0 && k < 0) || (p == 0 && m == 0 && l == 0 && k == 0) )
                    J_mod = J_mod + m_i(n)*sind(p*t*360 + m*phi + k*t*360/365.25 + l*nu*15);
                    n = n + 1;
                end
            end % l
        end % k
    end % m
end % p
J_mod_A = (1+R(k_QD)*EUVAC).*J_mod;

% model values SW-B
phi_B = mod(T_B/24 - t, 1)*360;
J_mod = zeros(size(t));
n = 1;
for p = 0:params.n_t
    for m = -params.n_phi:params.n_phi
        for k = -params.n_s:params.n_s
            if params.n_l == 1
                l_sel = 0;
            elseif params.n_l == 2
                l_sel = [0 -2];
            end
            for l = l_sel
                % coefficients a
                if ~( (p == 0 && l == 0 && m < 0) || (p == 0 && l == 0 && m == 0 && k < 0) )
                    J_mod = J_mod + m_i(n)*cosd(p*t*360 + m*phi_B + k*t*360/365.25 + l*nu*15);
                    n = n + 1;
                end
                % coefficients b
                if ~( (m < 0 && p == 0 && l == 0) || (p == 0 && m == 0 && l == 0 && k < 0) || (p == 0 && m == 0 && l == 0 && k == 0) )
                    J_mod = J_mod + m_i(n)*sind(p*t*360 + m*phi_B + k*t*360/365.25 + l*nu*15);
                    n = n + 1;
                end
            end % l
        end % k
    end % m
end % p
J_mod_B = (1+R(k_QD)*EUVAC).*J_mod;

%%
fac = J_B./J_mod_B;

LT = [0:24]';
rho = zeros(length(LT));
for i_LT_A = 1:length(LT)
    for i_LT_B = 1:length(LT)
        idx = abs(T_A - LT(i_LT_A)) < .5 & abs(T_B - LT(i_LT_B)) < .5;
        rho(i_LT_A, i_LT_B) = corr(J_A(idx), J_mod_A(idx));
    end
end

imagesc(LT, LT, rho)
caxis([-1 1])
colorbar
xlabel 'LT of SW-A [hours]'
ylabel 'LT of SW-B [hours]'
set(gca, 'XTick', [0:6:24], 'YTick', [0:6:24])
axis([0 24 0 24])
grid on

% plot(t(idx), J_A(idx), t(idx), J_B(idx))
cmap = brewermap(128, 'RdBu');
colormap(flipud(cmap));
axis square

set(gcf, 'PaperOrientation', 'Portrait', 'PaperType', 'A4', 'PaperUnits', ...
    'centimeters', 'PaperPosition', [0.2 0.2 [15 15]]);
filename_out = 'correlation_EEJ_A_B.png';
% print('-dpng', '-r200', filename_out); crop_png(filename_out);

%% correlation between observed and model values as a function of LT
clear all
rad = pi/180;

[ab, R, QD_lat, params] = read_EEJM('EEJ_SW-AB_01.txt');
k_QD = 41; % dip equator
m_i = ab(:,k_QD);

% define some colors
co = [        0    0.4470    0.7410
    0.8500    0.3250    0.0980
    0.9290    0.6940    0.1250
    0.4940    0.1840    0.5560
    0.4660    0.6740    0.1880
    0.3010    0.7450    0.9330
    0.6350    0.0780    0.1840];

load EEJ_A_extended.mat
% t = t(1:2:end); phi = phi(1:2:end); J = J(1:2:end,:); EUVAC = EUVAC(1:2:end); nu = nu(1:2:end); 
% only take every seciond value (no mixing of ascending and descending tracks)
t = [t(1:2:end); t(2:2:end)]; phi = [phi(1:2:end); phi(2:2:end)]; J = [J(1:2:end,:); J(2:2:end,:)]; 
EUVAC = [EUVAC(1:2:end); EUVAC(2:2:end)]; nu = [nu(1:2:end); nu(2:2:end)]; 


if size(J,2) > 1; J = J(:,41); end
idx = t > md2000(2023, 2, 11);
t(idx) = []; phi(idx) = []; J(idx) = []; nu(idx) = []; EUVAC(idx) = [];

% remove outliers
q_1 = prctile(J,1);
q_99 = prctile(J,99);
dq = q_99 - q_1;
idx = q_1-1.5*dq < J & J < q_99 + 1.5*dq & ~isnan(J);
t = t(idx); phi = phi(idx); J = J(idx); nu = nu(idx); EUVAC = EUVAC(idx);
disp([(length(idx)-sum(idx))/length(idx)*100])


T = mod(t+phi/360, 1)*24;

%
J_mod = zeros(size(t));
n = 1;
for p = 0:params.n_t
    for m = -params.n_phi:params.n_phi
        for k = -params.n_s:params.n_s
            if params.n_l == 1
                l_sel = 0;
            elseif params.n_l == 2
                l_sel = [-2 0];
            end
            for l = l_sel
                % coefficients a
                if ~( (p == 0 && l == 0 && m < 0) || (p == 0 && l == 0 && m == 0 && k < 0) )
                    J_mod = J_mod + m_i(n)*cosd(p*t*360 + m*phi + k*t*360/365.25 + l*nu*15);
                    n = n + 1;
                end
                % coefficients b
                if ~( (m < 0 && p == 0 && l == 0) || (p == 0 && m == 0 && l == 0 && k < 0) || (p == 0 && m == 0 && l == 0 && k == 0) )
                    J_mod = J_mod + m_i(n)*sind(p*t*360 + m*phi + k*t*360/365.25 + l*nu*15);
                    n = n + 1;
                end
            end % l
        end % k
    end % m
end % p

J_mod = (1+R(k_QD)*EUVAC).*J_mod;


%%
clf
i_col = 0;
% for N = [1 3 15 31 10*15+1]
for N = [1 3 7 15 31]
    i_col = i_col + 1;
    J_tmp = zp_filter(J, 'boxcar', N);
    J_mod_tmp = zp_filter(J_mod, 'boxcar', N);
    j_mod_tmp = J_mod;
    LT = [0.5:24]';
    clear rho MAD_LT;
    for i_LT = 1:length(LT)
        idx = ~isnan(J_tmp) & abs(T - LT(i_LT)) < .5; % & EUVAC < 81;
        % [b,stats] = robustfit(J(idx), J_mod(idx));
        rho(i_LT) = corr(J_tmp(idx), J_mod_tmp(idx));
        MAD_LT(i_LT) = mad(J_tmp(idx)-J_mod_tmp(idx));
    end
    subplot(1,2,1)
    hold on
    plot(LT, rho, 'LineWidth',2, 'Color', co(i_col,:))
    plot(LT, rho, '.', 'MarkerSize',30, 'Color', co(i_col,:))
grid on
xlabel 'LT [hours'
ylabel 'correlation \rho'
set(gca, 'XTick', [0:6:24])
axis([0 24 0 1])
    subplot(1,2,2)
    hold on
   h(i_col) = plot(LT, MAD_LT, 'LineWidth',2, 'Color', co(i_col,:))
    plot(LT, MAD_LT, '.', 'MarkerSize',30, 'Color', co(i_col,:))
end
hold off
grid on
xlabel 'LT [hours'
ylabel 'MAD [A/km]'
set(gca, 'XTick', [0:6:24])
axis([0 24 0 30 ])
legend([h(1) h(2) h(3) h(4) h(5)], 'N = 1', 'N = 3', 'N = 7', 'N = 15', 'N = 30', 'Location','northeast')

set(gcf, 'PaperOrientation', 'Portrait', 'PaperType', 'A4', 'PaperUnits', ...
    'centimeters', 'PaperPosition', [0.2 0.2 [25 15]]);
filename_out = 'correlation_MAD_EEJ_LT.png';
print('-dpng', '-r200', filename_out); crop_png(filename_out);

set(gcf, 'PaperOrientation', 'Portrait', 'PaperType', 'A4', 'PaperUnits', ...
    'centimeters', 'PaperPosition', [0.2 0.2 [35 12]]);
filename_out = 'correlation_MAD_EEJ_LT.eps';
print_pdf(filename_out)

%%
clear all


% define some colors
co = [        0    0.4470    0.7410
    0.8500    0.3250    0.0980
    0.9290    0.6940    0.1250
    0.4940    0.1840    0.5560
    0.4660    0.6740    0.1880
    0.3010    0.7450    0.9330
    0.6350    0.0780    0.1840];

[ab, R, QD_lat, params] = read_EEJM('EEJ_SW-AB_01.txt');
k_QD = 41;
m_i = ab(:,k_QD);

T_SWA = @(t) mod(471.24 -0.089947*t, 24); % Swarm Alpha
T_SWB = @(t) mod(448.86 -0.085625*t, 24); % Swarm Bravo

clf
for i_case = 1:4
    if i_case <= 2
        load EEJ_A_extended.mat
    else
        load EEJ_B_extended.mat
    end

    if size(J,2) > 1; J = J(:,41); end
    idx = t > md2000(2023, 2, 11);
    t(idx) = []; phi(idx) = []; J(idx) = []; nu(idx) = []; EUVAC(idx) = [];

    % remove outliers
    q_1 = prctile(J,1);
    q_99 = prctile(J,99);
    dq = q_99 - q_1;
    idx = q_1-1.5*dq < J & J < q_99 + 1.5*dq & ~isnan(J);
    t = t(idx); phi = phi(idx); J = J(idx); nu = nu(idx); EUVAC = EUVAC(idx);
    disp([(length(idx)-sum(idx))/length(idx)*100])

    %
    T = mod(t+phi/360,1)*24;
    if i_case == 1
        idx = abs(T - T_SWA(t)) < 1 | abs(T-T_SWA(t)) > 23;
    elseif i_case == 2
        idx = abs(T - T_SWA(t)) < 1 | abs(T-T_SWA(t)) > 23;
        idx = ~idx;
    elseif i_case == 3
        idx = abs(T - T_SWB(t)) < 1 | abs(T-T_SWB(t)) > 23;
    else
        idx = abs(T - T_SWB(t)) < 1 | abs(T-T_SWB(t)) > 23;
        idx = ~idx;
    end

    t = t(idx); phi = phi(idx); J = J(idx); nu = nu(idx); EUVAC = EUVAC(idx); T = T(idx);

    %
    J_mod = zeros(size(t));
    n = 1;
    for p = 0:params.n_t
        for m = -params.n_phi:params.n_phi
            for k = -params.n_s:params.n_s
                if params.n_l == 1
                    l_sel = 0;
                elseif params.n_l == 2
                    l_sel = [-2 0];
                end
                for l = l_sel
                    % coefficients a
                    if ~( (p == 0 && l == 0 && m < 0) || (p == 0 && l == 0 && m == 0 && k < 0) )
                        J_mod = J_mod + m_i(n)*cosd(p*t*360 + m*phi + k*t*360/365.25 + l*nu*15);
                        n = n + 1;
                    end
                    % coefficients b
                    if ~( (m < 0 && p == 0 && l == 0) || (p == 0 && m == 0 && l == 0 && k < 0) || (p == 0 && m == 0 && l == 0 && k == 0) )
                        J_mod = J_mod + m_i(n)*sind(p*t*360 + m*phi + k*t*360/365.25 + l*nu*15);
                        n = n + 1;
                    end
                end % l
            end % k
        end % m
    end % p

    J_mod = (1+R(k_QD)*EUVAC).*J_mod;

    i_plt = 1;
    for N = [1 3 15 31]
        J_tmp = zp_filter(J, 'boxcar', N);
        J_mod_tmp = zp_filter(J_mod, 'boxcar', N);
        LT = [0.5:24]';
        clear rho MAD_LT;
        for i_LT = 1:length(LT)
            idx = ~isnan(J_tmp) & abs(T - LT(i_LT)) < .5;
            rho(i_LT) = corr(J_tmp(idx), J_mod_tmp(idx),'Type', 'Spearman');
            MAD_LT(i_LT) = mad(J_tmp(idx)-J_mod_tmp(idx));
        end

        subplot(1,2,1)
        hold on
        plot(LT, rho, 'LineWidth',2, 'Color', co(i_plt,:))
        grid on
        xlabel 'LT [hours]'
        ylabel 'correlation \rho'
        set(gca, 'XTick', [0:6:24])
        legend('N = 1', 'N = 3', 'N = 15', 'N = 31', 'Location','southeast')
        axis([0 24 -Inf Inf])
        subplot(1,2,2)
        hold on
        plot(LT, MAD_LT, 'LineWidth',2, 'Color', co(i_plt,:))
        grid on
        xlabel 'LT [hours]'
        ylabel 'MAD [A/km]'
        set(gca, 'XTick', [0:6:24])
        legend('N = 1', 'N = 3', 'N = 15', 'N = 31', 'Location','northeast')
        axis([0 24 -Inf Inf])
        i_plt = i_plt + 1;
    end
end % i_case

%% Coherency between observed and model values
clear all

[ab, R, QD_lat, params] = read_EEJM('EEJ_SW-AB_01.txt');
k_QD = 41;
m_i = ab(:,k_QD);

T_SWA = @(t) mod(471.24 -0.089947*t, 24); % Swarm Alpha
T_SWB = @(t) mod(448.86 -0.085625*t, 24); % Swarm Bravo

load EEJ_A_extended.mat
if size(J,2) > 1; J = J(:,41); end
idx = t > md2000(2023, 2, 11);
t(idx) = []; phi(idx) = []; J(idx) = []; nu(idx) = []; EUVAC(idx) = [];

% remove outliers
q_1 = prctile(J,1);
q_99 = prctile(J,99);
dq = q_99 - q_1;
idx = q_1-1.5*dq < J & J < q_99 + 1.5*dq & ~isnan(J);
t = t(idx); phi = phi(idx); J = J(idx); nu = nu(idx); EUVAC = EUVAC(idx);
disp([(length(idx)-sum(idx))/length(idx)*100])

%
T = mod(t+phi/360,1)*24;
idx = abs(T - T_SWA(t)) < 1 | abs(T-T_SWA(t)) > 23;
% idx = ~idx;
t = t(idx); phi = phi(idx); J = J(idx); nu = nu(idx); EUVAC = EUVAC(idx); T = T(idx);

%
J_mod = zeros(size(t));
n = 1;
for p = 0:params.n_t
    for m = -params.n_phi:params.n_phi
        for k = -params.n_s:params.n_s
            if params.n_l == 1
                l_sel = 0;
            elseif params.n_l == 2
                l_sel = [-2 0];
            end
            for l = l_sel
                % coefficients a
                if ~( (p == 0 && l == 0 && m < 0) || (p == 0 && l == 0 && m == 0 && k < 0) )
                    J_mod = J_mod + m_i(n)*cosd(p*t*360 + m*phi + k*t*360/365.25 + l*nu*15);
                    n = n + 1;
                end
                % coefficients b
                if ~( (m < 0 && p == 0 && l == 0) || (p == 0 && m == 0 && l == 0 && k < 0) || (p == 0 && m == 0 && l == 0 && k == 0) )
                    J_mod = J_mod + m_i(n)*sind(p*t*360 + m*phi + k*t*360/365.25 + l*nu*15);
                    n = n + 1;
                end
            end % l
        end % k
    end % m
end % p

J_mod = (1+R(k_QD)*EUVAC).*J_mod;

%% continuous time series
% calculate coherency
[Qd, dQd, T, wsum, coh2d] = rob_tfe(detrend(J), detrend(J_mod));
T = T*93.4*60; % minutes to seconds

clf
semilogx(T/86400, coh2d, '-o', 'LineWidth',1);
grid on
text(12*3600/86400,   1.01, '12 hours', 'HorizontalAlignment', 'center')
text(86400/86400,   1.01, '1 day', 'HorizontalAlignment', 'center')
text(604800/86400,  1.01, '1 week', 'HorizontalAlignment', 'center')
text(2592000/86400, 1.01, '1 month', 'HorizontalAlignment', 'center')
text(7776000/86400, 1.01, '3 months', 'HorizontalAlignment', 'center')
hold off
% axis([2e4/86400 2.5e7/86400 0.6 1.001])
% set(gca,'XtickLabel',[]')
ylabel 'coh^2'
xlabel 'period [days]'

set(gcf, 'PaperOrientation', 'Portrait', 'PaperType', 'A4', 'PaperUnits', ...
    'centimeters', 'PaperPosition', [0.2 0.2 [20 12]]);
filename_out = 'coh2_Dst_RC_LT.png';
% print('-dpng', '-r200', filename_out); crop_png(filename_out);

%% rms statistics vs LT with EUVAC dependecies
clear all
rad = pi/180;

% [ab, R, QD_lat, params] = read_EEJM('EEJ_SW-AB_01.txt');
[ab, ab_EUVAC, QD_lat, params] = read_EEJM_EUVAC('EEJ_SW-AB_EUVAC_01.txt');

k_QD = 41;
m_i = ab(:,k_QD);
m_i_EUVAC = ab_EUVAC(:,k_QD);

% define some colors
co = [        0    0.4470    0.7410
    0.8500    0.3250    0.0980
    0.9290    0.6940    0.1250
    0.4940    0.1840    0.5560
    0.4660    0.6740    0.1880
    0.3010    0.7450    0.9330
    0.6350    0.0780    0.1840];

% read daily values of F10.7
filename = 'SW_OPER_AUX_F10_2__20060101T000000_20231015T000000_0001.DBL';
[t_F107_all, F107_all]  = textread(filename, '%f %f', 'commentstyle', 'shell');

load EEJ_A_extended.mat

if size(J,2) > 1; J = J(:,41); end
idx = t > md2000(2023, 2, 11);
t(idx) = []; phi(idx) = []; J(idx) = [];

% remove outliers
q_1 = prctile(J,1);
q_99 = prctile(J,99);
dq = q_99 - q_1;
idx = q_1-1.5*dq < J & J < q_99 + 1.5*dq & ~isnan(J);
t = t(idx); phi = phi(idx); J = J(idx); % nu = nu(idx); EUVAC = EUVAC(idx);

idx = ~isnan(F107_all);
F107 = interp1(t_F107_all(idx), zp_filter(F107_all(idx), 'boxcar', 81), t, 'linear'); % interpolate F107 to time instants of "synthetic" data
EUVAC = (F107 + interp1(t_F107_all(idx), F107_all(idx), t, 'linear'))/2;
nu = moon_phase(t);
T = mod(t+phi/360, 1)*24;

%
J_mod = zeros(size(t));
n = 1;
for p = 0:params.n_t
    for m = -params.n_phi:params.n_phi
        for k = -params.n_s:params.n_s
            if params.n_l == 1
                l_sel = 0;
            elseif params.n_l == 2
                l_sel = [-2 0];
            end
            for l = l_sel
                % coefficients a
                if ~( (p == 0 && l == 0 && m < 0) || (p == 0 && l == 0 && m == 0 && k < 0) )
                    J_mod = J_mod + m_i(n)*cosd(p*t*360 + m*phi + k*t*360/365.25 + l*nu*15);
                    J_mod = J_mod + m_i_EUVAC(n)*EUVAC.*cosd(p*t*360 + m*phi + k*t*360/365.25 + l*nu*15);
                    n = n + 1;
                end
                % coefficients b
                if ~( (m < 0 && p == 0 && l == 0) || (p == 0 && m == 0 && l == 0 && k < 0) || (p == 0 && m == 0 && l == 0 && k == 0) )
                    J_mod = J_mod + m_i(n)*sind(p*t*360 + m*phi + k*t*360/365.25 + l*nu*15);
                    J_mod = J_mod + m_i_EUVAC(n)*EUVAC.*sind(p*t*360 + m*phi + k*t*360/365.25 + l*nu*15);
                    n = n + 1;
                end
            end % l
        end % k
    end % m
end % p

% statistics
delta_d = [J - J_mod];


rms_LT_A = zeros(24,1);
for LT = 1:24
    idx = abs(LT - T - .5) < .5;
    rms_LT_A(LT) = mad(delta_d(idx));
end

%
clf
plot([.5:24], rms_LT_A, '.', 'MarkerSize', 30)
hold on
h1 = plot([.5:24], rms_LT_A, 'LineWidth', 2, 'Color', co(1,:))
hold off
axis([0 24 0 40])
grid on
set(gca, 'XTick', [0:3:24])
xlabel 'LT [hours]'
ylabel 'MAD [A/km]'

%%
tmp = load('EEJ_A_extended_EEJM2.txt');
idx = tmp(:,5) < -9e3;
tmp(idx,:) = [];
EEJM2.t = tmp(:,1);
EEJM2.phi = tmp(:,2);
EEJM2.EUVAC = tmp(:,3);
EEJM2.nu = tmp(:,4);
EEJM2.J = tmp(:,5)*1e3;
EEJM2.T = mod(EEJM2.t+EEJM2.phi/360,1)*24;

EEJM2.J_obs = interp1(t, J, EEJM2.t);

rms_LT = NaN*ones(24,1);
delta_d = EEJM2.J_obs - EEJM2.J;
for LT = 1:24
    idx = abs(LT - EEJM2.T - .5) < .5;
    if sum(idx) > 10
    rms_LT(LT) = mad(delta_d(idx));
    end
end

hold on
h2 = plot([.5:24], rms_LT, 'LineWidth', 2, 'Color', co(2,:))
plot([.5:24], rms_LT, '.', 'Markersize', 30, 'Color', co(2,:))
hold off

%% add rms with fixed R

[ab, R, QD_lat, params] = read_EEJM('EEJ_SW-AB_01.txt');
m_i = ab(:,k_QD);

%
J_mod = zeros(size(t));
n = 1;
for p = 0:params.n_t
    for m = -params.n_phi:params.n_phi
        for k = -params.n_s:params.n_s
            if params.n_l == 1
                l_sel = 0;
            elseif params.n_l == 2
                l_sel = [-2 0];
            end
            for l = l_sel
                % coefficients a
                if ~( (p == 0 && l == 0 && m < 0) || (p == 0 && l == 0 && m == 0 && k < 0) )
                    J_mod = J_mod + m_i(n)*cosd(p*t*360 + m*phi + k*t*360/365.25 + l*nu*15);
                    n = n + 1;
                end
                % coefficients b
                if ~( (m < 0 && p == 0 && l == 0) || (p == 0 && m == 0 && l == 0 && k < 0) || (p == 0 && m == 0 && l == 0 && k == 0) )
                    J_mod = J_mod + m_i(n)*sind(p*t*360 + m*phi + k*t*360/365.25 + l*nu*15);
                    n = n + 1;
                end
            end % l
        end % k
    end % m
end % p

J_mod = (1+R(k_QD)*EUVAC).*J_mod;

% statistics
delta_d = [J - J_mod];

rms_LT = zeros(24,1);
for LT = 1:24
    idx = abs(LT - T - .5) < .5;
    rms_LT(LT) = mad(delta_d(idx));
end

hold on
h3 = plot([.5:24], rms_LT, 'LineWidth', 2, 'Color', co(6,:))
plot([.5:24], rms_LT, '.', 'Markersize', 30, 'Color', co(6,:))
hold off


legend([h3 h2 h1], 'Swarm AB', 'EEJM-2.0', 'Swarm AB (individual EUVAC dependence)')

set(gcf, 'PaperOrientation', 'Portrait', 'PaperType', 'A4', 'PaperUnits', ...
    'centimeters', 'PaperPosition', [0.2 0.2 [20 12]]);
filename_out = 'MAD_LT.png';
print('-dpng', '-r200', filename_out); crop_png(filename_out);
