* Bug #13810 fixed - householder(v, k*v) returned column of %nan. Input parameters...
[scilab.git] / scilab / modules / linear_algebra / macros / householder.sce
1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) 2016 - Samuel GOUGEON
3 //
4 // Copyright (C) 2012 - 2016 - Scilab Enterprises
5 //
6 // This file is hereby licensed under the terms of the GNU GPL v2.0,
7 // pursuant to article 5.3.4 of the CeCILL v.2.1.
8 // This file was originally licensed under the terms of the CeCILL v2.1,
9 // and continues to be available under such terms.
10 // For more information, see the COPYING file which you should have received
11 // along with this program.
12
13 // Example for householder() : Reflection of an object using Householder matrix
14 // ----------------------------------------------------------------------------
15 // (OA) = [0 0 1] is reflected in O into (OB) = [ 1 1 0.3 ]:
16 [n, H] = householder([0 0 1]', [ 1 1 0.3 ]');
17 // "n" is the unit vector orthogonal to the reflecting plane
18
19 // Emitting object (feature from shell demo):
20 u = linspace(0,2*%pi,40);
21 v = linspace(0,2*%pi,20);
22 Xe = (cos(u).*u)'*(1+cos(v)/2)+10;
23 Ye = (u/2)'*sin(v);
24 Ze = (sin(u).*u)'*(1+cos(v)/2);
25
26 // Reflected object:
27 Pe = [ Xe(:)' ; Ye(:)' ; Ze(:)'];
28 Pr = H*Pe;
29 Xr = matrix(Pr(1,:),40,-1);
30 Yr = matrix(Pr(2,:),40,-1);
31 Zr = matrix(Pr(3,:),40,-1);
32
33 // Reflecting plane containing O: n(1).x + n(2).y + n(3).z = 0
34 //   Sampling space:
35 x = linspace(min([Xe(:);Xr(:)]), max([Xe(:);Xr(:)]),20);
36 y = linspace(min([Ye(:);Yr(:)]), max([Ye(:);Yr(:)]),20);
37 [X, Y] = meshgrid(x,y);
38 //   Generating the mirror:
39 deff("z = mirror(x,y,n)","z = -n(1)/n(3)*x - n(2)/n(3)*y")
40 Zm = mirror(X,Y,n);
41
42 // Plotting:
43 clf
44 drawlater
45 f = gcf();
46 f.color_map = [ 0.8 0.8 0.8 ; jetcolormap(100)];
47 surf(Xe,Ye,Ze)
48 surf(X,Y,Zm)
49 surf(Xr,Yr,Zr)
50 a = gca();
51 a.isoview = "on";
52 a.rotation_angles = [74 123];
53 a.children.color_flag = 0;
54 a.children.color_mode = 0;
55 a.children(1).foreground = color("red");
56 a.children(2).foreground = 1;
57 a.children(3).foreground = color("green");
58 drawnow
59 u = []
60 hm = []