view release on metacpan or  search on metacpan

lib/Math/PlanePath/ArchimedeanChords.pm  view on Meta::CPAN

  # loop with for(;;) since $n_lo..$n_hi limited to IV range
  for (my $n = $n_lo; ; $n += 1) {
    my ($nx,$ny) = $self->n_to_xy($n);
    # #### $n
    # #### $nx
    # #### $ny
    # #### hypot: hypot ($x-$nx,$y-$ny)
    if (hypot($x-$nx,$y-$ny) <= 0.5) {
      ### hypot in range ...
      return $n;
    }
    if (hypot($nx,$ny) >= $r_limit) {
      last;
    }
  }

  ### n not found ...
  return undef;
}

  # int (max (0, int(_PI*$r2) - 4*$r));
  #
  #  my $r2 = $r * $r;
  #  my $n_lo =  int (max (0, int(_PI*$r2) - 4*$r));
  #  my $n_hi = $n_lo + 7*$r + 2;
  # ### $r2
  # $n_lo == $n_lo-1 ||




# x,y has radius hypot(x,y), then want the next higher spiral arc which is r
# >= hypot(x,y)+0.5, with the 0.5 being the width of the circle figure on
# the spiral.
#
# The polar angle of x,y is a=atan2(y,x) and frac=a/2pi is the extra away
# from an integer radius for the spiral.  So seek integer k with k+a/2pi >=
# h with h=hypot(x,y)+0.5.
#
#     k + a/2pi >= h
#     k >= h-a/2pi
#     k = ceil(h-a/2pi)
#       = ceil(hypot(x,y) + 0.5 - atan2(y,x)/2pi)
#
#
# circle radius i has circumference 2*pi*i and at most that many N on it
# rectangle corner at radius Rcorner = hypot(x,y)
#
# sum i=1 to i=Rlimit of 2*pi*i = 2*pi/2 * Rlimit*(Rlimit+1)
#                               = pi * Rlimit*(Rlimit+1)
# is an upper bound, though a fairly slack one
#
#
# cf. arc length along the spiral r=a*theta with a=1/2pi
#     arclength = (1/2) * a * (theta*sqrt(1+theta^2) + asinh(theta))
#               = (1/4*pi) * (theta*sqrt(1+theta^2) + asinh(theta))
# and theta = 2*pi*r
#               = (1/4*pi) * (4*pi^2*r^2 * sqrt(1+1/theta^2) + asinh(theta))
#               = pi * (r^2 * sqrt(1+1/r^2) + asinh(theta)/(4*pi^2))
#
# and to compare to the circles formula
#
#               = pi * (r*(r+1) * r/(r+1) * sqrt(1+1/r^2)
#                       + asinh(theta)/(4*pi^2))
#
# so it's smaller hence better upper bound.  Only a little smaller than the
# squaring once get to 100 loops or so.
#
#
# not exact
sub rect_to_n_range {
  my ($self, $x1,$y1, $x2,$y2) = @_;
  ### rect_to_n_range() ...

  my $rhi = 0;
  foreach my $x ($x1, $x2) {
    foreach my $y ($y1, $y2) {
      my $frac = atan2($y,$x) / (2*_PI);  # -0.5 <= $frac <= 0.5
      $frac += ($frac < 0);                #    0 <= $frac < 1
      $rhi = max ($rhi, ceil(hypot($x,$y)+0.5 - $frac) + $frac);
    }
  }
  ### $rhi

  # arc length pi * (r^2 * sqrt(1+1/r^2) + asinh(theta)/(4*pi^2))
  #          = pi*r^2*sqrt(1+1/r^2) + asinh(theta)/4pi
  my $rhi2 = $rhi*$rhi;
  return (0,
          ceil (_PI * $rhi2 * sqrt(1+1/$rhi2)
                + asinh(2*_PI*$rhi) / (4*_PI)));


  # # each loop begins at N = pi*k^2 - 2 or thereabouts
  # return (0,
  #         int(_PI*$rhi*($rhi+1) + 1));
}

1;
__END__

    # my $slope = 2*($t + (-$c1-$s1*$t)*cos($t) + ($c1*$t-$s1)*sin($t));
    # my $dist = ( ($t*cos($t) - $c1) ** 2
    #           + ($t*sin($t) - $s1) ** 2
    #           - 4*_PI*_PI );
    # my $slope = (2*($t*cos($t)-$c1)*(cos($t) - $t*sin($t))
    #          + 2*($t*sin($t)-$s1)*(sin($t) + $t*cos($t)));
    # my $c1 = $t1 * cos($t1);
    # my $s1 = $t1 * sin($t1);
    # my $c1_2 = $c1*2;
    # my $s1_2 = $s1*2;
    # my $t = $t1 + 2*_PI/$t1; # estimate
    # my $ct = cos($t);
    # my $st = sin($t);
    # my $dist = (($t - $ct*$c1_2 - $st*$s1_2) * $t + $t1sqm);
    # my $slope = 2 * (($t*$ct - $c1) * ($ct - $t*$st)
    #                  + ($t*$st - $s1) * ($st + $t*$ct));
    #
    # my $sub = $dist/$slope;
    # $t -= $sub;

# use constant _A => 1 / (2*_PI);

lib/Math/PlanePath/ArchimedeanChords.pm  view on Meta::CPAN

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)] ]

=head2 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 C<atan2> is the same as the C library and gives -pi E<lt>= angle
E<lt>= pi, hence allowing for fracE<lt>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 C<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 C<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 C<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 C<xy_to_n()> and C<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?

=head2 Rectangle to N Range

For the C<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 E<lt>=
pi*(r^2+1) which is a tighter bound and an easy calculation, but it only
holds up to somewhere around r=10^7.)

=head1 SEE ALSO

L<Math::PlanePath>,
L<Math::PlanePath::TheodorusSpiral>,
L<Math::PlanePath::SacksSpiral>

=head1 HOME PAGE

L<http://user42.tuxfamily.org/math-planepath/index.html>

=head1 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 <http://www.gnu.org/licenses/>.

=cut