clear;
% clc;
% close all;

%% Options
run_optimisation = true;

%% Set up directories
current_dir = cd;
Idx_bactolash = find(current_dir == filesep);
main_dir = current_dir; %(1:Idx_bactolash(end));
data_dir = [main_dir '\data' filesep];
code_dir = [main_dir '\code' 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
% UNITS:
%     I = G_GHK * K [=] mA
%     G_GHK = I/K [=] Amp.litre/mol
%     G_PMR [=] mS [=] mA/V
    
V = (-120:1:60)/1000;

cKo = 4.5;
cKi = 104.75;
Cm = 60e-6; % Kernik  % 153400e-9; % Unit microF

E_K = R*T/F*log(cKo/cKi);

G_kr = 0.218025*sqrt(cKo/5.4)*Cm; % mA/V
I_lin = G_kr*(V-E_K); % Unit mA    don't consider gating variable: just finding conductance.

Vstart = 1;
Vend = 181;
error_func = @(G_GHK) square_error(I_lin(Vstart:Vend) - calc_IGHK(G_GHK,V(Vstart:Vend),cKi,cKo));

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

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

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

% G_GHK2 = 0.003;

I_GHK = calc_IGHK(G_GHK,V,cKi,cKo);
% I_GHK2 = calc_IGHK(G_GHK2,V,cKi,cKo);

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

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