%% construct matrix G that transform from dipole magnetisation (aligned with ambient magnetic field)
% to SH coefficients of magnetic field produced by this magnetisation
clear all
a = 6371.2;
rad = pi/180;

load coast

% region of interest
phi_region = [130 220]'; theta_region = 90 - [65 30]'; % Aleutian

h = 150; % altitude of "target" magnetic field about ellipsoid;

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

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

% magnetic field at regular grid at constant radius
Step = .25;
theta = [180:-Step:0]';
theta(end) = 1e-8; theta(1) = 180-1e-8;  % to avoid singularity at poles
Step = .5;
phi   = [-180+Step/2:Step:180]';
[~, ~, Z_LCS1, alpha] = synth_grid_ellipsoid(gh_L1, theta, phi, h);

% transform from phi = -180 ... 180 to 0 ... 360
idx = phi <= 0;
phi = [phi(~idx); phi(idx)+360];
Z_LCS1 = [Z_LCS1(:,~idx) Z_LCS1(:,idx)];

[phi_grid, theta_grid] = ndgrid(phi, theta);

% remove values outside region of interest
idx = phi < phi_region(1) | phi > phi_region(2); 
phi(idx) = []; phi_grid(idx,:) = []; theta_grid(idx,:) = []; Z_LCS1(:,idx) = [];
idx = theta < theta_region(1) | theta > theta_region(2); 
theta(idx) = []; alpha(idx) = []; phi_grid(:,idx) = []; theta_grid(:,idx) = []; Z_LCS1(idx,:) = [];

cos_theta = cosd(theta_grid(:));
sin_theta = sind(theta_grid(:));
phi_grid = phi_grid(:);
theta_grid = theta_grid(:);

% locations of dipoles (for whole Earth)
N_dipoles = 50000;
points_s = eq_point_set_polar(2, N_dipoles);
phi_s = points_s(1,:)'/rad;                 % point source longitude [deg]
theta_s = points_s(2,:)'/rad;               % point source co-latitude [deg]
theta_s(1) = 1e-8; theta_s(end) = 180-1e-8; % avoid singularities near poles

delta = a*sqrt(4*pi/N_dipoles); % mean separation of point sources in [km]
area = delta^2; % surface area of each prism, in [km^2]

% delete dipoles outside of region of interest
idx = phi_s > phi_region(1) & phi_s < phi_region(2) & theta_s > theta_region(1) & theta_s < theta_region(2);
theta_s = theta_s(idx);
phi_s = phi_s(idx);
N_dipoles = length(phi_s);

r_s = geod2geoc(pi/2-theta_s*rad, 0); % radius of dipoles 

cos_theta_s = cosd(theta_s);
sin_theta_s = sind(theta_s);

r_grid = geod2geoc(pi/2-theta_grid*rad, h); % radius of target points (h km above ellipsoid)

% background core field up to n = 13
B = synth_values(r_s, theta_s, phi_s, g_C7(1:195));
B_r = B(:,1); B_theta = B(:,2); B_phi = B(:,3);

F = sqrt(B_r.^2 + B_theta.^2 + B_phi.^2); % ambient magnetic field strength [nT]
H = sqrt(B_theta.^2 + B_phi.^2);
D = atan2(B_phi, -B_theta)/rad; % ambient magnetic field declination [deg]
I = atan2(-B_r, H)/rad; % ambient magnetic field inclination [deg]

ss = -sind(I); % Directional cosine 1 % -sin(I), Fox-Maule PhD eq. (7.24)
cc = -cosd(I).*cosd(D); % Directional cosine 2 % -cos(I)*cos(D)
cs = +cosd(I).*sind(D); % Directional cosine 3 % +cos(I)*sin(D)

% determine the matrix G that transforms from magnetisation (in dipole domain) 
% to magnetic field (on regular grid)
G = zeros(length(theta_grid(:)), N_dipoles);
for j = 1:N_dipoles
    H_r = zeros(size(theta_grid(:)));
    if mod(j, 50) == 0; disp(['dipole # ', num2str(j), ' out of ', num2str(N_dipoles)]); end
    mu_ij = cos_theta*cos_theta_s(j)+sin_theta*sin_theta_s(j).*cosd(phi_grid-phi_s(j));
    r_ij_2 = r_grid.^2 + r_s(j).^2 - 2*r_grid.*r_s(j).*mu_ij; % (r_ij)^2 in km^2
    r_ij_3 = (r_ij_2.^(3/2)); % (r_ij)^3 in [km^3]
    A = r_grid - r_s(j).*mu_ij;
    B = r_s(j)*(sin_theta.*sin_theta_s(j)+cos_theta.*cos_theta_s(j).*cosd(phi_grid-phi_s(j)));
    A1 = r_s(j) - r_grid.*mu_ij;
    B1 = r_grid.*(cos_theta.*sin_theta_s(j)-sin_theta.*cos_theta_s(j).*cosd(phi_grid-phi_s(j)));
    C1 = -r_grid.*sin_theta.*sind(phi_grid-phi_s(j));
    H_11 = 3*A.*A1./r_ij_2 + mu_ij; % in [1]
    H_12 = 3*A.*B1./r_ij_2 - B1./r_grid; % in [1]
    H_13 = 3*A.*C1./r_ij_2 - C1./r_grid; % in [1]
    H_r = -(H_11.*ss(j) + H_12.*cc(j) + H_13.*cs(j))./r_ij_3; % in [km^-3]
    G(:,j) = reshape(H_r, size(theta_grid));
end

%%  read subduction zone slabs from slab1.0 of Hayes et al 2012
chi = 0.04; % magnetic susceptibility [1]

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

slab_i = zeros(size(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(phi_s, 360), (90-theta_s), 'linear', 0);
end % loop over all 12 slab subduction zones

clf
subplot(1,2,1)
idx = slab_i > 0;
slabs.phi_s = mod(phi_s(idx), 360);
slabs.theta_s = 90-theta_s(idx);
plot(mod(phi_s, 360), 90-theta_s, '.k', 'Markersize', 2)
hold on
plot(slabs.phi_s, slabs.theta_s, '.b', 'MarkerSize', 7)
plot(mod(long, 360), lat, '-k')
k = boundary(slabs.phi_s, slabs.theta_s, .8);
P_positions = [slabs.phi_s(k) slabs.theta_s(k)];
P = images.roi.Polygon(gca,'Position',P_positions);
hold off
axis([phi_region(1) phi_region(2) 90-theta_region(2) 90-theta_region(1)])

%% run this cell again and again, after editing polygon in left-panel plot
% find positions of nonz-eros dipoles in polygon
P_Positions = P.Position;
tf = inROI(P, phi_s, 90-theta_s);
slab_i = zeros(size(slab_i));
slab_i(tf) = 15*chi; % slab depth integrated susceptibility of kappa = 15 km x chi SI = 0.6 SI x;

m_slab = (delta^2*F.*slab_i)/(4*pi); % dipole strengths in [nT*km^3]

B_r = G*m_slab;
B_r = reshape(B_r, length(phi), [])';


clf
% subplot(1,2,1)
axes('Position', [.05 .07 .45 .85])
imagesc(phi, 90-theta, B_r)
hold on
plot(phi_s, 90-theta_s,'.k', 'Markersize', 2)
plot(mod(long, 360), lat, '-k')
set(gca, 'YDir', 'normal')
P = images.roi.Polygon(gca,'Position',P_Positions);
hold off
axis([phi_region(1) phi_region(2) 90-theta_region(2) 90-theta_region(1)])
C_max = 120;
caxis(C_max*[-1 1])
cptcmap('tmp_70km.cpt')
cb = cptcbar(gca, 'tmp_70km', 'eastoutside', true);
title 'B_r due to slab model'

% subplot(1,2,2)
axes('Position', [.55 .07 .45 .85])
imagesc(phi, 90-theta, -Z_LCS1)
hold on
[~, hc] = contour(phi, 90-theta, B_r, [-10 0 10], '-k');
hc.LineColor = [0 .6 0];
hc.LineWidth = 1;
plot(phi_s, 90-theta_s,'.k', 'Markersize', 2)
plot(mod(long, 360), lat, '-k')
hold off
axis([phi_region(1) phi_region(2) 90-theta_region(2) 90-theta_region(1)])
caxis(C_max*[-1 1])
set(gca, 'YDir', 'normal')
cptcmap('tmp_70km.cpt')
cb = cptcbar(gca, 'tmp_70km', 'eastoutside', true);
title 'B_r from LCS-1 (green: contourlines of B_r of slab model)'
