source: ETALON/simulations/Smith_Purcell_infinite_width.m @ 1

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%  Computation of the SPM radiation using the surface current model
%           M. Labat - Decembre 2011

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Input parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%% Electron beam
Ee=100;                         % [MeV]
Qe=0.6*1e-9;                      % [C] Charge   
x_0=2e-3;                       % [m] Bunch height
sigma_x=0.5e-3;                   % [m] Bucnh size
sigma_y=1.1e-3;                   % [m] Bucnh size
sigma_t=3.3/2.35*1e-12;          % [s] Bunch length --> SOLEIL
%sigma_t=2/2.35*1e-15;           % [s] Bunch length --> LOA
prof=1;                         % [1: gaussian / 2: other] Bunch profile   

 %%% Gratings
N_gratings=1;               % Number of gratings
l=zeros(1,N_gratings);      % [m] Grating period
Ng1=zeros(1,N_gratings);    % [-] Number of periods
L1=zeros(1,N_gratings);     % [m] Grating length
h1=zeros(1,N_gratings);     % [m] Tooth depth
d1=zeros(1,N_gratings);     % [m] Tooth width
F1=zeros(1,N_gratings);     % [-] Number of facets
R1=zeros(1,N_gratings);

    %%% Grating 1 parameters
l(1)=1*1e-3;                % [m] Grating period --> SOLEIL
%l(1)=1*1e-6;                % [m] Grating period --> LOA
Ng(1)=40;                 % [-] Number of periods
alpha1(1)=30*pi/180;        % [rad] Blaze angle
    %%% Grating 2 parameters
%l(2)=2*1e-3;                % [m] Grating period --> SOLEIL
l(2)=5*1e-6;                % [m] Grating period --> LOA
Ng(2)=2000;                 % [-] Number of periods
alpha1(2)=20*pi/180;        % [rad] Blaze angle
    %%% Grating 3 parameters
l(3)=1*1e-3;               % [m] Grating period --> SOLEIL
%l(3)=5*1e-6;                % [m] Grating period --> LOA
Ng(3)=40;                 % [-] Number of periods
alpha1(3)=30*pi/180;        % [rad] Blaze angle

 %%% Detection
theta_min=1*pi/9;           % [rad]
theta_max=pi;               % [rad]
phi=0;                      % [rad]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Variable definition %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%% Constants
c=3e8;                  % [m/s]
q=1.6e-19;              % [C]
E0=0.511;               % [MeV] Electron rest energy
Ne=Qe/q                             % [-] Nb of electrons
gamma=Ee/E0;                        % [-]
beta=sqrt(1-1/power(gamma,2));      % [-]
v=beta*c;                           % [m]

 %%% Time coordinates
dt=sigma_t/20;          % Time sample
Dt=sigma_t*100;         % Time window
Nt=Dt/dt;
t_0=Dt/2;                           % [s] Beam position in time
for i=1:Nt
    t(i)=i*dt;
end

if (prof==1)
    rho = exp(-0.5*power((t-t_0)/sigma_t,2));     % gaussian profile
else
    rho=t*0;
end;

 %%% Electron beam Fourier transform
m = length(rho);                % Window length
n = pow2(nextpow2(m));          % Transform length
TF_rho_tot = fft(rho,n);            % DFT
for i=1:round(n/2)
    nu(i)=i/(n*dt);
    TF_rho(i)=TF_rho_tot(i);
    nu_tot(i)=i/(n*dt);
end
for i=(round(n/2)+1):n
    nu_tot(i)=(i-n)/(n*dt);
end
lambda_V=c./nu;
Np=length(lambda_V);
lambda_Vtot=c./nu_tot;
TF_rho_norm=abs(TF_rho);
TF_rho_norm=TF_rho_norm/max(TF_rho_norm);



for j=1:N_gratings
for ii=1:Np   
    
%%% Radiation parameters
nh=1;                   % [-] Order of the radiation
lambda=lambda_V(ii);    % [m]
theta=acos(1/beta-nh*lambda/l(j))
theta_V(j,ii)=acos(1/beta-nh*lambda/l(j));
omega=2*pi*c/lambda;    % [1/s]

%%% Calculation of lambda_e
lambda_e=lambda/(2*pi)*1/(sqrt(1+power(gamma*beta*sin(theta)*sin(phi),2))/(gamma*beta*c));

%%% Calculation of the grating parameters
L=Ng(j)*l(j);                   % [m] Grating length
l1=l(j)*power(cos(alpha1(j)),1);         % [m] Length of facet 1
    % ? % power(cos(alpha1(j)),2);
alpha2=pi/2+alpha1(j);          % [rad]
l2=l(j)-l1;                     % [m] Length of facet 2
h=abs(l1*tan(alpha1(j)));       % [m] Tooth heigth
x=h/lambda_e;

%%% 
kx=2*pi*sin(theta)*cos(phi)/lambda;
ky=2*pi*sin(theta)*sin(phi)/lambda;

tau1=2*pi*nh*cos(alpha1)-kx*l(j)*sin(alpha1(j));
tau2=2*pi*nh*cos(alpha2)-kx*l(j)*sin(alpha2);

eps1=2*pi*nh*l1/l(j)-kx*h;
eps2=2*pi*nh*l2/l(j);

a2t=(1/lambda_e*sin(alpha2)*(cos(eps2)*exp(-x)-cos(kx*h))+tau2*(sin(eps2)*exp(-x)+sin(kx*h)))/(power(1/lambda_e*sin(alpha2),2)+power(tau2,2));
b2t=(1/lambda_e*sin(alpha2)*(sin(eps2)*exp(-x)+sin(kx*h))-tau2*(cos(eps2)*exp(-x)-cos(kx*h)))/(power(1/lambda_e*sin(alpha2),2)+power(tau2,2));
    
a1=(1/lambda_e*sin(alpha1(j))*(cos(eps1)-exp(-x))+tau1*sin(eps1))/(power(1/lambda_e*sin(alpha1(j)),2)+power(tau1,2));
b1=(1/lambda_e*sin(alpha1(j))*sin(eps1)-tau1*(cos(eps1)-exp(-x)))/(power(1/lambda_e*sin(alpha1(j)),2)+power(tau1,2));

a2=a2t*cos(2*pi*l1/l(j))-b2t*sin(2*pi*l1/l(j));
b2=b2t*cos(2*pi*l1/l(j))+a2t*sin(2*pi*l1/l(j));

as=a1*sin(alpha1(j))+a2*sin(alpha2);
bs=b1*sin(alpha1(j))+b2*sin(alpha2);
ac=a1*cos(alpha1(j))+a2*cos(alpha2);
bc=b1*cos(alpha1(j))+b2*cos(alpha2);
a111=(power(as,2)+power(bs,2))*power(cos(theta)*cos(phi),2);
a112=(power(as,2)+power(bs,2))*power(ky*lambda_e,2)*power(cos(theta)*cos(phi),2);
a113=(power(ac,2)+power(bc,2))*power(sin(theta),2);
a114=-2*(as*ac+bs*bc)*sin(theta)*cos(theta)*cos(phi);
a115=2*(as*bc-bs*ac)*ky*lambda_e*sin(theta)*cos(theta)*sin(phi);

a211=power(as,2)+power(bs,2);
a212=power(sin(phi),2)*(1+power(kx*lambda_e,2));

aa1t=a111+a112+a113+a114+a115;
aa2t=a211*a212;

R2=aa1t+aa2t;
R2_V(ii)=R2;
%R2=0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%% Calculation of W1
W1=2*pi*power(q,2)*power(nh,2)*power(beta,3)/(power(l(j),2)*power(1-beta*cos(theta),3))*exp(-2*x_0/lambda_e)*L*R2;
    % ? % *1e-7;
    % ? % *power(beta,2);
    % ? % *power(q*Ne,2);

%%% Calculation of S_inc
S_inc=Ne*0.5*exp(2*power(sigma_x/lambda_e,2)-2*x_0/lambda_e)*erfc(-x_0/(sqrt(2)*sigma_x)+sqrt(2)*sigma_x/lambda_e);

%%% Calculation of S_coh
S_coh_X=power(0.5*exp(0.5*power(sigma_x/lambda_e,2)-x_0/lambda_e)*erfc(-x_0/(sqrt(2)*sigma_x)+sigma_x/(sqrt(2)*lambda_e)),2);
S_coh_Y=power(exp(-0.5*power(ky*sigma_y,2)),2);
S_coh_Z=power(TF_rho_norm,2);
S_coh=power(Ne,2)*S_coh_X*S_coh_Y*S_coh_Z;

%%% Calculation of W
if (theta>theta_min)&&(theta<theta_max)
    W_inc(j,ii)=W1*S_inc/(L*0.01)                 % [J/sr/cm]
    W(j,ii)=W1*(S_inc+S_coh(ii))/(L*0.01)   % [J/sr/cm]
else
    W_inc(j,ii)=0;
    W(j,ii)=0;
end

end

end

figure(1)
subplot(2,1,1)
plot(t*1e12,rho)
xlabel('Time (s)')
ylabel('Power')
title('{\bf Time sample}')
grid on

subplot(2,1,2)
semilogx(lambda_V,TF_rho_norm)
xlabel('Wavelength [m]')
ylabel('TF_rho')
title('{\bf Bunch profile Fourier Transform}')
grid on
%xlim([1e-6 1e-3])


figure(2)
subplot(2,1,1)
loglog(lambda_V,W(1,:),'*');%,lambda_V,W(2,:),'*',lambda_V,W(3,:),'*');
hold on
loglog(lambda_V,W_inc(1,:));%,lambda_V,W_inc(2,:),lambda_V,W_inc(3,:));
%xlim([200*1e-6 40000*1e-6])
xlabel('Wavelength [m]')
ylabel('W [J/sr/cm]')
%ylim ([1e-10 1e-0])
legend('g1','g2','g3',3)
grid on
hold off

subplot(2,1,2)
plot(theta_V(1,:)*180/pi,W(1,:),'*');%,theta_V(2,:),W(2,:),theta_V(3,:),W(3,:));
%xlim([0 3.15])
xlabel('Angle [degrees]')
ylabel('W [J/sr]')
grid on
%ylim ([1e-20 1e-3])

figure(3)
plot(theta_V(1,:)*180/pi,R2_V(:),'*');
xlabel('Angle [degres]')
ylabel('R2')
grid on