function reversal_potential_sim
% REVERSAL_POTENTIAL_SIM  Passive single-neuron simulator with Na, K, Cl, Ca
% conductances. Erev per ion is computed from inside/outside concentrations
% via the Nernst equation. No spikes; passive RC only. Press Update to
% recompute and animate the Vm trace.
%
% Layout: Plot on top. Bottom bar = 4 ion columns (Na,K,Cl,Ca) + a compact
% control column (V0, Speed, Update). The entire UI is inside a scrollable
% panel; if the figure gets too small, OS-style scrollbars appear for the
% whole app.

%% Figure with OS-style scrollbars for the WHOLE app
fig = uifigure('Name','Reversal Potential Simulator','Color','w');
fig.Position(3:4) = [700 420];  % smaller default size (compact)  % smaller default size (compact)

% Larger natural content canvas ensures scrollbars when window is smaller

% Scrollable container hosting all content
mainPanel = uipanel(fig,'Scrollable','on');
mainPanel.Units = 'normalized';
mainPanel.Position = [0 0 1 1];

% Larger content canvas ensures scrollbars appear when figure is smaller
content = uipanel(mainPanel,'Units','pixels','Position',[0 0 900 600]);

% Root grid inside the content canvas
root = uigridlayout(content,[2 1], ...
    'RowHeight',{'1x','fit'}, ...
    'ColumnWidth',{'1x'}, ...
    'Padding',10,'RowSpacing',8,'ColumnSpacing',8);

%% TOP: plot
ax = uiaxes(root); ax.Layout.Row = 1; ax.Layout.Column = 1; ax.Box = 'on';
ax.XLabel.String = 'Time (ms)'; ax.YLabel.String = 'V_m (mV)';
ax.Title.String = 'Membrane potential (passive)';

%% BOTTOM: ion slider row (4 columns) + control column
ions = {"Na","K","Cl","Ca"};
colors = [0.80 0.10 0.10; 0.10 0.60 0.10; 0.10 0.40 0.80; 0.55 0.10 0.55];
valence = struct('Na',1,'K',1,'Cl',-1,'Ca',2);
% Default concentrations (mM)
defCout = struct('Na',145,'K',5,'Cl',110,'Ca',2);
defCin  = struct('Na',15,'K',140,'Cl',10,'Ca',0.0001);

bottom = uigridlayout(root,[1 5], 'ColumnWidth',{'1x','1x','1x','1x','fit'}, 'RowHeight',{'fit'});
bottom.Layout.Row = 2; bottom.Layout.Column = 1;

% UI handle containers
S.g = struct(); S.Cout = struct(); S.Cin = struct();

for i = 1:numel(ions)
    ion = ions{i}; c = colors(i,:);
    col = uigridlayout(bottom,[4 1], 'RowHeight',{'fit','fit','fit','fit'}, 'ColumnWidth',{'1x'});
    col.Layout.Row = 1; col.Layout.Column = i;

    % Large, centered ion label
    uilabel(col,'Text',ion,'FontWeight','bold','FontSize',16,'HorizontalAlignment','center', ...
        'FontColor',c);

    % g slider
    gRow = uigridlayout(col,[1 2],'ColumnWidth',{'fit','1x'});
    uilabel(gRow,'Text','g');
    S.g.(ion) = uislider(gRow,'Limits',[0 2],'Value',1,'MajorTicks',[0 1 2], ...
        'MajorTickLabels',{'0','Norm','High'});

    % [out] slider + live value
    outRow = uigridlayout(col,[1 3],'ColumnWidth',{'fit','1x','fit'});
    uilabel(outRow,'Text','[out]');
    S.Cout.(ion) = uislider(outRow);
    lblOut = uilabel(outRow,'Text','');

    % [in] slider + live value
    inRow = uigridlayout(col,[1 3],'ColumnWidth',{'fit','1x','fit'});
    uilabel(inRow,'Text','[in]');
    S.Cin.(ion) = uislider(inRow);
    lblIn = uilabel(inRow,'Text',''); %#ok<NASGU>

    % Limits and defaults
    switch ion
        case {"Na","K","Cl"}
            set(S.Cout.(ion),'Limits',[0 200],'Value',defCout.(ion));
            set(S.Cin.(ion), 'Limits',[0 200],'Value',defCin.(ion));
        case "Ca"
            set(S.Cout.(ion),'Limits',[0 10],'Value',defCout.(ion));
            set(S.Cin.(ion), 'Limits',[0 5] ,'Value',defCin.(ion));
    end

    % Numeric labels update live while dragging + on release
    lblOut.Text = sprintf('%.4g', S.Cout.(ion).Value);
    lblIn.Text  = sprintf('%.4g', S.Cin.(ion).Value);
    S.Cout.(ion).ValueChangingFcn = @(s,e) set(lblOut,'Text',sprintf('%.4g',e.Value));
    S.Cout.(ion).ValueChangedFcn  = @(~,~) set(lblOut,'Text',sprintf('%.4g',S.Cout.(ion).Value));
    S.Cin.(ion).ValueChangingFcn  = @(s,e) set(lblIn,'Text', sprintf('%.4g',e.Value));
    S.Cin.(ion).ValueChangedFcn   = @(~,~) set(lblIn,'Text', sprintf('%.4g',S.Cin.(ion).Value));
end

% Rightmost column: V0, Speed, Update
ctrl = uigridlayout(bottom,[3 2],'ColumnWidth',{'fit','1x'}, 'RowHeight',{'fit','fit','fit'});
ctrl.Layout.Row = 1; ctrl.Layout.Column = 5;

uilabel(ctrl,'Text','V0 (mV)'); V0Edit = uieditfield(ctrl,'numeric','Limits',[-120 80],'Value',-65);
uilabel(ctrl,'Text','Speed');   speedEdit = uieditfield(ctrl,'numeric','Limits',[0.1 10],'Value',1);

btn = uibutton(ctrl,'Text','Update','ButtonPushedFcn',@(~,~) runSim());
btn.Layout.Row = 3; btn.Layout.Column = [1 2];

%% Core: simulation + animation
    function runSim()
        % Constants and fixed sim params (short plot height -> keep traces tight)
        Tdeg = 37;           % deg C (fixed)
        C_m  = 1;            % uF/cm^2 (fixed)
        Tms  = 1000;         % ms (fixed duration)
        dt   = 0.1;          % ms (fixed step)
        V0   = V0Edit.Value; % initial voltage
        spd  = speedEdit.Value; % playback speed

        % Conductances
        g.Na = S.g.Na.Value; g.K = S.g.K.Value; g.Cl = S.g.Cl.Value; g.Ca = S.g.Ca.Value;

        % Erev via Nernst (mV)
        E.Na = nernst_mV(S.Cout.Na.Value, S.Cin.Na.Value,  1, Tdeg);
        E.K  = nernst_mV(S.Cout.K.Value,  S.Cin.K.Value,   1, Tdeg);
        E.Cl = nernst_mV(S.Cout.Cl.Value, S.Cin.Cl.Value, -1, Tdeg);
        E.Ca = nernst_mV(S.Cout.Ca.Value, S.Cin.Ca.Value,  2, Tdeg);

        % Time vector and analytic solution
        t = 0:dt:Tms;
        gtot = g.Na + g.K + g.Cl + g.Ca;
        if gtot <= 0
            V = V0 + zeros(size(t));
        else
            E_inf = (g.Na*E.Na + g.K*E.K + g.Cl*E.Cl + g.Ca*E.Ca)/gtot;
            tau   = C_m/gtot; % ms
            V     = E_inf + (V0 - E_inf)*exp(-t./tau);
        end

        % Animate Vm
        cla(ax);
        h = animatedline(ax,'LineWidth',2);
        ax.XLim = [0 Tms];
        % Tighter y-lims (shorter plot height), with margin
        ylo = min(min(V), min([E.Na E.K E.Cl E.Ca])) - 6;
        yhi = max(max(V), max([E.Na E.K E.Cl E.Ca])) + 6;
        if ~isfinite(ylo) || ~isfinite(yhi)
            ylo = min(V0, -80)-6; yhi = max(V0, 60)+6;
        end
        ax.YLim = [ylo yhi];

        % Dotted Erev HORIZONTAL lines in ion colors (y = Erev)
        ev = [E.Na, E.K, E.Cl, E.Ca];
        for ii = 1:4
            yline(ax, ev(ii), ':', 'LineWidth',1.5, 'Color', colors(ii,:));
        end

        % Playback pacing (decimated for responsiveness)
        maxFrames = 400; frameStep = max(1, ceil(numel(t)/maxFrames));
        tic;
        for k = 1:frameStep:numel(t)
            addpoints(h, t(k), V(k));
            drawnow limitrate nocallbacks;
            if spd > 0
                target_dt = (dt*frameStep) / spd / 1000;  % seconds per frame
                while toc < target_dt, pause(0.001); end
                tic;
            end
        end
    end

end % reversal_potential_sim

%% Helpers
function E = nernst_mV(Cout,Cin,z,Tdeg)
% Nernst potential in mV; concentrations in mM, temp in C
R = 8.314462618; F = 96485.33212; T = Tdeg + 273.15;
E = (R*T/(z*F)) * log(Cout/Cin) * 1e3; % mV
end