I'm currently working with wavefront and psf.

I'm trying to convert a wavefront map extracted from opticstudio to a psf using a rather complicated but accurate numerical method, not a linear approximation method, but I'm not getting the right results. 

I am attaching the matlab code I am currently using.

Please let me know if there are any fixes.

Thanks.

CODE:

function psf_out=wavefront2psf(wvf, center_wvf, data_space_wvf, lambda, film_Dist, rect, data_space_psf)
    
% wvf: wavefront map of n x n size from opticsutdio
% center_wvf: center coordinate of the wavefront map
% data_space_wvf: data space of the wavefront map
% lambda: the wavelength
% film_Dist: distance between film and the last lens
% rect : point spread are estimated within the rectangle
% data_space_psf: data space of the point spread function

    dr=1;
    min_dA=2;

    cx=center_wvf(1);
    cy=center_wvf(2);
    
    rr=[];
    th=[];

    RR=0;
    TH=0;
    AA=pi;

    [h_wf, ~]=size(wvf);
    
    
    for p=1:dr:((h_wf/2)-3)
        r=p;

        r1=r;
        r2=r+dr;

        area_dR=pi*(r2^2-r1^2);
        n_th=ceil(area_dR/min_dA);

        n_th=ceil(n_th/4)*4;

        if mod(n_th,2)==1
            n_th=n_th+1;
        end

        ths=linspace(0, 2*pi, n_th+1);
        ths=ths(1:end-1);
        dth=ths(2)-ths(1);
        dA=dth/2*(r2^2-r1^2);

        r_=(r1+r2)/2*ones(size(ths));
        a_=dA*ones(size(ths));

        RR=[RR r_];
        TH=[TH ths];
        AA=[AA a_];

    end


    [xc, yc]=pol2cart(TH, RR);

    xc=xc+cx;
    yc=yc+cy;

    wsf_rth=interp2(wvf, xc, yc);


    x_=(xc-cx)*data_space_wvf;
    y_=(yc-cy)*data_space_wvf;
        
    hh=[rect(1):data_space_psf:rect(2)];
    ww=[rect(3):data_space_psf:rect(4)];
    
    phase_sum=zeros(size(hh,2), size(ww,2));
    [h_psf, w_psf]=size(phase_sum);
    
    
    for p=1:h_psf
        for q=1:w_psf

            px=(p-h_psf/2)*data_space_psf;
            py=(q-w_psf/2)*data_space_psf;
            ipx=x_-px;
            ipy=y_-py;

            path=sqrt(ipx.^2+ipy.^2+(film_Dist).^2);
            phase=(wsf_rth)+(path-film_Dist)/(lambda);
            phase_sum(p,q)=sum(AA.*exp(i*phase(:)'));
        end
    end
    
    psf_out=sqrt(abs(phase_sum).^2);
    psf_out=psf_out/sum(psf_out(:));

end

[Mod note: moved to more appropriate forum for ZOS-API-related discussions.]