diff --git a/ProgramFiles/v4.2/Utilities/GetFullPath.m b/ProgramFiles/v4.2/Utilities/GetFullPath.m new file mode 100755 index 0000000..79788af --- /dev/null +++ b/ProgramFiles/v4.2/Utilities/GetFullPath.m @@ -0,0 +1,326 @@ +function File = GetFullPath(File, Style) +% GetFullPath - Get absolute canonical path of a file or folder +% Absolute path names are safer than relative paths, when e.g. a GUI or TIMER +% callback changes the current directory. Only canonical paths without "." and +% ".." can be recognized uniquely. +% Long path names (>259 characters) require a magic initial key "\\?\" to be +% handled by Windows API functions, e.g. for Matlab's FOPEN, DIR and EXIST. +% +% FullName = GetFullPath(Name, Style) +% INPUT: +% Name: String or cell string, absolute or relative name of a file or +% folder. The path need not exist. Unicode strings, UNC paths and long +% names are supported. +% Style: Style of the output as string, optional, default: 'auto'. +% 'auto': Add '\\?\' or '\\?\UNC\' for long names on demand. +% 'lean': Magic string is not added. +% 'fat': Magic string is added for short names also. +% The Style is ignored when not running under Windows. +% +% OUTPUT: +% FullName: Absolute canonical path name as string or cell string. +% For empty strings the current directory is replied. +% '\\?\' or '\\?\UNC' is added on demand. +% +% NOTE: The M- and the MEX-version create the same results, the faster MEX +% function works under Windows only. +% Some functions of the Windows-API still do not support long file names. +% E.g. the Recycler and the Windows Explorer fail even with the magic '\\?\' +% prefix. Some functions of Matlab accept 260 characters (value of MAX_PATH), +% some at 259 already. Don't blame me. +% The 'fat' style is useful e.g. when Matlab's DIR command is called for a +% folder with les than 260 characters, but together with the file name this +% limit is exceeded. Then "dir(GetFullPath([folder, '\*.*], 'fat'))" helps. +% +% EXAMPLES: +% cd(tempdir); % Assumed as 'C:\Temp' here +% GetFullPath('File.Ext') % 'C:\Temp\File.Ext' +% GetFullPath('..\File.Ext') % 'C:\File.Ext' +% GetFullPath('..\..\File.Ext') % 'C:\File.Ext' +% GetFullPath('.\File.Ext') % 'C:\Temp\File.Ext' +% GetFullPath('*.txt') % 'C:\Temp\*.txt' +% GetFullPath('..') % 'C:\' +% GetFullPath('..\..\..') % 'C:\' +% GetFullPath('Folder\') % 'C:\Temp\Folder\' +% GetFullPath('D:\A\..\B') % 'D:\B' +% GetFullPath('\\Server\Folder\Sub\..\File.ext') +% % '\\Server\Folder\File.ext' +% GetFullPath({'..', 'new'}) % {'C:\', 'C:\Temp\new'} +% GetFullPath('.', 'fat') % '\\?\C:\Temp\File.Ext' +% +% COMPILE: +% Automatic: InstallMex GetFullPath.c uTest_GetFullPath +% Manual: mex -O GetFullPath.c +% Download: http://www.n-simon.de/mex +% Run the unit-test uTest_GetFullPath after compiling. +% +% Tested: Matlab 6.5, 7.7, 7.8, 7.13, WinXP/32, Win7/64 +% Compiler: LCC2.4/3.8, BCC5.5, OWC1.8, MSVC2008/2010 +% Assumed Compatibility: higher Matlab versions +% Author: Jan Simon, Heidelberg, (C) 2009-2013 matlab.THISYEAR(a)nMINUSsimon.de +% +% See also: CD, FULLFILE, FILEPARTS. + +% $JRev: R-G V:032 Sum:7Xd/JS0+yfax Date:15-Jan-2013 01:06:12 $ +% $License: BSD (use/copy/change/redistribute on own risk, mention the author) $ +% $UnitTest: uTest_GetFullPath $ +% $File: Tools\GLFile\GetFullPath.m $ +% History: +% 001: 20-Apr-2010 22:28, Successor of Rel2AbsPath. +% 010: 27-Jul-2008 21:59, Consider leading separator in M-version also. +% 011: 24-Jan-2011 12:11, Cell strings, '~File' under linux. +% Check of input types in the M-version. +% 015: 31-Mar-2011 10:48, BUGFIX: Accept [] as input as in the Mex version. +% Thanks to Jiro Doke, who found this bug by running the test function for +% the M-version. +% 020: 18-Oct-2011 00:57, BUGFIX: Linux version created bad results. +% Thanks to Daniel. +% 024: 10-Dec-2011 14:00, Care for long names under Windows in M-version. +% Improved the unittest function for Linux. Thanks to Paul Sexton. +% 025: 09-Aug-2012 14:00, In MEX: Paths starting with "\\" can be non-UNC. +% The former version treated "\\?\C:\\file" as UNC path and +% replied "\\?\UNC\?\C:\\file". +% 032: 12-Jan-2013 21:16, 'auto', 'lean' and 'fat' style. + +% Initialize: ================================================================== +% Do the work: ================================================================= + +% ############################################# +% ### USE THE MUCH FASTER MEX ON WINDOWS!!! ### +% ############################################# + +% Difference between M- and Mex-version: +% - Mex does not work under MacOS/Unix. +% - Mex calls Windows API function GetFullPath. +% - Mex is much faster. + +% Magix prefix for long Windows names: +if nargin < 2 + Style = 'auto'; +end + +% Handle cell strings: +% NOTE: It is faster to create a function @cell\GetFullPath.m under Linux, but +% under Windows this would shadow the fast C-Mex. +if isa(File, 'cell') + for iC = 1:numel(File) + File{iC} = GetFullPath(File{iC}, Style); + end + return; +end + +% Check this once only: +isWIN = strncmpi(computer, 'PC', 2); +MAX_PATH = 260; + +% Warn once per session (disable this under Linux/MacOS): +persistent hasDataRead +if isempty(hasDataRead) + % Test this once only - there is no relation to the existence of DATAREAD! + %if isWIN + % Show a warning, if the slower Matlab version is used - commented, because + % this is not a problem and it might be even useful when the MEX-folder is + % not inlcuded in the path yet. + % warning('JSimon:GetFullPath:NoMex', ... + % ['GetFullPath: Using slow Matlab-version instead of fast Mex.', ... + % char(10), 'Compile: InstallMex GetFullPath.c']); + %end + + % DATAREAD is deprecated in 2011b, but still available. In Matlab 6.5, REGEXP + % does not know the 'split' command, therefore DATAREAD is preferred: + hasDataRead = ~isempty(which('dataread')); +end + +if isempty(File) % Accept empty matrix as input: + if ischar(File) || isnumeric(File) + File = cd; + return; + else + error(['JSimon:', mfilename, ':BadTypeInput1'], ... + ['*** ', mfilename, ': Input must be a string or cell string']); + end +end + +if ischar(File) == 0 % Non-empty inputs must be strings + error(['JSimon:', mfilename, ':BadTypeInput1'], ... + ['*** ', mfilename, ': Input must be a string or cell string']); +end + +if isWIN % Windows: -------------------------------------------------------- + FSep = '\'; + File = strrep(File, '/', FSep); + + % Remove the magic key on demand, it is appended finally again: + if strncmp(File, '\\?\', 4) + if strncmpi(File, '\\?\UNC\', 8) + File = ['\', File(7:length(File))]; % Two leading backslashes! + else + File = File(5:length(File)); + end + end + + isUNC = strncmp(File, '\\', 2); + FileLen = length(File); + if isUNC == 0 % File is not a UNC path + % Leading file separator means relative to current drive or base folder: + ThePath = cd; + if File(1) == FSep + if strncmp(ThePath, '\\', 2) % Current directory is a UNC path + sepInd = strfind(ThePath, '\'); + ThePath = ThePath(1:sepInd(4)); + else + ThePath = ThePath(1:3); % Drive letter only + end + end + + if FileLen < 2 || File(2) ~= ':' % Does not start with drive letter + if ThePath(length(ThePath)) ~= FSep + if File(1) ~= FSep + File = [ThePath, FSep, File]; + else % File starts with separator: + File = [ThePath, File]; + end + else % Current path ends with separator: + if File(1) ~= FSep + File = [ThePath, File]; + else % File starts with separator: + ThePath(length(ThePath)) = []; + File = [ThePath, File]; + end + end + + elseif FileLen == 2 && File(2) == ':' % "C:" current directory on C! + % "C:" is the current directory on the C-disk, even if the current + % directory is on another disk! This was ignored in Matlab 6.5, but + % modern versions considers this strange behaviour. + if strncmpi(ThePath, File, 2) + File = ThePath; + else + try + File = cd(cd(File)); + catch % No MException to support Matlab6.5... + if exist(File, 'dir') % No idea what could cause an error then! + rethrow(lasterror); + else % Reply "K:\" for not existing disk: + File = [File, FSep]; + end + end + end + end + end + +else % Linux, MacOS: --------------------------------------------------- + FSep = '/'; + File = strrep(File, '\', FSep); + + if strcmp(File, '~') || strncmp(File, '~/', 2) % Home directory: + HomeDir = getenv('HOME'); + if ~isempty(HomeDir) + File(1) = []; + File = [HomeDir, File]; + end + + elseif strncmpi(File, FSep, 1) == 0 + % Append relative path to current folder: + ThePath = cd; + if ThePath(length(ThePath)) == FSep + File = [ThePath, File]; + else + File = [ThePath, FSep, File]; + end + end +end + +% Care for "\." and "\.." - no efficient algorithm, but the fast Mex is +% recommended at all! +if ~isempty(strfind(File, [FSep, '.'])) + if isWIN + if strncmp(File, '\\', 2) % UNC path + index = strfind(File, '\'); + if length(index) < 4 % UNC path without separator after the folder: + return; + end + Drive = File(1:index(4)); + File(1:index(4)) = []; + else + Drive = File(1:3); + File(1:3) = []; + end + else % Unix, MacOS: + isUNC = false; + Drive = FSep; + File(1) = []; + end + + hasTrailFSep = (File(length(File)) == FSep); + if hasTrailFSep + File(length(File)) = []; + end + + if hasDataRead + if isWIN % Need "\\" as separator: + C = dataread('string', File, '%s', 'delimiter', '\\'); %#ok + else + C = dataread('string', File, '%s', 'delimiter', FSep); %#ok + end + else % Use the slower REGEXP, when DATAREAD is not available anymore: + C = regexp(File, FSep, 'split'); + end + + % Remove '\.\' directly without side effects: + C(strcmp(C, '.')) = []; + + % Remove '\..' with the parent recursively: + R = 1:length(C); + for dd = reshape(find(strcmp(C, '..')), 1, []) + index = find(R == dd); + R(index) = []; + if index > 1 + R(index - 1) = []; + end + end + + if isempty(R) + File = Drive; + if isUNC && ~hasTrailFSep + File(length(File)) = []; + end + + elseif isWIN + % If you have CStr2String, use the faster: + % File = CStr2String(C(R), FSep, hasTrailFSep); + File = sprintf('%s\\', C{R}); + if hasTrailFSep + File = [Drive, File]; + else + File = [Drive, File(1:length(File) - 1)]; + end + + else % Unix: + File = [Drive, sprintf('%s/', C{R})]; + if ~hasTrailFSep + File(length(File)) = []; + end + end +end + +% "Very" long names under Windows: +if isWIN + if ~ischar(Style) + error(['JSimon:', mfilename, ':BadTypeInput2'], ... + ['*** ', mfilename, ': Input must be a string or cell string']); + end + + if (strncmpi(Style, 'a', 1) && length(File) >= MAX_PATH) || ... + strncmpi(Style, 'f', 1) + % Do not use [isUNC] here, because this concerns the input, which can + % '.\File', while the current directory is an UNC path. + if strncmp(File, '\\', 2) % UNC path + File = ['\\?\UNC', File(2:end)]; + else + File = ['\\?\', File]; + end + end +end + +% return; diff --git a/ProgramFiles/v4.2/Utilities/LowRe_toolbox/LowRE_dissipation_metric_from_curvature_bases.m b/ProgramFiles/v4.2/Utilities/LowRe_toolbox/LowRE_dissipation_metric_from_curvature_bases.m index a70b5f9..8198c1b 100644 --- a/ProgramFiles/v4.2/Utilities/LowRe_toolbox/LowRE_dissipation_metric_from_curvature_bases.m +++ b/ProgramFiles/v4.2/Utilities/LowRe_toolbox/LowRE_dissipation_metric_from_curvature_bases.m @@ -5,7 +5,7 @@ int_limit = L*[-0.5 0.5]; % Define the tangential, lateral drag matrix for unit/double drag - drag = [1 0; 0 2]*c; + drag = [1 0; 0 drag_ratio]*c; % Get the backbone locus, Jacobian, and Local Connection functions [A, h, J] = LowRE_local_connection_from_curvature_bases(kappa_basis_input,r,L,c,drag_ratio); diff --git a/ProgramFiles/v4.2/Utilities/LowRe_toolbox/LowRE_dissipation_metric_from_general_curvature.m b/ProgramFiles/v4.2/Utilities/LowRe_toolbox/LowRE_dissipation_metric_from_general_curvature.m new file mode 100644 index 0000000..130edd3 --- /dev/null +++ b/ProgramFiles/v4.2/Utilities/LowRe_toolbox/LowRE_dissipation_metric_from_general_curvature.m @@ -0,0 +1,42 @@ +function Mp = LowRE_dissipation_metric_from_general_curvature(curvdef,cparams,L,c,drag_ratio) +% Calculate the dissipation power metric for a set of curvature bases + + % Specified integration limits + int_limit = L*[-0.5 0.5]; + + % Define the tangential, lateral drag matrix for unit/double drag + drag = [1 0; 0 drag_ratio]*c; + + % Get the backbone locus, Jacobian, and Local Connection functions + [A, h, J,Omega] = LowRE_local_connection_from_general_curvature(curvdef,cparams,L,c,drag_ratio); + + % Integrate along the body for the power metric + Mp_sol = ode_multistart(@ode45,@(s,Mp) dMetric(s,Mp,A,h,J,drag),int_limit,int_limit(1),zeros(length(cparams)^2,1)); + + Mp = reshape(Mp_sol(int_limit(end)),length(cparams),[]); + +end + + +function localJ = local_body_velocity_J(A,h,J) +% Calculate the Jacobian from shape parameter velocity to local tangential +% and normal velocity + + R = [cos(h(3)) sin(h(3)); + -sin(h(3)) cos(h(3))]; + + %Velocity is sum of body velocity and velocity within body frame + localJ = R*(-A(1:2,:) + J(1:2,:) + [-h(2)*(-A(3,:)); h(1)*(-A(3,:))]); + +end + +function dMp = dMetric(s,Mp,A,h,J,drag) %#ok +% Calculate the local contribution to the power metric + + localJ = local_body_velocity_J(A,h(s),J(s)); + + dMp = localJ.'*drag*localJ; + + dMp = dMp(:); + +end \ No newline at end of file diff --git a/ProgramFiles/v4.2/Utilities/LowRe_toolbox/LowRE_local_connection_from_general_curvature.m b/ProgramFiles/v4.2/Utilities/LowRe_toolbox/LowRE_local_connection_from_general_curvature.m new file mode 100644 index 0000000..c8486cb --- /dev/null +++ b/ProgramFiles/v4.2/Utilities/LowRe_toolbox/LowRE_local_connection_from_general_curvature.m @@ -0,0 +1,122 @@ +function [A, h, J,Omega] = LowRE_local_connection_from_general_curvature(curvdef,cparams,L,c,drag_ratio) +% Calculate the local connection for a set of curvature bases +% +% Inputs: +% kappa_basis_input: cell vector of handles to basis functions +% r: coefficients for basis functions +% L: total length of swimmer +% c: drag per unit length +% drag_ratio: ratio of lateral to longitudinal drag + +% Specified integration limits +int_limit = L*[-0.5 0.5]; + +% First, get the backbone locus and Jacobian functions +[h, J] = backbone_from_general_curvature(curvdef,cparams,L); + +% Now integrate to get the pfaffian +%Omega_sol = ode45( @(s,Omega) connection_helper(s,h(s),J(s)),int_limit,zeros(3,3+length(r))); +Omega_sol = ode_multistart(@ode45,@(s,Omega) connection_helper(s,h(s),J(s),c,drag_ratio),int_limit,int_limit(1),zeros(3,3+length(cparams))); + +Omega = reshape(Omega_sol(int_limit(end)),3,[]); + +A = Omega(:,1:3)\Omega(:,4:end); + +end + +function dOmega = connection_helper(s,h,J,c,drag_ratio) %#ok +% Calculate the derivative of the local connection as it's built up along +% the backbone + + % Midpoint-tangent-frame force jacobian + gdot_to_xi_local = [cos(h(3)) sin(h(3)) 0; + -sin(h(3)) cos(h(3)) 0; + 0 0 1]; + + % Local drag, based on unit longitudinal drag, lateral according to the ratio, no local + % torsional drag, multiplied by drag coefficient + xi_local_to_F_local = [-1 0 0;0 -drag_ratio 0;0 0 0]*c; + + % Force local element applies at midpoint-tangent-frame + F_local_to_F_midpoint = [cos(h(3)) -sin(h(3)) 0; + sin(h(3)) cos(h(3)) 0; + 0 0 1]; + + % Moment around midpoint frame induced by force + F_midpoint_to_FM_midpoint = [1 0 0; 0 1 0; -h(2) h(1) 1]; + + % shape component of pfaffian + omega2 = F_midpoint_to_FM_midpoint * F_local_to_F_midpoint... + * xi_local_to_F_local * gdot_to_xi_local * J; + + % body velocity component of pfaffian + omega1 = F_midpoint_to_FM_midpoint * F_local_to_F_midpoint... + * xi_local_to_F_local * gdot_to_xi_local * [1 0 -h(2); 0 1 h(1); 0 0 1]; + + dOmega = [omega1 omega2]; + dOmega = dOmega(:); + +% % local contribution to the local connection +% dA = omega1\omega2; +% +% % reshape into a vector +% dA = dA(:); + + + +end + +function X = torow( X ) +% +% TOROW Converts a vector or a matrix into a row vector. +% If input is already a row vector, it is returned with no change. +% If input is a column vector, it is converted into a row vector and +% returned. +% If input is a matrix, each column is converted into a row, and all +% resulting rows are placed in series into a single row which is +% returned. +% +% Input: +% X - input vector or matrix +% +% Output: +% X - row vector +% +% Examples: +% torow([ 0 1 2 3 ]) +% returns [ 0 1 2 3 ] +% torow([ 0 1 2 3 ]') +% returns [ 0 1 2 3 ] +% torow([ 0 1; 2 3 ]) +% returns [ 0 2 1 3 ] +% torow([ 0 1; 2 3 ]') +% returns [ 0 1 2 3 ] +% +% Author: Tashi Ravach +% Version: 1.0 +% Date: 07/07/2010 +% + + % check if input is a vector + [ m, n ] = size(X); + if m == 1 + return % input is already a row vector with n columns + elseif n==1 + X = X'; % input is converted from column vector to row vector + elseif (m*n>n) || (m*n>m) + X = X(:)'; % input is converted from matrix to row vector by column + else + X = []; % input is unknown and an empty output is returned + end + +end + +function Y = toz( X ) + + %first, make X a row + X = torow(X); + + % then, make it a 3rd-dimension column + Y(1,1,:) = X; + +end \ No newline at end of file diff --git a/ProgramFiles/v4.2/Utilities/curvature_mode_toolbox/backbone_from_curvature_bases.m b/ProgramFiles/v4.2/Utilities/curvature_mode_toolbox/backbone_from_curvature_bases.m index 43ae91f..5884d3b 100644 --- a/ProgramFiles/v4.2/Utilities/curvature_mode_toolbox/backbone_from_curvature_bases.m +++ b/ProgramFiles/v4.2/Utilities/curvature_mode_toolbox/backbone_from_curvature_bases.m @@ -110,8 +110,6 @@ % delta-functions and the like %summed integral of the curvature functions -%theta_fun = @(s) sum(cell2mat(cellfun( @(k,r) [0 -%r]*(k(torow(s))-k(torow(zeros(size(s))))), kappa_basis,num2cell(r),'UniformOutput',false)),1); theta_fun = @(s) sum(cell2mat(cellfun( @(k,r) theta_fun_helper(k,r,s), kappa_basis,num2cell(r),'UniformOutput',false)),1); diff --git a/ProgramFiles/v4.2/Utilities/curvature_mode_toolbox/backbone_from_general_curvature.m b/ProgramFiles/v4.2/Utilities/curvature_mode_toolbox/backbone_from_general_curvature.m new file mode 100644 index 0000000..6be5d4a --- /dev/null +++ b/ProgramFiles/v4.2/Utilities/curvature_mode_toolbox/backbone_from_general_curvature.m @@ -0,0 +1,87 @@ +function [h, J] = backbone_from_general_curvature(curvdef,cparams,L) + +% Set the integration limits. Internally, everything uses [-0.5 0.5], and +% then scales by the length factor at the end. +all_limits = [-0.5 0 0.5]; + +% Get backbone theta as a function of s +theta_fun = curvdef(cparams,'angle'); + +% Get the x and y locations of points on the backbone +locus_sol = ode_multistart(@ode45, @(s,h) locus_helper(s,theta_fun),all_limits,0,[0;0]); + +% Return the locus data as 3xn matrix with x,y,theta at each of n points +h_norm = @(s) [locus_sol(torow(s)); theta_fun(torow(s))]; % in normalized internal coordinates +h = @(s) [L*locus_sol(torow(s)/L); theta_fun(torow(s)/L)]; % with actual length + +%%%%%%%%%%%%%% +% Get the jacobian to body point velocities + +dcurvdef = curvdef(cparams,'dcurvature_int'); + +%Jacobian of theta function is the derivative of the curvature with respect +%to each of the parameters, as definded in the curvdef function +J_theta_fun = dcurvdef; %ode_multistart(@ode45,@(s,J) dcurvdef(s),all_limits,0,zeros(length(dcurvdef(0)),1)); + +% SE(2) adjoint integration to get the x and y jacobians +jacobian_sol = ode_multistart(@ode45,@(s,J) J_helper(s,J,dcurvdef,h_norm(s)),all_limits,0,zeros(2*length(dcurvdef(0)),1)); + +% Concatenate xy and theta jacobians. +J = @(s) cat(1,reshape(L*jacobian_sol(toz(s/L)),2,[],length(s)),ipermute(J_theta_fun(s/L),[2,3,1])); + + + +end + + +% Get the derivative of the locus along the backbone length +function dh = locus_helper(s,theta_fun) + + % Build the rotation matrix + theta = theta_fun(s); + R = [cos(theta) -sin(theta);sin(theta) cos(theta)]; + + % Rotate the basic tangent vector to point along the backbone + dh = R * [1; 0]; + +end + +% get the derivative of the jacobian along the backbone length (note that J +% is only here for ODE45, and isn't actually used) +function dJ = J_helper(s,J,dcurvdef,h) %#ok + + nvars = length(dcurvdef(0)); + + dJ = zeros(2,nvars); + + dcurvdef_eval = dcurvdef(s); + + for i = 1:nvars + + int_basis = dcurvdef_eval(i,:); + + dJ(:,i) = [-sin(h(3)) -cos(h(3)); + cos(h(3)) -sin(h(3))] * [int_basis; 0]; + + end + + %Reshape into a vector for ODE + dJ = dJ(:); + +end + +function X = torow( X ) % Make the argument a row + + X = X(:).'; + +end + +function Y = toz( X ) + + %first, make X a row + X = torow(X); + + % then, make it a 3rd-dimension column + Y(1,1,:) = X; + +end \ No newline at end of file diff --git a/ProgramFiles/v4.2/Utilities/curvature_mode_toolbox/fatbackbone_from_general_curvature.m b/ProgramFiles/v4.2/Utilities/curvature_mode_toolbox/fatbackbone_from_general_curvature.m new file mode 100644 index 0000000..6cd3227 --- /dev/null +++ b/ProgramFiles/v4.2/Utilities/curvature_mode_toolbox/fatbackbone_from_general_curvature.m @@ -0,0 +1,9 @@ +function B = fatbackbone_from_general_curvature(curvdef,cparams,L,width) + + [h, J] = backbone_from_general_curvature(curvdef,cparams,L); + + + B = fatbackbone(h,L*[-.5 .5],width); + + +end \ No newline at end of file diff --git a/ProgramFiles/v4.2/Utilities/curvature_mode_toolbox/illustrate_backbone_shapespace.m b/ProgramFiles/v4.2/Utilities/curvature_mode_toolbox/illustrate_backbone_shapespace.m new file mode 100644 index 0000000..d16813c --- /dev/null +++ b/ProgramFiles/v4.2/Utilities/curvature_mode_toolbox/illustrate_backbone_shapespace.m @@ -0,0 +1,38 @@ +function illustrate_backbone_shapespace(curvdef,paramvalues) + +% Get the gradient of parameter values in the first two dimensions +[~,grid_spacing_x] = gradient(paramvalues{1}); +[grid_spacing_y,~] = gradient(paramvalues{2}); + +fh = figure(171); +close(fh); +fh = figure(171); + + +axh = axes('Parent',fh); +axis equal +hold on + +for idx = 1:numel(paramvalues{1}) + + % Extract the parameter values for this combination of parameters + p = zeros(size(paramvalues)); + for idx2 = 1:numel(paramvalues) + p(idx2) = paramvalues{idx2}(idx); + end + + % The backbone should be as long as the spacing between the nearby + % elements + blnth = .9 * max(grid_spacing_x(idx),grid_spacing_y(idx)); + + %Generate the backbone locus + B = fatbackbone_from_general_curvature(curvdef,p,blnth,blnth/20); + + for idx2 = 1:2 + B(:,idx2) = B(:,idx2) + p(idx2);% - mean(B(:,idx2)); + end + + % draw the backbone at the specified location + plot(B(:,1),B(:,2),'Parent',axh) + +end \ No newline at end of file diff --git a/ProgramFiles/v4.2/Utilities/curvature_mode_toolbox/make_curvdef.m b/ProgramFiles/v4.2/Utilities/curvature_mode_toolbox/make_curvdef.m new file mode 100644 index 0000000..4478d6b --- /dev/null +++ b/ProgramFiles/v4.2/Utilities/curvature_mode_toolbox/make_curvdef.m @@ -0,0 +1,300 @@ +function make_curvdef(curv_fun_string,paramlist,name,attempt_analytic) +% curv_fun_string can be either a string expression or a symbolic +% expression + +% Load the paths where files should be found +load sysplotter_config +flist = {}; + +% Symbolically calculate all of the requested functions unless told +% otherwise +if ~exist('attempt_analytic','var') + attempt_analytic = [1,1,1,1]; +end + +% Turn the curvature function from a string into a symbolic expression +if attempt_analytic(1) + curv_fun = sym(curv_fun_string); + flist = [flist,{curv_fun}]; +else + flist = [flist,{[]}]; +end + +% Integrate the curvature function along the backbone +if attempt_analytic(2) + int_curv_ds_fun = int(curv_fun,'s',0,'s'); + flist = [flist,{int_curv_ds_fun}]; +else + flist = [flist,{[]}]; +end + +% Take the derivative of the curvature with respect to each parameter +if attempt_analytic(3) + d_curv_dp_fun = jacobian(curv_fun,sym(paramlist)).'; + flist = [flist,{d_curv_dp_fun}]; +else + flist = [flist,{[]}]; +end + +% Integrate the parameter-derivative of the curvature function along the +% backbone +if attempt_analytic(4) + int_d_curv_dp_ds_fun = int(d_curv_dp_fun,'s',0,'s'); + flist = [flist,{int_d_curv_dp_ds_fun}]; +else + flist = [flist,{[]}]; +end + + + +% %%%%%% +% %%%%%% +%%% +%%% +% Print functions out to mfiles in a tempdir + +% List of all the function names +fnamelist = {'curv_fun','int_curv_ds_fun','d_curv_dp_fun','int_d_curv_dp_ds_fun'}; + +% Make a temporary directory to hold mfiles for curvature functions +curvdef_tempdir = fullfile(sysplotterpath,'Utilities','curvature_mode_toolbox',... +'curvdef_tempdir'); + +mkdir(curvdef_tempdir) + +% Test these functions to see if we can turn them into matlabFunctions +valid_mfile = ones(size(fnamelist)); +for idx = 1:numel(fnamelist) + + if attempt_analytic(idx) + + %%%%% + % See if we can turn this into a matlab function + ftest = matlabFunction(flist{idx},'vars',['s';paramlist(:)]); + testvars = num2cell(zeros(size(paramlist))); + + no_analytical = 0; + try + ftest(0,testvars{:}); % Attempt to call the function generated by matlabfunction + catch ME + % If matlab couldn't make an analytical form of the equation, mark + % this + if strcmp(ME.identifier,'MATLAB:UndefinedFunction') + no_analytical = 1; + end + end + + else + + % Mark as no analytical form if none was attempted + no_analytical = 1; + + end + + + %%%% + % If the matlabfunction passes the analytical-form test, save it to a + % file and then concatenate it into the curvdef file. If it fails the + % test, evaluate it at high density, save that to a file, and append a + % file that loads this data and interpolates it. + + if ~no_analytical + % Save the function to an mfile in the tempdir + + matlabFunction(flist{idx},'file',fullfile(curvdef_tempdir,fnamelist{idx}),'vars',['s';paramlist(:)]); + + switch fnamelist{idx} + + case 'curv_fun' + + curv_fun_call = {'output = @(s) curv_fun(s,params{:});'}; + + case 'int_curv_ds_fun' + + int_curv_ds_call = {'output = @(s) int_curv_ds_fun(s,params{:});'}; + + case 'd_curv_dp_fun' + + d_curv_dp_call = {'output = @(s) d_curv_dp_fun(s,params{:});'}; + + + case 'int_d_curv_dp_ds_fun' + + int_d_curv_dp_ds_call = {'output = @(s) int_d_curv_dp_ds_fun(s,params{:});'}; + + otherwise + + error('fnamelist entry is not handled by switch statement') + + end + + else + + % Mark that no mfile was created (so that later code does not + % attempt to copy it in + valid_mfile(idx) = 0; + + % Fill in the function calls + switch fnamelist{idx} + + + case 'curv_fun' + + error('Curvature could not be turned into a matlab function. Either implement alternate behavior in make_curvdef, or correct curvature definition') + + case 'int_curv_ds_fun' + + int_curv_ds_call = {'%% Padded length of unit backbone'... + ,'all_limits = [-.51 0 .51];'... + ,''... + ,'%% Make dummy integration function'... + ,['curv_fun_dummy = curv_' name '(cell2mat(params),''curvature'');']... + ,'curvature = @(s,~) curv_fun_dummy(s);'... + ,''... + ,'%% Integral of the integrand function along s'... + ,'output = ode_multistart(@ode45,curvature,all_limits,0,0);'}; + + case 'd_curv_dp_fun' + + d_curv_dp_call = {'%% Create a dummy function that takes in a vector of parameters'... + ,'%% including s, and deals them into individual function parameters'... + ,'curv_intermediate = @(all_params) vector_to_list_input(@curv_fun,all_params);'... + ,''... + ,'%% Create a function that takes the jacobian of the dummy function'... + ,'fulljacobian = @(s) jacobianest(curv_intermediate,[s,params{:}]);'... + ,''... + ,'%% Create a function that truncates the s-derivative from the full jacobian'... + ,'output = @(s) reshape_truncate_jacobian(fulljacobian(s));'}; + + + case 'int_d_curv_dp_ds_fun' + + int_d_curv_dp_ds_call = {'%% Padded length of unit backbone'... + ,'all_limits = [-.51 0 .51];'... + ,''... + ,'%% Make dummy integration function'... + ,['d_curv_dp_fun_dummy = curv_' name '(cell2mat(params),''dcurvature'');']... + ,'dcurvature = @(s,~) d_curv_dp_fun_dummy(s);'... + ,''... + ,'%% Integral of the integrand function along s'... + ,'output = ode_multistart(@ode45,dcurvature,all_limits,0,zeros(size(params(:).'')));'}; + + otherwise + + error('fnamelist entry is not handled by switch statement') + + end + + end + + + +end + + + +% Open the template file +fidi = fopen(fullfile(sysplotterpath,'Utilities','curvature_mode_toolbox',... + 'make_curvdef_template.txt')); + +% Create the output file +fido = fopen(fullfile(syspath,['curv_' name '.m']),'w'); + +% Place the input function string as a comment in the top line of the file +fprintf(fido,'%% %s\n',char(curv_fun_string)); + +while ~feof(fidi) + + % Read the next line from the input file + s = fgetl(fidi); + + + % Replace the function calls + switch strtrim(s) + + case 'AA_curv_fun_call' + + s = sprintf('\t\t%s\n',curv_fun_call{:}); + + case 'AA_int_curv_ds_call' + + s = sprintf('\t\t%s\n',int_curv_ds_call{:}); + + case 'AA_d_curv_dp_call' + + s = sprintf('\t\t%s\n',d_curv_dp_call{:}); + + case 'AA_int_d_curv_dp_ds_call' + + s = sprintf('\t\t%s\n',int_d_curv_dp_ds_call{:}); + + otherwise + + % Insert the filename + s = strrep(s,'AA_CURVFILENAME',['curv_' name]); + + end + + + % Copy put the line into the new file + fprintf(fido,'%s\n',s); + +end + + +% Close the template function +fclose(fidi); + +% Concatenate any mfile functions onto the end of the template + +for idx = 1:numel(flist) + + if valid_mfile(idx) %Only copy valid mfiles into the + + % Insert a blank line + fprintf(fido,'%s\n',[]); + firstline = 1; + + % Load the matlab function + fidi = fopen(fullfile(curvdef_tempdir,[fnamelist{idx} '.m'])); + + + + % Copy all lines of the function into the curvdef file + while ~feof(fidi) + + % Read the next line from the input file + s = fgetl(fidi); + + % Print the line unless it is a comment or empty + if ~strncmp(s,'%%',1) && ~strcmp(s,'') + + % If it is not the first line, indent + if firstline + firstline = 0; + else + fprintf(fido,'%s\t',[]); + end + + fprintf(fido,'%s\n',s); + end + + end + + % Insert an end command and an extra blank line + fprintf(fido,'%s\n\n','end'); + + fclose(fidi); + + end + +end + + + +% cleanup +rmdir(curvdef_tempdir,'s') + +fclose(fido); + +end \ No newline at end of file diff --git a/ProgramFiles/v4.2/Utilities/curvature_mode_toolbox/make_curvdef_template.txt b/ProgramFiles/v4.2/Utilities/curvature_mode_toolbox/make_curvdef_template.txt new file mode 100644 index 0000000..b5cb3a1 --- /dev/null +++ b/ProgramFiles/v4.2/Utilities/curvature_mode_toolbox/make_curvdef_template.txt @@ -0,0 +1,41 @@ +function output = AA_CURVFILENAME(params,mode) + +% Turn params into a cell matrix +params = num2cell(params); + +switch mode + + case 'curvature' + + AA_curv_fun_call + + case 'angle' + + AA_int_curv_ds_call + + case 'dcurvature' + + AA_d_curv_dp_call + + case 'dcurvature_int' + + AA_int_d_curv_dp_ds_call + +end + +end + +function output = vector_to_list_input(funhandle,all_params) + + all_params = num2cell(all_params); + + output = funhandle(all_params{:}); + +end + + +function output = reshape_truncate_jacobian(J) + + output = J(2:end)'; + +end \ No newline at end of file diff --git a/ProgramFiles/v4.2/Utilities/curvature_mode_toolbox/pinched_sin_scale.m b/ProgramFiles/v4.2/Utilities/curvature_mode_toolbox/pinched_sin_scale.m new file mode 100644 index 0000000..e5d93c2 --- /dev/null +++ b/ProgramFiles/v4.2/Utilities/curvature_mode_toolbox/pinched_sin_scale.m @@ -0,0 +1,65 @@ +function [pinch_wave_matched]= pinched_sin_scale + +% Generate a scale factor so that the maximum angle reached by the sin and +% pinched sin basis functions is the same +n = logspace(-1,1.5,100); + +x = linspace(0,pi/2); + +base_area = trapz(x,cos(x)); + +for idx = 1:numel(n) + + scale_factor(idx) = base_area/trapz(x,cos(x).^n(idx)); + +end + +P = polyfit(n,scale_factor,3); + +% figure(18) +% +% plot(n,scale_factor,n,polyval(P,n)) + + +% Now generate the pinched sin function + +syms s n A phi a1 a2 + +% Raw wave function +raw_wave = cos(2*pi*s+phi); + +% Pinched in wave +pinched_wave = A * (abs(raw_wave))^n * sign(raw_wave); + +% Substitute in the a1 and a2 values +amp = (a1^2+a2^2)^.5; +phase = atan2(a2,a1); +pinched_wave_restored = subs(pinched_wave,[A phi],[amp phase]); + + +% % Wave with sin and cosine amplitudes +% raw_wave = (a1*cos(2*pi*s)+a2*sin(2*pi*s)); +% +% % Amplitude of the wave +% amp_wave = (a1^2+a2^2)^(.5); +% +% % Normalize the wave +% norm_wave = raw_wave/amp_wave; +% +% % Pinch the wave +% pinched_wave = (abs(norm_wave)^n) * sign(norm_wave); +% +% % Restore the scale of the wave +% pinched_wave_restored = pinched_wave*amp_wave; + + + + +% Scale to match max angle with original wave +pinch_wave_matched = pinched_wave_restored * (P*[n^3;n^2;n;1]); + + + + + +end diff --git a/ProgramFiles/v4.2/Utilities/evaluate_displacement_and_cost.m b/ProgramFiles/v4.2/Utilities/evaluate_displacement_and_cost.m index 8b05d6f..7aaf7f3 100644 --- a/ProgramFiles/v4.2/Utilities/evaluate_displacement_and_cost.m +++ b/ProgramFiles/v4.2/Utilities/evaluate_displacement_and_cost.m @@ -6,7 +6,7 @@ % S should be a sysplotter 's' structure loaded from a file % sysf_FILENAME_calc.mat (note the _calc suffix) % -% P should be a structure with fields "phi" and "dphi", returning a +% P should be a structure with fields "phi_def" and "dphi_def", returning a % vector of shapes and shape velocities respectively. If it is not % convenient to analytically describe the shape velocity function, % gait.dphi should be defined as @@ -70,7 +70,7 @@ % Convert the final motion into its representation in optimal % coordinates - startshape = p.phi(0); + startshape = p.phi_def{1}{1}(0); startshapelist = num2cell(startshape); beta_theta = interpn(s.grid.eval{:},s.B_optimized.eval.Beta{3},startshapelist{:},'cubic'); net_disp_opt = [cos(beta_theta) sin(beta_theta) 0;... @@ -84,7 +84,7 @@ % Convert the final motion into its representation in optimal % coordinates - startshape = p.phi(0); + startshape = p.phi_def{1}{1}(0); startshapelist = num2cell(startshape); beta_theta = interpn(s.grid.eval{:},s.B_optimized.eval.Beta{3},startshapelist{:},'cubic'); net_disp_opt = [cos(beta_theta) sin(beta_theta) 0;... @@ -99,9 +99,9 @@ function [xi, dcost] = get_velocities(t,s,gait,ConnectionEval) % Get the shape and shape derivative at the current time - shape = gait.phi(t); + shape = gait.phi_def{1}{1}(t); shapelist = num2cell(shape); - dshape = gait.dphi(t); + dshape = gait.dphi_def{1}{1}(t); % Get the local connection and metric at the current time, in the new coordinates @@ -116,7 +116,7 @@ A = -cellfun(@(C) interpn(s.grid.eval{:},C,shapelist{:}),s.vecfield.eval.content.Avec); - M = cellfun(@(C) interpn(s.grid.eval{:},C,shapelist{:}),s.metricfield.eval.content.metric); + M = cellfun(@(C) interpn(s.grid.metric_eval{:},C,shapelist{:}),s.metricfield.metric_eval.content.metric); otherwise error('Unknown method for evaluating local connection'); diff --git a/ProgramFiles/v4.2/Utilities/gait_gui_draw/gait_gui_draw.m b/ProgramFiles/v4.2/Utilities/gait_gui_draw/gait_gui_draw.m index e683d65..4e4c015 100644 --- a/ProgramFiles/v4.2/Utilities/gait_gui_draw/gait_gui_draw.m +++ b/ProgramFiles/v4.2/Utilities/gait_gui_draw/gait_gui_draw.m @@ -1,13 +1,27 @@ -function gait_gui_draw +function gait_gui_draw(hAx) % Load the sysplotter configuration information load sysplotter_config -% Bring current axis to foreground -figure(gcf) +% Bring selected axis to foreground (use gca if this is not specified) +if ~exist('hAx','var') || isempty(hAx) + hAx = gca; +end + +found_figure = 0; +hSearch = hAx; +while ~found_figure + hSearch = get(hSearch,'Parent'); + if strcmpi(get(hSearch,'type'),'figure') + found_figure = 1; + hFig = hSearch; + end +end +figure(hFig) + % Use the mouse to select a series of points -[alpha1,alpha2,button] = ginputc('ShowPoints',true,'Color',Colorset.spot); +[alpha1,alpha2,button] = ginputc('AxHandle',hAx,'ShowPoints',true,'Color',Colorset.spot); % Unless a non-primary button was used for the last click, copy the first % point to the end to make a closed loop @@ -21,10 +35,14 @@ % Evenly space the points along one period t = linspace(0,period,numel(alpha1)); - % Generate spline structures for the selected points - spline_alpha1 = csape(t,alpha1','periodic'); - spline_alpha2 = csape(t,alpha2','periodic'); - + % Fit a periodic spline to the selected points; the endslopes are found + % by averaging the positions of the points before and after the join + endslope1 = (alpha1(2)-alpha1(end-1))/(t(end)-t(end-2)); + endslope2 = (alpha2(2)-alpha2(end-1))/(t(end)-t(end-2)); + spline_alpha1 = spline(t,[endslope1;alpha1(:);endslope1]); + spline_alpha2 = spline(t,[endslope2;alpha2(:);endslope2]); + + else % Set the period of the open motion as 1 @@ -33,20 +51,31 @@ % Evenly space the points along time from zero to 2pi t = linspace(0,period,numel(alpha1)); - % Generate spline structures for the selected points - spline_alpha1 = csape(t,alpha1'); - spline_alpha2 = csape(t,alpha2'); - + % Fit a non-periodic spline to the selected points + spline_alpha1 = spline(t,alpha1(:)); + spline_alpha2 = spline(t,alpha2(:)); end % Upsample and plot to show gait to user n_plot = 100; t_plot = linspace(0,period,n_plot); -alpha1_plot = fnval(spline_alpha1,t_plot); -alpha2_plot = fnval(spline_alpha2,t_plot); -gaitline = line(alpha1_plot,alpha2_plot,'Color',Colorset.spot,'LineWidth',5); +alpha1_plot = ppval(spline_alpha1,t_plot); +alpha2_plot = ppval(spline_alpha2,t_plot); + +% Provide zdata to line if necessary +maxZ = 0; +hAxChildren = get(hAx,'Children'); +if ~isempty(hAxChildren) + for idx = 1:numel(hAxChildren) + if ~isempty(hAxChildren(idx).ZData) + maxZ = max(maxZ,max(hAxChildren(idx).ZData(:))); + end + end +end + +gaitline = line('Parent',hAx,'XData',alpha1_plot,'YData',alpha2_plot,'ZData',maxZ*ones(size(alpha1_plot)),'Color',Colorset.spot,'LineWidth',5); %%%% Ask the user for a filename current_dir = pwd; % Remember where we started @@ -66,8 +95,8 @@ [~,paramfilenamebare,ext] = fileparts(paramfilename); % Break open the filename to get the extension end - if ( strncmp(paramfilename,'params_',7) && strcmp(ext,'.mat') )... - && strcmp(fullfile(pathname,paramfilename),fullfile(shchpath,paramfilename)) + if ( strncmpi(paramfilename,'params_',7) && strcmpi(ext,'.mat') )... + && strcmpi(fullfile(pathname,paramfilename),fullfile(shchpath,paramfilename)) usercancel = 0; break else @@ -97,15 +126,8 @@ end - % Refresh the sysplotter gui - try - set(groot,'ShowHiddenHandles','on'); % Crack open the sysplotter visibility - refresh_handle = findobj('tag','refresh_gui'); % Get the handle for the button - refresh_handle.Callback(refresh_handle,0) % Push the refresh button - set(groot,'ShowHiddenHandles','off'); % Reseal the sysplotter - catch - set(groot,'ShowHiddenHandles','off'); % Make sure the sysplotter goes back to hidden handles - end + refresh_handle = findall(0,'tag','refresh_gui'); % Get the handle for the button + refresh_handle.Callback(refresh_handle,0) % Push the refresh button end diff --git a/ProgramFiles/v4.2/Utilities/gait_gui_draw/gait_gui_draw_template.txt b/ProgramFiles/v4.2/Utilities/gait_gui_draw/gait_gui_draw_template.txt index 954d00e..04078a0 100644 --- a/ProgramFiles/v4.2/Utilities/gait_gui_draw/gait_gui_draw_template.txt +++ b/ProgramFiles/v4.2/Utilities/gait_gui_draw/gait_gui_draw_template.txt @@ -39,23 +39,39 @@ function output = AA_SHCHFILENAME(input_mode) end % Generate spline structures for the selected points - spline_alpha1 = csape(t,alpha1',splinemode); - spline_alpha2 = csape(t,alpha2',splinemode); + switch splinemode + + case 'periodic' + + % Fit a periodic spline to the selected points; the endslopes are found + % by averaging the positions of the points before and after the join + endslope1 = (alpha1(2)-alpha1(end-1))/(t(end)-t(end-2)); + endslope2 = (alpha2(2)-alpha2(end-1))/(t(end)-t(end-2)); + spline_alpha1 = spline(t,[endslope1;alpha1(:);endslope1]); + spline_alpha2 = spline(t,[endslope2;alpha2(:);endslope2]); + + case 'complete' + + % Fit a non-periodic spline to the selected points + spline_alpha1 = spline(t,alpha1(:)); + spline_alpha2 = spline(t,alpha2(:)); + + end % The gait path is now defined by evaluating the two splines at % specified times - p.phi_def = @(t) [fnval(spline_alpha1,t(:)),fnval(spline_alpha2,t(:))]; + p.phi_def = @(t) [ppval(spline_alpha1,t(:)),ppval(spline_alpha2,t(:))]; % Speed up execution by defining the gait velocity as well (so % the evaluator doesn't need to numerically take the derivative % at every time - spline_dalpha1 = fnder(spline_alpha1); - spline_dalpha2 = fnder(spline_alpha2); + spline_dalpha1 = ppdiff(spline_alpha1); + spline_dalpha2 = ppdiff(spline_alpha2); % The gait path is now defined by evaluating the two splin % derivatives at specified times - p.dphi_def = @(t) [fnval(spline_dalpha1,t(:)),fnval(spline_dalpha2,t(:))]; + p.dphi_def = @(t) [ppval(spline_dalpha1,t(:)),ppval(spline_dalpha2,t(:))]; %marker locations diff --git a/ProgramFiles/v4.2/Utilities/ginputc.m b/ProgramFiles/v4.2/Utilities/ginputc.m index 407a034..89825a1 100644 --- a/ProgramFiles/v4.2/Utilities/ginputc.m +++ b/ProgramFiles/v4.2/Utilities/ginputc.m @@ -59,6 +59,8 @@ % Jiro Doke % October 19, 2012 % Copyright 2012 The MathWorks, Inc. +% +% Small modifications by Ross Hatton, 2016 try if verLessThan('matlab', '7.5') @@ -75,21 +77,23 @@ addOptional(p, 'N', inf, @(x) validateattributes(x, {'numeric'}, ... {'scalar', 'integer', 'positive'})); -addParamValue(p, 'FigHandle', [], @(x) numel(x)==1 && ishandle(x)); -addParamValue(p, 'Color', 'k', @colorValidFcn); -addParamValue(p, 'LineWidth', 0.5 , @(x) validateattributes(x, ... +addParameter(p, 'FigHandle', [], @(x) numel(x)==1 && ishandle(x)); +addParameter(p, 'AxHandle', [], @(x) numel(x)==1 && ishandle(x)); +addParameter(p, 'Color', 'k', @colorValidFcn); +addParameter(p, 'LineWidth', 0.5 , @(x) validateattributes(x, ... {'numeric'}, {'scalar', 'positive'})); -addParamValue(p, 'LineStyle', '-' , @(x) validatestring(x, ... +addParameter(p, 'LineStyle', '-' , @(x) validatestring(x, ... {'-', '--', '-.', ':'})); -addParamValue(p, 'ShowPoints', false, @(x) validateattributes(x, ... +addParameter(p, 'ShowPoints', false, @(x) validateattributes(x, ... {'logical'}, {'scalar'})); -addParamValue(p, 'ConnectPoints', true, @(x) validateattributes(x, ... +addParameter(p, 'ConnectPoints', true, @(x) validateattributes(x, ... {'logical'}, {'scalar'})); parse(p, varargin{:}); N = p.Results.N; hFig = p.Results.FigHandle; +hAx = p.Results.AxHandle; color = p.Results.Color; linewidth = p.Results.LineWidth; linestyle = p.Results.LineStyle; @@ -118,13 +122,25 @@ end % Try to get the current axes even if it has a hidden handle. -hAx = get(hFig, 'CurrentAxes'); if isempty(hAx) - allAx = findall(hFig, 'Type', 'axes'); - if ~isempty(allAx) - hAx = allAx(1); - else - hAx = axes('Parent', hFig); + hAx = get(hFig, 'CurrentAxes'); + if isempty(hAx) + allAx = findall(hFig, 'Type', 'axes'); + if ~isempty(allAx) + hAx = allAx(1); + else + hAx = axes('Parent', hFig); + end + end +else + found_figure = 0; + hSearch = hAx; + while ~found_figure + hSearch = get(hSearch,'Parent'); + if strcmpi(get(hSearch,'type'),'figure') + found_figure = 1; + hFig = hSearch; + end end end @@ -166,14 +182,25 @@ % Create an invisible axes for displaying the full crosshair cursor hInvisibleAxes = axes(... - 'Parent', hFig, ... + 'Parent', get(hAx,'Parent'), ... 'Units', 'normalized', ... - 'Position', [0 0 1 1], ... + 'Position', [0 0 1 1],...get(hAx,'Position'), ... 'XLim', [0 1], ... 'YLim', [0 1], ... 'HitTest', 'off', ... 'HandleVisibility', 'off', ... 'Visible', 'off'); +%uistack(hInvisibleAxes,'top') + +% hInvisibleAxes = axes(... +% 'Parent', hFig, ... +% 'Units', 'normalized', ... +% 'Position', [0 0 1 1], ... +% 'XLim', [0 1], ... +% 'YLim', [0 1], ... +% 'HitTest', 'off', ... +% 'HandleVisibility', 'off', ... +% 'Visible', 'off'); % Create line object for the selected points if showpoints @@ -233,6 +260,25 @@ % Wait until enter is pressed. uiwait(hFig); + for idx = 1:length(propsToChange) + set(allHG(validObjects(:, idx)), propsToChange(idx), curCallbacks{idx}); + end + + % Restore window functions and pointer + % set(hFig, 'WindowButtonDownFcn', curWBDF); + % set(hFig, 'WindowButtonMotionFcn', curWBMF); + % set(hFig, 'WindowButtonUpFcn', curWBUF); + % set(hFig, 'KeyPressFcn', curKPF); + % set(hFig, 'KeyReleaseFcn', curKRF); + % set(hFig, 'ResizeFcn', curRF); + + % Restore pointer + set(hFig, 'Pointer', curPointer); + set(hFig, 'PointerShapeCData', curPointerShapeCData); + + % Delete invisible axes and return control + delete(hInvisibleAxes); + %-------------------------------------------------------------------------- function mouseMoveFcn(varargin) @@ -304,7 +350,6 @@ function keyPressFcn(obj, edata) %#ok function updatePoints(clickType) % This function captures the information for the selected point - hAx = gca; pt = get(hAx, 'CurrentPoint'); x = [x; pt(1)]; y = [y; pt(3)]; @@ -384,24 +429,24 @@ function resizeFcn(varargin) function exitFcn() % This function exits GINPUTC and restores previous figure settings - for idx = 1:length(propsToChange) - set(allHG(validObjects(:, idx)), propsToChange(idx), curCallbacks{idx}); - end - - % Restore window functions and pointer - % set(hFig, 'WindowButtonDownFcn', curWBDF); - % set(hFig, 'WindowButtonMotionFcn', curWBMF); - % set(hFig, 'WindowButtonUpFcn', curWBUF); - % set(hFig, 'KeyPressFcn', curKPF); - % set(hFig, 'KeyReleaseFcn', curKRF); - % set(hFig, 'ResizeFcn', curRF); - - % Restore pointer - set(hFig, 'Pointer', curPointer); - set(hFig, 'PointerShapeCData', curPointerShapeCData); - - % Delete invisible axes and return control - delete(hInvisibleAxes); +% for idx = 1:length(propsToChange) +% set(allHG(validObjects(:, idx)), propsToChange(idx), curCallbacks{idx}); +% end +% +% % Restore window functions and pointer +% % set(hFig, 'WindowButtonDownFcn', curWBDF); +% % set(hFig, 'WindowButtonMotionFcn', curWBMF); +% % set(hFig, 'WindowButtonUpFcn', curWBUF); +% % set(hFig, 'KeyPressFcn', curKPF); +% % set(hFig, 'KeyReleaseFcn', curKRF); +% % set(hFig, 'ResizeFcn', curRF); +% +% % Restore pointer +% set(hFig, 'Pointer', curPointer); +% set(hFig, 'PointerShapeCData', curPointerShapeCData); +% +% % Delete invisible axes and return control +% delete(hInvisibleAxes); uiresume(hFig); end %-------------------------------------------------------------------------- diff --git a/ProgramFiles/v4.2/Utilities/inpaint_nans.m b/ProgramFiles/v4.2/Utilities/inpaint_nans.m new file mode 100755 index 0000000..2460b51 --- /dev/null +++ b/ProgramFiles/v4.2/Utilities/inpaint_nans.m @@ -0,0 +1 @@ +function B=inpaint_nans(A,method) % INPAINT_NANS: in-paints over nans in an array % usage: B=INPAINT_NANS(A) % default method % usage: B=INPAINT_NANS(A,method) % specify method used % % Solves approximation to one of several pdes to % interpolate and extrapolate holes in an array % % arguments (input): % A - nxm array with some NaNs to be filled in % % method - (OPTIONAL) scalar numeric flag - specifies % which approach (or physical metaphor to use % for the interpolation.) All methods are capable % of extrapolation, some are better than others. % There are also speed differences, as well as % accuracy differences for smooth surfaces. % % methods {0,1,2} use a simple plate metaphor. % method 3 uses a better plate equation, % but may be much slower and uses % more memory. % method 4 uses a spring metaphor. % method 5 is an 8 neighbor average, with no % rationale behind it compared to the % other methods. I do not recommend % its use. % % method == 0 --> (DEFAULT) see method 1, but % this method does not build as large of a % linear system in the case of only a few % NaNs in a large array. % Extrapolation behavior is linear. % % method == 1 --> simple approach, applies del^2 % over the entire array, then drops those parts % of the array which do not have any contact with % NaNs. Uses a least squares approach, but it % does not modify known values. % In the case of small arrays, this method is % quite fast as it does very little extra work. % Extrapolation behavior is linear. % % method == 2 --> uses del^2, but solving a direct % linear system of equations for nan elements. % This method will be the fastest possible for % large systems since it uses the sparsest % possible system of equations. Not a least % squares approach, so it may be least robust % to noise on the boundaries of any holes. % This method will also be least able to % interpolate accurately for smooth surfaces. % Extrapolation behavior is linear. % % Note: method 2 has problems in 1-d, so this % method is disabled for vector inputs. % % method == 3 --+ See method 0, but uses del^4 for % the interpolating operator. This may result % in more accurate interpolations, at some cost % in speed. % % method == 4 --+ Uses a spring metaphor. Assumes % springs (with a nominal length of zero) % connect each node with every neighbor % (horizontally, vertically and diagonally) % Since each node tries to be like its neighbors, % extrapolation is as a constant function where % this is consistent with the neighboring nodes. % % method == 5 --+ See method 2, but use an average % of the 8 nearest neighbors to any element. % This method is NOT recommended for use. % % % arguments (output): % B - nxm array with NaNs replaced % % % Example: % [x,y] = meshgrid(0:.01:1); % z0 = exp(x+y); % znan = z0; % znan(20:50,40:70) = NaN; % znan(30:90,5:10) = NaN; % znan(70:75,40:90) = NaN; % % z = inpaint_nans(znan); % % % See also: griddata, interp1 % % Author: John D'Errico % e-mail address: woodchips@rochester.rr.com % Release: 2 % Release date: 4/15/06 % I always need to know which elements are NaN, % and what size the array is for any method [n,m]=size(A); A=A(:); nm=n*m; k=isnan(A(:)); % list the nodes which are known, and which will % be interpolated nan_list=find(k); known_list=find(~k); % how many nans overall nan_count=length(nan_list); % convert NaN indices to (r,c) form % nan_list==find(k) are the unrolled (linear) indices % (row,column) form [nr,nc]=ind2sub([n,m],nan_list); % both forms of index in one array: % column 1 == unrolled index % column 2 == row index % column 3 == column index nan_list=[nan_list,nr,nc]; % supply default method if (nargin<2) || isempty(method) method = 0; elseif ~ismember(method,0:5) error 'If supplied, method must be one of: {0,1,2,3,4,5}.' end % for different methods switch method case 0 % The same as method == 1, except only work on those % elements which are NaN, or at least touch a NaN. % is it 1-d or 2-d? if (m == 1) || (n == 1) % really a 1-d case work_list = nan_list(:,1); work_list = unique([work_list;work_list - 1;work_list + 1]); work_list(work_list <= 1) = []; work_list(work_list >= nm) = []; nw = numel(work_list); u = (1:nw)'; fda = sparse(repmat(u,1,3),bsxfun(@plus,work_list,-1:1), ... repmat([1 -2 1],nw,1),nw,nm); else % a 2-d case % horizontal and vertical neighbors only talks_to = [-1 0;0 -1;1 0;0 1]; neighbors_list=identify_neighbors(n,m,nan_list,talks_to); % list of all nodes we have identified all_list=[nan_list;neighbors_list]; % generate sparse array with second partials on row % variable for each element in either list, but only % for those nodes which have a row index > 1 or < n L = find((all_list(:,2) > 1) & (all_list(:,2) < n)); nl=length(L); if nl>0 fda=sparse(repmat(all_list(L,1),1,3), ... repmat(all_list(L,1),1,3)+repmat([-1 0 1],nl,1), ... repmat([1 -2 1],nl,1),nm,nm); else fda=spalloc(n*m,n*m,size(all_list,1)*5); end % 2nd partials on column index L = find((all_list(:,3) > 1) & (all_list(:,3) < m)); nl=length(L); if nl>0 fda=fda+sparse(repmat(all_list(L,1),1,3), ... repmat(all_list(L,1),1,3)+repmat([-n 0 n],nl,1), ... repmat([1 -2 1],nl,1),nm,nm); end end % eliminate knowns rhs=-fda(:,known_list)*A(known_list); k=find(any(fda(:,nan_list(:,1)),2)); % and solve... B=A; B(nan_list(:,1))=fda(k,nan_list(:,1))\rhs(k); case 1 % least squares approach with del^2. Build system % for every array element as an unknown, and then % eliminate those which are knowns. % Build sparse matrix approximating del^2 for % every element in A. % is it 1-d or 2-d? if (m == 1) || (n == 1) % a 1-d case u = (1:(nm-2))'; fda = sparse(repmat(u,1,3),bsxfun(@plus,u,0:2), ... repmat([1 -2 1],nm-2,1),nm-2,nm); else % a 2-d case % Compute finite difference for second partials % on row variable first [i,j]=ndgrid(2:(n-1),1:m); ind=i(:)+(j(:)-1)*n; np=(n-2)*m; fda=sparse(repmat(ind,1,3),[ind-1,ind,ind+1], ... repmat([1 -2 1],np,1),n*m,n*m); % now second partials on column variable [i,j]=ndgrid(1:n,2:(m-1)); ind=i(:)+(j(:)-1)*n; np=n*(m-2); fda=fda+sparse(repmat(ind,1,3),[ind-n,ind,ind+n], ... repmat([1 -2 1],np,1),nm,nm); end % eliminate knowns rhs=-fda(:,known_list)*A(known_list); k=find(any(fda(:,nan_list),2)); % and solve... B=A; B(nan_list(:,1))=fda(k,nan_list(:,1))\rhs(k); case 2 % Direct solve for del^2 BVP across holes % generate sparse array with second partials on row % variable for each nan element, only for those nodes % which have a row index > 1 or < n % is it 1-d or 2-d? if (m == 1) || (n == 1) % really just a 1-d case error('Method 2 has problems for vector input. Please use another method.') else % a 2-d case L = find((nan_list(:,2) > 1) & (nan_list(:,2) < n)); nl=length(L); if nl>0 fda=sparse(repmat(nan_list(L,1),1,3), ... repmat(nan_list(L,1),1,3)+repmat([-1 0 1],nl,1), ... repmat([1 -2 1],nl,1),n*m,n*m); else fda=spalloc(n*m,n*m,size(nan_list,1)*5); end % 2nd partials on column index L = find((nan_list(:,3) > 1) & (nan_list(:,3) < m)); nl=length(L); if nl>0 fda=fda+sparse(repmat(nan_list(L,1),1,3), ... repmat(nan_list(L,1),1,3)+repmat([-n 0 n],nl,1), ... repmat([1 -2 1],nl,1),n*m,n*m); end % fix boundary conditions at extreme corners % of the array in case there were nans there if ismember(1,nan_list(:,1)) fda(1,[1 2 n+1])=[-2 1 1]; end if ismember(n,nan_list(:,1)) fda(n,[n, n-1,n+n])=[-2 1 1]; end if ismember(nm-n+1,nan_list(:,1)) fda(nm-n+1,[nm-n+1,nm-n+2,nm-n])=[-2 1 1]; end if ismember(nm,nan_list(:,1)) fda(nm,[nm,nm-1,nm-n])=[-2 1 1]; end % eliminate knowns rhs=-fda(:,known_list)*A(known_list); % and solve... B=A; k=nan_list(:,1); B(k)=fda(k,k)\rhs(k); end case 3 % The same as method == 0, except uses del^4 as the % interpolating operator. % del^4 template of neighbors talks_to = [-2 0;-1 -1;-1 0;-1 1;0 -2;0 -1; ... 0 1;0 2;1 -1;1 0;1 1;2 0]; neighbors_list=identify_neighbors(n,m,nan_list,talks_to); % list of all nodes we have identified all_list=[nan_list;neighbors_list]; % generate sparse array with del^4, but only % for those nodes which have a row & column index % >= 3 or <= n-2 L = find( (all_list(:,2) >= 3) & ... (all_list(:,2) <= (n-2)) & ... (all_list(:,3) >= 3) & ... (all_list(:,3) <= (m-2))); nl=length(L); if nl>0 % do the entire template at once fda=sparse(repmat(all_list(L,1),1,13), ... repmat(all_list(L,1),1,13) + ... repmat([-2*n,-n-1,-n,-n+1,-2,-1,0,1,2,n-1,n,n+1,2*n],nl,1), ... repmat([1 2 -8 2 1 -8 20 -8 1 2 -8 2 1],nl,1),nm,nm); else fda=spalloc(n*m,n*m,size(all_list,1)*5); end % on the boundaries, reduce the order around the edges L = find((((all_list(:,2) == 2) | ... (all_list(:,2) == (n-1))) & ... (all_list(:,3) >= 2) & ... (all_list(:,3) <= (m-1))) | ... (((all_list(:,3) == 2) | ... (all_list(:,3) == (m-1))) & ... (all_list(:,2) >= 2) & ... (all_list(:,2) <= (n-1)))); nl=length(L); if nl>0 fda=fda+sparse(repmat(all_list(L,1),1,5), ... repmat(all_list(L,1),1,5) + ... repmat([-n,-1,0,+1,n],nl,1), ... repmat([1 1 -4 1 1],nl,1),nm,nm); end L = find( ((all_list(:,2) == 1) | ... (all_list(:,2) == n)) & ... (all_list(:,3) >= 2) & ... (all_list(:,3) <= (m-1))); nl=length(L); if nl>0 fda=fda+sparse(repmat(all_list(L,1),1,3), ... repmat(all_list(L,1),1,3) + ... repmat([-n,0,n],nl,1), ... repmat([1 -2 1],nl,1),nm,nm); end L = find( ((all_list(:,3) == 1) | ... (all_list(:,3) == m)) & ... (all_list(:,2) >= 2) & ... (all_list(:,2) <= (n-1))); nl=length(L); if nl>0 fda=fda+sparse(repmat(all_list(L,1),1,3), ... repmat(all_list(L,1),1,3) + ... repmat([-1,0,1],nl,1), ... repmat([1 -2 1],nl,1),nm,nm); end % eliminate knowns rhs=-fda(:,known_list)*A(known_list); k=find(any(fda(:,nan_list(:,1)),2)); % and solve... B=A; B(nan_list(:,1))=fda(k,nan_list(:,1))\rhs(k); case 4 % Spring analogy % interpolating operator. % list of all springs between a node and a horizontal % or vertical neighbor hv_list=[-1 -1 0;1 1 0;-n 0 -1;n 0 1]; hv_springs=[]; for i=1:4 hvs=nan_list+repmat(hv_list(i,:),nan_count,1); k=(hvs(:,2)>=1) & (hvs(:,2)<=n) & (hvs(:,3)>=1) & (hvs(:,3)<=m); hv_springs=[hv_springs;[nan_list(k,1),hvs(k,1)]]; end % delete replicate springs hv_springs=unique(sort(hv_springs,2),'rows'); % build sparse matrix of connections, springs % connecting diagonal neighbors are weaker than % the horizontal and vertical springs nhv=size(hv_springs,1); springs=sparse(repmat((1:nhv)',1,2),hv_springs, ... repmat([1 -1],nhv,1),nhv,nm); % eliminate knowns rhs=-springs(:,known_list)*A(known_list); % and solve... B=A; B(nan_list(:,1))=springs(:,nan_list(:,1))\rhs; case 5 % Average of 8 nearest neighbors % generate sparse array to average 8 nearest neighbors % for each nan element, be careful around edges fda=spalloc(n*m,n*m,size(nan_list,1)*9); % -1,-1 L = find((nan_list(:,2) > 1) & (nan_list(:,3) > 1)); nl=length(L); if nl>0 fda=fda+sparse(repmat(nan_list(L,1),1,2), ... repmat(nan_list(L,1),1,2)+repmat([-n-1, 0],nl,1), ... repmat([1 -1],nl,1),n*m,n*m); end % 0,-1 L = find(nan_list(:,3) > 1); nl=length(L); if nl>0 fda=fda+sparse(repmat(nan_list(L,1),1,2), ... repmat(nan_list(L,1),1,2)+repmat([-n, 0],nl,1), ... repmat([1 -1],nl,1),n*m,n*m); end % +1,-1 L = find((nan_list(:,2) < n) & (nan_list(:,3) > 1)); nl=length(L); if nl>0 fda=fda+sparse(repmat(nan_list(L,1),1,2), ... repmat(nan_list(L,1),1,2)+repmat([-n+1, 0],nl,1), ... repmat([1 -1],nl,1),n*m,n*m); end % -1,0 L = find(nan_list(:,2) > 1); nl=length(L); if nl>0 fda=fda+sparse(repmat(nan_list(L,1),1,2), ... repmat(nan_list(L,1),1,2)+repmat([-1, 0],nl,1), ... repmat([1 -1],nl,1),n*m,n*m); end % +1,0 L = find(nan_list(:,2) < n); nl=length(L); if nl>0 fda=fda+sparse(repmat(nan_list(L,1),1,2), ... repmat(nan_list(L,1),1,2)+repmat([1, 0],nl,1), ... repmat([1 -1],nl,1),n*m,n*m); end % -1,+1 L = find((nan_list(:,2) > 1) & (nan_list(:,3) < m)); nl=length(L); if nl>0 fda=fda+sparse(repmat(nan_list(L,1),1,2), ... repmat(nan_list(L,1),1,2)+repmat([n-1, 0],nl,1), ... repmat([1 -1],nl,1),n*m,n*m); end % 0,+1 L = find(nan_list(:,3) < m); nl=length(L); if nl>0 fda=fda+sparse(repmat(nan_list(L,1),1,2), ... repmat(nan_list(L,1),1,2)+repmat([n, 0],nl,1), ... repmat([1 -1],nl,1),n*m,n*m); end % +1,+1 L = find((nan_list(:,2) < n) & (nan_list(:,3) < m)); nl=length(L); if nl>0 fda=fda+sparse(repmat(nan_list(L,1),1,2), ... repmat(nan_list(L,1),1,2)+repmat([n+1, 0],nl,1), ... repmat([1 -1],nl,1),n*m,n*m); end % eliminate knowns rhs=-fda(:,known_list)*A(known_list); % and solve... B=A; k=nan_list(:,1); B(k)=fda(k,k)\rhs(k); end % all done, make sure that B is the same shape as % A was when we came in. B=reshape(B,n,m); % ==================================================== % end of main function % ==================================================== % ==================================================== % begin subfunctions % ==================================================== function neighbors_list=identify_neighbors(n,m,nan_list,talks_to) % identify_neighbors: identifies all the neighbors of % those nodes in nan_list, not including the nans % themselves % % arguments (input): % n,m - scalar - [n,m]=size(A), where A is the % array to be interpolated % nan_list - array - list of every nan element in A % nan_list(i,1) == linear index of i'th nan element % nan_list(i,2) == row index of i'th nan element % nan_list(i,3) == column index of i'th nan element % talks_to - px2 array - defines which nodes communicate % with each other, i.e., which nodes are neighbors. % % talks_to(i,1) - defines the offset in the row % dimension of a neighbor % talks_to(i,2) - defines the offset in the column % dimension of a neighbor % % For example, talks_to = [-1 0;0 -1;1 0;0 1] % means that each node talks only to its immediate % neighbors horizontally and vertically. % % arguments(output): % neighbors_list - array - list of all neighbors of % all the nodes in nan_list if ~isempty(nan_list) % use the definition of a neighbor in talks_to nan_count=size(nan_list,1); talk_count=size(talks_to,1); nn=zeros(nan_count*talk_count,2); j=[1,nan_count]; for i=1:talk_count nn(j(1):j(2),:)=nan_list(:,2:3) + ... repmat(talks_to(i,:),nan_count,1); j=j+nan_count; end % drop those nodes which fall outside the bounds of the % original array L = (nn(:,1)<1)|(nn(:,1)>n)|(nn(:,2)<1)|(nn(:,2)>m); nn(L,:)=[]; % form the same format 3 column array as nan_list neighbors_list=[sub2ind([n,m],nn(:,1),nn(:,2)),nn]; % delete replicates in the neighbors list neighbors_list=unique(neighbors_list,'rows'); % and delete those which are also in the list of NaNs. neighbors_list=setdiff(neighbors_list,nan_list,'rows'); else neighbors_list=[]; end \ No newline at end of file diff --git a/ProgramFiles/v4.2/Utilities/inpaint_nans_bc.m b/ProgramFiles/v4.2/Utilities/inpaint_nans_bc.m new file mode 100755 index 0000000..d6ae57e --- /dev/null +++ b/ProgramFiles/v4.2/Utilities/inpaint_nans_bc.m @@ -0,0 +1 @@ +function B=inpaint_nans_bc(A,method,bcclass) % INPAINT_NANS_BC: in-paints over nans in an array, with spherical or toroidal boundary conditions % usage: B=inpaint_nsns_bc(A) % default method % usage: B=inpaint_nsns_bc(A,method) % specify method used % usage: B=inpaint_nsns_bc(A,method,bcclass) % specify class of boundary conditions applied % % Solves approximation to one of several pdes to % interpolate and extrapolate holes in an array. % Depending upon the boundary conditions specified, % the array will effectively be treated as if it lies % on either the surface of a sphere or a toroid. % % arguments (input): % A - nxm array with some NaNs to be filled in % % method - (OPTIONAL) scalar numeric flag - specifies % which approach (or physical metaphor to use % for the interpolation.) All methods are capable % of extrapolation, some are better than others. % There are also speed differences, as well as % accuracy differences for smooth surfaces. % % The methods employed here are a subset of the % methods of the original inpaint_nans. % % methods {0,1} use a simple plate metaphor. % method 4 uses a spring metaphor. % % method == 0 --> (DEFAULT) see method 1, but % this method does not build as large of a % linear system in the case of only a few % NaNs in a large array. % Extrapolation behavior is linear. % % method == 1 --> simple approach, applies del^2 % over the entire array, then drops those parts % of the array which do not have any contact with % NaNs. Uses a least squares approach, but it % does not modify known values. % In the case of small arrays, this method is % quite fast as it does very little extra work. % Extrapolation behavior is linear. % % method == 4 --> Uses a spring metaphor. Assumes % springs (with a nominal length of zero) % connect each node with every neighbor % (horizontally, vertically and diagonally) % Since each node tries to be like its neighbors, % extrapolation is as a constant function where % this is consistent with the neighboring nodes. % % DEFAULT: 0 % % bcclass - (OPTIONAL) character flag, indicating how % the array boundaries will be treated in the % inpainting operation. bcclass may be either % 'sphere' or 'toroid', or any simple contraction % of these words. % % bcclass = 'sphere' --> The first and last rows % of the array will be treated as if they are % at the North and South poles of a sphere. % Adjacent to those rows will be singular % phantom nodes at each pole. % % bcclass = 'toroid' --> The first and last rows % of the array will be treated as if they are % adjacent to ech other. As well, the first and % last columns will be adjacent to each other. % % DEFAULT: 'sphere' % % arguments (output): % B - nxm array with NaNs replaced % % % Example: % [x,y] = meshgrid(0:.01:1); % z0 = exp(x+y); % znan = z0; % znan(20:50,40:70) = NaN; % znan(30:90,5:10) = NaN; % znan(70:75,40:90) = NaN; % % z = inpaint_nans(znan); % % % See also: griddata, interp1 % % Author: John D'Errico % e-mail address: woodchips@rochester.rr.com % Release: 2 % Release date: 4/15/06 % I always need to know which elements are NaN, % and what size the array is for any method [n,m]=size(A); A=A(:); nm=n*m; k=isnan(A(:)); % list those nodes which are known, and which will % be interpolated nan_list=find(k); known_list=find(~k); % how many nans overall nan_count=length(nan_list); % convert NaN indices to (r,c) form % nan_list==find(k) are the unrolled (linear) indices % (row,column) form [nr,nc]=ind2sub([n,m],nan_list); % both forms of index in one array: % column 1 == unrolled index % column 2 == row index % column 3 == column index nan_list=[nan_list,nr,nc]; % supply default method if (nargin<2) || isempty(method) method = 0; elseif ~ismember(method,[0 1 4]) error('INPAINT_NANS_BC:improperargument', ... 'If supplied, method must be one of: {0,1,4}.') end % supply default value for bcclass if (nargin < 3) || isempty(bcclass) bcclass = 'sphere'; elseif ~ischar(bcclass) error('INPAINT_NANS_BC:improperargument', ... 'If supplied, bcclass must be ''sphere'' or ''toroid''') else % it was a character string valid = {'sphere' 'toroid'}; % check to see if it is valid [bcclass,errorclass] = validstring(arg,valid); if ~isempty(errorclass) error('INPAINT_NANS_BC:improperargument', ... 'If supplied, bcclass must be ''sphere'' or ''toroid''') end end % choice of methods switch method case 0 % The same as method == 1, except only work on those % elements which are NaN, or at least touch a NaN. % horizontal and vertical neighbors only talks_to = [-1 0;0 -1;1 0;0 1]; neighbors_list=identify_neighbors(n,m,nan_list,talks_to); % list of all nodes we have identified all_list=[nan_list;neighbors_list]; % generate sparse array with second partials on row % variable for each element in either list, but only % for those nodes which have a row index > 1 or < n L = find((all_list(:,2) > 1) & (all_list(:,2) < n)); nl=length(L); if nl>0 fda=sparse(repmat(all_list(L,1),1,3), ... repmat(all_list(L,1),1,3)+repmat([-1 0 1],nl,1), ... repmat([1 -2 1],nl,1),nm,nm); else fda=spalloc(n*m,n*m,size(all_list,1)*5); end % 2nd partials on column index L = find((all_list(:,3) > 1) & (all_list(:,3) < m)); nl=length(L); if nl>0 fda=fda+sparse(repmat(all_list(L,1),1,3), ... repmat(all_list(L,1),1,3)+repmat([-n 0 n],nl,1), ... repmat([1 -2 1],nl,1),nm,nm); end % eliminate knowns rhs=-fda(:,known_list)*A(known_list); k=find(any(fda(:,nan_list(:,1)),2)); % and solve... B=A; B(nan_list(:,1))=fda(k,nan_list(:,1))\rhs(k); case 1 % least squares approach with del^2. Build system % for every array element as an unknown, and then % eliminate those which are knowns. % Build sparse matrix approximating del^2 for % every element in A. % Compute finite difference for second partials % on row variable first [i,j]=ndgrid(1:n,1:m); ind=i(:)+(j(:)-1)*n; np=n*m; switch bcclass case 'sphere' % we need to have two phantom nodes at the poles np = np + 2; end fda=sparse(repmat(ind,1,3),[ind-1,ind,ind+1], ... repmat([1 -2 1],np,1),n*m,n*m); % now second partials on column variable [i,j]=ndgrid(1:n,2:(m-1)); ind=i(:)+(j(:)-1)*n; np=n*(m-2); fda=fda+sparse(repmat(ind,1,3),[ind-n,ind,ind+n], ... repmat([1 -2 1],np,1),nm,nm); % eliminate knowns rhs=-fda(:,known_list)*A(known_list); k=find(any(fda(:,nan_list),2)); % and solve... B=A; B(nan_list(:,1))=fda(k,nan_list(:,1))\rhs(k); case 4 % Spring analogy % interpolating operator. % list of all springs between a node and a horizontal % or vertical neighbor hv_list=[-1 -1 0;1 1 0;-n 0 -1;n 0 1]; hv_springs=[]; for i=1:4 hvs=nan_list+repmat(hv_list(i,:),nan_count,1); k=(hvs(:,2)>=1) & (hvs(:,2)<=n) & (hvs(:,3)>=1) & (hvs(:,3)<=m); hv_springs=[hv_springs;[nan_list(k,1),hvs(k,1)]]; end % delete replicate springs hv_springs=unique(sort(hv_springs,2),'rows'); % build sparse matrix of connections, springs % connecting diagonal neighbors are weaker than % the horizontal and vertical springs nhv=size(hv_springs,1); springs=sparse(repmat((1:nhv)',1,2),hv_springs, ... repmat([1 -1],nhv,1),nhv,nm); % eliminate knowns rhs=-springs(:,known_list)*A(known_list); % and solve... B=A; B(nan_list(:,1))=springs(:,nan_list(:,1))\rhs; end % all done, make sure that B is the same shape as % A was when we came in. B=reshape(B,n,m); end % mainline % ==================================================== % end of main function % ==================================================== % ==================================================== % begin subfunctions % ==================================================== function neighbors_list=identify_neighbors(n,m,nan_list,talks_to) % identify_neighbors: identifies all the neighbors of % those nodes in nan_list, not including the nans % themselves % % arguments (input): % n,m - scalar - [n,m]=size(A), where A is the % array to be interpolated % nan_list - array - list of every nan element in A % nan_list(i,1) == linear index of i'th nan element % nan_list(i,2) == row index of i'th nan element % nan_list(i,3) == column index of i'th nan element % talks_to - px2 array - defines which nodes communicate % with each other, i.e., which nodes are neighbors. % % talks_to(i,1) - defines the offset in the row % dimension of a neighbor % talks_to(i,2) - defines the offset in the column % dimension of a neighbor % % For example, talks_to = [-1 0;0 -1;1 0;0 1] % means that each node talks only to its immediate % neighbors horizontally and vertically. % % arguments(output): % neighbors_list - array - list of all neighbors of % all the nodes in nan_list if ~isempty(nan_list) % use the definition of a neighbor in talks_to nan_count=size(nan_list,1); talk_count=size(talks_to,1); nn=zeros(nan_count*talk_count,2); j=[1,nan_count]; for i=1:talk_count nn(j(1):j(2),:)=nan_list(:,2:3) + ... repmat(talks_to(i,:),nan_count,1); j=j+nan_count; end % form the same format 3 column array as nan_list neighbors_list=[sub2ind([n,m],nn(:,1),nn(:,2)),nn]; % delete replicates in the neighbors list neighbors_list=unique(neighbors_list,'rows'); % and delete those which are also in the list of NaNs. neighbors_list=setdiff(neighbors_list,nan_list,'rows'); else neighbors_list=[]; end end % function identify_neighbors function [str,errorclass] = validstring(arg,valid) % validstring: compares a string against a set of valid options % usage: [str,errorclass] = validstring(arg,valid) % % If a direct hit, or any unambiguous shortening is found, that % string is returned. Capitalization is ignored. % % arguments: (input) % arg - character string, to be tested against a list % of valid choices. Capitalization is ignored. % % valid - cellstring array of alternative choices % % Arguments: (output) % str - string - resulting choice resolved from the % list of valid arguments. If no unambiguous % choice can be resolved, then str will be empty. % % errorclass - string - A string argument that explains % the error. It will be one of the following % possibilities: % % '' --> No error. An unambiguous match for arg % was found among the choices. % % 'No match found' --> No match was found among % the choices provided in valid. % % 'Ambiguous argument' --> At least two ambiguous % matches were found among those provided % in valid. % % % Example: % valid = {'off' 'on' 'The sky is falling'} % % % See also: parse_pv_pairs, strmatch, strcmpi % % Author: John D'Errico % e-mail: woodchips@rochester.rr.com % Release: 1.0 % Release date: 3/25/2010 ind = strmatch(lower(arg),lower(valid)); if isempty(ind) % No hit found errorclass = 'No match found'; str = ''; elseif (length(ind) > 1) % Ambiguous arg, hitting more than one of the valid options errorclass = 'Ambiguous argument'; str = ''; return else errorclass = ''; str = valid{ind}; end end % function validstring \ No newline at end of file diff --git a/ProgramFiles/v4.2/Utilities/ppdiff.m b/ProgramFiles/v4.2/Utilities/ppdiff.m new file mode 100644 index 0000000..991464c --- /dev/null +++ b/ProgramFiles/v4.2/Utilities/ppdiff.m @@ -0,0 +1,52 @@ +function qq = ppdiff(pp,j) +%PPDIFF Differentiate piecewise polynomial. +% QQ = PPDIFF(PP,J) returns the J:th derivative of a piecewise +% polynomial PP. PP must be on the form evaluated by PPVAL. QQ is a +% piecewise polynomial on the same form. Default value for J is 1. +% +% Example: +% x = linspace(-pi,pi,9); +% y = sin(x); +% pp = spline(x,y); +% qq = ppdiff(pp); +% xx = linspace(-pi,pi,201); +% plot(xx,cos(xx),'b',xx,ppval(qq,xx),'r') +% +% See also PPVAL, SPLINE, SPLINEFIT, PPINT + +% Author: Jonas Lundgren 2009 + +if nargin < 1, help ppdiff, return, end +if nargin < 2, j = 1; end + +% Check diff order +if ~isreal(j) || mod(j,1) || j < 0 + msgid = 'PPDIFF:DiffOrder'; + message = 'Order of derivative must be a non-negative integer!'; + error(msgid,message) +end + +% Get coefficients +coefs = pp.coefs; +[m n] = size(coefs); + +if j == 0 + % Do nothing +elseif j < n + % Derivative of order J + D = [n-j:-1:1; ones(j-1,n-j)]; + D = cumsum(D,1); + D = prod(D,1); + coefs = coefs(:,1:n-j); + for k = 1:n-j + coefs(:,k) = D(k)*coefs(:,k); + end +else + % Derivative kills PP + coefs = zeros(m,1); +end + +% Set output +qq = pp; +qq.coefs = coefs; +qq.order = size(coefs,2); \ No newline at end of file diff --git a/ProgramFiles/v4.2/Utilities/waypoints_along_gait.m b/ProgramFiles/v4.2/Utilities/waypoints_along_gait.m new file mode 100644 index 0000000..b84fe3e --- /dev/null +++ b/ProgramFiles/v4.2/Utilities/waypoints_along_gait.m @@ -0,0 +1,23 @@ +function [alpha,t] = waypoints_along_gait(gaitname,n) +% Get a set of n waypoints along the gait. Gaitname should be the filename +% for the gait, including the shchf_ heading, and the .mat extension can be +% included or neglected without affecting the function of the program + +% For now, this program assumes that there is only a single gait, with one +% segment (and works with segment 1 of gait 1 if more are defined in the +% file). + + % Get the paths for the sysplotter directories + load sysplotter_config + + % Load the .mat file for the gait + load(fullfile(datapath,gaitname)) + + % Generate n points evenly spaced along the gait + time_range = p.time_def{1}{1}; + t = linspace(time_range(1), time_range(2),n); + + % Find the gait points at these times + alpha = p.phi_def{1}{1}(t); + +end \ No newline at end of file diff --git a/ProgramFiles/v4.2/sys_calcpath.m b/ProgramFiles/v4.2/sys_calcpath.m index e280d25..4b5c6c2 100644 --- a/ProgramFiles/v4.2/sys_calcpath.m +++ b/ProgramFiles/v4.2/sys_calcpath.m @@ -45,6 +45,7 @@ %Vectorize the shape change equation if it's an inline function p = vectorize_shch(p); %#ok + save(infile2,'p') %get shape change loci p = find_loci(s,p); diff --git a/ProgramFiles/v4.2/sys_calcpath_fcns/vectorize_shch.m~ b/ProgramFiles/v4.2/sys_calcpath_fcns/vectorize_shch.m~ deleted file mode 100644 index 6bcb316..0000000 --- a/ProgramFiles/v4.2/sys_calcpath_fcns/vectorize_shch.m~ +++ /dev/null @@ -1,216 +0,0 @@ -%Vectorize the shape change equation -function p = vectorize_shch(p) - - pathspec = {'phi_def',; - 'dphi_def'; - 'time_def'; - 'phi_marker'; - 'phi_arrows' - 'phi_res' - 'cBVI_method'}; - - - default_values = {[]; % Placeholder - []; % Placeholder - [0 2*pi]; % default timespan - []; % no marker - 0; % no arrows - 30; % 30 points along the path - 'none'}; % default to no area integration for cBVI - %%%%%%%%% - % First, make sure that there is something in each of the pathspec - % fields - for i = 1:length(pathspec) - - % Make sure that phi_def is the first element in pathspec, - % otherwise the handling of the other elements could break - if i == 1 && ~strcmp(pathspec{i},'phi_def') - error('Please make sure that phi_def is the first element of the pathspec list at the top of this file') - end - - switch pathspec{i} - - - case 'phi_def' - - % check for absence emptiness of field - if ~isfield(p,pathspec{i}) || isempty(p.(pathspec{i})) - - % error if no shape change is actually specified - error('Shape change has no function specified for phi_def. Use the No Shapechange menu option to generate a plot with no shape change overlaid') - - else % make sure phi_def is a double-nested cell function - - p.phi_def = doublewrap(p.phi_def,'function_handle'); - - end - - % If shape change derivative is not defined, use a - % numerical approximation of it - case 'dphi_def' - - % check for absence emptiness of field - if ~isfield(p,pathspec{i}) || isempty(p.(pathspec{i})) - - % Create a double-wrapped set of numerical derivatives - % If creating all the derivatives, do not warn user - % about them - for j = 1:numel(p.phi_def) - for k = 1:numel(p.phi_def{j}) - p.dphi_def{j}{k} = @(T) jacobianest(@(TT) p.phi_def{j}{k}(TT),T); - end - end - - else - - % Fill in any unspecified derivatives with numerical - % derivatives - - % first, strokes and sections that were only partly - % addressed - for j = 1:min(numel(p.dphi_def),numel(p.phi_def)) % first fill in partially-defined derivatives - - if numel(p.dphi_def) < numel(p.phi_def) - warning('Some (but not all) shape changes *segments* have analytically specified derivatives. The remainder have been filled in numerically; fully specify the dphi_def analytically to improve performance, or specify an a numerical derivative function in the shape change file to suppress this warning') - end - - for k = (numel(p.dphi_def{j})+1):numel(p.phi_def{j}) - - - p.dphi_def{j}{k} = @(T) jacobianest(@(TT) p.phi_def{j}{k}(TT),T); - - end - - end - - % second, strokes that were completely unspecified - if numel(p.dphi_def) < numel(p.phi_def) - warning('Some (but not all) shape changes have analytically specified derivatives. The remainder have been filled in numerically; fully specify the dphi_def analytically to improve performance, or specify an a numerical derivative function in the shape change file to suppress this warning') - end - - for j = (numel(p.dphi_def)+1) : numel(p.phi_def) % first fill in partially-defined derivatives - - for k = 1:numel(p.phi_def{j}) - - p.dphi_def{j}{k} = @(T) jacobianest(@(TT) p.phi_def{j}{k}(TT),T); - - end - - end - - - - - end - - otherwise - - % check for absence emptiness of field - if ~isfield(p,pathspec{i}) || isempty(p.(pathspec{i})) - - % Create a single instance of the default object - p.(pathspec{i}) = default_values{i}; - - end - - % Ensure double-wrapping of field - p.(pathspec{i}) = doublewrap(p.(pathspec{i}),class(default_values{i})); - - %%%% - % Match the dimension of the field to that of the phi_def - % function - - % number of paths here and in phi_def - n_thisfield = numel(p.(pathspec{i})); - n_def = numel(p.phi_def); - - if n_thisfield < n_def; - - if n_thisfield ~= 1 - warning(['Fewer ' pathspec{i} ' entries than phi_defs. Replicating the last entry']) - end - - % Replicate last entry into needed paths - p.(pathspec{i})( (n_thisfield+1) : n_def ) ... - = repmat(p.(pathspec{i})(n_thisfield),[n_def-n_thisfield,1]); - - end - - % Loop over the paths, matching the segment counts - for j = 1:n_def - - n_seg_def = numel(p.phi_def{j}); - n_seg_this = numel(p.(pathspec{i}){j}); - - if n_seg_this < n_seg_def - - if n_seg_this ~= 1 - warning(['Fewer ' pathspec{i} ' segment entries than phi_def segments. Replicating the last entry']) - end - - % Replicate last entry into needed paths - p.(pathspec{i}){j}( (n_seg_this+1) : n_seg_def ) ... - = repmat(p.(pathspec{i}){j}(n_seg_this),[n_seg_def-n_seg_this,1]); - - end - - end - - - - end - - end - -end - -function B = doublewrap(A,classname) -% take object A, verify that its core is of class CLASSNAME, and nest it inside -% cell functions sufficient to make it doulbe-wrapped to return as B - - % test for one layer deep - if iscell(A) - - % test for two layers deep -- make sure all the contents of A are - % cells - if all(cellfun(@(x) iscell(x),A)) - - % if properly nested, do nothing - - else - - % if it's only one cell deep, add a cell layer around it - A = {A}; - - end - - else - - % if it's not a cell array, make it a double layered one - A = {{A}}; - - end - - - % Check to make sure all contents are of the right class - class_test = zeros(numel(A),1); - for i = 1:numel(A) - class_test(i) = all(cellfun(@(x) isa(x,classname),A{i})); - end - - if all(class_test) - - % If all meet the specified condition, do nothing - - else - - % If it's not the right class, throw an error - error(['The core of the input object should be of class ' classname ', but is of class ' class(A) '.']) - - end - - B = A; - - -end - diff --git a/ProgramFiles/v4.2/sys_calcsystem_fcns/four_way_symmetry.m b/ProgramFiles/v4.2/sys_calcsystem_fcns/four_way_symmetry.m index b6addae..f445db9 100644 --- a/ProgramFiles/v4.2/sys_calcsystem_fcns/four_way_symmetry.m +++ b/ProgramFiles/v4.2/sys_calcsystem_fcns/four_way_symmetry.m @@ -110,7 +110,57 @@ f.theta2.original = At2; f.theta2.reflect_even = At2; f.theta2.reflect_odd = At2; - f.theta2.rotate_180 = At2; + f.theta2.rotate_180 = At2; + + case 'None' + + % Shape reflections + f.a1.original = a1; + f.a1.reflect_even = a1; + f.a1.reflect_odd = a1; + f.a1.rotate_180 = a1; + + f.a2.original = a2; + f.a2.reflect_even = a2; + f.a2.reflect_odd = a2; + f.a2.rotate_180 = a2; + + % X field + f.x1.original = Axa1; + f.x1.reflect_even = Axa1; + f.x1.reflect_odd = Axa1; + f.x1.rotate_180 = Axa1; + + f.x2.original = Axa2; + f.x2.reflect_even = Axa2; + f.x2.reflect_odd = Axa2; + f.x2.rotate_180 = Axa2; + + % Y field + f.y1.original = Aya1; + f.y1.reflect_even = Aya1; + f.y1.reflect_odd = Aya1; + f.y1.rotate_180 = Aya1; + + f.y2.original = Aya2; + f.y2.reflect_even = Aya2; + f.y2.reflect_odd = Aya2; + f.y2.rotate_180 = Aya2; + + % Theta field + f.theta1.original = At1; + f.theta1.reflect_even = At1; + f.theta1.reflect_odd = At1; + f.theta1.rotate_180 = At1; + + f.theta2.original = At2; + f.theta2.reflect_even = At2; + f.theta2.reflect_odd = At2; + f.theta2.rotate_180 = At2; + + otherwise + + error('Invalid symmetry averaging mode') end diff --git a/ProgramFiles/v4.2/sys_calcsystem_fcns/optimize_coordinate_choice.m b/ProgramFiles/v4.2/sys_calcsystem_fcns/optimize_coordinate_choice.m index 7a30019..5930f74 100644 --- a/ProgramFiles/v4.2/sys_calcsystem_fcns/optimize_coordinate_choice.m +++ b/ProgramFiles/v4.2/sys_calcsystem_fcns/optimize_coordinate_choice.m @@ -74,33 +74,62 @@ %%%%% %Zero the potential function for theta at the center of the domain - % If the shapespace origin is within the domain, use that as E_theta=0. + % If the shapespace origin (zero or unity) is within the domain, use that as E_theta=0. % Otherwise, use the value at the center of the domain - - % Create a mask to turn all the lower-bounds negative - n_dim = numel(size(E_theta)); - mask = ones(1,2*n_dim); - mask(1:2:2*n_dim) = -1; - within = (mask .* s.grid_range)>0; - - if all(within) - - zcell = num2cell(zeros(n_dim,1)); - - E_theta_offset = interpn(s.grid.eval{:},E_theta,zcell{:},'cubic'); - - else +% +% % Create a mask to turn all the lower-bounds negative +% n_dim = numel(size(E_theta)); +% mask = ones(1,2*n_dim); +% mask(1:2:2*n_dim) = -1; +% within = (mask .* s.grid_range)>0; +% +% if all(within) +% +% zcell = num2cell(zeros(n_dim,1)); +% +% E_theta_offset = interpn(s.grid.eval{:},E_theta,zcell{:},'cubic'); +% +% else +% +% % Loop over the dimensions, getting the center values in each dimension +% centervals = zeros(n_dim,2); +% for i = 1:n_dim +% centervals(i,:) = E_theta([floor((size(E_theta,i)+1)/2) ceil((size(E_theta,i)+1)/2)]); +% end +% % Mean these value to get the value at the center of the space +% E_theta_offset = mean(centervals(:)); +% +% end + + n_dim = numel(size(E_theta)); + centervals = num2cell(zeros(n_dim,1)); + + for idx = 1:n_dim; + + dim_range = s.grid_range([(2*idx)-1,2*idx]); + + % If zero is in the range for a given dimension, center on that value + if (dim_range(1) <= 0) && (dim_range(2) >= 0) + + centervals{idx} = 0; + + % If one (unity) is in the range, center on that + elseif (dim_range(1) <= 1) && (dim_range(2) >= 1) + + centervals{idx} = 1; + + % Otherwise, put the center of the E_theta function at the middle + % of the grid range + else + + centervals{idx} = mean(dim_range); + + end + + end + + E_theta_offset = interpn(s.grid.eval{:},E_theta,centervals{:},'cubic'); - % Loop over the dimensions, getting the center values in each dimension - centervals = zeros(n_dim,2); - for i = 1:n_dim - centervals(i,:) = E_theta([floor((size(E_theta,i)+1)/2) ceil((size(E_theta,i)+1)/2)]); - end - % Mean these value to get the value at the center of the space - E_theta_offset = mean(centervals(:)); - - end - E_theta = E_theta-E_theta_offset; diff --git a/ProgramFiles/v4.2/sys_draw.m b/ProgramFiles/v4.2/sys_draw.m index c653262..9bc4410 100644 --- a/ProgramFiles/v4.2/sys_draw.m +++ b/ProgramFiles/v4.2/sys_draw.m @@ -1,4 +1,4 @@ -function sys_draw(plot_structure,sys,shch,progress,update,resolution,handles) +function plot_info = sys_draw(plot_structure,sys,shch,progress,update,resolution,handles) %make sure plot data file is up to date if update @@ -40,7 +40,7 @@ function sys_draw(plot_structure,sys,shch,progress,update,resolution,handles) eval(['plot_command = @' plot_structure(i).category '_draw;']); %call that plot command - plot_info = plot_command(s,p,plot_structure(i),sys,shch,convert,resolution); + plot_info(i,1) = plot_command(s,p,plot_structure(i),sys,shch,convert,resolution); end diff --git a/ProgramFiles/v4.2/sys_draw_fcns/overlay_shape_change_2d.m b/ProgramFiles/v4.2/sys_draw_fcns/overlay_shape_change_2d.m index b8b7f5a..e11a760 100644 --- a/ProgramFiles/v4.2/sys_draw_fcns/overlay_shape_change_2d.m +++ b/ProgramFiles/v4.2/sys_draw_fcns/overlay_shape_change_2d.m @@ -18,7 +18,12 @@ function overlay_shape_change_2d(ax,p,convert) %draw the path itself - line(p.phi_locus_full{i}.shape(:,1),p.phi_locus_full{i}.shape(:,2),'Color',Colorset.spot,'LineWidth',5,'Parent',ax); + + for idx = 1:size(p.phi_locus_full{i}.shape,2) + pathpoints{idx} = p.phi_locus_full{i}.shape(:,idx); + end + + line(pathpoints{:},'Color',Colorset.spot,'LineWidth',5,'Parent',ax); %draw the direction arrows on the line plot_dir_arrows(p.phi_locus_full{i}.shape(:,1),p.phi_locus_full{i}.shape(:,2),p.phi_arrows{i}{1},'Color',Colorset.spot,'LineWidth',4,'Parent',ax); diff --git a/ProgramFiles/v4.2/sys_draw_fcns/xy_draw_helper.m b/ProgramFiles/v4.2/sys_draw_fcns/xy_draw_helper.m index 0555490..01cd97c 100644 --- a/ProgramFiles/v4.2/sys_draw_fcns/xy_draw_helper.m +++ b/ProgramFiles/v4.2/sys_draw_fcns/xy_draw_helper.m @@ -147,9 +147,11 @@ % Set the axes with a nice buffering [x_min,x_max,y_min,y_max] = ... set_axis_limits(ax,x_collect,y_collect,.07,.07); - end + end + set(ax,'ZLim',10*[-1 1]); % prevent layering from going outside zlim axis(ax,'equal'); new_lim = [get(ax,'xlim') get(ax,'ylim')]; + %%%%%%% diff --git a/ProgramFiles/v4.2/sysplotter.fig b/ProgramFiles/v4.2/sysplotter.fig index 8b8d422..b2b86d1 100644 Binary files a/ProgramFiles/v4.2/sysplotter.fig and b/ProgramFiles/v4.2/sysplotter.fig differ diff --git a/ProgramFiles/v4.2/sysplotter.m b/ProgramFiles/v4.2/sysplotter.m index 89eb4c5..7ff0350 100644 --- a/ProgramFiles/v4.2/sysplotter.m +++ b/ProgramFiles/v4.2/sysplotter.m @@ -22,12 +22,14 @@ % Edit the above text to modify the response to help sysplotter -% Last Modified by GUIDE v2.5 02-Dec-2015 17:23:09 +% Last Modified by GUIDE v2.5 13-Jul-2016 14:49:31 + + addpath('./Utilities') %path to gui functions - addpath(genpath('sysplotter_gui_fcns'), genpath('sys_calcpath_fcns'),... - genpath('sys_calcsystem_fcns'), genpath('sys_draw_fcns'), ... - genpath('sys_update_fcns'),genpath('sysplotter_config_fcns'),genpath('Utilities')) + addpath(genpath(GetFullPath('sysplotter_gui_fcns')), genpath(GetFullPath('sys_calcpath_fcns')),... + genpath(GetFullPath('sys_calcsystem_fcns')), genpath(GetFullPath('sys_draw_fcns')), ... + genpath(GetFullPath('sys_update_fcns')),genpath(GetFullPath('sysplotter_config_fcns')),genpath(GetFullPath('Utilities'))) %%% % Ensure that system files are properly accessible diff --git a/ProgramFiles/v4.2/sysplotter_config.m b/ProgramFiles/v4.2/sysplotter_config.m index 37c4374..919c2b2 100644 --- a/ProgramFiles/v4.2/sysplotter_config.m +++ b/ProgramFiles/v4.2/sysplotter_config.m @@ -68,7 +68,7 @@ function sysplotter_config_OpeningFcn(hObject, eventdata, handles, varargin) %#o addpath('Utilities/gait_gui_draw') % Default values for HH and refpoint paths - set(handles.inputpathconfig,'String','../../UserFiles/') + set(handles.inputpathconfig,'String',GetFullPath('../../UserFiles/v4/GenericUser')) set(handles.HodgeHelmholtzconfig,'String','HodgeHelmholtz') set(handles.Refpointconfig,'String','RefPointOptimizer') set(handles.Colorconfig,'String','sys_draw_fcns/colorsets/color_Red.m') diff --git a/ProgramFiles/v4.2/sysplotter_gui_fcns/gait_gui_draw_Callback.m b/ProgramFiles/v4.2/sysplotter_gui_fcns/gait_gui_draw_Callback.m new file mode 100644 index 0000000..8ab89bd --- /dev/null +++ b/ProgramFiles/v4.2/sysplotter_gui_fcns/gait_gui_draw_Callback.m @@ -0,0 +1,13 @@ +% --- Executes on button press in gait_gui_draw. +function gait_gui_draw_Callback(hObject, eventdata, handles) +% hObject handle to gait_gui_draw (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% Plot whatever gaits are selected in the third checkbox column +plot_info = plotpushbutton_Callback(findall(0,'tag','plotpushbutton3'), eventdata, handles); + +% Execute the gait_gui_draw command + gait_gui_draw(plot_info(1).axes(1)); + +end diff --git a/ProgramFiles/v4.2/sysplotter_gui_fcns/plotpushbutton_Callback.m b/ProgramFiles/v4.2/sysplotter_gui_fcns/plotpushbutton_Callback.m index 828e9a9..fc7955b 100644 --- a/ProgramFiles/v4.2/sysplotter_gui_fcns/plotpushbutton_Callback.m +++ b/ProgramFiles/v4.2/sysplotter_gui_fcns/plotpushbutton_Callback.m @@ -1,5 +1,5 @@ % --- Executes on button press in any plotpushbutton. -function plotpushbutton_Callback(hObject, eventdata, handles) +function plot_info = plotpushbutton_Callback(hObject, eventdata, handles) % hObject handle to plotpushbutton2 (see GCBO) % eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) @@ -64,7 +64,7 @@ function plotpushbutton_Callback(hObject, eventdata, handles) %Call the draw function % test_plot(plots_to_make,current_system,current_shch) -sys_draw(plots_to_make,current_system,current_shch,handles.progresspanel,1,resolution,handles) +plot_info = sys_draw(plots_to_make,current_system,current_shch,handles.progresspanel,1,resolution,handles); %Show full progress bar waitbar2a(1,handles.progresspanel,'Finished plotting') diff --git a/UserFiles/v4/rlhatton/Shape_Changes/optimalgait_serpenoid_varfreq.mat b/UserFiles/v4/rlhatton/Shape_Changes/optimalgait_serpenoid_varfreq.mat new file mode 100644 index 0000000..bcd99ee Binary files /dev/null and b/UserFiles/v4/rlhatton/Shape_Changes/optimalgait_serpenoid_varfreq.mat differ diff --git a/UserFiles/v4/rlhatton/Shape_Changes/params_.mat b/UserFiles/v4/rlhatton/Shape_Changes/params_.mat new file mode 100644 index 0000000..13e5eae Binary files /dev/null and b/UserFiles/v4/rlhatton/Shape_Changes/params_.mat differ diff --git a/UserFiles/v4/rlhatton/Shape_Changes/params_10.mat b/UserFiles/v4/rlhatton/Shape_Changes/params_10.mat new file mode 100644 index 0000000..0484a2f Binary files /dev/null and b/UserFiles/v4/rlhatton/Shape_Changes/params_10.mat differ diff --git a/UserFiles/v4/rlhatton/Shape_Changes/params_107demo.mat b/UserFiles/v4/rlhatton/Shape_Changes/params_107demo.mat new file mode 100644 index 0000000..e890d1b Binary files /dev/null and b/UserFiles/v4/rlhatton/Shape_Changes/params_107demo.mat differ diff --git a/UserFiles/v4/rlhatton/Shape_Changes/params_11.mat b/UserFiles/v4/rlhatton/Shape_Changes/params_11.mat new file mode 100644 index 0000000..93d8c7a Binary files /dev/null and b/UserFiles/v4/rlhatton/Shape_Changes/params_11.mat differ diff --git a/UserFiles/v4/rlhatton/Shape_Changes/params_12.mat b/UserFiles/v4/rlhatton/Shape_Changes/params_12.mat new file mode 100644 index 0000000..8d0ca50 Binary files /dev/null and b/UserFiles/v4/rlhatton/Shape_Changes/params_12.mat differ diff --git a/UserFiles/v4/rlhatton/Shape_Changes/params_13.mat b/UserFiles/v4/rlhatton/Shape_Changes/params_13.mat new file mode 100644 index 0000000..e16d7ad Binary files /dev/null and b/UserFiles/v4/rlhatton/Shape_Changes/params_13.mat differ diff --git a/UserFiles/v4/rlhatton/Shape_Changes/params_6.mat b/UserFiles/v4/rlhatton/Shape_Changes/params_6.mat new file mode 100644 index 0000000..ccf7765 Binary files /dev/null and b/UserFiles/v4/rlhatton/Shape_Changes/params_6.mat differ diff --git a/UserFiles/v4/rlhatton/Shape_Changes/params_7.mat b/UserFiles/v4/rlhatton/Shape_Changes/params_7.mat new file mode 100644 index 0000000..f545d47 Binary files /dev/null and b/UserFiles/v4/rlhatton/Shape_Changes/params_7.mat differ diff --git a/UserFiles/v4/rlhatton/Shape_Changes/params_8.mat b/UserFiles/v4/rlhatton/Shape_Changes/params_8.mat new file mode 100644 index 0000000..6c32d01 Binary files /dev/null and b/UserFiles/v4/rlhatton/Shape_Changes/params_8.mat differ diff --git a/UserFiles/v4/rlhatton/Shape_Changes/params_9.mat b/UserFiles/v4/rlhatton/Shape_Changes/params_9.mat new file mode 100644 index 0000000..f04cf52 Binary files /dev/null and b/UserFiles/v4/rlhatton/Shape_Changes/params_9.mat differ diff --git a/UserFiles/v4/rlhatton/Shape_Changes/params_edgegait.mat b/UserFiles/v4/rlhatton/Shape_Changes/params_edgegait.mat new file mode 100644 index 0000000..a44c9ee Binary files /dev/null and b/UserFiles/v4/rlhatton/Shape_Changes/params_edgegait.mat differ diff --git a/UserFiles/v4/rlhatton/Shape_Changes/params_mudskipper_curved.mat b/UserFiles/v4/rlhatton/Shape_Changes/params_mudskipper_curved.mat new file mode 100644 index 0000000..3f15f2a Binary files /dev/null and b/UserFiles/v4/rlhatton/Shape_Changes/params_mudskipper_curved.mat differ diff --git a/UserFiles/v4/rlhatton/Shape_Changes/params_mudskipper_direct.mat b/UserFiles/v4/rlhatton/Shape_Changes/params_mudskipper_direct.mat new file mode 100644 index 0000000..07da0f7 Binary files /dev/null and b/UserFiles/v4/rlhatton/Shape_Changes/params_mudskipper_direct.mat differ diff --git a/UserFiles/v4/rlhatton/Shape_Changes/params_triangle_wavelength1_innercontour.mat b/UserFiles/v4/rlhatton/Shape_Changes/params_triangle_wavelength1_innercontour.mat new file mode 100644 index 0000000..f70184b Binary files /dev/null and b/UserFiles/v4/rlhatton/Shape_Changes/params_triangle_wavelength1_innercontour.mat differ diff --git a/UserFiles/v4/rlhatton/Shape_Changes/params_triangle_wavelength1_zerocontour.mat b/UserFiles/v4/rlhatton/Shape_Changes/params_triangle_wavelength1_zerocontour.mat new file mode 100644 index 0000000..c2514d5 Binary files /dev/null and b/UserFiles/v4/rlhatton/Shape_Changes/params_triangle_wavelength1_zerocontour.mat differ diff --git a/UserFiles/v4/rlhatton/Shape_Changes/shchf_LowRe_serpeff_forvarfeq.m b/UserFiles/v4/rlhatton/Shape_Changes/shchf_LowRe_serpeff_forvarfeq.m new file mode 100644 index 0000000..384540e --- /dev/null +++ b/UserFiles/v4/rlhatton/Shape_Changes/shchf_LowRe_serpeff_forvarfeq.m @@ -0,0 +1,67 @@ +function output = shchf_LowRe_serpeff_forvarfeq_freqfamily(input_mode) + + % Default argument + if ~exist('input_mode','var') + + input_mode = 'initialize'; + + end + + switch input_mode + + case 'name' + + output = 'Low Re serp eff for varfreq'; + + case 'dependency' + + output.dependency = {}; % the .mat file sharing the same name as this .m file + + case 'initialize' + + %%%% + % Filename to save to + + output = mfilename; + + + %% + %Path definitions + + %path definition + for idx = 1:5 + p.phi_def{idx} = {@(t) strokedef(t,.9+idx/20)}; + end + + %marker locations + p.phi_marker = []; + + %arrows to plot + p.phi_arrows{1} = {2}; + + %time to run path + p.time_def{1} = {[-pi pi]};%{[-1 1]}; + + %path resolution + p.phi_res{1} = {1000}; + + + %%%% + %Save the system properties + output = p; + + end + +end + +function [stroke] = strokedef(t,freq) +load('X_eff_serpenoid.mat') +XV=X; +len=length(X); + t = t(:); + stroke=[fourier_eval(t,XV(1:len/4),XV(len/4+1:len/2),2*pi); fourier_eval(t,XV(len/2+1:3*len/4),XV(3*len/4+1:len),2*pi);freq*ones(size(t'))]'; + + %stroke = [t t]*1.5; + + +end \ No newline at end of file diff --git a/UserFiles/v4/rlhatton/Shape_Changes/shchf_circle_6.m b/UserFiles/v4/rlhatton/Shape_Changes/shchf_circle_6.m new file mode 100644 index 0000000..055a94b --- /dev/null +++ b/UserFiles/v4/rlhatton/Shape_Changes/shchf_circle_6.m @@ -0,0 +1,63 @@ +function output = shchf_circle_6(input_mode) + + % Default argument + if ~exist('input_mode','var') + + input_mode = 'initialize'; + + end + + switch input_mode + + case 'name' + + output = 'Circle Stroke, 6 amplitude'; + + case 'dependency' + + output.dependency = {}; + + case 'initialize' + + %%%% + % Filename to save to + + %% + %Path definitions + + %path definition + p.phi_def{1} = {@strokedef}; + + %marker locations + p.phi_marker = []; + + %arrows to plot + p.phi_arrows{1} = {2}; + + %time to run path + p.time_def{1} = {[0 2*pi]}; + + %p.cBVI_method{1}{1} = 'simple'; + + %path resolution + p.phi_res{1} = {1000}; + + + %%%% + %Output the path properties + output = p; + end + +end + +function [stroke] = strokedef(t) + + t = -t(:)'; + + Rot=sqrt(2)/2*[1 -1;1 1]; + a=6; + + stroke=(Rot*[-a*cos(t);-a*sin(t)])'; + + +end \ No newline at end of file diff --git a/UserFiles/v4/rlhatton/Shape_Changes/shchf_optimalgait_serpenoid_varfreq.m b/UserFiles/v4/rlhatton/Shape_Changes/shchf_optimalgait_serpenoid_varfreq.m new file mode 100644 index 0000000..3feef66 --- /dev/null +++ b/UserFiles/v4/rlhatton/Shape_Changes/shchf_optimalgait_serpenoid_varfreq.m @@ -0,0 +1,84 @@ +function output = shchf_optimalgait_serpenoid_varfreq(input_mode) + + % Default argument + if ~exist('input_mode','var') + + input_mode = 'initialize'; + + end + + % Name the .mat file with the fourier coefficients + paramfile = 'optimalgait_serpenoid_varfreq'; + + + switch input_mode + + case 'name' + + output = 'serpenoid_varfreq'; + + case 'dependency' + + output.dependency = {fullfile(pwd,[paramfile '.mat'])}; + + case 'initialize' + + %%%% + %% + %Path definitions + + % Load the points that the user clicked and the time vector + load(paramfile) + + % Check if start and end are the same point, and set boundary + % conditions accordingly + if alpha1(1)==alpha1(end) && alpha2(1)==alpha2(end) && alpha3(1)==alpha3(end) + splinemode = 'periodic'; + else + splinemode = 'complete'; + end + + % Generate spline structures for the selected points + spline_alpha1 = csape(t,alpha1',splinemode); + spline_alpha2 = csape(t,alpha2',splinemode); + spline_alpha3 = csape(t,alpha3',splinemode); + + % The gait path is now defined by evaluating the two splines at + % specified times + p.phi_def = @(t) [fnval(spline_alpha1,t(:)),fnval(spline_alpha2,t(:)),fnval(spline_alpha3,t(:))]; + + + % Speed up execution by defining the gait velocity as well (so + % the evaluator doesn't need to numerically take the derivative + % at every time + spline_dalpha1 = fnder(spline_alpha1); + spline_dalpha2 = fnder(spline_alpha2); + spline_dalpha3 = fnder(spline_alpha3); + + % The gait path is now defined by evaluating the two splin + % derivatives at specified times + p.dphi_def = @(t) [fnval(spline_dalpha1,t(:)),fnval(spline_dalpha2,t(:)),fnval(spline_dalpha3,t(:))]; + + + %marker locations + p.phi_marker = []; + + %arrows to plot + p.phi_arrows = 2; + + %time to run path + p.time_def = t([1 end]); %#ok + + + %path resolution + p.phi_res = 50; + + + + + %%%% + %Save the shch properties + output = p; + end + +end diff --git a/UserFiles/v4/rlhatton/Shape_Changes/shchf_params_.m b/UserFiles/v4/rlhatton/Shape_Changes/shchf_params_.m new file mode 100644 index 0000000..de9ba6d --- /dev/null +++ b/UserFiles/v4/rlhatton/Shape_Changes/shchf_params_.m @@ -0,0 +1,82 @@ +function output = shchf_params_(input_mode) + + % Default argument + if ~exist('input_mode','var') + + input_mode = 'initialize'; + + end + + % Name the .mat file with the fourier coefficients + paramfile = 'params_'; + + + switch input_mode + + case 'name' + + output = 'params_'; + + case 'dependency' + + output.dependency = {fullfile(pwd,[paramfile '.mat'])}; + + case 'initialize' + + %%%% + %% + %Path definitions + + % Load the points that the user clicked and the time vector + load(paramfile) + + % Check if start and end are the same point, and set boundary + % conditions accordingly + if alpha1(1)==alpha1(end) && alpha2(1)==alpha2(end) + splinemode = 'periodic'; + else + splinemode = 'complete'; + end + + % Generate spline structures for the selected points + spline_alpha1 = csape(t,alpha1',splinemode); + spline_alpha2 = csape(t,alpha2',splinemode); + + % The gait path is now defined by evaluating the two splines at + % specified times + p.phi_def = @(t) [fnval(spline_alpha1,t(:)),fnval(spline_alpha2,t(:))]; + + + % Speed up execution by defining the gait velocity as well (so + % the evaluator doesn't need to numerically take the derivative + % at every time + spline_dalpha1 = fnder(spline_alpha1); + spline_dalpha2 = fnder(spline_alpha2); + + % The gait path is now defined by evaluating the two splin + % derivatives at specified times + p.dphi_def = @(t) [fnval(spline_dalpha1,t(:)),fnval(spline_dalpha2,t(:))]; + + + %marker locations + p.phi_marker = []; + + %arrows to plot + p.phi_arrows = 2; + + %time to run path + p.time_def = t([1 end]); %#ok + + + %path resolution + p.phi_res = 50; + + + + + %%%% + %Save the shch properties + output = p; + end + +end diff --git a/UserFiles/v4/rlhatton/Shape_Changes/shchf_params_10.m b/UserFiles/v4/rlhatton/Shape_Changes/shchf_params_10.m new file mode 100644 index 0000000..d0117df --- /dev/null +++ b/UserFiles/v4/rlhatton/Shape_Changes/shchf_params_10.m @@ -0,0 +1,82 @@ +function output = shchf_params_10(input_mode) + + % Default argument + if ~exist('input_mode','var') + + input_mode = 'initialize'; + + end + + % Name the .mat file with the fourier coefficients + paramfile = 'params_10'; + + + switch input_mode + + case 'name' + + output = 'params_10'; + + case 'dependency' + + output.dependency = {fullfile(pwd,[paramfile '.mat'])}; + + case 'initialize' + + %%%% + %% + %Path definitions + + % Load the points that the user clicked and the time vector + load(paramfile) + + % Check if start and end are the same point, and set boundary + % conditions accordingly + if alpha1(1)==alpha1(end) && alpha2(1)==alpha2(end) + splinemode = 'periodic'; + else + splinemode = 'complete'; + end + + % Generate spline structures for the selected points + spline_alpha1 = csape(t,alpha1',splinemode); + spline_alpha2 = csape(t,alpha2',splinemode); + + % The gait path is now defined by evaluating the two splines at + % specified times + p.phi_def = @(t) [fnval(spline_alpha1,t(:)),fnval(spline_alpha2,t(:))]; + + + % Speed up execution by defining the gait velocity as well (so + % the evaluator doesn't need to numerically take the derivative + % at every time + spline_dalpha1 = fnder(spline_alpha1); + spline_dalpha2 = fnder(spline_alpha2); + + % The gait path is now defined by evaluating the two splin + % derivatives at specified times + p.dphi_def = @(t) [fnval(spline_dalpha1,t(:)),fnval(spline_dalpha2,t(:))]; + + + %marker locations + p.phi_marker = []; + + %arrows to plot + p.phi_arrows = 2; + + %time to run path + p.time_def = t([1 end]); %#ok + + + %path resolution + p.phi_res = 50; + + + + + %%%% + %Save the shch properties + output = p; + end + +end diff --git a/UserFiles/v4/rlhatton/Shape_Changes/shchf_params_107demo.m b/UserFiles/v4/rlhatton/Shape_Changes/shchf_params_107demo.m new file mode 100644 index 0000000..13f4050 --- /dev/null +++ b/UserFiles/v4/rlhatton/Shape_Changes/shchf_params_107demo.m @@ -0,0 +1,98 @@ +function output = shchf_params_107demo(input_mode) + + % Default argument + if ~exist('input_mode','var') + + input_mode = 'initialize'; + + end + + % Name the .mat file with the fourier coefficients + paramfile = 'params_107demo'; + + + switch input_mode + + case 'name' + + output = 'params_107demo'; + + case 'dependency' + + output.dependency = {fullfile(pwd,[paramfile '.mat'])}; + + case 'initialize' + + %%%% + %% + %Path definitions + + % Load the points that the user clicked and the time vector + load(paramfile) + + % Check if start and end are the same point, and set boundary + % conditions accordingly + if alpha1(1)==alpha1(end) && alpha2(1)==alpha2(end) + splinemode = 'periodic'; + else + splinemode = 'complete'; + end + + % Generate spline structures for the selected points + switch splinemode + + case 'periodic' + + % Fit a periodic spline to the selected points; the endslopes are found + % by averaging the positions of the points before and after the join + endslope1 = (alpha1(2)-alpha1(end-1))/(t(end)-t(end-2)); + endslope2 = (alpha2(2)-alpha2(end-1))/(t(end)-t(end-2)); + spline_alpha1 = spline(t,[endslope1;alpha1(:);endslope1]); + spline_alpha2 = spline(t,[endslope2;alpha2(:);endslope2]); + + case 'complete' + + % Fit a non-periodic spline to the selected points + spline_alpha1 = spline(t,alpha1(:)); + spline_alpha2 = spline(t,alpha2(:)); + + end + + % The gait path is now defined by evaluating the two splines at + % specified times + p.phi_def = @(t) [ppval(spline_alpha1,t(:)),ppval(spline_alpha2,t(:))]; + + + % Speed up execution by defining the gait velocity as well (so + % the evaluator doesn't need to numerically take the derivative + % at every time + spline_dalpha1 = ppdiff(spline_alpha1); + spline_dalpha2 = ppdiff(spline_alpha2); + + % The gait path is now defined by evaluating the two splin + % derivatives at specified times + p.dphi_def = @(t) [ppval(spline_dalpha1,t(:)),ppval(spline_dalpha2,t(:))]; + + + %marker locations + p.phi_marker = []; + + %arrows to plot + p.phi_arrows = 2; + + %time to run path + p.time_def = t([1 end]); %#ok + + + %path resolution + p.phi_res = 50; + + + + + %%%% + %Save the shch properties + output = p; + end + +end diff --git a/UserFiles/v4/rlhatton/Shape_Changes/shchf_params_11.m b/UserFiles/v4/rlhatton/Shape_Changes/shchf_params_11.m new file mode 100644 index 0000000..700bba1 --- /dev/null +++ b/UserFiles/v4/rlhatton/Shape_Changes/shchf_params_11.m @@ -0,0 +1,98 @@ +function output = shchf_params_11(input_mode) + + % Default argument + if ~exist('input_mode','var') + + input_mode = 'initialize'; + + end + + % Name the .mat file with the fourier coefficients + paramfile = 'params_11'; + + + switch input_mode + + case 'name' + + output = 'params_11'; + + case 'dependency' + + output.dependency = {fullfile(pwd,[paramfile '.mat'])}; + + case 'initialize' + + %%%% + %% + %Path definitions + + % Load the points that the user clicked and the time vector + load(paramfile) + + % Check if start and end are the same point, and set boundary + % conditions accordingly + if alpha1(1)==alpha1(end) && alpha2(1)==alpha2(end) + splinemode = 'periodic'; + else + splinemode = 'complete'; + end + + % Generate spline structures for the selected points + switch splinemode + + case 'periodic' + + % Fit a periodic spline to the selected points; the endslopes are found + % by averaging the positions of the points before and after the join + endslope1 = (alpha1(2)-alpha1(end-1))/(t(end)-t(end-2)); + endslope2 = (alpha2(2)-alpha2(end-1))/(t(end)-t(end-2)); + spline_alpha1 = spline(t,[endslope1;alpha1(:);endslope1]); + spline_alpha2 = spline(t,[endslope2;alpha2(:);endslope2]); + + case 'complete' + + % Fit a non-periodic spline to the selected points + spline_alpha1 = spline(t,alpha1(:)); + spline_alpha2 = spline(t,alpha2(:)); + + end + + % The gait path is now defined by evaluating the two splines at + % specified times + p.phi_def = @(t) [ppval(spline_alpha1,t(:)),ppval(spline_alpha2,t(:))]; + + + % Speed up execution by defining the gait velocity as well (so + % the evaluator doesn't need to numerically take the derivative + % at every time + spline_dalpha1 = ppdiff(spline_alpha1); + spline_dalpha2 = ppdiff(spline_alpha2); + + % The gait path is now defined by evaluating the two splin + % derivatives at specified times + p.dphi_def = @(t) [ppval(spline_dalpha1,t(:)),ppval(spline_dalpha2,t(:))]; + + + %marker locations + p.phi_marker = []; + + %arrows to plot + p.phi_arrows = 2; + + %time to run path + p.time_def = t([1 end]); %#ok + + + %path resolution + p.phi_res = 50; + + + + + %%%% + %Save the shch properties + output = p; + end + +end diff --git a/UserFiles/v4/rlhatton/Shape_Changes/shchf_params_12.m b/UserFiles/v4/rlhatton/Shape_Changes/shchf_params_12.m new file mode 100644 index 0000000..91f2517 --- /dev/null +++ b/UserFiles/v4/rlhatton/Shape_Changes/shchf_params_12.m @@ -0,0 +1,98 @@ +function output = shchf_params_12(input_mode) + + % Default argument + if ~exist('input_mode','var') + + input_mode = 'initialize'; + + end + + % Name the .mat file with the fourier coefficients + paramfile = 'params_12'; + + + switch input_mode + + case 'name' + + output = 'params_12'; + + case 'dependency' + + output.dependency = {fullfile(pwd,[paramfile '.mat'])}; + + case 'initialize' + + %%%% + %% + %Path definitions + + % Load the points that the user clicked and the time vector + load(paramfile) + + % Check if start and end are the same point, and set boundary + % conditions accordingly + if alpha1(1)==alpha1(end) && alpha2(1)==alpha2(end) + splinemode = 'periodic'; + else + splinemode = 'complete'; + end + + % Generate spline structures for the selected points + switch splinemode + + case 'periodic' + + % Fit a periodic spline to the selected points; the endslopes are found + % by averaging the positions of the points before and after the join + endslope1 = (alpha1(2)-alpha1(end-1))/(t(end)-t(end-2)); + endslope2 = (alpha2(2)-alpha2(end-1))/(t(end)-t(end-2)); + spline_alpha1 = spline(t,[endslope1;alpha1(:);endslope1]); + spline_alpha2 = spline(t,[endslope2;alpha2(:);endslope2]); + + case 'complete' + + % Fit a non-periodic spline to the selected points + spline_alpha1 = spline(t,alpha1(:)); + spline_alpha2 = spline(t,alpha2(:)); + + end + + % The gait path is now defined by evaluating the two splines at + % specified times + p.phi_def = @(t) [ppval(spline_alpha1,t(:)),ppval(spline_alpha2,t(:))]; + + + % Speed up execution by defining the gait velocity as well (so + % the evaluator doesn't need to numerically take the derivative + % at every time + spline_dalpha1 = ppdiff(spline_alpha1); + spline_dalpha2 = ppdiff(spline_alpha2); + + % The gait path is now defined by evaluating the two splin + % derivatives at specified times + p.dphi_def = @(t) [ppval(spline_dalpha1,t(:)),ppval(spline_dalpha2,t(:))]; + + + %marker locations + p.phi_marker = []; + + %arrows to plot + p.phi_arrows = 2; + + %time to run path + p.time_def = t([1 end]); %#ok + + + %path resolution + p.phi_res = 50; + + + + + %%%% + %Save the shch properties + output = p; + end + +end diff --git a/UserFiles/v4/rlhatton/Shape_Changes/shchf_params_13.m b/UserFiles/v4/rlhatton/Shape_Changes/shchf_params_13.m new file mode 100644 index 0000000..8d16a8d --- /dev/null +++ b/UserFiles/v4/rlhatton/Shape_Changes/shchf_params_13.m @@ -0,0 +1,98 @@ +function output = shchf_params_13(input_mode) + + % Default argument + if ~exist('input_mode','var') + + input_mode = 'initialize'; + + end + + % Name the .mat file with the fourier coefficients + paramfile = 'params_13'; + + + switch input_mode + + case 'name' + + output = 'params_13'; + + case 'dependency' + + output.dependency = {fullfile(pwd,[paramfile '.mat'])}; + + case 'initialize' + + %%%% + %% + %Path definitions + + % Load the points that the user clicked and the time vector + load(paramfile) + + % Check if start and end are the same point, and set boundary + % conditions accordingly + if alpha1(1)==alpha1(end) && alpha2(1)==alpha2(end) + splinemode = 'periodic'; + else + splinemode = 'complete'; + end + + % Generate spline structures for the selected points + switch splinemode + + case 'periodic' + + % Fit a periodic spline to the selected points; the endslopes are found + % by averaging the positions of the points before and after the join + endslope1 = (alpha1(2)-alpha1(end-1))/(t(end)-t(end-2)); + endslope2 = (alpha2(2)-alpha2(end-1))/(t(end)-t(end-2)); + spline_alpha1 = spline(t,[endslope1;alpha1(:);endslope1]); + spline_alpha2 = spline(t,[endslope2;alpha2(:);endslope2]); + + case 'complete' + + % Fit a non-periodic spline to the selected points + spline_alpha1 = spline(t,alpha1(:)); + spline_alpha2 = spline(t,alpha2(:)); + + end + + % The gait path is now defined by evaluating the two splines at + % specified times + p.phi_def = @(t) [ppval(spline_alpha1,t(:)),ppval(spline_alpha2,t(:))]; + + + % Speed up execution by defining the gait velocity as well (so + % the evaluator doesn't need to numerically take the derivative + % at every time + spline_dalpha1 = ppdiff(spline_alpha1); + spline_dalpha2 = ppdiff(spline_alpha2); + + % The gait path is now defined by evaluating the two splin + % derivatives at specified times + p.dphi_def = @(t) [ppval(spline_dalpha1,t(:)),ppval(spline_dalpha2,t(:))]; + + + %marker locations + p.phi_marker = []; + + %arrows to plot + p.phi_arrows = 2; + + %time to run path + p.time_def = t([1 end]); %#ok + + + %path resolution + p.phi_res = 50; + + + + + %%%% + %Save the shch properties + output = p; + end + +end diff --git a/UserFiles/v4/rlhatton/Shape_Changes/shchf_params_6.m b/UserFiles/v4/rlhatton/Shape_Changes/shchf_params_6.m new file mode 100644 index 0000000..16424b4 --- /dev/null +++ b/UserFiles/v4/rlhatton/Shape_Changes/shchf_params_6.m @@ -0,0 +1,82 @@ +function output = shchf_params_6(input_mode) + + % Default argument + if ~exist('input_mode','var') + + input_mode = 'initialize'; + + end + + % Name the .mat file with the fourier coefficients + paramfile = 'params_6'; + + + switch input_mode + + case 'name' + + output = 'params_6'; + + case 'dependency' + + output.dependency = {fullfile(pwd,[paramfile '.mat'])}; + + case 'initialize' + + %%%% + %% + %Path definitions + + % Load the points that the user clicked and the time vector + load(paramfile) + + % Check if start and end are the same point, and set boundary + % conditions accordingly + if alpha1(1)==alpha1(end) && alpha2(1)==alpha2(end) + splinemode = 'periodic'; + else + splinemode = 'complete'; + end + + % Generate spline structures for the selected points + spline_alpha1 = csape(t,alpha1',splinemode); + spline_alpha2 = csape(t,alpha2',splinemode); + + % The gait path is now defined by evaluating the two splines at + % specified times + p.phi_def = @(t) [fnval(spline_alpha1,t(:)),fnval(spline_alpha2,t(:))]; + + + % Speed up execution by defining the gait velocity as well (so + % the evaluator doesn't need to numerically take the derivative + % at every time + spline_dalpha1 = fnder(spline_alpha1); + spline_dalpha2 = fnder(spline_alpha2); + + % The gait path is now defined by evaluating the two splin + % derivatives at specified times + p.dphi_def = @(t) [fnval(spline_dalpha1,t(:)),fnval(spline_dalpha2,t(:))]; + + + %marker locations + p.phi_marker = []; + + %arrows to plot + p.phi_arrows = 2; + + %time to run path + p.time_def = t([1 end]); %#ok + + + %path resolution + p.phi_res = 50; + + + + + %%%% + %Save the shch properties + output = p; + end + +end diff --git a/UserFiles/v4/rlhatton/Shape_Changes/shchf_params_7.m b/UserFiles/v4/rlhatton/Shape_Changes/shchf_params_7.m new file mode 100644 index 0000000..1f999dd --- /dev/null +++ b/UserFiles/v4/rlhatton/Shape_Changes/shchf_params_7.m @@ -0,0 +1,84 @@ +function output = shchf_params_7(input_mode) + + % Default argument + if ~exist('input_mode','var') + + input_mode = 'initialize'; + + end + + % Name the .mat file with the fourier coefficients + paramfile = 'params_7'; + + + switch input_mode + + case 'name' + + output = 'params_7'; + + case 'dependency' + + output.dependency = {fullfile(pwd,[paramfile '.mat'])}; + + case 'initialize' + + %%%% + %% + %Path definitions + + % Load the points that the user clicked and the time vector + load(paramfile) + + % Check if start and end are the same point, and set boundary + % conditions accordingly + if alpha1(1)==alpha1(end) && alpha2(1)==alpha2(end) + splinemode = 'periodic'; + else + splinemode = 'complete'; + end + + % Generate spline structures for the selected points + spline_alpha1 = csape(t,alpha1',splinemode); + spline_alpha2 = csape(t,alpha2',splinemode); + + % The gait path is now defined by evaluating the two splines at + % specified times + p.phi_def = @(t_val) [interp1(t,alpha1,t_val,'spline'), interp1(t,alpha2,t_val,'spline')]; + %p.phi_def = @(t) [fnval(spline_alpha1,t(:)),fnval(spline_alpha2,t(:))]; + + + % Speed up execution by defining the gait velocity as well (so + % the evaluator doesn't need to numerically take the derivative + % at every time + spline_dalpha1 = fnder(spline_alpha1); + spline_dalpha2 = fnder(spline_alpha2); + + % The gait path is now defined by evaluating the two splin + % derivatives at specified times + + %p.dphi_def = @(t) [fnval(spline_dalpha1,t(:)),fnval(spline_dalpha2,t(:))]; + + + %marker locations + p.phi_marker = []; + + %arrows to plot + p.phi_arrows = 2; + + %time to run path + p.time_def = t([1 end]); %#ok + + + %path resolution + p.phi_res = 50; + + + + + %%%% + %Save the shch properties + output = p; + end + +end diff --git a/UserFiles/v4/rlhatton/Shape_Changes/shchf_params_8.m b/UserFiles/v4/rlhatton/Shape_Changes/shchf_params_8.m new file mode 100644 index 0000000..7694c40 --- /dev/null +++ b/UserFiles/v4/rlhatton/Shape_Changes/shchf_params_8.m @@ -0,0 +1,82 @@ +function output = shchf_params_8(input_mode) + + % Default argument + if ~exist('input_mode','var') + + input_mode = 'initialize'; + + end + + % Name the .mat file with the fourier coefficients + paramfile = 'params_8'; + + + switch input_mode + + case 'name' + + output = 'params_8'; + + case 'dependency' + + output.dependency = {fullfile(pwd,[paramfile '.mat'])}; + + case 'initialize' + + %%%% + %% + %Path definitions + + % Load the points that the user clicked and the time vector + load(paramfile) + + % Check if start and end are the same point, and set boundary + % conditions accordingly + if alpha1(1)==alpha1(end) && alpha2(1)==alpha2(end) + splinemode = 'periodic'; + else + splinemode = 'complete'; + end + + % Generate spline structures for the selected points + spline_alpha1 = csape(t,alpha1',splinemode); + spline_alpha2 = csape(t,alpha2',splinemode); + + % The gait path is now defined by evaluating the two splines at + % specified times + p.phi_def = @(t) [fnval(spline_alpha1,t(:)),fnval(spline_alpha2,t(:))]; + + + % Speed up execution by defining the gait velocity as well (so + % the evaluator doesn't need to numerically take the derivative + % at every time + spline_dalpha1 = fnder(spline_alpha1); + spline_dalpha2 = fnder(spline_alpha2); + + % The gait path is now defined by evaluating the two splin + % derivatives at specified times + p.dphi_def = @(t) [fnval(spline_dalpha1,t(:)),fnval(spline_dalpha2,t(:))]; + + + %marker locations + p.phi_marker = []; + + %arrows to plot + p.phi_arrows = 2; + + %time to run path + p.time_def = t([1 end]); %#ok + + + %path resolution + p.phi_res = 50; + + + + + %%%% + %Save the shch properties + output = p; + end + +end diff --git a/UserFiles/v4/rlhatton/Shape_Changes/shchf_params_9.m b/UserFiles/v4/rlhatton/Shape_Changes/shchf_params_9.m new file mode 100644 index 0000000..50e25d2 --- /dev/null +++ b/UserFiles/v4/rlhatton/Shape_Changes/shchf_params_9.m @@ -0,0 +1,82 @@ +function output = shchf_params_9(input_mode) + + % Default argument + if ~exist('input_mode','var') + + input_mode = 'initialize'; + + end + + % Name the .mat file with the fourier coefficients + paramfile = 'params_9'; + + + switch input_mode + + case 'name' + + output = 'params_9'; + + case 'dependency' + + output.dependency = {fullfile(pwd,[paramfile '.mat'])}; + + case 'initialize' + + %%%% + %% + %Path definitions + + % Load the points that the user clicked and the time vector + load(paramfile) + + % Check if start and end are the same point, and set boundary + % conditions accordingly + if alpha1(1)==alpha1(end) && alpha2(1)==alpha2(end) + splinemode = 'periodic'; + else + splinemode = 'complete'; + end + + % Generate spline structures for the selected points + spline_alpha1 = csape(t,alpha1',splinemode); + spline_alpha2 = csape(t,alpha2',splinemode); + + % The gait path is now defined by evaluating the two splines at + % specified times + p.phi_def = @(t) [fnval(spline_alpha1,t(:)),fnval(spline_alpha2,t(:))]; + + + % Speed up execution by defining the gait velocity as well (so + % the evaluator doesn't need to numerically take the derivative + % at every time + spline_dalpha1 = fnder(spline_alpha1); + spline_dalpha2 = fnder(spline_alpha2); + + % The gait path is now defined by evaluating the two splin + % derivatives at specified times + p.dphi_def = @(t) [fnval(spline_dalpha1,t(:)),fnval(spline_dalpha2,t(:))]; + + + %marker locations + p.phi_marker = []; + + %arrows to plot + p.phi_arrows = 2; + + %time to run path + p.time_def = t([1 end]); %#ok + + + %path resolution + p.phi_res = 50; + + + + + %%%% + %Save the shch properties + output = p; + end + +end diff --git a/UserFiles/v4/rlhatton/Shape_Changes/shchf_params_edgegait.m b/UserFiles/v4/rlhatton/Shape_Changes/shchf_params_edgegait.m new file mode 100644 index 0000000..839008c --- /dev/null +++ b/UserFiles/v4/rlhatton/Shape_Changes/shchf_params_edgegait.m @@ -0,0 +1,98 @@ +function output = shchf_params_edgegait(input_mode) + + % Default argument + if ~exist('input_mode','var') + + input_mode = 'initialize'; + + end + + % Name the .mat file with the fourier coefficients + paramfile = 'params_edgegait'; + + + switch input_mode + + case 'name' + + output = 'params_edgegait'; + + case 'dependency' + + output.dependency = {fullfile(pwd,[paramfile '.mat'])}; + + case 'initialize' + + %%%% + %% + %Path definitions + + % Load the points that the user clicked and the time vector + load(paramfile) + + % Check if start and end are the same point, and set boundary + % conditions accordingly + if alpha1(1)==alpha1(end) && alpha2(1)==alpha2(end) + splinemode = 'periodic'; + else + splinemode = 'complete'; + end + + % Generate spline structures for the selected points + switch splinemode + + case 'periodic' + + % Fit a periodic spline to the selected points; the endslopes are found + % by averaging the positions of the points before and after the join + endslope1 = (alpha1(2)-alpha1(end-1))/(t(end)-t(end-2)); + endslope2 = (alpha2(2)-alpha2(end-1))/(t(end)-t(end-2)); + spline_alpha1 = spline(t,[endslope1;alpha1(:);endslope1]); + spline_alpha2 = spline(t,[endslope2;alpha2(:);endslope2]); + + case 'complete' + + % Fit a non-periodic spline to the selected points + spline_alpha1 = spline(t,alpha1(:)); + spline_alpha2 = spline(t,alpha2(:)); + + end + + % The gait path is now defined by evaluating the two splines at + % specified times + p.phi_def = @(t) [ppval(spline_alpha1,t(:)),ppval(spline_alpha2,t(:))]; + + + % Speed up execution by defining the gait velocity as well (so + % the evaluator doesn't need to numerically take the derivative + % at every time + spline_dalpha1 = ppdiff(spline_alpha1); + spline_dalpha2 = ppdiff(spline_alpha2); + + % The gait path is now defined by evaluating the two splin + % derivatives at specified times + p.dphi_def = @(t) [ppval(spline_dalpha1,t(:)),ppval(spline_dalpha2,t(:))]; + + + %marker locations + p.phi_marker = []; + + %arrows to plot + p.phi_arrows = 2; + + %time to run path + p.time_def = t([1 end]); %#ok + + + %path resolution + p.phi_res = 50; + + + + + %%%% + %Save the shch properties + output = p; + end + +end diff --git a/UserFiles/v4/rlhatton/Shape_Changes/shchf_params_mudskipper_curved.m b/UserFiles/v4/rlhatton/Shape_Changes/shchf_params_mudskipper_curved.m new file mode 100644 index 0000000..788ba6d --- /dev/null +++ b/UserFiles/v4/rlhatton/Shape_Changes/shchf_params_mudskipper_curved.m @@ -0,0 +1,82 @@ +function output = shchf_params_mudskipper_curved(input_mode) + + % Default argument + if ~exist('input_mode','var') + + input_mode = 'initialize'; + + end + + % Name the .mat file with the fourier coefficients + paramfile = 'params_mudskipper_curved'; + + + switch input_mode + + case 'name' + + output = 'params_mudskipper_curved'; + + case 'dependency' + + output.dependency = {fullfile(pwd,[paramfile '.mat'])}; + + case 'initialize' + + %%%% + %% + %Path definitions + + % Load the points that the user clicked and the time vector + load(paramfile) + + % Check if start and end are the same point, and set boundary + % conditions accordingly + if alpha1(1)==alpha1(end) && alpha2(1)==alpha2(end) + splinemode = 'periodic'; + else + splinemode = 'complete'; + end + + % Generate spline structures for the selected points + spline_alpha1 = csape(t,alpha1',splinemode); + spline_alpha2 = csape(t,alpha2',splinemode); + + % The gait path is now defined by evaluating the two splines at + % specified times + p.phi_def = @(t) [fnval(spline_alpha1,t(:)),fnval(spline_alpha2,t(:))]; + + + % Speed up execution by defining the gait velocity as well (so + % the evaluator doesn't need to numerically take the derivative + % at every time + spline_dalpha1 = fnder(spline_alpha1); + spline_dalpha2 = fnder(spline_alpha2); + + % The gait path is now defined by evaluating the two splin + % derivatives at specified times + p.dphi_def = @(t) [fnval(spline_dalpha1,t(:)),fnval(spline_dalpha2,t(:))]; + + + %marker locations + p.phi_marker = []; + + %arrows to plot + p.phi_arrows = 2; + + %time to run path + p.time_def = t([1 end]); %#ok + + + %path resolution + p.phi_res = 50; + + + + + %%%% + %Save the shch properties + output = p; + end + +end diff --git a/UserFiles/v4/rlhatton/Shape_Changes/shchf_params_mudskipper_direct.m b/UserFiles/v4/rlhatton/Shape_Changes/shchf_params_mudskipper_direct.m new file mode 100644 index 0000000..8e07ffa --- /dev/null +++ b/UserFiles/v4/rlhatton/Shape_Changes/shchf_params_mudskipper_direct.m @@ -0,0 +1,82 @@ +function output = shchf_params_mudskipper_direct(input_mode) + + % Default argument + if ~exist('input_mode','var') + + input_mode = 'initialize'; + + end + + % Name the .mat file with the fourier coefficients + paramfile = 'params_mudskipper_direct'; + + + switch input_mode + + case 'name' + + output = 'params_mudskipper_direct'; + + case 'dependency' + + output.dependency = {fullfile(pwd,[paramfile '.mat'])}; + + case 'initialize' + + %%%% + %% + %Path definitions + + % Load the points that the user clicked and the time vector + load(paramfile) + + % Check if start and end are the same point, and set boundary + % conditions accordingly + if alpha1(1)==alpha1(end) && alpha2(1)==alpha2(end) + splinemode = 'periodic'; + else + splinemode = 'complete'; + end + + % Generate spline structures for the selected points + spline_alpha1 = csape(t,alpha1',splinemode); + spline_alpha2 = csape(t,alpha2',splinemode); + + % The gait path is now defined by evaluating the two splines at + % specified times + p.phi_def = @(t) [fnval(spline_alpha1,t(:)),fnval(spline_alpha2,t(:))]; + + + % Speed up execution by defining the gait velocity as well (so + % the evaluator doesn't need to numerically take the derivative + % at every time + spline_dalpha1 = fnder(spline_alpha1); + spline_dalpha2 = fnder(spline_alpha2); + + % The gait path is now defined by evaluating the two splin + % derivatives at specified times + p.dphi_def = @(t) [fnval(spline_dalpha1,t(:)),fnval(spline_dalpha2,t(:))]; + + + %marker locations + p.phi_marker = []; + + %arrows to plot + p.phi_arrows = 2; + + %time to run path + p.time_def = t([1 end]); %#ok + + + %path resolution + p.phi_res = 50; + + + + + %%%% + %Save the shch properties + output = p; + end + +end diff --git a/UserFiles/v4/rlhatton/Shape_Changes/shchf_params_triangle_wavelength1_innercontour.m b/UserFiles/v4/rlhatton/Shape_Changes/shchf_params_triangle_wavelength1_innercontour.m new file mode 100644 index 0000000..53e9567 --- /dev/null +++ b/UserFiles/v4/rlhatton/Shape_Changes/shchf_params_triangle_wavelength1_innercontour.m @@ -0,0 +1,98 @@ +function output = shchf_params_triangle_wavelength1_innercontour(input_mode) + + % Default argument + if ~exist('input_mode','var') + + input_mode = 'initialize'; + + end + + % Name the .mat file with the fourier coefficients + paramfile = 'params_triangle_wavelength1_innercontour'; + + + switch input_mode + + case 'name' + + output = 'params_triangle_wavelength1_innercontour'; + + case 'dependency' + + output.dependency = {fullfile(pwd,[paramfile '.mat'])}; + + case 'initialize' + + %%%% + %% + %Path definitions + + % Load the points that the user clicked and the time vector + load(paramfile) + + % Check if start and end are the same point, and set boundary + % conditions accordingly + if alpha1(1)==alpha1(end) && alpha2(1)==alpha2(end) + splinemode = 'periodic'; + else + splinemode = 'complete'; + end + + % Generate spline structures for the selected points + switch splinemode + + case 'periodic' + + % Fit a periodic spline to the selected points; the endslopes are found + % by averaging the positions of the points before and after the join + endslope1 = (alpha1(2)-alpha1(end-1))/(t(end)-t(end-2)); + endslope2 = (alpha2(2)-alpha2(end-1))/(t(end)-t(end-2)); + spline_alpha1 = spline(t,[endslope1;alpha1(:);endslope1]); + spline_alpha2 = spline(t,[endslope2;alpha2(:);endslope2]); + + case 'complete' + + % Fit a non-periodic spline to the selected points + spline_alpha1 = spline(t,alpha1(:)); + spline_alpha2 = spline(t,alpha2(:)); + + end + + % The gait path is now defined by evaluating the two splines at + % specified times + p.phi_def = @(t) [ppval(spline_alpha1,t(:)),ppval(spline_alpha2,t(:))]; + + + % Speed up execution by defining the gait velocity as well (so + % the evaluator doesn't need to numerically take the derivative + % at every time + spline_dalpha1 = ppdiff(spline_alpha1); + spline_dalpha2 = ppdiff(spline_alpha2); + + % The gait path is now defined by evaluating the two splin + % derivatives at specified times + p.dphi_def = @(t) [ppval(spline_dalpha1,t(:)),ppval(spline_dalpha2,t(:))]; + + + %marker locations + p.phi_marker = []; + + %arrows to plot + p.phi_arrows = 2; + + %time to run path + p.time_def = t([1 end]); %#ok + + + %path resolution + p.phi_res = 50; + + + + + %%%% + %Save the shch properties + output = p; + end + +end diff --git a/UserFiles/v4/rlhatton/Shape_Changes/shchf_params_triangle_wavelength1_zerocontour.m b/UserFiles/v4/rlhatton/Shape_Changes/shchf_params_triangle_wavelength1_zerocontour.m new file mode 100644 index 0000000..380cd86 --- /dev/null +++ b/UserFiles/v4/rlhatton/Shape_Changes/shchf_params_triangle_wavelength1_zerocontour.m @@ -0,0 +1,98 @@ +function output = shchf_params_triangle_wavelength1_zerocontour(input_mode) + + % Default argument + if ~exist('input_mode','var') + + input_mode = 'initialize'; + + end + + % Name the .mat file with the fourier coefficients + paramfile = 'params_triangle_wavelength1_zerocontour'; + + + switch input_mode + + case 'name' + + output = 'params_triangle_wavelength1_zerocontour'; + + case 'dependency' + + output.dependency = {fullfile(pwd,[paramfile '.mat'])}; + + case 'initialize' + + %%%% + %% + %Path definitions + + % Load the points that the user clicked and the time vector + load(paramfile) + + % Check if start and end are the same point, and set boundary + % conditions accordingly + if alpha1(1)==alpha1(end) && alpha2(1)==alpha2(end) + splinemode = 'periodic'; + else + splinemode = 'complete'; + end + + % Generate spline structures for the selected points + switch splinemode + + case 'periodic' + + % Fit a periodic spline to the selected points; the endslopes are found + % by averaging the positions of the points before and after the join + endslope1 = (alpha1(2)-alpha1(end-1))/(t(end)-t(end-2)); + endslope2 = (alpha2(2)-alpha2(end-1))/(t(end)-t(end-2)); + spline_alpha1 = spline(t,[endslope1;alpha1(:);endslope1]); + spline_alpha2 = spline(t,[endslope2;alpha2(:);endslope2]); + + case 'complete' + + % Fit a non-periodic spline to the selected points + spline_alpha1 = spline(t,alpha1(:)); + spline_alpha2 = spline(t,alpha2(:)); + + end + + % The gait path is now defined by evaluating the two splines at + % specified times + p.phi_def = @(t) [ppval(spline_alpha1,t(:)),ppval(spline_alpha2,t(:))]; + + + % Speed up execution by defining the gait velocity as well (so + % the evaluator doesn't need to numerically take the derivative + % at every time + spline_dalpha1 = ppdiff(spline_alpha1); + spline_dalpha2 = ppdiff(spline_alpha2); + + % The gait path is now defined by evaluating the two splin + % derivatives at specified times + p.dphi_def = @(t) [ppval(spline_dalpha1,t(:)),ppval(spline_dalpha2,t(:))]; + + + %marker locations + p.phi_marker = []; + + %arrows to plot + p.phi_arrows = 2; + + %time to run path + p.time_def = t([1 end]); %#ok + + + %path resolution + p.phi_res = 50; + + + + + %%%% + %Save the shch properties + output = p; + end + +end diff --git a/UserFiles/v4/rlhatton/Shape_Changes/shchf_serpenoid_varfreq_test1.m b/UserFiles/v4/rlhatton/Shape_Changes/shchf_serpenoid_varfreq_test1.m new file mode 100644 index 0000000..3f0f220 --- /dev/null +++ b/UserFiles/v4/rlhatton/Shape_Changes/shchf_serpenoid_varfreq_test1.m @@ -0,0 +1,61 @@ +function output = shchf_serpenoid_varfreq_test1(input_mode) + + % Default argument + if ~exist('input_mode','var') + + input_mode = 'initialize'; + + end + + switch input_mode + + case 'name' + + output = 'Serpenoid varfreq test 1'; + + case 'dependency' + + output.dependency = {}; + + case 'initialize' + + %%%% + % Filename to save to + + %% + %Path definitions + + %path definition + p.phi_def{1} = {@strokedef}; + + %marker locations + p.phi_marker = []; + + %arrows to plot + p.phi_arrows{1} = {2}; + + %time to run path + p.time_def{1} = {[0 2]}; + + %p.cBVI_method{1}{1} = 'simple'; + + %path resolution + p.phi_res{1} = {100}; + + + %%%% + %Output the path properties + output = p; + end + +end + +function [stroke] = strokedef(t) + + t = t(:)'; + + stroke = [9*ones(size(t)); + 9*ones(size(t)); + .5 + 2*t].'; + +end \ No newline at end of file diff --git a/UserFiles/v4/rlhatton/Shape_Changes/shchf_serpenoid_varfreq_test2.m b/UserFiles/v4/rlhatton/Shape_Changes/shchf_serpenoid_varfreq_test2.m new file mode 100644 index 0000000..7df56d7 --- /dev/null +++ b/UserFiles/v4/rlhatton/Shape_Changes/shchf_serpenoid_varfreq_test2.m @@ -0,0 +1,61 @@ +function output = shchf_serpenoid_varfreq_test2(input_mode) + + % Default argument + if ~exist('input_mode','var') + + input_mode = 'initialize'; + + end + + switch input_mode + + case 'name' + + output = 'Serpenoid varfreq test 2'; + + case 'dependency' + + output.dependency = {}; + + case 'initialize' + + %%%% + % Filename to save to + + %% + %Path definitions + + %path definition + p.phi_def{1} = {@strokedef}; + + %marker locations + p.phi_marker = []; + + %arrows to plot + p.phi_arrows{1} = {2}; + + %time to run path + p.time_def{1} = {[0 1]}; + + %p.cBVI_method{1}{1} = 'simple'; + + %path resolution + p.phi_res{1} = {100}; + + + %%%% + %Output the path properties + output = p; + end + +end + +function [stroke] = strokedef(t) + + t = t(:)'; + + stroke = [zeros(size(t)); + 6*ones(size(t)); + .5 + 2*t].'; + +end \ No newline at end of file diff --git a/UserFiles/v4/rlhatton/Systems/Ross_vector_fields_00.mat b/UserFiles/v4/rlhatton/Systems/Ross_vector_fields_00.mat new file mode 100644 index 0000000..f7dda83 Binary files /dev/null and b/UserFiles/v4/rlhatton/Systems/Ross_vector_fields_00.mat differ diff --git a/UserFiles/v4/rlhatton/Systems/Ross_vector_fields_20.mat b/UserFiles/v4/rlhatton/Systems/Ross_vector_fields_20.mat new file mode 100644 index 0000000..d8b6b00 Binary files /dev/null and b/UserFiles/v4/rlhatton/Systems/Ross_vector_fields_20.mat differ diff --git a/UserFiles/v4/rlhatton/Systems/constrained_high_adduction_20.mat b/UserFiles/v4/rlhatton/Systems/constrained_high_adduction_20.mat new file mode 100644 index 0000000..8a79589 Binary files /dev/null and b/UserFiles/v4/rlhatton/Systems/constrained_high_adduction_20.mat differ diff --git a/UserFiles/v4/rlhatton/Systems/curv_serpenoid_def_compare.m b/UserFiles/v4/rlhatton/Systems/curv_serpenoid_def_compare.m new file mode 100644 index 0000000..54cd8ab --- /dev/null +++ b/UserFiles/v4/rlhatton/Systems/curv_serpenoid_def_compare.m @@ -0,0 +1,70 @@ +function output = curv_serpenoid_def_compare(params,mode) + +% Turn params into a cell matrix +params = num2cell(params); + +switch mode + + case 'curvature' + + output = @(s) curv_fun(s,params{:}); + + + case 'angle' + + output = @(s) int_curv_ds_fun(s,params{:}); + + + case 'dcurvature' + + output = @(s) d_curv_dp_fun(s,params{:}); + + + case 'dcurvature_int' + + output = @(s) int_d_curv_dp_ds_fun(s,params{:}); + + +end + +end + +function output = vector_to_list_input(funhandle,all_params) + + all_params = num2cell(all_params); + + output = funhandle(all_params{:}); + +end + + +function output = reshape_truncate_jacobian(J) + + output = J(2:end)'; + +end + +function out1 = curv_fun(s,a1,a2) + t2 = s.*pi.*2.0; + out1 = a1.*cos(t2)+a2.*sin(t2); +end + + +function out1 = int_curv_ds_fun(s,a1,a2) + t2 = sin(s.*pi); + out1 = (a1.*sin(s.*pi.*2.0).*(1.0./2.0)+a2.*t2.^2)./pi; +end + + +function out1 = d_curv_dp_fun(s,a1,a2) + t2 = s.*pi.*2.0; + out1 = [cos(t2);sin(t2)]; +end + + +function out1 = int_d_curv_dp_ds_fun(s,a1,a2) + t2 = 1.0./pi; + t3 = sin(s.*pi); + out1 = [t2.*sin(s.*pi.*2.0).*(1.0./2.0);t2.*t3.^2]; +end + diff --git a/UserFiles/v4/rlhatton/Systems/curv_serpenoid_def_test.m b/UserFiles/v4/rlhatton/Systems/curv_serpenoid_def_test.m new file mode 100644 index 0000000..e0d658f --- /dev/null +++ b/UserFiles/v4/rlhatton/Systems/curv_serpenoid_def_test.m @@ -0,0 +1,79 @@ +function output = curv_serpenoid_def_test(params,mode) + +% Turn params into a cell matrix +params = num2cell(params); + +switch mode + + case 'curvature' + + output = @(s) curv_fun(s,params{:}); + + + case 'angle' + + output = @(s) int_curv_ds_fun(s,params{:}); + + + case 'dcurvature' + + output = @(s) d_curv_dp_fun(s,params{:}); + + + case 'dcurvature_int' + + output = @(s) int_d_curv_dp_ds_fun(s,params{:}); + + +end + +end + +function output = vector_to_list_input(funhandle,all_params) + + all_params = num2cell(all_params); + + output = funhandle(all_params{:}); + +end + + +function output = reshape_truncate_jacobian(J) + + output = J(2:end)'; + +end + +function out1 = curv_fun(s,a1,a2,omega) + t2 = omega.*s.*pi.*2.0; + out1 = a1.*cos(t2)+a2.*sin(t2); +end + + +function out1 = int_curv_ds_fun(s,a1,a2,omega) + t2 = sin(omega.*s.*pi); + out1 = (a2.*t2.^2+a1.*sin(omega.*s.*pi.*2.0).*(1.0./2.0))./(omega.*pi); +end + + +function out1 = d_curv_dp_fun(s,a1,a2,omega) + t2 = omega.*s.*pi.*2.0; + t3 = cos(t2); + t4 = sin(t2); + out1 = [t3;t4;a1.*s.*t4.*pi.*-2.0+a2.*s.*t3.*pi.*2.0]; +end + + +function out1 = int_d_curv_dp_ds_fun(s,a1,a2,omega) + t2 = 1.0./omega; + t3 = 1.0./pi; + t7 = omega.*s.*pi; + t4 = sin(t7); + t5 = omega.*s.*pi.*2.0; + t6 = sin(t5); + t8 = t4.^2; + t9 = 1.0./omega.^2; + t10 = t8.*2.0; + out1 = [t2.*t3.*t6.*(1.0./2.0);t2.*t3.*t8;a2.*t3.*t9.*(t10-omega.*s.*t6.*pi.*2.0).*(-1.0./2.0)-a1.*t3.*t9.*(t6+omega.*s.*pi.*(t10-1.0).*2.0).*(1.0./2.0)]; +end + diff --git a/UserFiles/v4/rlhatton/Systems/curv_serpenoid_var_sharpness.m b/UserFiles/v4/rlhatton/Systems/curv_serpenoid_var_sharpness.m new file mode 100644 index 0000000..f588b65 --- /dev/null +++ b/UserFiles/v4/rlhatton/Systems/curv_serpenoid_var_sharpness.m @@ -0,0 +1,82 @@ +% abs(cos(atan2(a2, a1) + 2*s*pi))^n*sign(cos(atan2(a2, a1) + 2*s*pi))*(a1^2 + a2^2)^(1/2)*((153949471934815*n)/562949953421312 - (1370981160259359*n^2)/144115188075855872 + (174464881467265*n^3)/1152921504606846976 + 6368462077980939/9007199254740992) +function output = curv_serpenoid_var_sharpness(params,mode) + +% Turn params into a cell matrix +params = num2cell(params); + +switch mode + + case 'curvature' + + output = @(s) curv_fun(s,params{:}); + + + case 'angle' + + %% Padded length of unit backbone + all_limits = [-.51 0 .51]; + + %% Make dummy integration function + curv_fun_dummy = curv_serpenoid_var_sharpness(cell2mat(params),'curvature'); + curvature = @(s,~) curv_fun_dummy(s); + + %% Integral of the integrand function along s + output = ode_multistart(@ode45,curvature,all_limits,0,0); + + + case 'dcurvature' + + %% Create a dummy function that takes in a vector of parameters + %% including s, and deals them into individual function parameters + curv_intermediate = @(all_params) vector_to_list_input(@curv_fun,all_params); + + %% Create a function that takes the jacobian of the dummy function + fulljacobian = @(s) jacobianest(curv_intermediate,[s,params{:}]); + + %% Create a function that truncates the s-derivative from the full jacobian + output = @(s) reshape_truncate_jacobian(fulljacobian(s)); + + + case 'dcurvature_int' + + %% Padded length of unit backbone + all_limits = [-.51 0 .51]; + + %% Make dummy integration function + d_curv_dp_fun_dummy = curv_serpenoid_var_sharpness(cell2mat(params),'dcurvature'); + dcurvature = @(s,~) d_curv_dp_fun_dummy(s); + + %% Integral of the integrand function along s + output = ode_multistart(@ode45,dcurvature,all_limits,0,zeros(size(params(:).'))); + + +end + +end + +function output = vector_to_list_input(funhandle,all_params) + + all_params = num2cell(all_params); + + output = funhandle(all_params{:}); + +end + + +function output = reshape_truncate_jacobian(J) + + output = J(2:end)'; + +end + +function out1 = curv_fun(s,a1,a2,n) + t2 = a2.*1i; + t3 = a1+t2; + t4 = angle(t3); + t5 = s.*pi.*2.0; + t6 = t4+t5; + t7 = cos(t6); + t8 = n.^2; + out1 = abs(t7).^n.*sign(t7).*sqrt(a1.^2+a2.^2).*(n.*2.734691973934655e-1-t8.*9.513092815295325e-3+n.*t8.*1.513241628073878e-4+7.070413230426608e-1); +end + diff --git a/UserFiles/v4/rlhatton/Systems/curv_serpenoid_var_sharpness_symbolicd.m b/UserFiles/v4/rlhatton/Systems/curv_serpenoid_var_sharpness_symbolicd.m new file mode 100644 index 0000000..e68f0f1 --- /dev/null +++ b/UserFiles/v4/rlhatton/Systems/curv_serpenoid_var_sharpness_symbolicd.m @@ -0,0 +1,83 @@ +% cos(2*pi*s + atan2(a2, a1))^n*(a1^2 + a2^2) +function output = curv_serpenoid_var_sharpness_symbolicd(params,mode) + +% Turn params into a cell matrix +params = num2cell(params); + +switch mode + + case 'curvature' + + output = @(s) curv_fun(s,params{:}); + + + case 'angle' + + %% Padded length of unit backbone + all_limits = [-.51 0 .51]; + + %% Make dummy integration function + curv_fun_dummy = curv_serpenoid_var_sharpness_symbolicd(cell2mat(params),'curvature'); + curvature = @(s,~) curv_fun_dummy(s); + + %% Integral of the integrand function along s + output = ode_multistart(@ode45,curvature,all_limits,0,0); + + + case 'dcurvature' + + output = @(s) d_curv_dp_fun(s,params{:}); + + + case 'dcurvature_int' + + %% Padded length of unit backbone + all_limits = [-.51 0 .51]; + + %% Make dummy integration function + d_curv_dp_fun_dummy = curv_serpenoid_var_sharpness_symbolicd(cell2mat(params),'dcurvature'); + dcurvature = @(s,~) d_curv_dp_fun_dummy(s); + + %% Integral of the integrand function along s + output = ode_multistart(@ode45,dcurvature,all_limits,0,zeros(size(params(:).'))); + + +end + +end + +function output = vector_to_list_input(funhandle,all_params) + + all_params = num2cell(all_params); + + output = funhandle(all_params{:}); + +end + + +function output = reshape_truncate_jacobian(J) + + output = J(2:end)'; + +end + +function out1 = curv_fun(s,a1,a2,n) + out1 = cos(pi.*s.*2.0+atan2(a2,a1)).^n.*(a1.^2+a2.^2); +end + + +function out1 = d_curv_dp_fun(s,a1,a2,n) + t2 = pi.*s.*2.0; + t3 = atan2(a2,a1); + t4 = t2+t3; + t5 = cos(t4); + t6 = t5.^n; + t7 = n-1.0; + t8 = t5.^t7; + t9 = sin(t4); + t10 = a1.^2; + t11 = a2.^2; + t12 = t10+t11; + out1 = [a1.*t6.*2.0;a2.*t6.*2.0;t6.*t12.*log(t5)]; +end + diff --git a/UserFiles/v4/rlhatton/Systems/curv_serpenoid_varfreq.m b/UserFiles/v4/rlhatton/Systems/curv_serpenoid_varfreq.m new file mode 100644 index 0000000..80c7f57 --- /dev/null +++ b/UserFiles/v4/rlhatton/Systems/curv_serpenoid_varfreq.m @@ -0,0 +1,80 @@ +% (a1*cos(2*pi*s*omega)+a2*sin(2*pi*s*omega)) +function output = curv_serpenoid_varfreq(params,mode) + +% Turn params into a cell matrix +params = num2cell(params); + +switch mode + + case 'curvature' + + output = @(s) curv_fun(s,params{:}); + + + case 'angle' + + output = @(s) int_curv_ds_fun(s,params{:}); + + + case 'dcurvature' + + output = @(s) d_curv_dp_fun(s,params{:}); + + + case 'dcurvature_int' + + output = @(s) int_d_curv_dp_ds_fun(s,params{:}); + + +end + +end + +function output = vector_to_list_input(funhandle,all_params) + + all_params = num2cell(all_params); + + output = funhandle(all_params{:}); + +end + + +function output = reshape_truncate_jacobian(J) + + output = J(2:end)'; + +end + +function out1 = curv_fun(s,a1,a2,omega) + t2 = pi.*omega.*s.*2.0; + out1 = a1.*cos(t2)+a2.*sin(t2); +end + + +function out1 = int_curv_ds_fun(s,a1,a2,omega) + t2 = sin(pi.*omega.*s); + out1 = (a2.*t2.^2+a1.*sin(pi.*omega.*s.*2.0).*(1.0./2.0))./(pi.*omega); +end + + +function out1 = d_curv_dp_fun(s,a1,a2,omega) + t2 = pi.*omega.*s.*2.0; + t3 = cos(t2); + t4 = sin(t2); + out1 = [t3;t4;pi.*a1.*s.*t4.*-2.0+pi.*a2.*s.*t3.*2.0]; +end + + +function out1 = int_d_curv_dp_ds_fun(s,a1,a2,omega) + t2 = 1.0./pi; + t3 = 1.0./omega; + t7 = pi.*omega.*s; + t4 = sin(t7); + t5 = pi.*omega.*s.*2.0; + t6 = sin(t5); + t8 = t4.^2; + t9 = 1.0./omega.^2; + t10 = t8.*2.0; + out1 = [t2.*t3.*t6.*(1.0./2.0);t2.*t3.*t8;a2.*t2.*t9.*(t10-pi.*omega.*s.*t6.*2.0).*(-1.0./2.0)-a1.*t2.*t9.*(t6+pi.*omega.*s.*(t10-1.0).*2.0).*(1.0./2.0)]; +end + diff --git a/UserFiles/v4/rlhatton/Systems/curv_triangle_wave.m b/UserFiles/v4/rlhatton/Systems/curv_triangle_wave.m new file mode 100644 index 0000000..0624a55 --- /dev/null +++ b/UserFiles/v4/rlhatton/Systems/curv_triangle_wave.m @@ -0,0 +1,85 @@ +% (3665543810304250167*abs(cos(atan2(a2, a1) + 2*s*pi))^15*sign(cos(atan2(a2, a1) + 2*s*pi))*(a1^2 + a2^2)^(1/2))/1152921504606846976 +function output = curv_triangle_wave(params,mode) + +% Turn params into a cell matrix +params = num2cell(params); + +switch mode + + case 'curvature' + + output = @(s) curv_fun(s,params{:}); + + + case 'angle' + + %% Padded length of unit backbone + all_limits = [-.51 0 .51]; + + %% Make dummy integration function + curv_fun_dummy = curv_triangle_wave(cell2mat(params),'curvature'); + curvature = @(s,~) curv_fun_dummy(s); + + %% Integral of the integrand function along s + output = ode_multistart(@ode45,curvature,all_limits,0,0); + + + case 'dcurvature' + + %% Create a dummy function that takes in a vector of parameters + %% including s, and deals them into individual function parameters + curv_intermediate = @(all_params) vector_to_list_input(@curv_fun,all_params); + + %% Create a function that takes the jacobian of the dummy function + fulljacobian = @(s) jacobianest(curv_intermediate,[s,params{:}]); + + %% Create a function that truncates the s-derivative from the full jacobian + output = @(s) reshape_truncate_jacobian(fulljacobian(s)); + + + case 'dcurvature_int' + + %% Padded length of unit backbone + all_limits = [-.51 0 .51]; + + %% Make dummy integration function + d_curv_dp_fun_dummy = curv_triangle_wave(cell2mat(params),'dcurvature'); + dcurvature = @(s,~) d_curv_dp_fun_dummy(s); + + %% Integral of the integrand function along s + output = ode_multistart(@ode45,dcurvature,all_limits,0,zeros(size(params(:).'))); + + +end + +end + +function output = vector_to_list_input(funhandle,all_params) + + all_params = num2cell(all_params); + + output = funhandle(all_params{:}); + +end + + +function output = reshape_truncate_jacobian(J) + + output = J(2:end)'; + +end + +function out1 = curv_fun(s,a1,a2) + t6 = a2.*1i; + t7 = a1+t6; + t8 = angle(t7); + t9 = s.*pi.*2.0; + t10 = t8+t9; + t11 = cos(t10); + t2 = abs(t11); + t3 = t2.^2; + t4 = t3.^2; + t5 = t4.^2; + out1 = t2.*t3.*t4.*t5.*sign(t11).*sqrt(a1.^2+a2.^2).*3.179352449978129; +end + diff --git a/UserFiles/v4/rlhatton/Systems/floatsnake_connection.m b/UserFiles/v4/rlhatton/Systems/floatsnake_connection.m new file mode 100644 index 0000000..27998cc --- /dev/null +++ b/UserFiles/v4/rlhatton/Systems/floatsnake_connection.m @@ -0,0 +1,35 @@ +function A = floatsnake_connection(a1,a2) +%FLOATSNAKE_CONNECTION +% A = FLOATSNAKE_CONNECTION(A1,A2) + +% This function was generated by the Symbolic Math Toolbox version 6.3. +% 31-Aug-2016 15:49:32 + +t2 = a1+a2; +t3 = a1.*2.0; +t4 = cos(t2); +t5 = t4.*2.0; +t6 = cos(a1); +t7 = t6.*6.0; +t8 = cos(a2); +t9 = t8.*6.0; +t10 = t5+t7+t9+1.3e1; +t11 = 1.0./t10; +t12 = a2.*2.0; +t13 = sin(t2); +t14 = t13.*3.0; +t15 = sin(a1); +t16 = sin(a2); +t17 = a1-a2; +t18 = sin(t17); +t19 = t18.*9.0; +t20 = a1+t12; +t21 = sin(t20); +t22 = a2+t3; +t23 = sin(t22); +t24 = t4.*9.0; +t25 = cos(t17); +t26 = t25.*9.0; +t27 = cos(t20); +t28 = cos(t22); +A = reshape([t11.*(t14+t15.*2.1e1-t16.*7.0+t19-t21+t23+sin(t3).*3.0).*(1.0./1.2e1),t11.*(t6.*2.1e1+t8.*7.0+t24+t26+t27+t28+cos(t3).*3.0+3.0).*(1.0./1.2e1),-t11.*(t4+t6.*3.0+3.0),t11.*(t14-t15.*7.0+t16.*2.1e1-t19+t21-t23+sin(t12).*3.0).*(-1.0./1.2e1),t11.*(t6.*7.0+t8.*2.1e1+t24+t26+t27+t28+cos(t12).*3.0+3.0).*(1.0./1.2e1),t11.*(t4+t8.*3.0+3.0)],[3,2]); diff --git a/UserFiles/v4/rlhatton/Systems/floatsnake_connection_elie.m b/UserFiles/v4/rlhatton/Systems/floatsnake_connection_elie.m new file mode 100644 index 0000000..b5c3e5d --- /dev/null +++ b/UserFiles/v4/rlhatton/Systems/floatsnake_connection_elie.m @@ -0,0 +1,35 @@ +function A_elie = floatsnake_connection_elie(a1,a2) +%FLOATSNAKE_CONNECTION_ELIE +% A_ELIE = FLOATSNAKE_CONNECTION_ELIE(A1,A2) + +% This function was generated by the Symbolic Math Toolbox version 6.3. +% 31-Aug-2016 18:44:58 + +t2 = a1+a2; +t3 = a1.*2.0; +t4 = cos(t2); +t5 = t4.*2.0; +t6 = cos(a1); +t7 = t6.*6.0; +t8 = cos(a2); +t9 = t8.*6.0; +t10 = t5+t7+t9+1.9e1; +t11 = 1.0./t10; +t12 = a2.*2.0; +t13 = sin(t2); +t14 = t13.*3.0; +t15 = sin(a1); +t16 = sin(a2); +t17 = a1-a2; +t18 = sin(t17); +t19 = t18.*9.0; +t20 = a1+t12; +t21 = sin(t20); +t22 = a2+t3; +t23 = sin(t22); +t24 = t4.*9.0; +t25 = cos(t17); +t26 = t25.*9.0; +t27 = cos(t20); +t28 = cos(t22); +A_elie = reshape([t11.*(t14+t15.*2.9e1-t16.*1.1e1+t19-t21+t23+sin(t3).*3.0).*(1.0./6.0),t11.*(t6.*2.9e1+t8.*1.1e1+t24+t26+t27+t28+cos(t3).*3.0+3.0).*(1.0./6.0),-t11.*(t4+t6.*3.0+5.0),t11.*(t14-t15.*1.1e1+t16.*2.9e1-t19+t21-t23+sin(t12).*3.0).*(-1.0./6.0),t11.*(t6.*1.1e1+t8.*2.9e1+t24+t26+t27+t28+cos(t12).*3.0+3.0).*(1.0./6.0),t11.*(t4+t8.*3.0+5.0)],[3,2]); diff --git a/UserFiles/v4/rlhatton/Systems/sysf_float_snake.m b/UserFiles/v4/rlhatton/Systems/sysf_float_snake.m index 7a5a67e..ddecb37 100644 --- a/UserFiles/v4/rlhatton/Systems/sysf_float_snake.m +++ b/UserFiles/v4/rlhatton/Systems/sysf_float_snake.m @@ -22,8 +22,8 @@ case 'dependency' - output.dependency = {'honey_power_metric.m'}; - output.product = {[datapath mfilename '.mat']}; % the .mat file sharing the same name as this .m file + output.dependency = {}; +% output.product = {[datapath mfilename '.mat']}; % the .mat file sharing the same name as this .m file case 'initialize' @@ -46,7 +46,7 @@ s.singularity = 0; %Range over which to evaluate connection - s.grid_range = [-1,1,-1,1]*2.5; + s.grid_range = [-1,1,-1,1]*3; %densities for various operations s.density.vector = [11 11]; %density to display vector field @@ -59,8 +59,8 @@ %Display parameters %shape space tic locations - s.tic_locs.x = [-1 0 1]*1.5; - s.tic_locs.y = [-1 0 1]*1.5; + s.tic_locs.x = [-1 0 1]*2.5; + s.tic_locs.y = [-1 0 1]*2.5; %%%% @@ -77,11 +77,10 @@ function A_num = Conn_num(a1,a2) - R=1/6; - A_num = R * [zeros(size(a1)) zeros(size(a1)); + A_num =[zeros(size(a1)) zeros(size(a1)); zeros(size(a1)) zeros(size(a1)); - -(5+3*cos(a1)+cos(a1-a2)) (5+3*cos(a2)+cos(a1-a2))]; + -(5+3*cos(a1)+cos(a1+a2)) (5+3*cos(a2)+cos(a1+a2))]; @@ -93,6 +92,6 @@ A_den = [ ones(size(a1)), ones(size(a1)); ones(size(a1)), ones(size(a1)) ; - (19+6*cos(a1)+6*cos(a2)+2*cos(a1-a2)) (19+6*cos(a1)+6*cos(a2)+2*cos(a1-a2))]; + (19+6*cos(a1)+6*cos(a2)+2*cos(a1+a2)) (19+6*cos(a1)+6*cos(a2)+2*cos(a1+a2))]; end \ No newline at end of file diff --git a/UserFiles/v4/rlhatton/Systems/sysf_float_snake.mat b/UserFiles/v4/rlhatton/Systems/sysf_float_snake.mat index d938bdd..20428d6 100644 Binary files a/UserFiles/v4/rlhatton/Systems/sysf_float_snake.mat and b/UserFiles/v4/rlhatton/Systems/sysf_float_snake.mat differ diff --git a/UserFiles/v4/rlhatton/Systems/sysf_float_snake_generated.m b/UserFiles/v4/rlhatton/Systems/sysf_float_snake_generated.m new file mode 100644 index 0000000..995c71b --- /dev/null +++ b/UserFiles/v4/rlhatton/Systems/sysf_float_snake_generated.m @@ -0,0 +1,112 @@ +function output = sysf_float_snake_generated(input_mode,datapath) + + % Default arguments + if ~exist('input_mode','var') + + input_mode = 'name'; + + end + + if ~exist('datapath','var') + + datapath = ''; + end + + %%%%%%% + + switch input_mode + + case 'name' + + output = 'Floatsnake generated'; % Display name + + case 'dependency' + + output.dependency = {'floatsnake_connection'}; + + case 'initialize' + + %Initialize a kinematic snake with unit values for L and R + + %%%%% + % Filename to save to + output = mfilename; + + + %Functional representation of local connection + %Functional representation of local connection + s.A_num = @Conn_num; + s.A_den = @Conn_den; + + %%% + %Processing details + + %Mark that system does not have a singularity that needs to be accounted for + s.singularity = 0; + + %Range over which to evaluate connection + s.grid_range = [-1,1,-1,1]*3; + + %densities for various operations + s.density.vector = [11 11]; %density to display vector field + s.density.scalar = [31 31]; %density to display scalar functions + s.density.eval = [31 31]; %density for function evaluations + s.finite_element_density = 31; + + + %%% + %Display parameters + + %shape space tic locations + s.tic_locs.x = [-1 0 1]*1.5; + s.tic_locs.y = [-1 0 1]*1.5; + + + %%%% + %Save the system properties + save(fullfile(datapath, mfilename),'s') + + %%%% + %Output the system properties + output = s; + + end + +end + +function [Ar]=Conn_num(a1,a2) + + + Ar_cell = cell(size(a1)); + for i = 1:numel(a1); + Ar_cell{i} = A_num_helper(a1(i),a2(i)); + waitbar2a(i/numel(a1)); + end + + + %Pull Ar values into a matrix + Ar_woven = cell2mat(Ar_cell); + + %Rearrange the A matrix + Ar = [ Ar_woven(1:3:end,1:2:end) Ar_woven(1:3:end,2:2:end); + Ar_woven(2:3:end,1:2:end) Ar_woven(2:3:end,2:2:end); + Ar_woven(3:3:end,1:2:end) Ar_woven(3:3:end,2:2:end)]; + +end + +function A = A_num_helper(a1,a2) + +% % Add the path to the curvature functions +% addpath(genpath('/Users/rlhatton/Documents/MATLAB/curvature_mode_toolbox')) +% addpath(genpath('/Users/rlhatton/Documents/MATLAB/LowRE_toolbox/')) + % Get the local connection for the specified shape, with unit length + A = floatsnake_connection(a1,a2); + + +end + +function A_den = Conn_den(a1,a2) %#ok + + A_den = repmat(ones(size(a1)),[3,2]); + +end \ No newline at end of file diff --git a/UserFiles/v4/rlhatton/Systems/sysf_float_snake_generated.mat b/UserFiles/v4/rlhatton/Systems/sysf_float_snake_generated.mat new file mode 100644 index 0000000..6220007 Binary files /dev/null and b/UserFiles/v4/rlhatton/Systems/sysf_float_snake_generated.mat differ diff --git a/UserFiles/v4/rlhatton/Systems/sysf_float_snake_generated_elie.m b/UserFiles/v4/rlhatton/Systems/sysf_float_snake_generated_elie.m new file mode 100644 index 0000000..8ace52a --- /dev/null +++ b/UserFiles/v4/rlhatton/Systems/sysf_float_snake_generated_elie.m @@ -0,0 +1,112 @@ +function output = sysf_float_snake_generated_elie(input_mode,datapath) + + % Default arguments + if ~exist('input_mode','var') + + input_mode = 'name'; + + end + + if ~exist('datapath','var') + + datapath = ''; + end + + %%%%%%% + + switch input_mode + + case 'name' + + output = 'Floatsnake generated Elie'; % Display name + + case 'dependency' + + output.dependency = {'floatsnake_connection'}; + + case 'initialize' + + %Initialize a kinematic snake with unit values for L and R + + %%%%% + % Filename to save to + output = mfilename; + + + %Functional representation of local connection + %Functional representation of local connection + s.A_num = @Conn_num; + s.A_den = @Conn_den; + + %%% + %Processing details + + %Mark that system does not have a singularity that needs to be accounted for + s.singularity = 0; + + %Range over which to evaluate connection + s.grid_range = [-1,1,-1,1]*3; + + %densities for various operations + s.density.vector = [11 11]; %density to display vector field + s.density.scalar = [31 31]; %density to display scalar functions + s.density.eval = [31 31]; %density for function evaluations + s.finite_element_density = 31; + + + %%% + %Display parameters + + %shape space tic locations + s.tic_locs.x = [-1 0 1]*1.5; + s.tic_locs.y = [-1 0 1]*1.5; + + + %%%% + %Save the system properties + save(fullfile(datapath, mfilename),'s') + + %%%% + %Output the system properties + output = s; + + end + +end + +function [Ar]=Conn_num(a1,a2) + + + Ar_cell = cell(size(a1)); + for i = 1:numel(a1); + Ar_cell{i} = A_num_helper(a1(i),a2(i)); + waitbar2a(i/numel(a1)); + end + + + %Pull Ar values into a matrix + Ar_woven = cell2mat(Ar_cell); + + %Rearrange the A matrix + Ar = [ Ar_woven(1:3:end,1:2:end) Ar_woven(1:3:end,2:2:end); + Ar_woven(2:3:end,1:2:end) Ar_woven(2:3:end,2:2:end); + Ar_woven(3:3:end,1:2:end) Ar_woven(3:3:end,2:2:end)]; + +end + +function A = A_num_helper(a1,a2) + +% % Add the path to the curvature functions +% addpath(genpath('/Users/rlhatton/Documents/MATLAB/curvature_mode_toolbox')) +% addpath(genpath('/Users/rlhatton/Documents/MATLAB/LowRE_toolbox/')) + % Get the local connection for the specified shape, with unit length + A = floatsnake_connection_elie(a1,a2); + + +end + +function A_den = Conn_den(a1,a2) %#ok + + A_den = repmat(ones(size(a1)),[3,2]); + +end \ No newline at end of file diff --git a/UserFiles/v4/rlhatton/Systems/sysf_float_snake_generated_elie.mat b/UserFiles/v4/rlhatton/Systems/sysf_float_snake_generated_elie.mat new file mode 100644 index 0000000..2c37a10 Binary files /dev/null and b/UserFiles/v4/rlhatton/Systems/sysf_float_snake_generated_elie.mat differ diff --git a/UserFiles/v4/rlhatton/Systems/sysf_honey_swimmer_piecewise_2modes.m b/UserFiles/v4/rlhatton/Systems/sysf_honey_swimmer_piecewise_2modes.m new file mode 100644 index 0000000..f0e18dc --- /dev/null +++ b/UserFiles/v4/rlhatton/Systems/sysf_honey_swimmer_piecewise_2modes.m @@ -0,0 +1,120 @@ +function output = sysf_honey_swimmer_piecewise_2modes(input_mode) + + % Default arguments + if ~exist('input_mode','var') + + input_mode = 'initialize'; + + end + + if ~exist('datapath','var') + + datapath = ''; + end + + %%%%%%% + + switch input_mode + + case 'name' + + output = 'Low Re piecewise continuous two modes'; % Display name + + case 'dependency' + + output.dependency = {'Utilities/curvature_mode_toolbox/backbone_from_curvature_bases.m', + 'Utilities/curvature_mode_toolbox/curvatures/constant_curvature_1.m', + 'Utilities/curvature_mode_toolbox/curvatures/constant_curvature_2.m', + 'Utilities/LowRE_toolbox/LowRE_dissipation_metric_from_curvature_bases.m', + 'Utilities/LowRE_toolbox/LowRE_local_connection_from_curvature_bases.m'}; + + case 'initialize' + + %Initialize a kinematic snake with unit values for L and R + + %%%%% + % Filename to save to + output = mfilename; + + + %Functional representation of local connection + s.A_num = @Conn_num; + s.A_den = @Conn_den; + + + %%% + %Processing details + + %Mark that system has a singularity that needs to be accounted for + s.singularity = 0; + + %Range over which to evaluate connection + s.grid_range = [-1,1,-1,1]*12; + + %densities for various operations + s.density.vector = [11 11]; %density to display vector field + s.density.scalar = [21 21]; %density to display scalar functions + s.density.eval = [21 21]; %density for function evaluations + s.finite_element_density = 11; + % power metric + s.metric = @(x,y) LowRE_dissipation_metric_from_curvature_bases... + ({@(s) constant_curvature_1(s);@(s) constant_curvature_2(s)},[x;y],1,1,2); % @(x,y) eye(2);% + + %%% + %Display parameters + + %shape space tic locations + s.tic_locs.x = [-1 0 1]*6; + s.tic_locs.y = [-1 0 1]*6; + + + %%%% + %Save the system properties + output = s; + + + end + +end + +function [Ar]=Conn_num(a1,a2) + + + %Apply the inverse multiplication + wb = waitbar2a(0,['Building ' num2str(size(a1)) ' connection matrix']); + Ar_cell = cell(size(a1)); + for i = 1:numel(a1); + Ar_cell{i} = A_num_helper(a1(i),a2(i)); + waitbar2a(i/numel(a1)); + end + + close(wb) + + %Ar_cell = arrayfun(@(a1,a2) A_num_helper(a1,a2),a1,a2,'UniformOutput',false); + + %Pull Ar values into a matrix + Ar_woven = cell2mat(Ar_cell); + + %Rearrange the A matrix + Ar = [ Ar_woven(1:3:end,1:2:end) Ar_woven(1:3:end,2:2:end); + Ar_woven(2:3:end,1:2:end) Ar_woven(2:3:end,2:2:end); + Ar_woven(3:3:end,1:2:end) Ar_woven(3:3:end,2:2:end)]; + +end + +function A = A_num_helper(a1,a2) + +% % Add the path to the curvature functions +% addpath(genpath('/Users/rlhatton/Documents/MATLAB/curvature_mode_toolbox')) +% addpath(genpath('/Users/rlhatton/Documents/MATLAB/LowRE_toolbox/')) + % Get the local connection for the specified shape, with unit length + A = LowRE_local_connection_from_curvature_bases({@(s) constant_curvature_1(s);@(s) constant_curvature_2(s)},[a1;a2],1,1,2); + + +end + +function A_den = Conn_den(a1,a2) %#ok + + A_den = repmat(ones(size(a1)),[3,2]); + +end \ No newline at end of file diff --git a/UserFiles/v4/rlhatton/Systems/sysf_honey_swimmer_serpenoid.m b/UserFiles/v4/rlhatton/Systems/sysf_honey_swimmer_serpenoid.m index 8b7190b..e55a732 100644 --- a/UserFiles/v4/rlhatton/Systems/sysf_honey_swimmer_serpenoid.m +++ b/UserFiles/v4/rlhatton/Systems/sysf_honey_swimmer_serpenoid.m @@ -83,7 +83,7 @@ %Apply the inverse multiplication wb = waitbar2a(0,['Building ' num2str(size(a1)) ' connection matrix']); Ar_cell = cell(size(a1)); - for i = 1:numel(a1); + parfor i = 1:numel(a1); Ar_cell{i} = A_num_helper(a1(i),a2(i)); waitbar2a(i/numel(a1)); end diff --git a/UserFiles/v4/rlhatton/Systems/sysf_honey_swimmer_serpenoid_general.m b/UserFiles/v4/rlhatton/Systems/sysf_honey_swimmer_serpenoid_general.m new file mode 100644 index 0000000..8be0d65 --- /dev/null +++ b/UserFiles/v4/rlhatton/Systems/sysf_honey_swimmer_serpenoid_general.m @@ -0,0 +1,128 @@ +function output = sysf_honey_swimmer_serpenoid_general(input_mode) + + % Default arguments + if ~exist('input_mode','var') + + input_mode = 'initialize'; + + end + + if ~exist('datapath','var') + + datapath = ''; + end + + %%%%%%% + + switch input_mode + + case 'name' + + output = 'Low Re Serpenoid first modes generalized curvature'; % Display name + + case 'dependency' + + output.dependency = {'Utilities/curvature_mode_toolbox/backbone_from_general_curvature.m'; + 'Utilities/curvature_mode_toolbox/curvatures/curv_serpenoid_unit'; + 'Utilities/LowRE_toolbox/LowRE_dissipation_metric_from_general_curvature.m'; + 'Utilities/LowRE_toolbox/LowRE_local_connection_from_general_curvature.m'}; + + case 'initialize' + + %Initialize a kinematic snake with unit values for L and R + + %%%%% + % Filename to save to + output = mfilename; + + + %Functional representation of local connection + s.A_num = @Conn_num; + s.A_den = @Conn_den; + + + %%% + %Processing details + + %Mark that system has a singularity that needs to be accounted for + s.singularity = 0; + + %Range over which to evaluate connection + s.grid_range = [-1,1,-1,1]*12; + + %densities for various operations + s.density.vector = [11 11]; %density to display vector field + s.density.scalar = [21 21]; %density to display scalar functions + s.density.eval = [21 21]; %density for function evaluations + s.finite_element_density = 11; + % power metric + s.metric = @(x,y) LowRE_dissipation_metric_from_general_curvature... + (@curvelock,[x;y],1,1,2); % @(x,y) eye(2);% + + %%% + %Display parameters + + %shape space tic locations + s.tic_locs.x = [-1 0 1]*6; + s.tic_locs.y = [-1 0 1]*6; + + + %%%% + %Save the system properties + output = s; + + + end + +end + +function [Ar]=Conn_num(a1,a2) + + + %Apply the inverse multiplication + wb = waitbar2a(0,['Building ' num2str(size(a1)) ' connection matrix']); + Ar_cell = cell(size(a1)); + parfor i = 1:numel(a1); + Ar_cell{i} = A_num_helper(a1(i),a2(i)); + waitbar2a(i/numel(a1)); + end + + close(wb) + + %Ar_cell = arrayfun(@(a1,a2) A_num_helper(a1,a2),a1,a2,'UniformOutput',false); + + %Pull Ar values into a matrix + Ar_woven = cell2mat(Ar_cell); + + %Rearrange the A matrix + Ar = [ Ar_woven(1:3:end,1:2:end) Ar_woven(1:3:end,2:2:end); + Ar_woven(2:3:end,1:2:end) Ar_woven(2:3:end,2:2:end); + Ar_woven(3:3:end,1:2:end) Ar_woven(3:3:end,2:2:end)]; + +end + +function A = A_num_helper(a1,a2) + +% % Add the path to the curvature functions +% addpath(genpath('/Users/rlhatton/Documents/MATLAB/curvature_mode_toolbox')) +% addpath(genpath('/Users/rlhatton/Documents/MATLAB/LowRE_toolbox/')) + % Get the local connection for the specified shape, with unit length + A = LowRE_local_connection_from_general_curvature(@curvelock,[a1;a2],1,1,2); + + +end + +function A_den = Conn_den(a1,a2) %#ok + + A_den = repmat(ones(size(a1)),[3,2]); + +end + + +function output = curvelock(s,params,mode) + +% lock the variable frequency serpenoid curvature to frequency 1 + + output = curv_serpenoid_unit(s,params,mode); + +end \ No newline at end of file diff --git a/UserFiles/v4/rlhatton/Systems/sysf_honey_swimmer_serpenoid_var_sharpness.m b/UserFiles/v4/rlhatton/Systems/sysf_honey_swimmer_serpenoid_var_sharpness.m new file mode 100644 index 0000000..b7395b3 --- /dev/null +++ b/UserFiles/v4/rlhatton/Systems/sysf_honey_swimmer_serpenoid_var_sharpness.m @@ -0,0 +1,133 @@ +function output = sysf_honey_swimmer_serpenoid_var_sharpness(input_mode) + +% Name for the curvature functions +curv_name = 'serpenoid_var_sharpness'; +curv_fun = str2func(['curv_' curv_name]); + + % Default arguments + if ~exist('input_mode','var') + + input_mode = 'initialize'; + + end + + if ~exist('datapath','var') + + datapath = ''; + end + + %%%%%%% + + switch input_mode + + case 'name' + + output = 'Low Re Serpenoid first modes variable sharpness'; % Display name + + case 'dependency' + + + + output.dependency = {'Utilities/curvature_mode_toolbox/backbone_from_general_curvature.m'; + ['Utilities/curvature_mode_toolbox/curvatures/curv_' curv_name '.m']; + 'Utilities/LowRE_toolbox/LowRE_dissipation_metric_from_general_curvature.m'; + 'Utilities/LowRE_toolbox/LowRE_local_connection_from_general_curvature.m'}; + + case 'initialize' + + + %%%%% + % Filename to save to + output = mfilename; + + + %Functional representation of local connection + s.A_num = @(a1,a2,a3) Conn_num(a1,a2,a3,curv_fun); + s.A_den = @Conn_den; + + + %%% + %Processing details + + %Mark that system has a singularity that needs to be accounted for + s.singularity = 0; + + %Range over which to evaluate connection + s.grid_range = [[-1,1,-1,1]*10,[.1 20]]; + + %densities for various operations + s.density.vector = [1,1,1]*11; %density to display vector field + s.density.scalar = [1,1,1]*11; %density to display scalar functions + s.density.eval = [1,1,1]*21;%[21 21 21]; %density for function evaluations + s.finite_element_density = 11;%11; + % power metric + s.metric = @(x,y,z) LowRE_dissipation_metric_from_general_curvature... + (curv_fun,[x;y;z],1,1,2); % @(x,y) eye(2);% + + %%% + %Display parameters + + %shape space tic locations + s.tic_locs.x = [-1 0 1]*6; + s.tic_locs.y = [-1 0 1]*6; + + + %%%% + %Save the system properties + output = s; + + + end + +end + +function [Ar]=Conn_num(a1,a2,a3,curv_fun) + + + for i = 1:nargin + sizea1(i) = size(a1,i); + end + + + + %Apply the inverse multiplication + wb = waitbar2a(0,['Building ' num2str(size(a1)) ' connection matrix']); + Ar_cell = cell(sizea1(1),sizea1(2),sizea1(3)); + parfor i = 1:numel(a1); + Ar_cell{i} = A_num_helper(a1(i),a2(i),a3(i),curv_fun); + waitbar2a(i/numel(a1)); + end + + close(wb) + + Ar1=cell(3); + for i=1:1:9 + for j=1:sizea1(3) + for k=1:sizea1(2) + for l=1:sizea1(1) + Ar1{i}(l,k,j)=Ar_cell{l,k,j}(i); + end + end + end + end + + Ar=cell2mat(Ar1); + +end + +function A = A_num_helper(a1,a2,a3,curv_fun) + +% % Add the path to the curvature functions +% addpath(genpath('/Users/rlhatton/Documents/MATLAB/curvature_mode_toolbox')) +% addpath(genpath('/Users/rlhatton/Documents/MATLAB/LowRE_toolbox/')) + % Get the local connection for the specified shape, with unit length + A = LowRE_local_connection_from_general_curvature(curv_fun,[a1;a2;a3],1,1,2); + + +end + +function A_den = Conn_den(a1,a2,a3) %#ok + + A_den = repmat(ones(size(a1)),[3,3]); + +end \ No newline at end of file diff --git a/UserFiles/v4/rlhatton/Systems/sysf_honey_swimmer_serpenoid_varfreq.m b/UserFiles/v4/rlhatton/Systems/sysf_honey_swimmer_serpenoid_varfreq.m new file mode 100644 index 0000000..7e45a05 --- /dev/null +++ b/UserFiles/v4/rlhatton/Systems/sysf_honey_swimmer_serpenoid_varfreq.m @@ -0,0 +1,128 @@ +function output = sysf_honey_swimmer_serpenoid_varfreq(input_mode) + + % Default arguments + if ~exist('input_mode','var') + + input_mode = 'initialize'; + + end + + if ~exist('datapath','var') + + datapath = ''; + end + + %%%%%%% + + switch input_mode + + case 'name' + + output = 'Low Re Serpenoid first modes variable frequency'; % Display name + + case 'dependency' + + output.dependency = {'Utilities/curvature_mode_toolbox/backbone_from_general_curvature.m'; + 'Utilities/curvature_mode_toolbox/curvatures/curv_serpenoid_varfreq'; + 'Utilities/LowRE_toolbox/LowRE_dissipation_metric_from_general_curvature.m'; + 'Utilities/LowRE_toolbox/LowRE_local_connection_from_general_curvature.m'}; + + case 'initialize' + + %Initialize a kinematic snake with unit values for L and R + + %%%%% + % Filename to save to + output = mfilename; + + + %Functional representation of local connection + s.A_num = @Conn_num; + s.A_den = @Conn_den; + + + %%% + %Processing details + + %Mark that system has a singularity that needs to be accounted for + s.singularity = 0; + + %Range over which to evaluate connection + s.grid_range = [[-1,1,-1,1]*10,[.1 5]]; + + %densities for various operations + s.density.vector = [1,1,1]*11; %density to display vector field + s.density.scalar = [1,1,1]*11; %density to display scalar functions + s.density.eval = [1,1,1]*21;%[21 21 21]; %density for function evaluations + s.finite_element_density = 11;%11; + % power metric + s.metric = @(x,y,z) LowRE_dissipation_metric_from_general_curvature... + (@curv_serpenoid_def_test,[x;y;z],1,1,2); % @(x,y) eye(2);% + + %%% + %Display parameters + + %shape space tic locations + s.tic_locs.x = [-1 0 1]*6; + s.tic_locs.y = [-1 0 1]*6; + + + %%%% + %Save the system properties + output = s; + + + end + +end + +function [Ar]=Conn_num(a1,a2,a3) + + + for i = 1:nargin + sizea1(i) = size(a1,i); + end + + + + %Apply the inverse multiplication + wb = waitbar2a(0,['Building ' num2str(size(a1)) ' connection matrix']); + Ar_cell = cell(sizea1(1),sizea1(2),sizea1(3)); + parfor i = 1:numel(a1); + Ar_cell{i} = A_num_helper(a1(i),a2(i),a3(i)); + waitbar2a(i/numel(a1)); + end + + close(wb) + + Ar1=cell(3); + for i=1:1:9 + for j=1:sizea1(3) + for k=1:sizea1(2) + for l=1:sizea1(1) + Ar1{i}(l,k,j)=Ar_cell{l,k,j}(i); + end + end + end + end + + Ar=cell2mat(Ar1); + +end + +function A = A_num_helper(a1,a2,a3) + +% % Add the path to the curvature functions +% addpath(genpath('/Users/rlhatton/Documents/MATLAB/curvature_mode_toolbox')) +% addpath(genpath('/Users/rlhatton/Documents/MATLAB/LowRE_toolbox/')) + % Get the local connection for the specified shape, with unit length + A = LowRE_local_connection_from_general_curvature(@curv_serpenoid_def_test,[a1;a2;a3],1,1,2); + + +end + +function A_den = Conn_den(a1,a2,a3) %#ok + + A_den = repmat(ones(size(a1)),[3,3]); + +end \ No newline at end of file diff --git a/UserFiles/v4/rlhatton/Systems/sysf_honey_swimmer_triangle.m b/UserFiles/v4/rlhatton/Systems/sysf_honey_swimmer_triangle.m new file mode 100644 index 0000000..f7b79ce --- /dev/null +++ b/UserFiles/v4/rlhatton/Systems/sysf_honey_swimmer_triangle.m @@ -0,0 +1,139 @@ +function output = sysf_honey_swimmer_triangle(input_mode) + + % Default arguments + if ~exist('input_mode','var') + + input_mode = 'initialize'; + + end + + if ~exist('datapath','var') + + datapath = ''; + end + + %%%%%%% + + switch input_mode + + case 'name' + + output = 'Low Re triangular wave one period'; % Display name + + case 'dependency' + + output.dependency = {'Utilities/curvature_mode_toolbox/backbone_from_general_curvature.m'; + 'Utilities/curvature_mode_toolbox/curvatures/curv_serpenoid_var_sharness'; + 'Utilities/LowRE_toolbox/LowRE_dissipation_metric_from_general_curvature.m'; + 'Utilities/LowRE_toolbox/LowRE_local_connection_from_general_curvature.m'}; + + case 'initialize' + + %Initialize a kinematic snake with unit values for L and R + + %%%%% + % Filename to save to + output = mfilename; + + + %Functional representation of local connection + s.A_num = @Conn_num; + s.A_den = @Conn_den; + + + %%% + %Processing details + + %Mark that system has a singularity that needs to be accounted for + s.singularity = 1; + + %Range over which to evaluate connection + s.grid_range = [-1,1,-1,1]*12; + + %densities for various operations + s.density.vector = [1 1]*11; %density to display vector field + s.density.scalar = [1 1]*21; %density to display scalar functions + s.density.eval = [1 1]*21; %density for function evaluations + s.finite_element_density = 11; + % power metric + s.metric = @M_helper; %(x,y) LowRE_dissipation_metric_from_general_curvature... + %(@curv_serpenoid_var_sharpness,[x;y;10],1,1,2); % @(x,y) eye(2);% + + %%% + %Display parameters + + %shape space tic locations + s.tic_locs.x = [-1 0 1]*6; + s.tic_locs.y = [-1 0 1]*6; + + + %%%% + %Save the system properties + output = s; + + + end + +end + +function [Ar]=Conn_num(a1,a2) + + + %Apply the inverse multiplication + wb = waitbar2a(0,['Building ' num2str(size(a1)) ' connection matrix']); + Ar_cell = cell(size(a1)); + parfor i = 1:numel(a1); + Ar_cell{i} = A_num_helper(a1(i),a2(i)); + waitbar2a(i/numel(a1)); + end + + close(wb) + + %Ar_cell = arrayfun(@(a1,a2) A_num_helper(a1,a2),a1,a2,'UniformOutput',false); + + %Pull Ar values into a matrix + Ar_woven = cell2mat(Ar_cell); + + %Rearrange the A matrix + Ar = [ Ar_woven(1:3:end,1:2:end) Ar_woven(1:3:end,2:2:end); + Ar_woven(2:3:end,1:2:end) Ar_woven(2:3:end,2:2:end); + Ar_woven(3:3:end,1:2:end) Ar_woven(3:3:end,2:2:end)]; + +end + +function A = A_num_helper(a1,a2) + + % Get the local connection for the specified shape, with unit length + A = LowRE_local_connection_from_general_curvature(@curv_triangle_wave,[a1;a2],1,1,2); + + % Make any NaN values in the connection the mean of the neighboring + % values + A = inpaint_nans(A); + +end + +function A_den = Conn_den(a1,a2) %#ok + + A_den = repmat(ones(size(a1)),[3,2]); + +end + +function M = M_helper(a1,a2) + + % Get the local connection for the specified shape, with unit length + M = LowRE_dissipation_metric_from_general_curvature... + (@curv_triangle_wave,[a1;a2],1,1,2); + + % Make any NaN values in the metric the mean of the neighboring values + M = inpaint_nans(M); + +end + + +% function output = curvelock(params,mode) +% +% % lock the variable frequency serpenoid curvature to frequency 1 +% +% output = curv_serpenoid_var_sharpness(params,mode); +% +% end \ No newline at end of file diff --git a/UserFiles/v4/rlhatton/Systems/sysf_mudskipper.m b/UserFiles/v4/rlhatton/Systems/sysf_mudskipper.m new file mode 100644 index 0000000..e40038a --- /dev/null +++ b/UserFiles/v4/rlhatton/Systems/sysf_mudskipper.m @@ -0,0 +1,144 @@ +function output = sysf_mudskipper(input_mode) +% This file is an example of how to load the local connection as a data +% file, instead of generating it programmatically. + + + % Default arguments + if ~exist('input_mode','var') + + input_mode = 'initialize'; + + end + + + % File name where local connection is stored. This should be a file + % with: + % + % alpha1 and alpha2: ndgrid matrices (*not* meshgrid matrices) + % Ax1 through Atheta2: each component of the local connection, + % evaluated at all grid points + Local_connection = 'constrained_high_adduction_20'; + + % File name where power-usage data is stored + %Metric = 'Granular_Metric_Tensor.mat'; %Uncomment this line if loading a metric tensor + + + %%%%%%% + + switch input_mode + + case 'name' + + output = 'Mud Skipper'; % Display name + + case 'dependency' + + output.dependency = {Local_connection}; + + case 'initialize' + + %Initialize a kinematic snake with unit values for L and R + + %%%%% + % Filename to save to + output = mfilename; + + + % Load the connection and metric data + load(Local_connection,'A','alpha1','alpha2') + %load(Metric) + + %Functional representation of local connection + s.A_num = @(a1,a2) Conn_num(a1,a2,alpha1',alpha2',A{1}',A{4}',A{2}',A{5}',A{3}',A{6}'); + s.A_den = @Conn_den; + + + %%% + %Processing details + + %Mark that system has a singularity that needs to be accounted for + s.singularity = 0; + + %Range over which to evaluate connection + s.grid_range = [0,1,0,1.5]; + + %densities for various operations + s.density.vector = [11 11]; %density to display vector field + s.density.scalar = [21 21]; %density to display scalar functions + s.density.eval = [11 11]; %density for function evaluations + s.finite_element_density = 11; + s.density.metriceval = [11 11]; %density for metric evaluation + + % power metric + %s.metric = eye(2); %@(x,y) Granular_metric_calc(x,y,Metric_Tensor_raw,alpha1,alpha2); + + %%% + %Display parameters + + %shape space tic locations + s.tic_locs.x = [0 1]; + s.tic_locs.y = [0 1.5]; + + + %%%% + %Save the system properties + output = s; + + + end + +end + +function [Ar]=Conn_num(a1,a2,alpha1,alpha2,Ax1,Ax2,Ay1,Ay2,Atheta1,Atheta2) +% This function takes the components of the connection and arranges them +% into the internal structure used by the sysplotter back-end code. +% +% To improve the quality of the results, natural symmetries in the geometry +% of the system can be exploited to average the connection terms that +% should be the same as each other. +% +% There are three options for doing this, depending on the kind of symmetry +% expected in the data: +% +% 3link: If the system has three links with equal first and last links +% EvenOdd: If the system shape variables are for even and odd modes +% None: Don't average the connection across any points, and display exactly +% the data in the .mat file + +% Select averaging mode +% averaging_mode = '3link'; +% averaging_mode = 'EvenOdd'; + averaging_mode = 'None'; + + +Ar = four_way_symmetry(alpha1, alpha2, Ax1, Ax2, Ay1, Ay2, Atheta1,Atheta2, a1, a2, 1/3, averaging_mode); + +% else +% Ar_cell = A_num_helper(a1,a2); +% +% %Pull Ar values into a matrix +% Ar_woven = cell2mat(Ar_cell); +% +% %Rearrange the A matrix +% Ar = [ Ar_woven(1:3:end,1:2:end) Ar_woven(1:3:end,2:2:end); +% Ar_woven(2:3:end,1:2:end) Ar_woven(2:3:end,2:2:end); +% Ar_woven(3:3:end,1:2:end) Ar_woven(3:3:end,2:2:end)]; +% +% end + +end + +function A = A_num_helper(a1,a2) + + + % Get the local connection for the specified shape, with unit length + + A = feval(@Main_Granular_Swimmer,a2, a1); + +end + +function A_den = Conn_den(a1,a2) %#ok + + A_den = repmat(ones(size(a1)),[3,2]); + +end \ No newline at end of file diff --git a/UserFiles/v4/rlhatton/Systems/sysf_mudskipper_00.m b/UserFiles/v4/rlhatton/Systems/sysf_mudskipper_00.m new file mode 100644 index 0000000..e6cfca0 --- /dev/null +++ b/UserFiles/v4/rlhatton/Systems/sysf_mudskipper_00.m @@ -0,0 +1,144 @@ +function output = sysf_mudskipper_00(input_mode) +% This file is an example of how to load the local connection as a data +% file, instead of generating it programmatically. + + + % Default arguments + if ~exist('input_mode','var') + + input_mode = 'initialize'; + + end + + + % File name where local connection is stored. This should be a file + % with: + % + % alpha1 and alpha2: ndgrid matrices (*not* meshgrid matrices) + % Ax1 through Atheta2: each component of the local connection, + % evaluated at all grid points + Local_connection = 'Ross_vector_fields_00'; + + % File name where power-usage data is stored + %Metric = 'Granular_Metric_Tensor.mat'; %Uncomment this line if loading a metric tensor + + + %%%%%%% + + switch input_mode + + case 'name' + + output = 'Mud Skipper 00'; % Display name + + case 'dependency' + + output.dependency = {Local_connection}; + + case 'initialize' + + %Initialize a kinematic snake with unit values for L and R + + %%%%% + % Filename to save to + output = mfilename; + + + % Load the connection and metric data + load(Local_connection,'X','Y','vecX','vecY') + %load(Metric) + + %Functional representation of local connection + s.A_num = @(a1,a2) Conn_num(a1,a2,X,Y,vecX,vecY,zeros(size(X)),zeros(size(X)),zeros(size(X)),zeros(size(X))); + s.A_den = @Conn_den; + + + %%% + %Processing details + + %Mark that system has a singularity that needs to be accounted for + s.singularity = 0; + + %Range over which to evaluate connection + s.grid_range = [.1,1,0,1.5]; + + %densities for various operations + s.density.vector = [11 11]; %density to display vector field + s.density.scalar = [21 21]; %density to display scalar functions + s.density.eval = [11 11]; %density for function evaluations + s.finite_element_density = 11; + s.density.metriceval = [11 11]; %density for metric evaluation + + % power metric + %s.metric = eye(2); %@(x,y) Granular_metric_calc(x,y,Metric_Tensor_raw,alpha1,alpha2); + + %%% + %Display parameters + + %shape space tic locations + s.tic_locs.x = [0 1]; + s.tic_locs.y = [0 1.5]; + + + %%%% + %Save the system properties + output = s; + + + end + +end + +function [Ar]=Conn_num(a1,a2,alpha1,alpha2,Ax1,Ax2,Ay1,Ay2,Atheta1,Atheta2) +% This function takes the components of the connection and arranges them +% into the internal structure used by the sysplotter back-end code. +% +% To improve the quality of the results, natural symmetries in the geometry +% of the system can be exploited to average the connection terms that +% should be the same as each other. +% +% There are three options for doing this, depending on the kind of symmetry +% expected in the data: +% +% 3link: If the system has three links with equal first and last links +% EvenOdd: If the system shape variables are for even and odd modes +% None: Don't average the connection across any points, and display exactly +% the data in the .mat file + +% Select averaging mode +% averaging_mode = '3link'; +% averaging_mode = 'EvenOdd'; + averaging_mode = 'None'; + + +Ar = four_way_symmetry(alpha1, alpha2, Ax1, Ax2, Ay1, Ay2, Atheta1,Atheta2, a1, a2, 1/3, averaging_mode); + +% else +% Ar_cell = A_num_helper(a1,a2); +% +% %Pull Ar values into a matrix +% Ar_woven = cell2mat(Ar_cell); +% +% %Rearrange the A matrix +% Ar = [ Ar_woven(1:3:end,1:2:end) Ar_woven(1:3:end,2:2:end); +% Ar_woven(2:3:end,1:2:end) Ar_woven(2:3:end,2:2:end); +% Ar_woven(3:3:end,1:2:end) Ar_woven(3:3:end,2:2:end)]; +% +% end + +end + +function A = A_num_helper(a1,a2) + + + % Get the local connection for the specified shape, with unit length + + A = feval(@Main_Granular_Swimmer,a2, a1); + +end + +function A_den = Conn_den(a1,a2) %#ok + + A_den = repmat(ones(size(a1)),[3,2]); + +end \ No newline at end of file diff --git a/UserFiles/v4/rlhatton/Systems/sysf_mudskipper_20.m b/UserFiles/v4/rlhatton/Systems/sysf_mudskipper_20.m new file mode 100644 index 0000000..1cab12f --- /dev/null +++ b/UserFiles/v4/rlhatton/Systems/sysf_mudskipper_20.m @@ -0,0 +1,144 @@ +function output = sysf_mudskipper_20(input_mode) +% This file is an example of how to load the local connection as a data +% file, instead of generating it programmatically. + + + % Default arguments + if ~exist('input_mode','var') + + input_mode = 'initialize'; + + end + + + % File name where local connection is stored. This should be a file + % with: + % + % alpha1 and alpha2: ndgrid matrices (*not* meshgrid matrices) + % Ax1 through Atheta2: each component of the local connection, + % evaluated at all grid points + Local_connection = 'Ross_vector_fields_20'; + + % File name where power-usage data is stored + %Metric = 'Granular_Metric_Tensor.mat'; %Uncomment this line if loading a metric tensor + + + %%%%%%% + + switch input_mode + + case 'name' + + output = 'Mud Skipper 20'; % Display name + + case 'dependency' + + output.dependency = {Local_connection}; + + case 'initialize' + + %Initialize a kinematic snake with unit values for L and R + + %%%%% + % Filename to save to + output = mfilename; + + + % Load the connection and metric data + load(Local_connection,'X','Y','vecX','vecY') + %load(Metric) + + %Functional representation of local connection + s.A_num = @(a1,a2) Conn_num(a1,a2,X,Y,vecX,vecY,zeros(size(X)),zeros(size(X)),zeros(size(X)),zeros(size(X))); + s.A_den = @Conn_den; + + + %%% + %Processing details + + %Mark that system has a singularity that needs to be accounted for + s.singularity = 0; + + %Range over which to evaluate connection + s.grid_range = [0.1,1,0,1.5]; + + %densities for various operations + s.density.vector = [11 11]; %density to display vector field + s.density.scalar = [21 21]; %density to display scalar functions + s.density.eval = [11 11]; %density for function evaluations + s.finite_element_density = 11; + s.density.metriceval = [11 11]; %density for metric evaluation + + % power metric + %s.metric = eye(2); %@(x,y) Granular_metric_calc(x,y,Metric_Tensor_raw,alpha1,alpha2); + + %%% + %Display parameters + + %shape space tic locations + s.tic_locs.x = [0 1]; + s.tic_locs.y = [0 1.5]; + + + %%%% + %Save the system properties + output = s; + + + end + +end + +function [Ar]=Conn_num(a1,a2,alpha1,alpha2,Ax1,Ax2,Ay1,Ay2,Atheta1,Atheta2) +% This function takes the components of the connection and arranges them +% into the internal structure used by the sysplotter back-end code. +% +% To improve the quality of the results, natural symmetries in the geometry +% of the system can be exploited to average the connection terms that +% should be the same as each other. +% +% There are three options for doing this, depending on the kind of symmetry +% expected in the data: +% +% 3link: If the system has three links with equal first and last links +% EvenOdd: If the system shape variables are for even and odd modes +% None: Don't average the connection across any points, and display exactly +% the data in the .mat file + +% Select averaging mode +% averaging_mode = '3link'; +% averaging_mode = 'EvenOdd'; + averaging_mode = 'None'; + + +Ar = four_way_symmetry(alpha1, alpha2, Ax1, Ax2, Ay1, Ay2, Atheta1,Atheta2, a1, a2, 1/3, averaging_mode); + +% else +% Ar_cell = A_num_helper(a1,a2); +% +% %Pull Ar values into a matrix +% Ar_woven = cell2mat(Ar_cell); +% +% %Rearrange the A matrix +% Ar = [ Ar_woven(1:3:end,1:2:end) Ar_woven(1:3:end,2:2:end); +% Ar_woven(2:3:end,1:2:end) Ar_woven(2:3:end,2:2:end); +% Ar_woven(3:3:end,1:2:end) Ar_woven(3:3:end,2:2:end)]; +% +% end + +end + +function A = A_num_helper(a1,a2) + + + % Get the local connection for the specified shape, with unit length + + A = feval(@Main_Granular_Swimmer,a2, a1); + +end + +function A_den = Conn_den(a1,a2) %#ok + + A_den = repmat(ones(size(a1)),[3,2]); + +end \ No newline at end of file