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

% define synthetic orbit
t_0 = 0; % MD2000, t = 0 corresponds to 1 Jan 2000 00 UT
t = [0:1/24:30]' + t_0; % 30 days of hourly values starting on 1 Jan 2000

theta = 90*ones(size(t)); % co-latitude in deg
phi = 0*ones(size(t)); % longitude in deg
r = (a + 0)*ones(size(t)); % radius in km. Here at Earth's surface


% read MIO model
filename = 'SW_OPER_MIO_SHAi2C_20131201T000000_20170101T000000_0201.txt';
filename = 'MIO_CM6.DBL';
s_MIO = SwarmL2_MIO_SHA_Read(filename);

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


%% ionospheric field from MIO
F107 = interp1(t_F107_all, F107_all, t, 'linear'); % interpolate F107 to time instants of "synthetic" data
[B_ion_1, B_ion_2] = synth_values_MIO(r, theta, phi, t, F107, s_MIO);

%% plot 30 days time series
for k = 1:3
    subplot(3,1,k)
    plot(t, B_ion_1(:,k), t, B_ion_2(:,k))
    if k == 1
        ylabel '\DeltaB_r [nT]'
    elseif k == 2 
        ylabel '\DeltaB_\theta [nT]'
    elseif k == 3
            ylabel '\DeltaB_\phi [nT]'
    end
    legend('primary (ionospheric)', 'secondary (induced')
end
    
set(gcf, 'PaperOrientation', 'Portrait', 'PaperType', 'A4', 'PaperUnits', ...
    'centimeters', 'PaperPosition', [0.2 0.2 .5*[40 25]]);
% filename_out = ['Test_MIO.png']; print('-dpng', filename_out); crop_png(filename_out)
