%2/25/02
%Point is to integrate the HR equations with some constant current 

%Equations are from Rose and Hindmarsh (1989), Proc Royal Soc Lond B 237, 267-388
%Eqns. 9, p.271

%variables: q,v




%***************PARAMETERS************
global I C gbarNa ENa gbarK EK gbarL EL gA A Tb Tn

%I=4.9456;        %baseline; stim gets added in drift v below
I=4.95 ; %uA, microamp       %baseline; stim gets added in drift v below
C=1 ;  %uF/cm^2, microFarad/cm^2
gbarNa=120 ;  %(mS)/cm^2  mSimens/cm^2
ENa=55 ; %mV
gbarK=20 ;  %(mS)/cm^2  mSimens/cm^2
EK=-72 ;  %mV
gbarL=.3 ;  %(mS)/cm^2  mSimens/cm^2
EL=-17 ;  %mV
gA=47.7 ; %(mS)/cm^2  mSimens/cm^2

A=.21*gA/gbarK;
Tn=.52 ;
Tb=1 ;




%LOOP PARAMETERS
deltat=.01 %.005;

%num_iters=ceil(20500/deltat)
num_iters=ceil(3000/deltat)

%time is in ms

%%%%
timelist=linspace(0,num_iters*deltat,num_iters);
%*********INITIAL CONDITIONS**********
qlist=zeros(1,num_iters) ;         
vlist=zeros(1,num_iters) ;


q_ic=.106  ;      
v_ic=50.81  ;

qlist(1)=q_ic ;
vlist(1)=v_ic ; 

q=qlist(1) ;
v=vlist(1) ;

for iter=1:num_iters-1
    
    if mod(iter,round(num_iters/100))==0
        fracdone=iter/num_iters
    end
 
    stim=0 ;
    
    q=qlist(iter) ;
    v=vlist(iter) ;
    
       
           %rk O(4) method calculate newx1 and newx2 for the next timestep
           %---indexing on F's:  Fij, i=rk term #, j=var. #
           
           F1q=driftqHR(q,v) ;
           F1v=driftvHR(q,v,stim);
                             
           F2q=driftqHR(q+(1/2)*F1q*deltat, v+(1/2)*F1v*deltat);
           F2v=driftvHR(q+(1/2)*F1q*deltat, v+(1/2)*F1v*deltat, stim); 
           
           F3q=driftqHR(q+(1/2)*F2q*deltat, v+(1/2)*F2v*deltat);
           F3v=driftvHR(q+(1/2)*F2q*deltat, v+(1/2)*F2v*deltat,stim);
           
           F4q=driftqHR(q+F3q*deltat,v+F3v*deltat);
           F4v=driftvHR(q+F3q*deltat,v+F3v*deltat,stim);
           
           qlist(iter+1)=q +(1/6)*(F1q+2*(F2q+F3q)+F4q)*deltat ; 
           vlist(iter+1)=v +(1/6)*(F1v+2*(F2v+F3v)+F4v)*deltat ;
           
end 

figure
plot(qlist,vlist,'-')


figure
plot(timelist,vlist)
