function [outmat, outfile] = gdfwrite(x, xp, y, yp, z, energy, mass, outfile, enew)
%
% function [outmat, outfile] = gdfwrite(x, xp, y, yp, z, energy, mass, outfile, enew)
%
%    GDFWRITE - Write out ASCII file of particle data for input to GPT.
%
%    gdfwrite(x, xp, y, yp, z, energy, mass, outfile) - creates a text file
%    with the name 'outfile', containing columns of data that can be used
%    as an input file for GPT.
%
%    gdfwrite(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.
%
%    Although most particle data gives emittance information in X and X'
%    (position and angle), GPT requires X and Beta-X (position and
%    transverse velocity).  GDFWRITE converts position and angle
%    information into position and velocity information for each particle
%    to allow particle data sets to be initiated within GPT.
%
%    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 - either 1D array of length (N) of longitudinal particle
%            positions 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).
%
%        MASS - either 1D array of length (N) of particle masses OR single
%               value for all particle masses (AMU).
%
%        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 7) array that contains the data saved in
%                 the GDF 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   y   z   Bx  By  Bx  t
%
%    ...corresponding to the X, Y and Z positions and velocities.  The T
%    column, containing time information in seconds, is automatically set
%    to zero, ensuring that all particles are started simultaneously.  Each
%    row in the data file then corresponds to a single particle.  As such,
%    the lengths of all the input data arrays must match (except for the
%    ENERGY variable, which can be a single value for a mono-energetic
%    beam).
%
%    Note that the ASCII file created by GDFWRITE cannot be used directly
%    by GPT, but must be converted into a GDF data file with the GPT
%    command:
%
%            'asci2gdf -v -o outfile.gdf outfile.txt'

% File is the property of Simon Jolly, Imperial College London.  File
% released under the GNU public license.
%
% File history
%
%   11-Aug-2008 S. Jolly
%       Corrected error in output Energy to allow arbitrary Energy value to
%       be specified
%
%======================================================================

    if nargin < 8
        error('Too few input variables: syntax is [outmat, outfile] = gdfwrite(x, xp, y, yp, z, energy, mass, outfile)') ;
    elseif nargin > 9
        error('Too many input variables: syntax is [outmat, outfile] = gdfwrite(x, xp, y, yp, z, energy, mass, outfile, enew)') ;
    elseif nargout > 2
        error('Incorrect number of output variables: syntax is [outmat, outfile] = gdfwrite(x, xp, y, yp, z, energy, mass, 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(x) ~= size(z) ) & ( length(z) ~= 1 )
        error('Z variable must either be single number or 1D array of equal size to other arrays') ;
    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(mass) ~= size(x) ) & ( length(mass) ~= 1 ) & ( ~isempty(mass) )
        error('Mass 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,7) ;

    if ( length(z) == 1 )
        z = ones(numpart,1).*z ;
    end

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

%    Check for beam parameters; default is H-minus at 35keV

    c  = 299792458 ;            % Speed of light in m/s
    ma = 1.66053886e-27 ;       % Atomic mass in kg
    mzero = 1.0083736 * ma ;    % Mass of H-minus
    qe = -1.60217653e-19 ;      % Electron charge in C
    q0 = qe ;                   % H-minus is singularly charged
    Ezero = 35.0e3 ;            % Beam energy is 35 keV after ion source [eV]

%    Mass and energy of particles

    if isempty( mass )
        partmass = mzero ;
    else
%        partmass = mass.*ma ;
        partmass = mass ;
    end
    if isempty( energy )
        energy = Ezero ;
    end
    if nargin < 9 || 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

%    Calculate velocities and energy spreads

    Gam  = 1 - ((qe.*energy)./(partmass.*(c.^2))) ;   % Corresponding Lorentz factor Gamma []
    BG   = sqrt( (Gam.^2) - 1 ) ;               % Corresponding beta-gamma factor []
    Beta = sqrt( 1 - (Gam.^-2) ) ;              % Corresponding Normalized velocity []
    dE = (max(energy) - min(energy))./mean(energy) ;    % Energy spread [fraction of E0]
    dG = ( Gam-1 ).*dE ;                                % Gamma spread [fraction of Gam]

%    Convert angle data into tranverse velocities

    Bx = Beta.*sin(xp.*1e-3).*sqrt( ( 1 - (sin(yp.*1e-3).^2) )./ ...
        ( 1 - (sin(xp.*1e-3).^2).*(sin(yp.*1e-3).^2) ) ) ;
    By = Beta.*sin(yp.*1e-3).*sqrt( ( 1 - (sin(xp.*1e-3).^2) )./ ...
        ( 1 - (sin(xp.*1e-3).^2).*(sin(yp.*1e-3).^2) ) ) ;
    Bz = sqrt(Beta.^2 - Bx.^2 - By.^2) ;

%     Fill in data array with converted values

    datamat(:,1) = x.*1e-3 ;
    datamat(:,2) = y.*1e-3 ;
    datamat(:,3) = z ;
    datamat(:,4) = Bx ;
    datamat(:,5) = By ;
    datamat(:,6) = Bz ;
%    datamat(:,4) = -Bx ;
%    datamat(:,5) = -By ;
%    datamat(:,6) = -Bz ;

%    Write data to text file

    fid = fopen(outfile, 'w') ;
    fprintf(fid,'x\ty\tz\tBx\tBy\tBz\tt\r\n') ;
    fclose(fid) ;
    dlmwrite(outfile, datamat, '-append', 'delimiter','\t', 'newline','pc', 'precision',10) ;

    outmat = datamat ;

    return

