Merge pull request #26 from disordered-photonics/Hnf

compute H near field, update plotting functions

CELES_MAIN.m

@@ -14,9 +14,9 @@
 %   - 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), ...
+particles = celes_particles('positionArray',        data(:,2:4), ...
                             'refractiveIndexArray', data(:,5)+1i*data(:,6), ...
-                            'radiusArray',          data(:,4) ...
+                            'radiusArray',          data(:,1) ...
                             );
 
 % initialize initial field class instance
@@ -93,7 +93,7 @@
                           'solver',                 solver);
 
 % define a grid of points where the field will be evaulated
-[x,z] = meshgrid(-4000:50:4000, -3000:50:5000); y = zeros(size(x));
+[x,z] = meshgrid(-4000:40:4000, -3000:40:5000); y = zeros(size(x));
 % initialize output class instance
 %   - fieldPoints:          Nx3 array (float) points where to evaluate the
 %                           electric near field
@@ -131,9 +131,21 @@
                  simul.input.particles.refractiveIndexArray)
 
 % plot near field
+comp = ["real Ex", "real Ey", "real Ez", "abs E", "real Hx", "real Hy", "real Hz", "abs H"];
+fieldtype = 'total field';
+
 figure('Name','Near-field cross-cut','NumberTitle','off');
-plot_field(gca,simul,'abs E','Total field')
-caxis([0,2])
+for sp = 1:numel(comp)
+    subplot(2,4,sp)
+    plot_field(gca, simul, comp(sp), fieldtype);
+    caxis([-2*contains(comp(sp),'real'), 2])
+end
+linkaxes(findall(gcf,'type','axes'))
+
+% plot Poynting vector
+figure('Name','Pynting vector cross-cut','NumberTitle','off');
+S = plot_poynting(gca,simul,'Total field');
+ylim([min(z(:)), max(z(:))])
 
 % % export animated gif
 % figure('Name','Animated near-field cross-cut','NumberTitle','off');

examples/.gitignore

@@ -1,3 +1,2 @@
 /*
 !.gitignore
-!example_project/

examples/mstm_benchmark.inp

@@ -0,0 +1,60 @@
+begin_comment
+***********************************************************
+MSTM input file to replicate CELES_MAIN's near-field output
+***********************************************************
+end_comment
+number_spheres
+500
+sphere_position_file
+sphere_parameters.txt
+output_file
+mstm_benchmark_TE.dat
+write_sphere_data
+1
+length_scale_factor
+0.0114239733
+real_ref_index_scale_factor
+1.0d0
+imag_ref_index_scale_factor
+1.0d0
+mie_epsilon
+-3
+incident_or_target_frame
+0
+normalize_scattering_matrix
+1
+gaussian_beam_constant
+0.0437676d0
+gaussian_beam_focal_point
+0.0d0,0.0d0,0.0d0
+fixed_or_random_orientation
+0
+incident_azimuth_angle_deg
+0.000000d0
+incident_polar_angle_deg
+0.000000d0
+near_field_distance
+1d8
+calculate_scattering_coefficients
+1
+calculate_near_field
+1
+near_field_output_data
+2
+near_field_plane_coord
+2
+near_field_plane_position
+0.0d0
+near_field_plane_vertices
+-34.271920d0,-45.695893d0,57.119866d0,45.695893d0
+spacial_step_size
+0.45695893d0
+polarization_angle_deg
+90.d0
+near_field_output_file
+mstm_benchmark_TE_nf.dat
+near_field_output_data
+2
+track_iterations
+1
+end_of_options

examples/plot_mstm_nf.m

@@ -0,0 +1,146 @@
+% plot MSTM output
+clear; close all
+
+if exist('mstm_benchmark_nf.mat', 'file') == 2
+    load('mstm_benchmark_nf.mat');
+else
+    % get the number of lines to skip (# spheres cut by the near-field plane)
+    filename = 'mstm_benchmark_0deg_nf.dat';
+    fileID = fopen(filename,'r');
+    formatSpec = '%f%[^\n\r]';
+    dataArray = textscan(fileID, formatSpec, 1, 'HeaderLines', 1);
+    cutspheres = dataArray{:, 1};
+    fclose(fileID);
+    clearvars formatSpec fileID dataArray ans;
+
+    % read near-field data for 0 deg polarization
+    mstmOutput = dlmread(filename,'',cutspheres+2,0);
+    xvec=unique(mstmOutput(:,1));
+    zvec=unique(mstmOutput(:,2));
+    E0 = cat(3, rot90(reshape(mstmOutput(:,3) + 1i*mstmOutput(:,4), [length(zvec), length(xvec)])), ...
+                rot90(reshape(mstmOutput(:,5) + 1i*mstmOutput(:,6), [length(zvec), length(xvec)])), ...
+                rot90(reshape(mstmOutput(:,7) + 1i*mstmOutput(:,8), [length(zvec), length(xvec)])));
+    H0 = cat(3, rot90(reshape(mstmOutput(:,9) + 1i*mstmOutput(:,10),[length(zvec), length(xvec)])), ...
+                rot90(reshape(mstmOutput(:,11)+ 1i*mstmOutput(:,12),[length(zvec), length(xvec)])), ...
+                rot90(reshape(mstmOutput(:,13)+ 1i*mstmOutput(:,14),[length(zvec), length(xvec)])));
+    E0 = single(E0); % the ascii output has only a few digits...
+    H0 = single(H0);
+
+    % read near-field data for 90 deg polarization
+    filename = 'mstm_benchmark_90deg_nf.dat';
+    mstmOutput = dlmread(filename,'',cutspheres+2,0);
+    E90= cat(3, rot90(reshape(mstmOutput(:,3) + 1i*mstmOutput(:,4), [length(zvec), length(xvec)])), ...
+                rot90(reshape(mstmOutput(:,5) + 1i*mstmOutput(:,6), [length(zvec), length(xvec)])), ...
+                rot90(reshape(mstmOutput(:,7) + 1i*mstmOutput(:,8), [length(zvec), length(xvec)])));
+    H90= cat(3, rot90(reshape(mstmOutput(:,9) + 1i*mstmOutput(:,10),[length(zvec), length(xvec)])), ...
+                rot90(reshape(mstmOutput(:,11)+ 1i*mstmOutput(:,12),[length(zvec), length(xvec)])), ...
+                rot90(reshape(mstmOutput(:,13)+ 1i*mstmOutput(:,14),[length(zvec), length(xvec)])));
+    E90 = single(E90);
+    H90 = single(H90);
+    [xgrid,zgrid]=meshgrid(xvec,zvec);
+end
+
+%%
+
+cmap = interp1(1:3,[0 0 1; 1 1 1; 1 0 0],linspace(1,3,256)); % define a diverging colormap
+
+figure
+subplot(2,4,4)
+imagesc(xvec, zvec, sqrt(sum(abs(E0).^2, 3))), axis equal tight
+caxis([0,2])
+title('|E|, 0 deg'), xlabel('k*x'), ylabel('k*z')
+
+subplot(2,4,1)
+imagesc(xvec, zvec, real(E0(:,:,1))), axis equal tight
+caxis([-2,2])
+colormap(gca,cmap)
+title('E_x, 0 deg'), xlabel('k*x'), ylabel('k*z')
+
+subplot(2,4,2)
+imagesc(xvec, zvec, real(E0(:,:,2))), axis equal tight
+caxis([-2,2])
+colormap(gca,cmap)
+title('E_y, 0 deg'), xlabel('k*x'), ylabel('k*z')
+
+subplot(2,4,3)
+imagesc(xvec, zvec, real(E0(:,:,3))), axis equal tight
+caxis([-2,2])
+colormap(gca,cmap)
+title('E_z, 0 deg'), xlabel('k*x'), ylabel('k*z')
+
+subplot(2,4,8)
+imagesc(xvec, zvec, sqrt(sum(abs(H0).^2, 3))), axis equal tight
+caxis([0,2])
+title('|H|, 0 deg'), xlabel('k*x'), ylabel('k*z')
+
+subplot(2,4,5)
+imagesc(xvec, zvec, real(H0(:,:,1))), axis equal tight
+caxis([-2,2])
+colormap(gca,cmap)
+title('H_x, 0 deg'), xlabel('k*x'), ylabel('k*z')
+
+subplot(2,4,6)
+imagesc(xvec, zvec, real(H0(:,:,2))), axis equal tight
+caxis([-2,2])
+colormap(gca,cmap)
+title('H_y, 0 deg'), xlabel('k*x'), ylabel('k*z')
+
+subplot(2,4,7)
+imagesc(xvec, zvec, real(H0(:,:,3))), axis equal tight
+caxis([-2,2])
+colormap(gca,cmap)
+title('H_z, 0 deg'), xlabel('k*x'), ylabel('k*z')
+linkaxes(findall(gcf,'type','axes'))
+
+figure
+subplot(2,4,4)
+imagesc(xvec, zvec, sqrt(sum(abs(E90).^2, 3))), axis equal tight
+caxis([0,2])
+title('|E|, 90 deg'), xlabel('k*x'), ylabel('k*z')
+
+subplot(2,4,1)
+imagesc(xvec, zvec, real(E90(:,:,1))), axis equal tight
+caxis([-2,2])
+colormap(gca,cmap)
+title('E_x, 90 deg'), xlabel('k*x'), ylabel('k*z')
+
+subplot(2,4,2)
+imagesc(xvec, zvec, real(E90(:,:,2))), axis equal tight
+caxis([-2,2])
+colormap(gca,cmap)
+title('E_y, 90 deg'), xlabel('k*x'), ylabel('k*z')
+
+subplot(2,4,3)
+imagesc(xvec, zvec, real(E90(:,:,3))), axis equal tight
+caxis([-2,2])
+colormap(gca,cmap)
+title('E_z, 90 deg'), xlabel('k*x'), ylabel('k*z')
+
+subplot(2,4,8)
+imagesc(xvec, zvec, sqrt(sum(abs(H90).^2, 3))), axis equal tight
+caxis([0,2])
+title('|H|, 90 deg'), xlabel('k*x'), ylabel('k*z')
+
+subplot(2,4,5)
+imagesc(xvec, zvec, real(H90(:,:,1))), axis equal tight
+caxis([-2,2])
+colormap(gca,cmap)
+title('H_x, 90 deg'), xlabel('k*x'), ylabel('k*z')
+
+subplot(2,4,6)
+imagesc(xvec, zvec, real(H90(:,:,2))), axis equal tight
+caxis([-2,2])
+colormap(gca,cmap)
+title('H_y, 90 deg'), xlabel('k*x'), ylabel('k*z')
+
+subplot(2,4,7)
+imagesc(xvec, zvec, real(H90(:,:,3))), axis equal tight
+caxis([-2,2])
+colormap(gca,cmap)
+title('H_z, 90 deg'), xlabel('k*x'), ylabel('k*z')
+linkaxes(findall(gcf,'type','axes'))
+
+clearvars cutspheres filename fld_mstm mstmOutput;
+if ~(exist('mstm_benchmark_nf.mat', 'file') == 2)
+    save('mstm_benchmark_nf.mat')
+end