function Resistance_at_end_of_travel(varargin)
% Emulates (not simulates) the resitive force profile of a latching
% mechanism (e.g., a four-bar latch).  Most mechanical latches see a peak
% of that force near the point where they effect the latch.  Some
% implementations that are used to lock both the deployed and undeployed
% configuration of their parent mechanism can see a peak at (or near) both
% extremes of travel, with two-headed peaking as is shown here.
%
%This m-file is (obviously) just more "fun with numbers" (FWN), using the erfc
%function to create surges for the two peaks.  The variables of the
%"stroke" subfunction are used to tailor the sharpness and location of the
%peaks.  Although FWN, it is actually no less specific than a qualitative
%requirement on the subject of resistive load peaking.

%Inputs:
    %t: 1/stepsize
    %q: perturbation scaling constant

persistent h;

t_default = 10000;
q_default = 75;

switch nargin
    case 2
        t = varargin{1};
        q = varargin{2};
    case 1
        t = t_default;
        q = varargin{1};
    case 0
        t = t_default;
        q = q_default;
end;
        
x = 0:1/t:1;

if isempty(h) %Open new figure windows...
    h(1,1) = figure;
else %...or find the existing, re-opening with the same handle if they got closed.
    h(1,1) = figure(h(1,1));
end;

%Generate plot pairs: one with, and one without, randomized variation for
%each of four input value sets.
stroke(x,0.925,0,5,4,h,1,'a');
stroke(x,0.925,q,5,4,h,2,'b');

stroke(x,0.95,0,50,40,h,3,'c');
stroke(x,0.95,q,50,40,h,4,'d');

stroke(x,0.95,0,400,500,h,5,'e');
stroke(x,0.95,q,400,500,h,6,'f');

stroke(x,0.995,0,400,500,h,7,'g');
stroke(x,0.995,q,400,500,h,8,'h');

end

function y = stroke(x,oSet,q,r1,r2,fig,sp,ttl)

oSet1 = oSet;
oSet2 = 1 - oSet;

vSet = 0.5;

if q ~= 0
    oSet1 = oSet1 * (1 + (1 - oSet1) * vSet * randn(1,1));
    oSet2 = oSet2 * (1 + (1 - oSet2) * vSet * randn(1,1));
end;

y1 = erfc(-(x-oSet1)*r1);%one end of travel...
y2 = erfc(-(x-oSet2)*r2);%...the other end of travel

y = [0,diff(y1 + y2)];
y = y + min(y(y>0.01*max(y)));
y = abs(y);

if q ~= 0
    y = y .* (1 + randn(size(y))/q);
end;

y = y/max(y);

figure(fig);
subplot(4,2,sp),plot(x,y);

if ismember(sp,[1,3,5,7])
    ylabel 'Normalized Resistive Load (ndim)'; 
end;
title(ttl);
if ismember(sp,[7,8])
    xlabel 'Normalized Position (ndim)';
end;
end