%% read CDF files
clear all
a = 6371.2;
rad = pi/180;

filename_in = './CDF/CS_OPER_MAG_20110105T000000_20110105T235959_0101.cdf';


data = cdfread(filename_in, 'CombineRecords', 1, 'variable', {'Timestamp', 'Latitude', 'Longitude', ...
    'Radius', 'B_mod_NEC', 'B_FGM1', 'B_FGM2', 'B_FGM3', 'B_NEC1', 'B_NEC2', 'B_NEC3', 'q_NEC_CRF', 'q_error'});
t = todatenum(data{1}) - 730486;
r = double(data{4})/1e3; % radius in km
theta = 90-double(data{2});
phi = double(data{3});
B_mod_NEC = double(data{5});
B_FGM1 = double(data{6});
B_FGM2 = double(data{7});
B_FGM3 = double(data{8});
B_NEC1 = double(data{9});
B_NEC2 = double(data{10});
B_NEC3 = double(data{11});
q_NEC_CRF = double(data{12});
q_error = double(data{13});

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

%% plot in dependence on geographic latitude
B_NEC = B_NEC3; % data from FGM3

F = sqrt(B_NEC(:,1).^2 + B_NEC(:,2).^2 + B_NEC(:,3).^2);
F_mod = sqrt(B_mod_NEC(:,1).^2 + B_mod_NEC(:,2).^2 + B_mod_NEC(:,3).^2);

index_day = q_error < 40 & (T > 6 & T < 18);
index_night = q_error < 40 & (T > 18 | T < 6);
subplot(4,1,1)
plot(90-theta(index_day), B_NEC(index_day,1) - B_mod_NEC(index_day,1), '.r', ...
    90-theta(index_night), B_NEC(index_night,1) - B_mod_NEC(index_night,1), '.b')
ylabel '\DeltaB_N [nT]'
axis([-90 90 -500 500])
grid on
subplot(4,1,2)
plot(90-theta(index_day), B_NEC(index_day,2) - B_mod_NEC(index_day,2), '.r', ...
    90-theta(index_night), B_NEC(index_night,2) - B_mod_NEC(index_night,2), '.b')
ylabel '\DeltaB_E [nT]'
axis([-90 90 -500 500])
grid on
subplot(4,1,3)
plot(90-theta(index_day), B_NEC(index_day,3) - B_mod_NEC(index_day,3), '.r', ...
    90-theta(index_night), B_NEC(index_night,3) - B_mod_NEC(index_night,3), '.b')
ylabel '\DeltaB_C [nT]'
grid on
axis([-90 90 100*[-1 1]])
subplot(4,1,4)
plot(90-theta(index_day), F(index_day) - F_mod(index_day,1), '.r', ...
    90-theta(index_night), F(index_night) - F_mod(index_night), '.b')
ylabel '\DeltaF [nT]'
grid on
xlabel 'geographic latitude [deg]'
axis([-90 90 100*[-1 1]])
legend('dayside','nightside')

set(gcf, 'PaperOrientation', 'Portrait', 'PaperType', 'A4', 'PaperUnits', 'centimeters', 'PaperPosition', [0.2 0.2 [40 25]]);
print('-dpng', 'B_NEC_lat.png')



%% plot in dependence on QD latitude
QD_lat = qdipole_old(t, r/a, theta*rad, phi*rad, 'h:\m\tools\apexsh_1980-2020.txt');

subplot(4,1,1)
plot(QD_lat(index_day), B_NEC(index_day,1) - B_mod_NEC(index_day,1), '.r', ...
    QD_lat(index_night), B_NEC(index_night,1) - B_mod_NEC(index_night,1), '.b')
ylabel '\DeltaB_N [nT]'
axis([-90 90 -500 500])
grid on
subplot(4,1,2)
plot(QD_lat(index_day), B_NEC(index_day,2) - B_mod_NEC(index_day,2), '.r', ...
    QD_lat(index_night), B_NEC(index_night,2) - B_mod_NEC(index_night,2), '.b')
ylabel '\DeltaB_E [nT]'
axis([-90 90 -500 500])
grid on
subplot(4,1,3)
plot(QD_lat(index_day), B_NEC(index_day,3) - B_mod_NEC(index_day,3), '.r', ...
    QD_lat(index_night), B_NEC(index_night,3) - B_mod_NEC(index_night,3), '.b')
ylabel '\DeltaB_C [nT]'
grid on
axis([-90 90 100*[-1 1]])
subplot(4,1,4)
plot(QD_lat(index_day), F(index_day) - F_mod(index_day,1), '.r', ...
    QD_lat(index_night), F(index_night) - F_mod(index_night), '.b')
ylabel '\DeltaF [nT]'
grid on
xlabel 'QD latitude [deg]'
axis([-90 90 100*[-1 1]])
legend('dayside','nightside')

set(gcf, 'PaperOrientation', 'Portrait', 'PaperType', 'A4', 'PaperUnits', 'centimeters', 'PaperPosition', [0.2 0.2 [40 25]]);
print('-dpng', 'B_NEC_QD-lat.png')
