From 1286e46d21e8d570b11ea9c19956369896036575 Mon Sep 17 00:00:00 2001
From: Isaac Lenton <isaac@isuniversal.com>
Date: Wed, 20 Feb 2019 14:09:57 +1000
Subject: [PATCH] Added a load_stl_shape example

---
 +ott/+shapes/TriangularMesh.m | 33 +++++++++++--
 examples/load_stl_shape.m     | 88 +++++++++++++++++++++++++++++++++++
 2 files changed, 116 insertions(+), 5 deletions(-)
 create mode 100644 examples/load_stl_shape.m

diff --git a/+ott/+shapes/TriangularMesh.m b/+ott/+shapes/TriangularMesh.m
index d69aed10..8e33a831 100644
--- a/+ott/+shapes/TriangularMesh.m
+++ b/+ott/+shapes/TriangularMesh.m
@@ -121,17 +121,40 @@
     function surf(shape, varargin)
       % SURF generate a visualisation of the shape
       %
-      % SURF() displays a visualisation of the shape in the current figure.
+      % SURF(...) displays a visualisation of the shape in the current figure.
       %
-      % SURF(..., 'surfoptions', {varargin}) specifies the options to
-      % pass to the surf function.
+      % Optional named arguments:
+      %   offset   [x;y;z]   offset for location of surface
+      %   rotation   mat     rotation matrix to apply to surface
+      %   surfoptions   {varargin} options to be passed to surf.
 
       p = inputParser;
       p.addParameter('surfoptions', {});
+      p.addParameter('position', []);
+      p.addParameter('rotation', []);
       p.parse(varargin{:});
+      
+      X = shape.verts(1, :);
+      Y = shape.verts(2, :);
+      Z = shape.verts(3, :);
+      
+      % Apply a rotation to the surface
+      if ~isempty(p.Results.rotation)
+        XYZ = [X; Y; Z];
+        XYZ = p.Results.rotation * XYZ;
+        X = XYZ(1, :);
+        Y = XYZ(2, :);
+        Z = XYZ(3, :);
+      end
+
+      % Apply a translation to the surface
+      if ~isempty(p.Results.position)
+        X = X + p.Results.position(1);
+        Y = Y + p.Results.position(2);
+        Z = Z + p.Results.position(3);
+      end
 
-      trisurf(shape.faces.', shape.verts(1, :), ...
-          shape.verts(2, :), shape.verts(3, :), p.Results.surfoptions{:});
+      trisurf(shape.faces.', X, Y, Z, p.Results.surfoptions{:});
     end
 
     function writeWavefrontObj(shape, filename)
diff --git a/examples/load_stl_shape.m b/examples/load_stl_shape.m
new file mode 100644
index 00000000..1ff8e796
--- /dev/null
+++ b/examples/load_stl_shape.m
@@ -0,0 +1,88 @@
+% Load a STL file, generate a voxel grid and scatter a plane wave off the shape
+
+% Add ott to the path
+addpath('../');
+
+% Setup the figure for output
+figure();
+
+%% Load the STL object
+
+% STL reader only supports binary STL files for now
+shape = ott.shapes.StlLoader('mesh.stl');
+
+% Visualise the shape
+h = subplot(1, 3, 1);
+shape.surf('surfoptions', {'EdgeColor', 'none'});
+axis equal;
+view([-30, 25]);
+title('Mesh');
+xlabel('X');
+ylabel('Y');
+zlabel('Z');
+
+%% Convert the shape to a voxel grid
+
+spacing = 1.0;  % Adjust this for your mesh
+
+subplot(1, 3, 2);
+xyz = shape.voxels(spacing, 'visualise', true);
+axis equal;
+view([-30, 25]);
+title('Voxels');
+xlabel('X');
+ylabel('Y');
+zlabel('Z');
+
+%% Calculate the T-matrix for this shape
+
+scale = 0.1;   % Convert from mesh units to simulation units
+
+Tmatrix = ott.TmatrixDda(scale*xyz, 'index_particle', 1.4, 'index_medium', 1.33, ...
+  'wavelength0', 1.0);
+
+%% Calculate torque as a function of axial angle
+
+angles = linspace(0, 2*pi, 100);
+rotation = zeros(3, 3*length(angles));
+for ii = 1:length(angles)
+  rotation(:, (1:3) + (ii-1)*3) = roty(angles(ii)*180/pi);
+end
+
+beam = ott.BscPmGauss('NA', 1.02, 'power', 1.0, 'index_medium', 1.33, ...
+  'polarisation', [1, 1i], 'wavelength0', 1.0);
+
+[~, torque] = ott.forcetorque(beam, Tmatrix, 'rotation', rotation);
+
+subplot(1, 3, 3);
+plot(angles, torque);
+xlabel('Y-Angle [rad]');
+ylabel('Torque [a.u.]');
+legend({'X', 'Y', 'Z'}, 'Location', 'SouthEast');
+title('Torque');
+
+%% Animate the mesh rotating (no translation)
+
+beam = ott.BscPmGauss('NA', 1.02, 'power', 1.0, 'index_medium', 1.33, ...
+  'polarisation', [1, 1i], 'wavelength0', 1.0);
+
+x = [0;0;0];
+Rx = eye(3);
+dtheta = 2*pi/100;
+dx = 0.0;  % no translation
+
+figure();
+shape.surf('rotation', Rx, 'surfoptions', {'EdgeColor', 'none'});
+axis equal;
+axis([-10, 10, -10, 10, -2, 10]);
+
+for ii = 1:100
+  [f, t] = ott.forcetorque(beam, Tmatrix, 'rotation', Rx, 'position', x);
+  x = x + dx*f./vecnorm(f);
+  Rx = ott.utils.rotation_matrix(dtheta*t./vecnorm(t))*Rx;
+  
+  shape.surf('rotation', Rx, 'surfoptions', {'EdgeColor', 'none'});
+  axis equal;
+  axis([-10, 10, -10, 10, -2, 10]);
+  drawnow;
+end