#######  Script to generate radial doping profile  ####

# Inputs

ND = 1e19;  # donor concentration; unit /cm^3.  Use negative value for acceptors
P = 401;  # number of points along x and y direction
Pz = 21;  # number of points along z direction

# rectangular grid

x = linspace(-1.51,1.51,P)*1e-6;
y = linspace(-1.51,1.51,P)*1e-6;
z = linspace(-2.1,2.1,Pz)*1e-6;  # 

N = matrix(P,P,Pz);

N(1:P,1:P,1:Pz) = ND * 1e6;  # Set doping at all points to ND.  Make unit /m^3

R_bottom = 1.5e-6;  # cone radius at base 
R_top = 1e-6;  # cone radius at top
depth = 0.5e-6;  # dopant penetration

R = linspace(R_bottom,R_top,Pz) - depth;  # inner radius at base and at top

for (k=1:Pz) {
    for (i=1:P) {
        for (j=1:P) {
            if (sqrt(x(i)^2+y(j)^2)<R(k)) {N(i,j,k) = 0;}
        }
    }
}

# plot profile on XY plane

image(x,y,N(1:P,1:P,1),'x (m)','y (m)','N (/m^3)');
image(x,y,N(1:P,1:P,11),'x (m)','y (m)','N (/m^3)');
image(x,y,N(1:P,1:P,21),'x (m)','y (m)','N (/m^3)');

# save data in .mat file

matlabsave("charge",N,x,y,z);