function [trajmat,numderivs,cputime,cretime,ID,creator,destin,gdfver,gptver,destid,dummy1,dummy2] = ...
    gdf2trajmat(infile,outfile,gdfnamelen,gdfid,varargin)
% function [trajmat,numderivs,cputime,cretime,ID,creator,destin,gdfver,gptver,destid,dummy1,dummy2]
%              = gdf2trajmat(infile,outfile,gdfnamelen,gdfid,varargin)
%
%    GDF2TRAJMAT - Convert GPT binary data file into Matlab trajectory format
%
%    gdf2trajmat(infile,outfile)
%
%    gdf2trajmat(infile,outfile,gdfnamelen)
%
%    gdf2trajmat(infile,outfile,gdfnamelen,gdfid)
%
%    gdf2trajmat(infile,outfile,gdfnamelen,gdfid,'var1','var2','var3',...)
%
%    [trajmat] = gdf2trajmat(...)
%
%    [trajmat,numderivs,cputime] = gdf2trajmat(...)
%
%    [...,cretime,ID,creator,destin,gdfver,gptver,destid,dummy1,dummy2] = gdf2trajmat(...)
%
%    GDF2TRAJMAT is used to convert GPT particle data, from a GDF binary
%    file, into trajectory data that Matlab can read.  This requires
%    loading the binary data stored in time slices from the GDF file and
%    transposing the data into trajectory format.  To convert GPT data
%    into timeslice/screen data, without the matrix transpose required to
%    produce the trajectory data, use "GDF2MAT" instead.
%
%    gdf2trajmat(infile,outfile) - reads in GPT file INFILE, converts the
%    data into trajectory data in Matlab format, and saves it in the file
%    OUTFILE as the structure array TRAJMAT.  If the file OUTFILE already
%    exists, and contains partially written particle data, GDF2TRAJMAT
%    carries on from where it left off.  This works around Matlabs memory
%    limits.
%
%    gdf2trajmat(infile,outfile,gdfnamelen) - specify length of character
%    strings read from INFILE with GDFNAMELEN.  Default value is 
%    GDFNAMELEN = 16.
%
%    gdf2trajmat(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.
%    gdf2trajmat(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.
%
%    [trajmat] = gdf2trajmat(...) - reads in GPT file INFILE, converts the
%    data into Matlab format, tranposes the data to produce the trajectory
%    data and outputs the data as the structure array TRAJMAT.  It is
%    possible to specify both OUTFILE and TRAJMAT if one wishes to save the
%    data and keep it in the Matlab workspace.
%
%    [trajmat,numderivs,cputime] = gdf2trajmat(...) - 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] = gdf2trajmat(...)
%    - 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.
%
%    gdf2trajmat(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:
%
%        gdf2trajmat(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', 'time'.
%
%    The following variables should not be specified since they are
%    redundant (position data is not used to calculate trajectory data, and
%    particle ID specifies the array index TRAJMAT(ID)):
%
%        'ID', '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':
%
%        gdf2trajmat(infile,outfile,[],[],'gpt')
%
%    To extract all the particle data, plus calculating x', y' and particle
%    energy, use the variable 'all':
%
%        gdf2trajmat(infile,outfile,[],[],'all')
%
%    To extract only the most commonly used particle data - x, y, z, rxy,
%    time, x', y' and energy - use the variable 'simp':
%
%        gdf2trajmat(infile,outfile,[],[],'simp')
%

    if nargin < 1
        error('Too few input variables: syntax is [trajmat] = gdf2trajmat(infile)') ;
    elseif nargin > 24
        error('Too many input variables: syntax is [trajmat] = gdf2trajmat(infile, outfile, gdfnamelen, gdfid, varargin)') ;
    elseif nargout > 12
        error('GDF2TRAJMAT:tooManyNargout', ...
            ['Too many output variables: syntax is \n' ...
            '[trajmat,numderivs,cputime,cretime,ID,creator,destin,' ...
            'gdfver,gptver,destid,dummy1,dummy2] \n' ...
            '  = gdf2trajmat(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', ...
        'time'} ;

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

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

    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

            partset = varargin ;

            for i = 1:length(varargin)
                if sum(strcmp(fullpartset,varargin{i})) == 0
                    error('GDF2TRAJMAT: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 ;
                end
            end

        end

    end

    for i = 1:length(partset)
        eval(['trajmat.' partset{i} ' = [] ;']) ;
    end


%    Initialise variables

    trajmat = [] ;
    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 timeslices... ','Name','GDF File Progress') ;

%    tic ;
    if nargout < 5
        gdfreadmainhead(fid,gdfnamelen,gdfid) ;
    else
        [ID,cretime,creator,destin,gdfver,gptver,destid,dummy1,dummy2] = gdfreadmainhead(fid,gdfnamelen,gdfid) ;
    end

    timeloop = 0 ;
    trajinc = 0 ;
    endoffile = 0 ;
    loadtime = 0 ;
    calctime = 0 ;
    transtime = 0 ;

    while endoffile == 0

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

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

%    Read in data header and first data point

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

        if endoffile
            break
        end

        eval([datname ' = fread(fid, datlength, readstring)' extras ' ;']) ;
%        loadtime = loadtime + toc ;

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

        if dirstart

            if datname(1) == 't'
                trajinc = 1 ;
                timeloop = timeloop + 1 ;
            elseif datname(1) == 'p'
                trajinc = 0 ;
            else
                trajinc = 0 ;
                warning('GDF2TRAJMAT:unknownDirType', ...
                    ['Unknown directory type ' datname '...']) ;
            end

            dirend = 0 ;

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

            while dirend == 0

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

                if dirend
                    break
                else
                    eval([datname ' = fread(fid, datlength, readstring)' extras ' ;']) ;
                end
                loadtime = loadtime + toc ;

            end

            if dirend && trajinc

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

                tic ;
                if xpcalc
                    xp = atan2(Bx,Bz) ;
                end

                if ypcalc
                    yp = atan2(By,Bz) ;
                end

                if Ecalc
                    E = (m.*(G - 1).*(c.^2))./charge ;
                end
                calctime = calctime + toc ;

%    Transpose timeslice data into trajectory data

                tic ;
                for j = 1:length(ID)

                    partid = ID(j) ;

                    if length(trajmat) < partid || isempty(trajmat(partid).x)
                        for k = 1:length(partset)
                            if strcmp(partset{k},'time')
                                eval(['trajmat(partid).' partset{k} ' = ' partset{k} ' ;']) ;
                            else
                                eval(['trajmat(partid).' partset{k} ' = ' partset{k} '(j) ;']) ;
                            end
                        end
                    else
                        for k = 1:length(partset)
                            if strcmp(partset{k},'time')
                                eval(['trajmat(partid).' partset{k} '(end+1) = ' partset{k} ' ;']) ;
                            else
                                eval(['trajmat(partid).' partset{k} '(end+1) = ' partset{k} '(j) ;']) ;
                            end
                        end
                    end

                end
                transtime = transtime + toc ;

            end

        end

    end

%    Close file and waitbar

    fclose(fid) ;
    close(wbh) ;

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

    if nargout < 3
        disp(['CPU Time   : ' num2str(cputime,'%0.6f')]) ;
        disp(['Num. Derivs: ' num2str(numderivs,'%0.0f')]) ;
    end

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

%    Save data to .mat-file

    if nargin > 1 && ~isempty(outfile)
        tic ;
        disp(['Saving data file ' outfile ' ...']) ;
        save(outfile,'trajmat','-v7.3') ;
        toc ;
    end

%    Clear variables if they're not needed

    if nargout == 0
        clear variables ;
    end

    return
