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
e-REdING. Biblioteca de la Escuela Superior de Ingenieros de Sevilla.
