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

load coast.mat
load('LCS-1_model.mat')

gh_LCS_1 = read_shc('LCS-1.shc');

% map of strength of monopoles
clf
scatter(model.phi, 90-model.theta, 10, m, 'filled');
hold on
plot(long, lat, '-k')
hold off
caxis(.2*[-1 1])
colorbar
axis([-180 180 -90 90])
set(gca, 'Xtick', [-180:60:180], 'Ytick', [-90:30:90])
nio_colormap

%% calculate Gauss coefficients from monopoles and compare with LCS-1 Gauss coefficients
N = 185;

model_1.r = model.r/a;
model_1.theta = model.theta*rad;
model_1.phi = model.phi*rad;
gh_L1 = monopoles2gauss(model_1, m, N);

plot([1:N*(N+2)], gh_L1, [1:N*(N+2)], gh-gh_L1)

semilogy([1:185], powerspec(gh_LCS_1), [1:185], powerspec(gh_LCS_1- gh_L1))

%% calculate magnetic field Z from monopoles

Step = 1;
theta = [180:-Step:0];
phi   = [-180+Step/2:Step:180-Step/2];
theta(end) = 1e-8; theta(1) = 180-1e-8;  % to avoid singularity at poles

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

N_data = length(theta_grid(:));
N_monopoles = length(model.theta);

sin_theta = sind(theta_grid(:));
cos_theta = cosd(theta_grid(:));
sin_theta_model = sind(model.theta);
cos_theta_model = cosd(model.theta);
r_2 = r.^2;
r_model_2 = model.r.^2;

B_r = zeros(N_data, 1);
for i = 1:N_monopoles
    if mod(i, 1000) == 0; disp([i N_monopoles]); end
    cos_mu       = cos_theta.*cos_theta_model(i) + sin_theta.*sin_theta_model(i).*cosd(phi_grid(:)-model.phi(i));
    r_ik         = sqrt(r_model_2(i) + r_2 - 2.*r.*model.r(i).*cos_mu);
    rho_k_3      = (model.r(i) ./ r_ik).^3; 
    B_r     = B_r + (rho_k_3 .* (r./model.r(i) - cos_mu))*m(i);
end

Z = -reshape(B_r, length(theta), []);
imagesc(phi, 90-theta, Z)
colorbar
caxis(200*[-1 1])
hold on
plot(long, lat, '-k')
hold off
set(gca, 'YDir', 'normal')
nio_colormap

%% calculate magnetic field Z from SH expansion gh

Step = 1;
theta = [180:-Step:0];
phi   = [-180+Step/2:Step:180-Step/2];
theta(end) = 1e-8; theta(1) = 180-1e-8;  % to avoid singularity at poles

[X, Y, Z_SH] = synth_grid(gh, 1,  theta, phi);

imagesc(phi, 90-theta, Z_SH)
colorbar
caxis(200*[-1 1])
hold on
plot(long, lat, '-k')
hold off
set(gca, 'YDir', 'normal')
nio_colormap
