function [outmat, outfile] = lintrawrite(x, xp, y, yp, z, energy, outfile, enew)
%
% function [outmat, outfile] = lintrawrite(x, xp, y, yp, z, energy, outfile, enew)
%
%    LINTRAWRITE - Write out ASCII file of particle data for input to LINTRA.
%
%    lintrawrite(x, xp, y, yp, z, energy, outfile) - creates a text file
%    with the name 'outfile', containing columns of data that can be used
%    as an input file for LINTRA.
%
%    lintrawrite(x, xp, y, yp, z, energy, mass, outfile, enew) - scales up
%    the emittances from ENERGY to ENEW eg. a reading taken at 35keV can be
%    output as a 65keV beam by setting ENERGY = 35e3 and ENEW = 65e3.
%
%    To create an input file for N particles, the required input variables
%    are:
%
%        X - 1D array of length (N) of particle X position (mm).
%
%        XP - 1D array of length (N) of particle X angle (mrad).
%
%        Y - 1D array of length (N) of particle Y position (mm).
%
%        YP - 1D array of length (N) of particle Y angle (mrad).
%
%        Z - 1D array of length (N) of longitudinal particle position OR
%            single value for all particle positions (m).
%
%        ENERGY - either 1D array of length (N) of particle energies OR
%                 single value for all particle energies (eV).
%
%        OUTFILE - string of the full name of the file to be saved
%                  (including extension).
%
%    Note that the transverse positions are specified in millimetres, but
%    the longitudinal position is specified in metres.  All input variable
%    must be specified.
%
%    Optional input parameters:
%
%        ENEW - single value for all new particle energies (eV).
%
%    The 2 output variables are:
%
%        OUTMAT - returns the (N x 6) array that contains the data saved in
%                 the ASCII file.
%
%        OUTFILE - echoes the input variable of the same name.
%
%    The output file is created with columns of data in the following
%    order:
%
%        x   x'   y   y'   z   E
%
%    This is preceeded by a single number indicating the number of
%    particles in the dataset.
%
%    Each row in the data file corresponds to a single particle.  As such,
%    the lengths of all the input data arrays must match (except for the
%    ENERGY and Z variables, which can be single values for 2-D slices and
%    mono-energetic beams).
%

    if nargin < 7
        error('Too few input variables: syntax is [outmat, outfile] = lintrawrite(x, xp, y, yp, z, energy, outfile)') ;
    elseif nargin > 8
        error('Too few input variables: syntax is [outmat, outfile] = lintrawrite(x, xp, y, yp, z, energy, outfile, enew)') ;
    elseif nargout > 2
        error('Incorrect number of output variables: syntax is [outmat, outfile] = lintrawrite(x, xp, y, yp, z, energy, outfile)') ;
    elseif ( length(x) ~= length(xp) ) || ( length(y) ~= length(yp) ) || ...
            ( length(x) ~= length(y) )
        error('EmitCalc:incorrectSyntax', ...
            ['Input variables are different lengths: must contain same number '...
            '\n of data points ie. position and angle values for each particle']) ;
    elseif ( size(energy) ~= size(x) ) & ( length(energy) ~= 1 )
        error('Energy variable must either be single number or 1D array of equal size to other arrays') ;
    elseif ( size(z) ~= size(x) ) & ( length(z) ~= 1 )
        error('Z variable must either be single number or 1D array of equal size to other arrays') ;
    elseif ~ischar(outfile)
        error('Output filename must be a string array') ;
    end

%    Initialise variables for datafile

    numpart = length(x) ;
    datamat = zeros(numpart,6) ;

    if ( length(z) == 1 )
        zarr = z.*( ones(length(x),1) ) ;
    else
        zarr = z ;
    end

    disp(' ')
    disp(['Number of particles in dataset: ' num2str(numpart) '. '])
    disp(' ')

    if nargin < 8 || isempty(enew)
        enew = energy ;
    end

%    Rescale emittances for different output energy

    if enew ~= energy
        eold = mean( energy ) ;
        ediff = enew - eold ;
        xp = xp.*( sqrt( eold./enew ) ) ;
        yp = yp.*( sqrt( eold./enew ) ) ;
        energy = energy + ediff ;
    end

    if ( length(energy) == 1 )
        dEarr = zeros(numpart,1) ;
    else
        dEarr = energy - mean(energy) ;
    end

%     Fill in data array with converted values

    datamat(:,1) = x ;
    datamat(:,2) = xp ;
%    datamat(:,2) = -xp ;
    datamat(:,3) = y ;
    datamat(:,4) = yp ;
%    datamat(:,4) = -yp ;
%    datamat(:,5) = zarr ;
    datamat(:,6) = dEarr ;

%    Write data to text file

    fid = fopen(outfile, 'w') ;
    fprintf(fid,'% -7.0f\r\n',numpart) ;
    fprintf(fid,'% 15.7f % 15.7f % 15.7f % 15.7f % 15.7g % 15.7g\r\n',datamat') ;
    fclose(fid) ;
%    dlmwrite(outfile, datamat, '-append', 'delimiter','\t', 'newline','pc', 'precision',10) ;

    outmat = datamat ;

    return

