Dans le cadre d'un projet, je dois tracer la courbe de transit d'une exoplanète grace à python. J'ai rédigé ce programme et je ne comprends pas pourquoi il ne fonctionne pas. Si quelqu'un pouvait m'aider à comprendre mon erreur. Merci d'avance !!
Code:import numpy as np import matplotlib.pyplot as plt from math import cos from math import sin from math import sqrt from math import pi from math import atan #fonction de projection orthogonale de l'ellipse dans le plan d'observation def Proj(x,y,z,n1,n2,n3,position): i=0 Projecx=[] Projecy=[] Projecz=[] while i<position: Projecx.append(x[i]-(n1*x[i]+n2*y[i]+n3*z[i])*n1) Projecy.append(y[i]-(n1*x[i]+n2*y[i]+n3*z[i])*n2) Projecz.append(z[i]-(n1*x[i]+n2*y[i]+n3*z[i])*n3) #print(Projecx[i],Projecy[i],Projecz[i]) assert(-10**(-10)<=Projecx[i]*o1+Projecy[i]*o2+Projecz[i]*o3<=10**(-10)) i+=1 return(Projecx,Projecy,Projecz) #fonction de rotation de l'ellipse selon l'axe Z de notre repere def RZ(Matrice_X,Matrice_Y,Matrice_Z,position,angle): i=0 X_r=[] Y_r=[] Z_r=[] while i<position: X_r.append(Matrice_X[i]*cos(angle)-Matrice_Y[i]*sin(angle)) Y_r.append(Matrice_X[i]*sin(angle)+Matrice_Y[i]*cos(angle)) Z_r.append(Matrice_Z[i]) i+=1 Matrice_X=X_r Matrice_Y=Y_r Matrice_Z=Z_r return(Matrice_X,Matrice_Y,Matrice_Z) #fonction de rotation de l'ellipse selon l'axe X de notre repere def RX(Matrice_X,Matrice_Y,Matrice_Z,position,angle): i=0 X_r=[] Y_r=[] Z_r=[] while i<position: X_r.append(Matrice_X[i]) Y_r.append(Matrice_Y[i]*cos(angle)-Matrice_Z[i]*sin(angle)) Z_r.append(Matrice_Y[i]*sin(angle)+Matrice_Z[i]*cos(angle)) i+=1 Matrice_X=X_r Matrice_Y=Y_r Matrice_Z=Z_r return(Matrice_X,Matrice_Y,Matrice_Z) #donnée de notre étude #parametres de l'ellipse a=150 #demi grand axe de l'ellipse caracterisant la trajectoire de l'exo planete e=0.8 #excentricté de l'ellipse b=a*((1-(e)**2)**(1/2)) #demi petit axe de l'ellipse caracterisant la trajectoire de l'exo planete c=a*e P=(a*a-c*c)/a R=10 #Rayon de notre etoile Rp=1 #Rayon de notre planete G = 6.67*10**-11 # Constante gravitationnelle universelle en m3/kg/s2 M=10# Masse de l'étoile en kg T=sqrt((a**3)*4*pi**2/G*M) # Période orbitale à partir de la troisième loi de Kepler theta=2*pi pas=0.05 C=[c,0,0] i=0 k=0 p=0 Obs=[-100,-100,-20] q=[] x=[] y=[] z=[] X=[] Y=[] Z=[] #Tracé de l'ellipse t=0 p=0 pas=80000 pas2=7000000 # Calcul de la vitesse orbitale à partir des lois de Kepler v = np.sqrt(G * M * (2 / np.sqrt(a) - 1 / a)) #Ellipse while t<T: x.append( a * (np.cos(2 * np.pi *t/T) - e)) y.append(a * np.sqrt(1 - e**2) * np.sin(2 * np.pi * t/T)) z.append(0) t+=pas p+=1 #Sphere: i=0 while i<T: while k<T/2: X.append(R*cos(i)*sin(k)) Y.append(R*sin(i)*sin(k)) Z.append(R*cos(k)) k+=pas2 k=0 i+=pas2 #projection de l'ellipse o1=float(Obs[0]) o2=float(Obs[1]) o3=float(Obs[2]) print(o1,o2,o3) N=Obs print(N) N2= np. linalg. norm(N) print(N2) N3=N/N2 print(N3) n1=float(N3[0]) n2=float(N3[1]) n3=float(N3[2]) print(n1,n2,n3) vec_o1=[] vec_o2=[] vec_o3=[] i=0 while i<o1 and i<o2 and i<o3: vec_o1.append(o1*i) vec_o2.append(o2*i) vec_o3.append(o3*i) i+=0.05 #Rotation de l'ellipse x1,y1,z1=RZ(x,y,z,p,pi/2) x2,y2,z2=RX(x1,y1,z1,p,pi/2) #Porjection de l'ellipse xp,yp,zp=Proj(x2,y2,z2,n1,n2,n3,p) #Plot 2D Ellipse: plt.plot(xp,yp,'ro',X,Y,) plt.xlabel('x') plt.ylabel('y') plt.show() #Plot 3D Ellipse + Sphere: fig = plt.figure("truc") ax = plt.axes(projection='3d') ax.set_xlabel('x') ax.set_ylabel('y') ax.set_zlabel('z') # ax.set_xlim3d(-10,10) # ax.set_ylim3d(-10,10) # ax.set_zlim3d(-10,10) #ax.plot(P1,P2,P3,'green') ax.scatter(o1,o2,o3) ax.plot(xp, yp, zp,'red') ax.plot(X,Y,Z,'yellow') ax.quiver(0,0,0,o1,o2,o3) # plt.plot(courbe) ax.view_init(elev=0, azim=20) # Ajuster l'angle de vue plt.title("Courbe 3D") plt.tight_layout() plt.show() ypos=[] xpos=[] zpos=[] o=0 for i in range(p): if zp[i]>0: zpos.append(zp[i]) ypos.append(y[i]) xpos.append(x[i]) o+=1 ypos, xpos = zip(*sorted(zip(ypos, xpos), reverse=True)) # Convertir les tuples en listes ypos = list(ypos) xpos = list(xpos) Max_x=xpos[0] Min_x=xpos[0] for i in range(o): if xpos[i]>Max_x: Max_x=xpos[i] if xpos[i]<Min_x: Min_x=xpos[i] Max_y=ypos[0] Min_y=ypos[0] for i in range(o): if ypos[i]>Max_y: Max_y=ypos[i] if ypos[i]<Min_y: Min_y=ypos[i] print(Max_x,Min_x,Max_y,Min_y) if abs(Max_x)>abs(Max_y) and abs(Max_x)>abs(Min_x) and abs(Max_x)>abs(Min_y): t=abs(int(Max_x)) if abs(Max_y)>abs(Max_x) and abs(Max_y)>abs(Min_x) and abs(Max_y)>abs(Min_y): t=abs(int(Max_y)) if abs(Min_x)>abs(Max_x) and abs(Min_x)>abs(Min_y) and abs(Min_x)>abs(Max_y): t=abs(int(Min_x)) if abs(Min_y)>abs(Max_x) and abs(Min_y)>abs(Min_x) and abs(Min_y)>abs(Max_y): t=abs(int(Min_y)) image = [[0 for j in range(t)] for i in range(t)] for i in range(t): for j in range(t): r = sqrt((i - t/2)**2 + (j - t/2)**2) if r > R/2: image[i][j] = 0 else: teta = atan(r/(R*100)) omega =atan(R/2/(R*100)) image[i][j] = int(1000*(0.3 + 0.93*sqrt((cos(teta))**2 - (cos(omega))**2)/sin(omega) - 0.23*((cos(teta))**2 - (cos(omega))**2)/(sin(omega))**2)) for row in image: print("\t".join(str(x) for x in row)) inttot = 0 for i in range(t): for j in range(t): inttot += image[i][j] print(inttot) courbe = np.zeros(t) for c in range(t): courbe[c] = inttot a, b = np.polyfit(xpos, ypos, 1) print("f(x) = {0:.2f}x + {1:.2f}".format(a, b)) for x in np.arange(-t/2 + Rp/2 + 1, t/2 - Rp/2): y=a*x+b ix = int(x + t/2) jy = int(y + t/2) print("x=", x, "y=", y, "ix=", ix, "jy=", jy) intmoins = 0 for i in range(max(0, ix - int(Rp/2)), min(t, ix + int(Rp/2) - 1)): for j in range(max(0, jy - int(Rp/2)), min(t, jy + int(Rp/2) - 1)): if (i - ix)**2 + (j - jy)**2 <= (Rp/2)**2: intmoins += image[i][j] courbe[ix] = inttot - intmoins # k=o/t # o=1 # for x in np.arange(-t/2 + Rp/2 + 1, t/2 - Rp/2): # ix = int(x + t/2) # ymatrice=ypos[int(k*o)] # jy = int(ymatrice + t/2) # print(ix,jy) # intmoins = 0 # for i in range(ix - int(Rp/2), ix + int(Rp/2) - 1): # for j in range(jy - int(Rp/2), jy + int(Rp/2) - 1): # if (i - ix)**2 + (j - jy)**2 <= (Rp/2)**2: # intmoins += image[i][j] # courbe[ix] = inttot - intmoins # o+=1 print(courbe) Int=plt.figure() plt.plot(courbe) plt.show() print('Fin du programme \n')
-----
Envoyé par hormonut3 