function [Potential] = PotentialProfileWithPoisson(x, Steps, V, VShift, VBias, NDoping, eps1, eps2)
% [Potential] = PotentialProfileWithPoisson(x, Steps, V, VShift, VBias, NDoping, eps1, eps2)
% 
% Description:
% This function creates a potential profile given a set of potential steps.
% The Poisson equation is solved given the material parameters so that the
% potential includes space charging effects.  It is assumed that the first
% barrier begins at x = 0 and that enough room has been left in the domain
% for the bending due to VBias.
%
% Outputs:
% Potential - Spatial vector containing the potential at each point of 'x'
%
% Inputs:
% x       - Spatial vector containing the points of the x-axis where the
%           potential is defined.
% Steps   - Vector containing the locations of each of the potential steps
% V       - Vector containing the height of the potential at each location
%           of Steps
% VShift  - Shifts the entire potential up or down
% VBias   - Bias voltage
% NDoping - Doping concentration of leads in m^-3
% eps1    - Permittivity of the electrode in F/m
% eps2    - Permittivity of the barrier region in F/m

%*************************
% Constants and Matrices *
%*************************
% Universal Constants
echarge = 1.6021764e-19;	%electron charge (C)

% Calculated Parameters
Numx = max(size(x));                % Number of Samples in x-Axis
NumBarrier = max(size(Steps)) - 1;  % Number of Barriers

d2 = Steps(end);        %total width of potential profile (m) 
d1 = -(eps1*d2)/(2.0*eps2)+sqrt(VBias/echarge*eps1/NDoping+(eps1*d2)/(2.0*eps2)*(eps1*d2)/(2.0*eps2));
                            %width of charge accumulation/depletion region
q = VBias/(d1/eps1+d2/eps2);%precalculate parameter
V1 = q*d1/(2.0*eps1);       %depth of charge accumulation/depletion well
V2 = q*d2/eps2;             %potential decrease along potential profile

%***************************
% Create Potential Profile *
%***************************
for xIndex = 1:Numx
    if x(xIndex) < -d1
        Potential(xIndex) = 0.0;
    end
    if x(xIndex) >= -d1 
        if x(xIndex) < 0  
            Potential(xIndex) = -q*(x(xIndex)+d1)^2/(2.0*eps1*d1);   
            %form of potential profile in the negative charge accumulation 
            %region at the left electrode
        end
    end
    if x(xIndex) >= 0.0
        if x(xIndex) < d2
            for BarrIndex = 1:NumBarrier 
                if x(xIndex) >= Steps(BarrIndex)
                    if x(xIndex) < Steps(BarrIndex + 1)
                        Potential(xIndex) = V(BarrIndex); 
                    end 
                end 
            end
            Potential(xIndex) = Potential(xIndex) - V1 - V2*x(xIndex)/d2;
            %potential decrease along the barrier
        end 
    end
    if x(xIndex) >= d2
        if x(xIndex) < d2 + d1
            Potential(xIndex) = -VBias + q*(x(xIndex) - (d2 + d1))^2/(2.0*eps1*d1);
            %form of potential profile in the positive charge accumulation
            %region at the right electrode
        end 
    end
    if x(xIndex) >= d2 + d1
        Potential(xIndex) = -VBias;
    end
end

% Offset Potential Profile
Potential = Potential + VShift;
end