clear all
close all
NApmax = 1e17;
NDn = 1e15;
NDnmax = 1e17;
ju=1e-6;
ac = 2*10^-6;
y1 = 7*10^-6;
x1 = 5*10^-6;
ch = ju/sqrt(log(NApmax/NDn));
yy =  0:1e-7:y1;
xx = 0:1e-7:x1;
ni = 4.455e-11;
%ni=1.46e10;
i = 1;
j = 1;
q = 1.602e-19;
k = 1.38e-23;
T = 300;
c = q/(k*T);
for y = 0:1e-7:y1
    j=1;
    for x = 0:1e-7:x1
        if x < ac
            N(i,j) = NDn + NDnmax.*exp(-((y+y1)./ch).^2)-NApmax.*exp(-(y./ch).^2);
        else
           N(i,j) = NDn + NDnmax.*exp(-((y+y1)./ch).^2)-NApmax.*exp(-(y./ch).^2).*exp(-((x-ac)./ch).^2);
        end
        j = j +1;
    end
    i=i+1;
end
for y = 1:length(yy)
    for x = 1:length(xx)
        if(N(y,x)>=0)
            n_init(y,x)= (abs(N(y,x))/2+sqrt(N(y,x)^2/4+ni^2));
            p_init(y,x) = ni^2/(abs(N(y,x)/2+sqrt(N(y,x)^2/4+ni^2)));
            psi_init(y,x) = log(n_init(y,x)/ni);
        else
            n_init(y,x) = ni^2/(abs(N(y,x))/2+sqrt(N(y,x)^2/4+ni^2));
            p_init(y,x) = (abs(N(y,x))/2+sqrt(N(y,x)^2/4+ni^2));
            psi_init(y,x) = 1/c*(-log(p_init(y,x)/ni));
        end
    end
end
surf(xx,yy,N)
title 'N'
figure
surf(xx,yy,n_init)
title 'n_init'
figure
surf(xx,yy,p_init)
title 'p_init'
figure
surf(xx,yy,psi_init)
title 'psi_init'