IPSF_Vectorial_1Int(na, wl, n1, n2, r, z, d, ang,deltaZ)

function [r z I] = IPSF_Vectorial_1Int(na, wl, n1, n2, r, z, d, ang,deltaZ)

     [r z] = meshgrid(r,z);
    if n2 < n1
        a = asin(n2./n1);
    else
        a = asin(na/n1);
    end
    k = (2*pi)/wl;
    
    I0 = quadv(@E0field,0,a);
    I1 = quadv(@E1field,0,a);
    I2 = quadv(@E2field,0,a);
    
    I = abs(I0).^2 + 4 .* abs(I1).^2 .* (cos(ang)).^2 + abs(I2).^2 + 2 .* cos(2 .* ang) .* real(I0 .* conj(I2));
    
    function E0 = E0field(th1)
        % Theta 2 from snells law
        th2 = asin((n1.*sin(th1))./n2);
        % Abberation Function
%        Phi=0;
        %Phi =deltaZ.* sin(th1).^2;
         Phi = -d.*(n1.*cos(th1)-n2.*cos(th2))+deltaZ.*(1-cos(th1));
        % Fresnel coefficients
        [ts tp] = FresnelCoeff(th1, th2);
        % Electric Field
        E0 = (cos(th1).^(1/2)) .* sin(th1) .* (ts + tp .* cos(th2)) .* besselj(0, k .* r .* n1 .* sin(th1)) .* exp(-i .* k .* Phi - i.* k .* z .* n2 .* cos(th2));
    end

    function E1 = E1field(th1)
        % Theta 2 from snells law
        th2 = asin((n1.*sin(th1))./n2);
        
        % Abberation Function
%        Phi=0;
        %Phi =deltaZ.* sin(th1).^2;
          Phi = -d.*(n1.*cos(th1)-n2.*cos(th2))+deltaZ.*(1-cos(th1));
        % Fresnel coefficients
        [ts tp] = FresnelCoeff(th1, th2);
        % Electric Field
        E1 = (cos(th1).^(1/2)) .* sin(th1) .* (tp .* sin(th2)) .* besselj(1, k .* r .* n1 .* sin(th1)) .* exp(-i .* k .* Phi - i.* k .* z .* n2 .* cos(th2));
    end

    function E2 = E2field(th1)
        % Theta 2 from snells law
        th2 = asin((n1.*sin(th1))./n2);
        % Abberation Function
%        Phi=0;
        %   Phi =deltaZ.* sin(th1).^2;
         Phi = -d.*(n1.*cos(th1)-n2.*cos(th2))+deltaZ.*(1-cos(th1));
        % Fresnel coefficients
        [ts tp] = FresnelCoeff(th1, th2);
        % Electric Field
        E2 = (cos(th1).^(1/2)) .* sin(th1) .* (ts - tp .* cos(th2)) .* besselj(2, k .* r .* n1 .* sin(th1)) .* exp(-i .* k .* Phi - i.* k .* z .* n2 .* cos(th2));
    end
end