1 | function [nu_Hz,SP_spec_1,SP_spec_N,lam,E_total,sum_tr,thet]=GFWnew(var)
|
---|
2 |
|
---|
3 | % single_elec
|
---|
4 | % close all; clear all; clc
|
---|
5 | %% Configuration of calculation
|
---|
6 | %PrSelect == what profile to use: 1 --from file; 2 - assymmetric gaussian
|
---|
7 | %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)
|
---|
8 | %PlotFunc:
|
---|
9 | % 1 -- Form factor
|
---|
10 | % 2 -- Single electron yield from Theta
|
---|
11 | % 3 -- dI_N/dOmega (SP spetrum for whole bunch)
|
---|
12 | % 4 -- Bunch profile
|
---|
13 | % 5 -- dI(nu)/dnu (SP spetrum for signle elctron bunch)
|
---|
14 | % var='PrSelect=1; P_file=''../p20.mat''; PlotFunc=[1 1 1 1 1]; ';
|
---|
15 | PrSelect=2;
|
---|
16 | PlotFunc=[0 0 0 0 0];
|
---|
17 | P_file='tmp.mat';
|
---|
18 | %% Physical constants
|
---|
19 | c=299792458;%m/s speed of ligth
|
---|
20 | qe = 1.602176487e-19*c*10;%Electron charge, statColumns
|
---|
21 | %% Parameters block (by default)
|
---|
22 | T=1;%Transmissivity (from gun exit to experimental position)
|
---|
23 | m=1;%Order of radiation
|
---|
24 | lgr_mm=180;%Grating length, mm
|
---|
25 | wid_mm=40;%Grating width, mm
|
---|
26 | el_mm=8;%Grating period, mm
|
---|
27 | alpha1=30;%Blaze angle, degrees
|
---|
28 | gam=100;%relativistic gamma
|
---|
29 | Qbunch_nC=1.2;%Bunch charge, nC (if T=1,than chagre at experiment position, if not, than charge at the exit of gun)
|
---|
30 | x0_mm=6; %Position of beam centroid above grating, mm
|
---|
31 | fwhm_x_mm=2; %FWHM_X of Gaussian beam at the centre of grating, mm
|
---|
32 | fwhm_y_mm=2; %FWHM_Y of Gaussian beam at the centre of grating, mm
|
---|
33 | pulse_len_ps=2.5;%Nominal bunch duration (FWHM), ps
|
---|
34 | epsl=1;%Asymmetry factor for assymmetric Gaussian
|
---|
35 | BAd_mm=300;% beam - aperture distance , mm
|
---|
36 | Aed_mm=50;%Aperture entrance diameter, mm
|
---|
37 | %All angles are in radians
|
---|
38 | phiStop=atan(Aed_mm/2/BAd_mm);% phi go from 0 to phiStop
|
---|
39 | % phiStop=7/180*pi;
|
---|
40 | % ThetaStart=pi/2-atan(Aed_mm/2/BAd_mm);
|
---|
41 | % ThetaStop=pi/2+atan(Aed_mm/2/BAd_mm);
|
---|
42 | ThetaStart=alpha1/180*pi;%initial theta angle
|
---|
43 | ThetaStop=170/180*pi;%final theta angle
|
---|
44 | % ThetaStart=80/180*pi;%initial theta angle
|
---|
45 | % ThetaStop=140/180*pi;%final theta angle
|
---|
46 | ThetaDiv=10;%steps in theta
|
---|
47 | phiDiv=3;%steps in phi
|
---|
48 |
|
---|
49 | if (~isempty(var))
|
---|
50 | strread(var,'%s','delimiter',' ')
|
---|
51 | eval(var);
|
---|
52 | end %if
|
---|
53 |
|
---|
54 | %save all parameters
|
---|
55 | name=regexprep(var,'.mat','');
|
---|
56 | name=regexprep(name,'[../''[]]','');
|
---|
57 | name=regexprep(name,'[; ]','_');
|
---|
58 | % save([ name 'SP.mat'])
|
---|
59 | %% Profile selection
|
---|
60 | if PrSelect==1
|
---|
61 | FFTspec=1;
|
---|
62 | if exist(P_file,'file')==2;
|
---|
63 | load([P_file]);
|
---|
64 | else
|
---|
65 | error('Set correct file name!')
|
---|
66 | end
|
---|
67 | else
|
---|
68 | FFTspec=1;%how to calculate form factor: 1 -- FFT; 2 -- only fro assymetric gaussian
|
---|
69 | %Gaussian
|
---|
70 | sigma_t= pulse_len_ps/(sqrt(2*log(2))*(1+epsl))*1e-3; %The sigma of the leading edge ns
|
---|
71 |
|
---|
72 | nu=(1:1e+10:3e+16).*1e-9;%Ghz
|
---|
73 |
|
---|
74 | sigma_L= sigma_t*epsl;
|
---|
75 | sigma_H= sigma_t;
|
---|
76 | time=(1:(length(nu)*2))*(1/(length(nu)*2*(nu(2)-nu(1))));
|
---|
77 | xc=mean(time);
|
---|
78 | G1=exp(-(time-xc).^2./(2*sigma_L^2));
|
---|
79 | G2=exp(-(time-xc).^2./(2*sigma_H^2));
|
---|
80 | G1(time>xc)=0;
|
---|
81 | G2(time<xc)=0;
|
---|
82 | Profile=(G1+G2)./sum(G1+G2);
|
---|
83 | end
|
---|
84 |
|
---|
85 | %% Calculation block for additional parameters
|
---|
86 | n_e_t=round(Qbunch_nC*1e-9/1.6e-19*T); %Total number of electrons in the bunch
|
---|
87 | % n_e_t=1e+11;
|
---|
88 | SolAng=pi*((Aed_mm/2)^2)/(BAd_mm^2); %solid angle from beam to aperture
|
---|
89 | % Grating
|
---|
90 | alpha1=alpha1/180*pi;
|
---|
91 | el1=el_mm*(cos(alpha1)^2);
|
---|
92 | el2=el_mm-el1;
|
---|
93 | hi=abs(el1*tan(alpha1));
|
---|
94 | alpha2=-(pi/2-alpha1);
|
---|
95 | % Beam
|
---|
96 | beta=sqrt(gam^2 - 1)/gam;
|
---|
97 | sigma_x= fwhm_x_mm/(2*sqrt(2*log(2)));
|
---|
98 | sigma_y= fwhm_y_mm/(2*sqrt(2*log(2)));
|
---|
99 | % N_e above grating
|
---|
100 | N_sig=3;
|
---|
101 | lol=x0_mm-N_sig*sigma_x;
|
---|
102 | if (lol<hi) lol=1.001*hi; end
|
---|
103 | upl=x0_mm+N_sig*sigma_x;
|
---|
104 | GAU = @(x) exp(-(x-x0_mm).*(x-x0_mm)/(2*sigma_x*sigma_x));
|
---|
105 | res1 = integral(GAU,lol,upl)/(sigma_x*sqrt(2*pi));%part of the beam is above the grating
|
---|
106 | disp([num2str(res1) ' of the beam above the grating'])
|
---|
107 | n_e=n_e_t*res1; %Number of electrons above the grating
|
---|
108 | % save([ name '_additional_SP.mat'])
|
---|
109 | %% Estimation grating factor R^2, coherent & incoherent integrals
|
---|
110 | %%
|
---|
111 | %
|
---|
112 | % $$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$$
|
---|
113 | %
|
---|
114 | % In this block we calculate
|
---|
115 | %%
|
---|
116 | %
|
---|
117 | % $$| n \times n \times \sum G(...)|^2 $$
|
---|
118 | %
|
---|
119 | % where other terms (except L - grating length) are in f1, f2, f4. In general:
|
---|
120 | %%
|
---|
121 | %
|
---|
122 | % $$W_N=W_1 |_{x_0=0}(NS_{inc}+N^2S_{coh}) $$
|
---|
123 | %
|
---|
124 | % where s_coh and s_inc are:
|
---|
125 | %%
|
---|
126 | %
|
---|
127 | % $$S_{inc}=\int_h^{\infty}dxX(x)$$
|
---|
128 | %
|
---|
129 | % $$S_{coh}=|\int_h^{\infty}dxX(x)\tilde{Y}(ky)\tilde{T}(\omega)|^2 $$
|
---|
130 | %
|
---|
131 | % for now (in this block)
|
---|
132 | %%
|
---|
133 | %
|
---|
134 | % $$S^*_{inc}=\int_h^{\infty}dxX(x)*| n \times n \times \sum G(...)|^2\\ $$
|
---|
135 | %
|
---|
136 | % $$S^*_{coh}=\int_h^{\infty}dxX(x)*| n \times n \times \sum G(...)| $$
|
---|
137 | %
|
---|
138 | par = struct('Theta',NaN,'phi',NaN,'m',m,'alpha1',alpha1,'facet',NaN,'gamma',gam,'x0',x0_mm,'l',el_mm,'w',wid_mm);
|
---|
139 | for i=1:ThetaDiv%theta
|
---|
140 | thet(i)=(ThetaStart+(ThetaStop-ThetaStart)*(i-1)/(ThetaDiv-1))/pi*180;
|
---|
141 | theta=thet(i)*pi/180;
|
---|
142 | for ii=1:phiDiv;%phi
|
---|
143 | phi=(phiStop/(phiDiv-1))*(ii-1);
|
---|
144 |
|
---|
145 | for ind=1:21;%x0
|
---|
146 | xv=lol+ (ind-1)*(N_sig/10)*sigma_x;%height above grating
|
---|
147 |
|
---|
148 | par.phi=phi;
|
---|
149 | par.Theta=theta;
|
---|
150 | par.x0=xv;
|
---|
151 | par.facet=1;
|
---|
152 |
|
---|
153 | [res1,res2,res3,res4]=QUx(par,50,200);
|
---|
154 | G1x=(res1*tan(alpha1) +1i*res2*tan(alpha1))/(el_mm*pi);
|
---|
155 | G1z=(res1+1i*res2)/(el_mm*pi);
|
---|
156 | G1y=(res3*tan(alpha1)-1i*res4*tan(alpha1))/(el_mm*pi);
|
---|
157 |
|
---|
158 | par.facet=2;
|
---|
159 |
|
---|
160 | [res1,res2,res3,res4]=QUx(par,50,200);
|
---|
161 | G2x=(res1*tan(alpha2) +1i*res2*tan(alpha2))/(el_mm*pi);
|
---|
162 | G2z=(res1+1i*res2)/(el_mm*pi);
|
---|
163 | G2y=(res3*tan(alpha2)-1i*res4*tan(alpha2))/(el_mm*pi);
|
---|
164 |
|
---|
165 | Gx= G1x+ G2x;
|
---|
166 | Gy= G1y+G2y;
|
---|
167 | Gz= G1z+G2z;
|
---|
168 | G=[Gx Gy Gz];
|
---|
169 |
|
---|
170 | nv=complex([ sin(theta)*cos(phi) sin(theta)*sin(phi) cos(theta)]); %dir. unit vector
|
---|
171 | nvp=complex([cos(theta)*cos(phi) cos(theta)*sin(phi) -sin(theta)]);%1st pol. unit vector
|
---|
172 | nvv=complex([-sin(phi) cos(phi) 0]);%2nd pol. unit vector
|
---|
173 |
|
---|
174 | a1t=(G*nvp')*conj(G*nvp');% square of magn. in first polarisation
|
---|
175 | a2t=(G*nvv')*conj(G*nvv');% square of magn. in second polarisation
|
---|
176 |
|
---|
177 | f5c=a1t+a2t; % total
|
---|
178 |
|
---|
179 | if (ind==11)%
|
---|
180 |
|
---|
181 | single_elec(i,ii)= real(f5c);
|
---|
182 | single_elec_p1(i,ii)= real(a1t);
|
---|
183 | single_elec_p2(i,ii)= real(a2t) ;
|
---|
184 | end % /*store the single-electron R^2 */
|
---|
185 |
|
---|
186 | x0_bin(ind)= xv;
|
---|
187 | inc_bin(ind)= real(f5c)*exp(-(xv-x0_mm)*(xv-x0_mm)/(2*sigma_x*sigma_x));
|
---|
188 | inc_bin_p1(ind)= real(a1t)*exp(-(xv-x0_mm)*(xv-x0_mm)/(2*sigma_x*sigma_x));
|
---|
189 | inc_bin_p2(ind)= real(a2t)*exp(-(xv-x0_mm)*(xv-x0_mm)/(2*sigma_x*sigma_x));
|
---|
190 | coh_bin(ind)= sqrt(real(f5c))*exp(-(xv-x0_mm)*(xv-x0_mm)/(2*sigma_x*sigma_x));
|
---|
191 | coh_bin_p1(ind)= sqrt(real(a1t))*exp(-(xv-x0_mm)*(xv-x0_mm)/(2*sigma_x*sigma_x));
|
---|
192 | coh_bin_p2(ind)= sqrt(real(a2t))*exp(-(xv-x0_mm)*(xv-x0_mm)/(2*sigma_x*sigma_x));
|
---|
193 | end %of x0
|
---|
194 | % integral over hight from grating
|
---|
195 | toler=1000;
|
---|
196 | dist1=(x0_bin(end)-x0_bin(1))/toler;
|
---|
197 | dist=x0_bin(1):dist1:x0_bin(end);
|
---|
198 | incans=pchip(x0_bin,inc_bin,dist);
|
---|
199 | s_inc_x(i,ii)= trapz(dist,incans)/(sqrt(2*pi)*sigma_x);
|
---|
200 | cohans=pchip(x0_bin,coh_bin,dist);
|
---|
201 | s_coh_x(i,ii)= ( trapz(dist,cohans)/(sqrt(2*pi)*sigma_x))^2;
|
---|
202 |
|
---|
203 | incans=pchip(x0_bin,inc_bin_p1,dist);
|
---|
204 | s_inc_x_p1(i,ii)= trapz(dist,incans)/(sqrt(2*pi)*sigma_x);
|
---|
205 | cohans=pchip(x0_bin,coh_bin_p1,dist);
|
---|
206 | s_coh_x_p1(i,ii)= ( trapz(dist,cohans)/(sqrt(2*pi)*sigma_x))^2;
|
---|
207 |
|
---|
208 | incans=pchip(x0_bin,inc_bin_p2,dist);
|
---|
209 | s_inc_x_p2(i,ii)= trapz(dist,incans)/(sqrt(2*pi)*sigma_x);
|
---|
210 | cohans=pchip(x0_bin,coh_bin_p2,dist);
|
---|
211 | s_coh_x_p2(i,ii)= (trapz(dist,cohans)/(sqrt(2*pi)*sigma_x))^2;
|
---|
212 |
|
---|
213 | end%of phi
|
---|
214 | end%of theta
|
---|
215 | %% Form factor
|
---|
216 | % form factor is calculated numericaly by FFT
|
---|
217 | SpectrumFFT=abs(fft(Profile)).^2;%Fourier transform
|
---|
218 | SpectrumHalf=SpectrumFFT(1:round(length(SpectrumFFT)/2));%only half of FT spectrum hold usefull information
|
---|
219 | Frequency=(0:(length(SpectrumFFT)-1))./(time(2)-time(1))/length(SpectrumFFT)*1e-3;%THz
|
---|
220 | FreqHalf=Frequency(1:round(length(Frequency)/2));%half of frequency range
|
---|
221 | SpectrumHalf=SpectrumHalf./max(SpectrumHalf);%normalization of Spectrum
|
---|
222 | Nmax=find(FreqHalf>=60,1);%2.5 THz is more than enoght
|
---|
223 | % Interpolation is applyed to decrease error in calculation (freq step in FFT is to wide for this calculation )
|
---|
224 | FreqInter=FreqHalf(1):((FreqHalf(2)-FreqHalf(1))*1e-3):FreqHalf(Nmax);%interpolation diapason
|
---|
225 | SpecInter=pchip( FreqHalf(1:Nmax),SpectrumHalf(1:Nmax),FreqInter);%interpolation
|
---|
226 | if PlotFunc(1)==1
|
---|
227 | figure(1)
|
---|
228 | clf
|
---|
229 | plot(FreqInter,SpecInter,'-k','linew',2)
|
---|
230 | hold on
|
---|
231 | grid on
|
---|
232 | set(gca,'fontsize',16)
|
---|
233 | ylabel('F(\nu)')
|
---|
234 | xlabel('\nu, [THz]')
|
---|
235 | legend('Form factor')
|
---|
236 | end
|
---|
237 | %% I_N calculation
|
---|
238 | f1=2*pi*qe*qe*1e-7; % in Joules
|
---|
239 | f2=100/(el_mm^2);
|
---|
240 | f4=(m^2)*(beta./(1-beta*cos(thet*pi/180))).^3;
|
---|
241 | phi=(phiStop/(phiDiv-1))*((1:phiDiv)-1);
|
---|
242 | %lam_e=lam*beta*gam/(2*pi*sqrt(1+(beta*gam*sin(theta)*sin(phi0))^2));
|
---|
243 | lam= el_mm*(1/beta-cos(thet*pi/180))/m;%mm
|
---|
244 |
|
---|
245 | for i=1:ThetaDiv%theta
|
---|
246 | for ii=1:phiDiv;%phi
|
---|
247 | ky=(2*pi/lam(i))*sin(thet(i)*pi/180)*sin(phi(ii));
|
---|
248 | om=2*pi*0.3/lam(i); %c=0.3 mm/ps, hence om is in ps^-1
|
---|
249 | %calculate longitudinal coherence effects
|
---|
250 | %Form factor
|
---|
251 | if FFTspec==1
|
---|
252 | s_coh_long=SpecInter(find(2*pi*FreqInter>=om,1));
|
---|
253 | if isempty(s_coh_long)
|
---|
254 | s_coh_long=0;%ACHTUNG!!!
|
---|
255 | error('fail to find Form Factor value')
|
---|
256 | end
|
---|
257 | else
|
---|
258 | %this work with new versions of MATLAB (2015a for example)
|
---|
259 | u1= om*sigma_t;
|
---|
260 | u2= om*sigma_t*epsl;
|
---|
261 | f_z_c=(exp(-(u1*u1)/2)+ epsl*exp(-(u2*u2)/2))/(1+epsl);
|
---|
262 | f_z_s=-2/(sqrt(pi)*(1+epsl))*(dawson(u2/sqrt(2))*epsl-dawson(u1/sqrt(2)));
|
---|
263 | s_coh_long=f_z_c.^2+f_z_s.^2;
|
---|
264 | end
|
---|
265 | %%
|
---|
266 |
|
---|
267 | Y2=exp(-((ky*sigma_y)^2)); % Ysquare
|
---|
268 |
|
---|
269 | s_coh_xy= Y2*s_coh_x(i,ii);
|
---|
270 | s_coh_xy_p1= Y2*s_coh_x_p1(i,ii);
|
---|
271 | s_coh_xy_p2= Y2*s_coh_x_p2(i,ii);
|
---|
272 | FFact(i)=s_coh_long;
|
---|
273 | bin(i,1)=thet(i);
|
---|
274 | 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);
|
---|
275 | 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);
|
---|
276 | 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);
|
---|
277 | bin_single(i,1)= thet(i); % store for single electron yield
|
---|
278 | bin_single(i,ii+1)= f1*f2*f4(i)*single_elec(i,ii);
|
---|
279 | bin_single_p1(i,ii+1)= f1*f2*f4(i)*single_elec_p1(i,ii);
|
---|
280 | bin_single_p2(i,ii+1)= f1*f2*f4(i)*single_elec_p2(i,ii);
|
---|
281 | end% of ii (phi)
|
---|
282 | end % of i (theta)
|
---|
283 |
|
---|
284 | if PlotFunc(1)==1
|
---|
285 | figure(1)
|
---|
286 | hold on
|
---|
287 | plot(c./(lam*1e-3)*1e-12,FFact,'or','linew',2)
|
---|
288 | legend('Form factor','Used points')
|
---|
289 | end
|
---|
290 | %%
|
---|
291 | % At this point calculation of
|
---|
292 | %
|
---|
293 | % $$W_N=W_1 |_{x_0=0}(NS_{inc}+N^2S_{coh}) $$
|
---|
294 | %
|
---|
295 | % is complete. But there is dependence from angles. So next step is
|
---|
296 | % integration
|
---|
297 | %% Integration and averaging
|
---|
298 | philmt=phiDiv+1;%number of steps +1
|
---|
299 | %phiDiv + 2 -- average over phi
|
---|
300 | for i=1:ThetaDiv%theta
|
---|
301 | 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
|
---|
302 | bin_p1(i,phiDiv+2)=(bin_p1(i,2)+ 2*sum(bin_p1(i,3:philmt)))/(2*(philmt-1)-1); % ditto for pol. 1
|
---|
303 | bin_p2(i,phiDiv+2)=(bin_p2(i,2)+ 2*sum(bin_p2(i,3:philmt)))/(2*(philmt-1)-1); % and for pol. 2
|
---|
304 | bin_single(i,phiDiv+2)= (bin_single(i,2)+ 2*sum(bin_single(i,3:philmt)))/(2*(philmt-1)-1);
|
---|
305 | bin_single_p1(i,phiDiv+2)= (bin_single_p1(i,2)+ 2*sum(bin_single_p1(i,3:philmt)))/(2*(philmt-1)-1);
|
---|
306 | bin_single_p2(i,phiDiv+2)= (bin_single_p2(i,2)+ 2*sum(bin_single_p2(i,3:philmt)))/(2*(philmt-1)-1);
|
---|
307 | end
|
---|
308 |
|
---|
309 | %phiDiv + 3 -- integral
|
---|
310 | % integrate over phi
|
---|
311 | d_theta=(ThetaStop-ThetaStart)/(ThetaDiv-1); % d_theta expressed in rad
|
---|
312 | d_phi= (phiStop/(phiDiv-1)); %d_phi expressed in rad
|
---|
313 |
|
---|
314 | for i=1:ThetaDiv%theta
|
---|
315 |
|
---|
316 | theta=(ThetaStart+(ThetaStop-ThetaStart)*(i-1)/(ThetaDiv-1));
|
---|
317 |
|
---|
318 | sum0 = sin(theta)*d_theta*d_phi*sum(bin(i,3:(phiDiv+1)));
|
---|
319 | sum_p1 =sin(theta)*d_theta*d_phi*sum(bin_p1(i,3:(phiDiv+1))); %pol. 1 only
|
---|
320 | sum_p2 =sin(theta)*d_theta*d_phi*sum(bin_p2(i,3:(phiDiv+1))); %and pol. 2
|
---|
321 | sum_single =sin(theta)*d_theta*d_phi*sum(bin_single(i,3:(phiDiv+1)));
|
---|
322 | sum_single_p1 = sin(theta)*d_theta*d_phi*sum(bin_single_p1(i,3:(phiDiv+1)));
|
---|
323 | sum_single_p2 = sin(theta)*d_theta*d_phi*sum(bin_single_p2(i,3:(phiDiv+1)));
|
---|
324 | I3(i)=d_phi*(bin(i,2)+2*sum(bin(i,3:(phiDiv+1))));
|
---|
325 | bin(i,phiDiv+3)= sin(theta)*d_theta*d_phi*bin(i,2)+ 2*sum0;
|
---|
326 |
|
---|
327 | bin_p1(i,phiDiv+3) = sin(theta)*d_theta*d_phi*bin_p1(i,2)+ 2*sum_p1;% /* pol. 1 only */
|
---|
328 | bin_p2(i,phiDiv+3) = sin(theta)*d_theta*d_phi*bin_p2(i,2)+ 2*sum_p2; % /* and pol. 2 */
|
---|
329 | bin_single(i,phiDiv+3) = sin(theta)*d_theta*d_phi*bin_single(i,2) + 2*sum_single;
|
---|
330 | bin_single_p1(i,phiDiv+3) = sin(theta)*d_theta*d_phi*bin_single_p1(i,2) + 2*sum_single_p1;
|
---|
331 | bin_single_p2(i,phiDiv+3) = sin(theta)*d_theta*d_phi*bin_single_p2(i,2) + 2*sum_single_p2;
|
---|
332 | end
|
---|
333 |
|
---|
334 | % now integrate over theta
|
---|
335 | sum1=sum(bin(1:ThetaDiv,phiDiv+3));
|
---|
336 | sum3=sum(bin_p1(1:ThetaDiv,phiDiv+3));
|
---|
337 | sum4=sum(bin_p2(1:ThetaDiv,phiDiv+3));
|
---|
338 | sum2=sum(bin_single(1:ThetaDiv,phiDiv+3));
|
---|
339 |
|
---|
340 | 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*/
|
---|
341 | 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*/
|
---|
342 | percent=sum1*0.1*lgr_mm/(n_e_t*1.602177e-19*(gam-1)*0.511*1e+6);
|
---|
343 | disp(['This is equal to ' num2str(percent) ' of the energy in the bunch'])
|
---|
344 | disp(['Energy in polarization 1= ' num2str(sum3*0.1*lgr_mm) ' J'])
|
---|
345 | disp(['Energy in polarization 2= ' num2str(sum4*0.1*lgr_mm) ' J'])
|
---|
346 | lgr_e_mm = abs(Aed_mm./sin(thet/180*pi)); % variable effective grating length, determined basically by the cone entrance diameter of XXXmm*/
|
---|
347 | lgr_e_mm(lgr_e_mm>lgr_mm)=lgr_mm;
|
---|
348 | sum_tr = bin(:,phiDiv+2) .* 0.1 .* lgr_e_mm(:) .* SolAng; % /*variable grating length, now in J*/
|
---|
349 | sum_tr_single = bin_single(:,phiDiv+2) .* 0.1 .* lgr_e_mm(:) * SolAng;
|
---|
350 | sum_tr_p1 = bin_p1(:,phiDiv+2) .* 0.1 .* lgr_e_mm(:) * SolAng;
|
---|
351 | sum_tr_p2 = bin_p2(:,phiDiv+2) .* 0.1 .* lgr_e_mm(:) * SolAng;
|
---|
352 | sum_tr_single_p1 = bin_single_p1(:,phiDiv+2) .* 0.1 .* lgr_e_mm(:) * SolAng;
|
---|
353 | sum_tr_single_p2 = bin_single_p2(:,phiDiv+2) .* 0.1 .* lgr_e_mm(:) * SolAng;
|
---|
354 |
|
---|
355 | % KV. Precise calculation over phi
|
---|
356 | phiNew=phi(1):((phi(2)-phi(1))*1e-3):phi(end);
|
---|
357 | for i=1:ThetaDiv%theta
|
---|
358 | IntergalPhi=sum(pchip(phi,single_elec(i,:),phiNew))*(phiNew(2)-phiNew(1));
|
---|
359 | I1(i)=f1*f2*f4(i)*IntergalPhi;%interal over phi (0 to phiMax)
|
---|
360 | IntergalPhi=sum(pchip(phi,bin(i,2:(phiDiv+1)),phiNew))*(phiNew(2)-phiNew(1));
|
---|
361 | In(i)=IntergalPhi;%interal over phi (0 to phiMax)
|
---|
362 | end
|
---|
363 | E_total=sum(In)*d_theta*Aed_mm*0.1;%J
|
---|
364 | nu_Hz=c./(lam*1e-3);%Hz
|
---|
365 | SP_spec_1= I1(:).*0.1.*lgr_e_mm(:)./(el_mm/m*1e-3)*2.*(c./nu_Hz(:).^2)*1e+18;%muJ/THz
|
---|
366 | SP_spec_N=In(:) .*0.1.*lgr_e_mm(:)./(el_mm/m*1e-3)*2.*(c./nu_Hz(:).^2)*1e+18;%muJ/THz
|
---|
367 | %% Plots
|
---|
368 | if PlotFunc(2)==1
|
---|
369 | figure(2)
|
---|
370 | clf
|
---|
371 | plot(thet,bin_single(:,phiDiv+2),'k','linew',2)
|
---|
372 | set(gca,'fontsize',16)
|
---|
373 | ylabel('dI_1/d\Omega, [J sr^{-1}cm^{-1}]')
|
---|
374 | xlabel('\Theta, [deg]')
|
---|
375 | title('Singe electron yield')
|
---|
376 | grid on
|
---|
377 | end
|
---|
378 | if PlotFunc(3)==1
|
---|
379 | figure(3)
|
---|
380 | clf
|
---|
381 | plot(thet,sum_tr,'r','linew',2)
|
---|
382 | % thet_rad=thet/180*pi;
|
---|
383 | % polar(thet_rad(:),sum_tr(:));
|
---|
384 | % ylim([0 Inf])
|
---|
385 | % yl = get(gca,'YLim');
|
---|
386 | % set(gca,'YLim', [0 yl(2)]);
|
---|
387 | set(gca,'fontsize',16)
|
---|
388 | ylabel('dI/sin(\Theta)d\Theta, [J]')
|
---|
389 | xlabel('\Theta, [deg]')
|
---|
390 | title('SP spectrum')
|
---|
391 | grid on
|
---|
392 | end
|
---|
393 |
|
---|
394 | if PlotFunc(4)==1
|
---|
395 | figure(4)
|
---|
396 | clf
|
---|
397 | plot(time,Profile,'r','linew',2)
|
---|
398 | set(gca,'fontsize',16)
|
---|
399 | ylabel('Amplitude, [a.u.]')
|
---|
400 | xlabel('t, [ns]')
|
---|
401 | title('Temporal bunch profile')
|
---|
402 | grid on
|
---|
403 | end
|
---|
404 |
|
---|
405 | if PlotFunc(5)==1
|
---|
406 | figure(5)
|
---|
407 | clf
|
---|
408 | semilogy(nu_Hz*1e-12,SP_spec_1,'r','linew',2)
|
---|
409 | hold on
|
---|
410 | % load CTR.mat
|
---|
411 | % semilogy(nu*1e-12,IntergdOm*1e+18,'k','linew',2)
|
---|
412 | set(gca,'fontsize',16)
|
---|
413 | ylabel('Energy, [\muJ/THz]')
|
---|
414 | xlabel('\nu, [THz]')
|
---|
415 | title('SP spectrum (Single electron)')
|
---|
416 | % legend('SP','CTR')
|
---|
417 | grid on
|
---|
418 | xlim([min(nu_Hz) max(nu_Hz)].*1e-12)
|
---|
419 | end
|
---|
420 |
|
---|
421 |
|
---|
422 |
|
---|