ShaoChangk
commented
5 years ago 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');
endend
hold(ax,'off')
if exist('GIFoutputname','var') imwrite(imind,gifcmap,GIFoutputname,'DelayTime',0,'Loopcount',inf) end
end `
lpattelli
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') `