.\" Automatically generated by Pod::Man 4.14 (Pod::Simple 3.40) .\" .\" Standard preamble: .\" ======================================================================== .de Sp \" Vertical space (when we can't use .PP) .if t .sp .5v .if n .sp .. .de Vb \" Begin verbatim text .ft CW .nf .ne \\$1 .. .de Ve \" End verbatim text .ft R .fi .. .\" Set up some character translations and predefined strings. \*(-- will .\" give an unbreakable dash, \*(PI will give pi, \*(L" will give a left .\" double quote, and \*(R" will give a right double quote. \*(C+ will .\" give a nicer C++. Capital omega is used to do unbreakable dashes and .\" therefore won't be available. \*(C` and \*(C' expand to `' in nroff, .\" nothing in troff, for use with C<>. .tr \(*W- .ds C+ C\v'-.1v'\h'-1p'\s-2+\h'-1p'+\s0\v'.1v'\h'-1p' .ie n \{\ . ds -- \(*W- . ds PI pi . if (\n(.H=4u)&(1m=24u) .ds -- \(*W\h'-12u'\(*W\h'-12u'-\" diablo 10 pitch . if (\n(.H=4u)&(1m=20u) .ds -- \(*W\h'-12u'\(*W\h'-8u'-\" diablo 12 pitch . ds L" "" . ds R" "" . ds C` "" . ds C' "" 'br\} .el\{\ . ds -- \|\(em\| . ds PI \(*p . ds L" `` . ds R" '' . ds C` . ds C' 'br\} .\" .\" Escape single quotes in literal strings from groff's Unicode transform. .ie \n(.g .ds Aq \(aq .el .ds Aq ' .\" .\" If the F register is >0, we'll generate index entries on stderr for .\" titles (.TH), headers (.SH), subsections (.SS), items (.Ip), and index .\" entries marked with X<> in POD. Of course, you'll have to process the .\" output yourself in some meaningful fashion. .\" .\" Avoid warning from groff about undefined register 'F'. .de IX .. .nr rF 0 .if \n(.g .if rF .nr rF 1 .if (\n(rF:(\n(.g==0)) \{\ . if \nF \{\ . de IX . tm Index:\\$1\t\\n%\t"\\$2" .. . if !\nF==2 \{\ . nr % 0 . nr F 2 . \} . \} .\} .rr rF .\" ======================================================================== .\" .IX Title "Math::PlanePath::ArchimedeanChords 3pm" .TH Math::PlanePath::ArchimedeanChords 3pm "2021-01-23" "perl v5.32.0" "User Contributed Perl Documentation" .\" For nroff, turn off justification. Always turn off hyphenation; it makes .\" way too many mistakes in technical documents. .if n .ad l .nh .SH "NAME" Math::PlanePath::ArchimedeanChords \-\- radial spiral chords .SH "SYNOPSIS" .IX Header "SYNOPSIS" .Vb 3 \& use Math::PlanePath::ArchimedeanChords; \& my $path = Math::PlanePath::ArchimedeanChords\->new; \& my ($x, $y) = $path\->n_to_xy (123); .Ve .SH "DESCRIPTION" .IX Header "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. .PP .Vb 1 \& R = theta/(2*pi) .Ve .PP The result is roughly .PP .Vb 10 \& 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 .Ve .PP 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. .PP 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. .PP The unit chords here are a little like the \f(CW\*(C`TheodorusSpiral\*(C'\fR. But the \&\f(CW\*(C`TheodorusSpiral\*(C'\fR goes by unit steps at a fixed right-angle and approximates an Archimedean spiral (of 3.14 radial spacing). Whereas this \&\f(CW\*(C`ArchimedeanChords\*(C'\fR is an actual Archimedean spiral (of radial spacing 1), with unit steps angling along that. .SH "FUNCTIONS" .IX Header "FUNCTIONS" See \*(L"\s-1FUNCTIONS\*(R"\s0 in Math::PlanePath for behaviour common to all path classes. .ie n .IP """$path = Math::PlanePath::ArchimedeanChords\->new ()""" 4 .el .IP "\f(CW$path = Math::PlanePath::ArchimedeanChords\->new ()\fR" 4 .IX Item "$path = Math::PlanePath::ArchimedeanChords->new ()" Create and return a new path object. .ie n .IP """($x,$y) = $path\->n_to_xy ($n)""" 4 .el .IP "\f(CW($x,$y) = $path\->n_to_xy ($n)\fR" 4 .IX Item "($x,$y) = $path->n_to_xy ($n)" Return the X,Y coordinates of point number \f(CW$n\fR on the path. .Sp \&\f(CW$n\fR can be any value \f(CW\*(C`$n >= 0\*(C'\fR and fractions give positions on the chord between the integer points (ie. straight line between the points). \&\f(CW\*(C`$n==0\*(C'\fR is the origin 0,0. .Sp For \f(CW\*(C`$n < 0\*(C'\fR the return is an empty list, it being considered there are no negative points in the spiral. .ie n .IP """$n = $path\->xy_to_n ($x,$y)""" 4 .el .IP "\f(CW$n = $path\->xy_to_n ($x,$y)\fR" 4 .IX Item "$n = $path->xy_to_n ($x,$y)" Return an integer point number for coordinates \f(CW\*(C`$x,$y\*(C'\fR. Each integer N is considered the centre of a circle of diameter 1 and an \f(CW\*(C`$x,$y\*(C'\fR within that circle returns N. .Sp The unit spacing of the spiral means those circles don't overlap, but they also don't cover the plane and if \f(CW\*(C`$x,$y\*(C'\fR is not within one then the return is \f(CW\*(C`undef\*(C'\fR. .Sp The current implementation is a bit slow. .ie n .IP """$n = $path\->n_start ()""" 4 .el .IP "\f(CW$n = $path\->n_start ()\fR" 4 .IX Item "$n = $path->n_start ()" Return 0, the first \f(CW$n\fR on the path. .ie n .IP """$str = $path\->figure ()""" 4 .el .IP "\f(CW$str = $path\->figure ()\fR" 4 .IX Item "$str = $path->figure ()" Return \*(L"circle\*(R". .SH "FORMULAS" .IX Header "FORMULAS" .SS "N to X,Y" .IX Subsection "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 .PP .Vb 2 \& 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 .Ve .PP which simplifies to .PP .Vb 1 \& dist^2(u) = [ (t+u)^2 + t^2 \- 2*t*(t+u)*cos(u) ] / (4*pi^2) .Ve .PP 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, .PP .Vb 1 \& dist^2(u) = [ u^2 + 4*t*(t+u)*sin^2(u/2) ] / (4*pi^2) .Ve .PP 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 .PP .Vb 1 \& f(u) = u^2 + 4*t*(t+u)*sin^2(u/2) \- 4*pi^2 .Ve .PP Derivative f'(u) for the slope from the cos form is .PP .Vb 1 \& f\*(Aq(u) = 2*(t+u) \- 2*t*[ cos(u) \- (t+u)*sin(u) ] .Ve .PP And again switching from cos to sin in case u is small, .PP .Vb 1 \& f\*(Aq(u) = 2*[ u + t*[2*sin^2(u/2) + (t+u)*sin(u)] ] .Ve .SS "X,Y to N" .IX Subsection "X,Y to N" A given x,y point is at a fraction of a revolution .PP .Vb 2 \& frac = atan2(y,x) / 2pi # \-.5 <= frac <= .5 \& frac += (frac < 0) # 0 <= frac < 1 .Ve .PP And the nearest spiral arm measured radially from x,y is then .PP .Vb 1 \& r = int(hypot(x,y) + .5 \- frac) + frac .Ve .PP Perl's \f(CW\*(C`atan2\*(C'\fR is the same as the C library and gives \-pi <= angle <= pi, hence allowing for frac<0. It may also be \*(L"unspecified\*(R" 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 \f(CW\*(C`int\*(C'\fR rounds towards zero, so frac>.5 ends up as r=0. .PP 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. .PP The current code looks in saved \f(CW\*(C`n_to_xy()\*(C'\fR positions for an N below the target, and searches up from there until past the target and thus not covering x,y. With \f(CW\*(C`n_to_xy()\*(C'\fR points saved 500 apart this means searching somewhere between 1 and 500 points. .PP One possibility for calculating a lower bound for N, instead of the saved positions, and both for \f(CW\*(C`xy_to_n()\*(C'\fR and \f(CW\*(C`rect_to_n_range()\*(C'\fR, would be to add up chords in circles. A circle of radius k fits pi/arcsin(1/2k) many unit chords, so .PP .Vb 3 \& k=floor(r) pi \& total = sum \-\-\-\-\-\-\-\-\-\-\-\- \& k=0 arcsin(1/2k) .Ve .PP 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? .SS "Rectangle to N Range" .IX Subsection "Rectangle to N Range" For the \f(CW\*(C`rect_to_n_range()\*(C'\fR upper bound, the current code takes the arc length along with spiral with the usual formula .PP .Vb 1 \& arc = 1/4pi * (theta*sqrt(1+theta^2) + asinh(theta)) .Ve .PP Written in terms of the r radius (theta = 2pi*r) as calculated from the biggest of the rectangle x,y corners, .PP .Vb 1 \& arc = pi*r^2*sqrt(1+1/r^2) + asinh(2pi*r)/4pi .Ve .PP The arc length is longer than chords, so N=ceil(arc) is an upper bound for the N range. .PP 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. .PP .Vb 3 \& k=ceil(r) \& total = sum 2pi*k \& k=1 \& \& = pi*r*(r+1) .Ve .PP 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.) .SH "SEE ALSO" .IX Header "SEE ALSO" Math::PlanePath, Math::PlanePath::TheodorusSpiral, Math::PlanePath::SacksSpiral .SH "HOME PAGE" .IX Header "HOME PAGE" .SH "LICENSE" .IX Header "LICENSE" Copyright 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020 Kevin Ryde .PP This file is part of Math-PlanePath. .PP Math-PlanePath is free software; you can redistribute it and/or modify it under the terms of the \s-1GNU\s0 General Public License as published by the Free Software Foundation; either version 3, or (at your option) any later version. .PP Math-PlanePath is distributed in the hope that it will be useful, but \&\s-1WITHOUT ANY WARRANTY\s0; without even the implied warranty of \s-1MERCHANTABILITY\s0 or \s-1FITNESS FOR A PARTICULAR PURPOSE.\s0 See the \s-1GNU\s0 General Public License for more details. .PP You should have received a copy of the \s-1GNU\s0 General Public License along with Math-PlanePath. If not, see .