% Program to read the shell output (binary files) of BATS-R-US.
%
% The program can read the output of both the multi-species and the
% multi-fluid MHD codes.
%
% Input parameter: file name (.out)
%
% Output parameters: spherical grid (radius, longitude and latitude) and
% variables of the model (vars). Each column of the vars array corresponds
% to a specific variable as specified by the lines printed when running the
% program (together with the corresponding units).
%
% This program was prepared to read the outputs of the Mars version of
% BATS-R-US.
%
% Written by Leonardo Regoli, April 2018.
% Contact: lregoli@umich.edu

function [r,lon,lat,vars] = read_BATSRUS_shl(filename)
fid = fopen(filename);

header_size = fread(fid,1,'long');
header = fread(fid,header_size,'*char')';
header_size = fread(fid,1,'long');
length1 = fread(fid,1,'long');
iter=fread(fid,1,'long');
simtime=fread(fid,1,'float');
numdim=fread(fid,1,'long');
numpar = fread(fid,1,'long');
numvar = fread(fid,1,'long');
length1 = fread(fid,1,'int');
length1 = fread(fid,1,'int');
gridsizex = fread(fid,1,'long');
gridsizey = fread(fid,1,'long');
gridsizez = fread(fid,1,'long');
length1 = fread(fid,1,'int');
length1 = fread(fid,1,'int');
varspars = fread(fid,length1,'*char')';
varsparsaux=varspars;
disp(varspars)
disp(header)
count=size(varsparsaux,2);
while(varsparsaux(count) == ' ')
    varsparsaux(count) = [];
    count = size(varsparsaux,2);
end
blankspaces = strfind(varsparsaux,' ');
for count = 1:size(blankspaces,2);
    if count == 1
        varsparsaux(1:blankspaces(count)-1);
    elseif count == size(blankspaces,2)
        varsparsaux(blankspaces(count)+1:end);
    else
        varsparsaux(blankspaces(count-1)+1:blankspaces(count)-1);
    end
    
end
length1 = fread(fid,1,'int');
npoints = gridsizex * gridsizey * gridsizez;
length1 = fread(fid,1,'int');
r=fread(fid,npoints,'float');
lon=fread(fid,npoints,'float');
lat=fread(fid,npoints,'float');
length1 = fread(fid,1,'int');
dummylongs = zeros(2,size(blankspaces,2)-2);
vars = zeros(npoints,size(blankspaces,2)-2);
for count = 1:size(blankspaces,2)-2
    dummylongs(1,count) = fread(fid,1,'int');
    vars(:,count) = fread(fid,npoints,'float');
    dummylongs(2,count) = fread(fid,1,'int');
end
fclose(fid);
end