% adapted from Pan

clear;
% clc;
close all;

%% Options
run_optimisation = true;

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

%% Define constants
R = 8.314;
T = 310;
F = 96485;

%% Plot I-V curves
cKo = 6; % mM
cKi = 140;
C = 7.07; % pF
g_kinetic_per = 0.8; % nS_per_pF
g_kinetic = g_kinetic_per*C;
G_BK = g_kinetic*1e-6; % Unit mA/V = mS

V = (-120:1:60)/1000;
E_Ca = R*T/F*log(cKo/cKi);
I_lin = G_BK*(V-E_Ca); % Unit mA
fit_range = 10:round(length(I_lin)*0.85); % 31:91
error_func = @(G_GHK) square_error(I_lin(fit_range) - calc_IGHK(G_GHK,V(fit_range),cKi,cKo));

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

options_ps = optimoptions('particleswarm','UseParallel',true,'HybridFcn',@fminunc,'SwarmSize',2000);

if run_optimisation
    [G_GHK,fval,exitflag,output] = particleswarm(error_func,1,lb,ub,options_ps);
    save([storage_dir 'BK_G_GHK.mat'],'G_GHK');
else
    load([storage_dir 'BK_G_GHK.mat']);
end
I_GHK = calc_IGHK(G_GHK,V,cKi,cKo);

h = figure;
plot(1000*V,1e6*I_lin,'k--',1000*V,1e6*I_GHK,'k','LineWidth',2);
legend('LRd','BG','Location','southeast');
ylabel('Current (nA)');
xlabel('Voltage (mV)');
set(gca,'FontSize',16);

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