%% read all Swarm EEJ data of one satellite
clear all

d = dir('./Sat_A_extended/*.cdf');
% d = dir('./Sat_A/*.cdf');
% d = dir('./CSES/*.cdf');
phi = [];
t = [];
J = [];
for k = 1:length(d)
    filename_cdf = [d(k).folder '\' d(k).name];
    data = cdfread(filename_cdf, 'CombineRecords', 1, 'variable', {'Timestamp', 'Latitude', 'Longitude', 'EEJ'});
    t_tmp = todatenum(data{1}) - 730486;
    phi_tmp = double(data{3});
    J_tmp = double(data{4});

    if size(J_tmp, 1) == 81
        J_tmp = J_tmp';
    end

    J_tmp(J_tmp == 0) = NaN;

    t = [t; t_tmp];
    phi = [phi; phi_tmp];
    J = [J; J_tmp];
end

% sort in time
[t, idx] = sort(t);
phi = phi(idx); J = J(idx,:);

% remove double values
idx = diff(t) <= 0;
t(idx) = []; phi(idx) = []; J(idx,:) = [];

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

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); % lunar phase in hours (0 - 24 hours)

save EEJ_A_extended.mat t J phi EUVAC nu


