function [B_r, varargout] = synth_grid(m_i, r, theta, phi, varargin)
% B_r                      = synth_grid(m_i, r, theta, phi)
% [B_r, B_theta]           = synth_grid(m_i, r, theta, phi)
% [B_r, B_theta, B_phi]    = synth_grid(m_i, r, theta, phi)
% [B_r, B_theta, B_phi, F] = synth_grid(m_i, r, theta, phi)
% [B_r, B_theta, B_phi, F] = synth_grid(m_i, r, theta, phi, source)
% 
% Synthetic magnetic field values B_r, B_theta, B_phi, F (in nT) 
% on a regular grid given by radius r (in units of the Earth's radius),
% co-latitude theta(:) and longitude phi(:) (in degrees)
% m_i(:) contains the Gauss coefficients (in nT, Schmidt normalized)
% source: optional 5th argrument specifies if coefficients are for
% internal ('int', default), external ('ext') or toroidal ('tor') sources
%
% July 2003, Nils Olsen, DSRI

% Last change: NIO 080703: only the specified field components are calculated

theta = theta(:)';
phi = phi(:)';

rad = pi/180;
maxdeg_MF = sqrt(length(m_i)+1)-1;
cos_theta = cos(theta*rad);
sin_theta = sin(theta*rad);
N_theta   = length(theta);
N_phi     = length(phi);
n = [1:maxdeg_MF];

if nargin > 4
    source = varargin{1};
else
    source = 'int'; % internal sources by default
end
if strcmp(source, 'int')   
    r_n       = r.^(-(n+2));
    fac_r     = n+1;
elseif strcmp(source, 'ext')   
    r_n       = r.^(+(n-1));
    fac_r     = -n;
elseif strcmp(source, 'tor')   
    r_n       = repmat(r.^(1), size(n));
    fac_r     = -n.*(n+1);
else
    warning 'Unrecognized source: not ''int'', ''ext'', or ''tor'''
end

n_koeff = maxdeg_MF*(maxdeg_MF+2);
T_r     = zeros(N_theta, n_koeff);
cos_sin_m = zeros(N_phi, n_koeff);
if nargout > 1; % also calculation of B_theta
    T_theta = zeros(N_theta, n_koeff);
end;
if nargout > 2; % also calculation of B_phi
    T_phi   = zeros(N_theta, n_koeff);
    sin_cos_m = zeros(N_phi, n_koeff);
end

k=1;
for n = 1:maxdeg_MF  
    Pnm = legendre(n, cos_theta, 'sch')';
    
    if nargout > 1; % also calculation of B_theta
        dPnm(:,n+1) =  (sqrt(n/2).*Pnm(:,n));      % m=n
        dPnm(:,1) = -sqrt(n*(n+1)/2.).*Pnm(:,2);   % m=0
        if n > 1; dPnm(:,2)=(sqrt(2*(n+1)*n).*Pnm(:,1)-sqrt((n+2)*(n-1)).*Pnm(:,3))/2; end; % m=1
        for m = 2:n-1                              % m=2...n-1
            dPnm(:,m+1)=(sqrt((n+m)*(n-m+1)).*Pnm(:,m)-sqrt((n+m+1)*(n-m)).*Pnm(:,m+2))/2;
        end;
        if n == 1 dPnm(:,2) = sqrt(2)*dPnm(:,2); end
    end
    
    T_r(:,k)     = fac_r(n).*r_n(n).*Pnm(:,1);
    cos_sin_m(:,k) = ones(N_phi, 1);

    if nargout > 1; % also calculation of B_theta
        T_theta(:,k) = -r_n(n).*dPnm(:,1);
    end
    if nargout > 2; % also calculation of B_phi
        T_phi(:,k)   = zeros(N_theta, 1);
        sin_cos_m(:,k) = zeros(N_phi, 1);
    end
    k=k+1;
    for m = 1:n
        T_r(:,k:k+1)     =  fac_r(n).*r_n(n).*[Pnm(:,m+1) Pnm(:,m+1)];
        cos_sin_m(:,k:k+1) = [cos(m*phi'*rad)  sin(m*phi'*rad)];
        if nargout > 1; % also calculation of B_theta
            T_theta(:,k:k+1) = -r_n(n).*[dPnm(:,m+1) dPnm(:,m+1)];
        end
        if nargout > 2; % also calculation of B_phi
            T_phi(:,k:k+1)   = m.*r_n(n).*[Pnm(:,m+1)./sin_theta' Pnm(:,m+1)./sin_theta'];
            sin_cos_m(:,k:k+1) = [sin(m*phi'*rad) -cos(m*phi'*rad)];
        end
        k = k+2;
    end
end

m = spdiags(m_i, 0, n_koeff, n_koeff);
B_r = T_r*(cos_sin_m*m)';

if nargout > 1; % also calculation of B_theta
    varargout{1} = T_theta*(cos_sin_m*m)';
end;

if nargout > 2; % also calculation of B_phi
    varargout{2} = T_phi*(sin_cos_m*m)';
end;

if nargout == 4;  % also calculation of F
    varargout{3} = sqrt(B_r.^2 + varargout{1}.^2 + varargout{2}.^2); 
end;
