% plotItotVclamp.m
% Jessica Parker, October 18, 2024
%
% This Matlab script plots the total current of the data simulated in this directory only, so it does not do leak 
% substraction and therefore it also does not measure INaP or gNaP. This is mainly helpful for looking at the total 
% current in the leak measurement simulations where the magnitude of the voltage steps have been quartered.

clear all;
close all;

vrsn = 'A';

run1 = 5;
run2 = 0;
run3 = 0;
run4 = 2;
run5i = 1;
run5f = 3;
run6i = 1;
run6f = 3;
skpPltTrc = 1; % Number of run5 instances (+1) to skip between traces plotted. To skip none, set to 1.

dir1 = [num2str(run1) '_' num2str(run2) '_' num2str(run3)];
dir2 = [dir1 '_' num2str(run4)];

dataDir = 'data/';
prfx = 'NaiMnap'; % Prefix used in the name of the simulation data file
nvar = 2; % Number of variables in simulation data file
Naindx = 1; % Index of [Na+]i variable in simulation data file
mNaPindx = 2; % Index of mNaP variable in simulation data file
tint = 0.0001; % Time step used by integrator

fntnm = 'Lucida Sans';
fntsz = 12;
vsblty = 'on'; % Set to 'off' to make figures not pop up or set to 'on'
CloseFigs = 0; % Set to 1 to automatically close figures after printing them to a file, or set to 0 to leave them open.
ZoomInOnStepUp = 1; % Set to 1 to zoom in on the depolarization step or set to zero and set xmin and xmax below to the axis range you want
xmin = 75;
xmax = 100;

readpars; % Model parameters read from readpars.m, not all of these are actually used by this code

yyStrt = [];

for run5 = run5i:skpPltTrc:run5f
    Ipmpmx = 40.0*run5;

    yy = yyStrt;
    for run6=run6i:run6f   % Concatenating portions of time saved in separate files

        dir3 = [dir2 '_' num2str(run5)];

        VarN = [dataDir prfx dir3 '_' num2str(run6) '.dat'];
        fl = fopen(VarN);
        yy0 = fread(fl,[nvar, 10000000],'double')'; % Loading simulation data
        fclose(fl);

        yy = [yy(1:end-1,:); yy0]; % Remove last point here, because otherwise this time point will be repeating as the first time point of the next file
        indxs(run6) = length(yy(:,1));
        if run6 < run6f
            lngths(run6) = length(yy0(:,1))-1;
        else
            lngths(run6) = length(yy0(:,1));
        end
    end

    Vsteps = [-90, -80, -40];
    Vm = [];
    for aa = 1:length(Vsteps)
        Vm = [Vm; Vsteps(aa)*ones(lngths(aa),1)];
    end
    
    tt = 0.0:tint:tint*(indxs(end)-1);

    Inap = gNaP*yy(:,mNaPindx).*(Vm-rtf*log(Nae./yy(:,Naindx))); % INaP 
    IleakNa = gL*(EL-EK)/(ENa-EK)*(Vm-rtf*log(Nae./yy(:,Naindx))); % ILeakNa 
    IleakK = gL*(EL-ENa)/(EK-ENa)*(Vm-EK); % ILeakK 
    Ileak = IleakNa + IleakK;
    Itot = Inap + Ileak;

    if ZoomInOnStepUp
        xmin = tt(indxs(end-1)+1)-0.01;
        xmax = tt(indxs(end-1)+1) + 0.1;
        mrk = 'zm';
    else
        mrk = '';
    end

    f = figure('visible',vsblty);
    f.PaperPositionMode = 'manual';
    f.PaperUnits = 'inches';
    f.Units = 'inches';
    f.OuterPosition = [1 1 7.0 4.0];
    f.InnerPosition = [0.25 0.25 6.5 3.5];
    f.PaperPosition = [0.25 0.25 6.0 3.0];
    f.RendererMode = 'manual';

    axes('position',[0.15 0.6 0.8 0.35]);
    hold on;
    plot(tt, Itot, 'k', 'linewidth', 1.5);
    ylabel('I_{Tot} (pA)');
    xticks([]);
    xlim([xmin xmax]);
    box off;
    ax = gca;
    ax.FontName = fntnm;
    ax.FontSize = fntsz;
    ax.LineWidth = 2.0;

    axes('position',[0.15 0.2 0.8 0.35]);
    hold on;
    plot(tt, Vm, 'k', 'linewidth', 1.5);
    ylabel('V_{clamp} (mV)');
    xlabel('Time (s)');
    xlim([xmin xmax]);
    box off;
    ax = gca;
    ax.FontName = fntnm;
    ax.FontSize = fntsz;
    ax.LineWidth = 2.0;

    print(f,['plots/ItotVclamp' dir2 '_' num2str(run5) vrsn mrk '.eps'],'-depsc','-r0');

    if strcmp(vsblty,'off')
        close(f);
    elseif CloseFigs
        close(f);
    end
end