
function gh = monopoles2gauss(model, m_i, N)
% Calculates Schmidt-normalized spherical harmonic expansion coefficients
% for the monopole solution m_i and maximum SH degree and order N

% m_i = m_i./model.dS;

cos_theta = cos(model.theta);
gh = zeros(N*(N+2), 1);
k = 0;
for n = 1:N
%    if mod(n, 10) == 0; disp(['n = ' num2str(n)]); end
    rho_n2 = (model.r).^(n+2); % (r/a)^(n+2)
    Pnm = legendre(n, cos_theta, 'sch')';   % Schmidt normalized Legendre functions
    k = k+1;       % terms corresponding to g_n^0
    gh(k) = sum(rho_n2.*m_i.*Pnm(:,1));
    for m = 1:n
        k = k+1;   % terms corresponding to g_n^m
        gh(k) = sum(rho_n2.*m_i.*Pnm(:,m+1).*cos(m.*model.phi));
        k = k+1;   % terms corresponding to h_n^m
        gh(k) = sum(rho_n2.*m_i.*Pnm(:,m+1).*sin(m.*model.phi));
    end
end

