function [GtG, Gty, w_F, delta_d_F] = do_slave_scalar(t, r, theta, phi, Dst_e_scalar, Dst_i_scalar, F, ...
    options, model, params, m_i, N_coeff, N_Euler, N_para)

% $Rev: 6585 $ 
% $Date: 2014-03-01 23:02:33 +0000 (lø, 01 mar 2014) $

N_data_F = length(t);
B_r_scalar_synth = zeros(N_data_F, 1); B_theta_scalar_synth = zeros(N_data_F, 1); B_phi_scalar_synth = zeros(N_data_F, 1);
F_synth = zeros(N_data_F, 1); delta_d_F = zeros(N_data_F, 1);
w_F = zeros(N_data_F, 1); C_tilde = zeros(N_data_F, 1);

clear GtG Gty
GtG = zeros(N_para, N_para); Gty = zeros(N_para, 1);

for i=1:options.N_Data_Frac:N_data_F
    index = [i:min(i+options.N_Data_Frac-1, N_data_F)]';
    [A_r, A_theta, A_phi] = create_G_3(t(index), r(index), theta(index), phi(index), Dst_e_scalar(index), Dst_i_scalar(index), model, params);
    
    B_r_scalar_synth(index)      = A_r     * m_i(1:N_coeff);
    B_theta_scalar_synth(index)  = A_theta * m_i(1:N_coeff);
    B_phi_scalar_synth(index)    = A_phi   * m_i(1:N_coeff);
    F_synth(index)        = sqrt(B_r_scalar_synth(index).^2 + B_theta_scalar_synth(index).^2 + B_phi_scalar_synth(index).^2);
    delta_d_F(index)      = F(index) - F_synth(index);
    
    C_tilde(index) = params.sigma_ovh_2;
    
    residuals = 1./sqrt(C_tilde(index)).*delta_d_F(index);
    if strcmp(options.robust, 'Tuckey')
        w_F(index) = tuckey_w(residuals, params.c_Tuckey); % Tuckey weights
    elseif strcmp(options.robust, 'LS')
        w_F(index) = ones(size(index)); % unit weights        
    else 
        w_F(index) = min(1.5./abs(residuals),1); % Huber weights
    end
    
    if params.mu == 0
        C_y_1 = spdiags(w_F(index).*sin(theta(index))./C_tilde(index), 0, length(index), length(index)); % inverse of data covariance matrix
    else
        C_diag_tmp = C_tilde(index)./(w_F(index).*sin(theta(index))); % diagonal of data covariance matrix
        
        rad   = pi/180;
        theta_dip = 11.32*rad; phi_b = 289.59*rad;
        ST0   = sin(theta_dip); CT0 = cos(theta_dip);
        cos_theta_gg = cos(theta(index));
        sin_theta_gg = sin(theta(index));
        cos_theta_gm = CT0*cos_theta_gg + ST0*sin_theta_gg.*cos((phi(index) - phi_b));
        
        C_diag_RC = params.mu^2*((3*cos_theta_gm.^2-1).^2./(3*cos_theta_gm.^2+1));
        C_y_1 = spdiags(1./(C_diag_tmp+C_diag_RC), 0, length(index), length(index)); % inverse of data covariance matrix
    end
    
    G_tmp = repmat(B_r_scalar_synth(index)./F_synth(index), 1, N_coeff).*A_r + repmat(B_theta_scalar_synth(index)./F_synth(index), 1, N_coeff).*A_theta + ...
        repmat(B_phi_scalar_synth(index)./F_synth(index), 1, N_coeff).*A_phi;
    G_tmp = [G_tmp zeros(length(index), 3*sum(N_Euler))]; % add Euler angle parameters (not resolved by scalar data)
    GtG = GtG + G_tmp'*C_y_1*G_tmp;
    Gty = Gty + (C_y_1*G_tmp)'*delta_d_F(index);
end % summation over data segments of length not longer than options.N_Data_Frac


