%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% msis.m provides a background atmosphere 
% obtained from the MSIS atmospheric model. 
% Profiles of total number density and temperature are 
% obtained from an MSIS look-up-table in z, month and 
% latitude.
%
% input:  z - altitude in [km]
%         lat - latitude in [degree]
%         month - month
%         day - day
%
% output: T - temperature in[K]
%         Ntot - air density in [m-3]
%         T & Ntot are returned as function of altitude x month/day x latitude
%         in case of multidimensional input
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         

function [T,Ntot] = msis (z,lat,month,day);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         

% ##### DEFAULT OUTPUT #####

T=[];
Ntot=[];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         

% ##### LOAD DATA #####

% load look-up-table:
datafile='msisdata.mat';

test=exist(datafile,'file');
if test ~= 2
warning('Data file does not exist!')
return
end

load(datafile);   

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         

% ##### HANDLING #####

% look up as a function of time and latitude
Ntot = 1e6 * NMsis(:, round(month+day/30), round(lat/10)+10);
  
T = TMsis(:, round(month+day/30), round(lat/10)+10);

% interpolate to requested altitudes:
Ntot = exp (interp1 (zMsis, log(Ntot), z));  % logarithmic interpolation
T = interp1 (zMsis, T, z);                   % linear interpolation

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%