error about reshape when I change the number of spheres

[ShaoChangk]

I have modified the examples and change the distribution of spheres. I am prompted to use reshape incorrectly.
And my txt of spheres was shown below.
1.txt

And this is my main code.
`% =========================================================================
%> @brief Run this script to start the simulation
% =========================================================================

%% add all folders to the MATLAB search path
addpath(genpath('./src'))

%% import example file with sphere positions, radii and refractive indices
data = dlmread('examples/1.txt');

%% initialize the CELES class instances

% initialize particle class instance
% - positionArray: Nx3 array (float) in [x,y,z] format
% - refractiveIndexArray: Nx1 array (complex) of complex refractive indices
% - radiusArray: Nx1 array (float) of sphere radii
particles = celes_particles('positionArray', data(:,1:3), ...
'refractiveIndexArray', data(:,5)+1i*data(:,6), ...
'radiusArray', data(:,4) ...
);

% initialize initial field class instance
% - polarAngle: scalar (float) polar angle of incoming beam/wave,
% in radians. for Gaussian beams, only 0 or pi are
% currently possible
% - azimuthalAngle: scalar (float) azimuthal angle of incoming
% beam/wave, in radians
% - polarization: string (char) polarization of incoming beam/wave
% ('TE' or 'TM')
% - beamWidth: scalar (float) width of beam waist. use 0 or inf
% for plane wave
% - focalPoint: 1x3 array (float) focal point
initialField = celes_initialField('polarAngle', 0, ...
'azimuthalAngle', 0, ...
'polarization', 'TE', ...
'beamWidth', 2000, ...
'focalPoint', [0,0,0] ...
);

% initialize input class instance
% - wavelength: scalar (float) vacuum wavelength, same unit as
% particle positions and radii
% - mediumRefractiveIndex: scalar (complex) refractive index of environment
% - particles: valid instance of celes_particles class
% - initialField: valid instance of celes_initialField class
input = celes_input('wavelength', 550, ...
'mediumRefractiveIndex', 1.49, ...
'particles', particles, ...
'initialField', initialField ...
);

% initialize preconditioner class instance
% - type: string (char) type of preconditioner (currently
% only 'blockdiagonal' and 'none' available)
% - partitionEdgeSizes 1x3 array (float) edge size of partitioning cuboids
% (applies to 'blockdiagonal' type only)
precnd = celes_preconditioner('type', 'blockdiagonal', ...
'partitionEdgeSizes', [1200,1200,1200] ...
);

% initialize solver class instance
% - type: string (char) solver type (currently 'BiCGStab' or
% 'GMRES' are implemented)
% - tolerance: scalar (float) target relative accuracy of solution
% - maxIter: scalar (int) maximum number of iterations allowed
% - restart: scalar (int) restart parameter (applies only to
% GMRES solver)
% - preconditioner: valid instance of celes_preconditioner class
solver = celes_solver('type', 'GMRES', ...
'tolerance', 5e-4, ...
'maxIter', 1000, ...
'restart', 200, ...
'preconditioner', precnd);

% initialize numerics class instance
% - lmax: scalar (int) maximal expansion order of scattered
% fields around particle center
% - polarAnglesArray: 1xN array (float) sampling of polar angles in the
% plane wave patterns, in radians
% - azimuthalAnglesArray: sampling of azimuthal angles in the plane wave
% patterns, in radians
% - gpuFlag: scalar (bool) set to false if you experience GPU
% memory problems at evaluation time (translation
% operator always runs on GPU, though)
% - particleDistanceResolution: scalar (float) resolution of lookup table for
% spherical Hankel function (same unit as wavelength)
% - solver: valid instance of celes_solver class
numerics = celes_numerics('lmax', 3, ...
'polarAnglesArray', 0:pi/5e3:pi, ...
'azimuthalAnglesArray', 0:pi/1e2:2*pi, ...
'gpuFlag', true, ...
'particleDistanceResolution', 1, ...
'solver', solver);

% define a grid of points where the field will be evaulated
[x,z] = meshgrid(-5000:50:5000, -5000:50:5000); y = zeros(size(x));
% initialize output class instance
% - fieldPoints: Nx3 array (float) points where to evaluate the
% electric near field
% - fieldPointsArrayDims: 1x2 array (int) dimensions of the array, needed to
% recompose the computed field as a n-by-m image
output = celes_output('fieldPoints', [x(:),y(:),z(:)], ...
'fieldPointsArrayDims', size(x));

% initialize simulation class instance
% - input: valid instance of celes_input class
% - numerics: valid instance of celes_input class
% - output: valid instance of celes_output class
simul = celes_simulation('input', input, ...
'numerics', numerics, ...
'output', output);

%% run simulation
simul.run;

% evaluate transmitted and reflected power
simul.evaluatePower;
fprintf('transmitted power: \t%.4f %%\n', ...
simul.output.totalFieldForwardPower/simul.output.initialFieldPower100)
fprintf('reflected power: \t%.4f %%\n', ...
simul.output.totalFieldBackwardPower/simul.output.initialFieldPower100)

% evaluate field at output.fieldPoints
simul.evaluateFields;

%% plot results
% display particles
figure('Name','Particle positions','NumberTitle','off');
plot_spheres(gca,simul.input.particles.positionArray, ...
simul.input.particles.radiusArray, ...
simul.input.particles.refractiveIndexArray)

% plot near field
figure('Name','Near-field cross-cut','NumberTitle','off');
plot_field(gca,simul,'abs E','Total field')
colorbar
caxis([0,2])

% % export animated gif
% figure('Name','Animated near-field cross-cut','NumberTitle','off');
% plot_field(gca,simul,'real Ey','Total field','Ey_total.gif')
`

[ShaoChangk]

And this is the code of plot field
`function plot_field(ax,simulation,component,fieldType,GIFoutputname)

hold(ax,'on')

pArr = simulation.input.particles.positionArray;
rArr = simulation.input.particles.radiusArray;
dims = simulation.output.fieldPointsArrayDims;

switch lower(fieldType)
case 'total field'
E = simulation.output.totalField;
case 'scattered field'
E = simulation.output.scatteredField + simulation.output.internalField;
case 'initial field'
E = simulation.output.initialField;
end

switch lower(component)
case 'real ex'
fld = reshape(gather(E(:,1)), dims);
case 'real ey'
fld = reshape(gather(E(:,2)), dims);
case 'real ez'
fld = reshape(gather(E(:,3)), dims);
case 'abs e'
fld = reshape(gather(sqrt(sum(abs(E).^2,2))), dims);
end

if exist('GIFoutputname','var')
t = linspace(0,2*pi,26); t(end) = []; % 25 frames
cmap = interp1(1:3,[0 0 1; 1 1 1; 1 0 0],linspace(1,3,256-32));
gifcmap = [cmap; gray(32)]; % add some grayscale colors for the gif colormap
caxislim = [-max(abs(fld(:))), max(abs(fld(:)))];
else
t = 0;
cmap = parula(256);
caxislim = [0, max(abs(fld(:)))];
end

fldPnts = reshape([simulation.output.fieldPoints(:,1), ...
simulation.output.fieldPoints(:,2), ...
simulation.output.fieldPoints(:,3)], [dims,3]);

if all(fldPnts(:,:,1) == fldPnts(1,1,1)) % fldPoints are on the yz plane
perpdim = 1; % 1->x is the perp. direction
elseif all(fldPnts(:,:,2) == fldPnts(1,1,2))% fldPoints are on the xz plane
perpdim = 2; % 2->y is the perp. direction
elseif all(fldPnts(:,:,3) == fldPnts(1,1,3))% fldPoints are on the xy plane
perpdim = 3; % 3->z is the perp. direction
else
error('fieldPoint must define an xy, xz or yz-like plane')
end

xy = setdiff([1,2,3], perpdim); % here xy are the in-plane dimensions
x = fldPnts(:,:,xy(1));
y = fldPnts(:,:,xy(2));

dist = abs(pArr(:,perpdim) - fldPnts(1,1,perpdim));% particle distances from xy plane
idx = find(dist<rArr); % find particles intersecting the plane
rArr(idx) = sqrt(rArr(idx).^2 - dist(idx).^2); % overwrite radius of the intersection circle

if exist('GIFoutputname','var') % initialize imind array
f = getframe(gcf);
imind = rgb2ind(f.cdata,gifcmap,'nodither');
imind(1,1,1,numel(t)) = 0;
end

for ti=1:numel(t)
imagesc(x(1,:), y(:,1), real(fldexp(-1it(ti))))% plot field on a xy plane
colormap(cmap)
for i=1:length(idx)
rectangle(ax, ...
'Position', [pArr(idx(i),xy)-rArr(idx(i)), [2,2]*rArr(idx(i))], ...
'Curvature', [1 1], ...
'FaceColor', 'none', ...
'EdgeColor', [0,0,0], ...
'LineWidth', 0.75)
end
axis([min(x(:)),max(x(:)),min(y(:)),max(y(:))]) % set axis tight to fldPoints

labels = ['x'; 'y'; 'z'];
xlabel(labels(xy(1)))
ylabel(labels(xy(2)))

ax.DataAspectRatio = [1,1,1];
title([fieldType,', ',component])
caxis(caxislim)
colorbar
drawnow

if exist('GIFoutputname','var')
    f = getframe(gcf);
    imind(:,:,1,ti) = rgb2ind(f.cdata,gifcmap,'nodither');
end

end

hold(ax,'off')

if exist('GIFoutputname','var')
imwrite(imind,gifcmap,GIFoutputname,'DelayTime',0,'Loopcount',inf)
end

end
`

[lpattelli]

Hi,
the error you get is due to the fact that the simulation did not return a valid field array (simul.output.totalField is an empty array which cannot be reshaped into its intended shape). The reason why you get an invalid output is that your particle configuration comprises overlapping spheres, which is not allowed in the T-matrix approach used by CELES. Your script runs fine if you reduce the particle radii up to a point where they do not intersect each other.

Perhaps we should think about throwing an error if the user provides an input file with overlapping spheres, but that sounds computationally expensive as one should loop through several pairs of particles to check that no sphere is overlapping with any other sphere.

[ShaoChangk]

Thank you every much!!!

…

------------------ 原始邮件 ------------------ 发件人: "Lorenzo Pattelli"<notifications@github.com>; 发送时间: 2019年9月17日(星期二) 下午3:47 收件人: "disordered-photonics/celes"<celes@noreply.github.com>; 抄送: "邵常焜"<1224950293@qq.com>; "Manual"<manual@noreply.github.com>; 主题: Re: [disordered-photonics/celes] error about reshape when I changethe number of spheres (#25) Hi, the error you get is due to the fact that the simulation did not return a valid field array (simul.output.totalField is an empty array which cannot be reshaped into its intended shape). The reason why you get an invalid output is that your particle configuration comprises overlapping spheres, which is not allowed in the T-matrix approach used by CELES. Your script runs fine if you reduce the particle radii up to a point where they do not intersect each other. Perhaps we should think about throwing an error if the user provides an input file with overlapping spheres, but that sounds computationally expensive as one should loop through several pairs of particles to check that no sphere is overlapping with any other sphere. — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub, or mute the thread.

[ShaoChangk]

Hi,
the error you get is due to the fact that the simulation did not return a valid field array (simul.output.totalField is an empty array which cannot be reshaped into its intended shape). The reason why you get an invalid output is that your particle configuration comprises overlapping spheres, which is not allowed in the T-matrix approach used by CELES. Your script runs fine if you reduce the particle radii up to a point where they do not intersect each other.
Perhaps we should think about throwing an error if the user provides an input file with overlapping spheres, but that sounds computationally expensive as one should loop through several pairs of particles to check that no sphere is overlapping with any other sphere.

thank you very much!!!