Mercurial > repos > public > sbplib

changeset 104:ffd9e68f63cc feature/arclen-param

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Splined in the other direction! Works great. Runs faster. Removed ugly fix outside [0,1] in spline(). Camelcased numerical_derivative().
author Martin Almquist <martin.almquist@it.uu.se>
date Mon, 07 Dec 2015 18:46:59 +0100
parents bc5db54f9efd
children 1df4f3704b76
files +grid/Curve.m +util/fzero_vec.m +util/integral_vec.m
diffstat 3 files changed, 5 insertions(+), 68 deletions(-) [+]
line wrap: on
line diff
diff -r bc5db54f9efd -r ffd9e68f63cc +grid/Curve.m
--- a/+grid/Curve.m	Mon Dec 07 17:24:28 2015 +0100
+++ b/+grid/Curve.m	Mon Dec 07 18:46:59 2015 +0100
@@ -8,7 +8,6 @@
     methods
         %TODO:
         % Errors or FD if there is no derivative function added.
-        % -semi-done
 
         % Concatenation of curves
         % Subsections of curves
@@ -98,11 +97,9 @@
             % Construct arcLength function using splines
             tvec = linspace(0,1,N);
             arcVec = obj.arcLength(0,tvec);
-            arcLength = spline(tvec,arcVec);
-   
-            % Stretch the parameter, construct function with splines
-            arcParVec = fzero_vec(@(s)arcLength(s) - tvec*arcLength(1),[0-10*eps,1+10*eps]);
-            arcPar = spline(tvec,arcParVec);
+            tFunc = spline(arcVec,tvec); % t as a function of arcLength
+            L = obj.arcLength(0,1);
+            arcPar = @(s) tFunc(s*L);
             
             % New function and derivative
             g_new = @(t)obj.g(arcPar(t));
@@ -249,7 +246,8 @@
 
     methods (Static)
         
-        function gp_out = numerical_derivative(g,D1)
+        
+        function gp_out = numericalDerivative(g,D1)
             m = length(D1); L = 1; % Assume curve parameter from 0 to 1.
             t = linspace(0,L,m); 
             gVec = g(t)';
@@ -376,33 +374,12 @@
     end
 end
 
-function I = fzero_vec(f,lim)
-% Wrapper around the built-in function fzero that
-% handles multiple functions (vector-valued f).
-
-    fval = f(lim(1));
-    I = zeros(size(fval));
-    
-    for i = 1:length(fval)
-        e = zeros(size(fval));
-        e(i) = 1;
-        I(i) = fzero(@(t) sum(e.*f(t)),lim);
-    end
-    
-end
-
 % Returns a function handle to the spline.
 function f = spline(tval,fval,spline_order)
     default_arg('spline_order',4);
     [m,~] = size(tval);
     assert(m==1,'Need row vectors.');
 
-    % make vectors longer to be safe slightly beyond edges.
-    dt0 = tval(2)-tval(1); dt1 = tval(end)-tval(end-1);
-    df0 = fval(2)-fval(1); df1 = fval(end)-fval(end-1);
-    tval = [tval(1)-dt0, tval, tval(end)+dt1];
-    fval = [fval(1)-df0, fval, fval(end)+df1];
-
     f_spline = spapi( optknt(tval,spline_order), tval, fval );
     f = @(t) fnval(f_spline,t);
 end
diff -r bc5db54f9efd -r ffd9e68f63cc +util/fzero_vec.m
--- a/+util/fzero_vec.m	Mon Dec 07 17:24:28 2015 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,14 +0,0 @@
-function I = fzero_vec(f,lim)
-% Wrapper around the built-in function fzero that
-% handles multiple functions (vector-valued f).
-
-    fval = f(lim(1));
-    I = zeros(size(fval));
-    
-    for i = 1:length(fval)
-        e = zeros(size(fval));
-        e(i) = 1;
-        I(i) = fzero(@(t) sum(e.*f(t)),lim);
-    end
-    
-end
\ No newline at end of file
diff -r bc5db54f9efd -r ffd9e68f63cc +util/integral_vec.m
--- a/+util/integral_vec.m	Mon Dec 07 17:24:28 2015 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,26 +0,0 @@
-function I = integral_vec(f,a,b)
-% Wrapper around the built-in function integral that
-% handles multiple limits.
-
-    Na = length(a);
-    Nb = length(b);
-    assert(Na == 1 || Nb == 1 || Na==Nb,...
-        'a and b must have same length, unless one is a scalar.');
-    
-    if(Na>Nb);
-        I = zeros(size(a));
-        for i = 1:Na
-            I(i) = integral(f,a(i),b);
-        end
-    elseif(Nb>Na)
-        I = zeros(size(b));
-        for i = 1:Nb
-            I(i) = integral(f,a,b(i));
-        end
-    else
-        I = zeros(size(b));
-        for i = 1:Nb
-            I(i) = integral(f,a(i),b(i));
-        end
-    end
-end
\ No newline at end of file