Avis et conseils Simulation phénomène : Python

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()