%Calcolare e disegnare andamento di alpha beta n m h e delle conduttanze potassio e sodio al variare di Vm nell'intervallo [-100 : 40] mV
%condizioni: assone gigante di calamaro T=6.3 C

clear variables 
close all
clc

%% costanti

GNaM=120; %[mS/cm2]
GKM=36; %[mS/cm2]
GL=0.3; %[mS/cm2]
CM=1; %[uF/cm2]
Vr=-65; %[mv]
VNa=50; %[mv]
VK=-77; %[mv]
VL=-54.387; %[mv]



%% calcolo e grafico dei coefficienti di velocità di transizione [1/ms] e delle conduttanze al variare di Vm

Vm=-100:0.1:40; % intervallo variazione potenziale di membrana
dv=Vm-Vr; %Vm-Vr=V' variazione rispetto al potenziale di riposo 

%inizializzo le variabili alpha_nmh e beta_nmh
a_n=zeros(size(Vm));
b_n=zeros(size(Vm));
a_m=zeros(size(Vm));
b_m=zeros(size(Vm));
a_h=zeros(size(Vm));
b_h=zeros(size(Vm));


for i=1:length(Vm) % applico le formule di calcolo degli alpha e beta in funzione di V' (Vm-Vr)
a_n(i)=(0.01*(10-(dv(i))))/(exp((10-dv(i))/10)-1); 
b_n(i)=0.125/(exp(dv(i)/80));
a_m(i)=(0.1*(25-(dv(i))))/(exp((25-dv(i))/10)-1); 
b_m(i)=4/(exp(dv(i)/18));
a_h(i)=0.07/(exp(dv(i)/20));
b_h(i)=1/(exp((30-dv(i))/10)+1); 
end

%% potassio
%coefficienti alpha_n betha_n al variare di Vm
figure
plot(Vm,a_n,'b');
hold on
plot(Vm,b_n,'r');
xlabel('v [mv]');
legend('\alpha_n','\beta_n');


%calcolo valore di regime di n al variare di Vm
n_inf=a_n./(a_n+b_n);
figure
plot(Vm,n_inf,'b');
hold on
plot(Vm,n_inf.^4,'r');
xlabel('v [mv]');
legend('n_{inf}','n_{inf}^4');

%plotto Gk (regime) al variare di Vm
figure
plot(Vm,GKM*(n_inf.^4));
xlabel('v [mv]');
legend('G_K');

%% sodio

%coefficienti alpha_n betha_n al variare di Vm
figure
plot(Vm,a_m,'b');
hold on
plot(Vm,b_m,'r');
xlabel('v [mv]');
legend('\alpha_m','\beta_m');

%coefficienti alpha_h betha_h variare di Vm
figure
plot(Vm,a_h,'b');
hold on
plot(Vm,b_h,'r');
xlabel('v [mv]');
legend('\alpha_h','\beta_h');



%calcolo valore di regime al variare di Vm
m_inf=a_m./(a_m+b_m);
h_inf=a_h./(a_h+b_h);
figure
subplot(2,1,1)
plot(Vm,m_inf,'b');
hold on
plot(Vm,h_inf,'r');
xlabel('v [mv]');
legend('m_{inf}','h_{inf}');
subplot(2,1,2);
plot(Vm,(m_inf.^3).*h_inf);
xlabel('v [mv]');
legend('m_{inf}^3 h_{inf}');

%plotto GNa (regime) al variare di Vm
figure
plot(Vm,GNaM*((m_inf.^3).*h_inf));
xlabel('v [mv]');
legend('G_{Na}');

