import processing.core.*; 
import ORG.netlib.math.complex.*; 

public class MobiusDFS extends PApplet {

/**
	 * 
	 */
	private static final long serialVersionUID = 1L;
	
int windowsize=300;
int DEPTH=120;
float EPSILON=0.1f;
float SCALE=2.5f;
matrix[] word = new matrix[1000];
matrix[] gen = new matrix[5];
Complex[][] parabolics = new Complex[5][4];
int level=1;
int[] tag= new int[1000];
Complex oldPoint, newPoint, tempPoint;


public void setup()
{
  size(windowsize, windowsize);
  fill(255);
  rect(0,0,windowsize,windowsize);
  background(255); 
  Complex ta = new Complex(1.95f, 0.02f);
  Complex tb = new Complex(2.01f,0.0f);
  Complex tab = generateTab(ta, tb);
  Complex z0 = generatez0(ta, tb, tab);
  matrix Ma = generateMatrixA(ta, tb, tab, z0);
  matrix Mb = generateMatrixB(ta, tb, tab, z0);
  matrix MA = Ma.inv();
  matrix MB = Mb.inv();
  
  gen[0]=null;
  gen[1] = Ma;
  gen[2] = Mb;
  gen[3] = MA;
  gen[4] = MB;
  
  oldPoint = fixedPoint(gen[1].mul(gen[2]).mul(gen[3]).mul(gen[4]));  
  
  parabolics[1][1] = fixedPoint(Mb.mul(MA).mul(MB).mul(Ma));
  parabolics[1][2] = fixedPoint(Ma);
  parabolics[1][3] = fixedPoint(MB.mul(MA).mul(Mb).mul(Ma));
  parabolics[2][1] = fixedPoint(MA.mul(MB).mul(Ma).mul(Mb));
  parabolics[2][2] = fixedPoint(Mb);
  parabolics[2][3] = fixedPoint(Ma.mul(MB).mul(MA).mul(Mb));
  parabolics[3][1] = fixedPoint(MB.mul(Ma).mul(Mb).mul(MA));
  parabolics[3][2] = fixedPoint(MA);
  parabolics[3][3] = fixedPoint(Mb.mul(Ma).mul(MB).mul(MA));
  parabolics[4][1] = fixedPoint(Ma.mul(Mb).mul(MA).mul(MB));
  parabolics[4][2] = fixedPoint(MB);
  parabolics[4][3] = fixedPoint(MA.mul(Mb).mul(Ma).mul(MB));
  
  tag[1]=1; 
  word[1]=gen[1];
  
  frameRate(10000);
  //noFill();
  stroke(0,0,180);
  strokeWeight(0.2f);
  smooth();

}

public void draw()
{
  while(!(branchTermination()))
  {
    goForward();
  }
  goBackward();
  while( !((level==0)||(availableTurn())))
  {  
    goBackward();
  }  
  if(!((tag[1]==2)&&(level==0)))
  {
    turnAndGoForward();  
  }  
  else
  {
	// saveFrame("DFStest.tif");  
    noLoop();
  }
}

public void goBackward()
{
  level--;
}

public void goForward()
{
  level++;
  tag[level] = (tag[level-1])%4 +1;
  word[level] = word[level-1].mul(gen[tag[level]]);    
}

public boolean availableTurn()
{
  if ((tag[level+1])%4 == (tag[level]+3)%4)
  {
    return false;
  }
  else
  {
    return true;
  }
}

public void turnAndGoForward()
{
  tag[level+1] = (tag[level+1] +3)%4;
  if(tag[level+1]==0)
  {
    tag[level+1]=4;
  }
  if(level==0)
  {
    word[1] = gen[tag[1]];
  }
  else
  {
    word[level+1] = word[level].mul(gen[tag[level+1]]);
  }  
  level++;  
}

public boolean branchTermination()
{ 
  Complex[] z= new Complex[4];
  newPoint = mobius(oldPoint, word[level]);
  float dist=distance(oldPoint, newPoint); 
  
  for(int j=1; j<=3; j++)
  {  
	  z[j]=mobius(parabolics[tag[level]][j],word[level]);
  }
 
  if(level>=DEPTH || ((dist<EPSILON) && (distance(z[1],z[2])<EPSILON) && (distance(z[2],z[3])<EPSILON)))
  { 
	  stroke(0,4*level,250-5*level);
	  drawLine(oldPoint, newPoint);
	  oldPoint=newPoint;
      return true;
  }

  else
  {
    return false;
  }
}

public void drawLine(Complex z1, Complex z2)
{
  line((float) (windowsize/2+SCALE*z1.re()), (float) (windowsize/2+SCALE*z1.im()), (float) (windowsize/2+SCALE*z2.re()), (float) (windowsize/2+SCALE*z2.im()));
}

public float distance(Complex z1, Complex z2)
{
  return sqrt((float) ((z1.re()-z2.re())*(z1.re()-z2.re())+(z1.im()-z2.im())*(z1.im()-z2.im())));
}

public Complex mobius(Complex z, matrix M)
{
  Complex numerator = (M.a).mul(z).add(M.b);
  Complex denominator = (M.c).mul(z).add(M.d);
  return numerator.div(denominator);
} 

class matrix
{
  Complex a,b,c,d;
  matrix(Complex _a, Complex _b, Complex _c, Complex _d)
  {
    a=_a;
    b=_b;
    c=_c;
    d=_d;
  }   
  public Complex det()
  {
    return ((a).mul(d)).sub((b).mul(c));
  } 
  public matrix inv()
  {
    Complex det= det();
    matrix Minv=new matrix((d).div(det), (new Complex(-1,0).mul(b)).div(det), (new Complex(-1,0).mul(c)).div(det), a.div(det));
    return Minv;
  }  
  public matrix mul(matrix N)
  {
     Complex p1 = ((a).mul(N.a)).add((b).mul(N.c));
     Complex p2 = ((a).mul(N.b)).add((b).mul(N.d)); 
     Complex p3 = ((c).mul(N.a)).add((d).mul(N.c)); 
     Complex p4 = ((c).mul(N.b)).add((d).mul(N.d));
     return new matrix(p1, p2, p3, p4);
  }  
} 
  
public Complex fixedPoint(matrix M)
{
  Complex fixedPoint1;	
  Complex disc = ((M.a).sub(M.d)).mul((M.a).sub(M.d)).add(new Complex(4,0).mul(M.b).mul(M.c));  
  fixedPoint1 = ((M.a).sub(M.d).add(disc.sqrt())).div(new Complex(2,0).mul(M.c));
  // fixedPoint2 = ((M.a).sub(M.d).sub(disc.sqrt())).div(new Complex(2,0).mul(M.c));
  return fixedPoint1;
}
    
public Complex generateTab(Complex ta, Complex tb)
{
  return quadraticFormula(new Complex(1,0), ta.mul(tb).mul(new Complex(-1,0)), (ta.mul(ta)).add(tb.mul(tb)));
}

public Complex generatez0(Complex ta, Complex tb, Complex tab)
{
  Complex numerator = (tab.sub(new Complex(2,0))).mul(tb);
  Complex denominator = (tb.mul(tab)).add((new Complex(-2,0)).mul(ta)).add((new Complex(0,2)).mul(tab));
  return numerator.div(denominator);
} 

public Complex quadraticFormula(Complex a, Complex b, Complex c)
{
  Complex discriminant;
  discriminant = (b.mul(b)).sub(new Complex(4,0).mul(a).mul(c));
  return ((new Complex(-1,0).mul(b)).add(discriminant.sqrt())).div(new Complex(2,0).mul(a));
} 

public matrix generateMatrixA(Complex ta, Complex tb, Complex tab, Complex z0)
{
  
  Complex a = ta.div(new Complex(2,0));
  Complex b = ((ta.mul(tab)).add((new Complex(-2,0)).mul(tb)).add(new Complex(0,4))).div((((new Complex(2,0)).mul(tab)).add(new Complex(4,0))).mul(z0));
  Complex c = (ta.mul(tab).add((new Complex(-2,0)).mul(tb)).add(new Complex(0,-4))).mul(z0).div(new Complex(2,0).mul(tab).sub(new Complex(4,0)));  
  Complex d = a;
  return new matrix(a,b,c,d);
}

public matrix generateMatrixB(Complex ta, Complex tb, Complex tab, Complex z0)
{
  Complex a = (tb.sub(new Complex(0,2))).div(new Complex(2,0));
  Complex b = tb.div(new Complex(2,0));
  Complex c = b;
  Complex d = (tb.add(new Complex(0,2))).div(new Complex(2,0));
  return new matrix(a,b,c,d);
}  
  
static public void main(String args[]) {   PApplet.main(new String[] { "noCircleDFS2" });}}