Pb code calcul de la diffraction fraunhofer par methode numerique

Pb code calcul de la diffraction fraunhofer par methode numerique - C++ - Programmation

Marsh Posté le 08-01-2008 à 10:23:54    

Bonjour à tous,
 
Voila je post ce message car j'ai un code en langage c à réaliser qui calcule par une méthode numérique la diffraction de fraunhofer. Voici un sujet similaire au mien : http://www.impmc.jussieu.fr/impmc/ [...] 5/td5.html( attention monsujet ne comporte que la partie I).
J'ai utilisé la méthode de simpson pour faire mon calcul d'intégrale. Mon problème ce que mon code me rend une intensité nulle et je ne comprend pas pourquoi ??? :??:. J'ai utilisé 4 fonctions ( F1, F2, F3, F4) pour séparer la partie cosinus et la partie sinus et aussi pour séparer les x des y.
 
Voila si quelqu'un pourrait m'aider un peu car je patauge depuis un petit moment ?!?! :pt1cable: ^^
Merci d'avance !  
 
Voici la code que j'ai réalisé :

Code :
  1. #include <stdio.h>
  2. #include <math.h>
  3.    
  4. #define pi M_PI
  5. /* Longueur d'onde de la source lumineuse */
  6. #define lambda 500e-9
  7. /* Focale de la lentille */
  8. #define f 0.15
  11.        
  12. double F1(double x, double X)
  13. {
  14.     return cos((2*pi/(lambda*f))*(x*X));
  15. }
  17. double F2(double x, double X)
  18. {
  19.     return sin((2*pi/(lambda*f))*(x*X));
  20. } 
  21. double F3(double y, double Y)
  22. {
  23.     return cos((2*pi/(lambda*f))*(y*Y));
  24. }
  26. double F4(double y, double Y)
  27. {
  28.     return sin((2*pi/(lambda*f))*(y*Y));
  29. } 
  30. double simpson(int N, double x1, double x2,double y1, double y2, double X, double Y)
  32. {
  33.     double xx,yy, h1, IappF10, IappF11, IappF12, IappF1, h2, IappF20, IappF21, IappF22, IappF2, IappF30, IappF31, IappF32, IappF3, IappF40, IappF41, IappF42, IappF4,Iapp;
  34.     int NN, i;
  35.    
  36.     // Etape 1
  37.     h1 = (x2 - x1) / N;
  38.     h2 = (y2 - y1) / N;
  40.     // Etape 2
  41.     IappF10 = F1(x1,X) + F1(x2,X);
  42.     IappF11 = 0.0;
  43.     IappF12 = 0.0;
  44.     IappF20 = F2(x1,X) + F2(x2,X);
  45.     IappF21 = 0.0;
  46.     IappF22 = 0.0;
  47.     IappF30 = F3(y1,Y) + F3(y2,Y);
  48.     IappF31 = 0.0;
  49.     IappF32 = 0.0;
  50.     IappF40 = F4(y1,Y) + F4(y2,Y);
  51.     IappF41 = 0.0;
  52.     IappF42 = 0.0;
  53.     // Etape 3
  54.     NN = N -1;
  55.     for (i=1; i<=NN; i++)
  56.     {
  57. // Etape 4
  58. xx = x1 + i*h1;
  59. yy = y1 + i*h2;
  60.    
  61. // Etape 5
  62. if ((i%2) == 0)
  63. {
  64.     IappF12 = IappF12 + F1(xx,X);
  65.     IappF22 = IappF22 + F2(xx,X);
  66.     IappF32 = IappF32 + F3(yy,Y);
  67.     IappF42 = IappF42 + F4(yy,Y);
  68. }
  69. else
  70. {
  71.     IappF11 = IappF11 + F1(xx,X);
  72.     IappF21 = IappF21 + F2(xx,X);
  73.     IappF31 = IappF31 + F3(yy,Y);
  74.     IappF41 = IappF41 + F4(yy,Y);
  75. }
  76.     }
  77.    
  78.     // Etape 6
  79.     IappF1 = (IappF10 + 2.0 * IappF12 + 4.0 * IappF11) * h1 / 3.0;
  80.     IappF2 = (IappF20 + 2.0 * IappF22 + 4.0 * IappF21) * h1 / 3.0;
  81.     IappF3 = (IappF30 + 2.0 * IappF32 + 4.0 * IappF31) * h2 / 3.0;
  82.     IappF4 = (IappF40 + 2.0 * IappF42 + 4.0 * IappF41) * h2 / 3.0;
  83.     // Etape 7
  85.     Iapp = (IappF1*IappF2-IappF3*IappF4)*(IappF1*IappF2+IappF3*IappF4)+(IappF3*IappF2+IappF1*IappF4)*(IappF3*IappF2+IappF1*IappF4);
  87.     // Etape 8
  88.     return (Iapp);
  89. }
  91. int main(int argc, char *argv[])
  92. { 
  94.     int i,n;
  95.     double X,Y,x1,x2,y1,y2;
  96.    
  97.     x1=-0.000025;       /* Borne inferieure sur x*/
  98.     x2=+0.000025;       /* Borne superieure sur x*/
  99.     y1=-0.0025;         /* Borne inferieure sur y*/
  100.     y2=-0.0025;         /* Borne superieure sur y*/
  101.     n=100;              /* Nombre d'iteration */
  102.    
  103.     double I;
  105. /* Creation d'un fichier de donnée contenant l'intensité I(X,Y) demandé */
  106.     FILE *fp;
  107.     fp=fopen("intensite.data","w" );
  108.     for(X=0 ; X<=49 ; X=X+0.1)
  109.     {
  110. for(Y=0; Y<=49; Y=Y+0.1)
  111. {
  112.      I=simpson(n,x1,x2,y1,y2,X,Y);
  113.     fprintf(fp,"%lf %lf %lf\n",I,X,Y);
  115. }
  116.     }
  117.     fclose(fp);
  119. }


Reply

Marsh Posté le 08-01-2008 à 10:23:54   

Reply

Marsh Posté le 08-01-2008 à 10:31:21    

ta fonction simpson me parait douteuse à vue de nez.
Tu l'as sort de ou ?
 
 
Et accessoirement, c'est du C pas du C++


Message édité par Joel F le 08-01-2008 à 10:54:46
Reply

Marsh Posté le 08-01-2008 à 12:42:43    

voici le lien ou j'ai trouvé ma fonction simpson que j'ai ensuite adapté à mon pb :  
http://www.cppfrance.com/codes/INT [...] 27783.aspx

Reply

Marsh Posté le 08-01-2008 à 13:42:09    

Code :
  1. #     y1=-0.0025;         /* Borne inferieure sur y*/
  2. #     y2=-0.0025;         /* Borne superieure sur y*/


si c'est une intégration sur un plan, avec ces bornes ça fera 0 pour toute fonction.

Reply

Marsh Posté le 08-01-2008 à 15:04:22    

ah ouais merci erreur de frappe je suppose :p
je vais corriger ca !

Reply

Marsh Posté le 09-01-2008 à 02:24:20    

theoriquement ce genre de calcul est possible à coup de transformation de fourier, encore plus pour le cas fraunhofer car tu fais des hypothèses sur le fait que tu te trouves à r>>lambda. Passer par du simpson ne semble pas etre la meilleur idée! (apres j'ai pas vraiment tout lu en detail, je parle juste par experience)  [:spamafote]


---------------
--- WinSplit Revolution ---
Reply

Marsh Posté le 10-01-2008 à 09:54:33    

ouais va falloir que je modifie qque chose car malgré la correction sur les bornes de y ca me rend toujours une intensité nulle ...

Reply

Marsh Posté le 10-01-2008 à 10:46:13    

je viens d'installer la librairie fftw3. Quelqu'un aurait un site ou un document qui m'expliquerait comment l'utiliser ?? mci d'avance

Reply

Sujets relatifs:

Leave a Replay

Make sure you enter the(*)required information where indicate.HTML code is not allowed