%% modmclb.m is the equatorial magnetic field model from McIlwain 1986.
%%

function bf=modmclb(mlt,lsh)

% mlt = magnetic local time (fractional hours)
% lsh = L-shell
% bf = magnetic field strength (T)

%constant values
c1=6;
c2=-24;
c3=18;
c4=1728;
c5=31000;

mlt2=(pi/12)*mlt;             %MLT in rads

% B field strength
bf1=c1;
bf2=c2*cos(mlt2);
bf3=(c3*cos(mlt2).^2)./(1+c4./lsh.^3);
bf4=c5./lsh.^3;
bf=bf1+bf2+bf3+bf4;
bf=1e-9*bf;