function [B_x, B_y, B_z, varargout] = cm4(t, r, theta, phi, Dst, F107, varargin)
% [B_x, B_y, B_z]           = cm4(ut, r, theta, phi, Dst, F107)
% [B_x, B_y, B_z]           = cm4(ut, r, theta, phi, Dst, F107, options)
% [B_x, B_y, B_z, J_x, J_y] = cm4(ut, r, theta, phi, Dst, F107, options)
%
% Input:
%    t:         time in MJD2000 (days)
%    r:         radius (km)
%    theta:     geographic co-latitude (deg)
%    phi:       geographic longitude   (deg)
%    Dst:       Dst index (nT)
%    F107:      Solar flux F_10.7 (10^(-22)W/m^2/Hz)
% Output:
%    B_x(:,7)   B_x of 1) core,
%                      2) crustal,
%                      3) magnetospheric primary,
%                      4) magnetospheric secondary (induced),
%                      5) ionospheric primary,
%                      6) ionospheric secondary (induced),
%                      7) toroidal
%    J_x(:,4)   J_x of 1) Induced magnetospheric field [A/m]
%                      2) Primary ionospheric field
%                      3) Induced ionospheric field
%                      4) Poloidal field
%
% Nils Olsen, DTU Space, August 2010
%

if nargout < 5
    curr = 0;
else
    curr = 1; % calculate currents
end

a = 6371.2;
options.filename =  '/home/nio/m/tools/umdl.CM4';
if nargin == 6
    options.n_min = [ 1 14]; % n_min
    options.n_max = [13 15]; % n_max
    options.cord = 0;        % 0: geocentric coordinates
    options.pred = [1 1 1 1 1 1];
else
    options = varargin{1};
end

% determine size of input arrays
max_size = max([size(t); size(r); size(theta); size(phi); size(Dst); size(F107)]);
max_length = max_size(1)*max_size(2);
% convert to matrix if input parameter is scalar
if length(t)     == 1; t = t*ones(max_size); end;
if length(r)     == 1; r = r*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(Dst)   == 1; Dst = Dst*ones(max_size); end;
if length(F107)  == 1; F107 = F107*ones(max_size); end;
% check for equal length of all input parameters
if size(t) ~= size(r) | size(t) ~= size(theta) ...
        | size(t) ~= size(phi) | size(t) ~= size(Dst) | size(t) ~= size(F107);
    error('Variables must be of equal size (or scalars)');
    return
end
% convert to row vector
t = reshape(t, max_length, 1);
r = reshape(r, max_length, 1);
theta = reshape(theta, max_length, 1);
phi = reshape(phi, max_length, 1);
Dst = reshape(Dst, max_length, 1);
F107 = reshape(F107, max_length, 1);

% convert time from t=MJD2000 to UT in years
dof0 = [0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334]';
[iy, im, id, ih, min, sec] = datevec(t+730486); % t is MJD2000
ifleapyr = mod(iy, 4) == 0; % 1 if leap year, 0 otherwise
i2 = im > 2 & ifleapyr == 1;

% cm4 results in errorneous field values (also magnetosphere and Sq) after
% 2010.0, and therefore the time is shifted by 8 years if year >= 2010.0
if iy >= 2010
    ut = iy - 8 + (dof0(im) + id-1+i2 + (ih+min/60+sec/3600)/24)./(365+ifleapyr); % UT in yrs
else
    ut = iy + (dof0(im) + id-1+i2 + (ih+min/60+sec/3600)/24)./(365+ifleapyr); % UT in yrs
end
    
if exist(options.filename, 'file') > 0

    [B_x, B_y, B_z, J_x, J_y] = cm4_m(ut, r-a, theta, phi, Dst, F107, ...
        options.filename, options.n_max, options.n_min, options.cord, options.pred, curr);
    % convert output variables to matrices of same size as input variables
    B_x = squeeze(reshape(B_x, [max_size, 7]));
    B_y = squeeze(reshape(B_y, [max_size, 7]));
    B_z = squeeze(reshape(B_z, [max_size, 7]));

    % no proper core field determination after 2010.0 (since outside of
    % spline boundaries) and time has been shifted by 8 years
    if iy >= 2010 
        B_x(:,1) = NaN;
        B_y(:,1) = NaN;
        B_z(:,1) = NaN;
    end
    
    if curr == 1
        varargout{1} = squeeze(reshape(J_x, [max_size, 4]));
        varargout{2} = squeeze(reshape(J_y, [max_size, 4]));
    end
else
    error(['File ', options.filename, ' not found']);
end

