(* :Name: JuliaSet` *) (* :Author: Mark McClure, April 1998 *) (* :Copyright: Copyright Mark McClure, 1998 *) (* :Package Version: 1.1 *) (* :Mathematica Version: 2.2 *) (* :Summary: Generates images of Julia sets using several variants of the inverse iteration algorithm. *) BeginPackage["JuliaSet`"] Needs["Utilities`FilterOptions`"] JuliaSet::usage = "JuliaSet is a package containing several implementations of inverse iteration algorithms for the generation of Julia sets." Julia::usage = "Julia[formula,var] generates the Julia set of the rational function defined by formula, using the modified inverse iteration algorithm." JuliaSimple::usage = "JuliaSimple[c,depth] generates 2^depth points of the Julia set for z^2 + c using a simple deterministic algorithm." JuliaModified::usage = "JuliaModified[c] generates the Julia set for z^2 + c using a modified determinstic algorithm." JuliaStochastic::usage = "JuliaStochastic[c,numpoints] generates numpoints points of the Julia set for z^2 + c using a stochastic algorithm." Resolution::usage = "Resolution is a positive integer valued Option of Julia and JuliaModified which controls the level of detail of the generated set." Bound::usage = "Bound is a positive real valued Option of Julia which specifies the maximum absolute vale of points plotted in the Julia set of a rational function." OneBias::usage = "OneBias is an Option of JuliaStochastic which determines the probability of choosing the positive inverse of z^2 + c. OneBias must be a real number stricly between 0 and 1." Begin["`Private`"] (* Set the Options *) (* For JuliaSimple *) JuliaSimpleOptions = Join[{AspectRatio -> Automatic, Axes -> False, PlotStyle -> {AbsolutePointSize[0.4]}}, Options[ListPlot]]; keywords = Complement[Union[First /@ JuliaSimpleOptions], {DisplayFunction, DefaultFont}]; vals = keywords /. JuliaSimpleOptions; special = {DisplayFunction :> $DisplayFunction, DefaultFont :> $DefaultFont}; Options[JuliaSimple] = Union[Apply[Rule,Transpose[{keywords,vals}],{1}], special]; (* For JuliaModified *) Options[JuliaStochastic] = Join[{OneBias -> .5}, Options[JuliaSimple]]; Options[JuliaModified] = Join[{Resolution -> 200}, Options[JuliaSimple]]; Options[Julia] = Join[{Bound -> 4}, Options[JuliaModified]]; (* The functions *) (* Error Messages *) JuliaSet::argu = "`1` called with `2`; `3` expected." JuliaSet::badNumber = "Number expected at position 1 in `1`." JuliaSet::badInt = "Positive Integer expected at position 2 in `1`." JuliaSet::badOpt = "Options expected as optional arguments in `1`." JuliaSet::optx = "Unknown Option `1`." Julia::badFormula = "Rational function of degree at least two expected at position 1 in `1`." Julia::noInverses = "Can't find the the inverse functions." OneBias::outOfRange = "OneBias must be stricly between 0 and 1. Using 0.5." Resolution::outOfRange = "Resolution must be a positive integer. Using 200." Bound::outOfRange = "Bound must be a positive number. Using 4." JuliaSimple[c_, depth_, opts___] := Module[ {invImage, points, valid = First /@ Options[JuliaSimple]}, Scan[If[!MemberQ[valid, First[#]], Message[JuliaSet::optx, ToString[First[#]]]]&, Flatten[{opts}] ]; invImage := ( N[{1,-1} Sqrt[#-c]]& /@ # // Flatten ) &; points = {Re[#], Im[#]}& /@ Nest[invImage, {1}, depth]; ListPlot[points, FilterOptions[ListPlot,opts], FilterOptions[ListPlot,Sequence @@ Options[JuliaSimple]] ] ] /; ( NumberQ[N[c]] || Message[JuliaSet::badNumber, JuliaSimple, c] ) && ( (IntegerQ[depth] && depth > 0 ) || Message[JuliaSet::badInt, JuliaSimple, depth] ) && ( And @@ Map[OptionQ, {opts}] || Message[JuliaSet::badOpt,JuliaSimple] ) (* Special cases for JuliaSimple *) JuliaSimple[] := Message[JuliaSet::argu, "JuliaSimple", "0 arguments", "2 arguments are"] JuliaSimple[c_] := Message[JuliaSet::argu, "JuliaSimple", "1 argument", "2 arguments are"] JuliaModified[c_, opts___] := Module[ {invImage,reducedInvImage,pointsSoFar,res, valid = First /@ Options[JuliaModified]}, Scan[If[!MemberQ[valid, First[#]], Message[JuliaSet::optx, ToString[First[#]]]]&, {opts} ]; res = Resolution /. {opts} /. Options[JuliaModified]; If[!(IntegerQ[res] && res>0), {Message[Resolution::outOfRange]; res = Resolution /. Options[JuliaModified]}]; invImage := ( N[{1,-1} Floor[res Sqrt[#-c]]/res]& /@ # // Flatten )&; reducedInvImage[points_] := Module[{newPoints}, newPoints = Complement[invImage[points], pointsSoFar]; pointsSoFar = Union[newPoints, pointsSoFar]; newPoints]; pointsSoFar = Nest[invImage, {1}, 5]; FixedPoint[reducedInvImage, pointsSoFar]; ListPlot[{Re[#], Im[#]}& /@ pointsSoFar, FilterOptions[ListPlot,opts], FilterOptions[ListPlot,Sequence @@ Options[JuliaModified]] ] ] /; (NumberQ[c] || Message[JuliaSet::badNumber, JuliaModified, c]) && ( And @@ Map[OptionQ, {opts}] || Message[JuliaSet::badOpt,JuliaModified] ) (* Special case for JuliaModified *) JuliaModified[] := Message[JuliaSet::argu, "JuliaModified", "0 arguments", "1 argument is"] Julia[formula_, var_Symbol, opts___] := Module[ {z0,inverses,inverseRules,funcs,image, reducedImage,pointsSoFar, res,bound, valid = First /@ Options[Julia]}, Scan[If[!MemberQ[valid, First[#]], Message[JuliaSet::optx, ToString[First[#]]]]&, {opts} ]; res = Resolution /. {opts} /. Options[Julia]; If[!(IntegerQ[res] && res>0), {Message[Resolution::outOfRange]; res = Resolution /. Options[Julia]}]; bound = Bound /. {opts} /. Options[Julia]; If[!(bound>0), {Message[Bound::outOfRange]; bound = Bound /. Options[Julia]}]; inverseRules = NSolve[#==formula,var]; If[And @@ ListQ /@ inverseRules, If[Length[inverseRules] > 1, inverses = var /. inverseRules; funcs = Function[anInverse,N[Floor[anInverse res]/res] &] /@ inverses; image = Flatten[Through[funcs[#]]& /@ #,1]&; If[PolynomialQ[formula,var], reducedImage[points_] := Module[{newPoints}, newPoints = Complement[image[points], pointsSoFar]; pointsSoFar = Union[newPoints, pointsSoFar]; newPoints], reducedImage[points_] := Module[{newPoints}, newPoints = Complement[image[points], pointsSoFar]; newPoints = Select[newPoints,(N[Abs[#]]<=bound)&]; pointsSoFar = Union[newPoints, pointsSoFar]; newPoints] ]; z0 = N[Nest[funcs[[1]][funcs[[2]][#]]&, N[Pi],10]]; pointsSoFar = Select[Nest[image,{z0},3],(N[Abs[#]]<=bound)&]; FixedPoint[reducedImage, pointsSoFar]; ListPlot[{Re[#],Im[#]}& /@ pointsSoFar, FilterOptions[ListPlot,opts], FilterOptions[ListPlot,Sequence @@ Options[JuliaModified]] ], Message[Julia::badFormula, Julia]], Message[Julia::noInverses] ] ] /; (NumberQ[N[formula /. var -> Pi]] && PolynomialQ[Numerator[Together[formula]], var] && PolynomialQ[Denominator[Together[formula]], var]) || Message[Julia::badFormula, Julia] (* Special cases for Julia *) Julia[] := Message[JuliaSet::argu, "Julia", "0 arguments", "2 arguments are"] Julia[c_] := Message[JuliaSet::argu, "Julia", "1 argument", "2 arguments are"] JuliaStochastic[c_, numpoints_, opts___] := Module[ {inv1,inv2,points,z0,oneBias,chooser, valid = First /@ Options[JuliaStochastic]}, Scan[If[!MemberQ[valid, First[#]], Message[JuliaSet::optx, ToString[First[#]]]]&, {opts} ]; oneBias = OneBias /. {opts} /. Options[JuliaStochastic]; If[!(0 0 ) || Message[JuliaSet::badInt, JuliaStochastic, numpoints] ) && ( And @@ Map[OptionQ, {opts}] || Message[JuliaSet::badOpt,JuliaStochastic] ) (* Special cases for JuliaStochastic *) JuliaStochastic[] := Message[JuliaSet::argu, "JuliaStochastic", "0 arguments", "2 arguments are"] JuliaStochastic[c_] := Message[JuliaSet::argu, "JuliaStochastic", "1 argument", "2 arguments are"] End[] (* End Private Context *) Protect[JuliaSet, Julia, JuliaModified, JuliaSimple, JuliaStochastic, Bound, Resolution, OneBias] EndPackage[]