File Coverage

blib/lib/Math/PlanePath/VogelFloret.pm

Criterion	Covered	Total	%
statement	76	124	61.2
branch	10	32	31.2
condition	14	26	53.8
subroutine	20	28	71.4
pod	11	11	100.0
total	131	221	59.2

# Copyright 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019 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 .


# n_start=>0 to include N=0 at the origin, but that not a documented feature
# yet.

# http://algorithmicbotany.org/papers/#abop
#
# http://www.sciencedirect.com/science/article/pii/0025556479900804
# http://dx.doi.org/10.1016/0025-5564(79)90080-4 Helmut Vogel, "A Better Way
# to Construct the Sunflower Head", Volume 44, Issues 3-4, June 1979, Pages
# 179-189

# http://artemis.wszib.edu.pl/~sloot/2_1.html
#
# http://www.csse.monash.edu.au/publications/2003/tr-2003-149-full.pdf
#     on 3D surfaces of revolution or some such maybe
#     14 Mbytes (or preview with google)

# Count of Zeckendorf bits plotted on Vogel floret.
# Zeckendorf/Fibbinary with N bits makes radial spokes.  cf FibbinaryBitCount
# http://www.ms.unimelb.edu.au/~segerman/papers/sunflower_spiral_fibonacci_metric.pdf
# private copy ?

# closest two for phi are 1 and 4
#     n=1   r=sqrt(1) = 1
#           t=1/phi^2 = 0.381 around
#           x=-.72 y=.68
#     n=4   r=sqrt(4) = 2
#           t=4/phi^2 = 1.527 = .527 around
#           x=-1.97 y=-.337
#     diff angle=4/phi^2 - 1/phi^2 = 3/phi^2 = 3*(2-phi) = 1.14 = .14
#     diff dx=1.25 dy=1.017  hypot=1.61
#     dang = 2*PI()*(5-3*phi)
#     y = sin()
#     x = sin(2*PI()*(5-3*phi))

# Continued fraction
#               1
#     x = k + ------
#             k +  1
#                 ------
#                 k +  1
#                     ---
#                     k + ...
#
#     x = k + 1/x
#     (x-k/2)^2 = 1 + (k^2)/4
#
#         k + sqrt(4+k^2)
#     x = ---------------
#               2
#
#    k       x
#    1    (1+sqrt(5)) / 2
#    2    1 + sqrt(2)
#    3    (3+sqrt(13)) / 2
#    4    2 + sqrt(5)
#    5    (5 + sqrt(29)) / 2
#    6    3 + sqrt(10)
#   2e    e + sqrt(1+e^2)  even



package Math::PlanePath::VogelFloret;
use 5.004;
use strict;
use Carp 'croak';
use Math::Libm 'hypot';

use vars '$VERSION', '@ISA';
$VERSION = 128;
use Math::PlanePath;
@ISA = ('Math::PlanePath');

use Math::PlanePath::Base::Generic
  'is_infinite';
use Math::PlanePath::SacksSpiral;

# uncomment this to run the ### lines
#use Smart::Comments '###';

use constant figure => 'circle';

use constant 1.02; # for leading underscore
use constant _PHI => (1 + sqrt(5)) / 2;
use constant _TWO_PI => 4*atan2(1,0);

# not documented yet ...
use constant rotation_types =>
  { phi   => { rotation_factor => 2 - _PHI(),
               radius_factor   => 0.624239116809924,
               # closest_Ns      => [ 1,4 ],
               # continued_frac  => [ 1,1,1,1,1,... ],
             },
    sqrt2 => { rotation_factor => sqrt(2)-1,
               radius_factor   => 0.679984167849259,
               # closest_Ns      => [ 3,8 ],
               # continued_frac  => [ 2,2,2,2,2,... ],
             },
    sqrt3 => { rotation_factor => sqrt(3)-1,
               radius_factor   => 0.755560810248419,
               # closest_Ns      => [ 3,7 ],
               # continued_frac  => [ 1,2,1,2,1,2,1,2,... ],
             },
    sqrt5 => { rotation_factor => sqrt(5)-2,
               radius_factor   => 0.853488207169303,
               # closest_Ns      => [ 4,8 ],
               # continued_frac  => [ 4,4,4,4,4,4,... ],
             },
  };

use constant parameter_info_array =>
  [
   {
    name      => 'rotation_type',
    type      => 'enum',
    display   => 'Rotation Type',
    share_key => 'vogel_rotation_type',
    choices   => ['phi', 'sqrt2', 'sqrt3', 'sqrt5', 'custom'],
    default   => 'phi',
   },
   {
    name => 'rotation_factor',
    type => 'float',
    type_hint => 'expression',
    display   => 'Rotation Factor',
    description => 'Rotation factor.  If you have Math::Symbolic then this  can be an expression like pi+2*e-phi (constants phi,e,gam,pi), otherwise it should be a plain number.',
    default => - (1 + sqrt(5)) / 2,
    default_expression => '-phi',
    width => 10,
    when_name  => 'rotation_type',
    when_value => 'custom',
   },
   { name           => 'radius_factor',
     display        => 'Radius Factor',
     description    => 'Radius factor, spreading points out to make them non-overlapping.  0 means the default factor.',
     type           => 'float',
     minimum        => 0,
     maximum        => 999,
     page_increment => 1,
     step_increment => .1,
     decimals       => 2,
     default        => 1,
     when_name  => 'rotation_type',
     when_value => 'custom',
   },
  ];

sub x_negative_at_n {
  my ($self) = @_;
  return int(.25 / $self->{'rotation_factor'}) + 1;
}
sub y_negative_at_n {
  my ($self) = @_;
  return int(.5 / $self->{'rotation_factor'}) + 1;
}
sub sumabsxy_minimum {
  my ($self) = @_;
  my ($x,$y) = $self->n_to_xy($self->n_start);
  return abs($x)+abs($y);
}
sub rsquared_minimum {
  my ($self) = @_;
  # starting N=1 at R=radius_factor*sqrt(1), theta=something
  return $self->{'radius_factor'} ** 2;
}
use constant gcdxy_maximum => 0;

sub turn_any_left {  # always left if rot<=0.5
  my ($self) = @_;
  return ($self->{'rotation_factor'} <= 0.5);
}
sub turn_any_right {  # always left if rot<=0.5
  my ($self) = @_;
  return ($self->{'rotation_factor'} > 0.5);
}
use constant turn_any_straight => 0;  # never straight


#------------------------------------------------------------------------------

sub new {
  my $self = shift->SUPER::new (@_);
  ### $self

  my $rotation_type = ($self->{'rotation_type'} ||= 'phi');
  my $defaults = rotation_types()->{$rotation_type}
    || croak 'Unrecognised rotation_type: "',$rotation_type,'"';

  $self->{'radius_factor'} ||= ($self->{'rotation_factor'}
                                ? 1.0
                                : $defaults->{'radius_factor'});
  $self->{'rotation_factor'} ||= $defaults->{'rotation_factor'};
  return $self;
}

# R=radius_factor*sqrt($n)
# R^2 = radius_factor^2 * $n
# avoids sqrt and sin/cos in the main n_to_xy()
#
sub n_to_rsquared {
  my ($self, $n) = @_;
  ### VogelFloret RSquared: $i, $seq->{'planepath_object'}

  if ($n < 0) { return undef; }
  my $rf = $self->{'radius_factor'};
  $rf *= $rf;  # squared

  # don't round BigInt*flonum if radius_factor is not an integer, promote to
  # BigFloat instead
  if (ref $n && $n->isa('Math::BigInt') && $rf != int($rf)) {
    require Math::BigFloat;
    $n = Math::BigFloat->new($n);
  }
  return $n * $rf;
}


sub n_to_xy {
  my ($self, $n) = @_;
  if ($n < 0) { return; }

  my $two_pi = _TWO_PI();

  if (ref $n) {
    if ($n->isa('Math::BigInt')) {
      $n = Math::PlanePath::SacksSpiral::_bigfloat()->new($n);
    }
    if ($n->isa('Math::BigRat')) {
      $n = $n->as_float;
    }
    if ($n->isa('Math::BigFloat')) {
      $two_pi = 2 * Math::BigFloat->bpi;
    }
  }

  my $r = sqrt($n) * $self->{'radius_factor'};

  # take the frac part of 1==circle and then convert to radians, so as not
  # to lose precision in an fmod(...,2*pi)
  #
  my $theta = $n * $self->{'rotation_factor'};    # 1==full circle
  $theta = $two_pi * ($theta - int($theta));  # radians 0 to 2*pi
  return ($r * cos($theta),
          $r * sin($theta));

  # cylindrical_to_cartesian() is only perl code, so may as well sin/cos
  # here directly
  # return (Math::Trig::cylindrical_to_cartesian($r, $theta, 0))[0,1];
}

sub xy_to_n {
  my ($self, $x, $y) = @_;

  # Slack approach just trying all the N values between r-.5 and r+.5.
  #
  # r = sqrt(n)*FACTOR
  # n = (r/FACTOR)^2
  #
  # The target N satisfies N = K * phi + epsilon for integer K.  What's an
  # easy way to find the first integer N >= (r-.5)**2 satisfying -small <= N
  # mod .318 <= +small ?
  #
  my $r = sqrt($x*$x + $y*$y);  # hypot
  my $factor = $self->{'radius_factor'};
  my $n_lo = int( (($r-.6)/$factor)**2 );
  if ($n_lo < 0) { $n_lo = 0; }
  my $n_hi = int( (($r+.6)/$factor)**2 + 1 );
  #### $r
  #### xy: "$x,$y"
  #### $n_lo
  #### $n_hi

  if (is_infinite($n_lo) || is_infinite($n_hi)) {
    ### infinite range, r inf or too big
    return undef;
  }

  # for(;;) loop since "reverse $n_lo..$n_hi" limited to IV range
  for (my $n = $n_hi; $n >= $n_lo; $n--) {
    my ($nx, $ny) = $self->n_to_xy($n);
    ### hypot: "$n ".hypot($nx-$x,$ny-$y)
    if (hypot($nx-$x,$ny-$y) <= 0.5) {
      #### found: $n
      return $n;
    }
  }
  return undef;

  # my $theta_frac = Math::PlanePath::MultipleRings::_xy_to_angle_frac($x,$y);
  # ### assert: 0 <= $frac && $frac < 1
  #
  #   # seeking integer k where (k+theta)*PHIPHI == $r*$r == $n or nearby
  #   my $k = $r*$r / (PHI*PHI) - $theta;
  #
  #   ### $x
  #   ### $y
  #   ### $r
  #   ### $theta
  #   ### $k
  #
  #   foreach my $ki (POSIX::floor($k), POSIX::ceil($k)) {
  #     my $n = int (($ki+$theta)*PHI*PHI + 0.5);
  #
  #     # look for within 0.5 radius
  #     my ($nx, $ny) = $self->n_to_xy($n);
  #     ### $ki
  #     ### n frac: ($ki+$theta)*PHI*PHI
  #     ### $n
  #     ### hypot: hypot($nx-$x,$ny-$y)
  #     if (hypot($nx-$x,$ny-$y) <= 0.5) {
  #       return $n;
  #     }
  #   }
  #   return;
}

# max corner at R
# R+0.5 = sqrt(N) * radius_factor
# sqrt(N) = (R+0.5)/rfactor
# N = (R+0.5)^2 / rfactor^2
#   = (R^2 + R + 1/4) / rfactor^2
#   <= (X^2+Y^2 + X+Y + 1/4) / rfactor^2
#   <= (X(X+1) + Y(Y+1) + 1) / rfactor^2
#
# min corner at R
# R-0.5 = sqrt(N) * radius_factor
# sqrt(N) = (R-0.5)/rfactor
# N = (R-0.5)^2 / rfactor^2
#   = (R^2 - R + 1/4) / rfactor^2
#   >= (X^2+Y^2 - (X+Y)) / rfactor^2      because x+y >= r
#   = (X(X-1) + Y(Y-1)) / rfactor^2

# not exact
sub rect_to_n_range {
  my $self = shift;
  ### VogelFloret rect_to_n_range(): @_
  my ($n_lo, $n_hi) = Math::PlanePath::SacksSpiral->rect_to_n_range(@_);

  my $rf = $self->{'radius_factor'};
  $rf *= $rf; # squared

  # avoid BigInt/flonum if radius_factor is not an integer, promote to
  # BigFloat instead
  if ($rf == int($rf)) {
    $n_hi += $rf-1; # division round upwards
  } else {
    if (ref $n_lo && $n_lo->isa('Math::BigInt')) {
      require Math::BigFloat;
      $n_lo = Math::BigFloat->new($n_lo);
    }
    if (ref $n_hi && $n_lo->isa('Math::BigInt')) {
      require Math::BigFloat;
      $n_hi = Math::BigFloat->new($n_hi);
    }
  }

  $n_lo = int($n_lo / $rf);
  if ($n_lo < 1) { $n_lo = 1; }

  $n_hi = _ceil($n_hi / $rf);

  return ($n_lo, $n_hi);
}

sub _ceil {
  my ($x) = @_;
  my $int = int($x);
  return ($x > $int ? $int+1 : $int);
}

1;
__END__