(* Copyright (c) 1992 The Geometry Center; University of Minnesota
   1300 South Second Street;  Minneapolis, MN  55454, USA;
   
This file is part of CirclePack. CirclePack is free software; you can
redistribute it and/or modify it only under the terms given in the
file COPYING, which you should have received along with this file.
This and other software may be obtained via anonymous ftp from
geom.umn.edu; email: software@geom.umn.edu. *)

(* Author: Oliver Goodman *)
(* $Id: CirclePack.m,v 1.4 1992/07/17 16:09:52 oag Exp $ *)

(* Mathematica package for calculating circle packings. The input 
   is a triangulation of a surface and a list of required intersection
   angles. CirclePack can then be used to find the packing having 
   the given triangulation as its nerve and the specified intersection
   angles. The surface can be any closed orientable surface. 
*)

BeginPackage["CirclePack`",{"Triangulation`"}]

CirclePack::usage = "CirclePack[triang] finds a circle packing with
nerve equal to the triangulation. CirclePack[triang, opts] with the
option Intersects -> {c1, c2, ...} finds the packing with circles
intersecting at angle ArcCos[c1] along edges[triang][[1]] etc. 
CirclePack returns a CirclePacking object. The following functions 
may be applied to such an object: radii, lengths, angles."

Intersects::usage = "An optional parameter for CirclePack[]. See
the help on CirclePack for details."

Steps::usage = "An option for CirclePack[] which controls how many 
steps are taken in trying to find a circle packing numerically before 
giving up."

nerve::usage = "nerve[circlepacking] returns a Triangulation object
corresponding to the nerve of the packing."

radii::usage = "radii[circlepacking] returns a list of the radii 
of the circles in the packing. Circles are listed in the order
given by vertices[nerve[circlepacking]]."

lengths::usage = "lengths[circlepacking] returns the lengths of 
edges of the nerve of the packing. Edges are listed in the order
given by edges[nerve[circlepacking]]."

angles::usage = "angles[circlepacking] returns the angles at each
corner of the nerve of the packing. Angles are listed in the order
given by corners[nerve[circlepacking]]."

Format[_CirclePacking] = "CirclePacking"

Begin["`Private`"]

nerve[c_CirclePacking]:= c[[1]]
radii[c_CirclePacking]:= c[[2]]
lengths[c_CirclePacking]:= c[[3]]
angles[c_CirclePacking]:= c[[4]]

(* An obvious way of working with functions of several variables is
   to pass and return lists to and from the function. 
   A much neater method is to pass a function, x say, as a parameter 
   where x[i] is to be the ith input to the function. 
   We may nonetheless wish to apply such a function to a list so that 
   f[{a,b,c}] is equivalent to f[x] where x[1]=a, x[2]=b and x[3]=c. 
   "listForm[f]" overloads f with an appropriate definition of f[_List]. 
   "listForm[f,n]" assumes that f[x] is a function defined on 1,...,n
   and sets f[x_List] to {f[x][1], ... }. 
*)
listForm[f_Symbol]:= (f[l_List]:= f[l[[#]]&])
listForm[{f_Symbol,n_Integer}]:= (f[l_List]:= Array[f[l[[#]]&], n])
listForm[x_, y__]:= (listForm[x]; listForm[y])

(* Three versions of the cosine rule, one for each of the three geometries *)
cosrule[-1][a_,b_,c_,cos_] = (Cosh[a]==Cosh[b] Cosh[c] - Sinh[b] Sinh[c] cos)
cosrule[ 0][a_,b_,c_,cos_] = (a^2==b^2+c^2-2 b c cos)
cosrule[ 1][a_,b_,c_,cos_] = (Cos[a]==Cos[b] Cos[c] + Sin[b] Sin[c] cos)

(* define functions for getting lengths and angles of triangles *)
Off[Solve::ifun]
Do[
   gl[k,b_,c_,ct_]  = Solve[cosrule[k][a,b,c,-ct],{a}][[1,1,2]];
   ga[k,a_,b_,c_]   = ArcCos[Solve[cosrule[k][a,b,c,ct],{ct}][[1,1,2]]];
   Dgl[k,b_,c_,ct_] = D[gl[k,b,c,ct],b];
   D1ga[k,a_,b_,c_] = D[ga[k,a,b,c],a];
   D3ga[k,a_,b_,c_] = D[ga[k,a,b,c],c],
   {k,-1,1}
]
On[Solve::ifun]

(* give definitions for the symbolic derivatives of gl and ga *)
Derivative[0,1,0,0][gl][k_,a_,b_,ct_] = Dgl[k,a,b,ct]
Derivative[0,0,1,0][gl][k_,a_,b_,ct_] = Dgl[k,b,a,ct]
Derivative[0,1,0,0][ga][k_,a_,b_,c_] = D1ga[k,a,b,c]
Derivative[0,0,1,0][ga][k_,a_,b_,c_] = D3ga[k,a,c,b]
Derivative[0,0,0,1][ga][k_,a_,b_,c_] = D3ga[k,a,b,c]

pi = N[Pi]

(* Calling cpFunctions generates a set of functions associated with 
   finding circle packings on the given triangulated surface.
*)
cpFunctions[S_, opts_:{}]:= 
Module[{el, aa, kv, DLDR, DKDL},
Block[{gl, ga, Dgl, D1ga, D3ga, a, b, c, k, avsum, ints},
   ints = Intersects /. opts /. Intersects -> Table[1.,{ne[S]}];
   k = Max[-1,Min[1,euler[S]]];
   Do[
      {a,b} = vN[S,#]& /@ edges[S][[i]];
      el[r_][i] = gl[k,r[a],r[b],ints[[i]]], 
      {i,ne[S]}
   ];
   avsum = Transpose[
      Table[({a,b,c} = corners[S][[i]]; 
         aa[l_][i] = ga[k,l[eN[S,{b,c}]],l[eN[S,{a,c}]],l[eN[S,{a,b}]]];
         IdentityMatrix[nv[S]][[vN[S,a]]]), 
         {i,nc[S]}
      ]
   ];
   kv[l_] = 2 pi - avsum.Array[aa[l],{nc[S]}];
   DLDR[r_] = Array[D[el[r][#1],r[#2]]&, {ne[S],nv[S]}];
   DKDL[l_] = Transpose[Array[D[kv[l],l[#]]&, {ne[S]}]];
   listForm[{el,ne[S]}, {aa,nc[S]}, kv, DKDL, DLDR];
   {el, aa, kv, DLDR, DKDL}
]]

norm[v_]:= Sqrt[v.v]

(* finds a vector u such that mx.u is close to v but behaves well when
   mx is singular or nearly so. If mx is singular an exact solution is 
   impossible but nearlySolve finds a close vector. If mx is nearly 
   singular nearlySolve avoids jumping off towards infinty ! 
*)
nearlySolve[mx_,v_]:= Block[{o1, d, o2, dinv, m},
   {o1, d, o2} = SingularValues[mx];
   m = Max[Abs /@ d]/100;
   dinv = (If[#>m, 1/#, 0])& /@ d;
   Transpose[o2].DiagonalMatrix[dinv].o1.v
]

(* iter applies one step of a Newton-Raphson approximation to 
   find a solution to the circle packing problem. It uses nearlySolve.
   The effect should be to allow us to find a solution even when 
   the jacobian is singular or almost singular. 
*)
iter[r_List]:= Block[{l,k,d,s},
   k = kv[l = el[r]];
   If[norm[k] == 0., Return[r]];
   d = DKDL[l].DLDR[r];
   s = nearlySolve[d,k];
   While[norm[Im[kv[el[r-s]]]]> 0, s = s/2]; (* reduce step if out of bounds *)
   r - s
]

cpSolve::initl = "For the given triangulation InitialRadii should be a 
vector with `1` numerical components."

Options[CirclePack] = {Steps -> 20}
CirclePack[S_Triangulation, opts_:{}]:= 
Block[{el, aa, kv, DLDR, DKDL, rl},
   rl = InitialRadii /. opts /. InitialRadii -> Table[1.,{nv[S]}];
   If[Length[rl]=!=nv[S], Message[cpSolve::initl, nv[S]];Return[$Failed]];
   {el, aa, kv, DLDR, DKDL} = cpFunctions[S, opts];
   rl = FixedPoint[iter, rl, Steps /. opts /. Options[CirclePack]];
   CirclePacking[S, rl, el = el[rl], aa[el]]
]

End[]
EndPackage[]

