% Set format to long for higher precision
format long;

% Constants and Parameters
d_microns = 50;         % Thickness of glass in micrometers
d = d_microns * 1e-6;   % Thickness in meters
N = 500;                % Number of spatial points
L = d;                  % Length of the glass in meters
x = linspace(0, L, N);  % Spatial grid
dx = L / (N - 1);       % Grid spacing corrected

V0 = 20;                % Voltage applied across the glass (in volts)
e = 1.602e-19;          % Elementary charge in coulombs
epsilon_0 = 8.854e-12;  % Permittivity of free space (F/m)
er = 6;                 % Relative permittivity of the glass
k_B = 1.381e-23;        % Boltzmann constant 
E_a = 0.8 * e;          % Activation energy
mu_0 = 1e-5;            % Pre-exponential factor for mobility (m^2/V/s)
Tp = 250 + 273.15;      % Poling temperature in Kelvin

% Initial concentrations (in ions/m^3)
c1_0 = 4e27;
c2_0 = 2e27;
c3_0 = c1_0 + c2_0;

% Charges of the ionic species (z1, z2, z3)
z1 = 1;   % Charge of c1
z2 = 1;   % Charge of c2
z3 = -1;  % Charge of c3 (assuming it's negative for electroneutrality)

% Function to calculate mobility as a function of temperature
mu1_profile = mu_0 * exp(-E_a / (k_B * Tp));

% Diffusion coefficients using Einstein relationship
D1 = ((k_B * Tp) / e) * mu1_profile;

% Time parameters
dt = 1;              % Time step in seconds
T = 10;               % Total simulation time in seconds
Nt = ceil(T / dt);   % Number of time steps

% Initialize concentration profiles
c1_profile = zeros(N, Nt);
c2_profile = zeros(N, Nt);
c3_profile = zeros(N, Nt);

% Set initial concentrations
c1_profile(:, 1) = c1_0 * ones(N, 1);
c2_profile(:, 1) = c2_0 * ones(N, 1);
c3_profile(:, 1) = c3_0 * ones(N, 1);

% Initialize potential phi and electric field E
phi = zeros(N, Nt);
E_profile = zeros(N, Nt);

% Set initial potential profile
phi(:, 1) = linspace(V0, 0, N)';  % Linear gradient from V0 to 0

% Solver parameters
maxIter = 1000;       % Maximum iterations for potential solver
tol = 1e-4;           % Tolerance for potential convergence

% Time evolution loop
for t = 2:Nt
    % Solve Poisson equation in 1D using finite differences
    for iter = 1:maxIter
        % Update charge density rho based on current potential phi
        rho = z1 * e * c1_profile(:, t) + z2 * e * c2_profile(:, t) + z3 * e * c3_profile(:, t);
        
        % Save current phi for comparison
        phi_old = phi(:, t);
        
        % Update potential phi using finite difference method
        phi_new = phi(:, t);  % Start with the current phi values
        phi_new(2:end-1) = 0.5 * (phi(1:end-2, t) + phi(3:end, t) - rho(2:end-1) * dx^2 / (epsilon_0 * er));
        
        % Apply boundary conditions
        phi_new(1) = V0;
        phi_new(N) = 0;
        
        % Check for convergence
        if max(abs(phi_new - phi_old)) < tol
            phi(:, t + 1) = phi_new;
            break;
        end
    end
    
    % Update concentrations using drift-diffusion equations
    
    % Calculate electric field E as negative gradient of phi
    E_profile(:, t) = -gradient(phi(:, t), dx);
    
    % Calculate drift velocities (assuming linear drift)
    v1 = mu1_profile * E_profile(:, t);
    v2 = zeros(N, 1);  % Assume zero drift velocity for simplicity
    v3 = zeros(N, 1);  % Assume zero drift velocity for simplicity
    
    % Calculate spatial gradients for concentration
    dc1_dx = gradient(c1_profile(:, t), dx);
    dc2_dx = gradient(c2_profile(:, t), dx);
    dc3_dx = gradient(c3_profile(:, t), dx);

    % Diffusion flux based on Fick's law
    J_diffusion1 = -D1 .* dc1_dx;
    J_diffusion2 = zeros(N, 1);  % Assume zero diffusion for c2
    J_diffusion3 = zeros(N, 1);  % Assume zero diffusion for c3

    % Drift flux based on drift velocity
    J_drift1 = -z1 * e * v1 .* c1_profile(:, t);
    J_drift2 = zeros(N, 1);  % Assume no drift for c2
    J_drift3 = zeros(N, 1);  % Assume no drift for c3

    % Mass balance equation: d(c)/dt = -div(J)
    dc1_dt = -gradient(J_diffusion1 + J_drift1, dx);
    dc2_dt = zeros(N, 1);  % Assume no change for c2
    dc3_dt = zeros(N, 1);  % Assume no change for c3

    % Update concentrations using forward Euler method
    c1_profile(:, t + 1) = c1_profile(:, t) +dt .* dc1_dt;
    c2_profile(:, t + 1) = c2_profile(:, t) +dt .* dc2_dt;
    c3_profile(:, t + 1) = c3_profile(:, t) +dt .* dc3_dt;
end

% Plotting
figure;

subplot(3, 1, 1);
plot(x * 1e6, phi(:, :), 'LineWidth', 1.5);
xlabel('Position (\mum)');
ylabel('Electric Potential (V)');
title('Electric Potential Profile');
grid on;
legend(arrayfun(@(n) sprintf('t = %.2f s', n * dt), 0:dt:(T - dt), 'UniformOutput', false));

subplot(3, 1, 2);
plot(x * 1e6, E_profile(:, :), 'LineWidth', 1.5);
xlabel('Position (\mum)');
ylabel('Electric Field (V/m)');
title('Electric Field Profile');
grid on;

subplot(3, 1, 3);
plot(x * 1e6, c1_profile(:, :), 'r-', 'LineWidth', 2);
hold on;
plot(x * 1e6, c2_profile(:, :), 'b-', 'LineWidth', 2);
plot(x * 1e6, c3_profile(:, :), 'g-', 'LineWidth', 2);
hold off;
xlabel('Position (\mum)');
ylabel('Concentration (m^{-3})');
title('Final Concentration Profiles');
legend('c1', 'c2', 'c3', 'Location', 'best');
grid on;
