function [data1, data2] = igwread(fname, seq)
% Loads IGW binary output. The file "startup_graphics" that is written
% by IGW should be present with the output files.
%
% To load a specific file, use the one parameter form:
%      D  = igwread('bin_D5');
%   Jinv  = igwread('bin_inv_Jac');
%   [x,z] = igwread('bin_vgrid');
% [px,pz] = igwread('bin_pressure9');
% sgpbotx = igwread('bin_sgpbotx');
%
% You can also let this function build out the filename, such as inside
% of a loop, by using the two parameter form:
%       D = igwread('bin_D',5);
%       D = igwread('D',5);
% [px,pz] = igwread('pressure',9);
% for ii=1:10, u = igwread('U', ii); pcolor(u); end
%
% v0.13, Aug 2013, Michael Dunphy, mdunphy@math.uwaterloo.ca

if nargin==0, disp([mfilename ': must specify file to read!']); return; end
if nargin==2,
    % Build out filename
    fname=[fname num2str(seq)];
    if ~exist(fname,'file'), fname=['bin_' fname]; end
end
if ~exist(fname,'file')
    disp([mfilename ': unable to find ' fname ', exiting!']);
    data1=[]; data2=[]; return;
end;

% I only know how to handle these files so far
list1={'bin_sgpbotx','bin_sgpbotz'};
list2={'bin_D','bin_U','bin_V','bin_W','bin_inv_Jac'};
list3={'bin_pressure','bin_vgrid'};
% ToDo: handle 'bin_volvec_full_domain' and 'bin_bndvls*'

% Get I and J
fid=fopen('startup_graphics','r'); IJ=textscan(fid,'%f',2); fclose(fid);
IJ=IJ{1}; I=IJ(1); J=IJ(2);

% Handle (I+1)x(1) files
for ii=1:length(list1)
    if strncmp(list1{ii},fname,length(list1{ii}))
        data1 = igwread_grabdata(fname,I+1,1);
    end
end

% Handle 2D (I)x(J+2) files
for ii=1:length(list2)
    if strncmp(list2{ii},fname,length(list2{ii}))
        data1 = igwread_grabdata(fname,I,J+2);
    end
end

% Handle 2D (I)x(J+2) files with two fields
for ii=1:length(list3)
    if strncmp(list3{ii},fname,length(list3{ii}))
        [data1, data2] = igwread_grabdata(fname,I,J+2);
    end
end

% Permute to [x,z] ordering - comment these lines if you want [z,x] order
data1=permute(data1,[2 1]);
if exist('data2','var'), data2=permute(data2,[2 1]); else data2=[]; end

end

function [d1,d2] = igwread_grabdata(fname,I,J)
% Helper function to load the contents of the file and reshape as requested
fid=fopen(fname,'rb'); dd=fread(fid,inf,'double'); fclose(fid);
if nargout == 1
    d1=reshape(dd,J,I);
else
    d1=reshape(dd(1:length(dd)/2)    ,J,I); % first field
    d2=reshape(dd(1+length(dd)/2:end),J,I); % second field
end
end
