##  This script creates a radial diffusion doping profile.  The doping is applied
##  in the region between radius_out and radius_in (see definition below)
##  Inside radius_in the doping is set to a low value (Nref)


#### INPUTS  #####
#
#  Npeak = peak doping concentration in /m^3
#  Nref = reference (minimum) doping concentration in /m^3
#  radius_out = outer radius of the geometry.  Make this slightly larger than the actual geometry to ensure doping gets applied to the surface properly
#  radius_in = inner radius of the doped portion of the geometry
#  junction width = distance over which the concentration goes from peaf to ref.
#  zspan = length of the doping profile in z direction
#  x0, y0, z0 = x, y, and z coordinate of the center of the doping profile
#  P = number of points in rectangular grid.  Note a larger value of P will result in a large matrix and a long calculation time.
#
###################

# Radial doping
clear;

## Set Inputs

Npeak = 1e21;  
Nref = 1e12;  
radius_out = 1.1e-6; 
radius_in = 0.125e-6; 
junction_width = .05e-6; 
x0 = 0;  
y0 = 0;  
z0 = 0.125e-6;  
zspan = 5.25e-6; 
P = 400; 

##

#  Create rectangular grid

x = linspace(x0-radius_out,x0+radius_out,P);
y = linspace(y0-radius_out,y0+radius_out,P);
z = [z0-zspan/2, z0+zspan/2];

#  Initialize matrices

N2d = matrix(P,P) + Npeak;
N3d = matrix(P,P,2) + 1;
R = linspace(radius_out,0,P);

#  Locate junction

Rj = find( (R<(junction_width+radius_in)) & (R>radius_in) );
Nr_peak = matrix(Rj(1)-1,1);
R_erfc = R(Rj);
temp = length(Rj);
Nr_ref = matrix(P-Rj(temp),1) - (Npeak - Nref);
kern = length(R_erfc);
end = kern;

#  Create vector to model doping variation inside junction

R_erfc = R_erfc - R_erfc(1);
R_erfc = (R_erfc / R_erfc(end) * 5) - 2.5;
Nr_erfc = erfc(R_erfc);
Nr_erfc = Nr_erfc - Nr_erfc(1);
Nr_erfc = Nr_erfc / (-Nr_erfc(end));
Nr_erfc = Nr_erfc * (Npeak - Nref);

Nr = [Nr_peak; Nr_erfc; Nr_ref];
#R = radius - R;

#  Generate 2D doping profile

for (i=1:P) {
    for (j= 1:P) {
        Rcheck = sqrt( (x(i) - x0)^2 + (y(j) - y0)^2 );
        kern = find(R<=Rcheck);
        N2d(i,j) = N2d(i,j) + Nr(kern(1));
    }
}

# Create 3D doping profile

N3d(1:P,1:P,1) = N2d;
N3d(1:P,1:P,2) = N2d;
N = N3d;

#  Create rectilinear dataset

doping = rectilineardataset(x,y,z);
doping.addattribute("N",N3d);

#  Save Matlab file

matlabsave("rad_doping_dataset",doping);
matlabsave("rad_diff_doping",x,y,z,N);