% curve fitting for i_ns_Ca, which consists of 3 currents (each carrying
% Ca, Na, K)
% 3 GHK fits must be done to find 3 GHK conductances

clear;
% clc;
% close all;

%% Options
run_optimisation = true;

%% Set up directories
current_dir = pwd;
Idx_backslash = find(current_dir == filesep);
main_dir = current_dir; %(1:Idx_backslash(end));
output_dir = [main_dir '\output' filesep];
storage_dir = current_dir; %[main_dir '\storage' filesep];

%% Define constants - coonsistent with VOLTS.
R = 8.314;
T = 310;
F = 96485;

%% choose which ion to plot current of
ion = 'Ca'; %Ca, Na, K

%% Plot I-V curves
% UNITS:
%     I = G_GHK * [K] [=] mA
%     G_GHK = I/[K] [=] Amp.litre/mol
%     G_lin [=] mS [=] mA/V
    
V = (-110:1:65)/1000;
V(find(V == 0,1)) = [];
V = V';

cCao = 1.8;
cCai = 6e-5; %millimolar
cNao = 140;
cNai = 10;
cKo = 5.4;
cKi = 145;
Cm = 153400e-9; % Unit microF
SA = 11400e-6; % [=] mm2, using SA(um2) = 0.3*vol (in um3) (Bers)
P_ns_Ca = 1.75e-7; % [=] cm_per_second
K_m_ns_Ca = 0.0012; % [=] mM
gamma = 0.75; % [=] dimensionless

if true
    % offset, the same as LRd94 code
    EnsCa = R*T/F*log((cKo+cNao)/(cKi+cNai));
    Vns = V - EnsCa;
else
    Vns = V;
end

% the following currents in units [=] uA_per_mm2
I_ns_Na = P_ns_Ca*Vns*(F^2)/(R*T).*(gamma*cNai*exp(1*Vns*F/(R*T))-gamma*cNao)./(exp(Vns*F/(R*T))-1);
I_ns_K = P_ns_Ca*Vns*(F^2)/(R*T).*(gamma*cKi*exp(1*Vns*F/(R*T))-gamma*cKo)./(exp(Vns*F/(R*T))-1);
i_ns_Na = I_ns_Na./(1+(K_m_ns_Ca/cCai ^ 3));
i_ns_K = I_ns_K./(1+(K_m_ns_Ca/cCai ^ 3));
i_ns_Ca = i_ns_Na+i_ns_K;

% convert uA/mm2 to unit mA
% assume 1mm2 area == 1 uF
i_ns_Ca_lrd = 1e-3*SA*i_ns_Ca; % sign flip not required. this current in same dir as other currents e.g. fastNa
i_ns_Na_lrd = 1e-3*SA*i_ns_Na;
i_ns_K_lrd = 1e-3*SA*i_ns_K;

Vstart = 1;
Vend =  length(V);
switch ion
    case 'Ca'
        z = 2;
        error_func = @(G_GHK) square_error(i_ns_Ca_lrd(Vstart:Vend) - calc_IGHK(G_GHK,V(Vstart:Vend),cCai,cCao, z));
        mat_name = '\nsCa_G_GHK.mat';
    case 'Na'
        z = 1;
        error_func = @(G_GHK) square_error(i_ns_Na_lrd(Vstart:Vend) - calc_IGHK(G_GHK,V(Vstart:Vend),cNai,cNao, z));
        mat_name = '\nsNa_G_GHK.mat';
    case 'K'
        z = 1;
        error_func = @(G_GHK) square_error(i_ns_K_lrd(Vstart:Vend) - calc_IGHK(G_GHK,V(Vstart:Vend),cKi,cKo, z));
        mat_name = '\nsK_G_GHK.mat';
    otherwise
        error('not an ion')
end

A = [];
b = [];
Aeq = [];
beq = [];
lb = [-Inf];
ub = [Inf];

options_ps = optimoptions('particleswarm','UseParallel',false,'HybridFcn',@fminunc,'SwarmSize',1500, ...
    'FunctionTolerance', 1e-16);

if run_optimisation
    [G_GHK,fval,exitflag,output] = particleswarm(error_func,1,lb,ub,options_ps);
    save([storage_dir mat_name],'G_GHK');
else
    load([storage_dir mat_name]);
end


h = figure;
% subplot(1,2,1)
if strcmp(ion,'Ca')
    I_GHK = calc_IGHK(G_GHK,V,cCai,cCao,z);
    plot(1000*V,1e6*i_ns_Ca_lrd,'k--',1000*V,1e6*I_GHK);
elseif strcmp(ion,'Na')
    I_GHK = calc_IGHK(G_GHK,V,cNai,cNao,z);
    plot(1000*V,1e6*i_ns_Na_lrd,'k--',1000*V,1e6*I_GHK);
else
    I_GHK = calc_IGHK(G_GHK,V,cKi,cKo,z);
    plot(1000*V,1e6*i_ns_K_lrd,'k--',1000*V,1e6*I_GHK);
end
legend('LRd','GHK fit','Location','southeast');
ylabel('Current (nA)');
xlabel('Voltage (mV)');
title(ion)
set(gca,'FontSize',16);
grid on
% subplot(1,2,2)
% plot(1000*V,1e6*I_lin,'k--',1000*V, 1e6*I_GHK2);
% legend('LRd','BG_test','Location','southeast');
% ylabel('Current (nA)');
% xlabel('Voltage (mV)');
% set(gca,'FontSize',16);

% print_figure(h,output_dir,'tcc_IV_curve');