Code covered by the BSD License

# International Geomagnetic Reference Field (IGRF) Model

### Drew Compston (view profile)

30 Dec 2011 (Updated )

Compute the Earth's magnetic field at points in space according to the IGRF model.

igrfline(time, lat_start, lon_start, alt_start, coord, distance, nsteps)
```function [latitude, longitude, altitude] = igrfline(time, lat_start, ...
lon_start, alt_start, coord, distance, nsteps)

% IGRFLINE Trace IGRF magnetic field line.
%
% Usage: [LATITUDE, LONGITUDE, ALTITUDE] = IGRFLINE(TIME, LAT_START,
%           LON_START, ALT_START, COORD, DISTANCE, NSTEPS)
%     or LLA = IGRFLINE(...)
%
% Gives the coordinates of the magnetic field line starting at a given
% latitude, longitude, and altitude. A total of NSTEPS points on the field
% line over a distance DISTANCE are given. Note that the step length
% (DISTANCE/NSTEPS) should be small (the smaller it is, the more accurate
% the results will be). The output coordinates are NSTEPS+1 long column
% vectors with LAT_START, LON_START, and ALT_START being the first element
% in each vector. The input coordinates can either be in the geocentric or
% geodetic (default) system (specified by COORD), and the output will be in
% the same system as the input. If just one output argument is requested,
% LATITUDE, LONGITUDE, and ALTITUDE, are output as an NSTEPS+1-by-3 matrix
% as LLA = [LATITUDE, LONGITUDE, ALTITUDE].
%
% 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.
%
% Inputs:
%   -TIME: Time to get the magnetic field line coordinates 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
%   -LAT_START: Starting point geocentric or geodetic latitude in degrees.
%   -LON_START: Starting point geocentric or geodetic longitude in degrees.
%   -ALT_START: For geodetic coordiates, the starting height in km above
%   the Earth's surface. For geocentric coordiates, the starting 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).
%   -DISTANCE: Distance in km to trace along the magnetic field line.
%   Positive goes up the field line (towards the geographic north pole /
%   geomagnetic south pole), while negative goes the other way.
%   -NSTEPS: Number of steps to compute the magnetic field line.
%
% Outputs:
%   -LATITUDE: Column vector of the geocentric or geodetic latitudes in
%   degrees of the NSTEPS+1 points along the magnetic field line.
%   -LONGITUDE: Column vector of the geocentric or geodetic longitudes in
%   degrees of the NSTEPS+1 points along the magnetic field line.
%   -ALTITUDE: For geodetic coordiates, the heights in km above the Earth's
%   surface of the NSTEPS+1 points along the magnetic field line. For
%   geocentric coordiates, the radii in km from the center of the Earth of
%   the NSTEPS+1 points along the magnetic field line. Also a column
%   vector.
%   LLA: NSTEPS+1-by-3 matrix of [LATITUDE, LONGITUDE, ALTITUDE].
%

error(nargchk(7, 7, nargin));

% Length of each step.
steplen = distance/abs(nsteps);

% Convert from geodetic coordinates to geocentric coordinates if necessary.
% The coordinate system used here is spherical coordinates (r,phi,theta)
% corresponding to radius, azimuth, and elevation, respectively.
if isempty(coord) || strcmpi(coord, 'geodetic') || ...
strcmpi(coord, 'geod') || strcmpi(coord, 'gd')
% First convert to ECEF, then convert to spherical. The function
% geod2ecef assumes meters, but we use km here.
[x, y, z] = geod2ecef(lat_start, lon_start, alt_start*1e3);
[phi, theta, r] = cart2sph(x, y, z); r = r/1e3;
elseif strcmpi(coord, 'geocentric') || strcmpi(coord, 'geoc') || ...
strcmpi(coord, 'gc')
r = alt_start;
phi = lon_start*pi/180;
theta = lat_start*pi/180;
else
error('igrfline:coordInputInvalid', ['Unrecognized command ' coord ...
' for COORD input.']);
end

% Initialize for loop.
r = [r; zeros(nsteps, 1)];
phi = [phi; zeros(nsteps, 1)];
theta = [theta; zeros(nsteps, 1)];

for index = 1:nsteps

% Get magnetic field at this point. Note that IGRF outputs the
% Northward (x), Eastward (y), and Downward (z) components, but we want
% the radial (-z), azimuthal (y), and elevation (x) components
% corresponding to (r,phi,theta), respectively.
[Bt, Bp, Br] = igrf(gh, theta(index)*180/pi, phi(index)*180/pi, ...
r(index), 'geoc'); Br = -Br;
B = hypot(Br, hypot(Bp, Bt));

% Unit vector in the (r,phi,theta) direction:
dr = Br/B; dp = Bp/B; dt = Bt/B;

% The next point is steplen km from the previous point in the direction
% of the unit vector just found above.
r(index+1) = r(index) + steplen*dr;
theta(index+1) = theta(index) + steplen*dt/r(index);
phi(index+1) = phi(index) + steplen*dp/(r(index)*cos(theta(index)));

end

% Convert the field line to the proper coordinate system.
if isempty(coord) || strcmpi(coord, 'geodetic') || ...
strcmpi(coord, 'geod') || strcmpi(coord, 'gd')
% First convert to ECEF, then geodetic. The function ecef2geod assumes
% meters, but we want km here.
[x, y, z] = sph2cart(phi, theta, r*1e3);
[latitude, longitude, altitude] = ecef2geod(x, y, z);
altitude = altitude/1e3;
else
latitude = theta*180/pi;
longitude = phi*180/pi;
altitude = r;
end

if nargout <= 1
latitude = [latitude(:), longitude(:), altitude(:)];
end```