function [firstdat, lastdat, colours, rangemat, lossbins, binmids] = gptlosses(trajmat, ranges)
%
% function [firstdat, lastdat, colours, rangemat, lossbins, binmids] = gptlosses(trajmat, ranges)
%
%    GPTLOSSES - Analyse GPT trajectory data and mark lost particles.
%
%    GPTLOSSES is used to analyse GPT trajectory data and colour code each
%    particle according to the z-position at which the particle was lost.
%    The input data, TRAJMAT, needs to be a Matlab structure array
%    containing trajectory data created by the function GDFTRAJ2MAT.
%
%    [firstdat, lastdat] = gptlosses(trajmat)
%
%    [firstdat, lastdat, colours] = gptlosses(trajmat)
%
%    [firstdat, lastdat, colours] = gptlosses(trajmat, ranges)
%
%    [firstdat, lastdat, colours, rangemat] = gptlosses(...)
%
%    [firstdat, lastdat, colours, rangemat, lossbins] = gptlosses(...)
%
%    [firstdat, lastdat, colours, rangemat, lossbins, binmids] = gptlosses(...)
%
%    [firstdat, lastdat] = gptlosses(trajmat) - reads in trajectory data in
%    TRAJMAT and outputs 2 arrays of particle data:
%
%        FIRSTDAT - initial particle distribution ie. from TRAJMAT.NN(1)
%                   where NN is any of the data structures in TRAJMAT (eg.
%                   x, time, ID etc.).
%
%        LASTDAT - final particle distribution, containing the last data
%                  point at which each particle was lost ie. from
%                  TRAJMAT.NN(end).  This is not a real distribution.
%
%    The FIRSTDAT and LASTDAT arrays are both [ N x 10 ] arrays containing
%    10 columns of data for N particles.  The ordering of the columns is as
%    follows:
%
%        (N,1) - 'x' - particle x position
%        (N,2) - 'xp' - particle xp data
%        (N,3) - 'y' - particle y position
%        (N,4) - 'yp' - particle yp data
%        (N,5) - 'z' - particle z position
%        (N,6) - 'Bx' - particle Bx data
%        (N,7) - 'By' - particle By data
%        (N,8) - 'Bz' - particle Bz data
%        (N,9) - 'rxy' - particle rxy distance
%        (N,10) - 'time' - particle timeslice
%
%    The particle data is ordered by particle ID number ie. the order they
%    are stored in TRAJMAT.
%
%    [firstdat, lastdat, colours] = gptlosses(trajmat) - analyses particle
%    data in array  LASTDAT and creates an [ N x 3 ] array of colour
%    information of the form 'xk-'.  'Surviving' particles are determined
%    by using ZMAX = max(LASTDAT(:,5) - particles with Z==ZMAX are marked
%    as surviving and are coloured black ('xk-'), while those with Z<ZMAX
%    do not and are coloured red ('xr-').
%
%    [...] = gptlosses(lastdat, ranges) - additionally specifies the range
%    data in which to 'bin' the particle losses.
%
%    For example, to colour code particle losses in units of 10cm, from
%    0-50cm, set RANGES = [0:0.1:0.5].  As such, particles lost before 0cm
%    are coloured red, those between 0-10cm green, between 10-20cm blue
%    etc. and those that survive beyond 50cm are coloured black.  The
%    colours cycle through [R G B C M Y] for consecutive bins: when colours
%    start to repeat, different marker shapes are used.  Surviving
%    particles are always coloured black.
%
%    When specifying RANGES, ZMAX = RANGES(end): hence it is possible that
%    no particle survives (and none are coloured black).  To change the
%    choice of ZMAX, simply set RANGES = ZMAX for the ZMAX of your choice.
%
%    [firstdat, lastdat, colours, rangemat] = gptlosses(...) - return the
%    ranges used to assign colours.  RANGEMAT = [-Inf RANGES Inf].
%
%    [firstdat, lastdat, colours, rangemat, lossbins] = gptlosses(...) -
%    output the particles lost in each bin specified by RANGES to the
%    LOSSBINS array.
%
%    [firstdat, lastdat, colours, rangemat, lossbins, binmids] = gptlosses(...) 
%    - output the midpoint of each bin used to bin the particle data to the
%    BINMIDS array.  The first and last bins, which are nominally in the
%    range [-Inf:0] and [ZMAX:Inf] are given artificially given a midpoint
%    such that they are the same width as the adjacent bins.
%

    if nargin < 1
        error('Too few input variables: syntax is [firstdat, lastdat] = gptlosses(trajmat)') ;
    elseif nargin > 2
        error('Too many input variables: syntax is [firstdat, lastdat] = gptlosses(trajmat, ranges)') ;
    elseif nargout > 6
        error('Too many output variables: syntax is [firstdat, lastdat, colours, rangemat, lossbins, binmids] = gptlosses(trajmat, ranges)') ;
    end

    numpart = length(trajmat) ;
%    looppart = round( numpart./20 ) ;
    firstdat = zeros(numpart,10) ;
    lastdat = firstdat ;
    fields = {'x';'xp';'y';'yp';'z';'Bx';'By';'Bz';'rxy';'time'} ;

    for i = 1:numpart

        eval(['check = trajmat(i).' fields{1} ' ;']) ;

        if ~isempty(check)
        
            for j = 1:length(fields)

                eval(['firstdat(i,j) = trajmat(i).' fields{j} '(1) ;']) ;
                eval(['lastdat(i,j) = trajmat(i).' fields{j} '(end) ;']) ;

            end

        end

        if ( rem(i,1000) == 0 )
            disp(['Reading particle ' num2str(i) ' ... ']) ;
        end

    end

    if nargin < 2 || isempty(ranges)
        ranges = max(lastdat(:,5)) ;
    end

    if nargout >= 3

        [colours, rangemat, lossbins, binmids] = gptloscol(lastdat, ranges) ;

    end

    return