function [B_magn_1, B_magn_2] = synth_values_MMA(r, theta, phi, t, s_MMA)
rad = pi/180;
a = 6371.2;

s_MMA.nmax_e = max(s_MMA.nm_qs_2(:,1)); k = 1; s_MMA.terms_e_2 = [];
for n = 1:s_MMA.nmax_e
    for m = 0:n
        for l = 1:length(s_MMA.nm_qs_2)
            if n == s_MMA.nm_qs_2(l,1) && m == s_MMA.nm_qs_2(l,2); s_MMA.terms_e_2 = [s_MMA.terms_e_2; k]; end
        end
        k = k+1;
        if m > 0
            for l = 1:length(s_MMA.nm_qs_2)
                if n == s_MMA.nm_qs_2(l,1) && m == -s_MMA.nm_qs_2(l,2); s_MMA.terms_e_2 = [s_MMA.terms_e_2; k]; end
            end
            k = k+1;
        end
    end
end

s_MMA.nmax_i = max(s_MMA.nm_gh_2(:,1)); k = 1; s_MMA.terms_i_2 = [];
for n = 1:s_MMA.nmax_i
    for m = 0:n
        for l = 1:length(s_MMA.nm_gh_2)
            if n == s_MMA.nm_gh_2(l,1) && m == s_MMA.nm_gh_2(l,2); s_MMA.terms_i_2 = [s_MMA.terms_i_2; k]; end
        end
        k = k+1;
        if m > 0
            for l = 1:length(s_MMA.nm_gh_2)
                if n == s_MMA.nm_gh_2(l,1) && m == -s_MMA.nm_gh_2(l,2); s_MMA.terms_i_2 = [s_MMA.terms_i_2; k]; end
            end
            k = k+1;
        end
    end
end

% calculate dipolar coordinates + matrix R
[theta_d, phi_d, R] = gg2gm(theta, phi, +1, [], [], [s_MMA.theta_NGP, s_MMA.phi_NGP]); % in [deg]

% calculate design matrices for internal sources
[A_r_i, A_theta_i, A_phi_i] = design_SHA(r/a, theta_d*rad, phi_d*rad, s_MMA.nmax_i, 'int');
[A_r_e, A_theta_e, A_phi_e] = design_SHA(r/a, theta_d*rad, phi_d*rad, s_MMA.nmax_e, 'ext');

qs_tmp = interp1(s_MMA.t_qs_2, s_MMA.qs_2, t, 'linear');
% add g10 term
qs_tmp = [interp1(s_MMA.t_qs_1, s_MMA.qs_1, t, 'linear') qs_tmp];
terms_e = [1; s_MMA.terms_e_2];
B_magn_1 = zeros(length(t), 3);
for k = 1:length(terms_e)
    B_magn_1(:,1) = B_magn_1(:,1) + A_r_e(:,terms_e(k)).*qs_tmp(:,k);
    B_magn_1(:,2) = B_magn_1(:,2) + A_theta_e(:,terms_e(k)).*qs_tmp(:,k);
    B_magn_1(:,3) = B_magn_1(:,3) + A_phi_e(:,terms_e(k)).*qs_tmp(:,k);
end
% rotate from dipole to geographic frame
B_tmp = B_magn_1(:,2);
B_magn_1(:,2) = B_tmp.*squeeze(R(:,1,1)) + B_magn_1(:,3).*squeeze(R(:,2,1));
B_magn_1(:,3) = B_tmp.*squeeze(R(:,1,2)) + B_magn_1(:,3).*squeeze(R(:,2,2));

% secondary field
gh_tmp = interp1(s_MMA.t_gh_2, s_MMA.gh_2, t, 'linear');
% add q10 term
gh_tmp = [interp1(s_MMA.t_gh_1, s_MMA.gh_1, t, 'linear') gh_tmp];
terms_i = [1; s_MMA.terms_i_2];
B_magn_2 = zeros(length(t), 3);
for k = 1:length(terms_i) %%% terms_i_2 => terms_i 
    B_magn_2(:,1) = B_magn_2(:,1) + A_r_i(:,terms_i(k)).*gh_tmp(:,k);
    B_magn_2(:,2) = B_magn_2(:,2) + A_theta_i(:,terms_i(k)).*gh_tmp(:,k);
    B_magn_2(:,3) = B_magn_2(:,3) + A_phi_i(:,terms_i(k)).*gh_tmp(:,k);
end
% rotate from dipole to geographic frame
B_tmp = B_magn_2(:,2);
B_magn_2(:,2) = B_tmp.*squeeze(R(:,1,1)) + B_magn_2(:,3).*squeeze(R(:,2,1));
B_magn_2(:,3) = B_tmp.*squeeze(R(:,1,2)) + B_magn_2(:,3).*squeeze(R(:,2,2));


% % primary field
% qs_tmp = interp1(s_MMA.t_qs_2, s_MMA.qs_2, t, 'linear');
% % add g10 term
% qs_tmp = [interp1(s_MMA.t_qs_1, s_MMA.qs_1, t, 'linear') qs_tmp];
% terms_e = [1; s_MMA.terms_e_2];
% B_magn_1 = zeros(length(t), 3);
% for k = 1:length(terms_e)
%     B_magn_1(:,1) = B_magn_1(:,1) + A_r_e(:,terms_e(k)).*qs_tmp(:,k);
%     B_magn_1(:,2) = B_magn_1(:,2) + A_theta_e(:,terms_e(k)).*qs_tmp(:,k);
%     B_magn_1(:,3) = B_magn_1(:,3) + A_phi_e(:,terms_e(k)).*qs_tmp(:,k);
% end
% % rotate from dipole to geographic frame
% B_tmp = B_magn_1(:,2);
% B_magn_1(:,2) = B_tmp.*squeeze(R(:,1,1)) + B_magn_1(:,3).*squeeze(R(:,2,1));
% B_magn_1(:,3) = B_tmp.*squeeze(R(:,1,2)) + B_magn_1(:,3).*squeeze(R(:,2,2));
% 
% % secondary field
% gh_tmp = interp1(s_MMA.t_gh_2, s_MMA.gh_2, t, 'linear');
% % add q10 term
% gh_tmp = [interp1(s_MMA.t_gh_1, s_MMA.gh_1, t, 'linear') gh_tmp];
% terms_i = [1; s_MMA.terms_i_2];
% B_magn_2 = zeros(length(t), 3);
% for k = 1:length(terms_i)
%     B_magn_2(:,1) = B_magn_2(:,1) + A_r_i(:,terms_i(k)).*gh_tmp(:,k);
%     B_magn_2(:,2) = B_magn_2(:,2) + A_theta_i(:,terms_i(k)).*gh_tmp(:,k);
%     B_magn_2(:,3) = B_magn_2(:,3) + A_phi_i(:,terms_i(k)).*gh_tmp(:,k);
% end
% % rotate from dipole to geographic frame
% B_tmp = B_magn_2(:,2);
% B_magn_2(:,2) = B_tmp.*squeeze(R(:,1,1)) + B_magn_2(:,3).*squeeze(R(:,2,1));
% B_magn_2(:,3) = B_tmp.*squeeze(R(:,1,2)) + B_magn_2(:,3).*squeeze(R(:,2,2));

end

