forked from wb-bgs/m_IGRF
-
Notifications
You must be signed in to change notification settings - Fork 0
/
getigrfcoefs.m
129 lines (114 loc) · 4.7 KB
/
getigrfcoefs.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
function out = getigrfcoefs(datadirectory)
% GETIGRFCOEFS Extract IGRF coefficients from .dat files.
%
% Usage: GETIGRFCOEFS(DATADIRECTORY)
% or COEFS = GETIGRFCOEFS(DATADIRECTORY)
%
% Use this function to extract the IGRF coefficients from .dat files
% provided on this IGRF website here:
%
% https://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
% is returned to the user as the output from this function if requested. If
% no outputs are requested, the variable COEFS is saved to the file
% igrfcoefs.mat in the same directory as this file as 'coefs' (i.e., in
% lower case).
%
% 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.
%
% See also: IGRF, LOADIGRFCOEFS.
% Get the directory with the datfiles.
if nargin < 1 || isempty(datadirectory)
fpath = fileparts(mfilename('fullpath'));
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
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'));
end
% Sort coefs based on the year.
yrs = cell2mat({coefs.year});
[~, sortind] = sort(yrs);
coefs = coefs(sortind);
if nargout >= 1
out = coefs;
else
% Save coefs.
save(fullfile(fpath, 'igrfcoefs.mat'), 'coefs');
end