Использование quad2d

X = quad2d (@myXcul,-2,2,0,2*pi, par1,par2,par3);

xi=x1+m*t;
yi=y1+n*t;
zi=z1+p*t;

function E = myE(x,y,z,xi,yi,zi)
% x,y,z - координаты точки сферы, для которой считаем вектор Е
% E = [Ex, Ey, Ez] - вектор излучения
 
[Nx,Ny,Nz] = surfnorm(xi,yi,zi); % векторы нормали
Ex = 0; Ey = 0; Ez =0; 
for i = 1:size(xi,1)
    for j = 1:size(xi,2)
        Bxy = [x-xi(i,j), y-yi(i,j)];
        Bxz = [x-xi(i,j), z-zi(i,j)]; % вектор от точки поверхности к искомой точке
        Nxy = [Nx(i,j), Ny(i,j)];
        Nxz = [Nx(i,j), Nz(i,j)];
        % угол между нормалью излучателяи вектором к точке
        fi_xoy = acos( sum(Bxy.*Nxy)/(norm(Bxy)*norm(Nxy)) ); % вплостоксти XOY
        fi_xoz = acos( sum(Bxz.*Nxz)/(norm(Bxz)*norm(Nxz)) ); % вплостоксти XOZ
        % суммируем проекции векторов на оси
        Dx = norm(Nxy)*cos(fi_xoy) * cos(atan2(Bxy(2),Bxy(1)));
        Dy = norm(Nxy)*cos(fi_xoy) * sin(atan2(Bxy(2),Bxy(1))); 
        Dz = norm(Nxz)*sin(fi_xoz) * sin(atan2(Bxz(2),Bxz(1)));
        Ex = Ex + Dx*((fi_xoy>=-pi/2)&(fi_xoy<=pi/2));  
        Ey = Ey + Dy*((fi_xoy>=0)&(fi_xoy<=pi));
        Ez = Ez + Dz*((fi_xoz>=0)&(fi_xoz<=pi)); 
    end
end
E = [Ex, Ey, Ez]/(size(xi,1)*size(xi,2)); % нормируем вектор

clear
clc
 
% сфера
R=1;
u1 = linspace(0,2*pi,50)';
v1 = linspace(0,2*pi,20);
X = R*sin(u1)*cos(v1);
Y = R*sin(u1)*sin(v1);
Z = R*cos(u1)*ones(size(v1));
hold on
surf(X,Y,Z,'Facealpha',0.0);
%shading interp
%випромінювач
u = linspace(-0.5,0.5,20)';
v = linspace(0,2*pi,20);
rr = 0.1;
xi = rr*ones(size(u))*cos(v);
yi = rr*ones(size(u))*sin(v);
zi = u*ones(size(v));
surf(xi,yi,zi)
grid on
axis equal
xlabel('ось  X'),ylabel('ось Y'),zlabel('ось Z')
% plot3(x11,y11,z11)
for i=1:size(X,1)
    for j=1:size(X,2)
        E = myE(X(i,j),Y(i,j),Z(i,j), xi,yi,zi);
        XX(i,j) = E(1);
        YY(i,j) = E(2);
        ZZ(i,j) = E(3);
    end
end
quiver3 (X,Y,Z,XX,YY,ZZ,3)
view(40,20)
% mesh(X1+XX1.*0.004-0.05,Y1+YY1.*0.004-5,Z1+ZZ1.*0.004)

Ex = Ex + Dx*((fi_xoy>=-pi/2)&(fi_xoy<=pi/2));  
Ey = Ey + Dy*((fi_xoy>=-pi/2)&(fi_xoy<=pi/2));
Ez = Ez + Dz*((fi_xoz>=-pi/2)&(fi_xoz<=pi/2));