#
# bunch train statistics
#
function bunch_train_stat(beam,statfile,nbunches,nmacro,nslice)
format long;

fid=fopen(statfile,"wt");
if (fid!=-1)
  fprintf(fid,"\# bunch train data\n");
  fprintf(fid,"\# %18s %20s %20s %20s %20s %20s %20s %20s %20s %20s %20s %20s %20s %20s %20s %20s %20s\n",
              "semit nx [nm rad]","cemit nx [nm rad]","semit ny [nm rad]","cemit ny [nm rad]",
              "mean x [um]","rms x [um]","mean xp [urad]","rms xp [urad]",
              "mean y [um]","rms y [um]","mean yp [urad]","rms yp [urad]",
              "mean s [um]","rms s [um]","mean e [GeV]","rms e [GeV]","dee [1]");
end;

% full bunch train
bf=beam;
sliceemit_x=mean(bf(:,3),bf(:,2))/0.0005109989*sqrt(mean(bf(:,8),bf(:,2))*mean(bf(:,10),bf(:,2)) - mean(bf(:,9),bf(:,2))^2)/1000
sliceemit_y=mean(bf(:,3),bf(:,2))/0.0005109989*sqrt(mean(bf(:,11),bf(:,2))*mean(bf(:,13),bf(:,2)) - mean(bf(:,12),bf(:,2))^2)/1000
coreemit_x =mean(bf(:,3),bf(:,2))/0.0005109989*sqrt(mean(bf(:,4).^2,bf(:,2))*mean(bf(:,5).^2,bf(:,2)) - mean(bf(:,4).*bf(:,5),bf(:,2))^2)/1000
coreemit_y =mean(bf(:,3),bf(:,2))/0.0005109989*sqrt(mean(bf(:,6).^2,bf(:,2))*mean(bf(:,7).^2,bf(:,2)) - mean(bf(:,6).*bf(:,7),bf(:,2))^2)/1000

mean_x =mean(bf(:,4),bf(:,2));
mean_xp=mean(bf(:,5),bf(:,2));
mean_y =mean(bf(:,6),bf(:,2));
mean_yp=mean(bf(:,7),bf(:,2));

rms_x =rms(bf(:,4),bf(:,2));
rms_xp=rms(bf(:,5),bf(:,2));
rms_y =rms(bf(:,6),bf(:,2));
rms_yp=rms(bf(:,7),bf(:,2));

mean_s=mean(bf(:,1),bf(:,2));
mean_e=mean(bf(:,3),bf(:,2));
rms_s =rms(bf(:,1),bf(:,2));
rms_e =rms(bf(:,3),bf(:,2));

dee=rms_e/mean_e;

if (fid!=-1)
    fprintf(fid,"\# full bunch train\n");
    fprintf(fid,"%20.12e %20.12e %20.12e %20.12e %20.12e %20.12e %20.12e %20.12e %20.12e %20.12e %20.12e %20.12e %20.12e %20.12e %20.12e %20.12e %20.12e\n",
	        sliceemit_x,coreemit_x,sliceemit_y,coreemit_y,
	        mean_x,rms_x,mean_xp,rms_xp,
	        mean_y,rms_y,mean_yp,rms_yp,
	        mean_s,rms_s,mean_e,rms_e,dee);
    fprintf(fid,"\# single bunches\n");
end;

% single bunches
for i=1:nbunches
  is=nmacro*nslice*(i-1)+1;
  ie=nmacro*nslice*i;

  bf=beam(is:ie,:);

  printf("--------------------------------------------------\n");
  printf("Bunch %d of %d:\n",i,$nbunches);
  sliceemit_x=mean(bf(:,3),bf(:,2))/0.0005109989*sqrt(mean(bf(:,8),bf(:,2))*mean(bf(:,10),bf(:,2)) - mean(bf(:,9),bf(:,2))^2)/1000
  sliceemit_y=mean(bf(:,3),bf(:,2))/0.0005109989*sqrt(mean(bf(:,11),bf(:,2))*mean(bf(:,13),bf(:,2)) - mean(bf(:,12),bf(:,2))^2)/1000
  coreemit_x =mean(bf(:,3),bf(:,2))/0.0005109989*sqrt(mean(bf(:,4).^2,bf(:,2))*mean(bf(:,5).^2,bf(:,2)) - mean(bf(:,4).*bf(:,5),bf(:,2))^2)/1000
  coreemit_y =mean(bf(:,3),bf(:,2))/0.0005109989*sqrt(mean(bf(:,6).^2,bf(:,2))*mean(bf(:,7).^2,bf(:,2)) - mean(bf(:,6).*bf(:,7),bf(:,2))^2)/1000

  mean_x =mean(bf(:,4),bf(:,2));
  mean_xp=mean(bf(:,5),bf(:,2));
  mean_y =mean(bf(:,6),bf(:,2));
  mean_yp=mean(bf(:,7),bf(:,2));

  rms_x =rms(bf(:,4),bf(:,2));
  rms_xp=rms(bf(:,5),bf(:,2));
  rms_y =rms(bf(:,6),bf(:,2));
  rms_yp=rms(bf(:,7),bf(:,2));

  mean_s=mean(bf(:,1),bf(:,2));
  mean_e=mean(bf(:,3),bf(:,2));
  rms_s =rms(bf(:,1),bf(:,2));
  rms_e =rms(bf(:,3),bf(:,2));

  dee=rms_e/mean_e;

  if (fid!=-1)
  fprintf(fid,"%20.12e %20.12e %20.12e %20.12e %20.12e %20.12e %20.12e %20.12e %20.12e %20.12e %20.12e %20.12e %20.12e %20.12e %20.12e %20.12e %20.12e\n",
	      sliceemit_x,coreemit_x,sliceemit_y,coreemit_y,
	      mean_x,rms_x,mean_xp,rms_xp,
	      mean_y,rms_y,mean_yp,rms_yp,
	      mean_s,rms_s,mean_e,rms_e,dee);
  end;
end;

fclose(fid);

endfunction;
