#!/bin/bash /home/richmit/bin/ruby
# -*- Mode:Ruby; Coding:us-ascii-unix; fill-column:158 -*-
################################################################################################################################################################
##
# @file      gen.rb
# @author    Mitch Richling <https://www.mitchr.me>
# @date      2017-05-16
# @version   VERSION
# @brief     Generate geometry files (SVG or Pov-Ray) representing an Apollonian Gasket.@EOL
# @keywords  fractal geometry SVG Pov-Ray Apollonian Gasket
# @std       Ruby 2.0
# @copyright 
#  @parblock
#  Copyright (c) 2017, Mitchell Jay Richling <https://www.mitchr.me> All rights reserved.
#
#  Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
#
#  1. Redistributions of source code must retain the above copyright notice, this list of conditions, and the following disclaimer.
#
#  2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions, and the following disclaimer in the documentation
#     and/or other materials provided with the distribution.
#
#  3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without
#     specific prior written permission.
#
#  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
#  IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
#  FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
#  SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR
#  TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#  @endparblock
################################################################################################################################################################

#---------------------------------------------------------------------------------------------------------------------------------------------------------------
require 'cmath'
require 'optparse'
require 'optparse/time'

#---------------------------------------------------------------------------------------------------------------------------------------------------------------
# Command line parse and generation parameters

outputFileName = 'apollonianFractal'
$doSVG         = FALSE
$doPOV         = TRUE
r1             = 1.0
r2             = 1.0
r3             = 1.0
$curvLimit     = 20000
$povClipRad    = 0.5
$povMinFakeRad = 0.003

debug = 0
opts = OptionParser.new do |opts|
  opts.banner = "Usage: agen.rb [options]"
  opts.separator ""
  opts.separator "Options:"
  opts.on("-h",        "--help",             "Show this message")                      { puts opts; exit                 }
  opts.on(             "--r1 RADIUS",        "Radius for first circle")                { |v| r1=v.to_f                   }
  opts.on(             "--r2 RADIUS",        "Radius for second circle")               { |v| r2=v.to_f                   }
  opts.on(             "--r3 RADIUS",        "Radius for third circle")                { |v| r3=v.to_f                   }
  opts.on(             "--doSVG T/F",        "Dump SVG file")                          { |v| $doSVG=v.match(/^[YyTt1]/)  }
  opts.on(             "--doPOV T/F",        "Dump Povray include file")               { |v| $doPOV=v.match(/^[YyTt1]/)  }
  opts.on(             "--output BASE",      "Name for output file(s)")                { |v| outputFileName=v    }
  opts.on(             "--crvLim VALUE",     "Drop circles with large curvature")    { |v| $curvLimit=v.to_f           }
  opts.on(             "--povClipRad VALUE", "Drop spheres bigger than this")          { |v| $povClipRad=v.to_f          }
  opts.on(             "--povMinRad VALUE",  "Puff up small spheres to this size")      { |v| $povMinFakeRad=v.to_f       }
  opts.separator ""
  opts.separator "Generate geometry files (SVG or Pov-Ray) representing an Apollonian Gasket.  "
  opts.separator "Inputs are the starting radii of the initial interior circles.  The exterior "
  opts.separator "circle is automatically generated.  Note that not all radii will work, and   "
  opts.separator "we don't do any error checking.  Circle generation stops when radii become   "
  opts.separator "too small.  The Pov-Ray output is in the form of union object with a version "
  opts.separator "declaration of 3.7.  The SVG output is 1 pix wide, so use a converter that   "
  opts.separator "will do the right thing or adjust the width.                                 "
  opts.separator ""
  opts.separator "Examples:"
  opts.separator "   * Classic apollonian fractal starting with three circles of equal size."
  opts.separator "     * agen.rb --doSVG F --doPOV T --output apf2 --r1 1.0    --r1 1.0    --r1 1.0     --crvLim 20000.0 --povClipRad 0.5 --povMinRad 0.003"
  opts.separator "     * agen.rb --doSVG T --doPOV F --output apf2 --r1 1000.0 --r1 1000.0 --r1 1000.0  --crvLim 1.0"
  opts.separator "   * The tow eye apollonian fractal starting with two bit circles of equal size centered in the enclosing circle."
  opts.separator "     * agen.rb --doSVG T --doPOV F --output apf3 --r1 1000.0 --r1 1000.0 --r1 666.6666666666  --crvLim 1.0"
end
opts.parse!(ARGV)

if( !($doSVG || $doPOV)) then
  puts("ERROR: No output will be produced (see --doPOV & --doSVG)")
  exit
end

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

$circles = Hash.new   # The primary data structure with all the geometry
$quads   = Array.new  # Used for holding lists of tangent circles

#---------------------------------------------------------------------------------------------------------------------------------------------------------------
# Given three circles, produce the two tangent ones -- if possible.  No error checking, so be careful.  We only use this for the enclosing circle.  Method is
# that of Descartes.
#
# https://en.wikipedia.org/wiki/Descartes'_theorem

def newCircles (circle1, circle2, circle3)
  k1, k2, k3 = circle1[0], circle2[0], circle3[0]
  c1, c2, c3 = circle1[1], circle2[1], circle3[1]
  # Compute the inner and outer circle radius
  tmp1 = k1+k2+k3
  tmp2 = 2*Math::sqrt(k1*k2+k2*k3+k3*k1)
  iRad = tmp1+tmp2
  oRad = tmp1-tmp2
  if(iRad < oRad) then
    iRad, oRad = oRad, iRad
  end
  # Compute the coordinates
  tmp1 = c1*k1+c2*k2+c3*k3
  tmp2 = 2*CMath::sqrt(k1*k2*c1*c2+k2*k3*c2*c3+k3*k1*c1*c3)

  circles = [ [iRad, (tmp1+tmp2)/iRad ],
              [oRad, (tmp1-tmp2)/oRad ]
            ]
  return circles  
end

#---------------------------------------------------------------------------------------------------------------------------------------------------------------
# Given four circles, find an alternate solution for the first one.  This is a matter of finding a second root to a polynomial once you already know one.

def otherSol (circle1, circle2, circle3, circle4)
  k1, k2, k3, k4 = circle1[0], circle2[0], circle3[0], circle4[0]
  c1, c2, c3, c4 = circle1[1], circle2[1], circle3[1], circle4[1]
  newK1 = 2*(k2+k3+k4) - k1
  newC1 = (2*(k2*c2+k3*c3+k4*c4) - k1*c1) / newK1
  return [ newK1, newC1 ]
end

#---------------------------------------------------------------------------------------------------------------------------------------------------------------
# Add a circle to our master list.  We return a hash for the circle if the add was successful.  It won't be successful if the circle was already on the list
# or if the circle radius was too small.  Note the custom has key -- the idea is to call circles the same if radius and coordinates match to 6 digits.

def addAcircle(circle)
  theKey = sprintf("%0.5f/%0.5f/%0.5f", circle[0], circle[1].real, circle[1].imag);
  if( ($circles.member?(theKey)) || (circle[0].abs>$curvLimit) ) then
    return nil
  else
    $circles[theKey] = circle
    return theKey
  end
end

#---------------------------------------------------------------------------------------------------------------------------------------------------------------
# Create first quad element

tmp = (r1*r1+r1*r3+r1*r2-r2*r3)/(r1+r2)

$quads.push(Array.new)

$quads[0].push(addAcircle([ 1/r1, Complex(    0,                             0) ]))
$quads[0].push(addAcircle([ 1/r2, Complex(r1+r2,                             0) ]))
$quads[0].push(addAcircle([ 1/r3, Complex(  tmp, Math::sqrt((r1+r3)**2-tmp**2)) ]))

outsideCircle = addAcircle(newCircles(*($circles.values_at(*($quads[0]))))[1])

$quads[0].unshift(outsideCircle)

#---------------------------------------------------------------------------------------------------------------------------------------------------------------
# Compute circles

while(curQuad = $quads.shift) do
  a, b, c, d = curQuad
  tmp=addAcircle(otherSol($circles[a], $circles[b], $circles[c], $circles[d])); if(!tmp.nil?) then $quads.push([tmp, b, c, d]); end
  tmp=addAcircle(otherSol($circles[d], $circles[a], $circles[b], $circles[c])); if(!tmp.nil?) then $quads.push([tmp, a, b, c]); end
  tmp=addAcircle(otherSol($circles[b], $circles[a], $circles[c], $circles[d])); if(!tmp.nil?) then $quads.push([tmp, a, c, d]); end
  tmp=addAcircle(otherSol($circles[c], $circles[a], $circles[b], $circles[d])); if(!tmp.nil?) then $quads.push([tmp, a, b, d]); end
end

#---------------------------------------------------------------------------------------------------------------------------------------------------------------
# Dump SVG

if($doSVG) then
  open(outputFileName + '.svg', 'w') do |file|
    file.puts("<svg width='#{1}px' height='#{1}px'>")
    $circles.each do |key, circle|
      k, c = circle
      file.puts("    <circle cx='#{500+c.real}' cy='#{500+c.imag}' r='#{(1/k).abs}' fill='none' stroke='red'/>")
    end
    file.puts('</svg>')
  end
end

#---------------------------------------------------------------------------------------------------------------------------------------------------------------
# Dump Povray include file

if($doPOV) then
  open(outputFileName + '.inc', 'w') do |file|
    file.puts('#version 3.7;')
    file.puts('#declare apollonianfractal = union {')
    ox,oy = $circles[outsideCircle][1].rect
    $circles.each do |key, circle|
      k, c = circle
      r = (1/k).abs
      if(r < $povMinFakeRad) then
        r = $povMinFakeRad
      end
      if(r<=$povClipRad) then
        file.puts("  sphere {<#{c.real-ox},#{c.imag-oy},#{0}> #{r}}")
      end
    end
    file.puts('}')
  end
end