Operator '.^' is not supported for operands of type 'cell'

any idea to solve the error at the end of the code (visualisation part)?

close all
clear all
clc
%%

mu0 = 4*pi*1e-7; % Vs/Am
M0 = 1e3; % A/m
maxnum = 31;
rho1_min = 0;
rho1_max = 0.25;

xlimit = [-1, 1];
ylimit = xlimit;
zlimit = xlimit;
x = linspace(min(xlimit), max(xlimit), maxnum);
y = linspace(min(ylimit), max(ylimit), maxnum);
z = linspace(min(zlimit), max(zlimit), maxnum);
[Xg, Yg, Zg] = ndgrid(x, y, z);

rho = sqrt(Xg.^2 + Yg.^2 + Zg.^2);
phi = angle(Xg + 1i*Yg);
theta = angle(Zg + 1i*sqrt(Xg.^2 + Yg.^2));
%%

RHO = sqrt(x.^2 + y.^2 + z.^2);
THETA = linspace(0, pi, 31);                          % Trapz
PHI = linspace(0, 2*pi, 31);                          % Trapz

%% 
% Pre-allocate F_x , B1x and B1z as cell array
F_x = cell(numel(RHO), numel(THETA), numel(PHI));
B1x = cell(numel(RHO), numel(THETA), numel(PHI));
B1z = cell(numel(RHO), numel(THETA), numel(PHI));

for ii=1:numel(RHO)
    for jj=1:numel(THETA)
        for kk=1:numel(PHI)
        
        
F_x{ii,jj,kk} = (RHO(ii)>= rho1_max) .* 2/3*M0*mu0 .* sin(theta) .* (RHO(ii) .* (sin(THETA(jj)) .* cos(theta) .* cos(PHI(kk)-phi) - cos(THETA(jj)) .* sin(theta)) ./ ...
                                   (RHO(ii).^2 + rho1_max.^2 - 2.*RHO(ii) .* rho1_max .* (sin(THETA(jj)) .* sin(theta) .* cos(PHI(kk)-phi) + cos(THETA(jj)).* cos(theta))).^3/2) .* rho1_max.^2 .* sin(theta);
                                 
B1x{ii,jj,kk} = squeeze(-trapz(PHI,trapz(THETA,F_x{ii,jj,kk},2))) ; 



F_z{ii,jj,kk} = (RHO(ii)>= rho1_max) .* (2/3*M0*mu0 .* cos(theta) .* sin(phi) .* (RHO(ii) .* (sin(THETA(jj)) .* cos(theta) .* cos(PHI(kk)-phi) - cos(THETA(jj)) .* sin(theta)) ... 
                                  - sin(phi) .* RHO(ii).*sin(THETA(jj)).*sin(PHI(kk)-phi)) ./ ...
                                    (RHO(ii).^2 + rho1_max.^2 - 2.*RHO(ii) .* rho1_max .* (sin(THETA(jj)) .* sin(theta) .* cos(PHI(kk)-phi) + cos(THETA(jj)).* cos(theta))).^3/2) .* rho1_max.^2 .* sin(theta);
                                 
B1z{ii,jj,kk} = (rho < rho1_max).* 2/3*M0*mu0 - squeeze(trapz(PHI,trapz(THETA,F_z{ii,jj,kk},2))) ; 


        end
     end 
end

B1y =                                 (rho>= rho1_max).* 0 ; 

pcolor(squeeze(Xg(:, y==0, :)), squeeze(Zg(:, y==0, :)), squeeze(sqrt(B1x(:,y==0,:).^2 + B1y(:,y==0,:).^2 + B1z(:,y==0,:).^2)))
shading flat
colorbar
step_num = 5;
hold on
quiver(squeeze(Xg(1:step_num:end, y==0, 1:step_num:end)), squeeze(Zg(1:step_num:end, y==0, 1:step_num:end)), squeeze(B1x(1:step_num:end, y==0, 1:step_num:end)), squeeze(B1z(1:step_num:end, y==0, 1:step_num:end)), 'k')
axis equal

[KSSV]

If you have cells...you need to access them using {}. ...your code is a mess....you should re think on it.

[Image Analyst]

Why are you using cell arrays rather than double arrays?

[Nima]

in case of double array, there are problems with calculations of B1x and B1z!

Answer 1

[Image Analyst]

Bly is a double, and the other two are cell arrays. Try this:

y0Index = find(y==0)
xp = B1x{:,y0Index,:} % braces
yp = B1y(:,y0Index,:) % parentheses
zp = B1z{:,y0Index,:} % braces
xyz = sqrt(xp.^2 + yp.^2 + zp.^2)

Then use that in pcolor. You should pack so many things together into a single line. It makes the code easier to maintain if you break it down into smaller chunks.

[Image Analyst]

We have no idea what all this code does (not enough comments). Maybe pcolor is appropriate, maybe not. Maybe you should just call plot(xp) and imshow(squeeze(yp)). Not sure what zp is or how you want to visualize it. Why did you (at first) think that the square root of the sum of the squares of those things would be a 2-d matrix that could be visualized as an image with pcolor? Why don't you just think through this a lot more and figure out what you have and what needs to be visualized and how to do it? Because we're lost in your code. I have no idea what this is all about -- you didn't post the use case or any background context whatsoever. Good luck though.