Skip to content

Commit

Permalink
matlab demo [ci skip]
Browse files Browse the repository at this point in the history
  • Loading branch information
scivision committed Jul 24, 2018
1 parent 75225f5 commit f213554
Show file tree
Hide file tree
Showing 5 changed files with 63 additions and 1 deletion.
1 change: 1 addition & 0 deletions .gitignore
@@ -1,3 +1,4 @@
*.m~
*.svg
bin/
.mypy_cache/
Expand Down
2 changes: 1 addition & 1 deletion MANIFEST.in
@@ -1,2 +1,2 @@
include LICENSE tests/msise00_demo.gif
include LICENSE tests/msise00_demo.gif tests/msis_matlab.png
recursive-include tests *.py *.nc
4 changes: 4 additions & 0 deletions README.md
Expand Up @@ -17,6 +17,10 @@ The yellow ball represents the sun footprint on Earth.

![MSIS global time animation](tests/msise00_demo.gif)

This plot is from Matlab:

![MSISE00 Matlab](tests/msis_matlab.png)

## Install

- Mac: `brew install gcc`
Expand Down
Binary file added tests/msis_matlab.png
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
57 changes: 57 additions & 0 deletions tests/test_msise00_matlab.m
@@ -0,0 +1,57 @@
function test_msise00_matlab()
time = {'2013-03-31', '2013-04-01'};
altkm = 150.;
glat = 65;
glon = -148.;

%% run MSISE00
atmos = py.msise00.run(time, altkm, glat, glon);
%% Extract time
t = datetime(xarrayInd2vector(atmos, 'time') / 1e9, 'convertfrom', 'posixtime');

assert(t(1) == datetime(time{1}))
%% extract values
He = xarrayDataArray2mat(atmos{'He'});
N2 = xarrayDataArray2mat(atmos{'N2'});
O = xarrayDataArray2mat(atmos{'O'});

assert_allclose(N2(13), 3.051389580214272e16)
%% plot

figure()
plot(t, N2)
ylabel('density [m^{-3}]')
title('N_2 vs. time')
end


function V = xarrayDataArray2mat(V)
% convert xarray DataArray to Matlab matrix


V= V.values;
S = V.shape;
V = cell2mat(cell(V.ravel('F').tolist()));
V = reshape(V,[int64(S{1}), int64(S{2})]);

end

function I = xarrayInd2vector(V,key)

Iv = V{key}.indexes{key}.values.tolist();

I = cellfun(@double, cell(Iv));

end

function assert_allclose(actual, desired, atol,rtol)
narginchk(2,4)

if nargin<3 || isempty(atol), atol=0; end
if nargin<4 || isempty(rtol), rtol=1e-7; end

measdiff = abs(actual-desired);
tol = atol + rtol * abs(desired);
assert(all(measdiff(:) <= tol))

end

0 comments on commit f213554

Please sign in to comment.