function [ang1 ang2 R0 Rn]=angulo_difraccion(ptx,prx,difractor,pol)
%Función que calcula con que ángulo se difracta un rayo proveniente de un
%punto especificado en una esquina difractora concreta hacia otro punto
%              *v3 ptx *v4 prx
%               \     /                             |v1|
%                \   /                              |v2|
%                 \ /                      vectores=|v3|
%  difractor p_dif o------* v1 cara 0               |v4|
%                   \                      material=|mat cara 0|
%                    \                              |mat cara n|
%                     *v2 cara n

global MATERIALES FACETAS;
%Formamos primero los vectores unitarios v1, v2 y v3
ang1=[];
ang2=[];
cambio=0;
p_dif=difractor(1:2);
vectores(3,1:2)=(ptx-p_dif)/norm(ptx-p_dif);
vectores(4,1:2)=(prx-p_dif)/norm(prx-p_dif);
%Primera faceta de la difracción
for i=3:4 %Obtenemos los vectores que forman las facetas de difracción
    p_ini=FACETAS(difractor(i),1:2);
    p_fin=FACETAS(difractor(i),3:4);
    material(i-2)=MATERIALES(FACETAS(difractor(i),7),1);
    %Escogemos como punto inicial el que este en el p_dif
    if(abs(p_dif(1)-p_fin(1))<0.001 && abs(p_dif(2)-p_fin(2))<0.001) %En este caso, se intercambian el final con el inicial
        aux=p_ini;
        p_ini=p_fin;
        p_fin=aux;
    end    
    vectores(i-2,1:2)=(p_fin-p_ini)/norm(p_fin-p_ini); %Almacenamos los vectores normalizados
end
%Ahora los disponemos como v1 y v2, solo necesito colocar el v1 en abcisas
%y rotar el resto
mat_0=MATERIALES(FACETAS(difractor(3),7),1);
mat_n=MATERIALES(FACETAS(difractor(4),7),1);
x1=vectores(1,1);
y1=vectores(1,2);
x2=vectores(2,1);
y2=vectores(2,2);
if(x1==x2 && y1==y2) %Son el mismo, difracción en una sola faceta, v1=v2
    giro=atan2(y1,x1);
    if (giro<0)
        giro=2*pi+giro; %Para rotar en sentido antihorario    
    end
else
    giro1=atan2(y1,x1); %Ángulo que forma con abcisas
    giro2=atan2(y2,x2);
    if ((giro1>0) && (giro2)>0) %Cogemos el ángulo más grande
        if(giro1>giro2)
            giro=giro1;
        else
            giro=giro2;
            cambio=1;
        end
        giro=giro-2*pi; %Para rotar en sentido antihorario
    else
        if ((giro1<0) && (giro2)<0) %Cogemos el ángulo más pequeño en valor absoluto o el menos negativo
            if(giro1>giro2)
                giro=2*pi+giro1;
            else
                giro=2*pi+giro2;
                cambio=1;
            end
        else %Si estan en distintos cuadrantes (distinto signo ya que atan2 esta entre -pi y pi)
            if(giro1>giro2) %giro1 es positivo y giro2 negativo 
                giro1a2=(2*pi+giro2)-giro1;
                giro2a1=2*pi-giro1a2;
            else %giro2 es positivo y giro1 negativo
                giro2a1=(2*pi+giro1)-giro2;
                giro1a2=2*pi-giro2a1;
            end
            if(giro1a2>giro2a1)
                giro=giro1;
                if(giro<0)
                    giro=2*pi+giro;
                end
            else
                giro=giro2;
                if(giro<0)
                    giro=2*pi+giro;
                end
                cambio=1;
            end
        end
    end
end
%En este punto ya tenemos el ángulo que forma la cara 0 con el eje de
%abcisas en positivo
if(cambio==1) %Hay que intercambiar los vectores de la cara 0 y la cara n
    aux=mat_0;
    mat_0=mat_n;
    mat_n=aux;
    vectores(1,1:2)=[x2 y2]; %Antiguo v2, ahora es el v1
    vectores(2,1:2)=[x1 y1]; %Antiguo v1, ahora es el v2
end
            
M=[cos(giro) sin(giro);-sin(giro),cos(giro)]; %Rotamos hacia el eje de abcisas
%Realizamos la traslación de los ejes, para comprobar si v3 tiene la
%dirección correcta, ahora todos los vectores estan girados respecto al eje
%de abcisas

v1=(M*vectores(1,1:2)')'; %Deberia ser el eje de abcisas 
v2=(M*vectores(2,1:2)')';
v3=(M*vectores(3,1:2)')'; %Del punto tx
v4=(M*vectores(4,1:2)')'; %Del punto rx
%Si la dirección es correcta, el ángulo en sentido antihorario del vector
%de nuestros puntos será siempre más pequeño que entre el de v1 y v2

if(x1~=x2 && y1~=y2) %Cuando son paralelos falla
    angv1v2=atan2(v2(2),v2(1)); %Calculamos el ángulo de cada vector con v1
else
    angv1v2=2*pi; %Si no, no salia 0 exacto
end
angv1v3=atan2(v3(2),v3(1));
angv1v4=atan2(v4(2),v4(1));

if(abs(angv1v3)<1e-10) %Evita rayos que atraviesen la pared, ya que el angulo no sale 0 exacto   
    angv1v3=2*pi;
end
if(abs(angv1v4)<1e-10)
    angv1v4=2*pi;
end

if(angv1v2<0) %Se deberia cumplir siempre que no sean paralelos
    angv1v2=2*pi+atan2(v2(2),v2(1));
end
if(angv1v3<0)
    angv1v3=2*pi+atan2(v3(2),v3(1));
end
if(angv1v4<0)
    angv1v4=2*pi+atan2(v4(2),v4(1));
end
%PRUEBAS
% v2
% angv1v2
% angv1v3
% angv1v4
% if (prx(1)==6 && prx(2)==2)%PRUEBAS
%     angv1v3
%     angv1v4
% end

if(angv1v2>angv1v3) %Dirección correcta entre difractor y ptx
    ang1=angv1v3;
end
if(angv1v2>angv1v4) %Dirección correcta entre difractor y prx
    ang2=angv1v4;
end

%Aqui calculamos los valores de reflexión
n=difractor(5);
if(ang1<=ang2) %rayo incidente más próximo a la cara 0 de la cuña que el reflejado
    if(pol==1)
        R0=reflexion_horizontal_fresnel(mat_0,ang1);
        Rn=reflexion_horizontal_fresnel(mat_n,n*pi-ang2);
    else
        R0=reflexion_vertical_fresnel(mat_0,ang1);
        Rn=reflexion_vertical_fresnel(mat_n,n*pi-ang2);
    end
else %rayo reflejado más próximo a la cara 0 de la cuña que el incidente
    if(pol==1)
        R0=reflexion_horizontal_fresnel(mat_0,ang2);
        Rn=reflexion_horizontal_fresnel(mat_n,n*pi-ang1);
    else
        R0=reflexion_vertical_fresnel(mat_0,ang2);
        Rn=reflexion_vertical_fresnel(mat_n,n*pi-ang1);
    end
end
end