Avis et conseils Simulation phénomène : Python

Signaler
Messages postés
4
Date d'inscription
mardi 17 mars 2020
Statut
Membre
Dernière intervention
18 mars 2020
-
Jm2600
Messages postés
4
Date d'inscription
mardi 17 mars 2020
Statut
Membre
Dernière intervention
18 mars 2020
-
Bonjour,

Je programme depuis peu à programmer sur python. Dans le cadre d'un projet, j'ai réalisé une simulation du trajet des ondes sonores dans l'océan. Ce programme se base notament sur l'influence du gradient de vitesse et approxime le trajet des ondes par une succession d'arc de cercle.
J'ai pour l'instant ce résultat :

Mon but serait d'arriver à ce résultat :


J'ai commencé en parallèle l'implémentation d'un fond aléatoire avec différent sédiements avec un indice de réfraction propre:


Voici le programme principal:
# importation du module
import numpy as np
import matplotlib.pyplot as plt
import array as arr
from math import *
from matplotlib.collections import LineCollection
from matplotlib.colors import ListedColormap, BoundaryNorm
import matplotlib.colors as colors
import random
from mpl_toolkits.axes_grid1 import make_axes_locatable
import pylab

#Initialisation :
profondeur=-200
distanceMax=2000
pas=1
Nbr_courbure=120
Nbr_ondes=1500
Nbr_ondes_simplifiees=500
vitesse_surface=1520
derivation=5*pi/65
z=50
Poiss_x=[]
Poiss_y=[]
Xtot1,Ytot1=[],[]
l,L=20,20
l2,L2=100,1
L=[]
#tableaux
gradient=[0.5,-0.5,-0.05,0.05]
Ord_coude=[0,-50,-120,-250,-1000]


def tranche(a,A):
    n1=len(A)
    Rg_list=-1
    for i in range(n1-1):
        if(A[i]>=a):
            if(A[i+1]<a):
                Rg_list=i
                break
    return(Rg_list)

def Angle(V):
    if(V[0]!=0 and V[1]!=0):
        if(V[0]<=0):
            return(acos(V[1]/sqrt(V[0]**2+V[1]**2)))
        if(V[0]>0):
            return(acos(V[1]/sqrt(V[0]**2+V[1]**2)))
    else:
        return(1)

def Cercle(Rayon,Abscisses,Ordonnees,Agl_init,x_init,y_init,tab_prof,poiss,dist):
    Nbpt=int(abs(Rayon))
    stop,d=0,dist
    t_init=tranche(y_init,tab_prof)
    global L
    if(Rayon>=0):
        th=np.linspace(Agl_init,2*pi+Agl_init,Nbpt)
    elif(Rayon<0):
        th=np.linspace(2*pi+Agl_init,Agl_init,Nbpt)
    x_centre=x_init+Rayon*cos(-Agl_init)
    y_centre=y_init+Rayon*sin(-Agl_init)
    for i in range(Nbpt):
        Abscisses.append(x_centre+Rayon*cos(-pi-th[i]))
        Ordonnees.append(y_centre+Rayon*sin(-pi-th[i]))
        cond,indice=TestPoiss(poiss,Ordonnees,Abscisses)
        if(Abscisses[-1]>distanceMax):
            stop=1
            Abscisses[-1]=distanceMax
            break
        if(i>0):
            if(i==1):
                l=sqrt((Ordonnees[-1]-Ordonnees[-2])**2+(Abscisses[-1]-Abscisses[-2])**2)
                L.append(l)
            d=d+l
        if(cond):
            stop,poiss=3,1
            break
        if(Ordonnees[-1]<profondeur):
            stop = 1
            Ordonnees[-1]=profondeur
            break
        if(Ordonnees[-1]>0):
            stop = 2
            Ordonnees[-1]=0
            break
        if(tranche(Ordonnees[-1],tab_prof)!=t_init):
            break         
    return(Abscisses,Ordonnees,stop,poiss,indice,d)

def TestPoiss(poiss,Ordonnees,Abscisses):
    if(poiss==0):
        n=len(Poiss_y)
        for i in range(n):
            if(Ordonnees[-1]>-60 and Ordonnees[-1]<-40):
                if(Abscisses[-1]>Poiss_x[i]-L and Abscisses[-1]<Poiss_x[i]+L):
                    return(True,i)
    return(False,-1)

def vitesse(list_grad,list_coude,prof_pt):
    for h in range(len(list_coude)-1):
        if(prof_pt<=list_coude[h]):
            if(prof_pt>list_coude[h+1]):
                S=list_grad[0]
                for g in range(h):
                    S=S-list_grad[g]
                tamp=list_grad[h]*S+vitesse_surface
                vitesse_pt=list_grad[h]*prof_pt+tamp
                return(round(vitesse_pt,2))
    return(vitesse_surface)

def ligne(abs_trans,ord_trans,Alpha,Abscisses,Ordonnees):
    for i in range(0,7500):
        Abscisses.append(abs_trans+i)
        Ordonnees.append(ord_trans+i*cos(Alpha)/150)
        if(Ordonnees[-1]<Poiss_y-10):
            break
        if(Ordonnees[-1]>Poiss_y+10):
            break
        if(Abscisses[-1]>Poiss_x+200):
            break
    return(Abscisses,Ordonnees)

def Calcul(X1,Y1,cond):
    if(cond==0):
        angleux=Angle([X1[-1]-X1[-2],Y1[-1]-Y1[-2]])
        vites=vitesse(gradient,Ord_coude,Y1[-1])
        b=gradient[tranche(Y1[-1],Ord_coude)]
        return(angleux,vites,b)
    elif(cond==1):
        vites=vitesse(gradient,Ord_coude,Y1[-1])
        b=gradient[tranche(Y1[-1],Ord_coude)]
        return(vites,b)

def trace(foyer_x,foyer_z):
    X1,Y1=[],[]
    D=[0]*(Nbr_ondes-2*Nbr_ondes_simplifiees)
    for k in range(Nbr_ondes_simplifiees,Nbr_ondes-Nbr_ondes_simplifiees):
        theta0,x0,y0,X,Y=k*pi/Nbr_ondes,foyer_x,foyer_z,[],[]
        poiss,dist=0,D[k-Nbr_ondes_simplifiees]
        a0=gradient[tranche(y0,Ord_coude)]
        if(a0!=0):
            c0=vitesse(gradient,Ord_coude,y0)
            R0=c0/(sin(theta0)*a0)
            X,Y,stop,poiss,indice,dist=Cercle(R0,X,Y,theta0,x0,y0,Ord_coude,poiss,dist)
            while(X[-1]<distanceMax):
                if(len(X)>1 and len(Y)>1):
                    theta,c,a=Calcul(X,Y,0)
                else:
                    c,a=Calcul(X,Y,1)
                    theta=theta0
                if(a!=0):
                    R=c/(sin(theta)*a)
                    X,Y,stop,poiss,indice,dist=Cercle(R,X,Y,theta,X[-1],Y[-1],Ord_coude,poiss,dist)
                    if(stop==2):
                        c,a=Calcul(X,Y,1)
                        if(len(X)>=4):
                            theta1=-Angle([X[-1]-X[-4],Y[-1]-Y[-4]])+pi
                        else:
                            theta1=-Angle([X[-1]-X[-3],Y[-1]-Y[-3]])+pi
                        R=c/(sin(theta1)*a)
                        X,Y,stop,poiss,indice,dist=Cercle(R,X,Y,theta1,X[-1],-0.001,Ord_coude,poiss,dist)
                    if(stop==1):
                        break
                    if(stop==3):
                        theta,c,a=Calcul(X,Y,0)
                        theta=theta+derivation*random.randint(-2,2)
                        R=c/(sin(theta)*a)
                        if(poiss==1):
                            X,Y,stop,poiss,indice,dist=Cercle(R,X,Y,theta,Poiss_x[indice],Y[-1],Ord_coude,poiss,dist)
        X1,Y1=X1+X,Y1+Y
        Xtot1.append(X)
        Ytot1.append(Y)
        D[k-Nbr_ondes_simplifiees]=round(dist,2)
    #print(D)
    return(X1,Y1)

class Poisson():
    def __init__(self,x,y,orientation):
        Poiss_x.append(x)
        Poiss_y.append(y)
        self.T = np.linspace(0,2*np.pi,256)
        if(orientation==0):
            self.X = (l)*np.cos(self.T)+x 
            self.Y = L*np.sin(self.T)+y
        elif(orientation==1):
            print("er")
            self.X = l2*np.cos(self.T)+x 
            self.Y = L2*np.sin(self.T)+y

#Poissontest=Poisson(distanceMax/2,-50,0)

for i in range(z,z+1):
    Xa,Ya=trace(0,-i)

fig, (axes0,axes1,axes2)  = plt.subplots(nrows=1, ncols=3)
axes0.set_title('Densité')
axes2.set_title('Grad de vitesse')
axes1.set_title('Tracé réel')

axes0.axis([0,distanceMax,profondeur,0],"equal")
axes1.axis([0,distanceMax,profondeur,0],"equal")
axes2.axis([1400,1600,profondeur,0])

axes1.grid(True)
pas0=5
c =['cyan', 'skyblue', 'blue', 'navy', 'black']
for i in range(0,Nbr_ondes-2*Nbr_ondes_simplifiees,pas):
    longueur=len(Xtot1[i])
    for j in range(pas0,longueur,pas0):
        axes1.plot(Xtot1[i][j-pas0:j],Ytot1[i][j-pas0:j],c='navy',alpha=0.2-j/(10*longueur)) 
        #print(j-pas0,j+1)         #linestyle ='-'        #cmap='viridis'

pcm=axes0.hist2d(Xa,Ya,70,cmap=plt.get_cmap('hot'),cmin=0,cmax=120)       #cmap=plt.cm.jet
#fig.colorbar(pcm[0], ax=axes0, extend='max')


Y=[-i for i in range(-profondeur)]
avancement,p,X,cst=0,0,[],vitesse_surface
for j in range(-profondeur):
    X.append(cst+gradient[p]*(Y[j]-avancement))
    if(-j<=Ord_coude[p+1]):
        p=p+1
        avancement=-j
        cst=X[-1]
axes2.plot(X,Y)
#fig.colorbar(pcm[0], ax=axes0, shrink=0.5, aspect=5)

fig.set_size_inches(12.5, 4.5, forward=True)
plt.show()


et voici le second programme :

from mpl_toolkits.mplot3d import axes3d
from random import random, randint
import matplotlib
import matplotlib.pyplot as plt
from math import acos, sqrt ,pi , cos ,sin
import numpy as np

long_max,prof=1000,-200
nb_seg,solAngl,coef_sed,noise=100,[],4,5

pas=int(long_max/nb_seg)

def sign():
    if(int(randint(0,1))):
        signe=-1
    else:
        signe=1
    return(signe)

def sol():
    verif=[]
    global long_max
    global nb_seg
    global pas
    global coef_sed
    mat,drap,solAngl,composition=[],0,[],0
    color=['#B37426','#FFEA53','#747474']
    for i in range(0,long_max,pas*noise):
        if(drap):
            break
        if(i==0):
            sol=[randint(-185,-175)]
        else:
            sol=[sol[-1]]
        coef=sign()*randint(0,500)*0.02/(nb_seg)
        cst=sol[-1]
        if((i/(pas*noise))%coef_sed==0):
            composition=randint(0,2)
        mat.append(composition)
        for j in range(pas*noise):
            solAngl.append(coef)
            sol.append(cst+coef*(j))
            verif.append(cst+coef*(j))
            if(cst+coef*(j)<prof):
                for p in range(j,pas*noise):
                    sol.append(prof)
                    solAngl.append(0)
                    verif.append(prof)
                break      
    return(mat,solAngl,verif)

def norme(A,B,C):
    return((1/np.linalg.norm(A))*A,(1/np.linalg.norm(B))*B,(1/np.linalg.norm(C))*C)

def colision(X1,Y1,intensity):
    global verif
    indice =int(X1[-1])
    intensity= intensity * 0.75 #mat[int(X1[-1]//(pas*noise))]
    U=np.array([X1[-1]-X1[-2],Y1[-1]-Y1[-2]])
    Vt=np.array([1,verif[indice]-verif[indice-1]])
    V=np.array([-verif[indice]+verif[indice-1],1])
    V,Vt,V=norme(V,Vt,V)
    W=2*np.vdot(U,Vt)*Vt -U      #W=np.vdot(U,Vt)*Vt-np.vdot(U,V)*V
    return(X1[-1],Y1[-1],W,intensity)

def ligne(x0,y0,vect,intensity):
    Xl,Yl=[x0],[y0]
    for i in range(1,100):
        Xl.append(vect[0]*i+x0)
        Yl.append(vect[1]*i+y0)
        intensity=-1
        i=+1
    return(Xl,Yl,intensity)

mat,solAngl,verif=sol()

fig = plt.figure(figsize=(10,6))
ax1 = fig.add_subplot(111, projection='3d')

n=long_max/len(verif)
verif=np.asarray(verif)

x = np.arange(0,long_max,n)
y = np.arange(0,long_max,n)
X,Y = np.meshgrid(x,y)
Z = verif*np.ones((len(verif),len(verif)))

mycmap = plt.get_cmap('twilight_shifted')
ax1.set_title('sol marin')
surf1 = ax1.plot_surface(X, Y, Z, cmap=mycmap, antialiased=True)
fig.colorbar(surf1, ax=ax1, shrink=0.5, aspect=5)

for k in range(3,7):
    pastest=350
    X,Y=[],[]
    y=prof/2
    x,a=k*100,randint(1,5)
    for i in range(pastest):
        if(y<verif[int(x+i/a)]):
            break
        else:
            Y.append(y)
            X.append(x+i/a)
            y=y-(-prof/2)/pastest

    x,y,W,inten=colision(X,Y,100)
    X1,Y1,intensi=ligne(x,y,W,inten)
    n=len(X)
    ao=long_max/2
    ao1=randint(ao-100,ao+100)
    ax1.plot(X1,[ao1]*100,Y1,color='black')
    ax1.plot(X,[ao1]*n,Y,color='black')


ax1.set_xlabel('X en m')
ax1.set_xlim(0,long_max)
ax1.set_ylabel('Y en m')
ax1.set_ylim(0,long_max)
ax1.set_zlabel('Z en m')
ax1.set_zlim(-200,0)


plt.show()

2 réponses

Messages postés
3604
Date d'inscription
jeudi 16 juin 2005
Statut
Membre
Dernière intervention
7 avril 2020
944
Bonjour,

Je pense que ta question est beaucoup trop générale pour espérer obtenir une réponse ici.
On ne sait même pas ce qui te bloque : un point précis d'implémentation ? la formulation mathématique de ton diagramme cible ?

Xavier
Messages postés
4
Date d'inscription
mardi 17 mars 2020
Statut
Membre
Dernière intervention
18 mars 2020

Bonjour,
Merci de votre réponse. Effectivement mon sujet manque de clarté. Je précise. L’organisations générales, ce projet est mon premier « gros » code et je pense que sa complexité pourrait être diminué en optimisant sont fonctionnement mais je ne vois pas comment changer cela. Je n’arrive pas non plus à ajouter une colorbar pour le graphique hist2d. Autre point, j’aimerais me rapprocher du résultat souhaité en termes de représentation graphique.