inclusions_elliptiques.py 1.89 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 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