function [A_r, A_theta, A_phi] = design_SHA_Sq(rho,theta,phi,t,t_ut,...
    nmax,mmax,p_vec,s_vec, varargin)

% [A_r, A_theta, A_phi] = design_SHA_Sq(rho,theta,phi,t,t_ut,...
%    nmax,mmax,p_vec,s_vec, i_e_flag);
%
% Calculate design matrices A_i that connects the vector of REAL,
% Schmidt-normalized, spherical harmonic expansion coefficients,
%        x = g_{nsp}^{m} and h_{nsp}^{m}
% and the magnetic component B_r, B_theta, B_phi:
%        B_r     = x'*A_r
%        B_theta = x'*A_theta
%        B_phi   = x'*A_phi
% for the region r > a+h above (if i_e_flag = 'int') or below 
%             (if i_e_flag = 'int') the Sq currents (Sabaka et al., 2004).
%
% Inputs:   rho(:)                  radius [units of reference radius]
%           theta(:), phi(:)        co-latitude, longitude [deg]
%           t(:), t_ut(:)           time [yr]; UT [h]
%           nmax, mmax              maximum degree and order
%           p_vec(:)                diurnal wavenumbers
%           s_vec(:)                seasonal wavenumbers
%
% (Optimized version)
%
% A. Chulliat, 21/04/2011
% Nils Olsen (DTU Space), 17/10/2011: further optimisation
% Nils Olsen (DTU Space), 10/11/2011: change to new MIO format 
%                               (primary field is only given for r < r_Sq)

rad = pi/180;

w_s = 2*pi;
w_p = 2*pi/24;

N_data = length(theta);
if (length(phi) ~= N_data)
    error('design_SHA_Sq: theta and phi have different dimensions')
end

if nargin == 9;
    i_e_flag = 'int'; % internal sources by default
elseif nargin > 9
    i_e_flag = varargin{1};
else
    error('At least 9 inputs required')
end;


cos_theta = cosd(theta(:));
sin_theta = sind(theta(:));

% convert to row vector if input parameter is scalar

if length(rho) == 1; rho = rho*ones(N_data,1); end;
if length(t) == 1; t = t*ones(N_data,1); end;
if length(t_ut) == 1; t_ut = t_ut*ones(N_data,1); end;

% number of parameters

N_nm = mmax*(mmax+2) + (nmax-mmax)*(2*mmax+1);
N_sp    = length(p_vec)*length(s_vec);
N_coeff = 2*N_nm*N_sp;

% calculate sub-matrices for s=0, p=0

A_r_0     = zeros(N_nm, N_data);
A_theta_0 = zeros(N_nm, N_data);
A_phi_0   = zeros(N_nm, N_data);

k = 0;

if ~strcmp(i_e_flag, 'ext') % internal sources by default
    for n = 1:nmax
        rn1       = rho.^(-(n+2));
        
        Pnm = legendre(n, cos_theta, 'sch')';
        
        dPnm = zeros(size(Pnm,1),n+1);
        dPnm(:,n+1) = (sqrt(n/2).*Pnm(:,n));        % m=n
        dPnm(:,1) = -sqrt(n*(n+1)/2.).*Pnm(:,2);    % m=0
        if n > 1                                    % m=1
            dPnm(:,2) = (sqrt(2*(n+1)*n).*Pnm(:,1) - ...
                sqrt((n+2)*(n-1)).*Pnm(:,3))/2;
        end
        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
        
        for m = 0:min(n,mmax)
            
            alpha = m*phi*rad;
            cos_phi = cos(alpha);
            sin_phi = sin(alpha);
            
            if m==0
                k = k+1;
                A_r_0(k,:)     = (n+1).*rn1(:).*Pnm(:,1);
                A_theta_0(k,:) = -rn1(:).*dPnm(:,1);
                A_phi_0(k,:)   = 0;
            else
                k = k+1;
                A_r_0(k,:)     = (n+1).*rn1(:).*Pnm(:,m+1).*cos_phi(:);
                A_theta_0(k,:) = -rn1(:).*dPnm(:,m+1).*cos_phi(:);
                A_phi_0(k,:)   =  rn1(:).*Pnm(:,m+1)./sin_theta(:).*m.*sin_phi;
                
                k = k+1;
                A_r_0(k,:)     = (n+1).*rn1(:).*Pnm(:,m+1).*sin_phi(:);
                A_theta_0(k,:) = -rn1(:).*dPnm(:,m+1).*sin_phi(:);
                A_phi_0(k,:)   = -rn1(:).*Pnm(:,m+1)./sin_theta(:).*m.*cos_phi;
            end
            
        end
    end
else % external sources
    for n = 1:nmax
        rn1       = rho.^(n-1);
        
        Pnm = legendre(n, cos_theta, 'sch')';
        
        dPnm = zeros(size(Pnm,1),n+1);
        dPnm(:,n+1) = (sqrt(n/2).*Pnm(:,n));        % m=n
        dPnm(:,1) = -sqrt(n*(n+1)/2.).*Pnm(:,2);    % m=0
        if n > 1                                    % m=1
            dPnm(:,2) = (sqrt(2*(n+1)*n).*Pnm(:,1) - ...
                sqrt((n+2)*(n-1)).*Pnm(:,3))/2;
        end
        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
        
        for m = 0:min(n,mmax)
            
            alpha = m*phi*rad;
            cos_phi = cos(alpha);
            sin_phi = sin(alpha);
            
            if m==0
                k = k+1;
                A_r_0(k,:)     = -n.*rn1(:).*Pnm(:,1);
                A_theta_0(k,:) = -rn1(:).*dPnm(:,1);
                A_phi_0(k,:)   = 0;
            else
                k = k+1;
                A_r_0(k,:)     = -n.*rn1(:).*Pnm(:,m+1).*cos_phi(:);
                A_theta_0(k,:) = -rn1(:).*dPnm(:,m+1).*cos_phi(:);
                A_phi_0(k,:)   =  rn1(:).*Pnm(:,m+1)./sin_theta(:).*m.*sin_phi;
                
                k = k+1;
                A_r_0(k,:)     = -n.*rn1(:).*Pnm(:,m+1).*sin_phi(:);
                A_theta_0(k,:) = -rn1(:).*dPnm(:,m+1).*sin_phi(:);
                A_phi_0(k,:)   = -rn1(:).*Pnm(:,m+1)./sin_theta(:).*m.*cos_phi;
            end
            
        end
    end
end

% calculate full matrices

A_r     = zeros(N_coeff, N_data);
A_theta = zeros(N_coeff, N_data);
A_phi   = zeros(N_coeff, N_data);

i_sp = 0;

for s = s_vec;
    for p = p_vec;
        
        k_min_c = N_nm*(2*i_sp)+1;
        k_max_c = N_nm*(2*i_sp+1);
        k_min_s = N_nm*(2*i_sp+1)+1;
        k_max_s = N_nm*(2*i_sp+2);
        
        if (s==0 && p==0)
            A_r(k_min_c:k_max_c,:)     = A_r_0;
            A_theta(k_min_c:k_max_c,:) = A_theta_0;
            A_phi(k_min_c:k_max_c,:)   = A_phi_0;
        else
            beta = w_s*s*t + w_p*p*t_ut;
            cos_beta = cos(beta(:))';
            sin_beta = sin(beta(:))';
            A_r(k_min_c:k_max_c,:)     = bsxfun(@times, A_r_0, cos_beta);
            A_theta(k_min_c:k_max_c,:) = bsxfun(@times, A_theta_0, cos_beta);
            A_phi(k_min_c:k_max_c,:)   = bsxfun(@times, A_phi_0, cos_beta);
            
            A_r(k_min_s:k_max_s,:)     = bsxfun(@times, A_r_0, sin_beta);
            A_theta(k_min_s:k_max_s,:) = bsxfun(@times, A_theta_0, sin_beta);
            A_phi(k_min_s:k_max_s,:)   = bsxfun(@times, A_phi_0, sin_beta);
        end
        
        i_sp = i_sp+1;
        
    end
end
