import matplotlib.pyplot as plt
import random as rd
import numpy as np


## Exercice 1

xs=[-1,0,1]
ys=[0,1.732,0]

plt.scatter(xs,ys)

a,b,c=rd.random(),rd.random(),rd.random()

x0 = (a*xs[0]+b*xs[1]+c*xs[2])/(a+b+c)
y0 = (a*ys[0]+b*ys[1]+c*ys[2])/(a+b+c)

plt.scatter(x0,y0)

n=10000
x=x0
y=y0
X,Y=[],[]

for k in range(n):
    i=rd.randint(0,2)
    x=(x+xs[i])/2
    y=(y+ys[i])/2
    X.append(x)
    Y.append(y)

plt.scatter(X,Y,s=1)
    # pause(1e-10)

plt.axis("equal")
plt.show()

##

s=7 # nombre de sommets
f=3 # facteur d'échelle
n=5000 # nombre de points

# sommets
xs=[np.cos(2*k*np.pi/s+np.pi/2) for k in range(s)]
ys=[np.sin(2*k*np.pi/s+np.pi/2) for k in range(s)]

plt.scatter(xs,ys)

# point initial
L=[rd.random() for k in range(s)]
p=sum(L)
x0,y0=0,0
for k in range(s):
    x0=x0+L[k]*xs[k]/p
    y0=y0+L[k]*ys[k]/p

plt.scatter(x0,y0)
x=x0
y=y0
X,Y=[],[]


# récurrence
for k in range(n):
    i=rd.randint(0,s-1)
    x=(x+f*xs[i])/(1+f)
    y=(y+f*ys[i])/(1+f)
    X.append(x)
    Y.append(y)

plt.scatter(X,Y,s=0.1,color='r')

plt.axis("equal")
plt.show()

## Exercice 2

x,y=0,0
n=1000
X,Y=[],[]
for k in range(n):
    r=rd.random()
    if r<0.01:
        x=0
        y=0.16*y
    elif r<0.86:
        x,y=0.85*x+0.04*y,-0.04*x+0.85*y+1.6
    elif r<0.93:
        x,y=0.2*x-0.26*y,0.23*x+0.22*y+1.6
    else:
        x,y=-0.15*x+0.28*y,0.26*x+0.24*y+0.44
    X.append(x)
    Y.append(y)

plt.scatter(X,Y,s=0.1,color='g')
plt.axis("equal")
plt.show()

## Exercice 3 - Julia

# c=0.3+0.5*1j
# c=0.285+0.013*1j
# c=-0.8+0.156*1j
# c=-0.4+0.6*1j
c=-0.123+0.745*1j
N=50

from numba import jit
@jit(nopython=True)

def julia(x,y):
    z=x+y*1j
    n=0
    while abs(z)<2 and n<N:
        z=z**2+c
        n+=1
    return(n==N)

p=400
X=np.linspace(-2,2,p)
Y=np.linspace(-2,2,p)
Jx,Jy=[],[]

for x in X:
    for y in Y:
        if julia(x,y):
            Jx.append(x)
            Jy.append(y)

plt.scatter(Jx,Jy,s=1)
plt.axis("equal")
plt.show()

## Exercice 3 - Mandelbrot

N=50

from numba import jit
@jit(nopython=True)

def mandelbrot(x,y):
    z=0
    c=x+y*1j
    n=0
    while abs(z)<2 and n<N:
        z=z**2+c
        n+=1
    return(n==N)

p=100
X=np.linspace(-2,1,3*p)
Y=np.linspace(0,1,p)
Mx,My=[],[]

for x in X:
    for y in Y:
        if mandelbrot(x,y):
            Mx.append(x)
            My.append(y)
            if y!=0:
                Mx.append(x)
                My.append(-y)

plt.scatter(Mx,My,s=1)
plt.axis("equal")
plt.show()

#avec des couleurs :

# def mandelbrot(x,y):
#     z=0
#     c=x+y*1j
#     n=0
#     while abs(z)<2 and n<N:
#         z=z**2+c
#         n+=1
#     return(n)
#
# p=200
# X=np.linspace(-2,1,3*p)
# Y=np.linspace(0,1,p)
# Mx,My,C=[],[],[]
#
# for x in X:
#     for y in Y:
#         m=mandelbrot(x,y)
#         C.append(m/256)
#         Mx.append(x)
#         My.append(y)
#         if y!=0:
#             Mx.append(x)
#             My.append(-y)
#             C.append(m/256)
#
# plt.scatter(Mx,My,s=1,c=C)

## Exercice 4


import turtle as t
import math as m


def cote(l, a, n):
    x,y=t.pos()
    if n==0:
        t.goto(x+l*m.cos(a), y+l*m.sin(a))
    else:
        cote(l/3,a,n-1)
        cote(l/3,a+m.pi/3,n-1)
        cote(l/3,a-m.pi/3,n-1)
        cote(l/3,a,n-1)


def flocon(n):
    try:
        t.reset()          #réinitialisation de la figure
    except t.Terminator:
        pass
    cote(200,m.pi/3,n)   #tracé des cotés du triangle
    cote(200,-m.pi/3,n)
    cote(200,-m.pi,n)
    t.exitonclick()

for n in range(5):
    flocon(n)


## Exercice 5

def newton(f,df,u0,eps=2**(-51)):
    h=lambda x:x-f(x)/df(x)
    u=u0
    v=h(u)
    n=0
    while abs(u-v)>eps:
        n+=1
        u=v
        v=h(v)
    return(u,n)

f1=lambda x:x**2-2
df1=lambda x:2*x

print(newton(f1,df1,3,2**(-51)))
# print(newton(f2,df2,3,0.001))
# print(newton(f3,df3,0.5,0.001))

def bassin(z):
    f=lambda x:x**3-1
    df=lambda x:3*x**2
    l=newton(f,df,z)[0]
    c=l.imag
    if abs(c)<0.01:
        return 0
    elif abs(c-np.sqrt(3)/2)<0.01:
        return 1
    else:
        return 2

print(bassin(complex(2,4)))

a,b,c,d=-1.5,1.5,-1.5,1.5
N=200

X=np.linspace(a,b,N)
Y=np.linspace(c,d,N)

M=np.zeros([N,N])

for i in range(N):
    for j in range(N):
        M[i,j]=bassin(complex(X[i],Y[j]))
plt.imshow(M,interpolation="nearest")
plt.show()

##

def bassin5(z):
    f=lambda x:x**5-1
    df=lambda x:5*x**4
    l=newton(f,df,z)[0]
    c=l.imag
    if abs(c)<0.01:
        return 0
    elif abs(c-np.sin(2*np.pi/5))<0.01:
        return 1
    elif abs(c-np.sin(4*np.pi/5))<0.01:
        return 2
    elif abs(c-np.sin(6*np.pi/5))<0.01:
        return 3
    else:
        return 4

def bassin6(z):
    f=lambda x:x**6-1
    df=lambda x:6*x**5
    l=newton(f,df,z)[0]
    c=l.imag
    d=l.real
    if abs(c)<0.01:
        if abs(d-1)<0.01:
            return 0
        else:
            return 1
    elif abs(c-np.sqrt(3)/2)<0.01:
        if abs(d-1/2)<0.01:
            return 2
        else:
            return 3
    else:
        if abs(d-1/2)<0.01:
            return 4
        else:
            return 5

for i in range(N):
    for j in range(N):
        M[i,j]=bassin5(complex(X[i],Y[j]))
plt.imshow(M,interpolation="nearest")
plt.show()