source: ETALON/CLIO/CTR+gfw/GFWnew/GFWnew.m @ 383

 function  [nu_Hz,SP_spec]=GFWnew(var)
% 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 whole bunch)
% var='PrSelect=1; P_file=''../p20.mat''; PlotFunc=[1 1 1 1 1]; ';
 %% 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=8;%Grating period, mm
alpha1=30;%Blaze angle, degrees
gam=100;%relativistic gamma
Qbunch_nC=1.2e-9;%Bunch charge, nC (if T=1,than chagre at experiment position, if not, than charge at the exit of gun)
x0_mm=6; %Position of beam centroid above grating, mm
fwhm_x_mm=2; %FWHM_X  of Gaussian beam at the centre of grating, mm
fwhm_y_mm=2; %FWHM_Y  of Gaussian beam at the centre of grating, mm
pulse_len_ps=2.5;%Nominal bunch duration (FWHM), ps 
BAd_mm=300;% 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
% ThetaStart=pi/2-atan(Aed/2/BAd);
% ThetaStop=pi/2+atan(Aed/2/BAd);
ThetaStart=alpha1/180*pi;%initial theta angle
ThetaStop=pi-pi/10;%final theta angle
ThetaDiv=200;%steps in theta
phiDiv=10;%steps in phi

if (~isempty(var))
disp(var);
eval(var);
end %if

%save all parameters
name=regexprep(var,'.mat','');
name=regexprep(name,'../','');
name=regexprep(name,';','');
save([ name '.mat'])
%% Profile selection
if PrSelect==1
    FFTspec=1;
  if  exist(P_file,'file')==2;
    load([P_file]);
  else
      disp('Set correct file name!')
  end
else
    FFTspec=1;%how to calculate form factor: 1 -- FFT; 2 -- only fro assymetric gaussian
%Gaussian
epsl=0.5;                                                           %Asymmetry factor
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+13).*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/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
lol=x0_mm-5*sigma_x;                                       
if (lol<hi) lol=1.001*hi; end
upl=x0_mm+5*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

%% 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)*0.5*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==10)%
       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=(1:1:length(SpectrumFFT))./(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)
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!!!
     disp('fail to find Form Factro 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
%% Interration 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)));        
            
    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(i)/180*pi));           % variable effective grating length, determined basically by the cone entrance diameter of XXXmm*/

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)
end
nu_Hz=c./(lam*1e-3);%Hz
SP_spec=I1*0.1*lgr_e_mm*(el_mm/m*1e-3)*2.*(lam/nu_Hz)*1e+18;%muJ/THz

%% Plots
if PlotFunc(2)==1
 figure(2)
 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)
plot(thet,sum_tr,'r','linew',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)
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)
semilogy(nu_Hz*1e-12,SP_spec,'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')
% legend('SP','CTR')
 grid on
xlim([min(nu_Hz) max(nu_Hz)].*1e-12)
end