% 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
cCao = 2.5;
cCai_st = 0.000116;
cCa_SR = 0.2417;

g_kinetic = 0.1; % nS
% Original
G_SOC = g_kinetic*1e-6; % Unit mA/V = mS

V = (-120:1:60)/1000;
E_Ca = R*T/F*log(cCao/cCai_st);
I_lin = G_SOC*(V-E_Ca); % Unit mA
fit_range = 20:round(length(I_lin)*0.7); % 31:91
% fit_range = [fit_range, length(I_lin)];
error_func = @(G_GHK) square_error(I_lin(fit_range) - calc_IGHK(G_GHK,V(fit_range),cCai_st,cCao));

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);
    disp(G_GHK);
    save([storage_dir 'SOC_G_GHK.mat'],'G_GHK');
else
    load([storage_dir 'SOC_G_GHK.mat']);
end
I_GHK = calc_IGHK(G_GHK,V,cCai_st,cCao);

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');