%%  read subduction zone slabs from slab1.0 of Hayes et al 2012
% and calculate magnetic field produced by that slab model (forward modeling)
clear all
rad = pi/180;
a = 6371.2;

load coast

% load magnetic field model
load('/nio/m/mmod/CHAOS-7.4.mat')
g_C7 = g; g_C7(1:440) = ppval(pp, jd2000(2015, 1, 1)); % core fuield for epoch 2015

% load LCS-1
filename_in = '/nio/m/mmod/LCS-1.shc';
gh_L1 = read_shc(filename_in);
gh_L1(1:15*17,:) = 0;

% read dipole magnetisation <---> SH coefficient matrix
g = matfile('/nio/m/crust/G_crust_SHA_100km.mat'); % 70,000 dipoles
N_max = 160;
% g = matfile('/nio/m/crust/G_crust_SHA_70km.mat'); % 35,000 dipoles
% N_max = 185;

depth = a - g.r_s; % depth of dipoles
N_dipoles = length(g.theta_s); % number of dipoles
delta = a*sqrt(4*pi/N_dipoles); % mean separation of monopoles in [km]
chi = 0.04; % magnetic susceptibility [1]

d = dir('./slab1.0/*_slab1.0_clip.grd');

slab_i = zeros(size(g.phi_s));
for i = 1:length(d)
    filename_in = ['./slab1.0/' d(i).name];
    
    phi_slab = double(ncread(filename_in, 'x'));
    phi_slab = mod(phi_slab, 360);
    lat_slab = double(ncread(filename_in, 'y'));
    slab = -ncread(filename_in, 'z');
    
    slab(isnan(slab)) = 0;
    slab(slab > 69) = 0; % up to slab depth of 69 km
    slab(slab > 0) = 15*chi; % slab depth integrated susceptibility of kappa = 15 km x chi SI = 0.6 SI x km (if chi = 0.04)
    
    
    slab_i = slab_i + interp2(phi_slab, lat_slab, slab', mod(g.phi_s, 360), (90-g.theta_s), 'linear', 0);
    
end % loop over all 12 slab subduction zones


% write model of depth-integrated susceptibility at locations of dipoles
fid = fopen('h_slab1_0.txt','w');
fprintf(fid, '%6.2f %6.2f %10.4f\n', [g.phi_s 90-g.theta_s slab_i]');
fclose(fid);

% background core field at dipole locations from CHAOS-7
theta_s = g.theta_s; phi_s = g.phi_s;
B = synth_values(a-100, theta_s, phi_s, g_C7(1:195));
F = sqrt(B(:,1).^2 + B(:,2).^2 + B(:,3).^2); % ambient magnetic field strength [nT]

m_slab = (delta^2*F.*slab_i)/(4*pi); % dipole strengths in [nT*km^3]
m2h = 4*pi/(chi*delta^2); % factor transforming dipole strength to crustal thickness

gh_slab = g.G*m_slab; % determine Gauss coefficients of magnetic field produced by subduction model

factor = [];
for n = 1:N_max
    factor = [factor; repmat((a/(a+depth))^(n+2), 2*n+1, 1)];
end

gh_slab_surface = gh_slab./factor; % upward continuation from depth to surface

% magnetic field at regular grid at constant radius
Step = .25;
theta = [0:Step:180];
phi   = [-180+Step/2:Step:180];
theta(1) = 1e-8; theta(end) = 180-1e-8;  % to avoid singularity at poles

h = 0; % Earth' surface r = a;
B_r_L1 = synth_model(gh_L1, 1+h/6371, theta, phi); % LCS-1 model
B_r_slab = synth_model(gh_slab_surface, 1+h/6371, theta, phi); % from subduction model;

imagesc(phi, 90-theta, B_r_slab)
hold on
plot(long, lat, '-k')
hold off
set(gca, 'YDir', 'normal')
C_max = 300;
caxis(C_max*[-1 1])
blue_red_colormap
colorbar

% write in GMT format
grdwrite2(phi, 90-theta, B_r_L1, 'B_LCS1_0km.grd');
grdwrite2(phi, 90-theta, B_r_slab, 'B_slab1_0_0km.grd');

%% read and plot subduction zone slabs from slab1.0 of Hayes et al 2012
clear all
rad = pi/180;
a = 6371.2;

load coast

Step = .5;
theta = [0:Step:180]';
phi = [-180:Step:180]';
[phi_grid, theta_grid] = ndgrid(phi, theta);

d = dir('./slab1.0/*_slab1.0_clip.grd');

slab_i = zeros(size(phi_grid));
for i = 1:length(d)
    filename_in = ['./slab1.0/' d(i).name];
    
    phi_slab = double(ncread(filename_in, 'x'));
    % phi_slab = mod(phi_slab+180, 360) - 180;
    phi_slab = mod(phi_slab, 360);
    lat_slab = double(ncread(filename_in, 'y'));
    slab = -ncread(filename_in, 'z');
    
    slab(isnan(slab)) = 0;
    
    
    slab_i = slab_i + interp2(phi_slab, lat_slab, slab', mod(phi_grid, 360), (90-theta_grid), 'linear', 0);
    
end % loop over all 12 slab subduction zones

imagesc(phi, 90-theta, slab_i')
set(gca, 'YDir', 'normal')
hold on
plot(long, lat, '-k')
axis([-180 180 -90 90])
hold off
colorbar
caxis([0 200])


