-
Notifications
You must be signed in to change notification settings - Fork 1
/
calc_rvo.m
51 lines (42 loc) · 1.26 KB
/
calc_rvo.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
%==========================================================================
%
%
% input :
%
% output :
%
% Siqi Li, SMAST
% yyyy-mm-dd
%
% Updates:
%
%==========================================================================
function rvo = calc_rvo(u, v, dx, dy)
[nx, ny] = size(u);
% xc=zeros(nx-1,ny-1);
% yc=zeros(nx-1,ny-1);
rvo1=zeros(nx-1,ny-1);
for i=1:nx-1
for j=1:ny-1
% xc(i,j)=(x(i,j)+x(i+1,j)+x(i,j+1)+x(i+1,j+1))/4;
% yc(i,j)=(y(i,j)+y(i+1,j)+y(i,j+1)+y(i+1,j+1))/4;
% u_right=(u(i+1,j)+u(i+1,j+1))/2;
% u_left=(u(i,j+1)+u(i,j))/2;
% v_top=(v(i+1,j+1)+v(i,j+1))/2;
% v_bottom=(v(i,j)+v(i+1,j))/2;
% rvo(i,j)=(u_right-v_top-u_left+v_bottom)/dd;
u_bottom=(u(i,j)+u(i+1,j))/2;
u_top=(u(i+1,j+1)+u(i,j+1))/2;
v_right=(v(i+1,j)+v(i+1,j+1))/2;
v_left=(v(i,j+1)+v(i,j))/2;
% rvo1(i,j)=(u_bottom+v_right-u_top-v_left)/d;
rvo1(i,j)=(u_bottom-u_top)/dx+(v_right-v_left)/dy;
end
end
% Interpolate rvo back to the orginal grid points.
rvo = nan(nx, ny);
rvo(2:end-1, 2:end-1) = rvo1(1:end-1, 1:end-1) + ...
rvo1(2:end, 1:end-1) + ...
rvo1(1:end-1, 2:end) + ...
rvo1(2:end, 2:end);
rvo = rvo / 4;