% Program to read the full 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) and type ('msmhd' for multi-species and
% 'mfmhd' for multi-fluid).
%
% 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 [xmhd,ymhd,zmhd,vars] = read_BATSRUS_ful(filename,type)
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');
if(strcmp(type,'msmhd'))
    length1 = fread(fid,1,'int');
    param1 = fread(fid,1,'float');
    param2 = fread(fid,1,'float');
    param3 = fread(fid,1,'float');
    length1 = fread(fid,1,'int');
end
if(strcmp(type,'mfmhd'))
    length1 = fread(fid,1,'int');
    param1 = fread(fid,1,'float');
    param2 = fread(fid,1,'float');
    param3 = fread(fid,1,'float');
    param4 = fread(fid,1,'float');
    param5 = fread(fid,1,'float');
    param6 = fread(fid,1,'float');
    param7 = fread(fid,1,'float');
    length1 = fread(fid,1,'int');
end
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');
xmhd=fread(fid,npoints,'float');
ymhd=fread(fid,npoints,'float');
zmhd=fread(fid,npoints,'float');
length1 = fread(fid,1,'int');
if strcmp(type,'msmhd')
    dummy = 5;
elseif strcmp(type,'mfmhd')
    dummy = 9;
end
dummylongs = zeros(2,size(blankspaces,2)-dummy);
vars = zeros(npoints,size(blankspaces,2)-dummy);
for count = 1:size(blankspaces,2)-dummy
    dummylongs(1,count) = fread(fid,1,'int');
    vars(:,count) = fread(fid,npoints,'float');
    dummylongs(2,count) = fread(fid,1,'int');
end
fclose(fid);
end