Changeset 490 in ETALON for CLIO

Timestamp:

Apr 18, 2016, 12:36:31 PM (8 years ago)

Author:

hodnevuc

Message:

 

File:

• CLIO/CTR+gfw/GFWnew/GFWnew.m (modified) (5 diffs)

CLIO/CTR+gfw/GFWnew/GFWnew.m

@@ -1 @@
-  function  [nu_Hz,SP_spec_1,SP_spec_N,lam,E_total,sum_tr,thet]=GFWnew(var)
-
-%  single_elec
-% close all; clear all; clc
+ %SP space charge model, spectrum calculation
+ % close all; clear all; clc
+ % var='PrSelect=1; P_file=''../p20.mat''; PlotFunc=[1 1 1 1 1]; ';
+function  [nu_Hz,SP_spec_1,SP_spec_N,lam,E_total,sum_tr,thet]=GFWnew(var)
+%% Output parameters
+%nu_Hz -- frequency in Hz
+%SP_spec_1 -- single electron yield from frequency, J/THz
+%SP_spec_N -- whole spectrum from frequency, J/THz
+%lam -- wavelength in mm
+%E_total -- total energy emitted, J
+%sum_tr -- spectrum from Theta,  J 
+%thet -- observation angle, Thetadegrees
 %% Configuration of calculation
 %PrSelect == what profile to use: 1 --from file; 2 - assymmetric gaussian
@@ -12 +20 @@
 %   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];
@@ -46 +53 @@
 ThetaDiv=10;%steps in theta
 phiDiv=3;%steps in phi
-
+saveflag=0;%set to 1 if you want to save model parameters
 if (~isempty(var))
-strread(var,'%s','delimiter',' ')
+strread(var,'%s','delimiter',';')
 eval(var);
 end %if
 
 %save all parameters
+if (saveflag)
 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
-      disp('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+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);
+ save([ name 'SP.mat'])
 end
 
@@ -106 +90 @@
 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'])
+numax=c/(el_mm*(1/beta-cos(ThetaStart))/m*1e-3);%max freq in Hz
+
+ if (saveflag)
+save([ name '_additional_SP.mat'])
+ end
+ %% 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:(2*numax)).*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
+if isnan(sum(Profile))
+    error('Increase pulselength or decrease grating pitch')
+end
+%% 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>=numax*1e-12,1);%xxx 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
 %% Estimation grating factor R^2, coherent & incoherent integrals
 %%
@@ -213 +249 @@
 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)
-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