Voici: mon but est de simuler un système exoplanétaire avec la planète intérieure qui est une supergéante gazeuse ayant 2% de la masse du Soleil, située à un rayon orbital de 1 UA, L'on veut calculer le spectre de puissance (défini comme Ak^2+Bk^2, Ak et Bk étant des coefficients de Fourier) mais j'en arrive à une erreur de segmentation lorsque je veux le faire porter en graphique. (Veuillez m'excuser pour les paramètres de pllab; je n'ai pas encore ajusté le graphique mais je vais le faire une fois l'erreur de segmentation réglée)
Comment faire pour régler ce fâcheux problème?Code:#include <stdio.h> #include <math.h> #include <plplot.h> #define NPAS 41000 /* Nombre de pas de temps */ /* Ce programme integre les deux EDO couplees definisant */ /* le probleme du pendule lineaire par la methode d'Euler explicite */ int main(void) { /* Declarations ---------------------------------------------------- */ long c=0, k=0, rf, kx ; /* Compteurs d'iterations */ float h=1./1500 ; /* Taille du pas de temps */ float mplus, nplus, G, M, splus, wplus; double x, y, i, j, xmin, xmax, ymin, ymax, Mj, rp, rx, ry, sx, sy, tx, ty, fbase; float t[NPAS], moyenne1, moyenne2 ; /* Le temps */ double u, v, m, n, w, s, rayon1[NPAS], rayon2[NPAS] ; /* Les deux variables dependantes */ double Pk[NPAS], ak[NPAS], bk[NPAS]; xmin=-2.5;xmax=2.5; ymin=-2.5; ymax=2.5; plscolbg(255,255,255); plinit(); plenv(xmin, xmax, ymin, ymax, 1, 1); pllab("Orbite horizontale (UA)","Orbite verticale (UA)","Orbite en fonction du temps"); /* Executable ------------------------------------------------------ */ i=0. ; j=-1; m=0.9*6.2832; n=0.; /* Conditions initiales pour super-Jupiter */ x = 0. ; y = -1.6; s=6.2832/1.2649; w=0. ; /* Conditions initiales pour petite-planete */ G= 1.96*pow(10.,-29); M= 1.98*pow(10.,30); Mj=0.02*M; t[0]=0. ; for (k=0 ; k < NPAS ; k++ ) /* Boucle temporelle */ { rp =pow(((pow(x,2))+(pow(y,2))), -1.5); /* Variables intermediaires */ rx =pow((pow(((x-i)+h*(s-m)),2)+(pow(((y-j)+h*(w-n)),2))),1.5); ry =pow(((pow((x+h*s),2))+(pow((y+h*w),2))),1.5); tx =x-i+h*s-h*m; ty =y-j+h*w-h*n; sx =pow(((x-i)+h*s-h*m),2); sy =pow(((y-j)+h*w-h*n),2); splus =s - h*G*M*(x+h*s)/ry - h*G*Mj*((x-i)+h*(s-m))/rx; wplus =w - h*G*M*(y+h*w)/ry - h*G*Mj*((y-j)+h*(w-n))/rx; x =x+h/2 * ( s+splus ); /* petite-planete, x,y */ y =y+h/2 * ( w+wplus ); s =s - (h/2)*G*M*((x*rp) + ((x+h*s)/(pow(((pow((x+h*s),2))+(pow((y+h*w),2))),1.5))) +0.02*(((x-i)*pow(((pow((x-i),2))+(pow((y-j),2))), -1.5))+(tx*(pow((sx+sy),-1.5))))); w =w - (h/2)*G*M*((y*rp) + ((y+h*w)/(pow(((pow((x+h*s),2))+(pow((y+h*w),2))),1.5))) +0.02*(((y-j)*pow(((pow((x-i),2))+(pow((y-j),2))), -1.5))+(ty*(pow((sx+sy),-1.5))))); mplus =m - h*G*M*((i+h*m)/(pow(((pow((i+h*m),2))+(pow((j+h*n),2))),1.5))) ; nplus =n - h*G*M*((j+h*n)/(pow(((pow((i+h*m),2))+(pow((j+h*n),2))),1.5))) ; i =i+h/2 * ( m+mplus ); /* super-Jupiter, i,j */ j =j+h/2 * ( n+nplus ); m =m - (h/2)*G*M*(((i) *pow(((pow(i,2))+(pow(j,2))), -1.5)) + ((i+h*m)*(pow(((pow((i+h*m),2))+(pow((j+h*n),2))),-1.5)))); n =n - (h/2)*G*M*(((j) *pow(((pow(i,2))+(pow(j,2))), -1.5)) + ((j+h*n)*(pow(((pow((i+h*m),2))+(pow((j+h*n),2))),-1.5)))); rayon1[k]=pow(((pow(y,2))+(pow(x,2))), 0.5); rayon2[k]=pow(((pow(i,2))+(pow(j,2))), 0.5); moyenne1+=rayon1[k]/NPAS; moyenne2+=rayon2[k]/NPAS; t[k+1] =t[k] +h ; } fbase=1500*6.2832/NPAS; /*La frequence de base */ for (kx=1; kx<NPAS;kx++) /*Les nombres d'onde */ { ak[kx]=0.; bk[kx]=0.; for (k=0; k<fabs(NPAS/10); kx++) { ak[kx]+=(3000./NPAS)*(rayon1[k]-moyenne1)*cos(kx*fbase*t[k]); bk[kx]+=(3000./NPAS)*(rayon2[k]-moyenne1)*sin(kx*fbase*t[k]); } Pk[kx]=ak[kx]*ak[kx]+bk[kx]*bk[kx]; /* Spectre de puissance */ } for (kx=1; kx<fabs(NPAS/10);kx++) { plcol0(3); pljoin (kx*fbase,Pk[kx],(kx+1)*fbase,Pk[kx+1]); } plend(); }
-----
