%
% particle beam setup
%
% values are expected in Placet units and the output is done in Placet units.
% and intermediate conversion to standard units is done via pu_xyz, pu_xpyp, pu_e and pu_emit

function create_particle_distribution(outfilename,bx,ax,enx,by,ay,eny,rmsz,meanz,espread,echirp,energy,nsigma,npart,pu_xyz,pu_xpyp,pu_e,pu_emit)

% first we convert to some useful units ;-)
bx=bx;
ax=ax;
enx=enx/pu_emit;
by=by;
ay=ay;
eny=eny/pu_emit;
rmsz=rmsz/pu_xyz;
meanz=meanz/pu_xyz;
espread=espread;
echirp=echirp*pu_xyz;
energy=energy/pu_e;

ex=enx/energy*510998.9;
ey=eny/energy*510998.9;

sx=sqrt(ex*bx);
sxp=sqrt(ex/bx);
sy=sqrt(ey*by);
syp=sqrt(ey/by);
seunc=abs(espread)*energy;

xxpc=-ax/bx;
yypc=-ay/by;

rv=cutrandnv(npart,nsigma);

x =sx*rv(1,:);
x =x-mean(x,1);
x =x*sx/rms(x,1);

xp=sxp*rv(2,:);
xp=xp-mean(xp,1)+x*1e3;

sumxp=0;
for i=1:length(xp)
  if (x(i)>=0)
    sumxp=sumxp+xp(i);
  end;
end;
sign_ax=-sumxp/abs(sumxp);

ext=sqrt(rms(x,1)^2*rms(xp,1)^2-mean(x.*xp,1)^2);
dxpdx=sqrt(rms(xp,1)^2/rms(x,1)^2-ext^2/rms(x,1)^4);
xp=xp+x*sign_ax*dxpdx;
xp=xp*sxp/rms(xp,1);


y =sy*rv(3,:);
y =y-mean(y,1);
y =y*sy/rms(y,1);

yp=syp*rv(4,:);
yp=yp-mean(yp,1)+y*1e3;

sumyp=0;
for i=1:length(yp)
  if (y(i)>=0)
    sumyp=sumyp+yp(i);
  end;
end;
sign_ay=-sumyp/abs(sumyp);

eyt=sqrt(rms(y,1)^2*rms(yp,1)^2-mean(y.*yp,1)^2);
dypdy=sqrt(rms(yp,1)^2/rms(y,1)^2-eyt^2/rms(y,1)^4);
yp=yp+y*sign_ay*dypdy;
yp=yp*syp/rms(yp,1);


z =rmsz*rv(5,:);
z =z-mean(z,1);
z =z*rmsz/rms(z,1);

e =seunc*rv(6,:);
e =e-mean(e,1);
e =e*seunc/rms(e,1);

xp=xp+x*xxpc;
yp=yp+y*yypc;
e =e+z*echirp*energy+energy;

z=z+meanz;

% and now we convert back to Placet units
x=x*pu_xyz;
xp=xp*pu_xpyp;
y=y*pu_xyz;
yp=yp*pu_xpyp;
z=z*pu_xyz;
e=e*pu_e;

fid=fopen(outfilename,"wt");
outdata=[e',x',y',z',xp',yp'];
fprintf(fid,"%20.12e %20.12e %20.12e %20.12e %20.12e %20.12e\n",outdata');
fclose(fid);

endfunction;
