function reactiontime_later_reciprobit_fig5()
% Fig. 5-style LATER parameter changes + reciprobit plots
% Noorani & Carpenter (2016) 

close all; clc;

%% --- Baseline parameters
N       = 5000;
S0      = 0;
Delta0  = 1.0;
mu0     = 4.0;
sig0    = 0.6;
delay   = 0.00;

% schematic time axis
tmax = 0.35; dt = 0.001;
t = 0:dt:tmax;

%% --- Manipulations (3 conditions each)
% Row 1: vary Delta
Delta_set1 = Delta0 * [0.7 1.0 1.4];
mu_set1    = mu0;                 % scalar is OK now
sig_set1   = sig0 * [1 1 1];

% Row 2: vary mu
Delta_set2 = Delta0;              % scalar OK
mu_set2    = mu0   * [0.75 1.0 1.35];
sig_set2   = sig0;

% Row 3: vary sigma
Delta_set3 = Delta0;
mu_set3    = mu0;
sig_set3   = sig0  * [0.45 1.0 1.75];

%% --- Figure
figure(1);
clf

subplot(3,2,1);
plot_decision_schematics(t, S0, Delta_set1, mu_set1, ...
    'Vary \Delta (threshold range)');

subplot(3,2,2);
plot_reciprobit_from_sim(N, Delta_set1, mu_set1, sig_set1, delay, ...
    'Effect: swivel about intercept (k)');

subplot(3,2,3);
plot_decision_schematics(t, S0, Delta_set2, mu_set2, ...
    'Vary \mu (mean rate)');

subplot(3,2,4);
plot_reciprobit_from_sim(N, Delta_set2, mu_set2, sig_set2, delay, ...
    'Effect: horizontal, parallel shift');

subplot(3,2,5);
plot_sigma_schematics(t, S0, Delta_set3, mu_set3, sig_set3, ...
    'Vary \sigma (rate SD)');

subplot(3,2,6);
plot_reciprobit_from_sim(N, Delta_set3, mu_set3, sig_set3, delay, ...
    'Effect: rotate about median');




% optional, save figure
fname = fullfile('..','..','images',[mfilename '_ori.svg']);
savegraph(fname,'svg');
end

%% ===== Helper: decision schematic (supports scalar or vector mu/Delta) =====
function plot_decision_schematics(t, S0, Delta_in, mu_in, ttl)
hold on;

% Ensure both are vectors of the same length
n = max(numel(Delta_in), numel(mu_in));
Delta_vec = expand_to_len(Delta_in, n);
mu_vec    = expand_to_len(mu_in,    n);

for i = 1:n
    ST = S0 + Delta_vec(i);
    S  = min(ST, S0 + mu_vec(i).*t);
    plot(t*1000, S, 'LineWidth', 2);
end

% baseline (middle condition if >=3, otherwise first)
ibase = min(2, n);
yline(S0 + Delta_vec(ibase),'k--','Threshold','LabelVerticalAlignment','bottom');
xline(0,'k-','Stimulus onset');

xlabel('Time (ms)'); ylabel('Decision signal S');
title(ttl);
box off;
end

%% ===== Helper: sigma schematic (mean line + +/-2sigma rate bands) =====
function plot_sigma_schematics(t, S0, Delta_in, mu_in, sig_in, ttl)
hold on;

Delta0 = Delta_in(1);
mu0    = mu_in(1);
ST     = S0 + Delta0;

% mean line (rate = mu0)
Smean = min(ST, S0 + mu0.*t);
plot(t*1000, Smean, 'k-', 'LineWidth', 2);

% ensure sigma vector
sig_vec = sig_in(:)';

for i = 1:numel(sig_vec)
    rlo = max(1e-6, mu0 - 2*sig_vec(i));
    rhi = max(1e-6, mu0 + 2*sig_vec(i));
    Slo = min(ST, S0 + rlo.*t);
    Shi = min(ST, S0 + rhi.*t);

    patch([t*1000 fliplr(t*1000)], [Slo fliplr(Shi)], 0.85*[1 1 1], ...
        'EdgeColor','none','FaceAlpha',0.18);
end

yline(ST,'k--','Threshold','LabelVerticalAlignment','bottom');
xline(0,'k-','Stimulus onset');

xlabel('Time (ms)'); ylabel('Decision signal S');
title(ttl);
box off;
end

%% ===== Helper: reciprobit plot from simulation (supports scalar or vector params) =====
function plot_reciprobit_from_sim(N, Delta_in, mu_in, sig_in, delay, ttl)
hold on;

n = max([numel(Delta_in), numel(mu_in), numel(sig_in)]);
Delta_vec = expand_to_len(Delta_in, n);
mu_vec    = expand_to_len(mu_in,    n);
sig_vec   = expand_to_len(sig_in,   n);

for i = 1:n
    Delta = Delta_vec(i);
    mu    = mu_vec(i);
    sig   = sig_vec(i);

    r = mu + sig.*randn(N,1);
    r(r < 1e-6) = 1e-6;            % truncate negatives for clean curves

    T = delay + (Delta ./ r);      % seconds
    x = sort(-1 ./ T);              % promptness (1/s)

    p = ((1:N)' - 0.5) ./ N;
    z = probit(p);

    plot(x, z, 'LineWidth', 2);
end

xlabel('Promptness 1/T (s^{-1})');
ylabel('Probit (Z)');
title(ttl);
grid on; box off;
end

%% ===== Utilities =====
function out = expand_to_len(x, n)
% Replicate scalar to length n; if already length n, keep; if other length, error.
if isscalar(x)
    out = repmat(x, 1, n);
else
    x = x(:)'; % row
    if numel(x) ~= n
        error('Parameter vector has length %d but expected %d.', numel(x), n);
    end
    out = x;
end
end

function z = probit(p)
% Inverse normal CDF without toolboxes
p = min(max(p, 1e-12), 1-1e-12);
z = sqrt(2) * erfinv(2*p - 1);
end