designmatrixn.m

function A = designmatrixn(longitude,latitude,height,n_max,GM,R)
%
% Calculation of normalized Legendre Functions using stable recursion
% formulas. 
%
% --Input--
%   longitude     Longitudes (n x 1) in Radian
%   latitudes     Latitudes (n x 1)  in Radian
%   height        Heights (n x 1)
%   n_max         Maximum degree (1 x 1)
%   GM            Gravitational constant times mass of the Earth
%   R             Radius of the Earth
%
% --Output--
%   A             Design matrix (n x n_max)
%
% Author: Anne Springer
% Date:  2015-08-24
%

% Degree and order
idx = tril(ones(n_max+1))~=0;
degree = repmat((0:n_max)',1,n_max+1);
degree = degree(idx);
order = repmat(0:n_max,n_max+1,1);
order = order(idx);

% Design matrix
A = zeros(length(longitude),length(degree)*2-(n_max+1));
for i = 1:length(A)
    Pnm = calculatePnm((pi/2-latitude(i)), n_max);
    Pnm = Pnm(idx);
    factor = GM/R*(R/height(i)).^(degree+1);
    cos_ = cos(order*longitude(i));
    sin_ = sin(order*longitude(i));
    A(i,:) = [factor.*cos_.*Pnm; factor(n_max+2:end).*sin_(n_max+2:end).*Pnm(n_max+2:end)]';
end