## -*- texinfo -*-
## @deftypefn {Function File} {} placet_make_beam (@var{e0}, @var{e_spread}, @var{sigma_Z}, @var{twiss}, @var{N})
## Return a 6-column beam matrix, in guinea-pig format, with the particles.
##
## Input:
## @verbatim
##   - e0 : initial energy [GeV]
##   - e_spread : energy_spread [%]
##   - sigma_Z : bunch length [um]
##   - twiss : twiss parameters structure, twiss.beta_x [m], twiss.alpha_x, twiss.emitt_x [1e-7m], ...
##   - N : number of particles
## @end verbatim
## @end deftypefn

function B = placet_make_beam(e0, e_spread, sigma_Z, twiss, N)
  if nargin==5 && ischaruct(twiss)
    gamma = e0 / 0.000510998903076601;
    ex = twiss.emitt_x * 1e-7 / gamma;
    ey = twiss.emitt_y * 1e-7 / gamma;
    sx = sqrt(ex * twiss.beta_x) * 1e6;
    sxp = sqrt(ex / twiss.beta_x) * 1e6;
    sy = sqrt(ey * twiss.beta_y) * 1e6;
    syp = sqrt(ey / twiss.beta_y) * 1e6;
    se = 0.01 * e_spread;
    B = zeros(N,6);
    if e_spread < 0
      B(:,1) = (1. + se * (rand(N,1) - 0.5)) * e0;
    else
      B(:,1) = (1. + randn(N,1) * se) * e0;
    endif
    B(:,2) = randn(N,1) * sx;
    B(:,3) = randn(N,1) * sy;
    B(:,4) = randn(N,1) * sigma_Z; 
    B(:,5) = (randn(N,1) - twiss.alpha_x * B(:,2) / sx) * sxp; 
    B(:,6) = (randn(N,1) - twiss.alpha_y * B(:,3) / sy) * syp;
  else
    help placet_make_beam
  endif
endfunction
