%HW4
% APMA 0410
% Fall 2012

clear

dt = 0.0001;    % time step                             [sec]
R = 20e6;       % membrane's resistance                 [Ohm]
C = 1e-9;       % membrane's capacitance                [F] 
V_rest = -0.07; % resting potential                     [V]

tau = R*C;      % membrane time constant                [sec] 



%% PROBLEM 4.4
%  Constant injected current, starting at 200 msec

t_max = 0.4;    % end of recording period               [sec]
time_axis = 0:dt:t_max;
n_time_steps = length(time_axis);

I_injected = 2e-9;                     % amplitude of injected current    [A]
t_start = 0.2/dt;                       % current-injection starting time  [bins]

% We pre-allocate the array containing the injected current as a function of time,
%    setting it to 0 throughout the recording period: 
I_values = zeros(1,n_time_steps);     
% We set the injected current to its desired non-zero value from t_start to the end: 
I_values(t_start:n_time_steps) = I_injected;

V = zeros(1,n_time_steps);              % pre-allocate the membrane potential array
 
V_t = V_rest;                           % Initial value of V is resting potential
                                        %    V_t will be updated at each t

for t = 1 : n_time_steps
    
%     derivative = (V_rest - V_t)/tau + I_values(t)/C;    % derivative at time t
    derivative = (V_rest - V_t + R*I_values(t))/tau;    % equivalent form
    V_t = V_t + derivative*dt;                          % update the membrane potential
    V(t) = V_t;                                         % store the membrane potential
    
end

figure(1)
plot(time_axis,V)
title(['Constant Injected Current'])
xlabel('time   [ sec ]')
ylabel('membrane potential   [ V ]')
axis([0 t_max -.15 .15])
grid on


%% PROBLEM 4.5
%  Sinusoidal current

I_amplitude = 1*1e-9;   % amplitude of sinusoidal injected current [A]
                        %     Note: I_amplitude is positive.

I_frequency = 10;       % frequency of sinusoidal injected current [Hz]

t_max = .5;            % end of recording period               [sec]
time_axis = 0:dt:t_max;
n_time_steps = length(time_axis);

V = zeros(1,n_time_steps);              % pre-allocate the membrane potential array
 
V_t = V_rest;                           % Initial value of V is resting potential
                                        %    V_t will be updated at each t


% We create the array containing the injected current as a function of time:
I_values = I_amplitude*cos(2*pi*I_frequency*time_axis);
    

for t = 1 : n_time_steps
    
    derivative = (V_rest - V_t)/tau + I_values(t)/C;    % derivative at time t
    V_t = V_t + derivative*dt;                          % update the membrane potential
    V(t) = V_t;                                         % store the membrane potential
    
end

figure(2)
subplot(2,1,1)
plot(time_axis,V);
axis([0 t_max -.09 -.05])
title(['driving current:  amplitude =',num2str(I_amplitude*10^9),...
    ' nA    frequency = ',num2str(I_frequency),' Hz'])
grid on
ylabel('membrane potential [ V ]')
subplot(2,1,2)
plot(time_axis,I_values);
ylabel('injected current [ A ]')
xlabel('time [ sec ]')
grid on
