#
# 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]","std x [um]","mean xp [urad]","std xp [urad]",
              "mean y [um]","std y [um]","mean yp [urad]","std yp [urad]",
              "mean s [um]","std s [um]","mean e [GeV]","std 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));

std_x =std(bf(:,4),bf(:,2));
std_xp=std(bf(:,5),bf(:,2));
std_y =std(bf(:,6),bf(:,2));
std_yp=std(bf(:,7),bf(:,2));

mean_s=mean(bf(:,1),bf(:,2));
mean_e=mean(bf(:,3),bf(:,2));
std_s =std(bf(:,1),bf(:,2));
std_e =std(bf(:,3),bf(:,2));

dee=std_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,std_x,mean_xp,std_xp,
	        mean_y,std_y,mean_yp,std_yp,
	        mean_s,std_s,mean_e,std_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));

  std_x =std(bf(:,4),bf(:,2));
  std_xp=std(bf(:,5),bf(:,2));
  std_y =std(bf(:,6),bf(:,2));
  std_yp=std(bf(:,7),bf(:,2));

  mean_s=mean(bf(:,1),bf(:,2));
  mean_e=mean(bf(:,3),bf(:,2));
  std_s =std(bf(:,1),bf(:,2));
  std_e =std(bf(:,3),bf(:,2));

  dee=std_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,std_x,mean_xp,std_xp,
	      mean_y,std_y,mean_yp,std_yp,
	      mean_s,std_s,mean_e,std_e,dee);
  end;
end;

fclose(fid);

endfunction;
