-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathraio_unico.m
175 lines (133 loc) · 6.11 KB
/
raio_unico.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
%versão 1.1: encontrado erro no loop: fonte e receptor na mesma
%profundidade estava gerando wrong output.
% update: 23-05-17
function d = raio_unico(mod,S,R)
%=================REDEFININDO AS VARIAVEIS========================%
dx=mod(2,2); % Comprimento de cada célula na dimensão x
dz=mod(1,2); % comprimento de cada célula na dimensão z
nx = mod(2,3); % Número de células na dimensão x
nz = mod(1,3); % Número de células na dimensão z
o = mod(:,1); % Origem da malha
%%=====================COEF. ANGULAR E LINEAR%======================%
% CASO 1 %
% Caso em que a FONTE está MAIS a ESQUERDA e MAIS ACIMA que o RECEPTOR.
% xm e zm são representam os ultimos pontos multiplos de dx e dz,
% mais próximos da fonte, respectivamente.
% OBS:.Até o momento só valem para estes 2 casos
if S(1) <= R(1) && S(2) <= R(2)
m = (R(2)-S(2))/(R(1)-S(1)); % Coef. angular
b = S(2) - m*S(1); % Coef. linear
xm = (fix(S(1)/dx) +1)*dx; % xs<xr
zm = (fix(S(2)/dz) +1)*dz; % zs<zr
%%Essas expressões de xm e zm só valem para
%%essa configuração fonte-receptor.
%%é importante notar isto, pois um par de
%%expressão errado pode extrapolar o raio
%%para a direção errada, resultando em
%%pontos além da fonte ou do receptor.
points = [S;R]; % Inicializando a matriz points com as coordenadas
% da fonte e do receptor.
% No fim das interações abaixo, essa matriz conterá
% todos os pontos do receptor à fonte que interceptam a
% malha.
i=3; % Inicializando i com '3' porque a matriz points
% já possui duas linhas. O próximo ponto se encontrará
% na linha 3.
while xm < R(1) %
points(i,:) = [xm zrr(xm,m,b)];
xm = xm+dx; %Contador para sair do loop
i = i+1; %index da matriz de pontos
end
while zm < R(2) %
points(i,:) = [xrr(zm,m,b) zm];
zm = zm+dz; %Contador para sair do loop
i = i+1; %index da matriz de pontos
end
%=========================FINAL DO CASO 1=================================
%======================COEF. ANGULAR E LINEAR%======================%
%==============================CASO 2===============================%
% Caso em que a FONTE está MAIS a ESQUERDA e MAIS ABAIXO que o RECEPTOR.
elseif S(1)<=R(1) && S(2)>=R(2)
m = (R(2)-S(2))/(R(1)-S(1)); %Coef. angular
b = S(2) - m*S(1); %Coef. linear
xm = (fix(S(1)/dx) +1)*dx; % xs<xr
zm = (fix(S(2)/dz))*dz; % zs>zr
points = [S;R];
i=3;
while xm < R(1) %
p1 = [xm zrr(xm,m,b)];
points(i,:) = p1;
xm = xm+dx; %Contador para sair do loop
i = i+1; %index da matriz de pontos
end
while zm > R(2) %
p2 = [xrr(zm,m,b) zm];
points(i,:) = p2;
zm = zm-dz; %Contador para sair do loop
i = i+1; %index da matriz de pontos
end
%=====================COEF. ANGULAR E LINEAR%======================%
%===========================CASO 3=================================%
% Caso em que a FONTE está MAIS a DIREITA e MAIS ACIMA que o RECEPTOR.
elseif S(1)>=R(1) && S(2)<=R(2)
m = (R(2)-S(2))/(R(1)-S(1)); % Coef. angular
b = S(2) - m*S(1); % Coef. linear
xm = (fix(S(1)/dx))*dx; % xs>xr
zm = (fix(S(2)/dz) +1)*dz; % zs<zr
points = [S;R];
i=3;
while xm > R(1) %
p1 = [xm z(xm,m,b)];
points(i,:) = p1;
xm = xm-dx; %Contador para sair do loop
i = i+1; %index da matriz de pontos
end
while zm < R(2) %
p2 = [x(zm,m,b) zm];
points(i,:) = p2;
zm = zm+dz; %Contador para sair do loop
i = i+1; %index da matriz de pontos
end
else
%=====================COEF. ANGULAR E LINEAR%======================%
%%============================CASO 4===============================%
% Caso em que a FONTE está MAIS a DIREITA e MAIS ABAIXO que o RECEPTOR.
m = (R(2)-S(2))/(R(1)-S(1)); % Coef. angular
b = S(2) - m*S(1); % Coef. linear
xm = (fix(S(1)/dx))*dx; % xs>xr
zm = (fix(S(2)/dz))*dz; % zs>zr
points = [S;R];
i=3;
while xm > R(1) %
p1 = [xm z(xm,m,b)];
points(i,:) = p1;
xm = xm-dx; %Contador para sair do loop
i = i+1; %index da matriz de pontos
end
while zm > R(2) %
p2 = [x(zm,m,b) zm];
points(i,:) = p2;
zm = zm-dz; %Contador para sair do loop
i = i+1; %index da matriz de pontos
end
%=======================FINAL DO CASO 4 =======================%
end
% O cálculo começa a partir daqui, o código anterior foi apenas pra levar em conta
% as diferentes situações de posições relativas entre fonte e receptor.
points = sortrows(points,1); % Organizando em ordem crescente de x
% todos os pontos do raio que interceptam
% a malha.
%================ CRIANDO UM VETOR DE ZEROS ================================%
d= zeros(1,nx*nz); % Esse vetor tem 1 linha para cada raio de uma mesma fonte
% e o número de colunas é o número de células
% na malha
for i=1:length(points(:,1))-1
pm = (points(i,:)+points(i+1,:))/2; % Ponto médio auxiliar: utilizado
% para se obter o índice
% da célula na qual o raio está
% passando.
n = cell_numberrr(indexrr(pm,dx,dz),nx); % Índice da célula.
d(n) = distancerr(points(i,:),points(i+1,:)); % Cálculo do tamanho do raio
% dentro da n-ésima célula.
end
end