function [ header , dataheader, Data , Times ] = readflxhst( filename )
%READFLXHST Function designed to read PZFlex History files into MATLAB
%
%   [ header , dataheader, Data , Times ] = readflxhst( filename )
%
%   This code runs through the data array twice, initially to parse the
%   FORTRAN records to ensure the integrity of the data, and then to read
%   in the data. The function automatically determines the difference
%   between 32-bit and 64-bit data. Additionally, calling the function
%   without an input argument will open a dialogue box which prompts the
%   user to select a file.
%
%   The header structure contains general information about the
%   file as well as the time record header information. The dataheader
%   structure contains the specific information about each data record. The
%   Times array is an (n x 1) element array where the n is the number of time
%   steps and the Data array is an (n x m) array where m is the number of
%   data records.
%
%   The fields in the header structure are defined as:
%
%   lhedr - (C) 4-character label = 'hedr'
%   itype - (I) header type = 2
%   ftype - (C) 20-character variable defining file type = 'hist'
%   fform - (C) 20-character variable defining file form = 'bin'
%   date - (C) 80-character variable defining date file was written
%   code - (C) 20-character variable identifying the code = 'flex'
%   user - (C) 240-character variable defining user name
%   tag - (C) 280-character variable containing title and tag
%   iversion - (I) file type version number
%   nparts - (I) number of partitions for the model
%   ipart - (I) partition number for this file
%   intbyt - (I) number of bytes for standard integer
%   irlbyt - (I) number of bytes for standard real number
%   nch1 - (I) number of characters in string 1 variables = 20
%   nch2 - (I) number of characters in string 2 variables = 80
%   nch3 - (I) number of characters in string 3 variables = 200
%   iextra1 - (I) unused extra integer 1
%   iextra2 - (I) unused extra integer 2
%   int12 - (I) integer value = 12
%   nrecd - (I) number of data records on the file including the time record
%   ntim - (I) number of time values for which time history data was saved
%   ngrd - (I) number of 20-character words for grid name
%   ndes - (I) number of 20-character words for description
%   ntag - (I) number of 20-character words for curve tag
%   nshift - (I) 1 = time shift for each record; 0 = no timeshift
%   nlabl - (I) number of 20-character words for curve lable
%   nindx - (I) number of indicies for each record
%   ncord - (I) number of coordinates for each record
%   nchr - (I) number of bytes per character word = 20
%   nextra - (I) unused integer
%   label - (C) 140-characer label and description for time record
%   zero - (F) the number 0.0
%   tstep - (F) time step of calculation
%   i0 - (I) integer = 0
%   tbegin - (F) start time = time(1)
%   tshift - (F) time shift to add to time
%   tend - (F) end time = time(ntim)
%   tstart - (F) start time
%
%
%   Each individual data record has its own header sturcture which contains
%   the following fields:
%
%   label - (C) 20-characer label for data record
%   ltag - (C) 20-character curve tag
%   lgrid - (C) 20-character grid tag
%   ldesc - (C) 80-character description for data record
%   xcrd - (F) x-coordinates of this time history output
%   ycrd - (F) y-coordinates of this time history output
%   zcrd - (F) z-coordinates of this time history output
%   i - (I) i index for this time history record
%   j - (I) j index for this time history record
%   k - (I) k index for this time history record
%   datamin - (F) minimum data value
%   datamax - (F) maximum data value
%
%   Known Issues: The output array Data will always be a double precision
%   array as the preallocation is not varied depending on the bit
%   resolution of the file. Although this doubles the memory usage for the
%   variable uneccesarily it is considered that this is unlikely to present
%   a problem, while using an if statement could negatively impact the
%   speed of the code.
%
%

if nargin == 0
    
    [filename,pathname,FilterIndex] = uigetfile('*.flxhst','Please select a data file');
    
    if FilterIndex ~= 1
        
        error('Please select a single data file')
        
    end
    filename = [pathname filename];    
end

fid = fopen([filename]); % Opens the input file
 
% Loop checks the integrity of the file by inspecting the FORTRAN records
while feof(fid) ~= 1
    
    record1 = fread(fid,1,'*int32');
    
    if feof(fid) % If the end of the file is reached the integrity check is ended
        
        break
        
    end
    
    fseek(fid,record1,0);
    
    record2 = fread(fid,1,'*int32');
    
    if record1 ~= record2
        
        error('The file cannot be read and may be corrupted') % If the FORTRAN records do not match an error is returned
        
    end

end

frewind(fid); % The file is returned to the start to allow data to be read from it

bitrestest = fread(fid,1,'*int32'); % The first FORTRAN record is used to determine the bit resolution of the file

if bitrestest == 8
    
    % A value of 8 bytes indicates that the file is 32 bit
    inttype = '*int32';
    floattype = '*single';
    
elseif bitrestest == 12
    
    % A value of 12 bytes indicates that the file is 64 bit
    inttype = '*int64';
    floattype = '*double';
    
else
    
    % Any other value indicates that something is wrong
    error('The initial FORTRAN record is not an expected value. The file may be corrupted')
    
end

% The 1st header record is read in

header.lhedr = fread(fid,4,'*char');

header.itype = fread(fid,1,inttype);

fseek(fid,8,0);

% 2nd header record

header.ftype = fread(fid,20,'*char');

header.fform = fread(fid,20,'*char');

header.date = fread(fid,80,'*char');

header.code = fread(fid,20,'*char');

header.user = fread(fid,240,'*char');

fseek(fid,8,0);

% 3rd header record

header.tag = fread(fid,280,'*char');

fseek(fid,8,0);

% 4th header record

header.iversion = fread(fid,1,inttype);

header.nparts = fread(fid,1,inttype);

header.ipart = fread(fid,1,inttype);

header.intbyt = fread(fid,1,inttype);

header.irlbyt = fread(fid,1,inttype);

header.nch1 = fread(fid,1,inttype);

header.nch2 = fread(fid,1,inttype);

header.nch3 = fread(fid,1,inttype);

header.iextra1 = fread(fid,1,inttype);

header.iextra2 = fread(fid,1,inttype);

fseek(fid,8,0);

% Base record

header.int12 = fread(fid,1,inttype);

header.nrecd = fread(fid,1,inttype);

header.ntim = fread(fid,1,inttype);

header.ngrd = fread(fid,1,inttype);

header.ndes = fread(fid,1,inttype);

header.ntag = fread(fid,1,inttype);

header.nshift = fread(fid,1,inttype);

header.nlabl = fread(fid,1,inttype);

header.nindx = fread(fid,1,inttype);

header.ncord = fread(fid,1,inttype);

header.nchr = fread(fid,1,inttype);

header.nextra = fread(fid,1,inttype);

fseek(fid,8,0);

% Time header record rand

header.label = fread(fid,140,'*char');

fseek(fid,8,0);

header.zero1 = fread(fid,1,floattype);

header.zero2 = fread(fid,1,floattype);

header.tstep = fread(fid,1,floattype);

header.i0 = fread(fid,1,floattype);

header.zero3 = fread(fid,1,floattype);

header.tbegin = fread(fid,1,floattype);

header.tshift = fread(fid,1,floattype);

header.tend = fread(fid,1,floattype);

header.tstart = fread(fid,1,floattype);

% Time record

Times = fread(fid,header.ntim,floattype);

fseek(fid,8,0);

% Data record is preallocated

Data = zeros(header.ntim,header.nrecd - 1);

% 1st data header

% The header describes the first record but is put in the final header
% record in order to preallocate the structure array. It is then moved back
% to the first header record in order to preserve tho order of the headers

dataheader(header.nrecd - 1).label = fread(fid,20,'*char');
dataheader(1).label = dataheader(header.nrecd - 1).label;

dataheader(header.nrecd - 1).ltag = fread(fid,20,'*char');
dataheader(1).ltag = dataheader(header.nrecd - 1).ltag;

dataheader(header.nrecd - 1).lgrid = fread(fid,20,'*char');
dataheader(1).lgrid = dataheader(header.nrecd - 1).lgrid;

dataheader(header.nrecd - 1).ldesc = fread(fid,80,'*char');
dataheader(1).ldesc = dataheader(header.nrecd - 1).ldesc;

fseek(fid,8,0);

dataheader(header.nrecd - 1).xcrd = fread(fid,1,floattype);
dataheader(1).xcrd = dataheader(header.nrecd - 1).xcrd;

dataheader(header.nrecd - 1).ycrd = fread(fid,1,floattype);
dataheader(1).ycrd = dataheader(header.nrecd - 1).ycrd;

dataheader(header.nrecd - 1).zcrd = fread(fid,1,floattype);
dataheader(1).zcrd = dataheader(header.nrecd - 1).zcrd;

dataheader(header.nrecd - 1).icrd = fread(fid,1,inttype);
dataheader(1).icrd = dataheader(header.nrecd - 1).icrd;

dataheader(header.nrecd - 1).jcrd = fread(fid,1,inttype);
dataheader(1).jcrd = dataheader(header.nrecd - 1).jcrd;

dataheader(header.nrecd - 1).kcrd = fread(fid,1,inttype);
dataheader(1).kcrd = dataheader(header.nrecd - 1).kcrd;

dataheader(header.nrecd - 1).datamin = fread(fid,1,floattype);
dataheader(1).datamin = dataheader(header.nrecd - 1).datamin;

dataheader(header.nrecd - 1).datamax = fread(fid,1,floattype);
dataheader(1).datamax = dataheader(header.nrecd - 1).datamax;

dataheader(header.nrecd - 1).tshift = fread(fid,1,floattype);
dataheader(1).tshift = dataheader(header.nrecd - 1).tshift;

% 1st data record

Data(:,1) = fread(fid,header.ntim,floattype);

% The loop continues to read in data headers and records
for count = 2:header.nrecd - 1
       
    fseek(fid,8,0);
    
    % Data header
    
    dataheader(count).label = fread(fid,20,'*char');
    
    dataheader(count).ltag = fread(fid,20,'*char');
    
    dataheader(count).lgrid = fread(fid,20,'*char');
    
    dataheader(count).ldesc = fread(fid,80,'*char');
    
    fseek(fid,8,0);
    
    dataheader(count).xcrd = fread(fid,1,floattype);
    
    dataheader(count).ycrd = fread(fid,1,floattype);
    
    dataheader(count).zcrd = fread(fid,1,floattype);
    
    dataheader(count).icrd = fread(fid,1,inttype);
    
    dataheader(count).jcrd = fread(fid,1,inttype);
    
    dataheader(count).kcrd = fread(fid,1,inttype);
    
    dataheader(count).datamin = fread(fid,1,floattype);
    
    dataheader(count).datamax = fread(fid,1,floattype);
    
    dataheader(count).tshift = fread(fid,1,floattype);
    
    % Data record
    
    Data(:,count) = fread(fid,header.ntim,floattype);
    
end

fclose(fid);

end
