function hodghux

%the model used
% alpha_n=.01*(v+55)/(1-exp(-.1*(v+55)));
% beta_n=.125*exp(-.0125*(v+65));
% 
% alpha_m=.1*(v+40)/(1-exp(-.1*(v+40)));
% beta_m=4*exp(-.0556(v+65));
% 
% alpha_h=.07*exp(-.05*(v+65));
% beta_h=1/(1+exp(-.1*(v+35)));
%--------------------------------

close all
tic

%options -- set carefully  -- consider odephas2 as output
options = odeset('RelTol',1e-4,'AbsTol',[1e-4]);%,'stats','on','outputfcn','odeplot');


i_electrode=[0:10:100];
f_i_curve=[];
firstspike_times=[];
for iii=1:length(i_electrode)
    y0(1)=-65;%voltage
	y0(2)=.32;%n
	y0(3)=.05;%m
	y0(4)=.596;%h
	y0(5)=5;%tau a
	y0(6)=5;%tau b
	y0(7)=0.5557;%a
	y0(8)=0.228;%b
 
    [t,v] = ode15s(@connor_stevens_solve,[0 300],y0,options,i_electrode(iii));
    
%     figure(iii)
%     plot(t,v(:,1),'-o'), hold on
    
    v=v(:,1);
    thresh_=0;%mean(v);
	above_threshold=v>thresh_ ;
    if isempty(above_threshold),error('s'),end
	beginnings_iftrue=~above_threshold(1:end-1)&above_threshold(2:end);
	endings_iftrue=above_threshold(1:end-1)&~above_threshold(2:end);
	AP_start_indices=1+find(beginnings_iftrue);
    AP_end_indices=find(endings_iftrue);
    %if ~isempty(AP_end_indices) & ~isempty(AP_start_indices) & AP_end_indices(1) < AP_start_indices(1), AP_end_indices(1)=[];end
    if ~isempty(AP_start_indices) & AP_start_indices(end)>AP_end_indices(end),AP_start_indices(end)=[];end
	
	%find the max value within each group of suprathreshold values
	peak_values=[];
	peak_times=[];
	for ttt=1:length(AP_start_indices)
        [y,ind]=max(v(AP_start_indices(ttt):AP_end_indices(ttt)));
        ind=ind(1)+AP_start_indices(ttt)-1;
        peak_times=[peak_times,ind];
        peak_values=[peak_values,y(1)];
	end
	%remove any stim artifacts within stim_artifact_jitter of a previous stim artifact
	% z=length(AP_peaks);
	% if length(AP_peaks)>1
	% 	for bleb=2:length(AP_peaks)
	%         if AP_peaks(z)-AP_peaks(z-1)<stim_artifact_jitter
	%             AP_peaks(z)=[];
	%         end
	%         z=z-1;
	%     end
	% end
	%plot(t(peak_times),v(peak_times,1),'ms','markersize',5,'linewidth',3)
	c=t(peak_times);
	a=diff(c);
	b=mean(a);
    f_i_curve=[f_i_curve,b];
    if ~isempty(peak_times),first_time=t(peak_times(1));else first_time=0;end
    firstspike_times=[firstspike_times,first_time];
end
figure(50)
f_i_curve(isnan(f_i_curve))=0;
plot(i_electrode,f_i_curve,'ks','markersize',5,'linewidth',3)
title('F-I Relation'),xlabel('I_e_l_e_c_t_r_o_d_e'),ylabel('Firing rate')
figure(51)
plot(i_electrode(find(firstspike_times)),firstspike_times(find(firstspike_times)),'bo','markersize',5,'linewidth',3)
title('First Spike Times'),xlabel('I_e_l_e_c_t_r_o_d_e'),ylabel('Time')
toc