close all;clear all;clc;
height = 500000; % A scalar value, in meters
lat = -90:2.5:90;  % geodetic latitude, in degrees, where north latitude is
%positive, and south latitude is negative
lon = -180:2.5:180;  % geodetic longitude, in degrees,
%where east longitude is positive,and west longitude is negative
% dyearA scalar decimal year. Decimal year is the desired year in
%a decimal format to include any fraction of the year that has already passed.
sla = size(lat);
sln = size(lon);
for i = 1:sla(2)
    for j =1:sln(2)
        [Bxyz, h, dec, dip, f] = ...
        wrldmagm(height, lat(i), lon(j), decyear(2014,1,22),'2010');
    % Note that WMM-2010 isvalid from January 1, 2010, until December 31,
    % 2014
    % Bxyz Magnetic field vector in nanotesla (nT)
    % h Horizontal intensity in nanotesla (nT)
    % dec Declination in degrees
    % dip Inclination in degrees
    % f Total intensity in nanotesla (nT)
        Bx(i,j) = Bxyz(1,1);
        By(i,j) = Bxyz(2,1);
        Bz(i,j) = Bxyz(3,1);
        Bm(i,j) = f;
%          [Bxyz(i,j), h(i,j), dec(i,j), dip(i,j), f(i,j)] = ...
%             wrldmagm(height, lat(i), lon(j), decyear(2014,1,22),'2010');
    end
end

Z component

xwidth = 640;
ywidth = 320;
colormap jet
hFig = figure(1);
contourf(lon,lat,Bz)
hold on
set(gcf,'PaperPositionMode','auto')
set(hFig, 'Position', [100 100 xwidth ywidth])
axis([-180 180 -90 90])
colorbar
load('topo.mat','topo','topomap1')
contour(0:179,-90:89,topo(:,1:180),[0 0],'.k')
contour(-180:-1,-90:89,topo(:,181:360),[0 0],'.k')
ylabel('Latitude [deg]')
xlabel('Longitude [deg]')
title('World Magnetic Model, Bz(nT) at Alt = 500 km')

Component X

hFig = figure(2);
contourf(lon,lat,Bx)
hold on
set(gcf,'PaperPositionMode','auto')
set(hFig, 'Position', [100 100 xwidth ywidth])
axis([-180 180 -90 90])
colorbar
load('topo.mat','topo','topomap1')
contour(0:179,-90:89,topo(:,1:180),[0 0],'.k')
contour(-180:-1,-90:89,topo(:,181:360),[0 0],'.k')
ylabel('Latitude [deg]')
xlabel('Longitude [deg]')
title('World Magnetic Model, Bx(nT) at Alt = 500 km')

Component Y

hFig = figure(3);
contourf(lon,lat,By)
hold on
set(gcf,'PaperPositionMode','auto')
set(hFig, 'Position', [100 100 xwidth ywidth])
axis([-180 180 -90 90])
colorbar
load('topo.mat','topo','topomap1')
contour(0:179,-90:89,topo(:,1:180),[0 0],'.k')
contour(-180:-1,-90:89,topo(:,181:360),[0 0],'.k')
ylabel('Latitude [deg]')
xlabel('Longitude [deg]')
title('World Magnetic Model, By(nT) at Alt = 500 km')

Total intensity in nanotesla (nT)

hFig = figure(4);
contourf(lon,lat,Bm)
hold on
set(gcf,'PaperPositionMode','auto')
set(hFig, 'Position', [100 100 xwidth ywidth])
axis([-180 180 -90 90])
colorbar
load('topo.mat','topo','topomap1')
contour(0:179,-90:89,topo(:,1:180),[0 0],'.k')
contour(-180:-1,-90:89,topo(:,181:360),[0 0],'.k')
ylabel('Latitude [deg]')
xlabel('Longitude [deg]')
title('World Magnetic Model, Total intensity(nT) at Alt = 500 km')