function [y, xRecip, stats,h] = reciprobitplot(x, varargin)
% RECIPROBITPLOT Plots a reciprobit graph of reaction time data.
%
%   [Y, XRECIP, STATS] = RECIPROBITPLOT(X) plots the reciprobit graph
%   of input data X (reaction times). It returns:
%     Y       - probit-transformed cumulative probabilities
%     XRECIP  - negative reciprocal of sorted reaction times
%     STATS   - optional regression statistics if requested
%
%   Optional Name-Value Pairs:
%     'print'    - (logical) export figure as PDF (default: false)
%     'Color'    - (char) plot color (default: 'w')
%     'xlim'     - (1x2 array) x-axis limits
%     'xtick'    - (array) x-axis tick locations
%     'regress'  - (logical) perform linear regression (default: false)
%
% Example:
%   reciprobitplot(reactionTimes, 'Color', 'k', 'regress', true);

%% Input parser
p = inputParser;
p.addRequired('x', @(x) isnumeric(x) && isvector(x));

p.addParameter('print', false, @(x) islogical(x) || isnumeric(x));
p.addParameter('Color', 'k');
p.addParameter('xlim', [], @(x) isempty(x) || (isnumeric(x) && numel(x)==2));
p.addParameter('xtick', sort(-1000./(100+[0 calcfl2fu(50,0:6)])), @(x) isnumeric(x));
% p.addParameter('xtick', sort(-1000./calcfl2fu(100,0:6)), @(x) isnumeric(x));

p.addParameter('ytick', [1 2 5 10:20:90 95 98 99]/100, @(x) isnumeric(x));
p.addParameter('regress', false, @(x) islogical(x) || isnumeric(x));

p.parse(x, varargin{:});

x			= p.Results.x;
printFlag   = logical(p.Results.print);
plotColor   = p.Results.Color;
xLimits     = p.Results.xlim;
xTicks      = p.Results.xtick;
yTicks      = p.Results.ytick;
regressFlag = logical(p.Results.regress);

%% Ensure reaction times are in ms
if ~any(x > 100)
    x = x * 1000; % Convert from seconds to milliseconds
end

%% Compute reciprocal and probit
x		= sort(x(:));                    % Ensure x is a sorted column vector
xRecip	= -1000 ./ x;              % Negative reciprocal (promptness)
n		= numel(xRecip);
y		= probit((1:n).' / n);          % Probit of cumulative distribution



%% Plotting
plot(xRecip, y, 'ko', 'MarkerFaceColor', plotColor); hold on;

%% Linear regression
stats = struct();
if regressFlag
	b			= regstats(y,xRecip,'linear',{'beta','rsquare'});
	% h			= regline(b.beta,'k-');
	h = refline(b.beta(2),b.beta(1));
	set(h,'Color',plotColor,'LineWidth',1);
	stats.b		= b.beta;
	stats.rsq	= b.rsquare;
hold on
end

xlabel('reaction time (ms))');
ylabel('cumulative probability (%)');

%% Relabel ticks
set(gca,'XTick',xTicks(2:end-1),'XTickLabel',-1000./xTicks(2:end-1));
if isempty(xLimits),  xLimits=[min(xTicks) max(xTicks)]; end
xlim(xLimits);

ptick	= probit(yTicks);
set(gca,'YTick',ptick,'YTickLabel',yTicks*100);
ylim([probit(0.1/100) probit(99.9/100)]);

xlabel('reaction time (ms)');
ylabel('cumulative probability');
nicegraph;

%% Test for normality
H		= lillietest(x);
stats.H = H;

%% Save figure if requested
if printFlag
    savegraph('reciprobitplot','svg');
end
end



function F = calcfl2fu(F1,oct)
% FU = calcfl2fu(FL,OCT)
%
% Determine frequency FU that lies OCT octaves from frequency FL
%

F = F1 .* 2.^oct;
end


function phi = probit(x)
    % PROBIT computes the probit function, which is the quantile function 
    % associated with the standard normal distribution.
    %
    % Usage:
    %   phi = PROBIT(x)
    %
    % Input:
    %   x - A numeric array where each element is in the interval [0, 1]
    %
    % Output:
    %   phi - The corresponding probit values
    %
    % Example:
    %   PROBIT(0.5) returns 0
    %
    % See also ERFINV

    % Validate input
    if any(x < 0 | x > 1)
        error('Input x must be in the interval [0, 1]');
    end

    % Adjust extreme values to avoid Inf or -Inf in output
    x(x == 0) = 0.0001;
    x(x == 1) = 0.9999;

    % Compute the probit function
    phi = sqrt(2) * erfinv(2*x - 1);
end