1 /* makerotation - Construct rotation from x to y by alpha. */
2 /* XLISP-STAT 2.1 Copyright (c) 1990, by Luke Tierney */
3 /* Additions to Xlisp 2.1, Copyright (c) 1989 by David Michael Betz */
4 /* You may give out copies of this software; for conditions see the */
5 /* file COPYING included with this distribution. */
9 static double inner_product(n
, x
, y
)
15 for (; n
> 0; n
--, x
++, y
++) result
+= *x
* *y
;
19 #define NORM(n, x) (sqrt(inner_product(n, x, x)))
22 make_rotation(int n
, double **rot
, double *x
, double *y
, int use_alpha
, double alpha
)
24 double nx
, ny
, xy
, c
, s
;
27 for (i
= 0; i
< n
; i
++) {
28 for (j
= 0; j
< n
; j
++) rot
[i
][j
] = 0.0;
34 if (nx
== 0.0 || ny
== 0.0) return;
36 for (i
= 0; i
< n
; i
++) x
[i
] /= nx
;
37 for (i
= 0; i
< n
; i
++) y
[i
] /= ny
;
39 xy
= inner_product(n
, x
, y
);
40 c
= (use_alpha
) ? cos(alpha
) : xy
;
41 s
= (use_alpha
) ? sin(alpha
) : sqrt(1 - c
* c
);
43 for (i
= 0; i
< n
; i
++) y
[i
] -= xy
* x
[i
];
46 if (ny
== 0.0) return;
47 for (i
= 0; i
< n
; i
++) y
[i
] /= ny
;
49 for (i
= 0; i
< n
; i
++) {
50 for (j
= 0; j
< n
; j
++)
51 rot
[i
][j
] = x
[j
] * ( x
[i
] * (c
- 1.0) + y
[i
] * s
)
52 + y
[j
] * (- x
[i
] * s
+ y
[i
] * (c
- 1.0));