1 /* spanarc.i
2  * Compute equally spaced points along a circular arc.
3  */
4 /* Copyright (c) 2013, David H. Munro.
5  * All rights reserved.
6  * This file is part of yorick (http://yorick.sourceforge.net).
7  * Read the accompanying LICENSE file for details.
8  */
9 
10 func spanarc(a, b, c, n, cont=)
11 /* DOCUMENT spanarc(a, b, c, n)
12  *   return N points equally spaced along circular arc ABC.  Each of A,
13  *   B, and C is a point in M>=2 dimensional space.  That is, A = [x,y]
14  *   or [x,y,z], etc.  A, B, and C may also have trailing dimensions in
15  *   order to generate points along several arcs simultaneously.  The
16  *   result has dimensions dimsof(A,B,C)-by-N.
17  *   With non-nil, non-zero cont= keyword, the points will fall along the
18  *   circle from C to A continuing the arc ABC and completing the circle.
19  * SEE ALSO: span, spanl, bezier
20  */
21 {
22   if (cont) {  /* invert b through midpoint of ac */
23     m = 0.5*(a + c);
24     bm = m - b;
25     am = m - a;
26     b = ((am*am)(-,sum,..) / (bm*bm)(-,sum,..)) * bm + m;
27     m = a;  a = c;  c = m;  /* swap a and c */
28   }
29   abc = [a, b, c] - [b, c, a];
30   s = (abc*abc)(-,sum,..);
31   /* get p = vector perpendicular to ac toward b
32    * and q = exterior angle from ab to bc
33    */
34   ac = abc(..,3);
35   p = abc(..,2);
36   q = (abc(..,1)*p)(-,sum,..);
37   q = acos(q / sqrt(s(..,1)*s(..,2)));
38   s = s(..,3);
39   p -= (p*ac)(-,sum,..)/s * ac;
40   /* arc starts in direction q at a, finishes in direction -q at b */
41   q -= (((q+q)/(n-1.))(..,-:1:n-1))(..,cum)(..,zcen);  /* step angles */
42   cs = [cos(q), sin(q)](..,cum,);
43   s = sqrt(s);
44   cs *= s / cs(..,0,1);
45   ac *= 1./s;
46   p *= 1./sqrt((p*p)(-,sum,..));
47   abc = a + ac*cs(..,1) + p*cs(..,2);
48   abc(..,0) = c;  /* abc(..,1) = a already guaranteed */
49   return abc;
50 }
51