#-----------------------#
#                       #
#       Aufgabe 3       #
#                       #
#-----------------------#



calctime = Sys.time()             # let's monitor calculation time

maxiter = 100
c = 0 + 1i*0
plot(c,xlim=c(-2,1),ylim=c(-1.5,1.5))
for( x in seq(from=-2,to=1,by=0.005) )
{
  for(y in seq(from=-1.5,to=1.5,by=0.005) )
  {
    c = x+1i*y
    z = 0
    iter = 0
    while( abs(z)<1000 & iter < maxiter )
    {
      z = z^2 +c 
      iter = iter + 1
    }
    if(iter>= maxiter)
      points(c,pch=".",col="black")
  }
}

calctime = Sys.time() - calctime
calctime




# Damit ist die Aufgabe fertig. 
# Jetzt noch ein paar Bilder in Farbe:


# the whole set:

maxiter = 60
mycolor = rainbow(0.75*maxiter)
c = 0 + 1i*0
plot(c,xlim=c(-2,1),ylim=c(-1.5,1.5))
for( x in seq(from=-2,to=1,by=0.005) )
{
  for(y in seq(from=-1.5,to=1.5,by=0.005) )
  {
    c = x+1i*y
    z = 0
    iter = 0
    while( abs(z)<2 & iter < maxiter )      # falls abs(z)>2 folgt 
    {                                       # automatisch Divergenz, 
      z = z^2 +c                            # kann man beweisen.
      iter = iter + 1
    }
    points(c,pch=".",col=mycolor[iter])
  }
}


# zoom in at some point [c - Delta , c + Delta]
# npix = 600 per axis

c = -0.8 + 1i*0.2
Delta = 0.1
npix = 600
delta = (2*Delta)/npix
xmin = Re(c)-Delta
xmax = Re(c)+Delta
ymin = Im(c)-Delta
ymax = Im(c)+Delta

maxiter = 150
mycolor = rainbow(0.75*maxiter)

plot(c,xlim=c(xmin,xmax),ylim=c(ymin,ymax))
for( x in seq(from=xmin,to=xmax,by=delta) )
{
  for(y in seq(from=ymin,to=ymax,by=delta) )
  {
    c = x+1i*y
    z = 0
    iter = 0
    while( abs(z)<2 & iter < maxiter )
    {
      z = z^2 +c 
      iter = iter + 1
    }
    points(c,pch=".",col=mycolor[iter])
  }
}


# further zoomed in:

c = -0.73 + 1i*0.28
Delta = 0.01
npix = 600
delta = (2*Delta)/npix
xmin = Re(c)-Delta
xmax = Re(c)+Delta
ymin = Im(c)-Delta
ymax = Im(c)+Delta

maxiter = 150
mycolor = rainbow(0.75*maxiter)

plot(c,xlim=c(xmin,xmax),ylim=c(ymin,ymax))
for( x in seq(from=xmin,to=xmax,by=delta) )
{
  for(y in seq(from=ymin,to=ymax,by=delta) )
  {
    c = x+1i*y
    z = 0
    iter = 0
    while( abs(z)<2 & iter < maxiter )
    {
      z = z^2 +c 
      iter = iter + 1
    }
    points(c,pch=".",col=mycolor[iter])
  }
}


# further zoomed in:
# das dauert etwas:

c = -0.597214+1i*0.424206 
Delta = 0.0001
npix = 500
#delta = (2*Delta)/npix
delta = (Delta)/npix
xmin = Re(c)-Delta
xmax = Re(c)+Delta
ymin = Im(c)-Delta
ymax = Im(c)+Delta

maxiter = 500
mycolor = rainbow(0.75*maxiter)

plot(c,xlim=c(xmin,xmax),ylim=c(ymin,ymax))
for( x in seq(from=xmin,to=xmax,by=delta) )
{
  for(y in seq(from=ymin,to=ymax,by=delta) )
  {
    c = x+1i*y
    z = 0
    iter = 0
    while( abs(z)<2 & iter < maxiter )
    {
      z = z^2 +c 
      iter = iter + 1
    }
    points(c,pch=".",col=mycolor[iter])
  }
}



# Jetzt: Anstatt z^2+c -> z^p+c

p = 3
maxiter = 100
mycolor = rainbow(0.75*maxiter)
c = 0 + 1i*0
info = paste("z_{n+1} = z_n^p + c  mit  p =",p)
plot(c,xlim=c(-1.5,1.5),ylim=c(-1.5,1.5), main=info)
for( x in seq(from=-1.5,to=1.5,by=0.005) )
{
  for(y in seq(from=-1.5,to=1.5,by=0.005) )
  {
    c = x+1i*y
    z = 0
    iter = 0
    while( abs(z)<1000 & iter < maxiter )
    {
      z = z^p +c 
      iter = iter + 1
    }
    points(c,pch=".",col=mycolor[iter])
  }
}

p = 5
maxiter = 100
mycolor = rainbow(0.75*maxiter)
c = 0 + 1i*0
info = paste("z_{n+1} = z_n^p + c  mit  p =",p)
plot(c,xlim=c(-1.5,1.5),ylim=c(-1.5,1.5), main=info)
for( x in seq(from=-1.5,to=1.5,by=0.0025) )
{
  for(y in seq(from=-1.5,to=1.5,by=0.0025) )
  {
    c = x+1i*y
    z = 0
    iter = 0
    while( abs(z)<1000 & iter < maxiter )
    {
      z = z^p +c 
      iter = iter + 1
    }
    points(c,pch=".",col=mycolor[iter])
  }
}