Jm2600
Messages postés4Date d'inscriptionmardi 17 mars 2020StatutMembreDernière intervention18 mars 2020
-
17 mars 2020 à 23:08
Jm2600
Messages postés4Date d'inscriptionmardi 17 mars 2020StatutMembreDernière intervention18 mars 2020
-
18 mars 2020 à 09:28
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()
Reivax962
Messages postés3671Date d'inscriptionjeudi 16 juin 2005StatutMembreDernière intervention11 février 20211 011 18 mars 2020 à 08:45
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 ?
Jm2600
Messages postés4Date d'inscriptionmardi 17 mars 2020StatutMembreDernière intervention18 mars 2020 18 mars 2020 à 09:28
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.