function plot_alpha_beta

% Copyright 2007 John L Baker. All rights reserved.
% This software is provided AS IS under the terms of the Open Source
% MIT License. See http://www.opensource.org/licenses/mit-license.php.

% Plot an alpha beta table generated by VoltageDepTabChannel etc.

% Pick the directory where data results files are found
DATADIR='C:\BNSF_1_0\Data\Default\'; % adjust as needed

abfile='test-rallpack-h.txt'; % file to plot
disp(['Plotting alpha-beta values for ',abfile]);

ab=dlmread([DATADIR 'test-rallpack-m.txt']); % read alpha-beta values

n=1; % assumed exponent value in a Hodgkin-Huxley gating arrangement
nlabel=n; % exponent value to include in plot title

vmin=-80; % minimum voltage to plot
vmax= 20; % maximum voltage to plot

v=ab(:,1)*1e3;
alpha = ab(:,2);
beta = ab(:,3);
tau=1./(alpha+beta);
xinf=alpha.*tau;
xinfn = xinf.^n;

vhidx=find(sign(xinfn(1:end-1)-.5) ~= sign(xinfn(2:end)-.5));
if ~isempty(vhidx)
   vhidx=vhidx(1);
	vhalf=v(vhidx);
   k=(v(vhidx+1)-v(vhidx-1))./(xinfn(vhidx+1)-xinfn(vhidx-1))/4;
else
   vhalf=0;
   k=1;
end

if 1
   % Extract more parameters from the alpha-beta relationship (optional)
   maxTau = 1000*max(tau)
   minTau = 1000*min(tau(find(v>=vmin)))
   maxTauIdx = min(find(maxTau/1000==tau));
   vmaxtau = v(maxTauIdx)
   gamma=1/(exp(-(vmaxtau-vhalf)/k)+exp((vmaxtau-vhalf)/k) )
end

% Show estimated fitting parameter valued (not that accurate)
vhalf
k

if 1
   % Plot all table entries (optional)
   figure
	subplot(2,2,1)
	plot(v,alpha);
	title('alpha');
	subplot(2,2,2)
	plot(v,beta);
	title('beta')
	subplot(2,2,3)
	plot(v,xinfn)
	title('xinf^n')
	subplot(2,2,4)
	plot(v,tau*1000)
	title('tau')
end

if 0
   % Plot just stable state vs voltage (optional)
   figure
   vidx = find(v<=vmax & v>=vmin);
   plot(v(vidx),xinfn(vidx));
   title(['Xinf ^' num2str(n)])
   %grid on
   
   % Plot fit for comparison
   if 1
      hold on
      plot(v(vidx),1./(1+exp(-(v(vidx)-vhalf)/k)),'--')
   end
end

if 0
   % Plot time constant (optional)
	figure
	vidx = find(v<=vmax & v>=vmin);
	plot(v(vidx),1000*tau(vidx));
	title('Tau (msec)')
   %grid on
   
   % Plot fit for comparison
   if 1  
      z=1; % charge

      hold on  
      
      plot(v(vidx),minTau+(maxTau-minTau)* ...
         (exp(z*gamma*(vmaxtau-vhalf)/k)+exp(-z*(1-gamma)*(vmaxtau-vhalf)/k)) ./ ...
         (exp(z*gamma*(v(vidx)-vhalf)/k)+exp(-z*(1-gamma)*(v(vidx)-vhalf)/k)) ...
         ,'--')
   end
end

if 1
   % Plot Xinf^n and time constant together
   figure
   vidx = find(v<=vmax & v>=vmin);
   [ax,h1,h2]=plotyy(v(vidx),xinfn(vidx),v(vidx),1000*tau(vidx));
   xlabel('Vm (mV)')
   set(ax(1),'TickDir','out')
   set(ax(2),'TickDir','out')
   axes(ax(1))
   ylabel(['Xinf ^' num2str(nlabel)])
   axes(ax(2))
   ylabel('Tau (ms)')
   if max(tau(vidx))-min(tau(vidx))<1e-4
      axis([vmin vmax floor(min(1000*tau(vidx))) ceil(max(1000*tau(vidx)))])
   end
end