function D=coef_difraccion(k0,beta,R0,Rn,phi1,phi2,s,n,pol)
%Función que calcula el coeficiente de difracción según la UTD
% k0=constante de onda
% R0=coeficiente de reflexion de fresnel para la faceta 0
% Rn=coeficiente de reflexion de fresnel para la faceta n
% ph1=ángulo de difracción del transmisor en radianes
% phi2=ángulo de difracción del receptor en radianes
% s=distancia del punto de difracción al receptor
% n
% pol=polarización de la onda
L=s*sin(beta)^2;

X1=phi2-phi1; %X-
X2=phi2+phi1; %X+
N11=round((-pi+X1)/2*pi*n); %N- X- -> g-
N12=round((-pi+X2)/2*pi*n); %N- X+ -> g-
N21=round((pi+X1)/2*pi*n);  %N+ X- -> g+
N22=round((pi+X2)/2*pi*n);  %N+ X+ -> g+

g11=1+cos(X1-2*n*pi*N11); %g-(X-)
g12=1+cos(X2-2*n*pi*N12); %g-(X+)
g21=1+cos(X1-2*n*pi*N21); %g+(X-)
g22=1+cos(X2-2*n*pi*N22); %g+(X+)

if(phi1<=(n-1)*pi) %caso A   
   if(abs(X2-pi)<=pi/180) %singularidad en la RSB 
       D1=cot((pi+X1)/(2*n))*FresnelI(k0*L*g21); %Di=D1+D2
       D2=cot((pi-X1)/(2*n))*FresnelI(k0*L*g11);
       D3=-n*sqrt(2*pi*k0*L)*exp(i*pi/4); %Dr=D3+D4  %CAMBIADO DE SIGNO
       D4=cot((pi+X2)/(2*n))*FresnelI(k0*L*g22);
   elseif(abs(X1-pi)<=pi/180) %singularidad en la ISB 
       D1=cot((pi+X1)/(2*n))*FresnelI(k0*L*g21); %Di=D1+D2
       D2=-n*sqrt(2*pi*k0*L)*exp(i*pi/4); %CAMBIADO DE SIGNO
       D3=cot((pi-X2)/(2*n))*FresnelI(k0*L*g12); %Dr=D3+D4
       D4=cot((pi+X2)/(2*n))*FresnelI(k0*L*g22);
   else
       D1=cot((pi+X1)/(2*n))*FresnelI(k0*L*g21); %Di=D1+D2
       D2=cot((pi-X1)/(2*n))*FresnelI(k0*L*g11);
       D3=cot((pi-X2)/(2*n))*FresnelI(k0*L*g12); %Dr=D3+D4
       D4=cot((pi+X2)/(2*n))*FresnelI(k0*L*g22);
   end
else %caso B
   if(abs(X2-(2*n-1)*pi)<=pi/180) %singularidad en la RSB 
       D1=cot((pi+X1)/(2*n))*FresnelI(k0*L*g21); %Di=D1+D2
       D2=cot((pi-X1)/(2*n))*FresnelI(k0*L*g11);
       D3=cot((pi-X2)/(2*n))*FresnelI(k0*L*g12); %Dr=D3+D4
       D4=n*sqrt(2*pi*k0*L)*exp(i*pi/4); %CAMBIADO DE SIGNO
   elseif(abs(X2-pi)<=pi/180) 
       D1=cot((pi+X1)/(2*n))*FresnelI(k0*L*g21); %Di=D1+D2
       D2=cot((pi-X1)/(2*n))*FresnelI(k0*L*g11);
       D3=n*sqrt(2*pi*k0*L)*exp(i*pi/4); %Dr=D3+D4   %CAMBIADO DE SIGNO
       D4=cot((pi+X2)/(2*n))*FresnelI(k0*L*g22);
   elseif(abs(X1+pi)<=pi/180) %singularidad en la ISB
       D1=n*sqrt(2*pi*k0*L)*exp(i*pi/4); %Di=D1+D2  %CAMBIADO DE SIGNO
       D2=cot((pi-X1)/(2*n))*FresnelI(k0*L*g11);
       D3=cot((pi-X2)/(2*n))*FresnelI(k0*L*g12); %Dr=D3+D4
       D4=cot((pi+X2)/(2*n))*FresnelI(k0*L*g22);
   else
       D1=cot((pi+X1)/(2*n))*FresnelI(k0*L*g21); %Di=D1+D2
       D2=cot((pi-X1)/(2*n))*FresnelI(k0*L*g11);
       D3=cot((pi-X2)/(2*n))*FresnelI(k0*L*g12); %Dr=D3+D4
       D4=cot((pi+X2)/(2*n))*FresnelI(k0*L*g22);
   end
end
 
%coef de difraccion
D=-((exp(-i*pi/4))/(2*n*sqrt(2*pi*k0)))*(D1+D2+R0*D3+Rn*D4);
end