function [B_r, varargout] = synth_model(m_i, r, theta, phi, varargin)
% [B_r, B_theta, B_phi] = synth_model(m_i, r, theta, phi)
% 
% Synthetic magnetic field values B_r, B_theta, B_phi (in nT) 
% on a regular grid (at the Earth's surface) given by 
% co-latitude theta(:) and longitude phi(:) (in degrees)
% r is the radius in units of the Earth's radius a.
% m_i(:) contains the Gauss coefficients (in nT, Schmidt normalized)
%
% Note: phi = (0...360), not (0...360-dphi)!
%
% September 2003, Nils Olsen, DSRI, Copenhagen

% Latest change: 

rad = pi/180;

% check input values
theta = theta(:); % convert to column vector
r = r(:)'; % convert to column vector

N =  sqrt(size(m_i, 1)+1)-1;

% if nargin > 4
%     source = varargin{1};
% else
%     source = 'int'; % internal sources by default
% end

% source = 'int';
r_n = [];
% if strcmp(source, 'int')
    for n=1:N
        r_n       = [r_n; repmat(r.^(-(n+2)), 2*n+1, 1)];
    end
% end

m_i = repmat(m_i, 1, length(r)).*r_n;

factor = []; for n= 1:N; factor = [factor; repmat((n+1), 2*n+1, 1)]; end
B_r = SH_synthesis(theta, phi(1:end-1), m_i.*repmat(factor, 1, length(r)));
B_r(:,end+1,:) = B_r(:,1,:);

if nargout > 1 % B_theta and B_phi are also calculated
    [B_theta, B_phi] = SH_synthesis(theta, phi(1:end-1), -m_i);
    B_theta(:,end+1,:) = B_theta(:,1,:);
    B_phi(:,end+1,:)   = B_phi(:,1,:);
    varargout{1} = B_theta; varargout{2} = B_phi;
end


