function Sigmoid_noise(n_crit,scr_FOM_crit,singlePlot)
%This is a purely behavioral "model" (and a crude one at that!).  It
%emulates (not simulates!) the aggregate effect of things like motion,
%friction and strain energy, but does not directly represent any of them.
%It is (therefore) useful ONLY to demonstrate the appearance of how we want
%things to behave for discussion purposes, which includes parameterization
%of qualitative requirements.  It is NEVER the reality how we want them to
%"operate".
%
%I often refer this notion "fun with numbers".  When writing requirements
%for a widget to exhibit some specific behavior, I will rarely release any
%code such as this.  It makes it too easy for non-cognoscenti to think of
%it as a source of validation and verification data.

n_crit_default = 1;%Just a loop control criterion
scr_FOM_crit_default = 200;%Lower is better
scr_plot_default = true;%Convenience

switch nargin
    case 0
        n_crit = n_crit_default;
        scr_FOM_crit = scr_FOM_crit_default;
        singlePlot = scr_FOM_crit_default;
    case 1
        if floor(n_crit) == n_crit
            scr_FOM_crit = scr_FOM_crit_default;
        else
            error('number of passes must be an integer.');
        end;
        singlePlot = scr_FOM_crit_default;
    case 2
        if isempty(n_crit)
            n_crit = n_crit_default;
        end;
        if isempty(scr_FOM_crit)
            scr_FOM_crit = scr_FOM_crit_default;
        end;
        singlePlot = scr_FOM_crit_default;
    case 3
        if isempty(n_crit)
            n_crit = n_crit_default;
        end;
        if isempty(scr_FOM_crit)
            scr_FOM_crit = scr_FOM_crit_default;
        end;
        if isempty(singlePlot)
            singlePlot = scr_plot_default;
        end;
    otherwise
        error('Three inputs required: number of passes, Figure of Merit criterion.');
end;

%% Preliminaries

%A time step magnitude
scr_dt = 1/199;

%A simple time vector, as if "normalized" on the time it takes to actuate
%between stops.  Double-sized to start with; will adaptively truncate below.
scr_t = 0:scr_dt:2; 

%An arbitrary affine "stretch" factor to tailor the shape.  Lower values
%make the curve more like a straight line.  Higher values make the main
%rise steeper, but delay the transition into the rise.
scr_k = 4;

%Drag factor, as if dynamic drag (e.g., friction).  This will drive both an
%overall retardation and some randomness.
scr_d = 0.5;

%Ratio of as-if-static drag to as-if-dynamic drag.  Makes it harder to get
%moving again once the mechanism has stopped.
scr_st = 5;

%Windup factor, as if releasing strain energy accumulated when the actuator
%was applying load, but the mechanism wasn't moving (i.e., without a
%clutch).  Makes it easier to keep moving once started again after
%stopping.
scr_w = 0.4;

%% Generate a smooth trajectory, as if of a nominal kinematic stroke...
%...including the effect of any preload on the mechanism.
scr_s = 1./(1+exp(-scr_k*(2*scr_t-1)));

scr_s = scr_s - scr_s(1);%Start at position zero, as if at a hard stop
scr_s = scr_s / (scr_s(scr_t == 1));%Scale to finish at time == 1
scr_s(scr_s>1) = 1;%Trim, as if to a hard stop

scr_v = [0,diff(scr_s)]/scr_dt;%kinematic speed (starts at zero)

%Init loop
scr_FOM = scr_FOM_crit * 1.1;
scr_FOM_best = [];

scr_cutoff = -1;
scr_cutoff_crit = 0;
scr_cutoff_best = [];

n = 1;

%A crude loop, as if running a wee little Monte Carlo analysis
while ((scr_FOM > scr_FOM_crit || scr_cutoff < scr_cutoff_crit) && n <= n_crit)
    [scr_SN,scr_vnsw,scr_sn_f] = noisy(scr_dt,scr_k,scr_t,scr_d,scr_st,scr_w,scr_s,scr_v);
    [scr_pxx_kin,scr_freq_kin,scr_pxx_noise,scr_freq_noise,scr_diff,scr_cutoff,scr_cumsum,scr_FOM] = chk(scr_v,scr_vnsw);

    if isempty(scr_FOM_best) || isempty(scr_cutoff_best)
        scr_SN_best = scr_SN;
        scr_sn_f_best = scr_sn_f;
        scr_pxx_kin_best = scr_pxx_kin;
        scr_freq_kin_best = scr_freq_kin;
        scr_pxx_noise_best = scr_pxx_noise;
        scr_freq_noise_best = scr_freq_noise;
        scr_diff_best = scr_diff;
        scr_cutoff_best = scr_cutoff;
        scr_cumsum_best = scr_cumsum;
        scr_FOM_best = scr_FOM;
    elseif scr_FOM <= scr_FOM_best && scr_cutoff >= scr_cutoff_best
        scr_SN_best = scr_SN;
        scr_sn_f_best = scr_sn_f;
        scr_pxx_kin_best = scr_pxx_kin;
        scr_freq_kin_best = scr_freq_kin;
        scr_pxx_noise_best = scr_pxx_noise;
        scr_freq_noise_best = scr_freq_noise;
        scr_diff_best = scr_diff;
        scr_cutoff_best = scr_cutoff;
        scr_cumsum_best = scr_cumsum;
        scr_FOM_best = scr_FOM;
    end;
    n = n + 1;
end;

if n <= n_crit
    scr_SN_best = scr_SN;
    scr_sn_f_best = scr_sn_f;
    scr_pxx_kin_best = scr_pxx_kin;
    scr_freq_kin_best = scr_freq_kin;
    scr_pxx_noise_best = scr_pxx_noise;
    scr_freq_noise_best = scr_freq_noise;
    scr_diff_best = scr_diff;
    scr_cutoff_best = scr_cutoff;
    scr_cumsum_best = scr_cumsum;
    scr_FOM_best = scr_FOM;
else
    n = n-1;
end
%% See what I did there?

scr_t = scr_t(1:scr_sn_f_best);
scr_s = scr_s(1:scr_sn_f_best);
scr_SN_best = scr_SN_best(1:scr_sn_f_best);

if singlePlot
    figure;
    plot(scr_t,[scr_s;scr_SN_best]);
else
    subplot(4,1,1),plot(scr_t,[scr_s;scr_SN_best]);
end;

title(sprintf('Notional Trajectory \nMean Perturbed Velocity %1.3f\n(k = %2.1f; d = %3.1f; st = %2.0f; w = %2.1f)',1/scr_t(scr_sn_f_best),scr_k,scr_d,scr_st,scr_w));
legend('Nominal','Perturbed','Location','SouthEast');
ylabel 'Normalized Position (ndim)';
xlabel(sprintf('Normalized Time (ndim)'));
grid on;

if singlePlot
    figure;
    plot(scr_freq_kin_best/pi,10*log10(scr_pxx_kin_best));
else
    subplot(4,1,2),plot(scr_freq_kin_best/pi,10*log10(scr_pxx_kin_best));
end;

title 'Welch Power Spectral Density Estimates';
ylabel 'Pwr/freq (dB/rad/smpl)';
xlabel 'Normalized Frequency (\times\pi rad/sample)';

hold on;

scr_c_o1 = get(gca,'ColorOrder');
set(gca,'ColorOrder',[scr_c_o1(2:end,:);scr_c_o1(1,:)]);
if singlePlot
    plot(scr_freq_noise_best/pi,10*log10(scr_pxx_noise_best));
else
    subplot(4,1,2),plot(scr_freq_noise_best/pi,10*log10(scr_pxx_noise_best));
end;

hold off;

if singlePlot
    figure;
    plot(scr_freq_kin_best/pi,scr_diff_best);
else
    subplot(4,1,3),plot(scr_freq_kin_best/pi,scr_diff_best);
end;

title(sprintf('Impact on Welch Power Spectral Density Estimate'));
xlabel 'Normalized Frequency (\times\pi rad/sample)';
ylabel 'Pwr/freq (dB/rad/smpl)';
ylim([-20,20]);

if singlePlot
    figure;
    plot(scr_freq_kin_best/pi,cumsum(scr_cumsum_best));
else
    subplot(4,1,4),plot(scr_freq_kin_best/pi,cumsum(scr_cumsum_best));
end;

xlabel 'Normalized Frequency (\times\pi rad/sample)';
ylabel 'Pwr/freq (dB/rad/smpl)';
title(sprintf('Cumulative Power Discrepancy at Normalized Frequencies Below %0.2f (FOM = %3.0f dB)',scr_cutoff_best,scr_FOM_best));

if n_crit == 1
    display('Single-draw execution complete.');
else
    if n < n_crit
        display(sprintf('Success at %d iterations',n));
    else
        error('Criteria not met at n = %d.  Best results shown.',n);
    end;
end;

end

function [scr_SN,scr_vnsw,scr_sn_f] = noisy(scr_dt,scr_k,scr_t,scr_d,scr_st,scr_w,scr_s,scr_v)

%% Invent some noise, as if dynamic friction
scr_sn = 1./(1+exp(-scr_k*(2*scr_t-(1 + scr_d/(exp(1)^1.5)))));%Arbitrary
scr_sn = scr_sn - scr_sn(1);%Shift to start at zero position
%Will vertically scale after all noise sources are added (below)

scr_vn = [0,diff(scr_sn)]/scr_dt;

%A "noise" vector purturbing the velocity
scr_n =  1 - scr_d * abs(randn(size(scr_v)));
scr_vn = scr_vn .* scr_n;

%% Invent more noise, as if static drag ("stiction")

scr_stat = scr_st * rand(size(scr_vn));%Some arbitrary scaling.

scr_stat(scr_vn>scr_v(end)) = 0;%Applies only when stopped

%Impose the effect of static drag (scr_stat) on the retarded, randomized
%velocity.
scr_vns = scr_vn - scr_stat;
scr_vns(scr_vns<0) = 0;%But don't let it go backwards.

%% Invent some windup, as if the mechanism "popped"

%Windup is based on how far the mechanism would have gone with simple
%friction vs. the stoppage.
scr_windup = scr_w * (cumsum(scr_vn) - cumsum(scr_vns))*scr_dt;
scr_windup = [0,scr_windup(1:end-1)];%No windup at time 0

%While discharging windup, the mechanism can move faster than it would
%have if it had not stopped.
scr_vnsw = scr_vns + scr_windup;
scr_vnsw(scr_vnsw<0) = 0;

%% Accumulate the position based on the perturbed velocity

scr_SN = cumsum(scr_vnsw*scr_dt);%The perturbed position vector

scr_ds = scr_s(scr_t == 1) - scr_s(scr_t == (1 - scr_dt));%The last motion before hitting the stop
%scr_ds will be used as a "stopping" criterion below.

%Match (+/-) the last position change of the noisy trajectory to that of
%the original.  Could also have used velocity...doesn't much matter for
%demonstration purposes.
scr_sn_f = find(diff(scr_SN)>=scr_ds,1,'last')+1;
scr_SN = scr_SN/(scr_SN(scr_sn_f));%Scale to complete at position == 1

scr_SN(scr_SN>1) = 1;%Don't over-stroke
end

function [scr_pxx_kin,scr_freq_kin,scr_pxx_noise,scr_freq_noise,scr_diff,scr_cutoff,scr_cumsum,scr_FOM] = chk(scr_v,scr_vnsw)
%% As if analyzing the results of a test...
%Compare to the Kinematic Velocity (essentially, an error)

scr_wndw = 1:10;
[scr_pxx_kin,scr_freq_kin] = pwelch(scr_v,scr_wndw);
[scr_pxx_noise,scr_freq_noise] = pwelch(scr_vnsw,scr_wndw);

scr_diff = 10*(log10(scr_pxx_noise) - log10(scr_pxx_kin));
%Find the first cross-over between the two PSD estimates
%Because the kinematic trajectory includes the effect of any preload, the
%perturbed trajectory will always start more slowly, and will (therefore)
%have more power at low frequencies than the kinematic trajectory.  For
%mechanisms where this is NOT true, this algorithm will have to be
%modified.
scr_cutoff = scr_freq_noise(find(scr_diff > 0,1,'first'))/pi;
if isempty(scr_cutoff)
    scr_cutoff = 1;
end;

scr_cumsum = scr_diff;

scr_cumsum(scr_freq_kin/pi >  scr_cutoff) = 0;

scr_FOM = sum(abs(scr_cumsum))/scr_cutoff;
end