function [colours, rangemat, lossbins, binmids] = gptloscol(lastdat, ranges)
%
% function [colours, rangemat, lossbins, binmids] = gptloscol(lastdat, ranges)
%
%    GPTLOSCOL - Colour code particle loss data.
%
%    GPTLOSCOL provides colour-coding for particle trajectory data, based
%    on the z-position at which the particle was lost.  For N particles,
%    the input array LASTDAT needs to be an [ N x 10 ] array as produced by
%    GPTLOSSES - see HELP GPTLOSSES for more details.
%
%    [colours] = gptloscol(lastdat)
%
%    [colours] = gptloscol(lastdat, ranges)
%
%    [colours, rangemat] = gptloscol(...)
%
%    [colours, rangemat, lossbins] = gptloscol(...)
%
%    [colours, rangemat, lossbins, binmids] = gptloscol(...)
%
%    [colours] = gptloscol(lastdat) - 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-').
%
%    [...] = gptloscol(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.
%
%    [colours, rangemat] = gptloscol(...) - return the ranges used to
%    assign colours.  RANGEMAT = [-Inf RANGES Inf].
%
%    [colours, rangemat, lossbins] = gptloscol(...) - output the particles
%    lost in each bin specified by RANGES to the LOSSBINS array.
%
%    [colours, rangemat, lossbins, binmids] = gptloscol(...) - 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 [colours, rangemat, lossbins, binmids] = gptloscol(lastdat, ranges)') ;
    elseif nargin > 2
        error('Too many input variables: syntax is [colours, rangemat, lossbins, binmids] = gptloscol(lastdat, ranges)') ;
    elseif nargout > 4
        error('Too many output variables: syntax is [colours, rangemat, lossbins, binmids] = gptloscol(lastdat, ranges)') ;
    end

    endpart = max(lastdat(:,5)) ;

    if nargin < 2
        ranges = endpart ;
    end

    numpart = size(lastdat,1) ;
    colours = num2str(zeros(numpart,2)) ;
    colours(:,4:end) = [] ;
    colloop = ('rgbcmy') ;
    asloop = ('x+o*sd^v><ph') ;
    colnum = 0 ;
    zeds = lastdat(:,5) ;

    rangemat = [-Inf ranges Inf] ;
    lossbins = zeros(length(rangemat)-1,1) ;

    for i = 1:(length(rangemat) - 1)

        low = rangemat(i) ; up = rangemat(i+1) ;
        [pinds] = find( ( zeds >= low ) & ( zeds < up ) ) ;

        if ~isempty(pinds)

            colnum = colnum + 1 ;
            lossbins(i) = length(pinds) ;

            if up == Inf
                colchar = 'xk-' ;
            else
                colchar = [asloop(fix((colnum-1)./6)+1) colloop(rem(colnum-1,6)+1) '-'] ;
            end

            for j = 1:length(pinds)
                colours(pinds(j),:) = colchar ;
            end

        end

    end

    if length(ranges) == 1
        if ranges == 0 || ( ranges > endpart )
            binmids = ([(endpart./2), (endpart.*1.5)]) ;
        else
            rangediff = (endpart - ranges)./2 ;
            binmids = ([rangediff, ranges + rangediff]) ;
        end
    else
        rangelow = ranges(1:end-1) ;
        rangehigh = ranges(2:end) ;
        rangediff = ( rangehigh - rangelow )./2 ;
        binmids = ([rangelow(1) - rangediff(1), [rangelow + rangediff], rangehigh(end) + rangediff(end)]) ;
    end

    return


