function [A_r, A_theta, A_phi] = create_G_3(t, r, theta, phi, Dst_e, Dst_i, model, params)

% $Id: create_G_3.m 7614 2015-07-30 23:28:40Z nio $

rad = pi/180;

N_data = length(t);   % number of data points
N_t_segments_q10 = length(model.t_break_q10)-1; % breaks for q10
N_t_segments_qs11 = length(model.t_break_qs11)-1; % breaks for q11, s11
N_coeff_MF = model.maxdeg_MF*(model.maxdeg_MF+2); % number of static coefficients
K = model.maxdeg_spline*(model.maxdeg_spline+2);  % number of SHA coefficients per knot

% static (core and crustal) field, and secular variation
[A_r_i, A_theta_i, A_phi_i] = design_SHA(r, theta, phi, model.maxdeg_MF, 'int');

coll = spcol(model.knots, model.sp_order, [model.knots(1); t], 'sp', 'no');
coll = coll(2:end, :); % addition/removal of knots(1) is necesarry to get correct coll
L = size(coll, 2);    % number of knots
N_coeff_spline = K*L;

index_sort = []; for l=1:K; index_sort = [index_sort; [l:K:N_coeff_spline]']; end;
% the following two lines contain an accelerated version of
% G = []; for l=1:L; G = [G spdiags(repmat(coll(index, l)', l, K)', 0, N_data, N_data)*A_r_i(index,:)]; end
A_r = zeros(N_data, N_coeff_spline); A_theta = zeros(N_data, N_coeff_spline); A_phi = zeros(N_data, N_coeff_spline);
for l=1:L; tmp = spdiags(repmat(coll(:, l)', l, K)', 0, N_data, N_data)*A_r_i(:, 1:K);     A_r(:, (l-1)*K+1:l*K) = tmp;     end
for l=1:L; tmp = spdiags(repmat(coll(:, l)', l, K)', 0, N_data, N_data)*A_theta_i(:, 1:K); A_theta(:, (l-1)*K+1:l*K) = tmp; end
for l=1:L; tmp = spdiags(repmat(coll(:, l)', l, K)', 0, N_data, N_data)*A_phi_i(:, 1:K);   A_phi(:, (l-1)*K+1:l*K) = tmp;   end
A_r = A_r(:, index_sort); A_theta = A_theta(:, index_sort); A_phi = A_phi(:, index_sort);

% add static coefficients above spline truncation level
A_r =     [A_r     A_r_i(:,     K+1:N_coeff_MF)];
A_theta = [A_theta A_theta_i(:, K+1:N_coeff_MF)];
A_phi =   [A_phi   A_phi_i(:,   K+1:N_coeff_MF)];

% add external fields
if model.maxdeg_ext > 0 % q_n^m, q_{n,Dst}^m (resp. s_n^m and h_n^m)
    
%   [theta_sm, phi_sm, R_1]   = geo2sm(t, theta/rad, phi/rad);  %old version
    [theta_sm, phi_sm, R_1]   = geo2sm(t, theta/rad, phi/rad, params.igrf_dipole);  
    [A_r_ext, A_theta_sm_ext, A_phi_sm_ext] = design_SHA(r, theta_sm*rad, phi_sm*rad, model.maxdeg_ext, 'ext');
    [A_r_ind, A_theta_sm_ind, A_phi_sm_ind] = design_SHA(r, theta_sm*rad, phi_sm*rad, model.maxdeg_ext, 'int');
    N_coeff_sm = size(A_r_ext, 2);
    A_theta_ext = repmat(R_1(:,1,1), 1, N_coeff_sm).*A_theta_sm_ext + repmat(R_1(:,2,1), 1, N_coeff_sm).*A_phi_sm_ext;
    A_phi_ext =   repmat(R_1(:,1,2), 1, N_coeff_sm).*A_theta_sm_ext + repmat(R_1(:,2,2), 1, N_coeff_sm).*A_phi_sm_ext;
    A_theta_ind = repmat(R_1(:,1,1), 1, N_coeff_sm).*A_theta_sm_ind + repmat(R_1(:,2,1), 1, N_coeff_sm).*A_phi_sm_ind;
    A_phi_ind =   repmat(R_1(:,1,2), 1, N_coeff_sm).*A_theta_sm_ind + repmat(R_1(:,2,2), 1, N_coeff_sm).*A_phi_sm_ind;   
    
    [A_r_sm, A_theta_sm, A_phi_sm]     = design_SHA_sm(t, r, theta, phi, [1 1 1 1 1 1 1 1], params.filename_GSM_SM_SHA, params.SM_i_e_flag);
    [A_r_gsm, A_theta_gsm, A_phi_gsm]  = design_SHA_gsm(t, r, theta, phi, [1 1 1 1], params.filename_GSM_SM_SHA, params.GSM_i_e_flag);
      
    A_r_sm_gsm = [A_r_sm A_r_gsm(:, [1 4])];
    A_theta_sm_gsm = [A_theta_sm A_theta_gsm(:, [1 4])];
    A_phi_sm_gsm = [A_phi_sm A_phi_gsm(:, [1 4])];
    
    % add the terms that change explicitely in time (q10, ...)
    A_break_1 = zeros(N_data, N_t_segments_q10);
    [~, whichBin] = histc(t, model.t_break_q10);
    for it = 1:N_t_segments_q10
        A_break_1(whichBin == it, it) = 1;
    end
    A_break_2 = zeros(N_data, N_t_segments_qs11);
    [~, whichBin] = histc(t, model.t_break_qs11);
    for it = 1:N_t_segments_qs11
        A_break_2(whichBin == it, it) = 1;
    end
    
    % q10, q11, s11
    A_phi_q10   = [repmat(A_phi_sm_gsm(:,1),   1, N_t_segments_q10).*A_break_1 ...
        repmat(A_phi_sm_gsm(:,2),   1, N_t_segments_qs11).*A_break_2 ...
        repmat(A_phi_sm_gsm(:,3),   1, N_t_segments_qs11).*A_break_2];
    A_theta_q10 = [repmat(A_theta_sm_gsm(:,1), 1, N_t_segments_q10).*A_break_1 ...
        repmat(A_theta_sm_gsm(:,2), 1, N_t_segments_qs11).*A_break_2 ...
        repmat(A_theta_sm_gsm(:,3), 1, N_t_segments_qs11).*A_break_2];
    A_r_q10     = [repmat(A_r_sm_gsm(:,1),     1, N_t_segments_q10).*A_break_1 ...
        repmat(A_r_sm_gsm(:,2),     1, N_t_segments_qs11).*A_break_2 ...
        repmat(A_r_sm_gsm(:,3),     1, N_t_segments_qs11).*A_break_2];
     
    A_r = [A_r A_r_q10 A_r_sm_gsm(:,4:end) ...
        Dst_e.*A_r_ext(:,1) + Dst_i.*A_r_ind(:,1) ...
        Dst_e.*A_r_ext(:,2) + Dst_i.*A_r_ind(:,2) ...
        Dst_e.*A_r_ext(:,3) + Dst_i.*A_r_ind(:,3) ...
        ];
    A_theta = [A_theta A_theta_q10 A_theta_sm_gsm(:,4:end) ...
        Dst_e.*A_theta_ext(:,1) + Dst_i.*A_theta_ind(:,1) ...
        Dst_e.*A_theta_ext(:,2) + Dst_i.*A_theta_ind(:,2) ...
        Dst_e.*A_theta_ext(:,3) + Dst_i.*A_theta_ind(:,3) ...
        ];
    A_phi = [A_phi A_phi_q10 A_phi_sm_gsm(:,4:end) ...
        Dst_e.*A_phi_ext(:,1) + Dst_i.*A_phi_ind(:,1) ...
        Dst_e.*A_phi_ext(:,2) + Dst_i.*A_phi_ind(:,2) ...
        Dst_e.*A_phi_ext(:,3) + Dst_i.*A_phi_ind(:,3) ...
        ];
end
