%PotentialProfileWithPoisson.m
%calculates potential profile due to charge
%accumulation/depletion at the electrodes by solving Poisson equation
%in the depletion approximation
%Input: Potential profile (V(i), i = 1..NBarrierStep) build out of NBarrierStep barriers
%Output: Potential profile (Vpot(i), i = 1 .. xpoints), including Poisson

function [Vpot] = PotentialProfileWithPoisson(xpoints, x, dx, Lx, NBarrier, V, Vbias, NDoping, eps1, eps2);
%NDoping=NDoping*100;        %*100 to fix potential profile to look like Tsu
echarge = 1.6021764e-19;	%electron charge (C)

d2 = Lx(NBarrier+1);        %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
    

for i = 1 : 1 : xpoints
    if x(i) < -d1
        Vpot(i) = 0.0;
    end
    if x(i) >= -d1 
        if x(i) < 0  
            Vpot(i) = -q/(2.0*eps1*d1)*(x(i)+d1)*(x(i)+d1);   
            %form of potential profile in the negative charge accumulation 
            %region at the left electrode
        end
    end
    if x(i) >= 0.0
        if x(i) < d2
            for k = 1 : 1 : NBarrier 
                if x(i) >= Lx(k)
                    if x(i) < Lx(k+1)
                        Vpot(i) = V(k); 
                    end 
                end 
            end
            Vpot(i) = Vpot(i) - V1 - V2/d2*x(i);
            %potential decrease along the barrier
        end 
    end
    if x(i) >= d2
        if x(i) < d2+d1
            Vpot(i) = -Vbias + q/(2.0*eps1*d1)*(x(i)-(d2+d1))*(x(i)-(d2+d1));
            %form of potential profile in the positive charge accumulation
            %region at the right electrode
        end 
    end
    if x(i) >= d2+d1
        Vpot(i) = -Vbias;
    end
end