source: ETALON/CLIO/GUI_Spectrum_monitor/SEY/QUx.m @ 730

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));