function [Bx, By, Bz] = igrf(time, latitude, longitude, altitude, coord)
% IGRF Earth's magnetic field from IGRF model.
%
% Usage: [BX, BY, BZ] = IGRF(TIME, LATITUDE, LONGITUDE, ALTITUDE, COORD)
% or [BX, BY, BZ] = IGRF(COEFS, LATITUDE, LONGITUDE, ALTITUDE, COORD)
% or B = IGRF(TIME, LATITUDE, LONGITUDE, ALTITUDE, COORD)
% or B = IGRF(COEFS, LATITUDE, LONGITUDE, ALTITUDE, COORD)
%
% Calculates the components of the Earth's magnetic field using the
% International Geomagnetic Reference Field (IGRF) model. The inputs for
% the position can be scalars or vectors (in the latter case each should
% have the same number of elements or be a scalar), but TIME must be a
% scalar.
%
% When all the coordinate inputs are scalars, the function can be run more
% efficiently by providing the proper IGRF coefficient vector for a given
% time rather than the time itself. This mode is useful when making
% multiple calls to the function while keeping the time the same (meaning
% the coefficients will be the same for each run) as loading the
% coefficients can be the most time-consuming part of the function. The
% coefficient vector can be easily loaded using the function LOADIGRFCOEFS.
% This mode is assumed when all the coordinate inputs are scalars and the
% first input is a vector. In this case, the coefficient vector should be
% formatted as (LOADIGRFCOEFS provides this):
%
% [g(n=1,m=0) g(n=1,m=1) h(n=1,m=1) g(n=2,m=0) g(n=2,m=1) h(n=2,m=1) ...]
%
% Regardless of the size of the inputs, the outputs will be column vectors.
% If only one output is requested, B = [BX(:), BY(:), BZ(:)] is output.
% Note that the other parameters the IGRF gives can be computed from BX,
% BY, and BZ as:
%
% Horizonal intensity: hypot(BX, BY) (i.e., sqrt(BX.^2 + BY.^2) )
% Total intensity: hypot(BX, hypot(BY, BZ))
% Declination: atan2(BY, BX)
% Inclination: atan(BZ./hypot(BX, BY))
%
% This function relies on having the file igrfcoefs.mat in the MATLAB
% path to function properly when a time is input. If this file cannot be
% found, this function will try to create it by calling GETIGRFCOEFS.
%
% The IGRF is a spherical harmonic expansion of the Earth's internal
% magnetic field. Currently, the IGRF model is valid between the years 1900
% and 2015. See the health warning for the IGRF model here:
% http://www.ngdc.noaa.gov/IAGA/vmod/igrfhw.html
%
% Reference:
% International Association of Geomagnetism and Aeronomy, Working Group
% V-MOD (2010), International Geomagnetic Reference Field: the eleventh
% generation, _Geophys. J. Int._, _183_(3), 1216-1230,
% doi:10.1111/j.1365-246X.2010.04804.x.
%
% Inputs:
% -TIME: Time to get the magnetic field values either in MATLAB serial
% date number format or a string that can be converted into MATLAB serial
% date number format using DATENUM with no format specified (see
% documentation of DATENUM for more information).
% -COEFS: Instead of inputting a time, you can simply specify the proper
% coefficients for the time you want by inputting in the first argument
% the proper coefficient vector from igrfcoefs.mat.
% -LATITUDE: Geocentric or geodetic latitude in degrees.
% -LONGITUDE: Geocentric or geodetic longitude in degrees.
% -ALTITUDE: For geodetic coordiates, the height in km above the Earth's
% surface. For geocentric coordiates, the radius in km from the center of
% the Earth.
% -COORD: String specifying the coordinate system to use. Either
% 'geocentric' or 'geodetic' (optional, default is geodetic). Note that
% only geodetic coordinates have been verified.
%
% Outputs:
% -BX: Northward component of the magnetic field in nanoteslas (nT).
% -BY: Eastward component of the magnetic field in nT.
% -BZ: Downward component of the magnetic field in nT.
% -B: [BX(:), BY(:), BZ(:)].
%
% See also: LOADIGRFCOEFS, GETIGRFCOEFS, IGRFLINE, DATENUM, IGRF11MAGM.
% Run IGRFS if all position inputs are scalars.
if isscalar(latitude) && isscalar(longitude) && isscalar(altitude)
if nargin < 5
[Bx, By, Bz] = igrfs(time, latitude, longitude, altitude);
else
[Bx, By, Bz] = igrfs(time, latitude, longitude, altitude, coord);
end
% Otherwise run IGRFV.
else
if nargin < 5
[Bx, By, Bz] = igrfv(time, latitude, longitude, altitude);
else
[Bx, By, Bz] = igrfv(time, latitude, longitude, altitude, coord);
end
end
if nargout <= 1
Bx = [Bx(:), By(:), Bz(:)];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%% IGRF vector function. %%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Bx, By, Bz] = igrfv(time, latitude, longitude, altitude, coord)
% Fundamental constant.
Rearth_km = 6371.2;
%%% CHECK INPUT VALIDITY %%%
% Convert time to a datenumber if it is a string.
if ischar(time)
time = datenum(time);
end
% Make sure time has only one element.
if numel(time) > 1
error('igrf:timeInputInvalid', ['The input TIME can only have one ' ...
'element.']);
end
% Check that the inputs all have either one or the same number of elements.
numlat = numel(latitude);
numlon = numel(longitude);
numalt = numel(altitude);
if numlat > 1
if numlon == 1
longitude = repmat(longitude, size(latitude));
end
if numalt == 1
altitude = repmat(altitude, size(latitude));
end
elseif numlon > 1
latitude = repmat(latitude, size(longitude));
if numalt == 1
altitude = repmat(altitude, size(latitude));
end
elseif numalt > 1
latitude = repmat(latitude, size(altitude));
longitude = repmat(longitude, size(altitude));
end
numlat = numel(latitude);
numlon = numel(longitude);
numalt = numel(altitude);
if numlat ~= numlon || numlat ~= numalt || numlon ~= numalt
error('igrf:inputNotSameSize', ['The input coordinates must have ' ...
'the same number of elements.']);
end
%%% SPHERICAL COORDINATE CONVERSION %%%
% Convert the latitude, longitude, and altitude coordinates input into
% spherical coordinates r (radius), theta (inclination angle from +z axis),
% and phi (azimuth angle from +x axis). Also, make the coordinates go down
% the rows.
% We want cos(theta) and sin(theta) rather than theta itself.
costheta = cos((90 - latitude(:))*pi/180);
sintheta = sin((90 - latitude(:))*pi/180);
% Convert from geodetic coordinates to geocentric coordinates if necessary.
% This method was adapted from igrf11syn, which was a conversion of the
% IGRF subroutine written in FORTRAN.
if nargin < 5 || isempty(coord) || strcmpi(coord, 'geodetic') || ...
strcmpi(coord, 'geod') || strcmpi(coord, 'gd')
a = 6378.137; f = 1/298.257223563; b = a*(1 - f);
rho = hypot(a*sintheta, b*costheta);
r = sqrt( altitude(:).^2 + 2*altitude(:).*rho + ...
(a^4*sintheta.^2 + b^4*costheta.^2) ./ rho.^2 );
cd = (altitude(:) + rho) ./ r;
sd = (a^2 - b^2) ./ rho .* costheta.*sintheta./r;
oldcos = costheta;
costheta = costheta.*cd - sintheta.*sd;
sintheta = sintheta.*cd + oldcos.*sd;
elseif strcmpi(coord, 'geocentric') || strcmpi(coord, 'geoc') || ...
strcmpi(coord, 'gc')
r = altitude(:);
cd = 1;
sd = 0;
else
error('igrf:coordInputInvalid', ['Unrecognized command ' coord ...
' for COORD input.']);
end
% Special case when sin(theta) = 0.
sintheta0 = sintheta == 0;
anysintheta0 = any(sintheta0);
anysinthetanot0 = any(~sintheta0);
% Convert longitude to radians.
phi = longitude(:)*pi/180;
%%% GET PROPER IGRF COEFFICIENTS %%%
[g, h] = loadigrfcoefs(time);
nmax = size(g, 1);
% We need cos(m*phi) and sin(m*phi) multiple times, so precalculate into a
% matrix here:
cosphi = cos(bsxfun(@times, 0:nmax, phi));
sinphi = sin(bsxfun(@times, 0:nmax, phi));
%%% BEGIN MAGNETIC FIELD CALCULATION %%%
% Initialize variables used in for loop below.
Br = zeros(size(r));
Bt = zeros(size(r));
Bp = zeros(size(r));
lastP = 1;
lastdP_1 = 0;
lastdP_2 = 0;
% Sum for each n value.
for n = 1 : nmax
m = 0 : n;
% Calculate legendre values. The output of the function has each m
% value going down the rows, but since m goes along the columns
% (coordinates go down the rows, remember?), permute it.
P = legendre(n, costheta, 'sch').';
% We also need the derivative of the legendre with respect to theta. It
% is given by a recursive function of both the previous legendre values
% as well as the previous derivatives. Functionally, it is:
% dP(0, 0) = 0, dP(1, 1) = cos(theta)
% dP(n, n) = sqrt(1 - 1/(2n))*(sin(theta)*dP(n-1, n-1) +
% cos(theta)*P(n-1, n-1))
% dP(n, m) = (2n - 1)/sqrt(n^2 - m^2)*(cos(theta)*dP(n-1, m) -
% sin(theta)*P(n-1, m)) - sqrt(((n-1)^2 - m^2)/(n^2 - m^2))*
% dP(n-2, m)
dP = [bsxfun(@minus, bsxfun(@times, ...
(2*n - 1)./sqrt(n^2 - m(1:end-1).^2), ...
bsxfun(@times, costheta, lastdP_1) - bsxfun(@times, sintheta, ...
lastP)), bsxfun(@times, sqrt(((n - 1)^2 - m(1:end-1).^2)./...
(n^2 - m(1:end-1).^2)), lastdP_2)), zeros(size(costheta))];
if n > 1
dP(:, end) = sqrt(1 - 1/(2*n))*...
(sintheta*lastdP_1(end) + costheta*lastP(end));
lastdP_2 = [lastdP_1 zeros(size(costheta))];
else
dP(:, end) = costheta;
lastdP_2 = lastdP_1;
end
lastP = P;
lastdP_1 = dP;
% Multiply coefficients by proper longitude trigonemetric term.
gcos = bsxfun(@times, g(n, m + 1), cosphi(:, m + 1));
gsin = bsxfun(@times, g(n, m + 1), sinphi(:, m + 1));
hcos = bsxfun(@times, h(n, m + 1), cosphi(:, m + 1));
hsin = bsxfun(@times, h(n, m + 1), sinphi(:, m + 1));
% Calculate the magnetic field components as a running sum. Find
% explicit expressions for these in Global Earth Physics: a Handbook of
% Physical Constants by Thomas J. Aherns (1995), pg. 49. Link:
% http://books.google.com/books?id=aqjU_NHyre4C&lpg=PP1&dq=Global%20
% earth%20physics%3A%20a%20handbook%20of%20physical%20constants&pg=PA49
% #v=onepage&q&f=false
% (except equation 6 is missing a required 1/sin(theta) and m; correct
% equations on page 5 (equations 3a-3c) of:
% http://hanspeterschaub.info/Papers/UnderGradStudents/
% MagneticField.pdf)
a_r = (Rearth_km./r).^(n + 2);
Br = Br + a_r.*(n+1).*sum((gcos + hsin).*P, 2);
Bt = Bt - a_r.*sum((gcos + hsin).*dP, 2);
% Different case when sin(theta) == 0 for phi component.
if anysinthetanot0
Bp(~sintheta0) = Bp(~sintheta0) - 1./sintheta(~sintheta0).*...
a_r(~sintheta0).*sum(bsxfun(@times, m, ...
(-gsin(~sintheta0, :) + hcos(~sintheta0, :)).*...
P(~sintheta0, :)), 2);
end
if anysintheta0
Bp(sintheta0) = Bp(sintheta0) - costheta(sintheta0).*...
a_r(sintheta0).*sum((-gsin(sintheta0, :) ...
+ hcos(sintheta0, :)).*dP(sintheta0, :), 2);
end
end
% Convert from spherical to (x,y,z) = (North,East,Down).
Bx = -Bt;
By = Bp;
Bz = -Br;
% Convert back to geodetic coordinates if necessary.
Bx_old = Bx;
Bx = Bx.*cd + Bz.*sd;
Bz = Bz.*cd - Bx_old.*sd;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%% IGRF scalar function. %%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Bx, By, Bz] = igrfs(time, latitude, longitude, altitude, coord)
% Fundamental constant.
Rearth_km = 6371.2;
%%% CHECK INPUT VALIDITY %%%
% Convert time to a datenumber if it is a string.
if ischar(time)
time = datenum(time);
end
% Check that the input coordinates are scalars.
if ~isscalar(latitude) || ~isscalar(longitude) || ~isscalar(altitude)
error('igrf1:inputNotScalar', ...
'The input coordinates must be scalars.');
end
%%% SPHERICAL COORDINATE CONVERSION %%%
% Convert the latitude, longitude, and altitude coordinates input into
% spherical coordinates r (radius), theta (inclination angle from +z axis),
% and phi (azimuth angle from +x axis).
% We want cos(theta) and sin(theta) rather than theta itself.
costheta = cos((90 - latitude)*pi/180);
sintheta = sin((90 - latitude)*pi/180);
% Convert from geodetic coordinates to geocentric coordinates if necessary.
% This method was adapted from igrf11syn, which was a conversion of the
% IGRF subroutine written in FORTRAN.
if nargin < 5 || isempty(coord) || strcmpi(coord, 'geodetic') || ...
strcmpi(coord, 'geod') || strcmpi(coord, 'gd')
a = 6378.137; f = 1/298.257223563; b = a*(1 - f);
rho = hypot(a*sintheta, b*costheta);
r = sqrt( altitude.^2 + 2*altitude.*rho + ...
(a^4*sintheta.^2 + b^4*costheta.^2) ./ rho.^2 );
cd = (altitude + rho) ./ r;
sd = (a^2 - b^2) ./ rho .* costheta.*sintheta./r;
oldcos = costheta;
costheta = costheta.*cd - sintheta.*sd;
sintheta = sintheta.*cd + oldcos.*sd;
elseif strcmpi(coord, 'geocentric') || strcmpi(coord, 'geoc') || ...
strcmpi(coord, 'gc')
r = altitude;
cd = 1;
sd = 0;
else
error('igrf:coordInputInvalid', ['Unrecognized command ' coord ...
' for COORD input.']);
end
% Convert longitude to radians.
phi = longitude*pi/180;
%%% GET PROPER IGRF COEFFICIENTS %%%
if isscalar(time)
gh = loadigrfcoefs(time);
nmax = sqrt(numel(gh) + 1) - 1;
% Assume a vector input means the coefficients are the input.
else
gh = time;
nmax = sqrt(numel(gh) + 1) - 1;
% nmax should be an integer.
if nmax - round(nmax) ~= 0
error('igrf:timeInputInvalid', ['TIME input should either be ' ...
'a single date or a valid coefficient vector.']);
end
end
% We need cos(m*phi) and sin(m*phi) multiple times, so precalculate into a
% vector here:
cosphi = cos((1:nmax)*phi);
sinphi = sin((1:nmax)*phi);
Pmax = (nmax+1)*(nmax+2)/2;
%%% BEGIN MAGNETIC FIELD CALCULATION %%%
% Initialize variables used in for loop below.
Br = 0; Bt = 0; Bp = 0;
P = zeros(1, Pmax); P(1) = 1; P(3) = sintheta;
dP = zeros(1, Pmax); dP(1) = 0; dP(3) = costheta;
% For this initial condition, the first if will result in n = 1, m = 0.
m = 1; n = 0; coefindex = 1;
a_r = (Rearth_km/r)^2;
% Increment through all the n's and m's. gh will be a vector with g
% followed by h for incrementing through n and m except when h would be
% redundant (i.e., when m = 0).
for Pindex = 2:Pmax
% Increment to the next n when m becomes larger than n.
if n < m
m = 0;
n = n + 1;
a_r = a_r*(Rearth_km/r); % We need (Rearth_km./r)^(n+2)
end
% Calculate P and dP. They are given recursively according to:
%
% P(0, 0) = 1, P(1, 1) = sin(theta) <- Specified above
% P(n, n) = sqrt(1 - 1/(2n))*sin(theta)*P(n-1, n-1)
% P(n, m) = (2n - 1)/sqrt(n^2 - m^2)*cos(theta)*P(n-1, m) -
% sqrt(((n-1)^2 - m^2) / (n^2 - m^2)) * P(n-2, m)
%
% dP(0, 0) = 0, dP(1, 1) = cos(theta) <- Specified above
% dP(n, n) = sqrt(1 - 1/(2n))*(sin(theta)*dP(n-1, n-1) +
% cos(theta)*P(n-1, n-1))
% dP(n, m) = (2n - 1)/sqrt(n^2 - m^2)*(cos(theta)*dP(n-1, m) -
% sin(theta)*P(n-1, m)) - sqrt(((n-1)^2 - m^2)/(n^2 - m^2))*
% dP(n-2, m)
if m < n && Pindex ~= 3 % (Pindex=3 is n=1, m=1, initial cond. above)
last1n = Pindex - n;
last2n = Pindex - 2*n + 1;
P(Pindex) = (2*n - 1)/sqrt(n^2 - m^2)*costheta*P(last1n) - ...
sqrt(((n-1)^2 - m^2) / (n^2 - m^2)) * P(last2n);
dP(Pindex) = (2*n - 1)/sqrt(n^2 - m^2)*(costheta*dP(last1n) - ...
sintheta*P(last1n)) - sqrt(((n-1)^2 - m^2) / (n^2 - m^2)) * ...
dP(last2n);
elseif Pindex ~= 3
lastn = Pindex - n - 1;
P(Pindex) = sqrt(1 - 1/(2*m))*sintheta*P(lastn);
dP(Pindex) = sqrt(1 - 1/(2*m))*(sintheta*dP(lastn) + ...
costheta*P(lastn));
end
% Calculate the magnetic field components as a running sum. Find
% explicit expressions for these in Global Earth Physics: a Handbook of
% Physical Constants by Thomas J. Aherns (1995), pg. 49. Link:
% http://books.google.com/books?id=aqjU_NHyre4C&lpg=PP1&dq=Global%20
% earth%20physics%3A%20a%20handbook%20of%20physical%20constants&pg=PA49
% #v=onepage&q&f=false
% (except equation 6 is missing a required 1/sin(theta) and m; correct
% equations on page 5 (equations 3a-3c) of:
% http://hanspeterschaub.info/Papers/UnderGradStudents/
% MagneticField.pdf)
if m == 0 % Implies h = 0, so only coefficient in gh is g
coef = a_r*gh(coefindex); %*cos(0*phi) = 1
Br = Br + (n+1)*coef*P(Pindex);
Bt = Bt - coef*dP(Pindex);
% Bp is 0 for m = 0.
coefindex = coefindex + 1; % Only need to skip over g this time.
else
coef = a_r*(gh(coefindex)*cosphi(m) + gh(coefindex+1)*sinphi(m));
Br = Br + (n+1)*coef*P(Pindex);
Bt = Bt - coef*dP(Pindex);
if sintheta == 0 % Use different formula when dividing by 0.
Bp = Bp - costheta*a_r*(-gh(coefindex)*sinphi(m) + ...
gh(coefindex+1)*cosphi(m))*dP(Pindex);
else
Bp = Bp - 1/sintheta*a_r*m*(-gh(coefindex)*sinphi(m) + ...
gh(coefindex+1)*cosphi(m))*P(Pindex);
end
coefindex = coefindex + 2; % Skip over g and h this time.
end
% Increment m.
m = m + 1;
end
% Convert from spherical to (x,y,z) = (North,East,Down).
Bx = -Bt;
By = Bp;
Bz = -Br;
% Convert back to geodetic coordinates if necessary.
Bx_old = Bx;
Bx = Bx.*cd + Bz.*sd;
Bz = Bz.*cd - Bx_old.*sd;
