Mercurial > octave
changeset 28385:1888f07317a8
streamtube.m, ostreamtube.m: Minor code cleaning (patch #9916).
| author | Markus Meisinger <chloros2@gmx.de> |
|---|---|
| date | Thu, 02 Apr 2020 12:22:56 +0200 |
| parents | 23e6c897526a |
| children | 8a9a041db1dc |
| files | scripts/plot/draw/ostreamtube.m scripts/plot/draw/streamtube.m |
| diffstat | 2 files changed, 66 insertions(+), 71 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/plot/draw/ostreamtube.m Thu Apr 02 12:17:51 2020 +0200 +++ b/scripts/plot/draw/ostreamtube.m Thu Apr 02 12:22:56 2020 +0200 @@ -161,7 +161,7 @@ for i = 1 : length (xyz) sl = xyz{i}; if (! isempty (sl)) - slx = sl(:, 1); sly = sl(:, 2); slz = sl(:, 3); + slx = sl(:,1); sly = sl(:,2); slz = sl(:,3); mxx(j) = max (slx); mnx(j) = min (slx); mxy(j) = max (sly); mny(j) = min (sly); mxz(j) = max (slz); mnz(j) = min (slz); @@ -179,12 +179,12 @@ num_vertices = rows (sl); if (! isempty (sl) && num_vertices > 2) - usl = interp3 (x, y, z, u, sl(:, 1), sl(:, 2), sl(:, 3)); - vsl = interp3 (x, y, z, v, sl(:, 1), sl(:, 2), sl(:, 3)); - wsl = interp3 (x, y, z, w, sl(:, 1), sl(:, 2), sl(:, 3)); + usl = interp3 (x, y, z, u, sl(:,1), sl(:,2), sl(:,3)); + vsl = interp3 (x, y, z, v, sl(:,1), sl(:,2), sl(:,3)); + wsl = interp3 (x, y, z, w, sl(:,1), sl(:,2), sl(:,3)); vv = sqrt (usl.*usl + vsl.*vsl + wsl.*wsl); - div_sl = interp3 (x, y, z, div, sl(:, 1), sl(:, 2), sl(:, 3)); + div_sl = interp3 (x, y, z, div, sl(:,1), sl(:,2), sl(:,3)); is_singular_div = find (isnan (div_sl), 1, "first"); if (! isempty (is_singular_div)) @@ -211,10 +211,9 @@ cp = cos (phi); sp = sin (phi); - X0 = sl(1, :); - X1 = sl(2, :); - - ## 1st rotation axis + ## 1st streamline segment + X0 = sl(1,:); + X1 = sl(2,:); R = X1 - X0; RE = R / norm (R); @@ -224,10 +223,8 @@ XS0 = rotation (K, RE, cp, sp) + repmat (X0.', 1, num_circum); ## End of first segment - vold = vv(1); - vact = vv(2); - ract = rstart * exp (0.5 * div_sl(2) * norm (R) / vact) * sqrt (vold / vact); - vold = vact; + ract = rstart * exp (0.5 * div_sl(2) * norm (R) / vv(2)) * ... + sqrt (vv(1) / vv(2)); rold = ract; K = ract * KE; XS = rotation (K, RE, cp, sp) + repmat (X1.', 1, num_circum); @@ -249,22 +246,22 @@ for i = 3 : max_vertices - KK = K; + ## Next streamline segment X0 = X1; - X1 = sl(i, :); + X1 = sl(i,:); R = X1 - X0; RE = R / norm (R); ## Tube radius - vact = vv(i); - ract = rold * exp (0.5 * div_sl(i) * norm (R) / vact) * sqrt (vold / vact); - vold = vact; + ract = rold * exp (0.5 * div_sl(i) * norm (R) / vv(i)) * ... + sqrt (vv(i-1) / vv(i)); rold = ract; - ## Project KK onto RE and get the difference in order to calculate the next - ## guiding point - Kp = KK - RE * dot (KK, RE); - K = ract * Kp / norm (Kp); + ## Project KE onto RE and get the difference in order to transport + ## the normal vector KE along the vertex array + Kp = KE - RE * dot (KE, RE); + KE = Kp / norm (Kp); + K = ract * KE; ## Rotate around RE and collect surface patches XS = rotation (K, RE, cp, sp) + repmat (X1.', 1, num_circum); @@ -284,9 +281,9 @@ function N = get_normal1 (X) if ((X(3) == 0) && (X(1) == -X(2))) - N = [- X(2) - X(3), X(1), X(1)]; + N = [(- X(2) - X(3)), X(1), X(1)]; else - N = [X(3), X(3), - X(1) - X(2)]; + N = [X(3), X(3), (- X(1) - X(2))]; endif N /= norm (N); @@ -301,17 +298,17 @@ uy = U(2); uz = U(3); - Y(1, :) = X(1) * (cp + ux * ux * (1 - cp)) + ... - X(2) * (ux * uy * (1 - cp) - uz * sp) + ... - X(3) * (ux * uz * (1 - cp) + uy * sp); + Y(1,:) = X(1) * (cp + ux * ux * (1 - cp)) + ... + X(2) * (ux * uy * (1 - cp) - uz * sp) + ... + X(3) * (ux * uz * (1 - cp) + uy * sp); - Y(2, :) = X(1) * (uy * ux * (1 - cp) + uz * sp) + ... - X(2) * (cp + uy * uy * (1 - cp)) + ... - X(3) * (uy * uz * (1 - cp) - ux * sp); + Y(2,:) = X(1) * (uy * ux * (1 - cp) + uz * sp) + ... + X(2) * (cp + uy * uy * (1 - cp)) + ... + X(3) * (uy * uz * (1 - cp) - ux * sp); - Y(3, :) = X(1) * (uz * ux * (1 - cp) - uy * sp) + ... - X(2) * (uz * uy * (1 - cp) + ux * sp) + ... - X(3) * (cp + uz * uz * (1 - cp)); + Y(3,:) = X(1) * (uz * ux * (1 - cp) - uy * sp) + ... + X(2) * (uz * uy * (1 - cp) + ux * sp) + ... + X(3) * (cp + uz * uz * (1 - cp)); endfunction @@ -366,8 +363,8 @@ %!error <invalid number of inputs> ostreamtube (1,2,3) %!error <invalid number of inputs> ostreamtube (1,2,3,4) %!error <invalid number of inputs> ostreamtube (1,2,3,4,5) -%!error <invalid number of OPTIONS> ostreamtube (1,2,3,4,5,6, [1,2,3]) -%!error <SCALE must be a real scalar . 0> ostreamtube (1,2,3,4,5,6, [1i]) -%!error <SCALE must be a real scalar . 0> ostreamtube (1,2,3,4,5,6, [0]) -%!error <N must be greater than 2> ostreamtube (1,2,3,4,5,6, [1, 1i]) -%!error <N must be greater than 2> ostreamtube (1,2,3,4,5,6, [1, 2]) +%!error <invalid number of OPTIONS> ostreamtube (1,2,3,4,5,6,[1,2,3]) +%!error <SCALE must be a real scalar . 0> ostreamtube (1,2,3,4,5,6,[1i]) +%!error <SCALE must be a real scalar . 0> ostreamtube (1,2,3,4,5,6,[0]) +%!error <N must be greater than 2> ostreamtube (1,2,3,4,5,6,[1,1i]) +%!error <N must be greater than 2> ostreamtube (1,2,3,4,5,6,[1,2])
--- a/scripts/plot/draw/streamtube.m Thu Apr 02 12:17:51 2020 +0200 +++ b/scripts/plot/draw/streamtube.m Thu Apr 02 12:22:56 2020 +0200 @@ -169,7 +169,7 @@ for i = 1 : length (xyz) sl = xyz{i}; if (! isempty (sl)) - slx = sl(:, 1); sly = sl(:, 2); slz = sl(:, 3); + slx = sl(:,1); sly = sl(:,2); slz = sl(:,3); mxx(j) = max (slx); mnx(j) = min (slx); mxy(j) = max (sly); mny(j) = min (sly); mxz(j) = max (slz); mnz(j) = min (slz); @@ -226,10 +226,9 @@ cp = cos (phi); sp = sin (phi); - X0 = sl(1, :); - X1 = sl(2, :); - - ## 1st rotation axis + ## 1st streamline segment + X0 = sl(1,:); + X1 = sl(2,:); R = X1 - X0; RE = R / norm (R); @@ -244,35 +243,34 @@ py = zeros (num_circum, max_vertices); pz = zeros (num_circum, max_vertices); - px(:, 1) = XS0(1, :).'; - py(:, 1) = XS0(2, :).'; - pz(:, 1) = XS0(3, :).'; + px(:,1) = XS0(1,:).'; + py(:,1) = XS0(2,:).'; + pz(:,1) = XS0(3,:).'; - px(:, 2) = XS(1, :).'; - py(:, 2) = XS(2, :).'; - pz(:, 2) = XS(3, :).'; + px(:,2) = XS(1,:).'; + py(:,2) = XS(2,:).'; + pz(:,2) = XS(3,:).'; for i = 3 : max_vertices - KEold = KE; + ## Next streamline segment X0 = X1; - - X1 = sl(i, :); + X1 = sl(i,:); R = X1 - X0; RE = R / norm (R); - ## Project KE onto RE and get the difference in order to calculate the next - ## guiding point - Kp = KEold - RE * dot (KEold, RE); + ## Project KE onto RE and get the difference in order to transport + ## the normal vector KE along the vertex array + Kp = KE - RE * dot (KE, RE); KE = Kp / norm (Kp); K = radius_sl(i) * KE; ## Rotate around RE and collect surface patches XS = rotation (K, RE, cp, sp) + repmat (X1.', 1, num_circum); - px(:, i) = XS(1, :).'; - py(:, i) = XS(2, :).'; - pz(:, i) = XS(3, :).'; + px(:,i) = XS(1,:).'; + py(:,i) = XS(2,:).'; + pz(:,i) = XS(3,:).'; endfor @@ -285,7 +283,7 @@ ## the streamline vertex array "sl" ends function [div_sl_crop, max_vertices] = interp_sl (x, y, z, div, sl) - div_sl = interp3 (x, y, z, div, sl(:, 1), sl(:, 2), sl(:, 3)); + div_sl = interp3 (x, y, z, div, sl(:,1), sl(:,2), sl(:,3)); is_nan = find (isnan (div_sl), 1, "first"); is_inf = find (isinf (div_sl), 1, "first"); @@ -305,9 +303,9 @@ function N = get_normal1 (X) if ((X(3) == 0) && (X(1) == -X(2))) - N = [- X(2) - X(3), X(1), X(1)]; + N = [(- X(2) - X(3)), X(1), X(1)]; else - N = [X(3), X(3), - X(1) - X(2)]; + N = [X(3), X(3), (- X(1) - X(2))]; endif N /= norm (N); @@ -322,17 +320,17 @@ uy = U(2); uz = U(3); - Y(1, :) = X(1) * (cp + ux * ux * (1 - cp)) + ... - X(2) * (ux * uy * (1 - cp) - uz * sp) + ... - X(3) * (ux * uz * (1 - cp) + uy * sp); + Y(1,:) = X(1) * (cp + ux * ux * (1 - cp)) + ... + X(2) * (ux * uy * (1 - cp) - uz * sp) + ... + X(3) * (ux * uz * (1 - cp) + uy * sp); - Y(2, :) = X(1) * (uy * ux * (1 - cp) + uz * sp) + ... - X(2) * (cp + uy * uy * (1 - cp)) + ... - X(3) * (uy * uz * (1 - cp) - ux * sp); + Y(2,:) = X(1) * (uy * ux * (1 - cp) + uz * sp) + ... + X(2) * (cp + uy * uy * (1 - cp)) + ... + X(3) * (uy * uz * (1 - cp) - ux * sp); - Y(3, :) = X(1) * (uz * ux * (1 - cp) - uy * sp) + ... - X(2) * (uz * uy * (1 - cp) + ux * sp) + ... - X(3) * (cp + uz * uz * (1 - cp)); + Y(3,:) = X(1) * (uz * ux * (1 - cp) - uy * sp) + ... + X(2) * (uz * uy * (1 - cp) + ux * sp) + ... + X(3) * (cp + uz * uz * (1 - cp)); endfunction @@ -356,7 +354,7 @@ %!demo %! clf; %! t = 0:.15:15; -%! xyz{1} = [cos(t)' sin(t)' (t/3)']; +%! xyz{1} = [cos(t)', sin(t)', (t/3)']; %! dia{1} = cos(t)'; %! streamtube (xyz, dia); %! grid on;