Index: /CLIO/Analysis/Data19July/Analyse2.m
===================================================================
--- /CLIO/Analysis/Data19July/Analyse2.m	(revision 723)
+++ /CLIO/Analysis/Data19July/Analyse2.m	(revision 723)
@@ -0,0 +1,206 @@
+clear all; close all; clc;
+date=20170719;
+load('data1907.mat');
+%cut by electron signal
+%cut on signal on 12's detector (zero signal); and 11's detector (DAQ problem) 
+cut=(key_value(:,8)>3e+3);%&(abs(data(:,11,2))<100);%&(abs(data(:,12,2))<50);%&(sum(data(:,:,2)~=0,2));
+ %plot(rem(key_value(cut,4),date),data(cut,11,2),'or')
+% plot(rem(key_value(:,6),date),key_value(:,4),'*r')
+%plot(rem(key_value(cut,6),date),key_value(cut,8),'*r')
+
+position=key_value(:,4)*7.94e-6;
+hour=rem(key_value(:,6),date);
+
+ diap(1,:)=(hour>84900)&(hour<85750);
+ diap(2,:)=(hour>85750)&(hour<90400);
+ diap(3,:)=(hour>90400)&(hour<91100);
+ diap(4,:)=(hour>91100)&(hour<91800);
+ diap(5,:)=(hour>91800)&(hour<92500);
+ diap(6,:)=(hour>92500)&(hour<93300);
+ diap(7,:)=(hour>93300)&(hour<94000);
+ diap(8,:)=(hour>94000)&(hour<94700);
+
+
+fromd=-18;
+bkgsh=+4;
+for detector=1:11
+for dd=1:8
+cut2=cut&diap(dd,:)';
+x=position(cut2)-48.58-3.46;
+tmpE=key_value(cut2,8);
+Es(dd)=mean(tmpE(abs(tmpE-mean(tmpE))<std(tmpE)));
+Ess(dd)=std(tmpE(abs(tmpE-mean(tmpE))<std(tmpE)));
+
+y=data(cut2,detector,2)./Es(dd);
+
+fb(dd,detector) = polyfit(x(x<(fromd+bkgsh)),y(x<(fromd+bkgsh)),0);
+
+%[fs, gofs] = fit(x(x>fromd),y(x>fromd),'exp1');
+[fs, gofs] = fit(x((x>fromd)&y>fb(dd,detector)),y((x>fromd)&y>fb(dd,detector)),'exp1','Lower',[0 0.015]);
+chi(dd,detector) = gofs.adjrsquare;
+lamb(dd,detector)=fs.b;
+amp(dd,detector)=fs.a;
+% figure
+% plot(fs,x,y)
+tmp1=confint(fs);
+tmp2=(tmp1(1:2:3)-tmp1(2:2:4))/4;
+Ierrp(dd,detector,:)=tmp2(:);
+
+tmp1=confint(fs);
+tmp2=(tmp1(1:2:3)-tmp1(2:2:4))/4;
+Ierrp(dd,detector,:)=tmp2(:);
+
+
+end 
+end
+
+dd=5;
+cut2=cut&diap(dd,:)';
+x=position(cut2)-48.58-3.46;
+xfits=fromd:1:(-7);
+xfitb=(-22):1:fromd;
+
+figure 
+hold on
+detector=7;
+y=data(cut2,detector,2)./mean(key_value(cut2,8))*1e+3;
+semilogy(x,y,'^r','linewidth',2)
+semilogy(xfits,amp(dd,detector)*exp(xfits*lamb(dd,detector))*1e+3,'r','linewidth',2)
+semilogy(xfitb,repmat(fb(dd,detector),size(xfitb))*1e+3,'Color',[0.8 0.8 0.8],'linewidth',2)
+set(gca, 'YScale', 'log')
+detector=5;
+y=data(cut2,detector,2)./mean(key_value(cut2,8))*1e+3;
+semilogy(x,y,'>k','linewidth',2)
+semilogy(xfits,amp(dd,detector)*exp(xfits*lamb(dd,detector))*1e+3,'k','linewidth',2)
+semilogy(xfitb,repmat(fb(dd,detector),size(xfitb))*1e+3,'Color',[0.8 0.8 0.8],'linewidth',2)
+
+detector=9;
+y=data(cut2,detector,2)./mean(key_value(cut2,8))*1e+3;
+semilogy(x,y,'<b','linewidth',2)
+semilogy(xfits,amp(dd,detector)*exp(xfits*lamb(dd,detector))*1e+3,'b','linewidth',2)
+semilogy(xfitb,repmat(fb(dd,detector),size(xfitb))*1e+3,'Color',[0.8 0.8 0.8],'linewidth',2)
+
+set(gca,'fontsize',16)
+legend('76^o','Signal','Bkgnd','90^o','Signal','Bkgnd', '104^o','Signal','Bkgnd')
+xlabel('Position, [mm]')
+ylabel('Amplitude, [a.u.]')
+grid on 
+title('\phi_B=4.55')
+ xlim([-22 0])
+ ylim([10 120])
+
+
+Bphase=sort([4.3 4.4 4.45 4.5 4.55 4.6 4.65 4.7]);
+
+figure
+errorbar(Bphase,Es/1000,Ess/1000,'-*r','linewidth',2)
+set(gca,'fontsize',16)
+xlabel('\phi_B, [a.u.]')
+ylabel('Electron intensity, [a.u.]')
+grid on 
+
+
+
+
+figure
+chi(chi<0)=0;
+[hc hc]=contourf(Bphase(:), 48:7:118, chi',40)
+set(hc,'LineStyle','none');
+colorbar
+set(gca,'fontsize',16)
+ ylabel('Angle, [deg.]')
+ xlabel('\phi_B, [a.u.]')
+title('R^2')
+
+figure
+[hc hc]=contourf(Bphase(:), 48:7:118, 2./lamb',40)
+set(hc,'LineStyle','none');
+colorbar
+set(gca,'fontsize',16)
+ ylabel('Angle, [deg.]')
+ xlabel('\phi_B, [a.u.]')
+title('Evanescent wavelength, [mm]')
+
+%3d diap is bad
+
+figure
+hold on 
+errorbar(Bphase',2./lamb(:,4),2./lamb(:,4).^2.*Ierrp(:,4,2),'-^b','linewidth',2)
+errorbar(Bphase',2./lamb(:,7),2./lamb(:,7).^2.*Ierrp(:,7,2),'-*r','linewidth',2)
+errorbar(Bphase',2./lamb(:,10),2./lamb(:,10).^2.*Ierrp(:,10,2),'-ok','linewidth',2)
+set(gca, 'YScale', 'log')
+set(gca,'fontsize',16)
+legend('69^o','90^o', '111^o')
+xlabel('\phi_B, [a.u.]')
+ylabel('Evanescent wavelength, [mm]')
+grid on 
+
+
+figure
+hold on
+plot(Bphase',amp(:,4),'-^b','linewidth',2)
+plot(Bphase',amp(:,7),'-*r','linewidth',2)
+plot(Bphase',amp(:,10),'-ok','linewidth',2)
+
+set(gca,'fontsize',16)
+legend('69^o','90^o', '111^o')
+xlabel('\phi_B, [a.u.]')
+ylabel('Amplitude, [a.u.]')
+grid on 
+
+
+
+figure
+[hc hc]=contourf(Bphase(:), 48:7:118, fb',40)
+  set(hc,'LineStyle','none');
+colorbar
+set(gca,'fontsize',16)
+ ylabel('Angle, [deg.]')
+ xlabel('\phi_B, [a.u.]')
+title('Background, [a.u.]')
+%  (hf,'-depsc','file.eps')
+
+figure
+[hc hc]=contourf(Bphase(:), 48:7:118, amp',40)
+  set(hc,'LineStyle','none');
+colorbar
+set(gca,'fontsize',16)
+ ylabel('Angle, [deg.]')
+ xlabel('\phi_B, [a.u.]')
+title('Amplitude, [a.u.]')
+%  (hf,'-depsc','file.eps')
+
+
+
+
+figure
+hold on
+ 
+ 
+gam=44.2/0.511+1;
+bet=sqrt(1-1/gam^2);
+phi=2/180*pi;
+plot(48:7:118,6*(1/bet-cos((48:7:118)/180*pi))*bet*gam/(2*pi)./(sqrt(1+(bet*gam*sin((48:7:118)/180*pi)*sin(phi)).^2)),'m','linewidth',2)
+phi=1/180*pi;
+plot(48:7:118,6*(1/bet-cos((48:7:118)/180*pi))*bet*gam/(2*pi)./(sqrt(1+(bet*gam*sin((48:7:118)/180*pi)*sin(phi)).^2)),'g','linewidth',2)
+phi=0/180*pi;
+plot(48:7:118,6*(1/bet-cos((48:7:118)/180*pi))*bet*gam/(2*pi)./(sqrt(1+(bet*gam*sin((48:7:118)/180*pi)*sin(phi)).^2)),'b','linewidth',2)
+  
+W=1./(Ierrp(:,:,2)).^2;
+chilevel=0;
+W(isnan(W))=0;
+for c=1:11
+    LL(c)=sum(lamb((chi(:,c)>chilevel),c).*W((chi(:,c)>chilevel),c))/sum(W((chi(:,c)>chilevel),c));
+    LE(c)=sqrt(1/sum(W((chi(:,c)>chilevel),c)));
+end
+errorbar(48:7:118,2*(LL.^-1),2./(LL.^2).*LE,'-..k','linewidth',2)
+EWAW=2*(LL.^-1);
+EWAWERR=2./(LL.^2).*LE;
+save EW.mat EWAW EWAWERR
+ grid on
+ set(gca,'fontsize',16)
+ xlabel('Angle, [deg.]')
+ ylabel('Evanescent wavelength, [mm]')
+ legend('\phi=2^o','\phi=1^o','\phi=0^o','data','Location','northwest')
+ 
+ save analyze.mat
Index: /CLIO/Analysis/Data19July/Analyse4.m
===================================================================
--- /CLIO/Analysis/Data19July/Analyse4.m	(revision 723)
+++ /CLIO/Analysis/Data19July/Analyse4.m	(revision 723)
@@ -0,0 +1,78 @@
+close all; clear all; clc;
+load analyze.mat
+
+
+%% Spectrum
+
+
+separation=-10;
+
+figure
+hold on
+
+
+
+figure(9) 
+hold on
+Spectrum=amp(6,:).*exp(lamb(6,:)*separation);
+errorbar(48:7:118,Spectrum,Spectrum.*sqrt((Ierrp(6,:,1)./amp(6,:)).^2+(separation*Ierrp(6,:,2)).^2),'k','linewidth',2)
+ 
+%  plot(48:7:118, SP_spec_1(1:2:(end-1))*5e+16,'k','linewidth',2)
+load('GFW/Spec3.mat')
+ plot(48:7:118, Spectr(1:11)./max(Spectr(1:11))*max(Spectrum),'--or','linewidth',2)
+load('GFW/Spec.mat')
+ plot(48:7:118, Spectr(1:11)./max(Spectr(1:11))*max(Spectrum),'-^r','linewidth',2)
+load('GFW/Spec7.mat')
+ plot(48:7:118, Spectr(1:11)./max(Spectr(1:11))*max(Spectrum),'-.*r','linewidth',2)
+
+set(gca,'fontsize',16)
+
+legend('Data','3ps','5ps','7ps','location','northeast')
+xlabel('Angle, [deg.]')
+ylabel('Amplitude, [a.u.]')
+title('Spectrum')
+grid on 
+% ylim([0 0.02])
+% load('GFW/SEY.mat')
+
+
+figure(10) 
+  [hc hc]=contourf (48:7:118, Bphase, amp.*exp(lamb.*separation))
+    set(hc,'LineStyle','none');
+
+ set(gca,'fontsize',16)
+ xlabel('Angle, [deg.]')
+ ylabel('\phi_B, [a.u.]')
+ title('Spectrum map')
+ 
+ 
+ load('GFW/SEY.mat')
+
+figure(10341) 
+ [hc hc]=contourf(300./(6*(1-cos((48:7:118)/180*pi))), Bphase, amp.*exp(lamb.*(-5))./repmat((Integral(1:11)),[8 1])*1e-17,40)
+ colorbar
+  set(hc,'LineStyle','none');
+ set(gca,'fontsize',16)
+ xlabel('Frequency, [GHz]')
+ ylabel('\phi_B, [a.u.]')
+ title('Form Factor (BGS=5 mm)')
+
+figure
+hold on
+plot(300./(6*(1-cos((48:7:118)/180*pi))),  amp(1,:).*exp(lamb(1,:).*(-5))./(Integral(1:11))*1e-24,'-^b','linewidth',2)
+plot(300./(6*(1-cos((48:7:118)/180*pi))),  amp(3,:).*exp(lamb(3,:).*(-5))./(Integral(1:11))*1e-24,'-^m','linewidth',2)
+plot(300./(6*(1-cos((48:7:118)/180*pi))),  amp(4,:).*exp(lamb(4,:).*(-5))./(Integral(1:11))*1e-24,'-^k','linewidth',2)
+ plot(300./(6*(1-cos((48:7:118)/180*pi))),  amp(5,:).*exp(lamb(5,:).*(-5))./(Integral(1:11))*1e-24,'-^g','linewidth',2)
+plot(300./(6*(1-cos((48:7:118)/180*pi))),  amp(7,:).*exp(lamb(7,:).*(-5))./(Integral(1:11))*1e-24,'-^r','linewidth',2)
+
+ set(gca,'fontsize',16)
+ xlabel('Frequency, [GHz]')
+ legend('\phi_B=5.84','\phi_B=5.64','\phi_B=5.54','\phi_B=5.44','\phi_B=6.04')
+ title('Form Factor (BGS=10 mm)')
+ylabel('Amplitude, [a.u.]')
+grid on
+Freq=300./(6*(1-cos((48:7:118)/180*pi)));
+for c=1:7
+FF(c,:)=amp(c,:).*exp(lamb(c,:).*(-10))./(Integral(1:11))*1e-24;
+end
+save FFX.mat Freq FF
Index: /CLIO/Analysis/Data19July/GFW/Corrections.m
===================================================================
--- /CLIO/Analysis/Data19July/GFW/Corrections.m	(revision 723)
+++ /CLIO/Analysis/Data19July/GFW/Corrections.m	(revision 723)
@@ -0,0 +1,65 @@
+close all; clear all; clc
+load CLIOSpec.mat
+
+phi=0:0.5:6;
+for c=2:(size(bin,2)-2)
+   bin2(:,c-1)=bin_single(:,c);
+%    bin2(:,c-1)=bin(:,c);
+
+end
+ SEY=[flip(bin2(:,2:end)')',bin2]';
+
+% figure
+%  [hc hc]=contourf(thet,[ -flip(phi(2:end)),phi],SEY)
+%   set(hc,'LineStyle','none');
+phiNew=-6:0.1:6;
+ for c=1:size(SEY,2)
+     SEYi(:,c)=pchip([ -flip(phi(2:end)),phi],SEY(:,c),phiNew);
+ end
+ 
+%  figure
+%  [hc hc]=contourf(thet,phiNew,SEYi)
+%   set(hc,'LineStyle','none');
+%   
+ Thet=30:0.1:140;
+
+  for c=1:size(SEYi,1)
+     SEYij(c,:)=pchip(thet,SEYi(c,:), Thet);
+  end
+   figure
+ [hc hc]=contourf(Thet,phiNew,SEYij)
+  set(hc,'LineStyle','none');
+  TH=48:7:125
+  
+  Aed_mm=25;%Aperture entrance diameter,  mm 
+lgr_mm=180;%Grating length, mm
+el_mm=6;
+
+NormKV=[1.3458    1.2208    1.1325    1.0711    1.0305    1.0072    1.0002    1.0072 1.0305    1.0711    1.1325    1.2208];%180 length grating, 25 input mirror, 48:7:125
+lgr_e_mm = abs(Aed_mm./sin(Thet/180*pi));           % variable effective grating length, determined basically by the cone entrance diameter of XXXmm*/
+lgr_e_mm(lgr_e_mm>lgr_mm)=lgr_mm;
+CR=0.1.*lgr_e_mm(:);%J;
+
+
+SEYij=SEYij.*repmat(CR,[1 size(SEYij,1)])';
+
+ for  c=1:12
+  
+load(['map' num2str(TH(c)) '.mat']) 
+
+MB=Map(:,41:end-40)'./max(max(Map));
+
+CSEY=SEYij(:,[find(Thet>=TH(c)-10&Thet<=TH(c)+10)]);
+  figure
+  contourf(MB.*CSEY)
+  Integral(c)=sum(sum(MB.*CSEY))*(0.1/180*pi)^2.*NormKV(c);%*
+% SEYij(:,find(Thet>=38&Thet<=58))
+ end
+ figure
+ hold on
+ plot(TH,Integral)
+Spectr=Integral
+ 
+%   save SPRP15.mat Integral
+%    save SPRP10.mat  Spectr
+
Index: /CLIO/Analysis/Data19July/GFW/GFWnew.m
===================================================================
--- /CLIO/Analysis/Data19July/GFW/GFWnew.m	(revision 723)
+++ /CLIO/Analysis/Data19July/GFW/GFWnew.m	(revision 723)
@@ -0,0 +1,423 @@
+  function  [nu_Hz,SP_spec_1,SP_spec_N,lam,E_total,sum_tr,thet,bin_single, bin]=GFWnew(var)
+
+%  single_elec
+% close all; clear all; clc
+%% Configuration of calculation
+%PrSelect == what profile to use: 1 --from file; 2 - assymmetric gaussian
+%P_file -- name and path to file generated with ASTRA (use TakeData.m script on svn=> CLIO/Scripts to take data from ASTRA generated files)
+%PlotFunc:
+%   1 -- Form factor
+%   2 -- Single electron yield from Theta
+%   3 -- dI_N/dOmega (SP spetrum for whole bunch)
+%   4 -- Bunch profile
+%   5 -- dI(nu)/dnu (SP spetrum for signle elctron bunch)
+% var='PrSelect=1; P_file=''../p20.mat''; PlotFunc=[1 1 1 1 1]; ';
+PrSelect=2;
+PlotFunc=[0 0 0 0 0];
+P_file='tmp.mat';
+ %% Physical constants
+c=299792458;%m/s speed of ligth
+qe = 1.602176487e-19*c*10;%Electron charge, statColumns
+%% Parameters block (by default)
+T=1;%Transmissivity (from gun exit to experimental position)
+m=1;%Order of radiation
+lgr_mm=180;%Grating length, mm
+wid_mm=40;%Grating width, mm
+el_mm=6;%Grating period, mm
+alpha1=30;%Blaze angle, degrees
+gam=100;%relativistic gamma
+Qbunch_nC=1.2;%Bunch charge, nC (if T=1,than chagre at experiment position, if not, than charge at the exit of gun)
+x0_mm=8; %Position of beam centroid above grating, mm
+fwhm_x_mm=3.07; %FWHM_X  of Gaussian beam at the centre of grating, mm
+fwhm_y_mm=3.07; %FWHM_Y  of Gaussian beam at the centre of grating, mm
+pulse_len_ps=10;%Nominal bunch duration (FWHM), ps 
+epsl=1;%Asymmetry factor for assymmetric Gaussian
+BAd_mm=410;% beam - aperture distance ,  mm
+Aed_mm=25;%Aperture entrance diameter,  mm 
+%All angles are in radians
+% phiStop=atan(Aed_mm/2/BAd_mm);% phi go from 0 to phiStop
+ phiStop=6/180*pi;
+% ThetaStart=pi/2-atan(Aed_mm/2/BAd_mm);
+% ThetaStop=pi/2+atan(Aed_mm/2/BAd_mm);
+ThetaStart=alpha1/180*pi;%initial theta angle
+ThetaStop=170/180*pi;%final theta angle
+% ThetaStart=80/180*pi;%initial theta angle
+% ThetaStop=140/180*pi;%final theta angle
+ThetaDiv=10;%steps in theta
+phiDiv=3;%steps in phi
+
+if (~isempty(var))
+strread(var,'%s','delimiter',' ')
+eval(var);
+end %if
+
+%save all parameters
+name=regexprep(var,'.mat','');
+name=regexprep(name,'[../''[]]','');
+name=regexprep(name,'[; ]','_');
+% save([ name 'SP.mat'])
+%% Profile selection
+if PrSelect==1
+    FFTspec=1;
+  if  exist(P_file,'file')==2;
+    load([P_file]);
+  else
+      error('Set correct file name!')
+  end
+else
+    FFTspec=1;%how to calculate form factor: 1 -- FFT; 2 -- only fro assymetric gaussian
+%Gaussian
+sigma_t= pulse_len_ps/(sqrt(2*log(2))*(1+epsl))*1e-3;  %The sigma of the leading edge    ns
+
+nu=(1:1e+10:3e+16).*1e-9;%Ghz
+
+sigma_L= sigma_t*epsl;
+sigma_H= sigma_t;
+time=(1:(length(nu)*2))*(1/(length(nu)*2*(nu(2)-nu(1))));
+xc=mean(time);
+G1=exp(-(time-xc).^2./(2*sigma_L^2));
+G2=exp(-(time-xc).^2./(2*sigma_H^2));
+G1(time>xc)=0;
+G2(time<xc)=0;
+Profile=(G1+G2)./sum(G1+G2);
+end
+
+%% Calculation block for additional parameters
+ n_e_t=round(Qbunch_nC*1e-9/1.6e-19*T);              %Total number of electrons in the bunch
+% n_e_t=1e+11;
+SolAng=pi*((Aed_mm/2)^2)/(BAd_mm^2); %solid angle from beam to aperture
+% Grating
+alpha1=alpha1/180*pi;
+el1=el_mm*(cos(alpha1)^2);
+el2=el_mm-el1;
+hi=abs(el1*tan(alpha1));
+alpha2=-(pi/2-alpha1);
+% Beam 
+beta=sqrt(gam^2 - 1)/gam;
+sigma_x= fwhm_x_mm/(2*sqrt(2*log(2)));
+sigma_y= fwhm_y_mm/(2*sqrt(2*log(2)));
+% N_e above grating
+N_sig=3;
+lol=x0_mm-N_sig*sigma_x;                                       
+if (lol<hi) lol=1.001*hi; end
+upl=x0_mm+N_sig*sigma_x;
+GAU = @(x) exp(-(x-x0_mm).*(x-x0_mm)/(2*sigma_x*sigma_x));
+res1 = integral(GAU,lol,upl)/(sigma_x*sqrt(2*pi));%part of the beam is above the grating
+disp([num2str(res1) ' of the beam above the grating'])
+n_e=n_e_t*res1;                                                 %Number of electrons above the grating
+% save([ name '_additional_SP.mat'])
+%% Estimation grating factor R^2, coherent & incoherent integrals
+%%
+% 
+% $$W_1 |_{x_0=0}= \sum \frac{q^2\omega^3 L l }{4\pi^2|m|c^3}| n  \times n \times \sum G(...)|^2$$
+% 
+% In this block we calculate 
+%%
+%
+% $$| n  \times n \times \sum G(...)|^2 $$
+%
+% where other terms (except L - grating length) are in f1, f2, f4. In general:
+%%
+%
+% $$W_N=W_1 |_{x_0=0}(NS_{inc}+N^2S_{coh}) $$
+%
+% where s_coh and s_inc are:
+%%
+%
+% $$S_{inc}=\int_h^{\infty}dxX(x)$$
+%
+% $$S_{coh}=|\int_h^{\infty}dxX(x)\tilde{Y}(ky)\tilde{T}(\omega)|^2 $$
+%
+% for now (in this block)
+%%
+%
+% $$S^*_{inc}=\int_h^{\infty}dxX(x)*| n  \times n \times \sum G(...)|^2\\ $$
+%
+% $$S^*_{coh}=\int_h^{\infty}dxX(x)*| n  \times n \times \sum G(...)|  $$
+%
+par = struct('Theta',NaN,'phi',NaN,'m',m,'alpha1',alpha1,'facet',NaN,'gamma',gam,'x0',x0_mm,'l',el_mm,'w',wid_mm);
+for i=1:ThetaDiv%theta
+thet(i)=(ThetaStart+(ThetaStop-ThetaStart)*(i-1)/(ThetaDiv-1))/pi*180;
+theta=thet(i)*pi/180;
+for ii=1:phiDiv;%phi
+phi=(phiStop/(phiDiv-1))*(ii-1);
+
+for ind=1:21;%x0
+xv=lol+ (ind-1)*(N_sig/10)*sigma_x;%height above grating
+
+par.phi=phi;
+par.Theta=theta;
+par.x0=xv;
+par.facet=1;
+
+[res1,res2,res3,res4]=QUx(par,50,200);
+	G1x=(res1*tan(alpha1) +1i*res2*tan(alpha1))/(el_mm*pi); 
+	G1z=(res1+1i*res2)/(el_mm*pi);
+	G1y=(res3*tan(alpha1)-1i*res4*tan(alpha1))/(el_mm*pi);     
+
+par.facet=2;  
+
+[res1,res2,res3,res4]=QUx(par,50,200);
+	G2x=(res1*tan(alpha2) +1i*res2*tan(alpha2))/(el_mm*pi);
+	G2z=(res1+1i*res2)/(el_mm*pi);
+	G2y=(res3*tan(alpha2)-1i*res4*tan(alpha2))/(el_mm*pi);
+   
+      Gx= G1x+ G2x;
+      Gy= G1y+G2y;
+      Gz= G1z+G2z;
+      G=[Gx Gy Gz];
+      
+nv=complex([ sin(theta)*cos(phi) sin(theta)*sin(phi) cos(theta)]);   %dir. unit vector
+nvp=complex([cos(theta)*cos(phi) cos(theta)*sin(phi) -sin(theta)]);%1st pol. unit vector
+nvv=complex([-sin(phi) cos(phi) 0]);%2nd pol. unit vector
+
+    a1t=(G*nvp')*conj(G*nvp');% square of magn. in first polarisation
+    a2t=(G*nvv')*conj(G*nvv');% square of magn. in second polarisation 
+    
+      f5c=a1t+a2t;              %  total 
+
+       if (ind==11)%
+
+       single_elec(i,ii)= real(f5c);
+       single_elec_p1(i,ii)= real(a1t);
+       single_elec_p2(i,ii)= real(a2t) ;   
+       end         %   /*store the single-electron R^2 */  
+   
+      x0_bin(ind)= xv;
+      inc_bin(ind)= real(f5c)*exp(-(xv-x0_mm)*(xv-x0_mm)/(2*sigma_x*sigma_x));
+      inc_bin_p1(ind)= real(a1t)*exp(-(xv-x0_mm)*(xv-x0_mm)/(2*sigma_x*sigma_x));
+      inc_bin_p2(ind)= real(a2t)*exp(-(xv-x0_mm)*(xv-x0_mm)/(2*sigma_x*sigma_x));
+      coh_bin(ind)= sqrt(real(f5c))*exp(-(xv-x0_mm)*(xv-x0_mm)/(2*sigma_x*sigma_x));
+      coh_bin_p1(ind)= sqrt(real(a1t))*exp(-(xv-x0_mm)*(xv-x0_mm)/(2*sigma_x*sigma_x));
+      coh_bin_p2(ind)= sqrt(real(a2t))*exp(-(xv-x0_mm)*(xv-x0_mm)/(2*sigma_x*sigma_x));
+end %of x0
+%  integral over hight from grating
+toler=1000;
+dist1=(x0_bin(end)-x0_bin(1))/toler;
+dist=x0_bin(1):dist1:x0_bin(end);
+incans=pchip(x0_bin,inc_bin,dist);
+s_inc_x(i,ii)=   trapz(dist,incans)/(sqrt(2*pi)*sigma_x);
+cohans=pchip(x0_bin,coh_bin,dist);
+s_coh_x(i,ii)= ( trapz(dist,cohans)/(sqrt(2*pi)*sigma_x))^2;
+
+incans=pchip(x0_bin,inc_bin_p1,dist);
+s_inc_x_p1(i,ii)=    trapz(dist,incans)/(sqrt(2*pi)*sigma_x);
+cohans=pchip(x0_bin,coh_bin_p1,dist);
+s_coh_x_p1(i,ii)= ( trapz(dist,cohans)/(sqrt(2*pi)*sigma_x))^2;
+
+incans=pchip(x0_bin,inc_bin_p2,dist);
+s_inc_x_p2(i,ii)=   trapz(dist,incans)/(sqrt(2*pi)*sigma_x);
+cohans=pchip(x0_bin,coh_bin_p2,dist);
+s_coh_x_p2(i,ii)= (trapz(dist,cohans)/(sqrt(2*pi)*sigma_x))^2;
+
+end%of phi
+end%of theta
+%% Form factor
+% form factor is calculated numericaly by FFT
+SpectrumFFT=abs(fft(Profile)).^2;%Fourier transform
+SpectrumHalf=SpectrumFFT(1:round(length(SpectrumFFT)/2));%only half of FT spectrum hold usefull information 
+Frequency=(0:(length(SpectrumFFT)-1))./(time(2)-time(1))/length(SpectrumFFT)*1e-3;%THz
+FreqHalf=Frequency(1:round(length(Frequency)/2));%half of frequency range
+SpectrumHalf=SpectrumHalf./max(SpectrumHalf);%normalization of Spectrum
+Nmax=find(FreqHalf>=2.5,1);%2.5 THz is more than enoght
+% Interpolation is applyed to decrease error in calculation (freq step in FFT is to wide for this calculation ) 
+FreqInter=FreqHalf(1):((FreqHalf(2)-FreqHalf(1))*1e-3):FreqHalf(Nmax);%interpolation diapason
+
+SpecInter=pchip( FreqHalf(1:Nmax),SpectrumHalf(1:Nmax),FreqInter);%interpolation
+if PlotFunc(1)==1
+figure(1)
+clf
+plot(FreqInter,SpecInter,'-k','linew',2)
+hold on
+grid on
+set(gca,'fontsize',16)
+ylabel('F(\nu)')
+xlabel('\nu, [THz]')
+legend('Form factor')
+end
+ %% I_N calculation
+f1=2*pi*qe*qe*1e-7;                                           % in Joules
+f2=100/(el_mm^2); 
+f4=(m^2)*(beta./(1-beta*cos(thet*pi/180))).^3;
+phi=(phiStop/(phiDiv-1))*((1:phiDiv)-1);
+%lam_e=lam*beta*gam/(2*pi*sqrt(1+(beta*gam*sin(theta)*sin(phi0))^2));
+lam= el_mm*(1/beta-cos(thet*pi/180))/m;%mm
+
+for i=1:ThetaDiv%theta
+for ii=1:phiDiv;%phi
+ky=(2*pi/lam(i))*sin(thet(i)*pi/180)*sin(phi(ii));
+om=2*pi*0.3/lam(i);                          %c=0.3 mm/ps, hence om is in ps^-1                         
+%calculate longitudinal coherence effects
+%Form factor   
+if FFTspec==1
+ s_coh_long=SpecInter(find(2*pi*FreqInter>=om,1));
+ if isempty(s_coh_long)
+     s_coh_long=0;%ACHTUNG!!!
+     error('fail to find Form Factor value')
+ end
+else
+    %this work with new versions of MATLAB (2015a for example)
+u1= om*sigma_t;
+u2= om*sigma_t*epsl;
+ f_z_c=(exp(-(u1*u1)/2)+ epsl*exp(-(u2*u2)/2))/(1+epsl);
+ f_z_s=-2/(sqrt(pi)*(1+epsl))*(dawson(u2/sqrt(2))*epsl-dawson(u1/sqrt(2)));
+ s_coh_long=f_z_c.^2+f_z_s.^2;
+end
+%%
+
+Y2=exp(-((ky*sigma_y)^2));       % Ysquare
+
+s_coh_xy= Y2*s_coh_x(i,ii);
+s_coh_xy_p1= Y2*s_coh_x_p1(i,ii);
+s_coh_xy_p2= Y2*s_coh_x_p2(i,ii);
+FFact(i)=s_coh_long;
+bin(i,1)=thet(i);
+bin(i,ii+1)= f1*f2*f4(i)*(n_e*s_inc_x(i,ii)+ n_e*n_e*s_coh_xy*s_coh_long);
+bin_p1(i,ii+1)= f1*f2*f4(i)*(n_e*s_inc_x_p1(i,ii)+ n_e*n_e*s_coh_xy_p1*s_coh_long);
+bin_p2(i,ii+1)= f1*f2*f4(i)*(n_e*s_inc_x_p2(i,ii)+ n_e*n_e*s_coh_xy_p2*s_coh_long);
+bin_single(i,1)= thet(i);                     % store for single electron yield 
+bin_single(i,ii+1)= f1*f2*f4(i)*single_elec(i,ii);
+bin_single_p1(i,ii+1)= f1*f2*f4(i)*single_elec_p1(i,ii);
+bin_single_p2(i,ii+1)= f1*f2*f4(i)*single_elec_p2(i,ii);
+end% of ii (phi)           
+end % of i (theta)
+
+if PlotFunc(1)==1
+figure(1)
+hold on
+plot(c./(lam*1e-3)*1e-12,FFact,'or','linew',2)
+legend('Form factor','Used points')
+end
+%%
+% At this point calculation of
+% 
+% $$W_N=W_1 |_{x_0=0}(NS_{inc}+N^2S_{coh}) $$
+%
+% is complete. But there is dependence from angles. So next step is
+% integration
+%% Integration and averaging
+philmt=phiDiv+1;%number of steps +1
+%phiDiv + 2 -- average over phi
+for i=1:ThetaDiv%theta
+bin(i,phiDiv+2)=(bin(i,2)+ 2*sum(bin(i,3:philmt)))/(2*(philmt-1)-1);    %    col. 9 is the store for AVERAGE dE over all ACCEPTED phi's; in J/sr/cm 	 
+bin_p1(i,phiDiv+2)=(bin_p1(i,2)+ 2*sum(bin_p1(i,3:philmt)))/(2*(philmt-1)-1);     %   ditto for pol. 1 
+bin_p2(i,phiDiv+2)=(bin_p2(i,2)+ 2*sum(bin_p2(i,3:philmt)))/(2*(philmt-1)-1);     %   and for pol. 2 
+bin_single(i,phiDiv+2)= (bin_single(i,2)+ 2*sum(bin_single(i,3:philmt)))/(2*(philmt-1)-1);
+bin_single_p1(i,phiDiv+2)= (bin_single_p1(i,2)+ 2*sum(bin_single_p1(i,3:philmt)))/(2*(philmt-1)-1);
+bin_single_p2(i,phiDiv+2)= (bin_single_p2(i,2)+ 2*sum(bin_single_p2(i,3:philmt)))/(2*(philmt-1)-1);
+end
+
+%phiDiv + 3 -- integral
+% integrate over phi 
+d_theta=(ThetaStop-ThetaStart)/(ThetaDiv-1);                             % d_theta expressed in rad 
+d_phi= (phiStop/(phiDiv-1));                  %d_phi expressed in rad 
+
+for i=1:ThetaDiv%theta
+
+     theta=(ThetaStart+(ThetaStop-ThetaStart)*(i-1)/(ThetaDiv-1));
+     
+     sum0 = sin(theta)*d_theta*d_phi*sum(bin(i,3:(phiDiv+1)));   
+     sum_p1 =sin(theta)*d_theta*d_phi*sum(bin_p1(i,3:(phiDiv+1)));   %pol. 1 only
+     sum_p2 =sin(theta)*d_theta*d_phi*sum(bin_p2(i,3:(phiDiv+1)));   %and pol. 2
+     sum_single =sin(theta)*d_theta*d_phi*sum(bin_single(i,3:(phiDiv+1)));   
+     sum_single_p1 = sin(theta)*d_theta*d_phi*sum(bin_single_p1(i,3:(phiDiv+1)));   
+     sum_single_p2 =  sin(theta)*d_theta*d_phi*sum(bin_single_p2(i,3:(phiDiv+1)));        
+     I3(i)=d_phi*(bin(i,2)+2*sum(bin(i,3:(phiDiv+1))));
+    bin(i,phiDiv+3)= sin(theta)*d_theta*d_phi*bin(i,2)+ 2*sum0;
+    
+    bin_p1(i,phiDiv+3) = sin(theta)*d_theta*d_phi*bin_p1(i,2)+ 2*sum_p1;%  /* pol. 1 only */
+    bin_p2(i,phiDiv+3) = sin(theta)*d_theta*d_phi*bin_p2(i,2)+ 2*sum_p2; % /*  and pol. 2 */
+    bin_single(i,phiDiv+3) = sin(theta)*d_theta*d_phi*bin_single(i,2) + 2*sum_single;    
+    bin_single_p1(i,phiDiv+3) = sin(theta)*d_theta*d_phi*bin_single_p1(i,2) + 2*sum_single_p1;
+    bin_single_p2(i,phiDiv+3) = sin(theta)*d_theta*d_phi*bin_single_p2(i,2) + 2*sum_single_p2;
+end
+  
+% now integrate over theta 
+ sum1=sum(bin(1:ThetaDiv,phiDiv+3));       
+ sum3=sum(bin_p1(1:ThetaDiv,phiDiv+3));    
+ sum4=sum(bin_p2(1:ThetaDiv,phiDiv+3));    
+ sum2=sum(bin_single(1:ThetaDiv,phiDiv+3));    
+
+disp(['Rough estimate of total radiated energy (single electron) over hemisphere =' num2str(sum2*0.1*lgr_mm) ' J']);%/* no need for x2, since +-phi has been taken into account*/
+disp(['Rough estimate of total radiated energy (bunch) over hemisphere =' num2str(sum1*0.1*lgr_mm) ' J']);%/* no need for x2, since +-phi has been taken into account*/
+percent=sum1*0.1*lgr_mm/(n_e_t*1.602177e-19*(gam-1)*0.511*1e+6);
+disp(['This is equal to ' num2str(percent) ' of the energy in the bunch'])
+disp(['Energy in polarization 1= ' num2str(sum3*0.1*lgr_mm) ' J'])
+disp(['Energy in polarization 2= ' num2str(sum4*0.1*lgr_mm) ' J'])
+lgr_e_mm = abs(Aed_mm./sin(thet/180*pi));           % variable effective grating length, determined basically by the cone entrance diameter of XXXmm*/
+lgr_e_mm(lgr_e_mm>lgr_mm)=lgr_mm;
+sum_tr = bin(:,phiDiv+2) .* 0.1 .* lgr_e_mm(:) .* SolAng;   % /*variable grating length,  now in J*/
+sum_tr_single = bin_single(:,phiDiv+2) .* 0.1 .* lgr_e_mm(:) * SolAng;
+sum_tr_p1 = bin_p1(:,phiDiv+2) .* 0.1 .* lgr_e_mm(:) * SolAng;
+sum_tr_p2 = bin_p2(:,phiDiv+2) .* 0.1 .* lgr_e_mm(:) * SolAng;
+sum_tr_single_p1 = bin_single_p1(:,phiDiv+2) .* 0.1 .* lgr_e_mm(:) * SolAng;
+sum_tr_single_p2 = bin_single_p2(:,phiDiv+2) .* 0.1 .* lgr_e_mm(:) * SolAng;
+ 
+% KV. Precise calculation over phi
+phiNew=phi(1):((phi(2)-phi(1))*1e-3):phi(end);
+for i=1:ThetaDiv%theta
+        IntergalPhi=sum(pchip(phi,single_elec(i,:),phiNew))*(phiNew(2)-phiNew(1));
+        I1(i)=f1*f2*f4(i)*IntergalPhi;%interal over phi (0 to phiMax)
+         IntergalPhi=sum(pchip(phi,bin(i,2:(phiDiv+1)),phiNew))*(phiNew(2)-phiNew(1));
+        In(i)=IntergalPhi;%interal over phi (0 to phiMax)
+end
+E_total=sum(In)*d_theta*Aed_mm*0.1;%J
+nu_Hz=c./(lam*1e-3);%Hz
+SP_spec_1=  I1(:).*0.1.*lgr_e_mm(:)./(el_mm/m*1e-3)*2.*(c./nu_Hz(:).^2)*1e+18;%muJ/THz
+SP_spec_N=In(:) .*0.1.*lgr_e_mm(:)./(el_mm/m*1e-3)*2.*(c./nu_Hz(:).^2)*1e+18;%muJ/THz
+%% Plots
+if PlotFunc(2)==1
+ figure(2)
+ clf
+ plot(thet,bin_single(:,phiDiv+2),'k','linew',2)
+ set(gca,'fontsize',16)
+ylabel('dI_1/d\Omega, [J sr^{-1}cm^{-1}]')
+xlabel('\Theta, [deg]')
+title('Singe electron yield')
+ grid on
+end
+if PlotFunc(3)==1
+figure(3)
+clf
+ plot(thet,sum_tr,'r','linew',2)
+% thet_rad=thet/180*pi;
+% polar(thet_rad(:),sum_tr(:));
+% ylim([0 Inf])
+%  yl = get(gca,'YLim');
+% set(gca,'YLim', [0 yl(2)]);
+ set(gca,'fontsize',16)
+ylabel('dI/sin(\Theta)d\Theta, [J]')
+xlabel('\Theta, [deg]')
+title('SP spectrum')
+ grid on
+end
+
+if PlotFunc(4)==1
+figure(4)
+clf
+plot(time,Profile,'r','linew',2)
+ set(gca,'fontsize',16)
+ylabel('Amplitude, [a.u.]')
+xlabel('t, [ns]')
+title('Temporal bunch profile')
+ grid on
+end
+
+if PlotFunc(5)==1
+figure(5)
+clf
+semilogy(nu_Hz*1e-12,SP_spec_1,'r','linew',2)
+ hold on
+%  load CTR.mat
+% semilogy(nu*1e-12,IntergdOm*1e+18,'k','linew',2)
+set(gca,'fontsize',16)
+ylabel('Energy, [\muJ/THz]')
+xlabel('\nu, [THz]')
+title('SP spectrum (Single electron)')
+% legend('SP','CTR')
+ grid on
+xlim([min(nu_Hz) max(nu_Hz)].*1e-12)
+end
+
+
+
Index: /CLIO/Analysis/Data19July/GFW/GFWnewBKUP.m
===================================================================
--- /CLIO/Analysis/Data19July/GFW/GFWnewBKUP.m	(revision 723)
+++ /CLIO/Analysis/Data19July/GFW/GFWnewBKUP.m	(revision 723)
@@ -0,0 +1,422 @@
+  function  [nu_Hz,SP_spec_1,SP_spec_N,lam,E_total,sum_tr,thet,bin_single]=GFWnew(var)
+
+%  single_elec
+% close all; clear all; clc
+%% Configuration of calculation
+%PrSelect == what profile to use: 1 --from file; 2 - assymmetric gaussian
+%P_file -- name and path to file generated with ASTRA (use TakeData.m script on svn=> CLIO/Scripts to take data from ASTRA generated files)
+%PlotFunc:
+%   1 -- Form factor
+%   2 -- Single electron yield from Theta
+%   3 -- dI_N/dOmega (SP spetrum for whole bunch)
+%   4 -- Bunch profile
+%   5 -- dI(nu)/dnu (SP spetrum for signle elctron bunch)
+% var='PrSelect=1; P_file=''../p20.mat''; PlotFunc=[1 1 1 1 1]; ';
+PrSelect=2;
+PlotFunc=[0 0 0 0 0];
+P_file='tmp.mat';
+ %% Physical constants
+c=299792458;%m/s speed of ligth
+qe = 1.602176487e-19*c*10;%Electron charge, statColumns
+%% Parameters block (by default)
+T=1;%Transmissivity (from gun exit to experimental position)
+m=1;%Order of radiation
+lgr_mm=180;%Grating length, mm
+wid_mm=40;%Grating width, mm
+el_mm=6;%Grating period, mm
+alpha1=30;%Blaze angle, degrees
+gam=100;%relativistic gamma
+Qbunch_nC=1.2;%Bunch charge, nC (if T=1,than chagre at experiment position, if not, than charge at the exit of gun)
+x0_mm=8; %Position of beam centroid above grating, mm
+fwhm_x_mm=3.07; %FWHM_X  of Gaussian beam at the centre of grating, mm
+fwhm_y_mm=3.07; %FWHM_Y  of Gaussian beam at the centre of grating, mm
+pulse_len_ps=10;%Nominal bunch duration (FWHM), ps 
+epsl=1;%Asymmetry factor for assymmetric Gaussian
+BAd_mm=410;% beam - aperture distance ,  mm
+Aed_mm=25;%Aperture entrance diameter,  mm 
+%All angles are in radians
+phiStop=atan(Aed_mm/2/BAd_mm);% phi go from 0 to phiStop
+% phiStop=7/180*pi;
+% ThetaStart=pi/2-atan(Aed_mm/2/BAd_mm);
+% ThetaStop=pi/2+atan(Aed_mm/2/BAd_mm);
+ThetaStart=alpha1/180*pi;%initial theta angle
+ThetaStop=170/180*pi;%final theta angle
+% ThetaStart=80/180*pi;%initial theta angle
+% ThetaStop=140/180*pi;%final theta angle
+ThetaDiv=10;%steps in theta
+phiDiv=3;%steps in phi
+
+if (~isempty(var))
+strread(var,'%s','delimiter',' ')
+eval(var);
+end %if
+
+%save all parameters
+name=regexprep(var,'.mat','');
+name=regexprep(name,'[../''[]]','');
+name=regexprep(name,'[; ]','_');
+% save([ name 'SP.mat'])
+%% Profile selection
+if PrSelect==1
+    FFTspec=1;
+  if  exist(P_file,'file')==2;
+    load([P_file]);
+  else
+      error('Set correct file name!')
+  end
+else
+    FFTspec=1;%how to calculate form factor: 1 -- FFT; 2 -- only fro assymetric gaussian
+%Gaussian
+sigma_t= pulse_len_ps/(sqrt(2*log(2))*(1+epsl))*1e-3;  %The sigma of the leading edge    ns
+
+nu=(1:1e+10:3e+16).*1e-9;%Ghz
+
+sigma_L= sigma_t*epsl;
+sigma_H= sigma_t;
+time=(1:(length(nu)*2))*(1/(length(nu)*2*(nu(2)-nu(1))));
+xc=mean(time);
+G1=exp(-(time-xc).^2./(2*sigma_L^2));
+G2=exp(-(time-xc).^2./(2*sigma_H^2));
+G1(time>xc)=0;
+G2(time<xc)=0;
+Profile=(G1+G2)./sum(G1+G2);
+end
+
+%% Calculation block for additional parameters
+ n_e_t=round(Qbunch_nC*1e-9/1.6e-19*T);              %Total number of electrons in the bunch
+% n_e_t=1e+11;
+SolAng=pi*((Aed_mm/2)^2)/(BAd_mm^2); %solid angle from beam to aperture
+% Grating
+alpha1=alpha1/180*pi;
+el1=el_mm*(cos(alpha1)^2);
+el2=el_mm-el1;
+hi=abs(el1*tan(alpha1));
+alpha2=-(pi/2-alpha1);
+% Beam 
+beta=sqrt(gam^2 - 1)/gam;
+sigma_x= fwhm_x_mm/(2*sqrt(2*log(2)));
+sigma_y= fwhm_y_mm/(2*sqrt(2*log(2)));
+% N_e above grating
+N_sig=3;
+lol=x0_mm-N_sig*sigma_x;                                       
+if (lol<hi) lol=1.001*hi; end
+upl=x0_mm+N_sig*sigma_x;
+GAU = @(x) exp(-(x-x0_mm).*(x-x0_mm)/(2*sigma_x*sigma_x));
+res1 = integral(GAU,lol,upl)/(sigma_x*sqrt(2*pi));%part of the beam is above the grating
+disp([num2str(res1) ' of the beam above the grating'])
+n_e=n_e_t*res1;                                                 %Number of electrons above the grating
+% save([ name '_additional_SP.mat'])
+%% Estimation grating factor R^2, coherent & incoherent integrals
+%%
+% 
+% $$W_1 |_{x_0=0}= \sum \frac{q^2\omega^3 L l }{4\pi^2|m|c^3}| n  \times n \times \sum G(...)|^2$$
+% 
+% In this block we calculate 
+%%
+%
+% $$| n  \times n \times \sum G(...)|^2 $$
+%
+% where other terms (except L - grating length) are in f1, f2, f4. In general:
+%%
+%
+% $$W_N=W_1 |_{x_0=0}(NS_{inc}+N^2S_{coh}) $$
+%
+% where s_coh and s_inc are:
+%%
+%
+% $$S_{inc}=\int_h^{\infty}dxX(x)$$
+%
+% $$S_{coh}=|\int_h^{\infty}dxX(x)\tilde{Y}(ky)\tilde{T}(\omega)|^2 $$
+%
+% for now (in this block)
+%%
+%
+% $$S^*_{inc}=\int_h^{\infty}dxX(x)*| n  \times n \times \sum G(...)|^2\\ $$
+%
+% $$S^*_{coh}=\int_h^{\infty}dxX(x)*| n  \times n \times \sum G(...)|  $$
+%
+par = struct('Theta',NaN,'phi',NaN,'m',m,'alpha1',alpha1,'facet',NaN,'gamma',gam,'x0',x0_mm,'l',el_mm,'w',wid_mm);
+for i=1:ThetaDiv%theta
+thet(i)=(ThetaStart+(ThetaStop-ThetaStart)*(i-1)/(ThetaDiv-1))/pi*180;
+theta=thet(i)*pi/180;
+for ii=1:phiDiv;%phi
+phi=(phiStop/(phiDiv-1))*(ii-1);
+
+for ind=1:21;%x0
+xv=lol+ (ind-1)*(N_sig/10)*sigma_x;%height above grating
+
+par.phi=phi;
+par.Theta=theta;
+par.x0=xv;
+par.facet=1;
+
+[res1,res2,res3,res4]=QUx(par,50,200);
+	G1x=(res1*tan(alpha1) +1i*res2*tan(alpha1))/(el_mm*pi); 
+	G1z=(res1+1i*res2)/(el_mm*pi);
+	G1y=(res3*tan(alpha1)-1i*res4*tan(alpha1))/(el_mm*pi);     
+
+par.facet=2;  
+
+[res1,res2,res3,res4]=QUx(par,50,200);
+	G2x=(res1*tan(alpha2) +1i*res2*tan(alpha2))/(el_mm*pi);
+	G2z=(res1+1i*res2)/(el_mm*pi);
+	G2y=(res3*tan(alpha2)-1i*res4*tan(alpha2))/(el_mm*pi);
+   
+      Gx= G1x+ G2x;
+      Gy= G1y+G2y;
+      Gz= G1z+G2z;
+      G=[Gx Gy Gz];
+      
+nv=complex([ sin(theta)*cos(phi) sin(theta)*sin(phi) cos(theta)]);   %dir. unit vector
+nvp=complex([cos(theta)*cos(phi) cos(theta)*sin(phi) -sin(theta)]);%1st pol. unit vector
+nvv=complex([-sin(phi) cos(phi) 0]);%2nd pol. unit vector
+
+    a1t=(G*nvp')*conj(G*nvp');% square of magn. in first polarisation
+    a2t=(G*nvv')*conj(G*nvv');% square of magn. in second polarisation 
+    
+      f5c=a1t+a2t;              %  total 
+
+       if (ind==11)%
+
+       single_elec(i,ii)= real(f5c);
+       single_elec_p1(i,ii)= real(a1t);
+       single_elec_p2(i,ii)= real(a2t) ;   
+       end         %   /*store the single-electron R^2 */  
+   
+      x0_bin(ind)= xv;
+      inc_bin(ind)= real(f5c)*exp(-(xv-x0_mm)*(xv-x0_mm)/(2*sigma_x*sigma_x));
+      inc_bin_p1(ind)= real(a1t)*exp(-(xv-x0_mm)*(xv-x0_mm)/(2*sigma_x*sigma_x));
+      inc_bin_p2(ind)= real(a2t)*exp(-(xv-x0_mm)*(xv-x0_mm)/(2*sigma_x*sigma_x));
+      coh_bin(ind)= sqrt(real(f5c))*exp(-(xv-x0_mm)*(xv-x0_mm)/(2*sigma_x*sigma_x));
+      coh_bin_p1(ind)= sqrt(real(a1t))*exp(-(xv-x0_mm)*(xv-x0_mm)/(2*sigma_x*sigma_x));
+      coh_bin_p2(ind)= sqrt(real(a2t))*exp(-(xv-x0_mm)*(xv-x0_mm)/(2*sigma_x*sigma_x));
+end %of x0
+%  integral over hight from grating
+toler=1000;
+dist1=(x0_bin(end)-x0_bin(1))/toler;
+dist=x0_bin(1):dist1:x0_bin(end);
+incans=pchip(x0_bin,inc_bin,dist);
+s_inc_x(i,ii)=   trapz(dist,incans)/(sqrt(2*pi)*sigma_x);
+cohans=pchip(x0_bin,coh_bin,dist);
+s_coh_x(i,ii)= ( trapz(dist,cohans)/(sqrt(2*pi)*sigma_x))^2;
+
+incans=pchip(x0_bin,inc_bin_p1,dist);
+s_inc_x_p1(i,ii)=    trapz(dist,incans)/(sqrt(2*pi)*sigma_x);
+cohans=pchip(x0_bin,coh_bin_p1,dist);
+s_coh_x_p1(i,ii)= ( trapz(dist,cohans)/(sqrt(2*pi)*sigma_x))^2;
+
+incans=pchip(x0_bin,inc_bin_p2,dist);
+s_inc_x_p2(i,ii)=   trapz(dist,incans)/(sqrt(2*pi)*sigma_x);
+cohans=pchip(x0_bin,coh_bin_p2,dist);
+s_coh_x_p2(i,ii)= (trapz(dist,cohans)/(sqrt(2*pi)*sigma_x))^2;
+
+end%of phi
+end%of theta
+%% Form factor
+% form factor is calculated numericaly by FFT
+SpectrumFFT=abs(fft(Profile)).^2;%Fourier transform
+SpectrumHalf=SpectrumFFT(1:round(length(SpectrumFFT)/2));%only half of FT spectrum hold usefull information 
+Frequency=(0:(length(SpectrumFFT)-1))./(time(2)-time(1))/length(SpectrumFFT)*1e-3;%THz
+FreqHalf=Frequency(1:round(length(Frequency)/2));%half of frequency range
+SpectrumHalf=SpectrumHalf./max(SpectrumHalf);%normalization of Spectrum
+Nmax=find(FreqHalf>=60,1);%2.5 THz is more than enoght
+% Interpolation is applyed to decrease error in calculation (freq step in FFT is to wide for this calculation ) 
+FreqInter=FreqHalf(1):((FreqHalf(2)-FreqHalf(1))*1e-3):FreqHalf(Nmax);%interpolation diapason
+SpecInter=pchip( FreqHalf(1:Nmax),SpectrumHalf(1:Nmax),FreqInter);%interpolation
+if PlotFunc(1)==1
+figure(1)
+clf
+plot(FreqInter,SpecInter,'-k','linew',2)
+hold on
+grid on
+set(gca,'fontsize',16)
+ylabel('F(\nu)')
+xlabel('\nu, [THz]')
+legend('Form factor')
+end
+ %% I_N calculation
+f1=2*pi*qe*qe*1e-7;                                           % in Joules
+f2=100/(el_mm^2); 
+f4=(m^2)*(beta./(1-beta*cos(thet*pi/180))).^3;
+phi=(phiStop/(phiDiv-1))*((1:phiDiv)-1);
+%lam_e=lam*beta*gam/(2*pi*sqrt(1+(beta*gam*sin(theta)*sin(phi0))^2));
+lam= el_mm*(1/beta-cos(thet*pi/180))/m;%mm
+
+for i=1:ThetaDiv%theta
+for ii=1:phiDiv;%phi
+ky=(2*pi/lam(i))*sin(thet(i)*pi/180)*sin(phi(ii));
+om=2*pi*0.3/lam(i);                          %c=0.3 mm/ps, hence om is in ps^-1                         
+%calculate longitudinal coherence effects
+%Form factor   
+if FFTspec==1
+ s_coh_long=SpecInter(find(2*pi*FreqInter>=om,1));
+ if isempty(s_coh_long)
+     s_coh_long=0;%ACHTUNG!!!
+     error('fail to find Form Factor value')
+ end
+else
+    %this work with new versions of MATLAB (2015a for example)
+u1= om*sigma_t;
+u2= om*sigma_t*epsl;
+ f_z_c=(exp(-(u1*u1)/2)+ epsl*exp(-(u2*u2)/2))/(1+epsl);
+ f_z_s=-2/(sqrt(pi)*(1+epsl))*(dawson(u2/sqrt(2))*epsl-dawson(u1/sqrt(2)));
+ s_coh_long=f_z_c.^2+f_z_s.^2;
+end
+%%
+
+Y2=exp(-((ky*sigma_y)^2));       % Ysquare
+
+s_coh_xy= Y2*s_coh_x(i,ii);
+s_coh_xy_p1= Y2*s_coh_x_p1(i,ii);
+s_coh_xy_p2= Y2*s_coh_x_p2(i,ii);
+FFact(i)=s_coh_long;
+bin(i,1)=thet(i);
+bin(i,ii+1)= f1*f2*f4(i)*(n_e*s_inc_x(i,ii)+ n_e*n_e*s_coh_xy*s_coh_long);
+bin_p1(i,ii+1)= f1*f2*f4(i)*(n_e*s_inc_x_p1(i,ii)+ n_e*n_e*s_coh_xy_p1*s_coh_long);
+bin_p2(i,ii+1)= f1*f2*f4(i)*(n_e*s_inc_x_p2(i,ii)+ n_e*n_e*s_coh_xy_p2*s_coh_long);
+bin_single(i,1)= thet(i);                     % store for single electron yield 
+bin_single(i,ii+1)= f1*f2*f4(i)*single_elec(i,ii);
+bin_single_p1(i,ii+1)= f1*f2*f4(i)*single_elec_p1(i,ii);
+bin_single_p2(i,ii+1)= f1*f2*f4(i)*single_elec_p2(i,ii);
+end% of ii (phi)           
+end % of i (theta)
+
+if PlotFunc(1)==1
+figure(1)
+hold on
+plot(c./(lam*1e-3)*1e-12,FFact,'or','linew',2)
+legend('Form factor','Used points')
+end
+%%
+% At this point calculation of
+% 
+% $$W_N=W_1 |_{x_0=0}(NS_{inc}+N^2S_{coh}) $$
+%
+% is complete. But there is dependence from angles. So next step is
+% integration
+%% Integration and averaging
+philmt=phiDiv+1;%number of steps +1
+%phiDiv + 2 -- average over phi
+for i=1:ThetaDiv%theta
+bin(i,phiDiv+2)=(bin(i,2)+ 2*sum(bin(i,3:philmt)))/(2*(philmt-1)-1);    %    col. 9 is the store for AVERAGE dE over all ACCEPTED phi's; in J/sr/cm 	 
+bin_p1(i,phiDiv+2)=(bin_p1(i,2)+ 2*sum(bin_p1(i,3:philmt)))/(2*(philmt-1)-1);     %   ditto for pol. 1 
+bin_p2(i,phiDiv+2)=(bin_p2(i,2)+ 2*sum(bin_p2(i,3:philmt)))/(2*(philmt-1)-1);     %   and for pol. 2 
+bin_single(i,phiDiv+2)= (bin_single(i,2)+ 2*sum(bin_single(i,3:philmt)))/(2*(philmt-1)-1);
+bin_single_p1(i,phiDiv+2)= (bin_single_p1(i,2)+ 2*sum(bin_single_p1(i,3:philmt)))/(2*(philmt-1)-1);
+bin_single_p2(i,phiDiv+2)= (bin_single_p2(i,2)+ 2*sum(bin_single_p2(i,3:philmt)))/(2*(philmt-1)-1);
+end
+
+%phiDiv + 3 -- integral
+% integrate over phi 
+d_theta=(ThetaStop-ThetaStart)/(ThetaDiv-1);                             % d_theta expressed in rad 
+d_phi= (phiStop/(phiDiv-1));                  %d_phi expressed in rad 
+
+for i=1:ThetaDiv%theta
+
+     theta=(ThetaStart+(ThetaStop-ThetaStart)*(i-1)/(ThetaDiv-1));
+     
+     sum0 = sin(theta)*d_theta*d_phi*sum(bin(i,3:(phiDiv+1)));   
+     sum_p1 =sin(theta)*d_theta*d_phi*sum(bin_p1(i,3:(phiDiv+1)));   %pol. 1 only
+     sum_p2 =sin(theta)*d_theta*d_phi*sum(bin_p2(i,3:(phiDiv+1)));   %and pol. 2
+     sum_single =sin(theta)*d_theta*d_phi*sum(bin_single(i,3:(phiDiv+1)));   
+     sum_single_p1 = sin(theta)*d_theta*d_phi*sum(bin_single_p1(i,3:(phiDiv+1)));   
+     sum_single_p2 =  sin(theta)*d_theta*d_phi*sum(bin_single_p2(i,3:(phiDiv+1)));        
+     I3(i)=d_phi*(bin(i,2)+2*sum(bin(i,3:(phiDiv+1))));
+    bin(i,phiDiv+3)= sin(theta)*d_theta*d_phi*bin(i,2)+ 2*sum0;
+    
+    bin_p1(i,phiDiv+3) = sin(theta)*d_theta*d_phi*bin_p1(i,2)+ 2*sum_p1;%  /* pol. 1 only */
+    bin_p2(i,phiDiv+3) = sin(theta)*d_theta*d_phi*bin_p2(i,2)+ 2*sum_p2; % /*  and pol. 2 */
+    bin_single(i,phiDiv+3) = sin(theta)*d_theta*d_phi*bin_single(i,2) + 2*sum_single;    
+    bin_single_p1(i,phiDiv+3) = sin(theta)*d_theta*d_phi*bin_single_p1(i,2) + 2*sum_single_p1;
+    bin_single_p2(i,phiDiv+3) = sin(theta)*d_theta*d_phi*bin_single_p2(i,2) + 2*sum_single_p2;
+end
+  
+% now integrate over theta 
+ sum1=sum(bin(1:ThetaDiv,phiDiv+3));       
+ sum3=sum(bin_p1(1:ThetaDiv,phiDiv+3));    
+ sum4=sum(bin_p2(1:ThetaDiv,phiDiv+3));    
+ sum2=sum(bin_single(1:ThetaDiv,phiDiv+3));    
+
+disp(['Rough estimate of total radiated energy (single electron) over hemisphere =' num2str(sum2*0.1*lgr_mm) ' J']);%/* no need for x2, since +-phi has been taken into account*/
+disp(['Rough estimate of total radiated energy (bunch) over hemisphere =' num2str(sum1*0.1*lgr_mm) ' J']);%/* no need for x2, since +-phi has been taken into account*/
+percent=sum1*0.1*lgr_mm/(n_e_t*1.602177e-19*(gam-1)*0.511*1e+6);
+disp(['This is equal to ' num2str(percent) ' of the energy in the bunch'])
+disp(['Energy in polarization 1= ' num2str(sum3*0.1*lgr_mm) ' J'])
+disp(['Energy in polarization 2= ' num2str(sum4*0.1*lgr_mm) ' J'])
+lgr_e_mm = abs(Aed_mm./sin(thet/180*pi));           % variable effective grating length, determined basically by the cone entrance diameter of XXXmm*/
+lgr_e_mm(lgr_e_mm>lgr_mm)=lgr_mm;
+sum_tr = bin(:,phiDiv+2) .* 0.1 .* lgr_e_mm(:) .* SolAng;   % /*variable grating length,  now in J*/
+sum_tr_single = bin_single(:,phiDiv+2) .* 0.1 .* lgr_e_mm(:) * SolAng;
+sum_tr_p1 = bin_p1(:,phiDiv+2) .* 0.1 .* lgr_e_mm(:) * SolAng;
+sum_tr_p2 = bin_p2(:,phiDiv+2) .* 0.1 .* lgr_e_mm(:) * SolAng;
+sum_tr_single_p1 = bin_single_p1(:,phiDiv+2) .* 0.1 .* lgr_e_mm(:) * SolAng;
+sum_tr_single_p2 = bin_single_p2(:,phiDiv+2) .* 0.1 .* lgr_e_mm(:) * SolAng;
+ 
+% KV. Precise calculation over phi
+phiNew=phi(1):((phi(2)-phi(1))*1e-3):phi(end);
+for i=1:ThetaDiv%theta
+        IntergalPhi=sum(pchip(phi,single_elec(i,:),phiNew))*(phiNew(2)-phiNew(1));
+        I1(i)=f1*f2*f4(i)*IntergalPhi;%interal over phi (0 to phiMax)
+         IntergalPhi=sum(pchip(phi,bin(i,2:(phiDiv+1)),phiNew))*(phiNew(2)-phiNew(1));
+        In(i)=IntergalPhi;%interal over phi (0 to phiMax)
+end
+E_total=sum(In)*d_theta*Aed_mm*0.1;%J
+nu_Hz=c./(lam*1e-3);%Hz
+SP_spec_1=  I1(:).*0.1.*lgr_e_mm(:)./(el_mm/m*1e-3)*2.*(c./nu_Hz(:).^2)*1e+18;%muJ/THz
+SP_spec_N=In(:) .*0.1.*lgr_e_mm(:)./(el_mm/m*1e-3)*2.*(c./nu_Hz(:).^2)*1e+18;%muJ/THz
+%% Plots
+if PlotFunc(2)==1
+ figure(2)
+ clf
+ plot(thet,bin_single(:,phiDiv+2),'k','linew',2)
+ set(gca,'fontsize',16)
+ylabel('dI_1/d\Omega, [J sr^{-1}cm^{-1}]')
+xlabel('\Theta, [deg]')
+title('Singe electron yield')
+ grid on
+end
+if PlotFunc(3)==1
+figure(3)
+clf
+ plot(thet,sum_tr,'r','linew',2)
+% thet_rad=thet/180*pi;
+% polar(thet_rad(:),sum_tr(:));
+% ylim([0 Inf])
+%  yl = get(gca,'YLim');
+% set(gca,'YLim', [0 yl(2)]);
+ set(gca,'fontsize',16)
+ylabel('dI/sin(\Theta)d\Theta, [J]')
+xlabel('\Theta, [deg]')
+title('SP spectrum')
+ grid on
+end
+
+if PlotFunc(4)==1
+figure(4)
+clf
+plot(time,Profile,'r','linew',2)
+ set(gca,'fontsize',16)
+ylabel('Amplitude, [a.u.]')
+xlabel('t, [ns]')
+title('Temporal bunch profile')
+ grid on
+end
+
+if PlotFunc(5)==1
+figure(5)
+clf
+semilogy(nu_Hz*1e-12,SP_spec_1,'r','linew',2)
+ hold on
+%  load CTR.mat
+% semilogy(nu*1e-12,IntergdOm*1e+18,'k','linew',2)
+set(gca,'fontsize',16)
+ylabel('Energy, [\muJ/THz]')
+xlabel('\nu, [THz]')
+title('SP spectrum (Single electron)')
+% legend('SP','CTR')
+ grid on
+xlim([min(nu_Hz) max(nu_Hz)].*1e-12)
+end
+
+
+
Index: /CLIO/Analysis/Data19July/GFW/QUx.m
===================================================================
--- /CLIO/Analysis/Data19July/GFW/QUx.m	(revision 723)
+++ /CLIO/Analysis/Data19July/GFW/QUx.m	(revision 723)
@@ -0,0 +1,40 @@
+function [rez1,rez2,rez3,rez4]=QUx(par,Nstepsz,Nstepsy)
+Theta=par.Theta;
+phi=par.phi;
+m=par.m;
+alpha1=par.alpha1;
+facet=par.facet;
+gamma=par.gamma;
+x0=par.x0;
+l=par.l;
+w=par.w;
+%%%%%%%%%%%%%
+alpha2=-(pi/2-alpha1);%add minus
+l1=l*(cos(alpha1)^2);
+beta=sqrt(gamma^2 - 1)/gamma;
+lambda=l*(1/beta-cos(Theta))/m;
+k=2*pi/lambda;
+zs=[0 l1];
+zf=[l1 l];
+ky=k*sin(Theta)*sin(phi);
+kx=k*sin(Theta)*cos(phi);
+alpha=[alpha1 alpha2];
+D=k*(1/beta-cos(Theta))-kx*tan(alpha(facet));
+chi=[0 l*tan(alpha(2))];
+Psi=[0 2*pi*l1/l];%change phase
+dz=(zf(facet)-zs(facet))/(Nstepsz-1);
+dy=w/2/(Nstepsy-1);
+z=repmat(zs(facet):dz:zf(facet),Nstepsy,1)';
+y=repmat(0:dy:(w/2),Nstepsz,1);
+b=x0-z.*tan(alpha(facet))+chi(facet);
+
+QU1=2*dz*dy*cos(D.*z-kx*chi(facet)+Psi(facet)).*k.*b./(beta*gamma*sqrt(b.^2+y.^2)).*cos(ky.*y).*besselk(1,k.*sqrt(b.^2+y.^2)./(beta.*gamma));
+QU2=2*dz*dy*sin(D.*z-kx*chi(facet)+Psi(facet)).*k.*b./(beta*gamma*sqrt(b.^2+y.^2)).*cos(ky.*y).*besselk(1,k.*sqrt(b.^2+y.^2)./(beta.*gamma)); 
+QU3=2*dz*dy*sin(D.*z-kx*chi(facet)+Psi(facet)).*k.*y./(beta*gamma*sqrt(b.^2+y.^2)).*sin(ky.*y).*besselk(1,k.*sqrt(b.^2+y.^2)./(beta.*gamma));
+QU4=2*dz*dy*cos(D.*z-kx*chi(facet)+Psi(facet)).*k.*y./(beta*gamma*sqrt(b.^2+y.^2)).*sin(ky.*y).*besselk(1,k.*sqrt(b.^2+y.^2)./(beta.*gamma));  %minus of not????
+
+rez1=trapz(trapz(QU1));
+rez2=trapz(trapz(QU2));
+rez3=trapz(trapz(QU3));
+rez4=trapz(trapz(QU4));
+
Index: /CLIO/Analysis/Data19July/GFW/START_CLIO.m
===================================================================
--- /CLIO/Analysis/Data19July/GFW/START_CLIO.m	(revision 723)
+++ /CLIO/Analysis/Data19July/GFW/START_CLIO.m	(revision 723)
@@ -0,0 +1,6 @@
+clear all
+ gam=44.2/0.511;
+%  [nu_Hz,SP_spec_1,SP_spec_N,lam,E_total,sum_tr,thet,bin_single, bin]=GFWnew(['gam=' num2str(gam) '; x0_mm=15; pulse_len_ps=5;  PlotFunc=[0 0 0 0 0]; ThetaStart=30/180*pi; ThetaStop=140/180*pi; ThetaDiv=111; phiDiv=13;'])
+%   save CLIOP15.mat nu_Hz SP_spec_1 SP_spec_N lam E_total sum_tr thet bin_single bin
+ [nu_Hz,SP_spec_1,SP_spec_N,lam,E_total,sum_tr,thet,bin_single, bin]=GFWnew(['PrSelect=1; P_file=''../P6.mat''; gam=' num2str(gam) '; x0_mm=10; pulse_len_ps=7;  PlotFunc=[0 0 0 0 0]; ThetaStart=30/180*pi; ThetaStop=140/180*pi; ThetaDiv=111; phiDiv=13;'])
+  save CLIOSpec.mat nu_Hz SP_spec_1 SP_spec_N lam E_total sum_tr thet bin_single bin
Index: /CLIO/Analysis/Data19July/GlobalRead.m
===================================================================
--- /CLIO/Analysis/Data19July/GlobalRead.m	(revision 723)
+++ /CLIO/Analysis/Data19July/GlobalRead.m	(revision 723)
@@ -0,0 +1,22 @@
+clear all; close all;clc
+
+path='clio24/';
+ listing = dir([path '*.sig']);
+ Filelist={listing.name}';
+  name=char(Filelist(:));
+for c=1:(size(name,1)-1)
+    if (name(c,1:9)~='data_CLIO')
+        name(c,:)=[];
+    else
+    try
+    [key_name, key_value(c,:), data(c,:,:)]=readsig([path name(c,:)]);
+    
+    %Hour(c,:)=name(c,22:29);
+    catch
+        disp(['Bad try on ' name(c,:)])
+    end
+    disp(['Percent done=' num2str((c)/(size(name,1))*100)]);
+    end
+end
+
+save('data1907.mat', 'key_name',  'key_value', 'data')
Index: /CLIO/Analysis/Data19July/KK/GetProfile.m
===================================================================
--- /CLIO/Analysis/Data19July/KK/GetProfile.m	(revision 723)
+++ /CLIO/Analysis/Data19July/KK/GetProfile.m	(revision 723)
@@ -0,0 +1,13 @@
+function [RPH]=GetProfile(SP,SS)
+%Restore profiles
+[RSP(s,:),ISP]=SRf(SP,SS);
+%Kramers-Kronig restore
+% PhaseK=KK_phase(ISP);
+% RPK(s,:)=real((ifft(RSP(s,:).*exp(1i.*(PhaseK)))));
+% [ChiVk(s),Corrk(s),RPK(s,:)]=PositionDirection(RPK(s,:),RPH(s,:));
+%Hilbert restore
+PhaseH=-(imag(hilbert(log(abs(RSP(s,:))))));
+RPH(s,:)=real((ifft(RSP(s,:).*exp(1i.*(PhaseH)))));
+
+
+
Index: /CLIO/Analysis/Data19July/KK/KK_phase.m
===================================================================
--- /CLIO/Analysis/Data19July/KK/KK_phase.m	(revision 723)
+++ /CLIO/Analysis/Data19July/KK/KK_phase.m	(revision 723)
@@ -0,0 +1,24 @@
+ function theta=KK_phase(Pt)
+ len=length(Pt);
+ n=1: len;
+ 
+for m=1: len
+q=2*m/pi*log(Pt(n)/Pt(m));
+g=(m.*m-n.*n)';
+lm=logical(ones([1  len]));
+lm(m)=0;
+Th(m)=q(lm)*(1./g(lm));
+%     if (m~=1)&&(m~=len)
+%         Th(m)=Th(m)+(Th(m-1)+Th(m+1))/2;
+%     end
+end
+% lx=2:(len-1);
+% Th(lx)=Th(lx)+(Th(lx-1)+Th(lx+1))/2;
+
+ theta=[-Th,0,fliplr(Th(2:length(Th)))];
+ 
+% m=1: len;
+% q=2/pi*repmat(m,len,1)*log((Pt(n)')*(1./Pt(m)));
+% g=(repmat(m.*m,len,1)-repmat((n.*n)',1,len));
+% lm=~logical(diag(ones([1  len])));
+% Th(m)=q(lm)*(1./g(lm));%sum for some dimension
Index: /CLIO/Analysis/Data19July/KK/PositionDirection.m
===================================================================
--- /CLIO/Analysis/Data19July/KK/PositionDirection.m	(revision 723)
+++ /CLIO/Analysis/Data19July/KK/PositionDirection.m	(revision 723)
@@ -0,0 +1,172 @@
+function [ChiV,Corr,ExterProfile]=PositionDirection(prof2,prof)
+[ChiV1,Corr1,ExterProfile1]=Position(prof2,prof);
+[ChiV2,Corr2,ExterProfile2]=Position(fliplr(prof2),prof);
+if ChiV1<ChiV2
+    ChiV=ChiV1;
+    Corr=Corr1;
+    ExterProfile=ExterProfile1;
+else
+    ChiV=ChiV2;
+    Corr=Corr2;
+    ExterProfile=ExterProfile2;
+end
+end
+
+
+function [ChiV,Corr,ExterProfile]=Position(prof2,prof)
+c1=find(cumsum(prof2)>=(sum(prof2)/2),1);
+c2=find(cumsum(prof)>=(sum(prof)/2),1);
+fi=(1+sqrt(5))/2;
+  T=100;
+a=c2-T; b=c2+T;
+while(abs(b-a)>1)
+x1=b-(b-a)/fi;
+x2=a+(b-a)/fi;
+y1=Chi(circshift(prof2,[ 1 round(x1)-c1]),prof);
+y2=Chi(circshift(prof2,[ 1 round(x2)-c1]),prof);
+if(y1<=y2)
+b=x2; 
+else
+a=x1; 
+end
+end
+y1=Chi(circshift(prof2,[ 1 round(a)-c1]),prof);
+y2=Chi(circshift(prof2,[ 1 round(b)-c1]),prof);
+if ((y1-y2)>0)  rez=b; ChiV=y2;
+else rez=a;ChiV=y1
+end
+ExterProfile=circshift(prof2,[ 1 round(rez)-c1]);
+Corr=sum(prof.*ExterProfile)/(sqrt(sum((prof).^2))*sqrt(sum((ExterProfile).^2)));
+end
+
+function Chi=Chi(funk1,funk2)
+Chi=sum((funk1-funk2).^2./abs(funk1+funk2+1e-10)./(length(funk1)));
+end
+% %one side
+% c1=centerOFweight(abs(funk1));%find center of weight of function
+% c2=centerOFweight(funk2b);%base funk
+% c1=c1-(T+1);
+% for i=1:(2*T+1)
+% c1=c1+1;
+%     delta=c2-c1;
+%     N1=zeros(size(funk1));
+% if delta>0
+%     N1(delta:end)=funk1(1:length(funk1)-delta+1);   
+% end
+% if c1==c2
+%     N1=funk1;
+% end
+% if delta<0
+%     N1(1:end+delta+1)=funk1(-delta:end);   
+%  end
+% 
+
+% [chi1(i),cor1(i)]=ChiCor(N1,funk2b);
+% end
+
+% [chi1V,chi1P]=min(chi1);
+% [cor1V,cor1P]=max(cor1);
+% 
+% %Other side
+% funk1fliped=flip(funk1);
+% c1=centerOFweight(abs(funk1fliped));%find center of weight of function
+% c2=centerOFweight(funk2b);
+% c1=c1-(T+1);
+% for i=1:(2*T+1)
+% c1=c1+1;
+%     delta=c2-c1;
+%     N1=zeros(size(funk1));
+% if delta>0
+%     N1(delta:length(N1))=funk1fliped(1:length(funk1)-delta+1);   
+% end
+% if c1==c2
+%     N1=funk1;
+% end
+% if delta<0
+%     N1(1:length(N1)+delta+1)=funk1fliped(-delta:length(funk1));   
+%  end
+% 
+% [chi2(i),cor2(i)]=ChiCor(N1,funk2b);
+% end
+% [chi2V,chi2P]=min(chi2);
+% [cor2V,cor2P]=max(cor2);
+% 
+% %compare both sides
+% 
+% 
+% %by minimum of chi2
+% 
+% if chi1V<chi2V
+%     c1=centerOFweight(abs(funk1))-(T+1)+chi1P;
+%     c2=centerOFweight(funk2b);
+%      delta=c2-c1;
+%     N1=zeros(size(funk1));
+% if delta>0
+%     N1(delta:length(N1))=funk1(1:length(funk1)-delta+1);   
+% end
+% if c1==c2
+%     N1=funk1;
+% end
+% if delta<0
+%     N1(1:length(N1)+delta+1)=funk1(-delta:length(funk1));   
+% end
+% ChiVal1=chi1V;
+% CorVal1=cor1V;
+% else
+%     
+%     c1=centerOFweight(abs(funk1fliped))+chi2P-(T+1);
+%     c2=centerOFweight(funk2b);
+%         delta=c2-c1;
+%     N1=zeros(size(funk1));
+% if delta>0
+%     N1(delta:length(N1))=funk1fliped(1:length(funk1)-delta+1);   
+% end
+% if c1==c2
+%     N1=funk1;
+% end
+% if delta<0
+%     N1(1:length(N1)+delta+1)=funk1fliped(-delta:length(funk1));   
+% end
+%   ChiVal1=chi2V;
+%   CorVal1=cor2V;
+% end
+% 
+% % %by maximum of corr
+% % 
+% % if cor1V>cor2V
+% %     c1=centerOFweight(abs(funk1))-201+cor1P;
+% %     c2=centerOFweight(funk2b);
+% %      delta=c2-c1;
+% %     N1c=zeros(size(funk1));
+% % if delta>0
+% %     N1c(delta:length(N1c))=funk1(1:length(funk1)-delta+1);   
+% % end
+% % if c1==c2
+% %     N1c=funk1;
+% % end
+% % if delta<0
+% %     N1c(1:length(N1c)+delta+1)=funk1(-delta:length(funk1));   
+% % end
+% % ChiVal2=chi1V;
+% % CorVal2=cor1V;
+% % else
+% %     
+% %     c1=centerOFweight(abs(funk1fliped))+cor2P-201;
+% %     c2=centerOFweight(funk2b);
+% %         delta=c2-c1;
+% %     N1c=zeros(size(funk1));
+% % if delta>0
+% %     N1c(delta:length(N1c))=funk1fliped(1:length(funk1)-delta+1);   
+% % end
+% % if c1==c2
+% %     N1c=funk1;
+% % end
+% % if delta<0
+% %     N1c(1:length(N1c)+delta+1)=funk1fliped(-delta:length(funk1));   
+% % end
+% %   ChiVal2=chi2V;
+% %   CorVal2=cor2V;
+% % end
+% 
+% ExterProfile=N1;
+% %ExterProfileCor=N1c;
Index: /CLIO/Analysis/Data19July/KK/SRf.m
===================================================================
--- /CLIO/Analysis/Data19July/KK/SRf.m	(revision 723)
+++ /CLIO/Analysis/Data19July/KK/SRf.m	(revision 723)
@@ -0,0 +1,43 @@
+function [RSP,ISP,alp]=SRf(SP,SS)
+%RECOVERY PART
+w0=SP(1);
+%Make LF extrapolation
+DSP=diff(SS)./diff(SP);
+rs=sum(DSP(1:3));%must be negative
+r=SS(1);
+alp=rs/r/w0-log(r)/(w0*w0);
+bet=2*log(r)/w0-rs/r;
+x=0:(w0-1);
+% syms xx
+% w0
+% rs
+% vpasolve(-2*w0*xx*exp(-xx*w0^2)==rs, xx)
+%  LFpart=exp(alp.*x.*x+bet.*x);
+LFpart=exp(-abs(log(r)/(w0^2)).*x.^2);
+ LFpart( LFpart>1)= 1; 
+%  LFpart=ones(size(x));
+% LFpart=exp(-1.8262e-05*x.^2);
+%  LFpart=1;
+%Make interpolation
+% SS=SS.*exp(-1.8262e-05*w0.^2)
+ISP = pchip(SP,SS,w0:max(SP));
+%Make HF extrapolation
+w=1:(2^12-length(ISP)-length(LFpart));
+wmax=max(SP);
+rosh=DSP(end);
+
+  alp=(rosh/ISP(end))*wmax;
+        c=0;
+        while (alp>0)
+            c=c+1;
+          rosh=  (SS(end)-SS(end-c))/(SP(end)-SP(end-c));
+          alp=(rosh/ISP(end))*wmax;
+        end
+        if (alp>-2) alp=-2;end
+HFpart=ISP(end)* (1+w/wmax).^(alp);
+HFpart(HFpart<eps)=eps;
+
+%Connect all together
+ISP=[LFpart,ISP,HFpart];
+%Mirror image of the spectrum
+RSP=[ISP,fliplr(ISP)];
Index: /CLIO/Analysis/Data19July/ProfileX.m
===================================================================
--- /CLIO/Analysis/Data19July/ProfileX.m	(revision 723)
+++ /CLIO/Analysis/Data19July/ProfileX.m	(revision 723)
@@ -0,0 +1,154 @@
+close all; clear all; clc
+load FFX.mat
+cd('./KK')
+
+for c=[1:7]%, 14:18
+    c
+% for c=1:20
+% SP=flip(Freq([3:4, 6:end-1]));
+% SS(c,:)=flip(FF(c,[3:4, 6:end-1]))/9/5/3;
+ SP=flip(Freq([3:6 8:end-1]));
+ SS(c,:)=flip(FF(c,[3:6 8:end-1]))/9/5/1.5;
+% SSM(c)=max(SS);
+[RSP,ISP(c,:)]=SRf(SP,SS(c,:));
+% 
+
+% figure (1)
+% plot(SP,SS(c,:),'or')
+% hold on
+%  plot(RSP)
+
+ PhaseH=-(imag(hilbert(log(abs(RSP)))));
+RPH=real((ifft(RSP.*exp(1i.*(PhaseH)))));
+
+
+% figure(2)
+[mx,ps]=max(RPH)
+ RPH= [zeros(1,230-ps), RPH(1:end-(230-ps))]
+%  plot((1:length(RPH))*1e-9/length(RPH),RPH)
+
+Gauss=@(mu,sig,A,x)A.*exp(-(x-mu).^2/(2*sig^2))
+% hold on
+%2.1e-12
+Profile5=Gauss(ps*1e-9/length(RPH),2.1e-12,mx,(1:length(RPH))*1e-9/length(RPH));
+% plot((1:length(RPH))*1e-9/length(RPH),Profile5,'r')
+ FWHM(c)=sum(RPH>max(RPH)/2)*1e-9/length(RPH)
+  FW01M(c)=sum(RPH>max(RPH)*0.1)*1e-9/length(RPH)
+ FW09M(c)=sum(RPH>max(RPH)*0.9)*1e-9/length(RPH)
+
+ 
+PR(c,:)=RPH;
+sum(Profile5>max(Profile5)/2)*1e-9/length(Profile5)
+% figure(1)
+% plot(abs(fft(RPH)),'r')
+% plot(abs(fft(Profile5./sum(Profile5))),'k')
+
+ end
+Bphase=sort([5.84 5.74 5.64 5.54 5.44 5.94 6.04]);
+
+% ds=1:20;
+ figure
+ hold on
+plot(Bphase,FW01M(:)*1e+12,'->k','linewidth',2)
+plot(Bphase,FWHM(:)*1e+12,'->b','linewidth',2)
+plot(Bphase,FW09M(:)*1e+12,'->r','linewidth',2)
+
+ set(gca,'fontsize',16)
+  xlabel('\phi_B, [a.u.]')
+ ylabel('Pulsewidth, [ps]')
+ legend('FW0.1M','FWHM','FW0.9M','Orientation','horizontal')
+ grid on
+ title('Pulsewidth vs. Buncher phase')
+ 
+figure
+[hc hc]=contourf((1:length(RPH))*1e-9/length(RPH)*1e+12,Bphase,PR*500)
+colorbar
+  set(hc,'LineStyle','none');
+xlim([10 50])
+ set(gca,'fontsize',16)
+
+ xlabel('Time, [ps]')
+ ylabel('\phi_B, [a.u.]')
+ title('Bunch profile')
+ 
+ 
+ figure
+ 
+ hold on
+plot(SP,  SS(2,:),'ob','linewidth',2)
+plot(0:(length(ISP)-1),  ISP(2,:),'-.b','linewidth',2)
+plot(SP,  SS(4,:),'om','linewidth',2)
+plot(0:(length(ISP)-1),  ISP(4,:),'-.m','linewidth',2)
+plot(SP,  SS(6,:),'ok','linewidth',2)
+plot(0:(length(ISP)-1),  ISP(6,:),'-.k','linewidth',2)
+plot(SP,  SS(1,:),'og','linewidth',2)
+plot(0:(length(ISP)-1),  ISP(1,:),'-.g','linewidth',2)
+plot(SP,  SS(3,:),'or','linewidth',2)
+plot(0:(length(ISP)-1),  ISP(3,:),'-.r','linewidth',2)
+
+ set(gca,'fontsize',16)
+ xlabel('Frequency, [GHz]')
+ legend('\phi_B=5.84','\phi_B=5.74','\phi_B=5.54','\phi_B=5.44','\phi_B=6.04')
+ % 2 4 6 8 10
+ title('Form Factor (BGS=10 mm)')
+ylabel('Amplitude, [a.u.]')
+xlim([0, 300])
+grid on
+ 
+  figure
+ 
+ hold on
+plot((0:(length(RPH)-1))*1e-9/length(RPH)*1e+12,  PR(1,:)*500,'-.b','linewidth',2)
+plot((0:(length(RPH)-1))*1e-9/length(RPH)*1e+12,  PR(3,:)*500,'-.m','linewidth',2)
+plot((0:(length(RPH)-1))*1e-9/length(RPH)*1e+12,  PR(5,:)*500,'-.k','linewidth',2)
+plot((0:(length(RPH)-1))*1e-9/length(RPH)*1e+12,  PR(2,:)*500,'-.g','linewidth',2)
+plot((0:(length(RPH)-1))*1e-9/length(RPH)*1e+12,  PR(7,:)*500,'-.r','linewidth',2)
+
+ set(gca,'fontsize',16)
+ xlabel('Time, [ps]')
+ legend('\phi_B=5.84','\phi_B=5.64','\phi_B=5.44','\phi_B=5.74','\phi_B=6.04')
+ % 2 4 6 8 10
+ title('Reconstructed profile')
+ylabel('Amplitude, [a.u.]')
+xlim([0, 50])
+grid on
+ 
+%   
+%   figure
+%  
+%  hold on
+% plot((0:(length(RPH)-1))*1e-9/length(RPH)*1e+12,  PR(5,:)*500,'-.b','linewidth',2)
+% plot((0:(length(RPH)-1))*1e-9/length(RPH)*1e+12,  PR(6,:)*500,'-.m','linewidth',2)
+% plot((0:(length(RPH)-1))*1e-9/length(RPH)*1e+12,  PR(7,:)*500,'-.k','linewidth',2)
+% plot((0:(length(RPH)-1))*1e-9/length(RPH)*1e+12,  PR(2,:)*500,'-.g','linewidth',2)
+% plot((0:(length(RPH)-1))*1e-9/length(RPH)*1e+12,  PR(1,:)*500,'-.r','linewidth',2)
+% 
+%  set(gca,'fontsize',16)
+%  xlabel('Time, [ps]')
+%  legend('\phi_B=8.14','\phi_B=8.19','\phi_B=8.24','\phi_B=5.54','\phi_B=5.44')
+%  % 2 4 6 8 10
+%  title('Reconstructed profile')
+% ylabel('Amplitude, [a.u.]')
+% xlim([0, 50])
+% grid on
+%  
+%  Bphase=sort([5.84 5.74 5.64 5.54 5.44 5.94 6.04]);
+% 
+%  time =(0:(length(RPH)-1))*1e-9/length(RPH)*1e+9;%ns
+%  Profile=PR(6,:)*500;
+%  
+% %  save P6.mat time Profile
+ 
+ 
+ 
+ 
+ 
+ 
+ 
+ 
+ 
+ 
+ 
+ 
+ 
+ cd('..')
Index: /CLIO/Analysis/Data19July/readsig.m
===================================================================
--- /CLIO/Analysis/Data19July/readsig.m	(revision 723)
+++ /CLIO/Analysis/Data19July/readsig.m	(revision 723)
@@ -0,0 +1,25 @@
+function [key_name, key_value, data]=readsig(file)
+fid = fopen(file);
+tline = fgets(fid);
+c=1;
+key_name=[];
+key_value=[];
+while ischar(tline)
+    
+    if (tline(1)=='#')
+    %disp(tline)
+    tline=regexprep(tline,'-','');
+    tline=regexprep(tline,'_','');
+    C = strsplit(tline(2:end),'=');
+    key_name=[key_name,num2str(C{1,1})];
+    key_value=[key_value;str2double(C{1,2})];
+    else 
+    if c<13
+    data(c,:)=sscanf(tline,'%f\t%f\t%f\n',[3 Inf]);
+    c=c+1;
+    end
+    end
+    tline = fgets(fid);
+end
+
+fclose(fid);
