Main program(test_tri_contour.m)
%--------------------------------------------
% Test of Equipotential contour
% for triangular grids
%--------------------------------------------
clear all
clf
%---tri_da
tri_da=[...
1 1 6 2
2 2 6 7
3 3 2 7
4 3 7 8
5 3 8 9
6 3 9 4
7 4 9 5
8 6 11 12
9 7 6 12
10 7 12 13
11 7 13 8
12 8 13 14
13 8 14 10
14 8 10 9
15 11 15 16
16 11 16 12
17 12 16 17
18 12 17 13
19 13 17 18
20 13 18 14]
%----xy_da
xy_da=[...
1 0.0 0.0
2 1.0 0.0
3 2.0 0.0
4 3.0 0.0
5 4.0 0.0
6 0.0 1.0
7 1.0 1.0
8 1.7 1.3
9 3.0 1.0
10 2.4 1.6
11 0.0 1.75
12 0.65 1.8
13 1.25 2.0
14 1.66 2.34
15 0.0 2.0
16 0.518 2.069
17 1.0 2.268
18 1.414 2.586]
%-----phi_da
phi_da =[...
337.44
349.24
373.28
390.21
401.31
314.53
326.45
334.44
382.00
346.79
253.59
259.90
268.44
268.89
223.56
224.39
224.99
225.77]
%---------------------------
% y_scale = 1.8
% tri_num_prnt = 0;
% pnt_num_prnt = 0;
% %
% tri_grid(tri_da, xy_da, y_scale, tri_num_prnt, pnt_num_prnt)
%----------------------------------------
% call subroutine tri_cont
%----------------------------------------
y_scale = 1.8
tri_num_prnt = 1;
pnt_num_prnt = 1;
%
tri_cont(tri_da, xy_da, phi_da, y_scale, tri_num_prnt, pnt_num_prnt)
%==============================================
Subroutine tri_cont.m
function dummy=tri_cont(tri_da, xy_da, phi_da, y_scale, tri_num_prnt, pnt_num_prnt)
%
[n_tr, n] = size(tri_da);
[n_pt, n] = size(xy_da);
nmax = tri_da(1,1);
x = xy_da(:,2);
y = xy_da(:,3);
phi = phi_da;
%
%tri_num_prnt = 0;
%pnt_num_prnt = 0;
%
xmin = min(x);
xmax = max(x);
xcen = 0.5 * (xmin + xmax);
ymin = min(y*y_scale);
ymax = max(y*y_scale);
ycen = 0.5 * (ymin + ymax);
%
phimin = min(phi) + 0.01;
phimax = max(phi) - 0.01;
dx = xmax - xmin;
dy = ymax - ymin;
n_cont = 20;
dphi = (phimax - phimin)/n_cont;
kmax = n_cont;
s = phimin:dphi:phimax;
if dx < dy, xmin = xcen - dy / 2;
xmax = xcen + dy / 2;
end
if dx > dy, ymin = ycen - dx / 2;
ymax = ycen + dx / 2;
end
clf;
hold off;
clc;
axis([xmin, xmax, ymin, ymax]);
%
m = 0;
title('contour plot');
hold on;
%------- draw the border line
plot([xmin, xmax], [ymin, ymin])
plot([xmin, xmax], [ymax, ymax])
plot([xmin, xmin], [ymin, ymax])
plot([xmax, xmax], [ymin, ymax])
%
deltax = 0.1;
deltay = 0.1;
for k = 1:n_tr
for j = 1:3
p = tri_da(k, j+1);
xx(j) = x(p);
yy(j) = y(p);
ff(j) = phi(p);
end
xx(4) = xx(1);
yy(4) = yy(1);
ff(4) = ff(1);
%------ plot the grid mesh
plot(xx, yy*y_scale, 'k-', 'linewidth', 1.5)
%
phi_min = min([ff(1), ff(2), ff(3)]);
phi_max = max([ff(1), ff(2), ff(3)]);
for kv = 1:kmax
if phi_min <= s(kv) & s(kv) <= phi_max,
m = 0;
for i = 1:3
if ( s(kv) - ff(i) ) * (s(kv) - ff(i+1) ) <= 0,
m = m + 1;
if phi(i+1) == phi(i),
alpha = 0.5;
end
if phi(i+1) ~= phi(i),
alpha = (s(kv) - ff(i))/(ff(i+1) - ff(i));
end
xp(m) = alpha*xx(i+1) + (1-alpha)*xx(i);
yp(m) = alpha*yy(i+1) + (1-alpha)*yy(i);
end
if m == 2,
plot([xp(1),xp(2)], [yp(1)*y_scale, yp(2)*y_scale], ':', 'linewidth', 3);
break
end
end
end
end
%
xcen = sum(xx(1:3))/3;
ycen = sum(yy(1:3))/3;
if tri_num_prnt == 1
text(xcen-deltax, (ycen-deltay)*y_scale, int2str(k))
end
end
%
if pnt_num_prnt == 1
for n = 1:n_pt
text( x(n), y(n)*y_scale, ['(', int2str(n), ')'])
end
end
axis('off')
Trying to push the limits of App Inventor!
Snippets and
Tutorials from
Pura Vida Apps by
Taifun.