%% read SV and flow coefficients and save in *.mat format
clear all
a = 6371.2;
rad = pi/180;

% read SV predictions
gh_1 = NaN*ones(224, 187);
t_1 = NaN*ones(187, 1);
for k = 1:187
    filename_in = ['./anisotropic_errors/svpred.rmm.e-4.anis.ep' num2str(k, '%4.4d') '.mon.1'];
    disp([k exist(filename_in, 'file')])
    if exist(filename_in, 'file')
        fid = fopen(filename_in, 'r');
        t_1(k) = sscanf(fgetl(fid), '%f'); maxdeg = fgetl(fid);
        tmp = []; while ~feof(fid); tmp = [tmp fgetl(fid)]; end
        gh_1(:,k) = sscanf(tmp, '%f');
        fclose(fid);
    end
end

gh_2 = NaN*ones(224, 187);
t_2 = NaN*ones(187, 1);
for k = 1:187
    filename_in = ['./anisotropic_errors/svpred.rmm.e-4.anis.ep' num2str(k-1, '%4.4d') '.mon.2'];
    disp([k exist(filename_in, 'file')])
    if exist(filename_in, 'file')
        fid = fopen(filename_in, 'r');
        t_2(k) = sscanf(fgetl(fid), '%f'); maxdeg = fgetl(fid);
        tmp = []; while ~feof(fid); tmp = [tmp fgetl(fid)]; end
        gh_2(:,k) = sscanf(tmp, '%f');
        fclose(fid);
    end
end

gh_3 = NaN*ones(224, 187);
t_3 = NaN*ones(187, 1);
for k = 1:187
    filename_in = ['./anisotropic_errors/svpred.rmm.e-4.anis.ep' num2str(k-2, '%4.4d') '.mon.3'];
    disp([k exist(filename_in, 'file')])
    if exist(filename_in, 'file')
        fid = fopen(filename_in, 'r');
        t_3(k) = sscanf(fgetl(fid), '%f'); maxdeg = fgetl(fid);
        tmp = []; while ~feof(fid); tmp = [tmp fgetl(fid)]; end
        gh_3(:,k) = sscanf(tmp, '%f');
        fclose(fid);
    end
end

% read flow coefficients
ts = NaN*ones(448, 187);
for k = 1:187
    filename_in = ['./anisotropic_errors/flow.e-4.rmm.epoch' num2str(k, '%4.4d')];
    if exist(filename_in, 'file')
        fid = fopen(filename_in, 'r');
        tmp = fgetl(fid); tmp = fgetl(fid);
        tmp = []; while ~feof(fid); tmp = [tmp fgetl(fid)]; end
        ts(:,k) = sscanf(tmp, '%f');
        fclose(fid);
    end
end
t_coeff = ts(1:224,:);
s_coeff = ts(225:end,:);

save 'svpred_rmm_e-4_anis_ep.mat' t_1 t_2 t_3 gh_1 gh_2 gh_3 t_coeff s_coeff

%% Plot SV predictions at observatory locations and compare with observations and CHAOS-4
clear all
a = 6371.2;
rad = pi/180;

% load observed RMM
load \nio\m\Swarm\SCARF\WP16000\dB_obs_10a.mat

% list of observatories to be plotted
index = [1; find(abs(diff(theta_obs_all)) > 0)];
obs_mne = obs_all(index,:);

% list of observatories to be plotted
obs_mne = ['NGK'; 'HER'; 'CNB'; 'KAK'; 'MBO'; 'CLF'];

% load CHAOS-4 model
load \nio\m\mmod\CHAOS-01a_14.mat
d_pp_CHAOS = fnder(pp,1);

% load flow and SV predictions
load 'svpred_rmm_e-4_anis_ep.mat'

filename_plot = 'dB_rmm_anis.eps';
delete(filename_plot)

t_2(1) = t_2(2)-1/12;
t_2(end) = t_2(end-1)+1/12;

% loop over all observatories
for i_obs = 1:size(obs_mne, 1)
    obs = obs_mne(i_obs,1:3);
    
    index = find(obs_all(:,1) == obs(1) & obs_all(:,2) == obs(2) & obs_all(:,3) == obs(3));
    t_tmp = t_obs_all(index)/365.25+2000;
    
    % observations
    dX = -dB_obs_all(index,2); dY = dB_obs_all(index,3); dZ = -dB_obs_all(index,1);
    
    % CHAOS-4 model prediction
    dB_mod = synth_values(a, theta_obs_all(index)/rad, phi_obs_all(index)/rad, d_pp_CHAOS, (t_tmp-2000)*365.25)*365.25;
    dB_mod = [-dB_mod(:,2) dB_mod(:,3) -dB_mod(:,1)];
    
    % prediction from flow model
    for k = 1:length(t_1)
        dB_flow(k,:) = synth_values(a, theta_obs_all(index(1))/rad, phi_obs_all(index(1))/rad, gh_1(:,k));
    end
    dB_flow_1 = [-dB_flow(:,2) dB_flow(:,3) -dB_flow(:,1)];
    
    for k = 1:length(t_2)
        dB_flow(k,:) = synth_values(a, theta_obs_all(index(1))/rad, phi_obs_all(index(1))/rad, gh_2(:,k));
    end
    dB_flow_2 = [-dB_flow(:,2) dB_flow(:,3) -dB_flow(:,1)];
    
    for k = 1:length(t_3)
        dB_flow(k,:) = synth_values(a, theta_obs_all(index(1))/rad, phi_obs_all(index(1))/rad, gh_3(:,k));
    end
    dB_flow_3 = [-dB_flow(:,2) dB_flow(:,3) -dB_flow(:,1)];
    
    
    index = find(isnan(dX) == 1);
    dX(index) = []; dY(index) = []; dZ(index) = []; t_tmp(index) = [];
    
    for iplt = 1:3
        subplot(1,3,iplt)
        
        for k = 1:1
            if iplt == 1
                ylabel_text = '{\itdX/dt} [nT/yr]';
                b = dX;
            elseif iplt == 2
                ylabel_text = '{\itdY/dt} [nT/yr]';
                b = dY;
            elseif iplt == 3
                ylabel_text = '{\itdZ/dt} [nT/yr]';
                b = dZ;
            end
            
            % plot CHAOS-4 model values
            plot(t_tmp, dB_mod(:,iplt), '.b', 'Linewidth', 2) % CHAOS-4
            hold on
            % plot observations
            plot(t_tmp, b, '.k', 'Markersize', 10) % RMM observation
            % plot flow predictions
            plot(t_1, dB_flow_1(:,iplt), '.r', 'Markersize', 10) % flow prediction
            plot(t_2, dB_flow_2(:,iplt), '.y', 'Markersize', 10) % flow prediction
            plot(t_3, dB_flow_3(:,iplt), '.g', 'Markersize', 10) % flow prediction
            hold off
            grid on
            xlabel 'year'
            ylabel(ylabel_text)
            if iplt == 1
                text(.80, .965, [obs ], 'Fontsize', 16, 'Color', 'r', ...
                    'EdgeColor','w', 'Backgroundcolor', 'w', 'VerticalAlignment', 'top', 'Units', 'normalized', 'Fontangle', 'italic')
            end
            if k == 1;
                % y-axis scaling based on CHAOS-4 model values
                b_min(iplt) = min(b);
                b_max(iplt) = max(b);
            end
            axis([1997 2012 b_min(iplt) b_max(iplt)])
            
            % calculate statistics
            db_CHAOS = b - mean(b) - dB_mod(:,iplt);
            index = ~isnan(dB_flow_3(:,iplt));
            db_flow_3 = b - mean(b) - interp1(t_3(index), dB_flow_3(index,iplt), t_tmp);
            
            text(.01, .87, ['\sigma = ' num2str(std(db_CHAOS), '%5.1f') ' nT/yr'], 'Fontsize', 10, 'Color', 'b', ...
                'VerticalAlignment', 'top', 'Units', 'normalized', 'Fontangle', 'italic')
            text(.01, .92, ['\sigma = ' num2str(nanstd(db_flow_3), '%5.1f') ' nT/yr'], 'Fontsize', 10, 'Color', 'g', ...
                'VerticalAlignment', 'top', 'Units', 'normalized', 'Fontangle', 'italic')
        end
    end
    set(gcf, 'PaperOrientation', 'Landscape', 'PaperType', 'A4', 'PaperUnits', ...
        'centimeters', 'PaperPosition', [0.2 0.2 [30 8]]);
    
    print('-dpsc','-append',filename_plot)
end

%% make movie of SV map in Br
clear all
c = 3480.0;
a = 6371.2;
rad = pi/180;

filename_plot=' ';
filename_out='dBr.png';

% load flow and SV predictions
load 'svpred_rmm_e-4_anis_ep.mat'

Step = 1;
theta = [180-Step/2:-Step:0];
phi   = [-180+Step/2:Step:180-Step/2];

C_max = 25000;
dC = C_max/5;
load coast;

pos_vector = [50 150 [700 600]];
figure('Position', pos_vector, 'Renderer','zbuffer');


% colourmap and bar
cptcmap('cbcPuOr', 'flip', true);
map = colormap;

writerObj = VideoWriter('dB_CMB.mp4', 'MPEG-4');
writerObj.FrameRate = 10;
open(writerObj);



 for i_plt = 2:186
% for i_plt = 2:10
B_r_s = synth_grid(gh_2(:,i_plt), c/a, theta, phi); %Core surface SV
B_ref =  makerefmat('RasterSize', size(B_r_s), ...
    'Latlim', [-90 90], 'Lonlim', [-180 180]);

clf
% North pole
pole_width = .3;
axes('Position', [0  0.45 pole_width pole_width], 'Box', 'off')
axesm('MapProjection', 'ortho', 'Origin', [90 0], 'FLatLimit',[-Inf 45.001])
meshm(B_r_s, B_ref)
hold on
plotm(lat, long, '-k', 'Linewidth', 1)
hold off
gridm
tightmap
caxis(C_max*[-1 1])
set(gca, 'Box', 'off', 'Visible', 'off')

% South pole
axes('Position', [1-pole_width 0.45 pole_width pole_width ], 'Box', 'off')
axesm('MapProjection', 'ortho', 'Origin', [-90 0], 'FLatLimit',[-Inf 45.001])
meshm(B_r_s, B_ref)
hold on
plotm(lat, long, '-k','Linewidth', 1)
hold off
gridm
tightmap
caxis(C_max*[-1 1])
set(gca, 'Box', 'off', 'Visible', 'off')

% Hammer
axes('Position', [0 0 1 .5], 'Box', 'off')
axesm('MapProjection', 'hammer', 'Frame', 'off')
meshm(B_r_s, B_ref)
hold on
plotm(lat, long, '-k','Linewidth', 1)
hold off
gridm
set(gca, 'Box', 'off', 'Visible', 'off')
caxis(C_max*[-1 1])
tightmap

% add text
title_text = num2str(t_2(i_plt), '%6.2f')
axes('Position', [0 0 1 .5], 'Box', 'off')
text(.05, .05, title_text, 'Interpreter', 'none')
set(gca, 'Box', 'off', 'Visible', 'off')

axes('Position', [0.333 0.6 .333 .015], 'Box', 'off')
image([1:length(map)], 'XData', [-C_max C_max]);
set(gca,'XTick', [-C_max:dC:C_max], 'XtickLabel', [-C_max/1000:dC/1000:C_max/1000], 'Xdir', 'Normal', 'FontWeight', 'normal', 'Xcolor', 'k', 'FontSize', 8, 'YTickLabel', {}, 'Ydir', 'Normal')
text(0.52, 2,' mT ', 'Color', 'k', 'FontWeight','normal', 'FontSize', 10, 'Units', 'normalized', 'HorizontalAlignment', 'Center')
 
f = getframe(gcf);              %Capture screen shot
    writeVideo(writerObj,f);
end

close(writerObj);


 


%%
clear all

% load flow and SV predictions
load 'svpred_rmm_e-4_anis_ep.mat'

clf
plot(t_1, t_coeff(1,:), 'r', t_1, t_coeff(9,:), 'b', ...
    t_1, zp_filter(1.138*(t_coeff(1,:) + 12/7*t_coeff(9,:)), 'boxcar', 7), 'g', 'Linewidth', 2)
grid on
axis([1997.5 2013 -16 2])
ylabel [km/yr]
legend('t_1^0','t_3^0', '1.138(t_1^0 + 12/7t_3^0)','Location','best')
legend('t_1^0','t_3^0', 'Location','best')

set(gcf, 'PaperOrientation', 'Portrait', 'PaperType', 'A4', 'PaperUnits', ...
    'centimeters', 'PaperPosition', [0.2 0.2 .6*[30 20]]);
filename_pdf = 't_10.eps';
print('-depsc', filename_pdf); system(['epstopdf ' filename_pdf]); delete(filename_pdf)




%%
cmap = [         0         0    0.5625
    0         0    0.6250
    0         0    0.6875
    0         0    0.7500
    0         0    0.8125
    0         0    0.8750
    0         0    0.9375
    0         0    1.0000
    0.0208         0    0.9792
    0.0417         0    0.9583
    0.0625         0    0.9375
    0.0833         0    0.9167
    0.1042         0    0.8958
    0.1250         0    0.8750
    0.1458         0    0.8542
    0.1667         0    0.8333
    0.1875         0    0.8125
    0.2083         0    0.7917
    0.2292         0    0.7708
    0.2500         0    0.7500
    0.2708         0    0.7292
    0.2917         0    0.7083
    0.3125         0    0.6875
    0.3333         0    0.6667
    0.3542         0    0.6458
    0.3750         0    0.6250
    0.3958         0    0.6042
    0.4167         0    0.5833
    0.4375         0    0.5625
    0.4583         0    0.5417
    0.4792         0    0.5208
    0.5000         0    0.5000
    0.5208         0    0.4792
    0.5417         0    0.4583
    0.5625         0    0.4375
    0.5833         0    0.4167
    0.6042         0    0.3958
    0.6250         0    0.3750
    0.6458         0    0.3542
    0.6667         0    0.3333
    0.6875         0    0.3125
    0.7083         0    0.2917
    0.7292         0    0.2708
    0.7500         0    0.2500
    0.7708         0    0.2292
    0.7917         0    0.2083
    0.8125         0    0.1875
    0.8333         0    0.1667
    0.8542         0    0.1458
    0.8750         0    0.1250
    0.8958         0    0.1042
    0.9167         0    0.0833
    0.9375         0    0.0625
    0.9583         0    0.0417
    0.9792         0    0.0208
    1.0000         0         0
    0.9375         0         0
    0.8750         0         0
    0.8125         0         0
    0.7500         0         0
    0.6875         0         0
    0.6250         0         0
    0.5625         0         0
    0.5000         0         0];

N = 14;
rho = 3485/6371;
% rho = 1;
clf
for it = 1:6:length(t_1)
    h5 = semilogy([1:N], powerspec(gh_1(1:N*(N+2),it), rho, 'int'), '-b', 'Linewidth', 2);
    hold on
    icolor = round(1+(it-1)*63/(length(t_1)-1));
    set(h5, 'Color', cmap(icolor, :))
end
hold off
xlabel 'degree \itn'
ylabel 'Secular variation power [(nT/yr)^2]'
grid on
axis([0 N 3e4 2e7])
% axis([0 N 4e-2 4e3])
cax = colorbar;
colormap(cmap);
set(cax, 'Ytick',[1:64/8:65], 'YtickLabel', [1997:2:2013])
set(gcf, 'PaperOrientation', 'Portrait', 'PaperType', 'A4', 'PaperUnits', ...
    'centimeters', 'PaperPosition', [0.2 0.2 .6*[30 20]]);
filename_pdf = 'Spectra_SV_CMB.eps';
print('-depsc', filename_pdf); system(['epstopdf ' filename_pdf]); delete(filename_pdf)





