function [V2,A2] = HH1step(V1,N,i,A1)
% function [V2,A2] = HH1step(V1,N,i,A1)
% 
% INPUT
% V1 is the initial membrane potential
% N is the noise
% i is the injected current
% A1 is the vector of initial activation states 
%
% OUTPUT
% V2 is the final membrane potential
% A2 is the vector of final activation states
%
% ACTIVATION STATES
%   [A1(1), A1(2), A1(3)] = (Na, K activation, Na inactivation) or (m,n,h)
%   m,n,h are between 0 and 1
%   fraction of open sodium channels is m^3 h
%   fraction of open potassium channels is n^4

%%% DEFAULT VALUES FOR INPUTS %%%

% If no noise is passed, set noise to zero
if nargin < 2
  q = [0 0 0 0];
end

% if injected current not set, set to zero
if nargin < 3
  i = 0;
end

% If activation variables not set, set them to resting state
if nargin < 4
  A1(1) = 0.05293421762087;
  A1(2) = 0.31768116757978;
  A1(3) = 0.59611104634676;
end

%%% DEFAULT VALUES FOR PARMETERS %%%

% Maximal conductances
gNa = 120;   % mS/cm^2
gK  = 36;    % mS/cm^2
gL  = 0.3;   % mS/cm^2

% Reversal Potentials
ENa = 115;   % mV
EK  = -12;   % mV
EL  = 10.6;  % mV

% Membrane Capacitance
C = 1;       % uF/cm^2

% Time step
dt = 0.01;    % msec

%%% STATE VARIABLES

V = V1;
m = A1(1);
n = A1(2);
h = A1(3);

%%% COMPUTATIONS %%%

% Compute activation factors
alpha_m = (2.5 - 0.1.*V)./(exp(2.5 - 0.1.*V) - 1);
alpha_n = (0.1 - 0.01.*V)./(exp(1 - 0.1.*V) - 1);
alpha_h = 0.07.*exp(-V./20);
beta_m = 4.*exp(-V./18);
beta_n = 0.125.*exp(-V./80);
beta_h = 1./(exp(3-0.1.*V)+1);

% Compute ionic current
Ik = gNa.*m.^3.*h.*(V-ENa) + gK.*n^4.*(V-EK) + gL.*(V-EL);

% Scale noise for duration of time step
noise = sqrt(dt)*N; 

% Compute new membrane potential
V2 = V+dt./C.*(-Ik + i) + noise;

% Compute new activation variables
A2 = [m+dt.*(alpha_m.*(1-m) - beta_m.*m);
       n+dt.*(alpha_n.*(1-n) - beta_n.*n);
       h+dt.*(alpha_h.*(1-h) - beta_h.*h)]; 
