Code covered by the BSD License  

Highlights from
MSIS-E-90

image thumbnail

MSIS-E-90

by

 

NRL MSIS-E-90 atmosphere model.

msistest.m
% MSISTEST Test MSIS function.

%% Clear variables and close all figures.
tic;
clear;
close all;

%% Initial values.
font = 'Times New Roman';
axis_font = 12;
title_font = 12;

time = datenum([2000 7 17 12 0 0]); % Must be a scalar or string.
latitude = -90:90; % Degrees.
longitude = -180:180; % Degrees.
altitude = 100; % km. Must be a scalar.

% %% Compute B-field (needs function IGRF from MATLAB file exchange).
% [LAT, LON] = meshgrid(latitude, longitude); s = size(LAT);
% [Bx, By, Bz] = igrf(time, LAT, LON, altitude, 'geod');
% Bx = reshape(Bx, s); By = reshape(By, s); Bz = reshape(Bz, s);
% decl = atan(By./Bx)*180/pi;
% B = hypot(Bx, hypot(By, Bz));

%% Compute atomic oxygen number density for all values above.
O = zeros(numel(latitude), numel(longitude));
if numel(latitude) >= numel(longitude) || altitude ~= 100
    for index = 1:length(longitude)
        [rho, T] = msis(time, latitude(:), longitude(index), altitude);
        O(:, index) = rho(:, 1);
        fprintf('%4.4g%% Complete\n', index/numel(longitude)*100);
    end
else
    for index = 1:numel(latitude)
        [rho, T] = msis(time, latitude(index), longitude(:), altitude);
        O(index, :) = rho(:, 1).';
        fprintf('%4.4g%% Complete\n', index/numel(latitude)*100);
    end
end

%% Plot data.
figure;
if ~license('test', 'MAP_Toolbox')
    hold on;
    % surf(LON.', LAT.', B.'); %decl
    surf(longitude, latitude, O.');
    shading flat;
    clim = get(gca, 'CLim'); zlim = get(gca, 'ZLim');
    load('topo.mat', 'topo'); topo = [topo(:, 181:360), topo(:, 1:180)];
    [C, h] = contour3(-180:179, -90:+89, topo + zlim(2), [0 0] + zlim(2));
    set(h, 'EdgeColor', 0.25*[1 1 1]);
    set(gca, 'XLim', [-180 180], 'XTick', [], ...
        'YLim', [-90 90], 'YTick', [], 'CLim', clim);
else
    axesm miller; axis fill;
    hold on;
    % surfm(LAT, LON, B); %decl
    surfm(latitude, longitude, O);
    load coast;
    plotm(lat, long, 'Color', 0.25*[1 1 1]);
end
hc = colorbar;
title(hc, 'cm^{-3}', 'FontName', font, 'FontSize', axis_font);
set(gca, 'FontName', font, 'FontSize', axis_font);
% title(['Magnetic Field Declination in degrees at \ith\rm = ' ...
%     sprintf('%g km', altitude) ' at ' datestr(time) ' UTC'], ...
%     'FontName', font, 'FontSize', title_font);
title(['Atomic Oxygen (O) Density at  \ith\rm = ' ...
    sprintf('%g km', altitude) ' at ' datestr(time) ' UTC'], ...
    'FontName', font, 'FontSize', title_font);
print -dpng -r100 msis.png

%%
toc;

Contact us