1、The Single Helix225The Single HelixR. Robert Hentzel Scott Williams Iowa State University Ames, IA 50011 Advisor: Stephen J. WillsonIntroductionWe present an iterative algorithm that finds all points of intersection between a finite, infinitely thin helix a

2、nd an infinite plane, ordered along the helix. The algorithm has constant space complexity and has time complexity proportional to the number of intersections and the logarithm of the desired precision.Assumptions The helix is of finite length. The helix is infinitely thin. The plane is infinite in

3、extent. The helix has constant radius. No more than one helix needs to be considered simultaneously; if multiple helices must be modeled, the routine may be run sequentially oneach.Input and OutputInput to the program consists of: the endpoints of the central axis of the helix, one point which is on

4、 the helix, the winding number of the helix, the handedness of the helix,更多數(shù)學(xué)建模資料請關(guān)注微店店鋪“數(shù)學(xué)建模學(xué)習(xí)交流”/RHO6PSpA226The UMAP Journal16.3 the specification of a plane, and the desired precision of the solutions.The winding number of a helix is the number of times that a point traveling

5、on the helix makes a complete circle while its projection on the helix axis advances by one unit. In other words, it is the number of coils of the helix contained in one distance unit parallel to its axis. Figure 1 illustrates a helix with a winding number of 21.510.5010.510.500-0.5-0.5-1-1Fi

6、gure 1. A helix with winding number of 2.0 (note the axes!).Output consists of an ordered list of intersections, giving for each the distance along the axis of the helix, the coordinates of the intersection point, and the distance from the helix point to the plane (which will be less than the prescr

7、ibed precision.Helicoidal Normal FormWe will do all of our computations with geometric figures (helices and planes) which have been put into helicoidal normal form (HNF), described below. Doing so takes advantage of the symmetries of helices and of the infiniteness of planes to simplify the resultin

8、g calculations and increase the speed of calculating intersection points. The transformation is invertible,The Single Helix227so the coordinates of intersection points obtained may be transformed back into the original system.Definition of Helicoidal Normal FormIn HNF, helices and planes have the fo

9、llowing properties: The axis extends from (0, 0, 0) to (0, 0, 2L), where L is the length of the original helix and is the original winding number of the helix. The point (1, 0, 0) is on the helix; that is, the coordinate system is oriented so that the lowest point on the helix is on the x-axis and t

10、he radius isunity. The winding number of the helix is 1/2. The planes normal vector is normalized to unit length. The helix is right-handed. The z-component of the planes normal vector is nonnegative.Thus, a normalized helix may be parameterized asr(t) = 1,(t) = t,z(t) = tin radial coordinates or as

11、x(t) = cos t,y(t) = sin t,z(t) = tin Cartesian coordinates. The equation of a plane is given byax + by + cz + d = 0,with the triple (a, b, c) representing a unit normal vector, so thata2 + b2 + c2 = 1andc 0,as per the definition of the helicoidal normal form.Transformation to Helicoidal Normal FormT

12、he helix and plane are specified by the user with the following data: the endpoints of the symmetry axis of the helix:x0 = (x0 , y0 , z0 ), any point on the helix, p = (xp, yp, zp);x1 = (x1 , y1 , z1 );228The UMAP Journal16.3the winding number of the helix, 0;any (a, b, c, d) quadruplet representing

13、 the plane, so that the normal vectorn is (a, b, c); andthe handedness of the helix.The transformation to HNF consists of seven steps: translation of one axis endpoint to the origin,rotation of coordinates to bring the other endpoint to the z-axis,rotation of coordinates to eliminate any initial pha

14、se of thehelix, space inversion to ensure right-handedness of helix,scaling of coordinates to normalize the helix radius,scaling of coordinates to normalize the helix winding number, and normalizing the normal vector of the plane.The details of each step follow..5.6.7.Translation of one axis

15、endpoint to the originWe translate the coordinate system to bring x0 (one end of the helixs symmetry axis) to (0, 0, 0):x0x0 x0 = 0x1x1 x0pp x0 .We can find the effect on the representation of the plane as follows. If x is originally on the plane, thenn x + d = 0,so the translated point x x0 satisfi

16、es the conditionn (x x0 ) + (d + n x0 ) = 0.Accordingly, we transformd d + n x0 .The Single Helix229Rotation of coordinates to bring the other endpoint to the z-axisTakingx2+ y2 = arctan(x1 /y1 ) , = arctan 11z1,we can construct the rotation matrixcos 0sincos sin 0sincos 0001R = RyzRxy =010cos ,sin

17、0which can be applied to bring x1 to coincide with (0, 0, L), where L is the original length of the helix, given byL =x1 x1 .That is, we transformx0 Rx0 = 0,x1 Rx1 = (0, 0, L),p Rp.We can find the effect on the representation of the plane as follows. If x were originally on the plane, thenn x + d =

18、0.Now the rotated point Rx satisfies the conditionRn Rx + d = 0,so we make the transformationn Rn.Rotation of coordinates to eliminate any initial phaseWe now look at the point p on the helix and note its parametric repre- sentation asrp =x2+ y 2,p = arctan(yp/xp),zp = zp.ppThus, in going downward a

19、long the helixs axis from z = zp to z = 0 (the bottom end), we will go through zp rotations, bringing the initial angular position of the helix (0) to0 = p 2zp.We wish this to coincide with the x-axis, so we apply an additional counter- clockwise rotation of 0 about the z-axis to x1 , x2 , p and the

20、 representationof the plane, using the same techniques as in the previous section.230The UMAP Journal16.3Scaling of coordinates to normalize helix radiusWe can now parameterize our helix asx(t) = rcos t,y(t) = r sin t,z(t) = t,with r = rp =plane, thenx2 + y2. We note that if any point x(t) coincides

21、 with theppa r cos t + b r sin t + ct + d = 0;so the radius-normalized point (cos t, sin t, t) satisfies(ar) cos t + (br) sin t + ct + d = 0,and we can make the transformationa ar,b br,r 1.Space inversion to ensure right-handedness of helixIf the helix is originally left-handed, it can be made right

22、-handed by effecting a spatial inversion of the coordinate system about the xz-plane. This is accomplished by negating the y-coordinate of the planes normal vector:handedness right.b b,Scaling of coordinates to normalize helix winding numberThe final parameter to normalize is the helix winding numbe

23、r. We do this by forcing the helix to advance one rotation per 2 advance along its axis. To compensate, we scale our helix.If the helix originally advanced turns per unit axis advance and had a length of L, then the same number of turns will be made in a length of 2L with a winding number of 1/2. Th

24、us:1 2.L 2L,We can find the effect on the representation of the plane as follows. If a pointx = (x, y, z) is originally on the plane, thenax + by + cz + d = 0.Then, since the transformed point (x, y, z) fulfills the conditioncax + by +(2z)+ d = 0,2we can make the transformation:c c2.The Single Helix

25、231Normalizing the normal vector of the planeWe make the transformationsaa2 + b2 + c2b b,a ,a2 + b2 + c2 da2 + b2 + c2c c,d .a2 + b2 + c2If c 2/2, then the initial arcsin will have an argument larger thanunity and will be nonexistent. The arcsin function will yield two distinct values in 0, 2), both

26、 of which must be checked. If t is a solution of f (t) = 0, then t + n2 must also be a solution for integers n.232The UMAP Journal16.3Difference FunctionThe phrase “the (helix) point t” will be taken to mean the helix pointparameterized by t, namely (cos t, sin t, t) for t 0, L. Since all points x o

27、n the plane satisfyax + by + cz + d = 0,we see that finding the points of intersection between the helix and the plane amounts to finding t such thatf (t) = a cos(t)+ b sin(t)+ ct + d = 0for t 0, L.We name f our difference function since it provides a measure of thedistance between the helix and the

28、 plane for a given t. In fact, f (t) is the perpendicular distance from the helix point t to the plane. Wewill eventually seek intersections by attempting to minimize the absolute value of f(t).Slope of the Difference FunctionThe slope of the difference function at the helix point tis given by:df(x)

29、 dxf t() = a sin t + b cos t + cxtWe note: a2 + b2 (which is equivalent to c 2/2), the slope is every- For c where positive, so the difference function is monotonically increasing. The form of the slope function is that studied in the Transcendental Lemma above, so we know that we can find zeros of

30、the slope function,and hence local extrema of the difference function, analytically. Here the two values for the arcsin correspond to a local maximum and a local minimum of the distance function.Upper and Lower BoundsIt is easy to find bounds on t that contain all of the possible intersection points

31、. At an intersection point, we havea cos t + b sin t + ct + d = 0,so a cos t b sin t d ct =.The Single Helix233Simple calculus, plus using c =a2 + b2, | cos t| 1, and | sin t| 1, showsthat a2+ b 2 a cos t b sin t a2+ b 2,which provides bounds on t: a2 + b2 da2 + b2 d t ub.bccAll intersection points

32、t must lie within these bounds.If c = 0 (the discontinuity in the above formulae), the plane runs exactly parallel to the axis of the helix and would intersect an infinite helix either infinitely often (inside the radius) or never (outside the radius). In thiscase, b may be taken to be and ub to be

33、+, since both values will be truncated by the finiteness of the helix, as mentioned in the next subsection.Intervals and SubintervalsWe define the search interval lower bound B to be the larger of 0 (one end of the helix) and b. We define the search interval upper bound uB to be the smaller of L (th

34、e other end) and ub. All intersections must lie in the search interval B, uB.We divide the interval B, uB into subintervals, broken by the local extrema, both maxima and minima. Each subinterval, except the leftmost and rightmost, is bounded by adjacent local extrema. The leftmost and rightmost are

35、bounded on their “inner” side by a local extremum and on their “outer” side by either B or uB. We can do this because we can apply the Transcendental Lemma to the slope function to calculate two local extrema, one minimum and one maximum. Since the spacings of minima and maxima are both 2, we can no

36、w locate all local extrema within the search interval by repeated addition/subtraction of 2 from the original maximum and minimum (see Figure 2). It may be the case, if c 2/2, that there are no local extrema. Then there is only one subinterval, B, uB, which contains the only intersection. (There can

37、 be only one intersection, since f is monotonically increasing in this case.)Searching the SubintervalsWe now consider each subinterval separately. We let t and ut represent the initial left and right endpoints of the subinterval. We evaluate f ( t) and f (ut) and compare the signs. If the signs agr

38、ee, then there cannot be an intersection in the subinterval; if there were, there would have to be a local extremum within the subinterval, but all local extrema fall on subinterval boundaries by construction. However, if the signs are different, there must234The UMAP Journal16.343210-1-2-3-40510152

39、0t25303540Figure 2. Difference function f (t) with search interval limits shown as the outermost dotted vertical lines, zero shown as a dashed horizontal line, and subinterval boundaries shown as dotted vertical lines. Note the small(er) subintervals at the ends.be an intersection point, by applicat

40、ion of the Intermediate Value Theorem to the continuous function f , which is positive at one endpoint and negative at the other.If endpoints show that no intersection can exist (i.e., both have the same sign), then the subinterval is immediately discarded (this can occur only with the leftmost and

41、rightmost subintervals). Otherwise, we begin a binary search for the intersection point, provided that the endpoint itself is not a point of intersection within the desired precision.We evaluate f at the midpoint of the subinterval ( t + ut)/2 and com- pare its sign with the endpoints. By the Interm

42、ediate Value Theorem, the intersection must lie between the midpoint and whichever endpoint differs in sign from it (they cannot both differ, since they have different signs). We can therefore define a new subinterval equal to either the left or right half of t, ut as appropriate and repeat this pro

43、cess. Eventually, evaluation of f at the midpoint must yield a number within the prescribed accuracy from zero. At this point, the binary search terminates, having located the intersection with sufficient precision.We note that we apply the binary search only in subintervals that are guaranteed to c

44、ontain an intersection by virtue of the Intermediate Valuef(t)The Single Helix235Theorem.When intersections are found, the proper inverse transformations are performed to return the point to its original coordinate system, and it is displayed.After completing a subinterval, processing proceeds with

45、the next and terminates following the rightmost. Because of the left-to-right processing of subintervals, the intersections are generated in order of increasing t along the helix axis.Space and Time ComplexitySpace ComplexityIn our implementation, each subinterval is calculated and searched sub- seq

46、uently, obviating the need to store large arrays of points found or of edges of subintervals. Thus, the storage space needed is independent of the input. For our implementation, it is approximately 300 bytes.Time ComplexityThe time to put the helix and plane into HNF is independent of the input and,

47、 in the current implementation, involves approximately 118 floating- point multiplies/divides and 36 floating-point trigonometric functions and square roots. Floating-point adds, subtracts, and comparisons take negligi- ble time relative to multiplication, division, square roots, and trigonometric f

48、unction calls.The time per intersection is 40 floating-point multiplications/divisions and 8 complicated functions, plus 4 multiplications and 3 complicated func- tions per iteration required. The largest that a subinterval may be is 2, so to divide this by halving into a slice as small as the desired accuracyrequires log2(2/ ) steps. Thus, for each intersection, 40 + 4 log2(2/ ) mul- tiplications and 8 + 3 log2(2/ ) complicated functions are required. For an accuracy of one part in 106, this amounts to 131 mul