Пытался смоделировать дифракцию на круглом отверстии, но верный ответ получается только в центре. В чем проблема?

Question

[Polkovn1k]

Использовал интеграл Кирхгофа в виде суммы

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<complex.h>
#include<pthread.h>
#define a 0.5
#define PI 3.1415926535897932
#define P 300
#define P2 300
struct MatrixAndFile
{
    FILE*f3;
    struct Matrix** T;
    int m;

};
struct Matrix
{
    double x;
    double y;
    double R;
    double cosA;
    double z;

};

void* sumir (void* arg)
{
    int i,k,j,p,M,N0;
    double x1,y1,h,h2,k0,lm,s,cosB,dS,b,P0;
    double complex Sum;
    M=(((struct MatrixAndFile*)arg)->m);
    N0=2*P+1;
    h=a/P;
    h2=0.5/P2;;
    lm=632/pow(10,6);
    k0=2*PI/lm;
    b=0.25/(3*lm);
    dS=h*h;
    for(p=M;p<(M+1);p++)
    {
        y1=0.5-p*h2;
        for(k=0;k<P2+1;k++)
        {

            x1=-0.5+k*h2;
            Sum=0+0*I;
            for (i=0;i<N0;i++)
            {
                for (j=0;j<N0;j++)
                {
                s=pow(((((struct MatrixAndFile*)arg)->T[i][j].x)-x1),2)+pow(((((struct MatrixAndFile*)arg)->T[i][j].y)-y1),2)+pow(b,2);
                s=sqrt(s);
                cosB=b/s;
                P0=(((struct MatrixAndFile*)arg)->T[i][j].z)*(((((struct MatrixAndFile*)arg)->T[i][j].cosA)+cosB)*dS)/((((struct MatrixAndFile*)arg)->T[i][j].R)*s);
                Sum+=(double complex)(P0)*cexp(I*k0*(s+(((struct MatrixAndFile*)arg)->T[i][j].R)));
                }
            }
            Sum=(-I*Sum)/(2*lm);
            Sum=Sum*conj(Sum);
            Sum=cabs(Sum);
            fprintf(((struct MatrixAndFile*)arg)->f3,"%.d\t  %.d\t  %.10lf\n",p,k,creal(Sum));
            fprintf(((struct MatrixAndFile*)arg)->f3,"%.d\t  %.d\t  %.10lf\n",2*P2-p,k,creal(Sum));
            fprintf(((struct MatrixAndFile*)arg)->f3,"%.d\t  %.d\t  %.10lf\n",p,2*P2-k,creal(Sum));
            fprintf(((struct MatrixAndFile*)arg)->f3,"%.d\t  %.d\t  %.10lf\n",2*P2-p,2*P2-k,creal(Sum));
        }
    }
    return NULL;
}

int main()
{
    FILE *file;
    file=fopen("Kvadrat parallel2.txt","w");
    FILE *file3;
    file3=fopen("Kvadrat difrakcia parallel2.txt","w");

    int i,j,k,N,N0,N3,rc,NUMTHREADS;
    double x0,y0,h,R0;
    struct Matrix **A;
    NUMTHREADS=P2+1;

    h=a/P;
    N=P+1;
    N0=2*P+1;
    N3=2*P;
    R0=10000;


    A=malloc(N0*sizeof(struct Matrix*));
    for(i=0;i<N0;i++)
    {
        A[i]=(struct Matrix*)malloc(N0*sizeof(struct Matrix));
    }


    for (i=0;i<N;i++)
    {
        y0=0.5-i*h;
        for (j=0;j<N;j++)
        {
            x0=-0.5+j*h;
            A[i][j].x=x0;
            A[i][j].y=y0;
            A[i][j].R=(A[i][j].x)*(A[i][j].x)+(A[i][j].y)*(A[i][j].y)+R0*R0;
            A[i][j].R=sqrt(A[i][j].R);
            A[i][j].cosA=R0/A[i][j].R;
            if((A[i][j].x*A[i][j].x+A[i][j].y*A[i][j].y)<=0.25)
                A[i][j].z=100000;
            else
                A[i][j].z=0;
            fprintf(file,"%d\t  %d\t  %.10lf\n",i,j,A[i][j].z);

            A[N3-i][j].x=-x0;
            A[N3-i][j].y=y0;
            A[N3-i][j].R=A[i][j].R;
            A[N3-i][j].cosA=A[i][j].cosA;
            A[N3-i][j].z=A[i][j].z;
            fprintf(file,"%d\t  %d\t  %.10lf\n",N3-i,j,A[N3-i][j].z);

            A[i][N3-j].x=x0;
            A[i][N3-j].y=-y0;
            A[i][N3-j].R=A[i][j].R;
            A[i][N3-j].cosA=A[i][j].cosA;
            A[i][N3-j].z=A[i][j].z;
            fprintf(file,"%d\t  %d\t  %.10lf\n",i,N3-j,A[i][N3-j].z);

            A[N3-i][N3-j].x=-x0;
            A[N3-i][N3-j].y=-y0;
            A[N3-i][N3-j].R=A[i][j].R;
            A[N3-i][N3-j].cosA=A[i][j].cosA;
            A[N3-i][N3-j].z=A[i][j].z;
            fprintf(file,"%d\t  %d\t  %.10lf\n",N3-i,N3-j,A[N3-i][N3-j].z);
        }
    }

    struct MatrixAndFile *B = malloc((NUMTHREADS)*sizeof(struct MatrixAndFile));
    pthread_t threads[NUMTHREADS];

    for (k=0;k<NUMTHREADS;k++)
    {
        B[k].T=A;
        B[k].f3=file3;
        B[k].m=k;
        rc=pthread_create(&threads[k],NULL,sumir,(void*)&B[k]);
        if(rc)
       {
           printf("Error: unable to create thread %d\n",rc);
            exit(-1);
        }
    }
    for (k=0;k<NUMTHREADS;k++)
    {
        rc=pthread_join(threads[k],NULL);
    }


    fclose(file);


    for (i=0;i<N;i++)
        free(A[i]);
    free(A);
    free(B);
    fclose(file3);
    return 0;
}

Comments

[jcmvbkbc]

Пытался смоделировать дифракцию на круглом отверстии, но верный ответ получается только в центре.
Использовал интеграл Кирхгофа в виде суммы

Я предлагаю вам провести дополнительную работу и упростить вопрос до такого: "я пытаюсь посчитать такой-то интеграл чесленным методом, ожидаю такой-то результат, получаю такой-то, почему". Это позволит сконцентрироваться на коде вашей реализации, оставив переход от физики к математике на вас.