--- name: "Math::PlanePath::ArchimedeanChords" section: "3pm" platform: "linux" package: "libmath-planepath-perl" version: "129-1" release: "Debian 13 (Trixie)" platforms: ["linux"] url: "https://new.carta.tech/Math::PlanePath::ArchimedeanChords" --- # Math::PlanePath::ArchimedeanChords(3pm) ## NAME Math::PlanePath::ArchimedeanChords -- radial spiral chords ## SYNOPSIS ``` use Math::PlanePath::ArchimedeanChords; my $path = Math::PlanePath::ArchimedeanChords->new; my ($x, $y) = $path->n_to_xy (123); ``` ## DESCRIPTION This path puts points at unit chord steps along an Archimedean spiral. The spiral goes outwards by 1 unit each revolution and the points are spaced 1 apart. ``` R = theta/(2*pi) ``` The result is roughly ``` 31 32 30 ... 3 33 29 14 34 15 13 28 50 2 12 16 3 35 2 27 49 1 4 11 17 36 5 0 1 26 48 <- Y=0 10 18 37 6 25 47 -1 9 19 7 8 24 46 38 -2 20 23 39 21 22 45 -3 40 44 41 42 43 ^ -3 -2 -1 X=0 1 2 3 4 ``` X,Y positions returned are fractional. Each revolution is about 2\*pi longer than the previous, so the effect is a kind of 6.28 increment looping. Because the spacing is by unit chords, adjacent unit circles centred on each N position touch but don't overlap. The spiral spacing of 1 unit per revolution means they don't overlap radially either. The unit chords here are a little like the `"TheodorusSpiral"`. But the `"TheodorusSpiral"` goes by unit steps at a fixed right-angle and approximates an Archimedean spiral (of 3.14 radial spacing). Whereas this `"ArchimedeanChords"` is an actual Archimedean spiral (of radial spacing 1), with unit steps angling along that. ## FUNCTIONS See "FUNCTIONS" in Math::PlanePath for behaviour common to all path classes. - "$path = Math::PlanePath::ArchimedeanChords-\>new ()" Create and return a new path object. - "($x,$y) = $path-\>n\_to\_xy ($n)" Return the X,Y coordinates of point number `$n` on the path. `$n` can be any value `"$n >= 0"` and fractions give positions on the chord between the integer points (ie. straight line between the points). `"$n==0"` is the origin 0,0. For `"$n < 0"` the return is an empty list, it being considered there are no negative points in the spiral. - "$n = $path-\>xy\_to\_n ($x,$y)" Return an integer point number for coordinates `"$x,$y"`. Each integer N is considered the centre of a circle of diameter 1 and an `"$x,$y"` within that circle returns N. The unit spacing of the spiral means those circles don't overlap, but they also don't cover the plane and if `"$x,$y"` is not within one then the return is `"undef"`. The current implementation is a bit slow. - "$n = $path-\>n\_start ()" Return 0, the first `$n` on the path. - "$str = $path-\>figure ()" Return "circle". ## FORMULAS ### N to X,Y The current code keeps a position as a polar angle t and calculates an increment u needed to move along by a unit chord. If dist(u) is the straight-line distance between t and t+u, then squared is the hypotenuse ``` dist^2(u) = ((t+u)/2pi*cos(t+u) - t/2pi*cos(t))^2 # X + ((t+u)/2pi*sin(t+u) - t/2pi*sin(t))^2 # Y ``` which simplifies to ``` dist^2(u) = [ (t+u)^2 + t^2 - 2*t*(t+u)*cos(u) ] / (4*pi^2) ``` Switch from cos to sin using the half angle cos(u) = 1 - 2\*sin^2(u/2) in case if u is small then the cos(u) near 1.0 might lose floating point accuracy, and also as a slight simplification, ``` dist^2(u) = [ u^2 + 4*t*(t+u)*sin^2(u/2) ] / (4*pi^2) ``` Then want the u which has dist(u)=1 for a unit chord. The u\*sin(u) part probably doesn't have a good closed form inverse, so the current code is a Newton/Raphson iteration on f(u) = dist^2(u)-1, seeking f(u)=0 ``` f(u) = u^2 + 4*t*(t+u)*sin^2(u/2) - 4*pi^2 ``` Derivative f'(u) for the slope from the cos form is ``` f'(u) = 2*(t+u) - 2*t*[ cos(u) - (t+u)*sin(u) ] ``` And again switching from cos to sin in case u is small, ``` f'(u) = 2*[ u + t*[2*sin^2(u/2) + (t+u)*sin(u)] ] ``` ### X,Y to N A given x,y point is at a fraction of a revolution ``` frac = atan2(y,x) / 2pi # -.5 <= frac <= .5 frac += (frac < 0) # 0 <= frac < 1 ``` And the nearest spiral arm measured radially from x,y is then ``` r = int(hypot(x,y) + .5 - frac) + frac ``` Perl's `"atan2"` is the same as the C library and gives -pi \<= angle \<= pi, hence allowing for frac\<0. It may also be "unspecified" for x=0,y=0, and give +/-pi for x=negzero, which has to be a special case so 0,0 gives r=0. The `"int"` rounds towards zero, so frac\>.5 ends up as r=0. So the N point just before or after that spiral position may cover the x,y, but how many N chords it takes to get around to there is 's not so easily calculated. The current code looks in saved `"n_to_xy()"` positions for an N below the target, and searches up from there until past the target and thus not covering x,y. With `"n_to_xy()"` points saved 500 apart this means searching somewhere between 1 and 500 points. One possibility for calculating a lower bound for N, instead of the saved positions, and both for `"xy_to_n()"` and `"rect_to_n_range()"`, would be to add up chords in circles. A circle of radius k fits pi/arcsin(1/2k) many unit chords, so ``` k=floor(r) pi total = sum ------------ k=0 arcsin(1/2k) ``` and this is less than the chords along the spiral. Is there a good polynomial over-estimate of arcsin, to become an under-estimate total, without giving away so much? ### Rectangle to N Range For the `"rect_to_n_range()"` upper bound, the current code takes the arc length along with spiral with the usual formula ``` arc = 1/4pi * (theta*sqrt(1+theta^2) + asinh(theta)) ``` Written in terms of the r radius (theta = 2pi\*r) as calculated from the biggest of the rectangle x,y corners, ``` arc = pi*r^2*sqrt(1+1/r^2) + asinh(2pi*r)/4pi ``` The arc length is longer than chords, so N=ceil(arc) is an upper bound for the N range. An upper bound can also be calculated simply from the circumferences of circles 1 to r, since a spiral loop from radius k to k+1 is shorter than a circle of radius k. ``` k=ceil(r) total = sum 2pi*k k=1 = pi*r*(r+1) ``` This is bigger than the arc length, thus a poorer upper bound, but an easier calculation. (Incidentally, for smallish r have arc length \<= pi\*(r^2+1) which is a tighter bound and an easy calculation, but it only holds up to somewhere around r=10^7.) ## SEE ALSO Math::PlanePath, Math::PlanePath::TheodorusSpiral, Math::PlanePath::SacksSpiral ## HOME PAGE \ ## LICENSE Copyright 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020 Kevin Ryde This file is part of Math-PlanePath. Math-PlanePath is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 3, or (at your option) any later version. Math-PlanePath is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with Math-PlanePath. If not, see \.