function [B_1, B_2] = synth_values_MIO(r,theta,phi,mjd2000,F107,s)
% B_1 = synth_values_MIO(r, theta, phi, t, F107, s);
% [B_1,B_2] = synth_values_MIO(r, theta, phi, t, F107, s);
%
% Calculate the primary (1) and induced (2) Sq magnetic field from a model
% in DIPOLE coordinates (parameterization of Sabaka et al., 2002).
%
% Inputs:   r(:)                    radius [km]
%           theta(:), phi(:)        geographic co-latitude, longitude [deg]
%           t(:)                    time MD2000 [days]
%           F107(:)                 F10.7 index
%           s is structure with elements
%             m_e_b(:,:)            primary Sq model for r < rho_Sq (REAL coefficients)
%             m_i(:,:)              secondary Sq model (REAL coefficients)
%             nmax, mmax            max degree and order in dipole coord.
%             p_vec(:)              diurnal wavenumbers
%             s_vec(:)              seasonal wavenumbers
%             theta_NGP             co-latitude of Northern dipole (in deg)
%             phi_NGP               longitude of Northern dipole (in deg)
%             h                     altitude of Sq current layer (in km)
%             N                     coefficient of the F10.7 dependence
%
% Dependendies :
% - "gg2gm" where geomagnetic coordinates are defined at epoch 2010
% - "getmut"
% - "design_SHA_Sq"
%
% A. Chulliat, 23/04/2011
% Nils Olsen (DTU Space), 17/10/2011: further optimisation and change of I/O
%       variables

a = 6371.2;
rho = r/a;

rho_Sq = (a+s.h)/a;

% return empty output if at least one of the first four inputs is empty
if (size(mjd2000) == 0 | size(r) == 0 | size(theta) == 0 | size(phi) == 0);
    B_1 = zeros(0,3);
    B_2 = zeros(0,3);
    return
end

max_size = max([size(rho); size(theta); size(phi); size(mjd2000)]);
max_length = max_size(1)*max_size(2);

% convert to matrix if input parameter is scalar
if length(rho)     == 1; rho = rho*ones(max_size); end;
if length(theta)   == 1; theta = theta*ones(max_size); end;
if length(phi)     == 1; phi = phi*ones(max_size); end;
if length(mjd2000) == 1; mjd2000 = mjd2000*ones(max_size); end;
if length(F107) == 1; F107 = F107*ones(max_size); end;
% check for equal length of all input parameters
if size(mjd2000) ~= size(r) | size(mjd2000) ~= size(theta) ...
        | size(mjd2000) ~= size(phi) | size(mjd2000) ~= size(F107);
    error('Variables must be of equal size (or scalars)');
end

mjd2000 = mjd2000(:); rho = rho(:); theta = theta(:); phi = phi(:);
F107 = F107(:);

t_mut = getmut(mjd2000, s.theta_NGP, s.phi_NGP);  % t in MJD2000

% calculate time in year (season)
t_1 = datevec(mjd2000);
year = t_1(:,1) + 2000;
ndays = zeros(size(year));
for i = 1:length(year)
    ndays(i) = sum(eomday(year(i),1:12),2);
end
t = (mjd2000 - md2000(year,1,1,0))./ndays;

% calculate number of coefficients and data
N_nm       = s.mmax*(s.mmax+2) + (s.nmax-s.mmax)*(2*s.mmax+1);
N_models   = length(s.p_vec)*length(s.s_vec);

if (size(s.m_e_b, 1) ~= N_nm)
    error('wrong number of model coefficients')
end

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

if isfield(s, 'N_Data_Frac') == 0;
    N_Data_Frac = 20; % process data in pieces of N_Data_Frac data points
else
    N_Data_Frac = s.N_Data_Frac;
end

N_data = length(rho);
B_1 = zeros(N_data, 3);
B_2 = zeros(N_data, 3);
for i=1:N_Data_Frac:N_data
    index = [i:min(i+N_Data_Frac-1, N_data)]';
    
    % calculate design matrices for internal sources
    [A_r_i_d,A_theta_i_d,A_phi_i_d] = design_SHA_Sq(...
        rho(index),theta_d(index),phi_d(index),t(index),t_mut(index),s.nmax,s.mmax,s.p_vec,s.s_vec, 'int');
    
    % F10.7 dependence
    w = (1 + s.N*F107(index));
    % induced fields
    if nargout > 1
        B_r_2_tmp     = A_r_i_d'*s.m_i(:);
        B_tmp         = A_theta_i_d'*s.m_i(:);
        B_phi_2_tmp   = A_phi_i_d'*s.m_i(:);
        B_theta_2_tmp = B_tmp.*squeeze(R(index,1,1)) + B_phi_2_tmp.*squeeze(R(index,2,1));
        B_phi_2_tmp   = B_tmp.*squeeze(R(index,1,2)) + B_phi_2_tmp.*squeeze(R(index,2,2));
        % correct for F10.7 dependence
        B_2(index,:)  = bsxfun(@times, [B_r_2_tmp B_theta_2_tmp B_phi_2_tmp], w);
    end
    
    % CASE #1: above Sq currents
    if (min(rho(index)) > rho_Sq)
        
        % (primary) coefficients above ionosphere
        f = [];
        for n = 1:s.nmax
            f = [f; repmat(-n/(n+1)*rho_Sq.^(2*n+1), min(2*n+1, 2*s.mmax+1), 1)];
        end
        m_e_a = bsxfun(@times, s.m_e_b, f);
        
        % calculate magnetic field and rotate from dipole to geographic frame
        B_r_1_tmp     = A_r_i_d'*m_e_a(:);
        B_tmp         = A_theta_i_d'*m_e_a(:);
        B_phi_1_tmp   = A_phi_i_d'*m_e_a(:);
        B_theta_1_tmp = B_tmp.*squeeze(R(index,1,1)) + B_phi_1_tmp.*squeeze(R(index,2,1));
        B_phi_1_tmp   = B_tmp.*squeeze(R(index,1,2)) + B_phi_1_tmp.*squeeze(R(index,2,2));
        
        % CASE #2: below Sq currents
    elseif (min(rho(index)) >=1 && max(rho(index)) < rho_Sq)
        
        % calculate design matrices for external sources
        [A_r_e_d,A_theta_e_d,A_phi_e_d] = design_SHA_Sq(...
            rho(index),theta_d(index),phi_d(index),t(index),t_mut(index),s.nmax,s.mmax,s.p_vec,s.s_vec, 'ext');
        
        % calculate magnetic field and rotate from dipole to geographic frame
        B_r_1_tmp     = A_r_e_d'*s.m_e_b(:);
        B_tmp         = A_theta_e_d'*s.m_e_b(:);
        B_phi_1_tmp   = A_phi_e_d'*s.m_e_b(:);
        B_theta_1_tmp = B_tmp.*squeeze(R(index,1,1)) + B_phi_1_tmp.*squeeze(R(index,2,1));
        B_phi_1_tmp   = B_tmp.*squeeze(R(index,1,2)) + B_phi_1_tmp.*squeeze(R(index,2,2));
        
        % CASE #3: error
    else
        error('wrong region (r <= a+h)')
    end
    
    % correct for F10.7 dependence
    B_1(index,:)     = bsxfun(@times, [B_r_1_tmp B_theta_1_tmp B_phi_1_tmp], w);
end
