diff --git a/inclusions_elliptiques.py b/inclusions_elliptiques.py
new file mode 100644
index 0000000000000000000000000000000000000000..8f7d671fc155d3549d4c89eba53f236c594468c7
--- /dev/null
+++ b/inclusions_elliptiques.py
@@ -0,0 +1,93 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+from random import uniform
+
+
+def compteur(im):
+    nn,mm=im.shape
+    cpt=0
+    for i in range(nn):
+        for j in range(mm):
+            if im[i,j]==1:
+                cpt=cpt+1
+    return cpt/nn/mm
+    
+    
+def ajoute(im,a,b,th,iso,nb):
+    nn,mm=im.shape
+    centers=[[0.,0.] for i in range(nb)]
+    theta=[0. for i in range(nb)]
+    
+    if iso==True:
+        for i in range(nb):
+            theta[i]=uniform(0,np.pi)
+    else:
+        for i in range(nb):
+            theta[i]=th
+        
+    for i in range(nb):
+        centers[i]=[uniform(0,1),uniform(0,1)]
+    
+    
+
+    for i in range(nb):
+        c=max(a,b)
+        x=centers[i][0]
+        y=centers[i][1]  
+        the=theta[i]  
+        for j in range(int((x-c)*nn)-1,int((x+c)*nn)+1):
+            for k in range(int((y-c)*mm)-1,int((y+c)*mm)+1):                
+                x1=(j+0.5)/nn
+                y1=(k+0.5)/mm
+                x2=x1*np.cos(the)+y1*np.sin(the)
+                y2=y1*np.cos(the)-x1*np.sin(the)
+                if (x2-x)**2/a**2+(y2-y)**2/b**2<=1:
+                    if j>=nn:
+                        l=nn
+                    else:
+                        l=0
+                    if k>=mm:
+                        m=mm
+                    else:
+                        m=0
+                    im[j-l,k-m]=1
+    
+    
+
+def image_init(N,a,b,th,iso,f):
+
+    im=np.zeros((N,N))
+    
+    n=int(f/np.pi/a/b)
+    
+    c=0
+    while c<f:
+        ajoute(im,a,b,th,iso,n)
+        c=compteur(im)
+        n=int((f-c)/np.pi/a/b)+1
+
+    return im,c
+    
+
+def image_augmente(im,a,b,th,iso,f):
+    f0=compteur(im)
+    if f<=f0:
+        print('fraction volumique actuelle :'+str(f0))
+    else:
+        c=f0
+        n=int((f-c)/np.pi/a/b)+1
+        while c<f:
+            ajoute(im,a,b,th,iso,n)
+            c=compteur(im)
+            n=int((f-c)/np.pi/a/b)+1
+        
+        return c
+            
+
+
+    
+
+
+