par (mar=c(4,4,1,1))

HV <- 1500	# V
Th <- 1	# cm
Na <- 1e10	# cm-3
e <- 1.6e-19
eps <- 16 *8.8e-12 #F/m

mu.e <- 3900 #cm2/V/s
mu.h <- 1900 #cm2/V/s

#returns E in V/cm
E <- function (x) {
	ifelse (x < 0, 0, ifelse (x > Th, 0, HV/Th+e*Na/eps*(Th/2-x)*100))
}
v.e <- function (x) mu.e*E(x)
v.h <- function (x) mu.h*E(x)

dt <- 0.1	# ns

t <- seq(0, 1000, by=dt)

pos.vs.time <- function (x0, q) {
	x <- rep(0,length(t))
	x[1] <- x0
	i <- 1
	repeat {
		if (q=="e") v <- -v.e(x[i])
		else v <- v.h(x[i])
		x[i+1] <- x[i] + v * (dt * 1e-9)
		i <- i + 1
		if (v == 0) break
	}
	if (q=="e") x[i:length(x)] <- -0.001
	else x[i:length(x)] <- Th + 0.001
	x
} 

current.vs.time <- function (x, q) {
	if (q=="e") v <- v.e(x)
	else v <- v.h(x)

	i <- v * 1/Th * 40666 * 1.6e-19
}


x0 <- 0.1	# cm

##
## position vs time
##

x.e <- pos.vs.time (x0, "e")
x.h <- pos.vs.time (x0, "h")

# position
plot (c(0,400), c(0,Th), type="n",
	xlab="time (ns)", ylab="position(cm)")
grid()
lines (t, x.e, col="blue", lwd=2)
lines (t, x.h, col="red", lwd=2)

#velocity
plot (c(0,400), c(0,80), type="n",
	xlab="time (ns)", ylab="velocity (mm/µs)")
grid()
lines (t, v.e(x.e) * 1e-6*10, col="blue", lwd=2)
lines (t, v.h(x.h) * 1e-6*10, col="red", lwd=2)

i.e <- current.vs.time (x.e, "e")
i.h <- current.vs.time (x.h, "h")

sum(i.e) * dt * 1e-9
sum(i.h) * dt * 1e-9

#current (A)
plot (c(0,400), c(0,10e-8), type="n",
	xlab="time (ns)", ylab="current (A)")
grid()
lines (t, i.e, col="blue", lwd=2)
lines (t, i.h, col="red", lwd=2)
lines (t, i.e+i.h+0.04e-8, col="green", lwd=2)

#charge (C)
plot (c(0,400), c(0,7e-15), type="n",
	xlab = "time(ns)", ylab="charge")
grid()
lines (t, cumsum(i.e + i.h) * dt * 1e-9, col="green", lwd=2)