function PlotK(x, t, sol, PAR) % Plotting the results % Unwrapping parameters l = PAR(1); A_i = PAR(2); A_o = PAR(3); O_m = PAR(4); g_K0 = PAR(5); g_Na = PAR(6); g_Cl = PAR(7); C_m = PAR(8); c_Ko0 = PAR(9); c_Ki0 = PAR(10); c_Nao0 = PAR(11); c_Nai0 = PAR(12); c_Clo0 = PAR(13); c_Cli0 = PAR(14); v_m0 = PAR(15); X_iZ_i = PAR(16); X_oZ_o = PAR(17); R = PAR(18); F = PAR(19); T = PAR(20); psifac = PAR(21); D_K = PAR(22); D_Na = PAR(23); D_Cl = PAR(24); lambda_i = PAR(25); lambda_o = PAR(26); tstart = PAR(27); tstop = PAR(28); xstart = PAR(29); xstop = PAR(30); Imax = PAR(31); Kdec = PAR(32); R = 8.315; % Gas constant J*mol^-1*K^-1; F = 96500; % Faraday's constant (C/mol); psifac = R*T/F; % Standard potential RT/F(V); c_Ko = sol(:,:,1); c_Ki = sol(:,:,2); c_Nao = sol(:,:,3); c_Nai = sol(:,:,4); c_Clo = sol(:,:,5); c_Cli = sol(:,:,6); v_mo = sol(:,:,7); v_mi = sol(:,:,8); for j = 1:length(t) %Derivative of c_Ko v/ at all t given by duoutdx1 (V/m); [uout1(j,:),duoutdx1(j,:)] = pdeval(0,x,sol(j,:,1),x); [uout2(j,:),duoutdx2(j,:)] = pdeval(0,x,sol(j,:,2),x); [uout3(j,:),duoutdx3(j,:)] = pdeval(0,x,sol(j,:,3),x); [uout4(j,:),duoutdx4(j,:)] = pdeval(0,x,sol(j,:,4),x); [uout5(j,:),duoutdx5(j,:)] = pdeval(0,x,sol(j,:,5),x); [uout6(j,:),duoutdx6(j,:)] = pdeval(0,x,sol(j,:,6),x); [uout7(j,:),duoutdx7(j,:)] = pdeval(0,x,sol(j,:,7),x); end % Input/Output-vectors Kinput = zeros(size(sol(:,:,1))); Koutput = zeros(size(sol(:,:,1))); I_tstart = min(find(t>=tstart)); I_tstop = max(find(t=xstart)); I_xstop = max(find(x<=xstop)); Kinput(I_tstart:I_tstop, I_xstart:I_xstop) = Imax*10^6; %Convert to micromol/(m^2s) Koutput = Kdec*(c_Ko - c_Ko0)*10^6; %Convert to micromol/(m^2s) %Flux densities due to diffusion (mol/(m^2*s)); j_KoD = -(D_K/lambda_o^2)*duoutdx1; j_KiD = -(D_K/lambda_i^2)*duoutdx2; j_NaoD = -(D_Na/lambda_o^2)*duoutdx3; j_NaiD = -(D_Na/lambda_i^2)*duoutdx4; j_CloD = -(D_Cl/lambda_o^2)*duoutdx5; j_CliD = -(D_Cl/lambda_i^2)*duoutdx6; %Current densities due to diffusion (A/m^2); i_odiff = F*(j_KoD + j_NaoD - j_CloD); i_idiff = F*(j_KiD + j_NaiD - j_CliD); %Resistivities (ohm*m) r_o = psifac*lambda_o^2./(F*(D_Na*c_Nao+D_K*c_Ko+D_Cl*c_Clo)); r_i = psifac*lambda_i^2./(F*(D_Na*c_Nai+D_K*c_Ki+D_Cl*c_Cli)); % Gradients of v_o og v_i (V/m) dv_odx = (-duoutdx7 + r_i.*i_idiff+r_i*(A_o/A_i).*i_odiff).*r_o./(r_o+r_i*A_o/A_i); dv_idx = (duoutdx7 + r_o.*(A_i/A_o).*i_idiff+r_o.*i_odiff).*(r_i./(r_i+r_o*A_i/A_o)); %Flux densities due to electrical migration(mol/(m^2*s)); j_KoV = -(D_K/lambda_o^2).*(1/psifac).*c_Ko.*dv_odx; j_KiV = -(D_K/lambda_i^2).*(1/psifac).*c_Ki.*dv_idx; j_NaoV = -(D_Na/lambda_o^2).*(1/psifac).*c_Nao.*dv_odx; j_NaiV = -(D_Na/lambda_i^2).*(1/psifac).*c_Nai.*dv_idx; j_CloV = (D_Cl/lambda_o^2).*(1/psifac).*c_Clo.*dv_odx; j_CliV = (D_Cl/lambda_i^2).*(1/psifac).*c_Cli.*dv_idx; %Current densities due to electrical migration (A/m^2); i_ofield = F*(j_KoV + j_NaoV - j_CloV); i_ifield = F*(j_KiV + j_NaiV - j_CliV); % Total axial current densities (A/m^2); i_i = i_idiff + i_ifield; i_o = i_odiff + i_ofield; %%% Membrane mechanisms P = NaKpumprate(c_Ko, c_Nai); % Na/K pumprate (mol/(m^2*s)); f_Kir = Kirf(c_Ko,c_Ki,c_Ko0,c_Ki0, v_mo, v_m0, psifac); % Kir-factor (unitless) % Membrane-flux densities (mol/(m^2*s)); j_KM = (g_K0.*f_Kir/F).*(v_mo-psifac.*log(c_Ko./c_Ki)) - 2*P; j_NaM = (g_Na/F).*(v_mo-psifac.*log(c_Nao./c_Nai)) + 3*P; j_ClM = -(g_Cl/F).*(v_mo+psifac.*log(c_Clo./c_Cli)); %Membrane current density (A/m^2); i_m = F*(j_KM + j_NaM - j_ClM); %%%% PLOTS: % Defining locations for plots Itstart = min(find(t>tstart)); Itend = max(find(txstart)); Ixend = max(find(x