(* File: LSystem.m *)
(* Date: /12/08/97 *)
(* Author: J.-P. Kuska,
   email: jpk@osf2.mpae.gwdg.de
 *)

(* Copyright 1996 Jens-Peer Kuska
   this file may be freely distributet als long as the copyright notice
   is left intact
   *)

(* Version: 1.1             *)
(* Mathematica Version: 3.0 *)

(*
Sources:
  ``The Science of Fractal Images''
  Heinz-Otto Peitgen and Dietmar Saupe, Ed.
  Springer-Verlag 1988
  page 273, Appendix C
  
  fractint documentation
 *)

(* Requirements: None *)

(* Limitations: None *)

BeginPackage["LSystem`","Graphics`Colors`"]

LSystem::usage=
"LSystem[what,stage] will produce the 2D lines of the stage times iteratet
 L-sytem due to the parameters LSystemN[what], LSystemAxiom[what] and
 LSystemRules[what] of the database."

LSystemStringToTurtle::usage=
"LSystemStringToTurtle[lsys,divsions] puts the result of the
iteration lsys and will convert the lsys string to moves of the
turtle. The meaning of the characters is:\n
F      inc position and draw line\n
f      inc position and move\n
+      turn right by 2Pi/divisions\n
-      turn left by 2Pi/divisions\n
|      turn back\n
[      push the turtle state to stack\n
]      pop the turtle state from stack\n
@pnum  scale the width of the turtle step by the positive number pnum,\n \
       @Qnum means square root num @Inum means scale by 1/num\n
!      flip menaing of + and -\n
cnum   set color number num\n
<num   decrement color by num\n
>num   increment color by num"

LSystemRules::usage=
"LSystemRules[fractalname] gives the transformation rules applied to
 the string on every stage."

LSystemAxiom::usage=
"LSystemAxiom[fractalname] gives the initial state of the l-system for
the fractal fractal name."

LSystemN::usage=
"LSystemN[fractalname] gives the angle delta=2Pi/LSystem[fractalname] for
the given fractal."


(#::usage=ToString[#] <> "is a data base entry for LSystem["<>ToString[#]<>" ,int].") & /@
  {Sierpinsky, SierpinskySquare, HilbertCurve, KochSnowFlake, KochIsland, DragonCurve,
   Bush1, Bush2, Penrose1, Penrose2, Penrose3,
   Island1,Island2,Quartet,Fass1,Fass2,Peano2,
   SnowFlake1,SnowFlake2,SnowFlake3   
} 

Begin["`Private`"]

$LSDigit={"0","1","2","3","4","5","6","7","8","9",".","I","Q"};
$LSNumberFollows={"@","\\","/","c","<",">"};
$LSOrdinary={"F","f","+","-","]","[","|","!"};

$LSValidChar=Join[$LSOrdinary,$LSNumberFollows,$LSDigit]

$LSColors=Sort[ToExpression /@ AllColors];
$LSColorNo=Length[$LSColors]

Clear[LSTokenScan]


LSToNumber[a_String]:=
  Module[{val,lst=Characters[a]},

    Switch[First[lst],
      "I", val=LSToNumber[StringJoin @@ Rest[lst]]; val=1.0/val,
      "Q", val=LSToNumber[StringJoin @@ Rest[lst]]; val=N[Sqrt[val]],
      _, val=ToExpression[a]
      ]; 
     val
    ]  
				
LSTokenScan[{numfollow_  /;MemberQ[$LSNumberFollows,numfollow],an__}]:=
  Module[{nstr="",alst={an}},
    While[alst!={} &&MemberQ[$LSDigit,First[alst]],
      nstr=nstr<>First[alst];
      alst=Rest[alst]
     ];
    alst=LSTokenScan[alst]; 
    {{numfollow,LSToNumber[nstr]},Sequence @@ alst}
   ]

LSNumToken[{numfollow_  /;MemberQ[$LSNumberFollows,numfollow],an__}]:=
  Module[{nstr="",alst={an}},
    While[alst!={} &&MemberQ[$LSDigit,First[alst]],
      nstr=nstr<>First[alst];
      alst=Rest[alst]
     ];
    {{numfollow,LSToNumber[nstr]}, alst}
   ]

(* Iteration is used here to avoid to deep recursion and save memory *)

LSTokenScan[l_List]:=
  Module[{top={},toks=l,nftok},
    While[toks=!={},
      If[!MemberQ[$LSNumberFollows,First[toks]], (* Then *)
        AppendTo[top,First[toks]];
        toks=Rest[toks],
      (* Else *)
        {nftok,toks}=LSNumToken[toks];
        AppendTo[top,nftok]
       ]  
     ];
    top
   ]
   
LSTokenScan[{any_}]:={any}
LSTokenScan[{}]:={}


LSValidQ[tok_]:=MemberQ[$LSOrdinary,tok] || (Head[tok]===List && Lenght[tok]==2)

LSToTokens[s_String]:=
  Module[{clst,toklst},
    clst=Select[Characters[s],MemberQ[$LSValidChar,#]&];
    (* 
       Since LSTokenScan[] is recursive to the 
       depth of the Lenght of the list the $RecursionLimit 
       must be increased local. 
    *)
    Block[{$RecursionLimit=Max[21,2*Length[clst]]},
      LSTokenScan[clst]
     ]
   ]

Clear[LSValid]

LSValidQ[tok_String]:=MemberQ[$LSOrdinary,tok]
LSValidQ[{s_String,f_?NumericQ}]:=True
LSValidQ[_]:=False

IncDirection[]:=If[turtleDir+1<nn,turtleDir++,turtleDir=0]
DecDirection[]:=If[turtleDir>0, turtleDir--,turtleDir=nn-1]

SetColor[n_Integer]:=(If[n>=1 && n<=$LSColorNo,
                       AppendTo[lines,{$LSColors[[n]],Line[tmoves]}];colorIndex=n,
                       AppendTo[lines,{Black,Line[tmoves]}];colorIndex=1
                      ];
                      tmoves={{turtleX,turtleY}};)
                      

DecColor[n_]:=(If[colorIndex-n >=1,
                 colorIndex-=n,
                 colorIndex=$LSColorNo-n
                ];
                AppendTo[lines,{$LSColors[[colorIndex]],Line[tmoves]}];
                tmoves={{turtleX,turtleY}};)   

IncColor[n_]:=(If[colorIndex+n <=$LSColorNo,
                 colorIndex+=n,
                 colorIndex=colorIndex+n-$LSColorNo
                ];
                AppendTo[lines,{$LSColors[[colorIndex]],Line[tmoves]}];
                tmoves={{turtleX,turtleY}};)   

$LSystemTurtle=Dispatch[
       {"F" :> (AppendTo[tmoves,{turtleX,turtleY}+=turtleStep*scangle[[turtleDir+1]]];),
        "f" :> (AppendTo[lines,Line[tmoves]];
                tmoves={{turtleX,turtleY}+=turtleStep*scangle[[turtleDir+1]]};),
        "+" :>(If[!revertDir,IncDirection[],DecDirection[]];),
        "[" :>(PrependTo[pstack,{turtleX,turtleY}];
               PrependTo[dstack,{turtleDir,revertDir}];),
        "]" :>(If[pstack=!={},
                  {turtleX,turtleY}=First[pstack];
                   pstack=Rest[pstack];
                   {turtleDir,revertDir}=First[dstack];
                   dstack=Rest[dstack]];
                   AppendTo[lines,Line[tmoves]];
                   tmoves={{turtleX,turtleY}};),
        "-" :>(If[!revertDir,DecDirection[],IncDirection[]];),
        "|" :>(If[turtleDir<Round[nn/2],
                 turtleDir+=Round[nn/2],
                 turtleDir-=Round[nn/2]
                ];),
         {"@",scale_} :> (turtleStep*=scale ;),
         "!" :> (revertDir=!revertDir;) ,
         {"c",cno_Integer} :> (SetColor[cno];),
         {"<",cno_Integer} :> (DecColor[cno];),
         {">",cno_Integer} :> (IncColor[cno];)     

      }];
 
LSystemStringToTurtle[lstring_String,n_Integer]:=
  Block[{turtleX,turtleY,turtleDir,turtleStep,lsystemTurtle,nn=n,
         llst,scangle,i,revertDir=False,colorIndex,pnts,cno,
         lines={},pstack={},dstack={},
         tmoves={{0.0,0.0}}},  
     llst=LSToTokens[lstring];       
     llst=Select[llst,LSValidQ];
     colorIndex=1;
     turtleX=0.0;
     turtleY=0.0;
     turtleDir=0;
     turtleStep=1.0;
     scangle=Table[
               N[{Cos[i*2Pi/n],Sin[i*2Pi/n]}],
               {i,0,n-1}
              ];      
     Scan[# /. $LSystemTurtle &,llst];
     Append[lines,Line[tmoves]] /. {Line[{{_,_}}]:>Sequence[],Line[{}] :> Sequence[]}
   ]
   
FoldStringReplace[a_String,rules_]:=Fold[StringReplace,a,rules]   


LSystem[axiom_String,repRules_,angleN_Integer,iter_Integer]:=
  LSystemStringToTurtle[
   Nest[
      StringReplace[
        #,
        repRules] &,
      axiom,
      iter
     ] 
     ,
   angleN
  ]

LSystem[sym_Symbol,Iter_Integer]:=LSystem[
                                    LSystemAxiom[sym],
                                    LSystemRules[sym],
                                    LSystemN[sym],
                                    Iter]
 

(* the data base *)

Sierpinsky /:
  LSystemRules[Sierpinsky]:={"X"->"--FXF++FXF++FXF--","F"->"FF"}
Sierpinsky /:
  LSystemAxiom[Sierpinsky]:="FXF--FF--FF"
Sierpinsky /:
  LSystemN[Sierpinsky]:=6

HilbertCurve /:
  LSystemRules[HilbertCurve]:={"X"->"-YF+XFX+FY-", "Y"->"+XF-YFY-FX+"}
HilbertCurve /:
  LSystemAxiom[HilbertCurve]:="X"
HilbertCurve /:
  LSystemN[HilbertCurve]:=4

KochSnowFlake /:
  LSystemRules[KochSnowFlake]:={"F"->"F-F++F-F"}
KochSnowFlake /:
  LSystemAxiom[KochSnowFlake]:="F"
KochSnowFlake /:
  LSystemN[KochSnowFlake]:=6

KochIsland /:
  LSystemRules[KochIsland]:={"F"->"F+F-F-FF+F+F-F"}
KochIsland /:
  LSystemAxiom[KochIsland]:="F+F+F+F"
KochIsland /:
  LSystemN[KochIsland]:=4

DragonCurve /:
  LSystemRules[DragonCurve]:={"X"->"X+YF+","Y"->"-FX-Y"}
DragonCurve /:
  LSystemAxiom[DragonCurve]:="X"
DragonCurve /:
  LSystemN[DragonCurve]:=4
SierpinskySquare /:
  LSystemRules[SierpinskySquare]:={"F"->"FF+F+F+F+FF"}
SierpinskySquare /:
  LSystemAxiom[SierpinskySquare]:="F+F+F+F"
SierpinskySquare /:
  LSystemN[SierpinskySquare]:=4

Bush1 /:
  LSystemRules[Bush1]:={"F"->"F[+F]F[-F]F"}
Bush1 /:
  LSystemAxiom[Bush1]:="F"
Bush1 /:
  LSystemN[Bush1]:=14
Bush2 /:
  LSystemRules[Bush2]:={"G"->"GFX[+G][-G]",
                        "X"->"X[-FFF][+FFF]FX"}
Bush2 /:
  LSystemAxiom[Bush2]:="G"
Bush2 /:
  LSystemN[Bush2]:=14


Island1 /: LSystemN[Island1]:= 4
Island1 /:LSystemAxiom[Island1]:="F+F+F+F"
Island1 /:LSystemRules[Island1]:={"F"->"FFFF-F+F+F-F[-fFF+F+FF+F]FF",
                                  "f"->"@8f@I8"}


Island2 /: LSystemN[Island2]:= 4
Island2 /:LSystemAxiom[Island2]:="F+F+F+F"
Island2 /:LSystemRules[Island2]:={"F"->"F+fF-FF-F-FF+f+FF-fF+FF+F+FF-f-FFF" ,
                                  "f"->"@6f@I6"}
LSystemN[Quartet]^:=4
LSystemAxiom[Quartet]^:="FB"
LSystemRules[Quartet]:={
    "A"->"FBFA+HFA+FB-FA",
    "B"->"FB+FA-FB-JFBFA",
    "F"->"",
    "H"->"-",
    "J"->"+"}

LSystemN[Fass1]^:=4
LSystemAxiom[Fass1]^:="-L"
LSystemRules[Fass1]^:={"L"->"LF+RFR+FL-F-LFLFL-FRFR+",
		       "R"->"-LFLF+RFRFR+F+RF-LFL-FR"}    


LSystemN[Fass2]^:=4
LSystemAxiom[Fass2]^:="-L"
LSystemRules[Fass2]^:={"L"->"LFLF+RFR+FLFL-FRF-LFL-FR+F+RF-LFL-FRFRFR+",
		       "R"->"-LFLFLF+RFR+FL-F-LF+RFR+FLF+RFRF-LFL-FRFR"}


LSystemN[Peano2]^:=8
LSystemAxiom[Peano2]^:="FXY++F++FXY++F"
LSystemRules[Peano2]^:={"X"->"XY@Q2-F@IQ2-FXY++F++FXY",
                        "Y"->"-@Q2F-@IQ2FXY"}


LSystemN[SnowFlake1]^:=12
LSystemAxiom[SnowFlake1]^:="FR"
LSystemRules[SnowFlake1]^:={"R"->"++!FRFU++FU++FU!---@Q3FU|-@IQ3!FRFU!",
		           "U"->"!FRFU!|+@Q3FR@IQ3+++!FR--FR--FRFU!--",
		           "F"->""}
LSystemN[SnowFlake2]^:=12
LSystemAxiom[SnowFlake2]^:="F"
LSystemRules[SnowFlake2]^:={ "F"->"++!F!F--F--F@IQ3|+F!F--F--F!+++@Q3F@QI3|+F!F@Q3|+F!F"}
    
LSystemN[SnowFlake3]^:=12
LSystemAxiom[SnowFlake3]^:="FX"
LSystemRules[SnowFlake3]^:={
    "X"->"++F!X!FY--FX--FY|+@IQ3FYF!X!++F!Y!++F!Y!FX@Q3+++F!Y!FX",
    "Y"->"FYF!X!+++@IQ3FYF!X!++F!X!++F!Y!FX@Q3|+FX--FY--FXF!Y!++",
    "F"->""}    
    
Penrose1 /: LSystemN[Penrose1]:=10
Penrose1 /: LSystemAxiom[Penrose1]:="+WF--XF---YF--ZF"
Penrose1 /: LSystemRules[Penrose1]:={ "W"->"YF++ZF----XF[-YF----WF]++",
                                      "X"->"+YF--ZF[---WF--XF]+",
                                      "Y"->"-WF++XF[+++YF++ZF]-",
                                      "Z"->"--YF++++WF[+ZF++++XF]--XF",
                                      "F"->""}

Penrose2 /: LSystemN[Penrose2]:=10
Penrose2 /: LSystemAxiom[Penrose2]:="++ZF----XF-YF----WF"
Penrose2 /: LSystemRules[Penrose2]:={"W"->"YF++ZF----XF[-YF----WF]++",
                                     "X"->"+YF--ZF[---WF--XF]+",
                                     "Y"->"-WF++XF[+++YF++ZF]-",
                                     "Z"->"--YF++++WF[+ZF++++XF]--XF",
                                     "F"->""}

Penrose3 /: LSystemN[Penrose3]:=10
Penrose3 /: LSystemAxiom[Penrose3]:="[X]++[X]++[X]++[X]++[X]"
Penrose3 /: LSystemRules[Penrose3]:={"W"->"YF++ZF----XF[-YF----WF]++",
                                     "X"->"+YF--ZF[---WF--XF]+","Y"->"-WF++XF[+++YF++ZF]-",
                                     "Z"->"--YF++++WF[+ZF++++XF]--XF","F"->""}

End[]
EndPackage[]
Null


(* Examples:

Show[Graphics[LSystem[Bush1,4]],AspectRatio->Automatic];

Show[Graphics[LSystem[Bush2,6]],AspectRatio->Automatic];

Show[Graphics[LSystem[Sierpinsky,5]],AspectRatio->Automatic];

Show[Graphics[LSystem[KochIsland,2]],AspectRatio->Automatic];

Show[Graphics[LSystem[SierpinskySquare,4]],AspectRatio->Automatic];

Show[Graphics[LSystem[DragonCurve,12]],AspectRatio->Automatic];
*)