function [timemat, posmat, trajmat] = gdfa2mat(infile, outfile, trajfile)
%
% function [timemat, posmat, trajmat] = gdfa2mat(infile, outfile, trajfile)
%
%    GDFA2MAT - Convert GPT ASCII data file into Matlab format
%
%    GDFA2MAT is used to convert GPT particle data, from a GDF ASCII file,
%    into a format that Matlab can read.  The input file must have been
%    created with the GPT command 'gdf2a' and must be in ASCII format.
%
%    gdfa2mat(infile, outfile) - reads in GPT file INFILE, converts the
%    data into Matlab format, and saves it in the file OUTFILE.  If the
%    file OUTFILE already exists, and contains partially written particle
%    data, GDFA2MAT carries on from where it left off.  This works
%    around Matlabs memory limits.
%
%    gdfa2mat(infile, outfile, trajfile) - as above, but also calculates the
%    trajectory data (see GDFTRANS) and saves it to the .MAT-file TRAJFILE
%    as the structure array TRAJMAT.
%
%    [timemat] = gdfa2mat(infile) - reads in GPT file INFILE, converts
%    the data into Matlab format, and outputs the time data as the
%    structure array TIMEMAT.
%
%    [timemat, posmat] = gdfa2mat(infile) - reads in GPT file INFILE,
%    converts the data into Matlab format, and outputs the time data as the
%    structure array TIMEMAT and the position data as the structure array
%    POSMAT.
%
%    [timemat, posmat, trajmat] = gdfa2mat(infile) - as above, but
%    also calculates the trajectory data (see GDFTRANS) and outputs it as
%    the structure array TRAJMAT.

    if nargin < 1
        error('Too few input variables: syntax is [timemat, posmat, trajmat] = gdfa2mat(infile, outfile, trajfile)') ;
    elseif nargin > 3
        error('Too many input variables: syntax is [timemat, posmat, trajmat] = gdfa2mat(infile, outfile, trajfile)') ;
    elseif nargout > 3
        error('Too many output variables: syntax is [timemat, posmat, trajmat] = gdfa2mat(infile, outfile, trajfile)') ;
    elseif ~ischar(infile)
        error('Input filename must be a string') ;
    elseif ( nargin > 1 ) && ~ischar(outfile)
        error('Output filename must be a string') ;
    elseif ( nargin > 2 ) && ~ischar(trajfile)
        error('Trajectory filename must be a string') ;
    elseif ( exist(infile,'file') ~= 2 )
        error('Input file cannot be found')
    end

%    Set up initial variables

    partnum = [] ;
    type = [] ;
    colheads = [] ;
    goodline = 0 ;
    timeloop = 0 ;
    posloop = 0 ;
    timemat = [] ;
    posmat = [] ;
    time = [] ;
    position = [] ;

%    Set up 'trajmat' array for trajectory data if required

    if nargin > 2 || nargout > 2
        trajmat.x = [] ;
        trajmat.xp = [] ;
        trajmat.y = [] ;
        trajmat.yp = [] ;
        trajmat.z = [] ;
        trajmat.Bx = [] ;
        trajmat.By = [] ;
        trajmat.Bz = [] ;
        trajmat.rxy = [] ;
        trajmat.time = [] ;
    end

%    Open text file first time and read in header data

    fid = fopen(infile) ;

    logdat = textscan(fid, '%s',1,'delimiter','EOL') ;
    numderivstext = textscan(fid, '%s %n',1) ;
    cputext = textscan(fid, '%s %n',1) ;

    type = textscan(fid, '%s %n',1) ;

%    Read in column headings for first data block to find out how many columns there are

    colheads = textscan(fid, '%s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s',1) ;

    for i = 1:length(colheads)

        if ~isempty(char(colheads{i}))
            colhead{i} = colheads{i} ;
        else
            break
        end

    end

%    Create arrays to read in data using correct number of column entries

    colheadp = colhead ;
    colheadp{end+1} = 'ID' ;

    numcols = length(colhead) ;
    hstrt = [] ;
    nstrt = [] ;

    for i = 1:numcols

        hstrt = [hstrt '%s '] ;
        nstrt = [nstrt '%n '] ;

    end

    hstrt(end) = [] ;
    nstrt(end) = [] ;
    hstrp = [hstrt ' %s'] ;
    nstrp = [nstrt ' %n'] ;

    fclose(fid) ;

%    Open file a second time to read in actual data

    fid = fopen(infile) ;

    disp(['fid = ' num2str(fid)]) ;

%    Read in header data

    logdat = textscan(fid, '%s',1,'delimiter','EOL') ;
    numderivstext = textscan(fid, '%s %n',1) ;
    cputext = textscan(fid, '%s %n',1) ;

%    Loop through time/position data entries and read in data from file

    while goodline == 0

%    Work out whether data is time or position

        [type, filepos] = textscan(fid, '%s %n',1) ;
        typesize = size(type{2}) ;
        typetext = char(type{1}) ;

%    If data is empty, end of file reached

        if typesize(1) == 0

            goodline = 1 ;
            break

%    Read in time data

        elseif typetext(1) == 't'

            goodline = 0 ;
            timeloop = timeloop + 1 ;
            [headings, filepos] = textscan(fid, hstrt, 1) ;
            [dat, filepos] = textscan(fid, nstrt) ;
            datsize = size(dat) ;

            if ( rem(timeloop,50) == 0 )
%                disp(['Timestep ' num2str(timeloop) ' ... ']) ;
            end

%    Save time data to 'timemat' structure array

            for j = 1:datsize(2)
%                timemat(:,j,timeloop) = dat{j} ;
                eval(['timemat(' num2str(timeloop) ').' char(colhead{j}) ' = dat{j} ;']) ;
            end

%    Calculate angles and store timeslice data in array

            eval(['timemat(' num2str(timeloop) ').xp = atan2(timemat(' num2str(timeloop) ...
                ').Bx,timemat(' num2str(timeloop) ').Bz) ;']) ;
            eval(['timemat(' num2str(timeloop) ').yp = atan2(timemat(' num2str(timeloop) ...
                ').By,timemat(' num2str(timeloop) ').Bz) ;']) ;
            eval(['timemat(' num2str(timeloop) ').time = type{2} ;']) ;

%    Convert timeslice data to particle trajectories, one particle at a time

            if nargin > 2 || nargout > 2

                partids = timemat(timeloop).ID ;

                for j = 1:length(partids)

                    partid = partids(j) ;

                    if length(trajmat) < partid || isempty(trajmat(partid).x)
                        trajmat(partid).x = timemat(timeloop).x(j) ;
                        trajmat(partid).xp = timemat(timeloop).xp(j) ;
                        trajmat(partid).y = timemat(timeloop).y(j) ;
                        trajmat(partid).yp = timemat(timeloop).yp(j) ;
                        trajmat(partid).z = timemat(timeloop).z(j) ;
                        trajmat(partid).Bx = timemat(timeloop).Bx(j) ;
                        trajmat(partid).By = timemat(timeloop).By(j) ;
                        trajmat(partid).Bz = timemat(timeloop).Bz(j) ;
                        trajmat(partid).rxy = timemat(timeloop).rxy(j) ;
                        trajmat(partid).time = timemat(timeloop).time ;
                    else
                        trajmat(partid).x(end+1) = timemat(timeloop).x(j) ;
                        trajmat(partid).xp(end+1) = timemat(timeloop).xp(j) ;
                        trajmat(partid).y(end+1) = timemat(timeloop).y(j) ;
                        trajmat(partid).yp(end+1) = timemat(timeloop).yp(j) ;
                        trajmat(partid).z(end+1) = timemat(timeloop).z(j) ;
                        trajmat(partid).Bx(end+1) = timemat(timeloop).Bx(j) ;
                        trajmat(partid).By(end+1) = timemat(timeloop).By(j) ;
                        trajmat(partid).Bz(end+1) = timemat(timeloop).Bz(j) ;
                        trajmat(partid).rxy(end+1) = timemat(timeloop).rxy(j) ;
                        trajmat(partid).time(end+1) = timemat(timeloop).time ;
                    end

                end

            end

        elseif typetext(1) == 'p'

%    Read in position data

            goodline = 0 ;
            posloop = posloop + 1 ;
            [headings, filepos] = textscan(fid, hstrp, 1) ;
            [dat, filepos] = textscan(fid, nstrp) ;
            datsize = size(dat) ;

%    Save position data to 'posmat' structure array

            for j = 1:datsize(2)
%                posmat(:,j,posloop) = dat{j} ;
                eval(['posmat(' num2str(posloop) ').' char(colheadp{j}) ' = dat{j} ;']) ;
            end

%    Calculate angles and store screen position data in array

            eval(['posmat(' num2str(posloop) ').xp = atan2(posmat(' num2str(posloop) ...
                ').Bx,posmat(' num2str(posloop) ').Bz) ;']) ;
            eval(['posmat(' num2str(posloop) ').yp = atan2(posmat(' num2str(posloop) ...
                ').By,posmat(' num2str(posloop) ').Bz) ;']) ;
            eval(['posmat(' num2str(posloop) ').position = type{2} ;']) ;

        end

    end

%    Close file

    fclose(fid) ;

%    Save data to .mat-file

    if nargin > 1

        save(outfile,'timemat','posmat') ;

    end

%    If required, save trajectory data to separate .mat-file

    if nargin > 2

        save(trajfile,'trajmat') ;

    end

    return