function [timemat,posmat,trajmat,numderivs,cputime,cretime,ID,creator,destin,gdfver,gptver,destid,dummy1,dummy2] = ...
    gdf2mat(infile,outfile,gdfnamelen,gdfid,varargin)
% function [timemat,posmat,trajmat,numderivs,cputime,cretime,ID,creator,destin,gdfver,gptver,destid,dummy1,dummy2]
%              = gdf2mat(infile,outfile,gdfnamelen,gdfid,varargin)
%
%    GDF2MAT - Convert GPT binary data file into Matlab format
%
%    gdf2mat(infile,outfile)
%
%    gdf2mat(infile,outfile,gdfnamelen)
%
%    gdf2mat(infile,outfile,gdfnamelen,gdfid)
%
%    gdf2mat(infile,outfile,gdfnamelen,gdfid,'var1','var2','var3',...)
%
%    [timemat] = gdf2mat(...)
%
%    [timemat,posmat,trajmat] = gdf2mat(...)
%
%    [timemat,posmat,trajmat,numderivs,cputime] = gdf2mat(...)
%
%    [...,cretime,ID,creator,destin,gdfver,gptver,destid,dummy1,dummy2] = gdf2mat(...)
%
%    GDF2MAT is used to convert GPT particle data, from a GDF binary file,
%    into a format that Matlab can read.  To convert an ASCII data file,
%    use the command "GDFA2MAT" instead.
%
%    GDF2MAT can read both time/position slice data AND trajectory data,
%    created using the GPT function GDFTRANS.  Since a data file will only
%    contain one of these two kinds of data, arrays corresponding to the
%    other type will be empty.
%
%    gdf2mat(infile,outfile) - reads in GPT file INFILE, converts the
%    data into Matlab format, and saves it in the file OUTFILE.  Only
%    non-empty data variables (TIMEMAT, POSMAT and TRAJMAT - see below) are
%    saved.
%
%    gdf2mat(infile,outfile,gdfnamelen) - specify length of character
%    strings read from INFILE with GDFNAMELEN.  Default value is 
%    GDFNAMELEN = 16.
%
%    gdf2mat(infile,outfile,gdfnamelen,gdfid) - also specify GDF file
%    ID.  Default value is ID = 94325877.
%
%    Note: it is possible to leave any of the last 3 variables empty ie.
%    gdf2mat(infile,[],[],[]) is an acceptable syntax.  In this case, the
%    data will not be saved and the default values for GDFNAMELEN and GDFID
%    will be used.
%
%    [timemat] = gdf2mat(...) - reads in GPT file INFILE, converts the data
%    into Matlab format, and outputs the time data as the structure array
%    TIMEMAT.  It is possible to specify both OUTFILE and TIMEMAT if one
%    wishes to save the data and keep it in the Matlab workspace.
%
%    [timemat,posmat] = gdf2mat(...) - as above, but also outputs the
%    position data as the structure array POSMAT.
%
%    [timemat,posmat,trajmat] = gdf2mat(...) - as above, but also outputs
%    trajectory data as the structure array TRAJMAT.  This implies that the
%    GDF file being read contains trajectory data.
%
%    **Note** - by using this syntax: [timemat,posmat,trajmat] = gdf2mat,
%    either trajmat or timemat and posmat will be empty, since a GDF file
%    can either contain time/position slice data or trajectory data, but
%    not both.
%
%    [timemat,posmat,trajmat,numderivs,cputime] = gdf2mat(...) - return the
%    number of derivations required to produce the data, NUMDERIVS, and the
%    total CPU time required, CPUTIME.
%
%    [...,cretime,ID,creator,destin,gdfver,gptver,destid,dummy1,dummy2] = gdf2mat(...)
%    - read out the main header data at the start of the GDF file and
%    output each piece of data as a separate variable.  The output
%    variables are:
%
%        cretime - file creation date/time string
%        ID - GDF file ID; should be ID = 94325877 if the file is a valid
%            GDF file.
%        creator - app that created the file: should be GPT.
%        destin - destination.
%        gdfver - version of GDF-software.
%        gptver - version of creator (GPT).
%        destver - version of destination.
%        dummy1 - alignment (on 32-bit boundary).
%        dummy2 - empty dummy variable.
%
%    gdf2mat(infile,outfile,gdfnamelen,gdfid,'var1','var2','var3',...) -
%    select the variables to extract from the GDF file.  This allows a
%    limited amount of particle data to be extracted from a larger file.
%    For example, to only analyse X, Y and Z data, specify:
%
%        gdf2mat(infile,outfile,[],[],'x','y','z')
%
%    Each of the variables MUST be a string, and can take any one of the
%    following values (direct from the GDF file): 
%
%        'x', 'y', 'z', 'Bx', 'By', 'Bz', 'G', 'rxy', 'fEx', 'fEy', 'fEz',
%        'fBx', 'fBy', 'fBz', 'm', 'q', 'nmacro', 'rmacro'.
%
%    The following variables can be specified, but are redundant as they
%    are included automatically:
%
%        'ID', 'time', 'position'.
%
%    In addition, since GPT does not output x', y' or Energy, specify any
%    of the following variables to calculate them:
%
%        'xp', 'yp', 'E'.
%
%    To extract EVERY particle coordinate, use the variable 'gpt':
%
%        gdf2mat(infile,outfile,[],[],'gpt')
%
%    To extract all the particle data, plus calculating x', y' and particle
%    energy, use the variable 'all':
%
%        gdf2mat(infile,outfile,[],[],'all')
%
%    To extract only the most commonly used particle data - x, y, z, rxy,
%    time, position, ID, x', y' and energy - use the variable 'simp':
%
%        gdf2mat(infile,outfile,[],[],'simp')
%

% File is the property of Simon Jolly, Imperial College London.  File
% released under the GNU public license.
%
% File history
%
%   01-Apr-2008 S. Jolly
%       Full read-in of GPT .gdf file and conversion to Matlab format.
%
%   21-Apr-2008 S. Jolly
%       Added ability to specify which variables are required using
%       'varargin'.
%
%   26-Jun-2008 S. Jolly
%       Added ability to read in trajectory files created using GDFTRANS
%       and analysed files created using GDFA.
%
%======================================================================

    if nargin < 1
        error('Too few input variables: syntax is [timemat, posmat] = gdf2mat(infile)') ;
    elseif nargin > 24
        error('Too many input variables: syntax is [timemat, posmat] = gdf2mat(infile, outfile, gdfnamelen, gdfid, varargin)') ;
    elseif nargout > 14
        error('GDF2MAT:tooManyNargout', ...
            ['Too many output variables: syntax is \n' ...
            '[timemat,posmat,trajmat,numderivs,cputime,cretime,ID,creator,destin,' ...
            'gdfver,gptver,destid,dummy1,dummy2] \n' ...
            '  = gdf2mat(infile,outfile,gdfnamelen,gdfid)']) ;
    elseif ~ischar(infile)
        error('Input filename must be a string') ;
    elseif ( nargin > 1 ) && ~isempty(outfile) && ~ischar(outfile)
        error('Output filename must be a string') ;
    elseif ( exist(infile,'file') ~= 2 )
        error('Input file cannot be found')
    elseif ~strcmp( infile(end-2:end),'gdf')
        error('Input file is not a GDF file')
    end

    if nargin < 4 || isempty(gdfid)
        gdfid = 94325877 ;
    end

    if nargin < 3 || isempty(gdfnamelen)
        gdfnamelen = 16 ;
    end

%    Select particle coordinates to read out

    standpartset = {'x', 'y', 'z', 'Bx', 'By', 'Bz', 'G', 'rxy', 'fEx', ...
        'fEy', 'fEz', 'fBx', 'fBy', 'fBz', 'm', 'q', 'nmacro', 'rmacro', ...
        'ID', 'time', 'position'} ;

    fullpartset = {'x', 'y', 'z', 'Bx', 'By', 'Bz', 'G', 'rxy', 'fEx', ...
        'fEy', 'fEz', 'fBx', 'fBy', 'fBz', 'm', 'q', 'nmacro', 'rmacro', ...
        'ID', 'time', 'position', 'xp', 'yp', 'E'} ;

    simppartset = {'x', 'y', 'z', 'rxy', 'ID', 'time', 'position', 'xp', 'yp', 'E'} ;

    xpcalc = 0 ; ypcalc = 0 ; Ecalc = 0 ;

    partset = {} ;

    if nargin < 5

        partset = standpartset ;

    else

        if sum(strcmp(varargin,'gpt')) > 0
            partset = standpartset ;
        elseif sum(strcmp(varargin,'all')) > 0
            partset = fullpartset ;
            xpcalc = 1 ; ypcalc = 1 ; Ecalc = 1 ;
        elseif sum(strcmp(varargin,'simp')) > 0
            partset = simppartset ;
            xpcalc = 1 ; ypcalc = 1 ; Ecalc = 1 ;
        else

            for i = 1:length(varargin)
                if sum(strcmp(fullpartset,varargin{i})) == 0
                    error('GDF2MAT:unrecognisedVarargin', ...
                        ['Input variable ' varargin{i} ' not recognised: see the help file']) ;
                elseif sum(strcmp(varargin{i},'xp')) > 0
                    xpcalc = 1 ;
                elseif sum(strcmp(varargin{i},'yp')) > 0
                    ypcalc = 1 ;
                elseif sum(strcmp(varargin{i},'E')) > 0
                    Ecalc = 1 ;
                else
                    partset{end+1} = varargin{i} ;
                end
            end

            partset{end+1} = ('ID') ;
            partset{end+1} = ('time') ;
            partset{end+1} = ('position') ;

        end

    end

    Bxloc = [] ; Byloc = [] ; Bzloc = [] ;
    mloc = [] ; Gloc = [] ;

%    Initialise variables

    timemat = [] ;
    posmat = [] ;
    trajmat = [] ;
    tempmat = [] ;
    numderivs = 0 ;
    cputime = 0 ;
    charge = 1.602176487e-19 ;
    c = 299792458 ;

%    Open file and read in main data header

    fid = fopen(infile,'r') ;
    fseekstat = fseek(fid,0,'eof') ;
    endpos = ftell(fid) ;
    fseekstat = fseek(fid,0,'bof') ;

    wbh = waitbar(0,'Loading data slice... ','Name','GDF File Progress') ;

    tic ;

    [ID,cretime,creator,destin,gdfver,gptver,destid,dummy1,dummy2] = gdfreadmainhead(fid,gdfnamelen,gdfid) ;

    timeloop = 0 ;
    posloop = 0 ;
    trajloop = 0 ;
    endoffile = 0 ;

    while endoffile == 0

%    Find out how large the file is to set the length of the waitbar

        curpos = ftell(fid) ;
        percfile = curpos./endpos ;
        if trajloop > timeloop
            wbh = waitbar(percfile,wbh,['Loading trajectory slice ' num2str(trajloop+1) ' (' num2str(percfile.*100,'%12.0f') '%)']) ;
        else
            wbh = waitbar(percfile,wbh,['Loading timeslice ' num2str(timeloop+1) ' (' num2str(percfile.*100,'%12.0f') '%)']) ;
        end

%    Read in data header and first data point

        [datname, datlength, readstring, extras, blocktype, dirstart, dirend, endoffile] = gdfreadhead(fid) ;

        if endoffile
            break
        end

        eval([datname ' = fread(fid, datlength, readstring)' extras ' ;']) ;

        if sum(strcmp(creator,{'GDFA';'ASCI2GDF'})) && ~dirstart
            eval(['tempmat.' datname ' = ' datname ' ;']) ;
        end

%    If this is the start of a directory entry, choose the relevant
%    structure

        if dirstart
            switch datname(1:2)
                case 'ti'
                    fillarr = 'timemat' ;
                    fillloop = 'timeloop' ;
                    timeloop = timeloop + 1 ;
                case 'po'
                    fillarr = 'posmat' ;
                    fillloop = 'posloop' ;
                    posloop = posloop + 1 ;
                case 'ID'
                    fillarr = 'trajmat' ;
                    fillloop = 'trajloop' ;
                    trajloop = trajloop + 1 ;
                otherwise
                    fillarr = [] ;
                    fillloop = [] ;
                    warning('GDF2MAT:unknownDirType', ...
                        ['Unknown directory type ' datname '...']) ;
            end

            eval([fillarr '(' fillloop ').' datname ' = ' datname ' ;']) ;

            dirend = 0 ;

%    Loop through till the end of the directory entry, and read out the
%    directory to the relevant structure

            while dirend == 0

                [datname, datlength, readstring, extras, blocktype, dirstart, dirend, endoffile] = gdfreadhead(fid) ;

                if dirend
                    break
                elseif sum(strcmp(partset,datname)) > 0 || strcmp(creator,'GDFA')
                    eval([fillarr '(' fillloop ').' datname ' = fread(fid, datlength, readstring)' extras ' ;']) ;
                else
                    eval([datname ' = fread(fid, datlength, readstring)' extras ' ;']) ;
                end

            end

%    Set up some variables that allow xp, yp and E to be calculated
%    correctly

            if dirend

                if sum(strcmp(partset,'Bx')) > 0
                    Bxloc = [ fillarr '(' fillloop ').' ] ;
                end

                if sum(strcmp(partset,'By')) > 0
                    Byloc = [ fillarr '(' fillloop ').' ] ;
                end

                if sum(strcmp(partset,'Bz')) > 0
                    Bzloc = [ fillarr '(' fillloop ').' ] ;
                end

                if sum(strcmp(partset,'m')) > 0
                    mloc = [ fillarr '(' fillloop ').' ] ;
                end

                if sum(strcmp(partset,'G')) > 0
                    Gloc = [ fillarr '(' fillloop ').' ] ;
                end

%    If required, calculate xp, yp and E for current time/position slice

                eval(['Bxexists = isfield(' fillarr ', ''Bx'') ;']) ;
                eval(['Byexists = isfield(' fillarr ', ''By'') ;']) ;
                eval(['Bzexists = isfield(' fillarr ', ''Bz'') ;']) ;
                eval(['mexists = isfield(' fillarr ', ''m'') ;']) ;

                if xpcalc && Bxexists && Bzexists
                    eval([fillarr '(' fillloop ').xp = atan2(' Bxloc 'Bx,' ...
                        Bzloc 'Bz) ;']) ;
                end

                if ypcalc && Bxexists && Bzexists
                    eval([fillarr '(' fillloop ').yp = atan2(' Byloc 'By,' ...
                        Bzloc 'Bz) ;']) ;
                end

                if Ecalc && mexists
                    eval([fillarr '(' fillloop ').E = (' mloc 'm.*(' Gloc 'G - 1).*(c.^2))./charge ;']) ;
                end

            end

        end

    end

%    Close file and waitbar

    fclose(fid) ;
    close(wbh) ;

%    If the GDF file was created using GDFA or ASCI2GDF, some arrays may
%    not be created correctly, so reassign arrays

    if isempty(timemat) && isempty(posmat) && ~isempty(tempmat)
        if xpcalc && isfield(tempmat,'Bx') && isfield(tempmat,'Bz')
            tempmat.xp = atan2(tempmat.Bx,tempmat.Bz) ;
        end
        if ypcalc && isfield(tempmat,'By') && isfield(tempmat,'Bz')
            tempmat.yp = atan2(tempmat.By,tempmat.Bz) ;
        end
        if Ecalc && isfield(tempmat,'m') && isfield(tempmat,'G')
            tempmat.E = (tempmat.m.*(tempmat.G - 1).*(c.^2))./charge ;
        end
    end

    if strcmp(creator,'GDFA')
        if isempty(timemat) && isfield(tempmat, 'time')
            timemat = tempmat ;
        elseif isempty(posmat) && isfield(tempmat, 'position') && ~isfield(tempmat, 'time')
            posmat = tempmat ;
        end
    elseif strcmp(creator,'ASCI2GDF')
        if isempty(posmat)
            posmat = tempmat ;
        elseif isempty(timemat) && ~isfield(tempmat, 'position')
            timemat = tempmat ;
        end
    end

%    If there are less than 4 output arguments, display the CPUTIME and
%    NUMDERIVS data

    if nargout < 4
        if exist('cputime','var') && ( length(cputime) == 1 )
            disp(['CPU Time   : ' num2str(cputime,'%0.6f')]) ;
        end
        if exist('numderivs','var') && ( length(numderivs) == 1 )
            disp(['Num. Derivs: ' num2str(numderivs,'%0.0f')]) ;
        end
    end

    loadtime = toc ;
    disp(['Load time  : ' num2str(loadtime,'%0.6f') ' secs']) ;

%    Save data to .mat-file

    if nargin > 1 && ~isempty(outfile)

        savecom = ['save ' outfile ' '] ;

        if ~isempty(timemat)
            savecom = [savecom 'timemat '] ;
        end
        if ~isempty(posmat)
            savecom = [savecom 'posmat '] ;
        end
        if ~isempty(trajmat)
            savecom = [savecom 'trajmat '] ;
        end

        savecom = [savecom '-V7.3'] ;

        tic ;
        disp(['Saving data file ' outfile ' ...']) ;
        eval(savecom) ;
        savetime = toc ;
        disp(['Save time  : ' num2str(savetime,'%0.6f') ' secs']) ;
    end

%    Clear variables if they're not needed

    if nargout == 0
        clear variables ;
    end

    return
