Hi to all,
The problem was in the program where I use quadgk. It was just lacking the index into variable phi(p) and now all is good:
Sorry for the inconvenience.
Carolina R
clear all
close all
clc
%w = warning('query','last')
id1='MATLAB:quadgk:MaxIntervalCountReached';
id2='MATLAB:quadgk:MinStepSize';
warning('off',id1)
warning('off',id2)
%%%Optical parameters
lambda=532e-9;
k=2*pi/lambda;
n1=1.518;
k1=n1.*k;
NA=1*n1;
f=0.01;
alfa=asin(NA/n1);
%%%%Pupil radius (m)
R=f.*sin(alfa);
%%%%Axis (m)
N=64;
x=linspace(-N/2,N/2-1,N);
x=2/N.*x.*1.25e-6;
x=repmat(x,N,1);
y=x';
phi=angle(x+1i.*y);
r=sqrt(x.^2+y.^2);
%%%%In tight focusing work we use theta variable
theta=asin(r./f);
%%%%Vortex radii
rho1=0.00055.*R;
rho2=0.91.*rho1;
%%%%Topological charge, m=0 inclusive
m=0;
%%%%Polarization order, n=0 inclusive
n=0;
%%%%Integration order function I0_nu
nu=m+n;
%%%%Initializing
I0_nu1=zeros(N);
%%%%Focal point distance z=0
z=0;
for p=1:N*N
I0_nu1(p) =quadgk(@(theta)I0_nu(theta,nu,z,r(p),k1,f,m,rho1,rho2,R,phi(p)),0,alfa);
if mod(p,8*N)==0;
%%%%Percentage tellboard
fprintf('completo %3.0f%%\n',round(p.*100./(N*N)))
end
end