Skip to content
Permalink
Branch: master
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
313 lines (254 sloc) 8.28 KB
function lmcosi=igrf10(yr,yir)
% lmcosi=IGRF10(yr,yir)
%
% Interface to load the International Geomagnetic Reference Field
% and pass it on to other subroutines.
%
% INPUT:
%
% yr The year of interest (out of 1900:5:2005) [default: 2005]
% OR: a string with the demo numbre
% yir The year of interest if the first argument is a demo string
%
% OUTPUT:
%
% lmcosi The tradition matrix with the ordered real coefficients
% for the potential - however, in units of nT (nanoTesla)
%
% EXAMPLE:
%
% igrf10('demo1')
% igrf10('demo2') % The radial non-dipolar field, in nanoTesla
% igrf10('demo3') % The radial field, contoured, in nanoTesla
% igrf10('demo4') % The radial non-dipolar field, contoured, in nanoTesla
% igrf10('demo5') % The radial non-dipolar field
%
% SEE ALSO:
%
% PLM2MAG
%
% Last modified by fjsimons-at-alum.mit.edu, 10/11/2010
% See plates at: http://pubs.usgs.gov/sim/2007/2964/
defval('yr',2005)
if ~isstr(yr)
% Make it a single year - perhaps fix later
if prod(size(yr))~=1
error('Only a single year at the time can be requested for now')
end
% Open file
fid=fopen(fullfile(getenv('IFILES'),...
'EARTHMODELS','IGRF-10','igrf10coeffs.txt'));
% Define formats
fmt1=['%s %s %s' repmat('%n',1,22) '%s'];
fmt2=['%s' repmat('%n',1,25)];
% The maximum expansion
lmax=[repmat(10,1,20) 13 13];
% Read the first line % TEXTSCAN better than TEXTREAD
d=textscan(fid,fmt1,1);
% Read the rest - supply zeroes where unavailable
e=textscan(fid,fmt2,'emptyvalue',0);
% Close file
fclose(fid);
% Available years
years=[d{4:25}];
% Spherical harmonic degrees
EL=e{2};
% Spherical harmonic orders
EM=e{3};
% Secular variation
SV=e{26}(~isnan(e{26}));
% The actual field expansion coefficients
% Work from the back - you've got to know what's going on
cosi=[e{4:25}];
% Work-around the known knowns and the known unknowns
ilmax=addmoff(lmax)-1;
% This belongs with the last two
cosi(ilmax(1)+1:end,21)=cosi(ilmax(1)+1:end,1);
cosi(ilmax(1)+1:end,22)=cosi(ilmax(2)+1:end,2);
cosi(ilmax(1)+1:end,1)=0;
cosi(ilmax(2)+1:end,2)=0;
% Now extract the data
[C,iy]=intersect(years,yr);
prepar=cosi(:,iy);
% Stick in the non-existing zeros for degree and order zero
prepar=[zeros(1,size(prepar,2)) ; prepar];
% Reordering sequence
[dems,dels,mz,lmcosi,mzi,mzo,bigm,bigl,rinm,ronm,demin]=...
addmon(max(lmax));
% What we have from the file is, effectively
% [bigl(2:end) bigm(2:end)]
% and what we want is lmcosi with the coefficients in the right position
% This is the output
lmcosi(mzo+2*size(lmcosi,1))=prepar;
elseif strcmp(yr,'demo1')
clf
yir=2005;
h=igrf10(yir);
% Change Schmidt to full normalization for use in PLOTPLM
% This converts the COEFFICIENTS to be multiplied with Schmidt to the
% COEFFICIENTS to be multiplied with 4pi-normalized harmonics, which
% are Schmidt*sqrt(2l+1), see PLM2XYZ and note the TYPO in Blakely.
h(:,3:4)=h(:,3:4)./repmat(sqrt(2*h(:,1)+1),1,2);
% Make sure it is the RADIAL component of this at the surface
h(:,3:4)=repmat(h(:,1)+1,1,2).*h(:,3:4);
d=plotplm(h,[],[],4,1);
longticks(gca,2)
t(1)=title(sprintf('IGRF-10 magnetic field, year %i, degrees %i-%i',yir,...
h(min(find(h(:,3))),1),h(end,1)));
movev(t,5)
cb=colorbar('hor');
shrink(cb,2,2)
axes(cb)
longticks(cb,2)
xlabel('radial component (nT)')
movev(cb,-.1)
fig2print(gcf,'portrait')
figna=figdisp([],sprintf('%s-%i',yr,yir),[],1);
system(sprintf('degs %s.eps',figna));
system(sprintf('epstopdf %s.eps',figna));
elseif strcmp(yr,'demo2')
defval('yir',2005);
h=igrf10(yir);
% Change Schmidt to full normalization for use in PLOTPLM
h(:,3:4)=h(:,3:4)./repmat(sqrt(2*h(:,1)+1),1,2);
% The nondipole field, as Blakely p170 and eq. (8.20)
h(1:3,3:4)=0;
% Make sure it is the RADIAL component of this at the surface
h(:,3:4)=repmat(h(:,1)+1,1,2).*h(:,3:4);
clf
d=plotplm(h,[],[],4,1);
axis image
longticks(gca,2)
t(1)=title(sprintf('IGRF-10 magnetic field, year %i, degrees %i-%i',yir,...
h(min(find(h(:,3))),1),h(end,1)));
movev(t,5)
cb=colorbar('hor');
shrink(cb,2,2)
axes(cb)
longticks(cb,2)
xlabel('non-dipolar radial component (nT)');
movev(cb,-.1)
fig2print(gcf,'portrait')
figna=figdisp([],sprintf('%s-%i',yr,yir),[],1);
system(sprintf('degs %s.eps',figna));
system(sprintf('epstopdf %s.eps',figna));
elseif strcmp(yr,'demo3')
clf
yir=2005;
h=igrf10(yir);
% Change Schmidt to full normalization for use in PLOTPLM
h(:,3:4)=h(:,3:4)./repmat(sqrt(2*h(:,1)+1),1,2);
% Make sure it is the RADIAL component of this at the surface
h(:,3:4)=repmat(h(:,1)+1,1,2).*h(:,3:4);
d=plotplm(h,[],[],4,1); clf
lons=linspace(0,360,size(d,2));
lats=linspace(-90,90,size(d,1));
% Don't forget to flip up down for contouring!
[c,hh]=contour(lons,lats,...
flipud(d),[-65000:5000:-5000]);
set(hh,'EdgeC','r')
hold on
[c,hh]=contour(lons,lats,...
flipud(d),[5000:5000:65000]);
set(hh,'EdgeC','b')
[c,hh]=contour(lons,lats,...
flipud(d),[0 0]);
set(hh,'EdgeC','k','LineW',2)
plotcont; axis image; ylim([-90 90])
defval('dlat',45)
set(gca,'ytick',[-90:dlat:90])
set(gca,'xtick',[0:90:360])
deggies(gca)
longticks(gca,2)
t(1)=title(sprintf('IGRF-10 magnetic field, year %i, degrees %i-%i',yir,...
h(min(find(h(:,3))),1),h(end,1)));
movev(t,5)
xl=xlabel(sprintf('minimum %i nT ; maximum %i nT',round(min(d(:))), ...
round(max(d(:)))));
movev(xl,-10)
fig2print(gcf,'portrait')
figna=figdisp([],sprintf('%s-%i',yr,yir),[],1);
system(sprintf('degs %s.eps',figna));
system(sprintf('epstopdf %s.eps',figna));
elseif strcmp(yr,'demo4')
clf
yir=2005;
h=igrf10(yir);
% The nondipole field, as Blakely p170 and eq. (8.20)
h(1:3,3:4)=0;
% Change Schmidt to full normalization for use in PLOTPLM
h(:,3:4)=h(:,3:4)./repmat(sqrt(2*h(:,1)+1),1,2);
% Make sure it is the RADIAL component of this at the surface
h(:,3:4)=repmat(h(:,1)+1,1,2).*h(:,3:4);
d=plotplm(h,[],[],4,1); clf
lons=linspace(0,360,size(d,2));
lats=linspace(-90,90,size(d,1));
clf
% Don't forget to flip up down for contouring!
[c,hh]=contour(lons,lats,...
flipud(d),[-20000:1000:-1000]);
set(hh,'EdgeC','r')
hold on
[c,hh]=contour(lons,lats,...
flipud(d),[1000:1000:20000]);
set(hh,'EdgeC','b')
[c,hh]=contour(lons,lats,...
flipud(d),[0 0]);
set(hh,'EdgeC','k','LineW',2)
plotcont; axis image; ylim([-90 90])
defval('dlat',45)
set(gca,'ytick',[-90:dlat:90])
set(gca,'xtick',[0:90:360])
deggies(gca)
longticks(gca,2)
t(1)=title(sprintf('IGRF-10 magnetic field, year %i, degrees %i-%i',yir,...
h(min(find(h(:,3))),1),h(end,1)));
movev(t,5)
xl=xlabel(sprintf('minimum %i nT ; maximum %i nT',round(min(d(:))), ...
round(max(d(:)))));
movev(xl,-10)
fig2print(gcf,'portrait')
figna=figdisp([],sprintf('%s-%i',yr,yir),[],1);
system(sprintf('degs %s.eps',figna));
system(sprintf('epstopdf %s.eps',figna));
elseif strcmp(yr,'demo5')
clf
defval('yir',2005);
h=igrf10(yir);
% The nondipole field, as Blakely p170 and eq. (8.20)
h(1:3,3:4)=0;
% Change Schmidt to full normalization for use in PLOTPLM
h(:,3:4)=h(:,3:4)./repmat(sqrt(2*h(:,1)+1),1,2);
% Make sure it is the RADIAL component of this at the surface
h(:,3:4)=repmat(h(:,1)+1,1,2).*h(:,3:4);
d=plotplm(h,[],[],4,1);
lons=linspace(0,360,size(d,2));
lats=linspace(-90,90,size(d,1));
% Don't forget to flip up down for contouring!
[c,hh]=contour(lons,lats,...
flipud(d),[-20000:2000:-2000]);
set(hh,'EdgeC','r')
hold on
[c,hh]=contour(lons,lats,...
flipud(d),[2000:2000:20000]);
set(hh,'EdgeC','b')
[c,hh]=contour(lons,lats,...
flipud(d),[0 0]);
set(hh,'EdgeC','k','LineW',2)
plotcont; axis image; ylim([-90 90])
defval('dlat',45)
set(gca,'ytick',[-90:dlat:90])
set(gca,'xtick',[0:90:360])
deggies(gca)
longticks(gca,2)
t(1)=title(sprintf('IGRF-10 magnetic field, year %i, degrees %i-%i',yir,...
h(min(find(h(:,3))),1),h(end,1)));
movev(t,5)
xl=xlabel(sprintf('minimum %i nT ; maximum %i nT',round(min(d(:))), ...
round(max(d(:)))));
movev(xl,-10)
fig2print(gcf,'portrait')
figna=figdisp([],sprintf('%s-%i',yr,yir),[],1);
system(sprintf('degs %s.eps',figna));
system(sprintf('epstopdf %s.eps',figna));
end
You can’t perform that action at this time.
You signed in with another tab or window. Reload to refresh your session. You signed out in another tab or window. Reload to refresh your session.