function [Rtotal, Rseries, Rmem, Cmem, Ileak, current] = testpulse(Fs, scaling, gain, amplitude, pulsewidth, nrepeat)
%
% [Rtotal, Rseries, Rmem, Cmem, Ileak, current] = 
%            testpulse(Fs, scaling, gain, amplitude, pulsewidth, nrepeat)
%                  testpulse -- meant to be called from axontp
%
% Returned values:
% All resistances are returned in Ohms. Cmem is in Farads. Ileak
% and current are in pA.
%
% If nrepeat is > 1, then all returned values are averaged over 
% nrepeat trials.
%
% Revision history:
% Implemented Rm, Cm calculation from Axon diagram 5/2003 (Neville), 
% Cleaned up, added 2nd channel 11/2003 (Naveen), 
% Fixed Ra mislabeling, fixed bugs with nrepeat 2/2004 (Neville)

	% The STIM voltage array will be what we send through vstim.
	% We make it here.  We will send the same pulse to both the
	% output channels.
	
	pulse=amplitude*[zeros(pulsewidth,1); -ones(pulsewidth,1); zeros(pulsewidth,1)];
	stim=[pulse pulse];
	
	dV = amplitude * 10^-3;
	
	avgCurrent = zeros(length(pulse),2);
	Rtotal=zeros(1,2);
	Rseries=zeros(1,2);
	Rmem=zeros(1,2);
	Cmem=zeros(1,2);
	Ileak=zeros(1,2);
 	
	% In this loop, we send out the voltage pulse, and calculate
	% the different parameters based on the output current.
	
	for irepeat = 1:nrepeat

		current=vstim(stim, scaling, gain, Fs, 1);  % current in pA
		
		% The bulk of the calculations below are based on figures and
		% info in Axon Knowledge base article #652, found at: 
		% http://www.axon.com/mr_Axon_KB_Article.cfm?ArticleID=652
		
		% As in Axon, I1 is steady state current during pulse and I2 is
		% the steady state current before pulse
		I1=mean(current(350:390,:));
		I2=mean(current(80:160,:));
		
		% We will calculate four variables: rt, rm, rs, and cm, for each
		% of the two channels and then assign them to our returning variables
		
		for ichannel = 1:2

			Ileak(ichannel) = Ileak(ichannel) + mean(current(80:160,ichannel));
			
			% We want to prevent "divide-by-zero" kind of errors, so we make sure
			% the denominator isn't zero.
			denom = mean(current(80:160,ichannel)) - mean(current(280:360, ichannel));
			if denom == 0 
				rt = 0;
			else
				rt = amplitude / denom;
			end

			% Because amplitude is in mV, current is in pA, this returns
			% Rtotal in gOhms.  We will convert it to Ohms below.
			
			% p will have a polynomial fit to the first little bit of the
			% reponse to the downward step (which started at element 201
			% of the current array)
			
			p=polyfit(2:4,current(202:204,ichannel)',1);  % linear extrapolation
			p=polyfit(3:8,current(203:208,ichannel)',2);  % quadratic extrapolation
			
			% Convert measured Rt to Ohms
			rt = rt * 10 ^9;
			
			% Q1, Q2 and Qt are from the Axon diagram for Membrane Test
			
			% Trying to find all points on the capacitive transient to
			% calculate Q1
			%
			% Not sure if this is really the principled way to do this??
			% We are trying to get all points where current > the steady
			% steady state current during the pulse
			% ASSUMPTION: Capacitive transient lasts until at least 202
			%             (true for all cells)
			Q1_pts1 = find(current(202:400,ichannel) >= I1(ichannel));
			numPointsAfterQ1 = length(Q1_pts1);
			
			%  By taking the second point, we allow some error... this ensures
			% that we don't "stop including" points in Q1 too early.
			% (e.g. some digitization mistake, etc.)
			
			bestPointToEndTransient = 2;
			
			% If there's no transient, this IF statement prevents errors
			% (e.g. with a tight, tight seal in cell-attached mode)
			if (numPointsAfterQ1 == 0)
				Qt = 0;
			else
				if (numPointsAfterQ1 > bestPointToEndTransient)
					Q1_lastpoint1 = Q1_pts1(bestPointToEndTransient);
				else
					Q1_lastpoint1 = Q1_pts1(numPointsAfterQ1);
				end

				% Calculate the area under the capacitive transient ABOVE I1
				% This area is Q1. We are doing numerical integration here
				% OLD NUMERICAL INTEGRATION (Rectangle method)
				% Q1 = sum(abs(current((201:201+Q1_lastpoint1),ichannel)-I1(ichannel))/Fs);
				
				% The formula above gives the area under the curve with rectangles
				% at every timepoint.  To convert this to trapezoids at every point,
				% we have to subtract the areas of all the little traingles on top
				% every rectangle.  There are Q1_lastpoint1 such triangles, each of
				% same width.  Putting these together two at a time, we get Q1_lastpoint1/2
				% rectangles, which cover (roughly) half the height of current(201).
				% Because the curve stops at current(201+Q1_lastpoint1), we subtract the
				% value of current(201+Q1_lastpoint1)/2.
				
				rect_integral = sum(abs(current((201:201+Q1_lastpoint1),ichannel)-I1(ichannel)));
				trap_subtraction = abs ( ((current(201,ichannel) - I1(ichannel)) / 2 ) - ...
					((current(201+Q1_lastpoint1,ichannel) - I1(ichannel)) / 2 ) );

				Q1 = (rect_integral - trap_subtraction) / Fs;


				% Find the time constant, tau.
				% At t= tau, the current is .368 of its final value.
				% OLD, INCORRECT WAY: tauCurrent = .368 * p(end);
				% (Some recordings were taken this way.)
				%
				% NEW: We take the value 63.2% of the way to the SS current
				% from the peak of the transient spike as the current at tau.
				tauCurrent = .632*(p(end) - I1(ichannel));
				
				% Find indicies on the curve where current is greater than the
				% current at tau (i.e. closer to the steady state current).
				% Comparison sign is > because transient is a downward (-)
				% deflection.
				% +2 is because we started searching at 203 (so
				% currentsAterTau=1 implies tau=(1+2)/Fs=3/Fs msec)
				%
				% Arbitrarily, we've chosen to take tau to be at the first
				% recorded time point where the current *exceeded* the 63.2%
				% current value. We could have just as well used the last time
				% point where the current was *below* the 63.2% value.
				currentsAfterTau = find(current(203:280,ichannel) > tauCurrent);
				
				if length(currentsAfterTau > 0)
					tau = (currentsAfterTau(1)+2)/Fs;
				else
					tau = 0;
				end

				% From Axon calculation
				% abs measures the dI from the two steady-state currents (I2 & I1)
				% Note: Including Q2 in Qt seems to make almost no difference
				% in the final value of Cm.
				Q2 = (abs(I2 - I1)) .* tau;
				% Convert Q since I's are picoamps
				Qt=(Q1+Q2(1)) *10^-12;
				% Qt=(Q1) *10^-12;
				
			end
			
			% Derivations from Axon/Garcia-Diaz Membrane Test
			% rs is the series resistance (pipette + access)
			% rm is the membrane resistance
			% cm is the membrane capacitance
			%
			% The Axon Rs seems less accurate than the transient-based rs
			% calculation: rs2 = tau * dV / Qt;
			rs = dV / ( (mean(current(160:199,ichannel)) - p(end)) * 10^-12);
			rm = rt - rs;
			cm = (Qt*rt)/(dV*rm);
 
			% Another way to calculate rm. Gives similar results to rm above.
			% rm2 = dV / (abs(p(end) - I1(1))*10^-12);
 			
			% An alternate way to estimate cm that seems less accurate
			% on the model cell, true cap is 47 pF. cm=51pF while cm2=62pF
			% cm2 = ((ra+rm)/(ra*rm))*tau;
			
			Rtotal(ichannel) = Rtotal(ichannel) + rt;
			Rseries(ichannel) = Rseries(ichannel) + rs;
			Rmem(ichannel) = Rmem(ichannel) + rm;
			Cmem(ichannel) = Cmem(ichannel) + cm;
			
			% Old way of calculating Rseries (mislabeled as
			% Raccess, which is really Rser - Rpip). 
			% Raccess(ichannel)=Raccess(ichannel) + amplitude / (mean(current(195:199,ichannel))-p(end));
			
		end  % for ichannel            
				
		avgCurrent = avgCurrent + current; 
    
	end  % for irepeat
	
	% Now we calculate the final values that we will send back
	
	% Old calculation for Rseries
	% Raccess = (10 ^9) * Raccess / nrepeat;  % convert from GOhm to Ohms
	Rtotal = Rtotal/nrepeat;                % Rtotal is already in Ohms
	Ileak = Ileak/nrepeat;
	Rseries = Rseries/nrepeat;
	Rmem = Rmem/nrepeat;
	Cmem = Cmem/nrepeat;
	current = avgCurrent/nrepeat;


%---------------------------------------------------------------------------