clear all; close all; clf;
% Constants,Volumes, Patch Area
F=96485; VOL=(5*10^-5)^3; S=(5*10^-4)^2;

%MAX G---->> g(1)=Gk_Max,g(2)=G_Na_Max,g(3)=G_leak_Max
g(1)=36*S; g(2)=120*S; g(3)=0*0.3*S;

%Initial Value For E---->> E(1)=Ek,E(2)=Ena,E(3)=Ecl
E(1)=-59;E(2)=34.66;E(3)=-49.187;

%Loop Variables
%top=top of patch,bot=Bottom of the patch
dt=0.01; tt=2;  t_rec=1; t_final=1800;
I_ext=0; SumIbot=0; SumItop=0;
ZItop=zeros(3,(t_final)*(1/dt)); ZVtop=zeros(1,(t_final)*(1/dt)); INa_Ik_Newtop=zeros(1,(t_final)*(1/dt));
Integrate_Itop=zeros(1,(t_final)*(1/dt)); ten_step_sum = 0;
ZIbot=zeros(3,(t_final)*(1/dt)); ZVbot=zeros(1,(t_final)*(1/dt)); INa_Ik_Newbot =zeros(1,(t_final)*(1/dt));
Integrate_Ibot=zeros(1,(t_final)*(1/dt)); ten_step_sum = 0;

%vm---->>vm=Vm+60,Vm=Membrane voltage!sumI=Ibot+Itop,Itot=Itop+Ibot
vm=0; dV_SumI = 0;
INTEGRATE_Itot=zeros(1,(t_final)*(1/dt));

%Calculate rest voltage & Initial Value For Gn,m,h 
    alpha_top(1)=0.01*(10-vm)/((exp((10-vm)/10)-1));
    alpha_top(2)=0.1*(25-vm)/((exp((25-vm)/10)-1));
    alpha_top(3)=0.07*exp(-vm/20);
    beta_top(1)=0.125*exp(-vm/80);
    beta_top(2)=4*exp(-vm/18);
    beta_top(3)=1/(exp((30-vm)/10)+1);
    tau_top=1./(alpha_top+beta_top);
    x_inf=alpha_top.*tau_top;
    x=x_inf; x1=x;
    gnmh_top(1)=g(1)*x_inf(1)^4; gnmh_top(2)=g(2)*x_inf(2)^3*x_inf(3);gnmh_top(3)=g(3);
    Vrest =(gnmh_top*E')/(sum (gnmh_top));
    Vtop=Vrest;
   
for t=0:dt:t_final
    
    
    %dV_sin = 5*(1+ 0*sin(2*pi*0.01*t)).*sin(2*pi*1*t).*(t<900).*(t>100);
     dV_sin = 5*sin(2*pi*1*t).*(t<900).*(t>100);
    
    %-----------------top patch-----------------%
    Vtop=Vrest+dV_sin+dV_SumI;
    vm=Vtop-Vrest; 
    
    %Calculate Gn,m,h & I channle_na_k for top of patch
    alpha_top(1)=0.01*(10-vm)/((exp((10-vm)/10)-1));
    alpha_top(2)=0.1*(25-vm)/((exp((25-vm)/10)-1));
    alpha_top(3)=0.07*exp(-vm/20);
    beta_top(1)=0.125*exp(-vm/80);
    beta_top(2)=4*exp(-vm/18);
    beta_top(3)=1/(exp((30-vm)/10)+1);
    tau_top=1./(alpha_top+beta_top);
    x_inf=alpha_top.*tau_top;
    x=(1-dt./tau_top).*x+dt./tau_top.*x_inf;
    gnmh_top(1)=g(1)*x(1)^4; gnmh_top(2)=g(2)*x(2)^3*x(3);gnmh_top(3)=g(3);
    I_top=gnmh_top.*(Vtop-E);  ZItop(:,tt) =I_top' ;  ZVtop(:,tt) =Vtop ;%1=k, 2=na, 3=cl
    INa_Ik_lasttop(1,tt) = ZItop(1,tt-1)+ZItop(2,tt-1)+ZItop(3,tt-1);
    INa_Ik_Newtop(1,tt) = ZItop(1,tt)+ZItop(2,tt)+ZItop(3,tt);
    SumItop = SumItop+(INa_Ik_Newtop(1,tt)+INa_Ik_lasttop(1,tt))*(dt/2);
    Integrate_Itop(1,tt)= SumItop ;
    
    %-----------------bottom patch-----------------%
    Vbot=Vrest-dV_sin+dV_SumI;
    vmbot=Vbot-Vrest; 
    
    %Calculate Gn,m,h & I channle_na_k for bottom of patch
    alpha_bot(1)=0.01*(10-vmbot)/((exp((10-vmbot)/10)-1));
    alpha_bot(2)=0.1*(25-vmbot)/((exp((25-vmbot)/10)-1));
    alpha_bot(3)=0.07*exp(-vmbot/20);
    beta_bot(1)=0.125*exp(-vmbot/80);
    beta_bot(2)=4*exp(-vmbot/18);
    beta_bot(3)=1/(exp((30-vmbot)/10)+1);
    tau_bot=1./(alpha_bot+beta_bot);
    x_inf1=alpha_bot.*tau_bot;
    x1=(1-dt./tau_bot).*x1+dt./tau_bot.*x_inf1;
    gnmh_bot(1)=g(1)*x1(1)^4; gnmh_bot(2)=g(2)*x1(2)^3*x1(3);gnmh_bot(3)=g(3);
    I_bot=gnmh_bot.*(Vbot-E);  ZIbot(:,tt) =I_bot' ;  ZVbot(:,tt) =Vbot ;
    INaIk_lastbot(1,tt) = ZIbot(1,tt-1)+ZIbot(2,tt-1)+ZIbot(3,tt-1);
    INa_Ik_Newbot(1,tt) = ZIbot(1,tt)+ZIbot(2,tt)+ZIbot(3,tt);
    SumIbot = SumIbot+(INa_Ik_Newbot(1,tt)+INaIk_lastbot(1,tt))*(dt/2);
    Integrate_Ibot(1,tt)= SumIbot ;

     %-----------------Voltage-change of the Membrane-----------------%
    C=1;
    dV_SumI = -1*(((Integrate_Itop(1,tt))/(C*6*S))+((Integrate_Ibot(1,tt))/(C*6*S)));
    INTEGRATE_Itot(1,tt)=Integrate_Itop(1,tt)+Integrate_Ibot(1,tt);  
     
	%-----------------Ion Concentrations-----------------%   

	%	(You write it!)




     %-----------------------Registeration for later visualisation-----------------------%
    if t>=0;
        t_rec=t_rec+1; x_plot(t_rec)=t; y_plot(t_rec)=Vtop;
        G(t_rec,1)=gnmh_top(1); % GK
        G(t_rec,2)=gnmh_top(2); % GNa
        n_particle(t_rec)=x(1); m_particle(t_rec)=x(2);h_particle(t_rec)=x(3);
        
    end
    tt=tt+1;
end

%-----------------Visualization- Voltage Membrane Curve-----------------%
plot(x_plot,y_plot); xlabel('Time (ms)'); ylabel('Membrane Voltage (mV)');
Fs = 1/(dt*1e-3);% Sampling frequency                    
L = length(x_plot);% Length of signal
% t = (0:L-1)*dt;% Time vector
VMFFT=fft(y_plot);
P2 = abs(VMFFT/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = Fs*(0:(L/2))/L;

%---------------------Visualization-dV_SumI Curve--------------------%
figure;
plot(f,P1) 
title('Single-Sided Amplitude Spectrum of X(t)')
xlabel('f (Hz)')
ylabel('|P1(f)|')
figure;
plot(x_plot,-1/(C*6*S)*INTEGRATE_Itot); xlabel('Time (ms)'); ylabel('Vmem_Demodulated');