% Ion transport noise model dynamics for nano porous membranes

clear all;
clc;

% Input parameters

% Simulation parameters
D = 1.33e-9 ; % Diffusivity of Sodium
K_b = 1.38064e-23; % in kg m^2/s^2K
T = 300; % Temperature in K
lambda = (K_b*T)/D ; % Mass rate in (kg/s)
omega = D/(K_b*T) ; % Mobility in (s/kg)
epsilon_0 = 8.854e-12; % Dielectric permittivity of free space in F/m
epsilon_r = 80; % Dielectric permittivity of water (no units)
e_charge = 1.60217e-19; % electronic Charge

% Geometry of the micro-nano gating membrane

r_nano = 1.78e-9; % Radius of nanopore in nm
L_nano = 4.64e-9; % Length of nanopore in nm
r_micro = 3.095e-9; % Radius of micropore in nm
L_micro = 6.67e-9; % Length of micropore in nm
r_macro = 3.095e-9; % Radius of micropore in nm
L_macro = 6.67e-9; % Length of micropore in nm

A_micro = (pi*(r_micro)^2);
A_nano = (pi*(r_nano)^2);
A_macro = (pi*(r_macro)^2);
dS_nano = (2*pi*r_nano*L_nano);
dV_nano = A_nano*L_nano;
dV_micro = A_micro*L_micro;
dV_macro = A_macro*L_macro;

% Micro-nano-gating membrane parameters
sigma_micro = 0e-3; % surface chage of micropore in C/m^3
sigma_nano = -200e-3; % surface charge of nanopore in C/m^3
sigma_macro = 0; % surface charge of macropore in C/m^3
q_s_micro = sigma_micro*(2*pi*r_micro*L_micro); % Charge of micropore
q_s_nano = sigma_nano*(2*pi*r_nano*L_nano); % charge of nanopore
q_s_macro = sigma_macro*(2*pi*r_macro*L_macro); % charge of macropore

c_0 = 1; % concentration of NaHPO4/Ha2HPO4 = 1 mM
NA = 6.023e23; % Avagadro number
if (q_s_nano ~= 0)
    c_nano = (q_s_nano*-1)/(e_charge*dV_nano*NA); % ideal counter ion concentration in nanopore 
else if (q_s_nano == 0)
        c_nano = c_0;
    end 
end
n_V_micro = (c_0*NA); % number of molecules per m^3 in micropore
n_V_nano = (c_nano*NA); % Number of molecules per m^3 in nanopore
n_V_macro = (c_0*NA); % Number of molecules per m^3 in macropore


% Total mass of the Charged particle
% mass_Sodium = 22.9898e-3; % Kg/mol
% mass_micro = ((mass_Sodium*n_V_micro*dV_micro)/NA); % mass of one charged ion
% mass_nano = ((mass_Sodium*n_V_nano*dV_nano)/NA); % mass of one charged ion
% mass_macro = ((mass_Sodium*n_V_macro*dV_macro)/NA); % mass of one charged ion

% More than concentration mass of ion is used as a parameter to obtain the
% noise dynamics and I-V characteristics in nano porous membrane transport 

mass_micro = 3.848e-6; % mass of one charged ion
mass_nano = 3.848e-6; % mass of one charged ion
mass_macro = 3.848e-6; % mass of one charged ion


% Simulation time step
t_new = [1e-10:3e-10:100e-9];
N_size = size(t_new,2);

% Applied voltage
Amplitude = 3; % Voltage in (V)
V_0 = 3; % Normalized voltage (V)
% Voltage at each time
Voltage = ones(N_size,1)*Amplitude;

% Derived parameters
p = 4*pi*epsilon_0*epsilon_r; % constant
Gamma_micro = (1+((lambda*(t_new(1)-0))/(2*mass_micro)));
zeta_micro = (1-((lambda*(t_new(1)-0))/(2*mass_micro)));
Gamma_nano = (1+((lambda*(t_new(1)-0))/(2*mass_nano)));
zeta_nano = (1-((lambda*(t_new(1)-0))/(2*mass_nano)));
Gamma_macro = (1+((lambda*(t_new(1)-0))/(2*mass_macro)));
zeta_macro = (1-((lambda*(t_new(1)-0))/(2*mass_macro)));
x_0 = L_micro; % Length of micropore
x_b1 = 0;
x_b2 = L_micro+L_nano+L_macro;


% Initial guess of each parameter
q_new(1) = (q_s_nano*-1); % Initial total counter ionic charge (GUESS 1)
%q_new(1) = (1e-21)*(+1); % Counter ionic charge (+ ve)
x_variable(1) = 1e-9; % initial position of the charged particle
v_new(1) = 0.001; % initial velocity = 1e-15/1e-12 (electron_radius/time)
acc_new(1) = (q_s_micro*q_new(1))/(mass_micro*p*abs(x_variable(1)-x_0)^2) + (q_new(1)*Voltage(1))/(mass_micro*(x_b2-x_b1)) - (lambda*v_new(1))/mass_micro ; 

% Non-Dimensional numbers
Du_number = (abs(q_s_nano)/(c_0*NA*dV_nano*e_charge))
Nandigana_number = (abs(sigma_nano)*dS_nano)/(4*pi*epsilon_r*epsilon_0*1e-9*V_0)

for i = 2:N_size
% position of the particle
  x_variable(i) = x_variable(i-1) + (v_new(i-1)*(t_new(i)-t_new(i-1))) + (0.5*(acc_new(i-1)*(t_new(i)-t_new(i-1))^2));
 % IF PARTICLE IN MICROPORE
 if (x_variable(i) >=0 && x_variable(i) <= L_micro) 
  
 B1(i) = (q_s_micro*(t_new(i)-t_new(i-1)))/(2*mass_micro*p*Gamma_micro*L_micro*abs(x_variable(i)-x_0)^2);
 B2(i) = (Voltage(i)*(t_new(i)-t_new(i-1)))/(2*mass_micro*Gamma_micro*L_micro*abs(x_b2-x_b1));
 A_new(i) = B1(i)+B2(i);

 M(i) = (zeta_micro*v_new(i-1))/(Gamma_micro*L_micro);
 C1(i) = (q_s_micro*q_new(i-1)*(t_new(i)-t_new(i-1)))/(2*mass_micro*p*Gamma_micro*L_micro*abs(x_variable(i-1)-x_0)^2);
 C2(i) = (q_new(i-1)*Voltage(i-1)*(t_new(i)-t_new(i-1)))/(2*mass_micro*Gamma_micro*abs(x_b2-x_b1)*L_micro);
 E(i) = -1/(t_new(i)-t_new(i-1));
 
 B_new(i) = M(i)+C1(i)+C2(i)+E(i) ; 
 
 C_new(i) = q_new(i-1)/(t_new(i)-t_new(i-1)) ;
 
 % Charge of the particle in micropore
  q_new(i) = (-B_new(i) - sqrt((B_new(i)^2)-(4*A_new(i)*C_new(i))))/(2*A_new(i));
 
  % position of the particle
  x_variable(i) = x_variable(i-1) + (v_new(i-1)*(t_new(i)-t_new(i-1))) + (0.5*(acc_new(i-1)*(t_new(i)-t_new(i-1))^2));
 % velocity of particle
 v_new(i) = (zeta_micro*v_new(i-1))/Gamma_micro + (q_s_micro*q_new(i-1)*(t_new(i)-t_new(i-1)))/(2*mass_micro*p*Gamma_micro*abs(x_variable(i-1)-x_0)^2) + (q_new(i-1)*Voltage(i-1)*(t_new(i)-t_new(i-1)))/(2*mass_micro*Gamma_micro*abs(x_b2-x_b1)) + (q_s_micro*q_new(i)*(t_new(i)-t_new(i-1)))/(2*mass_micro*p*Gamma_micro*abs(x_variable(i)-x_0)^2) + (q_new(i)*Voltage(i)*(t_new(i)-t_new(i-1)))/(2*mass_micro*Gamma_micro*abs(x_b2-x_b1));
 % acceleration of the particle
 acc_new(i-1) = (q_s_micro*q_new(i-1))/(mass_micro*p*abs(x_variable(i-1)-x_0)^2) + (q_new(i-1)*Voltage(i-1))/(mass_micro*abs(x_b2-x_b1)) - (lambda*v_new(i-1))/(mass_micro);
 acc_new(i) = (q_s_micro*q_new(i))/(mass_micro*p*abs(x_variable(i)-x_0)^2) + (q_new(i)*Voltage(i))/(mass_micro*abs(x_b2-x_b1)) - (lambda*v_new(i))/(mass_micro);

 % Current in the micropore
   current(i) = (v_new(i)*q_new(i))/L_micro;
  
         
    
       IN1 = wiener2(current,[i i])-2875*wiener2(current,[i i])+wiener2(current,[i i])-3490*wiener2(current,[i i])-2840*wiener2(current,[i-20 i-20])+wiener2(current,[i-20 i-20])+500*wiener2(current,[i-20 i-20])-3470*wiener2(current,[i-80 i-80])-wiener2(current,[i-80 i-80])-2850*wiener2(current,[i-80 i-80])-2190*wiener2(current,[i-80 i-80])-wiener2(current,[i-80 i-80])-7480*wiener2(current,[i-180 i-180])-6870*wiener2(current,[i-180 i-180])-wiener2(current,[i-180 i-180])-8540*wiener2(current,[i-180 i-180])-7940*wiener2(current,[i-180 i-180])-17870*wiener2(current,[i-230 i-230])-wiener2(current,[i-230 i-230])-18540*wiener2(current,[i-230 i-230])-12350*wiener2(current,[i-230 i-230])-21870*wiener2(current,[i-50 i-50])-wiener2(current,[i-50 i-50])-25540*wiener2(current,[i-50 i-50])-18650*wiener2(current,[i-50 i-50])-7684*wiener2(current,[i-20 i-20])-wiener2(current,[i-20 i-20])-8165*wiener2(current,[i-20 i-20])-6576*wiener2(current,[i-20 i-20]);
       IN2 = wiener2(current,[i-40 i-40])+8800*wiener2(current,[i-40 i-40])+wiener2(current,[i-40 i-40])+9500*wiener2(current,[i-40 i-40])+wiener2(current,[i-40 i-40])+9500*wiener2(current,[i-40 i-40])+88*wiener2(current,[i-80 i-80])+wiener2(current,[i-80 i-80])+95*wiener2(current,[i-80 i-80])+wiener2(current,[i-80 i-80])+95*wiener2(current,[i-80 i-80])+wiener2(current,[i-40 i-40])+8800*wiener2(current,[i-40 i-40])+wiener2(current,[i-40 i-40])+9500*wiener2(current,[i-40 i-40])+wiener2(current,[i-40 i-40])+9500*wiener2(current,[i-40 i-40])+95*wiener2(current,[i-80 i-80])+wiener2(current,[i-40 i-40])+88*wiener2(current,[i-40 i-40])+wiener2(current,[i-180 i-180])+6800*wiener2(current,[i-180 i-180])+wiener2(current,[i-180 i-180])+9400*wiener2(current,[i-180 i-180])+21300*wiener2(current,[i-230 i-230])+wiener2(current,[i-230 i-230])+24800*wiener2(current,[i-230 i-230])+25700*wiener2(current,[i-40 i-40])+wiener2(current,[i-40 i-40])+27900*wiener2(current,[i-40 i-40])+4520*wiener2(current,[i-50 i-50])+2850*wiener2(current,[i-50 i-50]);
       
       INoise = IN1+IN2+IN1+IN2+IN1+IN2+IN1+IN2+IN1+IN2+IN1+IN2+IN1+IN2+IN1+IN2+IN1+IN2+IN1+IN2+IN1+IN2+IN1+IN2+IN1+IN2+IN1+IN2+IN1+IN2+IN1+IN2;

 
 
 
 
 
 
 
 
 end
    % IF PARTICLE IN NANOPORE      
 if (x_variable(i) > L_micro && x_variable(i) <= L_micro+L_nano ) 
  
 B1(i) = (q_s_nano*(t_new(i)-t_new(i-1)))/(2*mass_nano*p*Gamma_nano*L_nano*abs(x_variable(i)-x_0)^2);
 B2(i) = (Voltage(i)*(t_new(i)-t_new(i-1)))/(2*mass_nano*Gamma_nano*L_nano*abs(x_b2-x_b1));
 A_new(i) = B1(i)+B2(i);

 M(i) = (zeta_nano*v_new(i-1))/(Gamma_nano*L_nano);
 C1(i) = (q_s_nano*q_new(i-1)*(t_new(i)-t_new(i-1)))/(2*mass_nano*p*Gamma_nano*L_nano*abs(x_variable(i-1)-x_0)^2);
 C2(i) = (q_new(i-1)*Voltage(i-1)*(t_new(i)-t_new(i-1)))/(2*mass_nano*Gamma_nano*abs(x_b2-x_b1)*L_nano);
 E(i) = -1/(t_new(i)-t_new(i-1));
 
 B_new(i) = M(i)+C1(i)+C2(i)+E(i) ; 
 
 C_new(i) = q_new(i-1)/(t_new(i)-t_new(i-1)) ;
 
 % Charge of the particle in nanopore
  q_new(i) = (-B_new(i) - sqrt((B_new(i)^2)-(4*A_new(i)*C_new(i))))/(2*A_new(i));
 
  % position of the particle
  x_variable(i) = x_variable(i-1) + (v_new(i-1)*(t_new(i)-t_new(i-1))) + (0.5*(acc_new(i-1)*(t_new(i)-t_new(i-1))^2));
 % velocity of particle
 v_new(i) = (zeta_nano*v_new(i-1))/Gamma_nano + (q_s_nano*q_new(i-1)*(t_new(i)-t_new(i-1)))/(2*mass_nano*p*Gamma_nano*abs(x_variable(i-1)-x_0)^2) + (q_new(i-1)*Voltage(i-1)*(t_new(i)-t_new(i-1)))/(2*mass_nano*Gamma_nano*abs(x_b2-x_b1)) + (q_s_nano*q_new(i)*(t_new(i)-t_new(i-1)))/(2*mass_nano*p*Gamma_nano*abs(x_variable(i)-x_0)^2) + (q_new(i)*Voltage(i)*(t_new(i)-t_new(i-1)))/(2*mass_nano*Gamma_nano*abs(x_b2-x_b1));
 % acceleration of the particle
 acc_new(i-1) = (q_s_nano*q_new(i-1))/(mass_nano*p*abs(x_variable(i-1)-x_0)^2) + (q_new(i-1)*Voltage(i-1))/(mass_nano*abs(x_b2-x_b1)) - (lambda*v_new(i-1))/(mass_nano);
 acc_new(i) = (q_s_nano*q_new(i))/(mass_nano*p*abs(x_variable(i)-x_0)^2) + (q_new(i)*Voltage(i))/(mass_nano*abs(x_b2-x_b1)) - (lambda*v_new(i))/(mass_nano);

 % Current in the nanopore
   current(i) = (v_new(i)*q_new(i))/L_nano;
   
       
        IN1 = wiener2(current,[i i])-2875*wiener2(current,[i i])+wiener2(current,[i i])-3490*wiener2(current,[i i])-2840*wiener2(current,[i-20 i-20])+wiener2(current,[i-20 i-20])+500*wiener2(current,[i-20 i-20])-3470*wiener2(current,[i-80 i-80])-wiener2(current,[i-80 i-80])-2850*wiener2(current,[i-80 i-80])-2190*wiener2(current,[i-80 i-80])-wiener2(current,[i-80 i-80])-7480*wiener2(current,[i-180 i-180])-6870*wiener2(current,[i-180 i-180])-wiener2(current,[i-180 i-180])-8540*wiener2(current,[i-180 i-180])-7940*wiener2(current,[i-180 i-180])-17870*wiener2(current,[i-230 i-230])-wiener2(current,[i-230 i-230])-18540*wiener2(current,[i-230 i-230])-12350*wiener2(current,[i-230 i-230])-21870*wiener2(current,[i-50 i-50])-wiener2(current,[i-50 i-50])-25540*wiener2(current,[i-50 i-50])-18650*wiener2(current,[i-50 i-50])-7684*wiener2(current,[i-20 i-20])-wiener2(current,[i-20 i-20])-8165*wiener2(current,[i-20 i-20])-6576*wiener2(current,[i-20 i-20]);
       IN2 = wiener2(current,[i-40 i-40])+8800*wiener2(current,[i-40 i-40])+wiener2(current,[i-40 i-40])+9500*wiener2(current,[i-40 i-40])+wiener2(current,[i-40 i-40])+9500*wiener2(current,[i-40 i-40])+88*wiener2(current,[i-80 i-80])+wiener2(current,[i-80 i-80])+95*wiener2(current,[i-80 i-80])+wiener2(current,[i-80 i-80])+95*wiener2(current,[i-80 i-80])+wiener2(current,[i-40 i-40])+8800*wiener2(current,[i-40 i-40])+wiener2(current,[i-40 i-40])+9500*wiener2(current,[i-40 i-40])+wiener2(current,[i-40 i-40])+9500*wiener2(current,[i-40 i-40])+95*wiener2(current,[i-80 i-80])+wiener2(current,[i-40 i-40])+88*wiener2(current,[i-40 i-40])+wiener2(current,[i-180 i-180])+6800*wiener2(current,[i-180 i-180])+wiener2(current,[i-180 i-180])+9400*wiener2(current,[i-180 i-180])+21300*wiener2(current,[i-230 i-230])+wiener2(current,[i-230 i-230])+24800*wiener2(current,[i-230 i-230])+25700*wiener2(current,[i-40 i-40])+wiener2(current,[i-40 i-40])+27900*wiener2(current,[i-40 i-40])+4520*wiener2(current,[i-50 i-50])+2850*wiener2(current,[i-50 i-50]);
       
       INoise = IN1+IN2+IN1+IN2+IN1+IN2+IN1+IN2+IN1+IN2+IN1+IN2+IN1+IN2+IN1+IN2+IN1+IN2+IN1+IN2+IN1+IN2+IN1+IN2+IN1+IN2+IN1+IN2+IN1+IN2+IN1+IN2;

 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
     
 end    
   
 % IF Particle in the macropore
 
 if (x_variable(i) > L_micro+L_nano && x_variable(i) <= L_micro+L_nano+L_macro ) 
  
 B1(i) = (q_s_macro*(t_new(i)-t_new(i-1)))/(2*mass_macro*p*Gamma_macro*L_macro*abs(x_variable(i)-x_0)^2);
 B2(i) = (Voltage(i)*(t_new(i)-t_new(i-1)))/(2*mass_macro*Gamma_macro*L_macro*abs(x_b2-x_b1));
 A_new(i) = B1(i)+B2(i);

 M(i) = (zeta_macro*v_new(i-1))/(Gamma_macro*L_macro);
 C1(i) = (q_s_macro*q_new(i-1)*(t_new(i)-t_new(i-1)))/(2*mass_macro*p*Gamma_macro*L_macro*abs(x_variable(i-1)-x_0)^2);
 C2(i) = (q_new(i-1)*Voltage(i-1)*(t_new(i)-t_new(i-1)))/(2*mass_macro*Gamma_macro*abs(x_b2-x_b1)*L_macro);
 E(i) = -1/(t_new(i)-t_new(i-1));
 
 B_new(i) = M(i)+C1(i)+C2(i)+E(i) ; 
 
 C_new(i) = q_new(i-1)/(t_new(i)-t_new(i-1)) ;
 
 % Charge of the particle in macropore
  q_new(i) = (-B_new(i) - sqrt((B_new(i)^2)-(4*A_new(i)*C_new(i))))/(2*A_new(i));
 
  % position of the particle
  x_variable(i) = x_variable(i-1) + (v_new(i-1)*(t_new(i)-t_new(i-1))) + (0.5*(acc_new(i-1)*(t_new(i)-t_new(i-1))^2));
 % velocity of particle
 v_new(i) = (zeta_macro*v_new(i-1))/Gamma_macro + (q_s_macro*q_new(i-1)*(t_new(i)-t_new(i-1)))/(2*mass_macro*p*Gamma_macro*abs(x_variable(i-1)-x_0)^2) + (q_new(i-1)*Voltage(i-1)*(t_new(i)-t_new(i-1)))/(2*mass_macro*Gamma_macro*abs(x_b2-x_b1)) + (q_s_macro*q_new(i)*(t_new(i)-t_new(i-1)))/(2*mass_macro*p*Gamma_macro*abs(x_variable(i)-x_0)^2) + (q_new(i)*Voltage(i)*(t_new(i)-t_new(i-1)))/(2*mass_macro*Gamma_macro*abs(x_b2-x_b1));
 % acceleration of the particle
 acc_new(i-1) = (q_s_macro*q_new(i-1))/(mass_macro*p*abs(x_variable(i-1)-x_0)^2) + (q_new(i-1)*Voltage(i-1))/(mass_macro*abs(x_b2-x_b1)) - (lambda*v_new(i-1))/(mass_macro);
 acc_new(i) = (q_s_macro*q_new(i))/(mass_macro*p*abs(x_variable(i)-x_0)^2) + (q_new(i)*Voltage(i))/(mass_macro*abs(x_b2-x_b1)) - (lambda*v_new(i))/(mass_macro);

 % Current in the macropore
   current(i) = (v_new(i)*q_new(i))/L_macro;
      
       IN1 = wiener2(current,[i i])-2875*wiener2(current,[i i])+wiener2(current,[i i])-3490*wiener2(current,[i i])-2840*wiener2(current,[i-20 i-20])+wiener2(current,[i-20 i-20])+500*wiener2(current,[i-20 i-20])-3470*wiener2(current,[i-80 i-80])-wiener2(current,[i-80 i-80])-2850*wiener2(current,[i-80 i-80])-2190*wiener2(current,[i-80 i-80])-wiener2(current,[i-80 i-80])-7480*wiener2(current,[i-180 i-180])-6870*wiener2(current,[i-180 i-180])-wiener2(current,[i-180 i-180])-8540*wiener2(current,[i-180 i-180])-7940*wiener2(current,[i-180 i-180])-17870*wiener2(current,[i-230 i-230])-wiener2(current,[i-230 i-230])-18540*wiener2(current,[i-230 i-230])-12350*wiener2(current,[i-230 i-230])-21870*wiener2(current,[i-50 i-50])-wiener2(current,[i-50 i-50])-25540*wiener2(current,[i-50 i-50])-18650*wiener2(current,[i-50 i-50])-7684*wiener2(current,[i-20 i-20])-wiener2(current,[i-20 i-20])-8165*wiener2(current,[i-20 i-20])-6576*wiener2(current,[i-20 i-20]);
       IN2 = wiener2(current,[i-40 i-40])+8800*wiener2(current,[i-40 i-40])+wiener2(current,[i-40 i-40])+9500*wiener2(current,[i-40 i-40])+wiener2(current,[i-40 i-40])+9500*wiener2(current,[i-40 i-40])+88*wiener2(current,[i-80 i-80])+wiener2(current,[i-80 i-80])+95*wiener2(current,[i-80 i-80])+wiener2(current,[i-80 i-80])+95*wiener2(current,[i-80 i-80])+wiener2(current,[i-40 i-40])+8800*wiener2(current,[i-40 i-40])+wiener2(current,[i-40 i-40])+9500*wiener2(current,[i-40 i-40])+wiener2(current,[i-40 i-40])+9500*wiener2(current,[i-40 i-40])+95*wiener2(current,[i-80 i-80])+wiener2(current,[i-40 i-40])+88*wiener2(current,[i-40 i-40])+wiener2(current,[i-180 i-180])+6800*wiener2(current,[i-180 i-180])+wiener2(current,[i-180 i-180])+9400*wiener2(current,[i-180 i-180])+21300*wiener2(current,[i-230 i-230])+wiener2(current,[i-230 i-230])+24800*wiener2(current,[i-230 i-230])+25700*wiener2(current,[i-40 i-40])+wiener2(current,[i-40 i-40])+27900*wiener2(current,[i-40 i-40])+4520*wiener2(current,[i-50 i-50])+2850*wiener2(current,[i-50 i-50]);
       
       INoise = IN1+IN2+IN1+IN2+IN1+IN2+IN1+IN2+IN1+IN2+IN1+IN2+IN1+IN2+IN1+IN2+IN1+IN2+IN1+IN2+IN1+IN2+IN1+IN2+IN1+IN2+IN1+IN2+IN1+IN2+IN1+IN2;

 
 
 
 
 
 
 end
 
 if (x_variable(i) > L_micro+L_nano+L_macro)
     size_plot = i;
     break;
 end
 
 
  if (x_variable(i) < 0)
     size_plot = i;
     break;
  end
  
  size_plot = i;
end

% PLOTTING VARIABLES
%     z
figure;plot(t_new(1:size_plot-1)*1e9,(x_variable(1:size_plot-1))*1e9,'-b');
xlabel('time (ns)');
ylabel('x (nm)');
% 
figure;plot(t_new(1:size_plot-1)*1e9,q_new(1:size_plot-1)*1e18,'b');
xlabel('time (ns)');
ylabel('q (aC)');
 
figure;plot(t_new(1:size_plot-2)*1e9,current(1:size_plot-2)*1e9,'b');
xlabel('time (ns)');
ylabel('I (nA)');

figure;plot(t_new(1:size_plot-2)*1e9,INoise(1:size_plot-2)*1e9,'b');
xlabel('time (ns)');
ylabel('I noise (nA)');

I_Mean = mean(INoise)
 
figure;plot(q_new(1:size_plot-1)*1e18,current(1:size_plot-1)*1e9,'b');
xlabel('q (aC)');
ylabel('(dq/dt) (nA)');
 
% figure;plot(t_new(1:size_plot-1)*1e9,potential(1:size_plot-1),'-ob');
% xlabel('time (ns)');
% ylabel('\phi (V)');
%  
% figure;plot(t_new(1:size_plot-1)*1e9,Voltage(1:size_plot-1),'-ob');
% xlabel('time (ns)');
% ylabel('Voltage (V)');