% plotPumpPaperFig1A1.m
% Jessica Parker, October 18, 2024
%
% Makes Figure 1A1 and Figure 1A2. Plots Vm, [Na+]i, Ipump, and INaPn over time.

clear all;
close all;

mrk = 'sc';

run1 = 1;
run2 = 0;
run3 = 0;
run4i = 0;
run4f = 0;
run5i = 0;
run5f = 0;
run6i = 1;
run6f = 1;
datadir = 'data/'; % Directory where input data files are located
prfx = 'Vall'; % Prefix of input data file name
nvar = 5; % Number of variables saved in the input data file
vsblty = 'on';
CloseFigs = 0;
SingleBurst = 0; % Set to 1 to zoom in on a single burst, set to 0 to use xmin0 and xmax0 below to define x-axis range. 

% Parameters
readpars;

nFigRows = 4;
titles = {'V_m','[Na^+]_i','I_{Pump}','I_{NaPn}'};
xmin0 = 5.5; % Set the time range that you want to plot
xmax0 = 21.5;
ymins = [-80, 18, 5, -130];
ymaxs = [12, 30, 28, 0];
scalrefs = [-80, 20, 5, 0];
scalmins = [-40, 25, 15, -120];
scalmaxs = [0, 30, 25, -60];
pbp = 0.05; % If you set SingleBurst to 1, this is the proportion of burst period to plot on either side of the burst duration.
tint = 0.0001; %Time step between data points
fntsz = 16; %Font size
fntnm = 'Lucida Sans';
fntwght = 'bold';
rowUnits = {'mV','mM','pA','pA'};

% Set widths and heights of subplots and the spaces between and around them
pgBttm = 0.45; % Space below bottom figure in inches
pgTop = 0.2; % Space above top figure in inches, the title will be placed within this space if you have a title
figSep = 0.4;
figHght = 0.7; % Height of each subplot in inches
scl = 1; % Scaling factor if you want to make everything bigger or smaller without changing any proportions.

if contains(prfx,'Vall')
  nNai = 4;
  mNaPi = 5;
  hNaPi = 6;
  hNaFi = 3;
  mKi = 2;
elseif contains(prfx,'VNa')
  nNai = 2;
end


aa = 1;
for run4=run4i:run4f
  for run5=run5i:run5f

    yy = [];
    for run6=run6i:run6f   %Concatenating portions of time saved in separate files
      dir2 = [num2str(run1) '_' num2str(run2) '_' num2str(run3)];
      dir3 = [dir2 '_' num2str(run4) '_' num2str(run5)];

      VarN=[datadir prfx dir3 '_' num2str(run6) '.dat'];
      f1=fopen(VarN);
      yy1=fread(f1,[nvar, 10000000],'double')';
      fclose(f1);

      yy = [yy; yy1(1:end-1,:)];  %Remove last point here, because otherwise you will repeat the same time point
    end                           %twice between consecutive run6 files   
				
    
    lnt = length(yy(:,1));
    tt = 0:tint:tint*(lnt-1);

    plotDat(1,:) = yy(:,1);
    plotDat(2,:) = yy(:,nNai);
    plotDat(3,:) = Ipmpmx./(1.0+(Naih*exp(alpha*yy(:,1)/rtf)./yy(:,nNai)).^nhillNa);
    plotDat(4,:) = gNaP*yy(:,mNaPi).*(yy(:,1)-rtf*log(Nae./yy(:,nNai)));

    bc = [];
    if SingleBurst
      if exist([datadir 'bd1j' dir3 '_' num2str(run6f) mrk '.txt'])
        bc = load([datadir 'bd1j' dir3 '_' num2str(run6f) mrk '.txt']);
        xmin = bc(1,2)-pbp*bc(2,2);
        xmax = bc(1,2)+bc(3,2)+pbp*bc(2,2);
      elseif exist([datadir 'plts' dir3 '_' num2str(run6f) mrk '.txt'])
        bc = load([datadir 'plts' dir3 '_' num2str(run6f) mrk '.txt']);
        xmin = bc(1,2)-pbp*bc(2,2);
        xmax = bc(1,2)+bc(3,2)+pbp*bc(2,2);
      else
        xmin = xmin0;
        xmax = xmax0;
      end
    else
      xmin = xmin0;
      xmax = xmax0;
    end
    xl = xmax-xmin;
    xmini = round(xmin/tint)+1;
    xmaxi = round(xmax/tint)+1;


    pgHght = nFigRows*figHght + pgBttm + pgTop + (nFigRows-1)*figSep;

    f = figure('visible',vsblty);
    f.PaperPositionMode = 'manual';
    f.PaperUnits = 'inches';
    f.Units = 'inches';
    f.OuterPosition = [1.0 1.0 7.5 pgHght+1.0];
    f.InnerPosition = [0.25 0.25 7 pgHght+0.5];
    f.PaperPosition = [0.25 0.25 6.5 pgHght];
    f.RendererMode = 'manual';

    pgBttm = scl*pgBttm;
    pgTop = scl*pgTop;
    figSep = scl*figSep;
    figHght = scl*figHght;

    topPerc = 1 - pgTop/pgHght;
    bttmPerc = pgBttm/pgHght;
    figPerc = figHght/pgHght;
    restPerc = topPerc - bttmPerc;
    sepPerc = figSep/pgHght;

    for bb = 1:nFigRows
      ymin = ymins(bb);
      ymax = ymaxs(bb);
      yl = ymax - ymin;
      minmaxs(bb,1) = min(plotDat(bb,xmini:xmaxi));
      minmaxs(bb,2) = max(plotDat(bb,xmini:xmaxi));

      axes('position',[0.04 topPerc+sepPerc-(figPerc+sepPerc)*bb 0.89 figPerc]);
      plot(tt(xmini:xmaxi),plotDat(bb,xmini:xmaxi),'k-','linewidth',1.5);
      hold on;
      xlim([xmin-0.06*xl xmax+0.17*xl]);
      ylim([ymin ymax]);
      plot([xmax-0.03*xl xmax-0.03*xl]+0.05*xl,[scalmins(bb) scalmaxs(bb)],'k','linewidth',3.5);
      scalLbl = [num2str(scalmaxs(bb)-scalmins(bb)) ' ' rowUnits{bb}];
      text(xmax-0.01*xl+0.05*xl,0.5*(scalmins(bb)+scalmaxs(bb)),scalLbl,'fontname',fntnm,'fontsize',fntsz,'fontweight',fntwght);
      plot([xmax-0.03*xl xmax]+0.05*xl,scalrefs(bb)*[1 1],'k','linewidth',3.5);
      scalRefLbl = [num2str(scalrefs(bb)) ' ' rowUnits{bb}];
      text(xmax+0.02*xl+0.05*xl,scalrefs(bb),scalRefLbl,'fontname',fntnm,'fontsize',fntsz,'fontweight',fntwght);
      text(xmin-0.07*xl,ymin+0.7*yl,titles{bb},'fontname',fntnm,'fontsize',fntsz+1,'fontweight',fntwght);
      axis off;
    end

    axes('position',[0.04 0 0.89 bttmPerc+0.05]);
    if SingleBurst
      plot([xmax-0.73 xmax-0.23],[0.45 0.45],'k','linewidth',3.5);
      text(xmax-0.56,0.23,'500 ms','fontname',fntnm,'fontsize',fntsz,'fontweight',fntwght);
    else
      plot([xmax-6.3 xmax-2.3],[0.45 0.45],'k','linewidth',3.5);
      text(xmax-4.8,0.23,'4 s','fontname',fntnm,'fontsize',fntsz,'fontweight',fntwght);
    end
    xlim([xmin-0.06*xl xmax+0.17*xl]);
    ylim([0 1]);
    axis off;

    if SingleBurst && length(bc)>0
      print(f,['plots/sbVNaiIpumpInapFrml' dir3 '_' num2str(run6) '.eps'],'-depsc','-r0');
    else
      print(f,['plots/VNaiIpumpInapFrml' dir3 '_' num2str(run6) '.eps'],'-depsc','-r0');
    end
    
    if strcmp(vsblty,'off')
      close(f);
    elseif CloseFigs
      close(f);
    end
    aa = aa + 1;
  end
end