from NanoTCAD_ViDES import *

def charge_SC(self):
    vt=self.Temp*kboltz/q;
    dist=avervect(self.x)
    self.Ei=-self.Phi;
    #Electrons
    for i in range(0,size(self.charge)):
        self.charge[i]=-6*dist[i]*1e-9*(2/sqrt(pi))*2*(vt/(2*pi)*(self.massl[i]*m0/hbar)*(q/hbar))**1.5*fphalf(-(self.Ei[i]+self.Egap[i]*0.5-self.Ef)/vt)
    #Holes
    for i in range(0,size(self.charge)):
        self.charge[i]=self.charge[i]+dist[i]*1e-9*(2/sqrt(pi))*2*(vt/(2*pi)*(self.massh[i]*m0/hbar)*(q/hbar))**1.5*fphalf((self.Ei[i]-self.Egap[i]*0.5-self.Ef)/vt)

# I create the grid 1D
x=nonuniformgrid(array([0,0.1,2,0.05,102,10]))

grid=grid1D(x,x)  

top_gate=gate("hex",grid.xmin,grid.xmin+0.05);
top_gate.Ef=-2.0; 
bottom_gate=gate("hex",grid.xmax-1,grid.xmax);

SiO2=region("hex",grid.xmin,2);
SiO2.set_material("SiO2")

Si=region("hex",2,grid.xmax);
Si.rho=-1e24;
Si.set_material("Si")

p=interface1D(grid,top_gate,bottom_gate,SiO2,Si);
p.normd=1e-7

nx=size(x);
# Uncomment line 37 for semiclassical approximation 
# both for holes and electrons and comment
# line 39
#M=QM1D(nx,3,x,p,charge_SC);

M=QM1D(nx,3,x,p);

solve_self_consistent(grid,p,M);

b=[x,M.Psi[:,0]];
savetxt("Psi1.out",transpose(b));

dist=avervect(x)*1e-9;
a=[x,M.charge/dist];
savetxt("free_charge.out",transpose(a));

a=[x,p.Phi];
savetxt("Phi.out",transpose(a));

a=[x,M.Ei+M.Egap*0.5];
savetxt("Ec.out",transpose(a));

a=[x,M.Ei];
savetxt("Ei.out",transpose(a));

a=[x,p.fixed_charge];
savetxt("rho2.out",transpose(a));