Bonjour,
Me voici devant un probleme tres certainement numerique que je n'arrive pas a resoudre.. si quelqu'un pouvait m'aider ce serait le pied!
Je suis censee integrer une enorme fonction comprenant deux integrales triples et une simple.. Pour cela j'utilise les librairies de Numerical Recipes mais j'obtiens de mauvais resultats dans certains cas et je ne sais pas d'ou vient l'erreur..
je joins donc mon programme .. si quelqu'un avait la gentillesse d'y jeter un coup d'oeil...
Mon probleme est que si je choisis un kTe <2.0 je n'obtient plus le meme type de courbe qu'en kTe=5.0 par exemple... (soit une partie negative et une autre positive) ..
Il suffit de lancer mon programme et de ploter le fichier szre.dat..
De plus il semble avoir des distortions pour toutes les valeurs de kTe qui ne devrait pas exister!!!..
Voili...
Je sais que c'est un peu dur de rentrer comme ca dans un programme mais c'est aussi difficile de mieux expliquer...
En revanche n'hesitez pas pour les questions...
A bientot j'espere!
Annoka
#include<stdio.h>
#include<math.h>
#include"/home/anna/recipes_c-ansi/include/nrutil.h"
#include"/home/anna/recipes_c-ansi/include/nr.h"
#define mec2 511.0
#define kTe 1.0
#define k 1.38066
#define h 6.62608
#define c 2.99792
#define T 2.726
float xx;
float eta=mec2/kTe;
int e;
float Io()
{
float fact,io,tau,fxx;
fact=-69+68-16;
tau=0.001;
io=2.*pow(k*T,3.)*pow(10.,fact)/pow(h*c,2.);
fxx=pow(xx,3.)/(exp(xx)-1.);
return(fxx*tau*io);
}
/*fonctions a integrer*/
float func(float t,float beta,float mu)
{
float f1,f2;
if(e==1)
{
f1=((t*(exp(xx)-1.))/(exp(xx*t)-1.))*(beta*(1./sqrt(1.-pow(beta,2.)))*exp(-eta*((1./sqrt(1.-pow(beta,2.)))-1.)))*(((((3.*pow(mu,2.))-1.)*pow(((1.-(beta*mu))/t)-1.,2.)/pow(beta,2.))+3.-pow(mu,2.))/pow(1.-(beta*mu),2.));
return(f1);
}
else
{
f2=(((exp(xx)-1.))/(pow(t,3.)*(exp(xx/t)-1.)))*(beta*(1./sqrt(1.-pow(beta,2.)))*exp(-eta*((1./sqrt(1.-pow(beta,2.)))-1.)))*(((((3.*pow(mu,2.))-1.)*pow(((1.-(beta*mu))/(1./t))-1.,2.)/pow(beta,2.))+3.-pow(mu,2.))/pow(1.-(beta*mu),2.));
return(f2);
}
}
/*bornes*/
float z1(float t,float beta)
{
float mum;
if(e==1)
{
mum=(1.-t-(t*beta))/beta;
return(mum);
}
else
{
return(-1.0);
}
}
float z2(float t,float beta)
{
float muM;
if(e==1)
{
return(1.0);
}
else
{
muM=(t-1.+beta)/(t*beta);
return(muM);
}
}
float yy1(float t)
{
float betam;
betam=(1.-t)/(1.+t);
return(betam);
}
float yy2(float t)
{
return(0.99);
}
/*fonction a intergrer*/
float A(float beta)
{
float AA;
AA=pow(beta,2.)*pow((1./sqrt(1.-pow(beta,2.))),5.)*exp(-eta*((1./sqrt(1.-pow(beta,2.)))-1.));
return(AA);
}
int main(void)
{
FILE *fich;
float fi1,fi2,AA,fi;
int i;
xx=0.0000001;
fich=fopen("szre.dat","w");
for(i=0;i<5000;i++)
{
e=1;
fi1=quad3d(func,0.0,1.0);
e=2;
fi2=quad3d(func,0.0,1.0);
AA=3./(32.*qromb(A,1.e-5,0.9999999));
fi=AA*(fi1+fi2);
printf("fi-1=%f x=%f\n",Io()*(fi-1.)*1.e18,xx);
fprintf(fich,"%f %f\n",xx,Io()*(fi-1.)*1.e18);
xx=xx+0.004;
}
fclose(fich);
return(0);
}