clear all
rad = pi/180;
a = 6371.2;

fclose 'all'

% output files
filename_out_A = 'SWA_B_01.quiet';
filename_out_F_A = 'SWA_F_01.quiet';
delete(filename_out_A);
delete(filename_out_F_A);
i_sat = 3; % Swarm-A

% where are the CDF files?
filedir_CDF = '/hpc2/users/Swarm/CalVal6/';

filename_apex_grid = '/home/nio/m/tools/apexsh_1980-2025.txt';

% read Kp (file availabe at ftp.gfz-potsdam.de/pub/incoming/champ_payload/Nils/Kp_MJD_1998_QL.dat )
filename_Kp = '/r2/solsys/indices/Kp_MJD_1998_QL.dat';
[t_Kp, Kp_all]  = textread(filename_Kp, '%f %d %*d  %*s', 'commentstyle', 'shell');

% read Dst index (file available at http://www.spacecenter.dk/files/magnetic-models/RC/current/RC_1997-2022_augmented.dat )
filename_Dst = '/nfs/wwwfiles/projects/oersted/models/RC/current/RC_1997-2022_augmented.dat';
[t_Dst, Dst_all, Est_all, Ist_all]  = textread(filename_Dst, '%f %f %f %f %*c', 'commentstyle', 'shell'); % new version, with flag in 4th column
index = isnan(Dst_all);
t_Dst(index) = []; Dst_all(index) = [];
dDst_all = diff(Dst_all);
dDst_all = [dDst_all(1); dDst_all];


% read Swarm Alpha
year = 2021;
Step = 60; % downsample data to 60 sec
for month = 1:12
    for day = 1:31
        d_A = dir([filedir_CDF 'SW_OPER_MAGA_LR_1B_' num2str(year, '%04d') num2str(month, '%02d') num2str(day, '%02d') '*.cdf']);
        if ~isempty(d_A) % does file exist?
            filename_in = [filedir_CDF d_A(1).name];
            disp(filename_in)
            data = cdfread(filename_in, 'CombineRecords', 1, 'variable', {'Timestamp', 'Latitude', 'Longitude', ...
                'Radius', 'F', 'B_VFM', 'q_NEC_CRF', 'Flags_B' 'Flags_q'});
            t = todatenum(data{1}(1:Step:end)) - 730486;
            r = double(data{4}(1:Step:end))/1e3; % radius in km
            theta = 90-double(data{2}(1:Step:end));
            phi = double(data{3}(1:Step:end));
            F = double(data{5}(1:Step:end,:));
            B_VFM = double(data{6}(1:Step:end,:));
            q_NEC_CRF = double(data{7}(1:Step:end,:));
            Flags_B = double(data{8}(1:Step:end));
            Flags_q = double(data{9}(1:Step:end));
            source = i_sat*ones(size(t));
            
            % interpolate Kp and dDst to present time instants
            dDst = interp1(t_Dst, dDst_all, t, 'linear');
            Kp = interp1(t_Kp, Kp_all, t, 'linear'); % Kp of present time instant
            
            % dDst and Kp selection criteria
            index = find((abs(dDst) > 3) | (Kp > 30));
            t(index) = []; r(index) = []; theta(index) = []; phi(index) = []; B_VFM(index,:) = []; F(index) = []; q_NEC_CRF(index,:) = []; Flags_B(index) = []; Flags_q(index) = []; source(index) = [];
            
            % remove errortneous data
            index = F < 10000 | sqrt(B_VFM(:,1).^2 + B_VFM(:,2).^2 + B_VFM(:,3).^2) < 10000;
            t(index) = []; r(index) = []; theta(index) = []; phi(index) = []; B_VFM(index,:) = []; F(index) = []; q_NEC_CRF(index,:) = []; Flags_B(index) = []; Flags_q(index) = []; source(index) = [];
            
            % calculate (cosine of) solar zenith angle
            [~, decl] = sun_md2000(t);
            cos_zeta = cos(theta*rad).*sin(decl)+sin(theta*rad).*cos(decl).*cos(mod(t+.5,1)*2*pi+phi*rad);
            cos_zeta_0 = cosd(90-(-10)); % Threshold for defining sunlit or dark data, according to solar inclination angle
            % remove data from sunlit regions
            index = find(cos_zeta > cos_zeta_0);
            t(index) = []; r(index) = []; theta(index) = []; phi(index) = []; B_VFM(index,:) = []; F(index) = []; q_NEC_CRF(index,:) = []; Flags_B(index) = []; Flags_q(index) = []; source(index) = [];
            
            
            % calculate QD latitude
            QD_lat = qdipole(t, r/a, theta*rad, phi*rad, filename_apex_grid);
            
            % save scalar data
            fid_F_out = fopen(filename_out_F_A, 'a');
            
            fprintf(fid_F_out, '%s%s\n', '% Swarm scalar data, file created on: ', datestr(now));
            fprintf(fid_F_out, '%s%s\n', '% filename_Dst: ', filename_Dst);
            fprintf(fid_F_out, '%s%s\n', '% filename_Kp: ', filename_Kp);
            fprintf(fid_F_out, '%s%s\n', '% CDF input file: ', filename_in);
            fprintf(fid_F_out, '%s\n', '% ');
            
            index = find(abs(QD_lat) > 60);
            fprintf(fid_F_out,'%15.8f %9.3f %11.5f %11.5f %10.3f %2d\n',...
                [t(index) r(index) theta(index) phi(index) F(index) source(index)]');
            % fprintf(fid_F_out,'%15.8f %9.3f %11.5f %11.5f %10.3f %2d\n',...
            %     [t r theta phi F source]');
            fclose(fid_F_out);
            
            
            % save vector data
            q_nec = dcm2q(shiftdim(ecef_nec(theta, phi), 1));	% = q_nec_itrf
            q_ecef = zeros(size(q_nec)); q_ecef(:,4) = 1;
            q_j2000 = qmult(qmult(qconj(q_ecef), q_nec), qconj(q_NEC_CRF));	% = q_crf_icrf
            
            % remove high-latitude data and data with NaN
            index = isnan(q_j2000(:,1)) | abs(QD_lat) > 60;
            t(index) = []; r(index) = []; theta(index) = []; phi(index) = [];
            B_VFM(index,:) = []; q_j2000(index,:) = []; q_ecef(index,:) = [];
            Flags_B(index,:) = []; Flags_q(index,:) = [];
            
            fid_out = fopen(filename_out_A, 'a');
            
            fprintf(fid_out, '%s%s\n', '% Swarm vector data, file created on: ', datestr(now));
            fprintf(fid_out, '%s%s\n', '% filename_Dst: ', filename_Dst);
            fprintf(fid_out, '%s%s\n', '% filename_Kp: ', filename_Kp);
            fprintf(fid_out, '%s%s\n', '% CDF input file: ', filename_in);
            fprintf(fid_out, '%s\n', '% ');
            fprintf(fid_out,'%15.8f %9.3f %11.5f %11.5f %10.3f %10.3f %10.3f %10.6f %10.6f %10.6f %10.6f 0 0 0 1 %4.4d %4.4d\n',...
                [t r theta phi B_VFM q_j2000 Flags_B Flags_q]');
            fclose(fid_out);
            
        end % CDF file exist
    end % day
end % month
