Code covered by the BSD License  

Highlights from
International Geomagnetic Reference Field (IGRF) Model

image thumbnail

International Geomagnetic Reference Field (IGRF) Model

by

 

30 Dec 2011 (Updated )

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

getigrfcoefs(datadirectory)
function [out, years] = getigrfcoefs(datadirectory)

% GETIGRFCOEFS Extract IGRF coefficients from .dat files.
% 
% Usage: GETIGRFCOEFS(DATADIRECTORY)
%     or [COEFS, YEARS] = GETIGRFCOEFS(DATADIRECTORY)
% 
% Use this function to extract the IGRF coefficients from .dat files
% provided on the IGRF ftp here:
% 
% ftp://hanna.ccmc.gsfc.nasa.gov/pub/modelweb/geomagnetic/igrf/fortran_code/
% 
% The .dat files should be located in the directory DATADIRECTORY. If
% DATADIRECTORY is unspecified, this function will look for a directory
% called "datfiles" within the directory that this file is in or, if that
% does not exist, will just use the directory that this file is in. The
% coefficients are stored as a structure array with the following fields
% (case-sensitive):
% 
%    -year: Year for the coefficients of that index in the array.
%    -gh: g and h coefficients in a vector.
%    -g: g coefficients in a matrix with each n going down the rows and
%    all of the m's for each n going along the columns.
%    -h: h coefficients in a matrix with each n going down the rows and
%    all of the m's for each n going along the columns.
%    -slope: True if the coefficients are slopes, false otherwise.
% 
% The elements in the structure array will be sorted by year (i.e., the
% first structure in the array will have the earliest year and the last
% structure in the array will have the latest year). This structure array
% and a row vector corresponding to the year of each element in the
% structure array are returned to the user as COEFS and YEARS and no files
% are created if an output is requested from this function. If no outputs
% are requested, the variables COEFS and YEARS are saved to the file 
% igrfcoefs.mat in the same directory as this file as 'coefs' and 'years', 
% respectively.
% 
% Inputs:
%   -DATADIRECTORY: directory where the .dat IGRF coefficients can be
%   found (optional, default is the directory datfiles within the same
%   directory that this function is located or, if that does not exist, the
%   directory where this function is located.
% 
% Outputs:
%   -COEFS: Structure array with fields as explained above.
%   -YEARS: Year of each element in COEFS.
% 
% See also: IGRF, LOADIGRFCOEFS.

% Get the directory with the datfiles.
if nargin < 1 || isempty(datadirectory)
    fpath = mfilename('fullpath');
    fpath = fpath(1:end-length(mfilename));
    if exist(fullfile(fpath, 'datfiles'), 'dir')
        datadirectory = fullfile(fpath, 'datfiles');
    else
        datadirectory = fpath;
    end
end

% Get all the valid coefficient files.
datfiles = dir(fullfile(datadirectory, '*grf*.dat'));

if isempty(datfiles)
    error('getigrfcoef:filesNotFound', ['No valid .dat files found in ' ...
        datadirectory]);
end

years = zeros(1, numel(datfiles));
coefs = struct('year', [], 'g', [], 'h', [], 'gh', [], 'slope', []);
for index = 1:numel(datfiles)
    
    % Read the current file.
    fid = fopen(fullfile(datadirectory, datfiles(index).name));
    if fid == -1
        error('getigrfcoef:fileOpenError', ['Cannot open ' ...
            datfiles(index).name '.']);
    end
    line1 = fgetl(fid);
    line2 = fscanf(fid, '%d %g %g', [1 3]);
    thisnmax = line2(1); thisr = line2(2); thisyear = line2(3);
    thiscoefs = fscanf(fid, '%g', Inf);
    
    % Generate a vector that has all the m indices. It will go
    % [0 1 0 1 2 0 1 2 3 ...].
    mmat = tril(repmat(1:thisnmax, thisnmax, 1))';
    mmat(mmat == 0) = -1;
    mmat = [zeros(1, thisnmax); mmat];
    m = mmat(mmat ~= -1);
    
    % Also generate a vector that has all the n indices. It will go
    % [1 1 2 2 2 3 3 3 3 ...].
    nmat = [1:thisnmax; triu(repmat(1:thisnmax, thisnmax, 1))];
    n = nmat(nmat ~= 0);
    
    % thiscoefs omits h when m = 0. However, it is easier to just set h = 0
    % in those places. Find where those places should be in a full vector:
    m0 = m; m0(m0 == 0) = -1;
    h0 = [m'; m0'] == -1; h0 = h0(:);
    
    % Extract the g and h coefficients from thiscoefs.
    ncoefs = 2*sum((1:thisnmax) + 1);
    coefsvec = zeros(ncoefs, 1);
    coefsvec(~h0) = thiscoefs;
    gvec = coefsvec(1:2:end); hvec = coefsvec(2:2:end);
    
    % Put the coefficients in a matrix with n going down the rows and m
    % down the columns.
    thisg = zeros(thisnmax, thisnmax+1);
    thisg(sub2ind(size(thisg), n, m+1)) = gvec;
    thish = zeros(thisnmax, thisnmax+1);
    thish(sub2ind(size(thisg), n, m+1)) = hvec;
    
    % Save the data to coefs.
    coefs(index).year = thisyear;
    coefs(index).g = thisg;
    coefs(index).h = thish;
    coefs(index).gh = thiscoefs;
    coefs(index).slope = ~isempty(regexpi(datfiles(index).name, ...
        '\wgrf\d\d\d\ds.dat', 'once'));
    
    % Save the year for sorting later.
    years(index) = thisyear;
    
end

% Sort coefs based on the year.
[tmp, sortind] = sort(years);
coefs = coefs(sortind); years = years(sortind);

if nargout >= 1
    out = coefs;
else
    % Save coefs.
    save(fullfile(fpath, 'igrfcoefs.mat'), 'years', 'coefs');
end

Contact us